C++ TH2I::Write方法代码示例

//.........这里部分代码省略.........
	  if (verbose>5) {
	    std::cout<<"MISMATCH"<<std::endl ;
	    std::cout<<"(phi,eta, Nbxtals)="<<std::dec<<iphi<<" "<<ieta<<" "<<nbXtals<<std::endl ;
	    std::cout<<"Data Et, TTF: "<<tp<<" "<<ttf<<std::endl ;
	    std::cout<<"Emulator: " ;
	    for (int i=0 ; i<5 ; i++) std::cout<<emul[i]<<" " ;
	    std::cout<<std::endl ;
	  }
	}
      }
      if (maxOfTPEmul>0) TPEmulMaxIndex->Fill(indexOfTPEmulMax+1) ;


      // Fill TTF mismatch
      if ((ttf==1 || ttf==3) && nbXtals != 25) ttfMismatch->Fill(iphi, ieta) ;


    } // end loop towers


  } // endloop entries

  

  ///////////////////////
  // Format & write histos
  ///////////////////////


  // 1. TP Spectrum  
  TProfile2D * TPspectrumMap = TPspectrumMap3D->Project3DProfile("yx") ;
  TPspectrumMap->SetName("TPspectrumMap") ;

  // 2. TP Timing
  TH2F * TPMatchEmul2D = new TH2F("TPMatchEmul2D", "TP data matching Emulator", 72, 1, 73, 38, -19, 19) ;
  TH2F * TPMatchFraction2D = new TH2F("TPMatchFraction2D", "TP data: fraction of non-single timing", 72, 1, 73, 38, -19, 19) ;
  TPMatchEmul2D->GetYaxis()->SetTitle("eta index") ; 
  TPMatchEmul2D->GetXaxis()->SetTitle("phi index") ;
  TPMatchEmul2D->GetZaxis()->SetRangeUser(-1,6) ;
  TPMatchFraction2D->GetYaxis()->SetTitle("eta index") ; 
  TPMatchFraction2D->GetXaxis()->SetTitle("phi index") ;
  for (int binx=1 ; binx<=72 ; binx++)    
    for (int biny=1 ; biny<=38 ; biny++) {
      int maxBinz = 5 ;
      double maxCell = TPMatchEmul3D->GetBinContent(binx, biny, maxBinz) ;
      double totalCell(0) ;
      for (int binz=1; binz<=7 ; binz++) {
	double content = TPMatchEmul3D->GetBinContent(binx, biny, binz) ;
	if (content>maxCell) {
	  maxCell = content ;
	  maxBinz = binz ;
	}
	totalCell += content ;
      }
      if (maxCell <=0) maxBinz = 2 ; // empty cell
      TPMatchEmul2D->SetBinContent(binx, biny, float(maxBinz)-2.) ; //z must be in [-1,5] 
      double fraction = 0 ;
      if (totalCell>0) fraction = 1.- maxCell/totalCell ;
      TPMatchFraction2D->SetBinContent(binx, biny, fraction) ;
      if (totalCell > maxCell && verbose>9) {
	std::cout<<"--->"<<std::endl ;	
	for (int binz=1; binz<=7 ; binz++) {	  
	  std::cout<< "(phi,eta, z): (" 
		   << TPMatchEmul3D->GetXaxis()->GetBinLowEdge(binx) 
		   << ", " << TPMatchEmul3D->GetYaxis()->GetBinLowEdge(biny) 
		   << ", " << TPMatchEmul3D->GetZaxis()->GetBinLowEdge(binz)		   
		   << ") Content="<<TPMatchEmul3D->GetBinContent(binx, biny, binz)		   
		   << ", erro="<<TPMatchEmul3D->GetBinContent(binx, biny, binz)	   
		   << std::endl ;	
	}
      }
    }



  TFile saving (outputRootName.c_str (),"recreate") ;
  saving.cd () ;
  
  occupancyTP->Write() ;
  occupancyTPEmul->Write() ;
  
  TP->Write() ;
  TPEmul->Write() ;
  TPEmulMax->Write() ;
  TPspectrumMap->Write() ;

  TPMatchEmul->Write() ; 
  TPMatchEmul3D->Write() ; 
  TPEmulMaxIndex->Write() ;
  TPMatchEmul2D->Write() ; 
  TPMatchFraction2D->Write() ; 

  ttfMismatch->Write() ; 

     
  saving.Close () ;
  delete chain ;

  return 0 ;
}

//.........这里部分代码省略.........
    gStyle->SetPalette(1);
    while(s.Next()){                     //read the next raw digit; returns kFALSE if there is no digit left
      if(s.IsCompletedModule()) continue;
      if(s.IsCompletedDDL()) continue;
      iddl=rd->GetDDLID()-2; // -2 is temporary for test raw data

	isddmod=s.GetModuleNumber(iddl,s.GetCarlosId());        //this is the FEE Carlos
       	//cout<<"DDLID= "<<iddl <<"; Module number= " <<isddmod  <<endl;
	if(isddmod >= imodoffset) { 
	  hModulePattern->Fill(isddmod-imodoffset);             // 0 to 259    so  240 to 499
	  moduleSDD=2*(isddmod-imodoffset)+s.GetChannel();
          hModuleSidePattern->Fill(moduleSDD);                  // 0 to 519
	  hModuleCarlos->Fill(isddmod-imodoffset,s.GetCarlosId());
          //cout << "anode " << s.GetCoord1() << ", time bin: " << s.GetCoord2() << ", charge: " << s.GetSignal() << endl;	  
	  Int_t coord1 = s.GetCoord1();
 	  Int_t coord2 = s.GetCoord2();
 	  Int_t signal = s.GetSignal();
	  hModuleChargeMap[moduleSDD]->Fill(coord2, coord1,signal);
          hModuleCountsMap[moduleSDD]->Fill(coord2, coord1 );   
	  //hModuleSideBL[moduleSDD][coord1]->Fill(signal);             // BL  !!!
	} else {
	  nSkip++;
	}
    }    
    cout << "End of Event " << evCounter+FirstEvt-1 << ", " << nSkip << " wrong module numbers" << endl;
  } while(rd->NextEvent()); // end loop on events
  delete rd;

  cout << "end after " << evCounter << " events" << endl;
  /*
  TNtuple *Baseline = new TNtuple("Baseline","Baseline","HalfModule:Anode:Mean:RMS");      // BL!!!
  Float_t meanBL;
  Float_t rmsBL;  
  */
  for(Int_t i=0; i<modtotSDD; i++){   
    if(hModuleSidePattern->GetBinContent(i+1)){              //check if they're not empty
      hModuleChargeMap[i]->GetXaxis()->SetTitle("Time Bin");
      hModuleChargeMap[i]->GetYaxis()->SetTitle("Anode"); 
      hModuleCountsMap[i]->GetXaxis()->SetTitle("Time Bin");
      hModuleCountsMap[i]->GetYaxis()->SetTitle("Anode");  
      /*
      for(Int_t ianode=0; ianode<anode; ianode++ ){                                      // BL!!!
	hModuleSideBL[i][ianode]->GetXaxis()->SetTitle("ADC counts");
	hModuleSideBL[i][ianode]->GetYaxis()->SetTitle("#"); 
	meanBL = hModuleSideBL[i][ianode]->GetMean();
	rmsBL = hModuleSideBL[i][ianode]->GetRMS();
	gaussfitBL = hModuleSideBL[i][ianode]->Fit("gaus");	
	Baseline->Fill(i,ianode,meanBL,rmsBL);
      }
      */
    }
  }  
  
  hModuleSidePattern->GetXaxis()->SetTitle("2*(Module Number-1)+Side");
  hModuleSidePattern->GetYaxis()->SetTitle("Counts");  
  hModulePattern->GetXaxis()->SetTitle("Module Number");  
  hModulePattern->GetYaxis()->SetTitle("Counts");  


  //-------store Histograms
  cout << "Store Histograms" << endl;
  TString oFileName(iFile);
  oFileName.Append(".root");
  TFile *oFile = TFile::Open(oFileName,"recreate");
  hModulePattern->Write();
  hModuleSidePattern->Write();
  hModuleCarlos->Write();
  for(Int_t i=0; i<modtotSDD; i++){ 
      if(hModuleSidePattern->GetBinContent(i+1)){            //check if they're not empty
	hModuleChargeMap[i]->Write();
	hModuleCountsMap[i]->Write();     
	/* 
	for(Int_t ianode=0; ianode<anode; ianode++ ){                           // BL!!!
 	  hModuleSideBL[i][ianode]->Write();
	  Baseline->Write();
 	}
	*/
      }
  }
  
  oFile->Close();
  cout << "Clear memory" << endl;
  for(Int_t imod=0; imod<nSDDmodules;imod++){
    for(Int_t isid=0;isid<2;isid++){
      Int_t index=2*imod+isid;       //260*2 position
      delete hModuleChargeMap[index];
      delete hModuleCountsMap[index];       
      /*
      for(Int_t ianode=0; ianode<anode; ianode++ ){                              // BL!!!
	delete hModuleSideBL[index][ianode]; 
	delete Baseline;
      }
      */
    }
  }
  delete hModulePattern;
  delete hModuleSidePattern;
  delete hModuleCarlos;

}

//.........这里部分代码省略.........

  // choose loop over l1 tree first (smaller)
  //std::cout << "Begin making map." << std::endl;
  Long64_t l_entries = l1Tree->GetEntries();
  for(Long64_t j = 0; j < l_entries; ++j)
  {
    l1Tree->GetEntry(j);
    Long64_t key = makeKey(l1_run, l1_event);

    std::pair<Long64_t,Long64_t> p(key,j);
    kmap.insert(p);
  }
  //std::cout << "Finished making map." << std::endl;

  outFile->cd();

  TH1I *hists[22][18]; // [eta][pu bin], arbitrary value of 18 for # bins in pu
  for(int i = 0; i < 22; ++i)
    for(int j = 0; j < 18; ++j)
    {
      hists[i][j] = new TH1I(Form("hist_%d_%d",i,j),"", 1024,0,1024);
    }

  TH2I *centPUM = new TH2I("cenPUM","",200,0,200,396,0,396);
  int count = 0;

  Long64_t entries = fEvtTree->GetEntries();
  for(Long64_t j = 0; j < entries; ++j)
  {
    //if(j % 10000 == 0)
    //  printf("%lld / %lld\n",j,entries);

    fEvtTree->GetEntry(j);
    Long64_t key = makeKey(f_run, f_evt);

    std::map<Long64_t,Long64_t>::const_iterator got = kmap.find(key);
    if(got == kmap.end() ) {
      continue;
    } else {
      l1Tree->GetEntry(got->second);
      kmap.erase(key);
      count++;

      fSkimTree->GetEntry(j);
      if((pcollisionEventSelection == 1) && (pHBHENoiseFilter == 1))
      {

	//int pubin = (int) ( (double)hiBin * (18.0/200.0));
	int PUM0 = 0;
	for(int i = 0; i < NREG; ++i)
	{
	  if(region_hwPt[i] > 0)
	    ++PUM0;
	}
	int pubin = PUM0/22;
	if(pubin == 18) pubin = 17; //special case for every region firing
	for(int i = 0; i < NREG; ++i)
	{
	  hists[region_hwEta[i]][pubin]->Fill(region_hwPt[i]);
	}
	centPUM->Fill(hiBin,PUM0);

      }
    }
  }

  std::cout << "cms.vdouble(";
  TH1D *hists_eta[22];
  for(int i = 0; i < 22; ++i)
  {
    hists_eta[i] = new TH1D(Form("hists_eta_%d",i),"",18,0,17);
    for(int j = 0; j < 18; ++j)
    {
      double Mean = hists[i][j]->GetMean();
      double MeanError = hists[i][j]->GetMeanError();
      hists_eta[i]->SetBinContent(j,Mean);
      hists_eta[i]->SetBinError(j,MeanError);
      std::cout << Mean*0.5;
      if(!((i == 21) && (j == 17))) std::cout << ", ";
    }
  }
  std::cout << ")" << std::endl;

  for(int i = 0; i < 22; ++i)
  {
    hists_eta[i]->Write();
    for(int j = 0; j < 18; ++j)
    {
      hists[i][j]->Write();
    }
  }
  centPUM->Write();


  //std::cout << "Matching entries: " << count << std::endl;

  lFile->Close();
  fFile->Close();
  outFile->Close();
}