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

//.........这里部分代码省略.........
			chamber2->GetYaxis()->SetTitleOffset(1.4);
			chamber2->GetYaxis()->CenterTitle();

			pc2->cd(3);
			chamber3->Draw();
			chamber3->GetXaxis()->SetTitle("X position (px)");
			chamber3->GetXaxis()->CenterTitle();

			chamber3->GetYaxis()->SetTitle("Intensity");
			chamber3->GetYaxis()->SetTitleOffset(1.4);
			chamber3->GetYaxis()->CenterTitle();

			pc2->cd(4);
			chamber4->Draw();
			chamber4->GetXaxis()->SetTitle("X position (px)");
			chamber4->GetXaxis()->CenterTitle();

			chamber4->GetYaxis()->SetTitle("Intensity");
			chamber4->GetYaxis()->SetTitleOffset(1.4);
			chamber4->GetYaxis()->CenterTitle();

			//TFile myF("trackTree.root","RECREATE");
			

			//loop over all data in chain
			Int_t nevent = chain->GetEntries();	//get the number of entries in the TChain
			for (Int_t i=0;i<nevent;i++) {
				chain->GetEntry(i);

				//put the frame data into the histogram for this event
				for(int x=0;x<width;x++){
					for(int y=0;y<height;y++){
						if(frame[x][y]>0){
							frameHisto->Fill(x,y,frame[x][y]);
							if(y>580 && y<610){
								chamber1->Fill(x,frame[x][y]);
								chamber1y->Fill(y,frame[x][y]);
							}
							else if(y>400 && y<440){
								chamber2->Fill(x,frame[x][y]);
								chamber2y->Fill(y,frame[x][y]);
							}
							else if(y>240 && y<280){
								chamber3->Fill(x,frame[x][y]);
								chamber3y->Fill(y,frame[x][y]);
							}
							else if(y>50 && y<100){
								chamber4->Fill(x,frame[x][y]);
								chamber4y->Fill(y,frame[x][y]);
							}
						}	
					}
				}

				double x12[2];
				double y12[2];
				double x34[2];
				double y34[2];

				x12[0] = chamber1->GetMean();
				y12[0] = chamber1y->GetMean();			//593.3;
				x12[1] = chamber2->GetMean();
				y12[1] = chamber2y->GetMean();			//424.7;
				x34[0] = chamber3->GetMean();
				y34[0] = chamber3y->GetMean();			//262.5;
				x34[1] = chamber4->GetMean();

void ITSSDDQAMaker(char *iFile, Int_t MaxEvts=1000000, Int_t FirstEvt=0) {

  //To have Baseline Histos uncomment parts with " // BL!!! " comment

cout << "SDD Quality Assurance Prototype Macro" << endl; 

const Int_t nSDDmodules= 260;
const Int_t imodoffset = 240;
const Int_t modtotSDD  = nSDDmodules*2;
const Int_t anode = 256;

Float_t xi = -0.5;
Float_t xf = xi + nSDDmodules;
TH1F *hModulePattern = new TH1F("hModulePattern","Modules pattern",nSDDmodules,xi,xf);
xf = xi + modtotSDD;
TH1F *hModuleSidePattern = new TH1F("hModuleSidePattern","Modules/Side pattern",modtotSDD,xi,xf);

TH2F *hModuleChargeMap[modtotSDD];  //260 dx e 260 sx  with A, T, Q
TH2F *hModuleCountsMap[modtotSDD];  //260 dx e 260 sx  with A, T, Ncounts
TH2I *hModuleCarlos = new TH2I("hModuleCarlos","hModuleCarlos",modtotSDD,xi,xf,101,-0.5,100.5);
/*
  TH1F *hModuleSideBL[modtotSDD][anode];                                // BL!!!
*/ 


//-------histograms definition 
Char_t *hisnam = new Char_t[50];
Char_t *histit = new Char_t[50];  
Char_t *hisnam2 = new Char_t[50];
Char_t *histit2 = new Char_t[50];    
Char_t *hisnam3 = new Char_t[50];
 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

     sprintf(hisnam,"chargeMap%d",index);
     sprintf(histit,"Total Charge, module number %d",index);
     hModuleChargeMap[index]=new TH2F(hisnam,histit,256,-0.5,255.5,256,-0.5,255.5);  

     sprintf(hisnam2,"countsMap%d",index);
     sprintf(histit2,"Number of Counts, module number %d",index);
     hModuleCountsMap[index] = new TH2F(hisnam2,histit2,256,-0.5,255.5,256,-0.5,255.5);
     /*
       for(Int_t ianode=0; ianode<anode; ianode++){                         // BL!!!
       sprintf(hisnam3,"BL_module_%d_%d",index,ianode); 
       //cout<<hisnam3 <<endl;
       hModuleSideBL[index][ianode] = new TH1F(hisnam3,hisnam3,256,0.,1024.);
     }
     */
   }
 }

  TString strFile = iFile;
  strFile += "?EventType=7";
  AliRawReader *rd = new AliRawReaderDate(strFile.Data(),FirstEvt);  // open run
  Int_t evCounter = 0;
  do{                       // start loop on events
    if(++evCounter > MaxEvts) { cout << MaxEvts << " events read, stop" << endl; evCounter--; break; }  
    cout << "Read Event: " << evCounter+FirstEvt-1 << endl;

    rd->RequireHeader(kFALSE);             
    rd->Reset();                           // reset the current position to the beginning of the event
 
    Int_t nSkip = 0;                     // number of skipped signals
    AliITSRawStreamSDD s(rd);            //This class provides access to ITS SDD digits in raw data.
    Int_t iddl;
    Int_t isddmod;
    Int_t moduleSDD;
    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;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    
             
    // loop on towers
    for (unsigned int tower = 0 ; tower < treeVars.nbOfTowers ; tower++) {

      int tp = getEt(treeVars.rawTPData[tower]) ;
      int emul[5] = {getEt(treeVars.rawTPEmul1[tower]),
		     getEt(treeVars.rawTPEmul2[tower]),
		     getEt(treeVars.rawTPEmul3[tower]),
		     getEt(treeVars.rawTPEmul4[tower]),
		     getEt(treeVars.rawTPEmul5[tower])} ;
      int maxOfTPEmul = 0 ;
      int indexOfTPEmulMax = -1 ;
      for (int i=0 ; i<5 ; i++) if (emul[i]>maxOfTPEmul) {
	maxOfTPEmul = emul[i] ; 
	indexOfTPEmulMax = i ;
      }
      int ieta = treeVars.ieta[tower] ;
      int iphi = treeVars.iphi[tower] ;
      int nbXtals = treeVars.nbOfXtals[tower] ;
      int ttf = getTtf(treeVars.rawTPData[tower]) ;


      if (verbose>9 && (tp>0 || maxOfTPEmul>0)) {
	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 ;
      }


      // Fill TP spctrum
      TP->Fill(tp) ;
      TPEmul->Fill(emul[ref]) ;
      TPEmulMax->Fill(maxOfTPEmul) ;
      TPspectrumMap3D->Fill(iphi, ieta, tp) ;


      // Fill TP occupancy
      if (tp>occupancyCut) occupancyTP->Fill(iphi, ieta) ;
      if (emul[ref]>occupancyCut) occupancyTPEmul->Fill(iphi, ieta) ;


      // Fill TP-Emulator matching
      // comparison is meaningful when:
      if (tp>0 && nbXtals == 25) {
	bool match(false) ;
	for (int i=0 ; i<5 ; i++) {
	  if (tp == emul[i]) {
	    TPMatchEmul->Fill(i+1) ;
	    TPMatchEmul3D->Fill(iphi, ieta, i+1) ;
	    match = true ;
	  }
	}
	if (!match) {
	  TPMatchEmul->Fill(-1) ;
	  TPMatchEmul3D->Fill(iphi, ieta, -1) ;
	  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 ;
	  }

int frameStack2_Mall(char* arg){
	//Take the arguments and save them into respective strings
	std::string infileName, outfileName0, outAllfileName0, outfileName1, outAllfileName1;
	std::string inF, outF0, outF1, outAll0, outAll1;
	std::string inPrefix, outPrefix;
	std::string runs, layers;
	std::string runCount;
	std::istringstream stm(arg);

	inPrefix = "/home/p180f/Do_Jo_Ol_Ma/Analysis/MainProcedure/testMain/rawRoot/";
	outPrefix = "/home/p180f/Do_Jo_Ol_Ma/Analysis/MainProcedure/testMain/images/";

	outAll0 = "sliceXCuts_allLayers.png";
	outAllfileName0 = outPrefix + outAll0;
	std::cout << outAll0 << " created\n";
	outAll1 = "projYCuts_allLayers.png";
	outAllfileName1 = outPrefix + outAll1;
	std::cout << outAll1 << " created\n";

	const int width=480;	//width of the raw image
	const int height=640;	//height of the raw image
	
	TH2I *frameHistoAll = new TH2I("frameHistoAll","Stacked Frames After Edge Cuts",width/4,0,width,height/4,0,height);	//histogram for the stacked images
	TH1I *chamber1All = new TH1I("chamber1All","Chamber 1 After Edge Cuts",width/4,0,width);//histogram for chamber 1 (the top one)
	TH1I *chamber2All = new TH1I("chamber2All","Chamber 2 After Edge Cuts",width/4,0,width);//histogram for chamber 2
	TH1I *chamber3All = new TH1I("chamber3All","Chamber 3 After Edge Cuts",width/4,0,width);//histogram for chamber 3
	TH1I *chamber4All = new TH1I("chamber4All","Chamber 4 After Edge Cuts",width/4,0,width);//histogram for chamber 4 (the bottom one)

	TCanvas *projCAll = new TCanvas("projCAll","",0,0,800,600);
	TCanvas *pc2All = new TCanvas("pc2All", "Stack of 4 Layer Runs", 0, 0, 800, 600);
	
	while (true) {	
		if (std::getline(stm, layers, ' ')) {

			//create the output root file
			outF0 = "sliceXCuts_" + layers + "layers.png";
			outfileName0 = outPrefix + outF0;
			std::cout << outF0 << " created\n";
			outF1 = "projYCuts_" + layers + "layers.png";
			outfileName1 = outPrefix + outF1;
			std::cout << outF1 << " created\n";
			
			//load the input root files
			TChain *chain = new TChain("fourChamTree");

			for (int i=0; ; i++) {
				runCount = std::to_string(i);
				inF = "run" + runCount + "_" + layers + "layers.root";
				infileName = inPrefix + inF;
				ifstream fin;
				fin.open(infileName.c_str());
				if (!fin.fail()) {
					fin.close();
					chain->Add(infileName.c_str());
					std::cout << "Got " << inF << std::endl;
				} else break;
			}


			int x=-10;	//x from file
			int y=-10;	//y from file
			int intensity=-10;	//pixle intensity from file
			int pNum=0;//the order in which the frame was processed

			//the 2d array which will store each frame of image data.
			int frame[480][640]={0};

			//variables
			int UNIXtime=0;
			float tdc[2]={-10,-10};

			//TTree *T = new TTree("T","TTree of muplus data");
			//add the 'branches' to the tree we will now read in
			chain->SetBranchAddress("pNum",&pNum);	//branch for the frame number
			chain->SetBranchAddress("frame",&frame);	//branch for frame data


			TH2I *frameHisto = new TH2I("frameHisto","Stacked Frames After Edge Cuts",width/4,0,width,height/4,0,height);	//histogram for the stacked images
			TH1I *chamber1 = new TH1I("chamber1","Chamber 1 After Edge Cuts",width/4,0,width);//histogram for chamber 1 (the top one)
			TH1I *chamber2 = new TH1I("chamber2","Chamber 2 After Edge Cuts",width/4,0,width);//histogram for chamber 2
			TH1I *chamber3 = new TH1I("chamber3","Chamber 3 After Edge Cuts",width/4,0,width);//histogram for chamber 3
			TH1I *chamber4 = new TH1I("chamber4","Chamber 4 After Edge Cuts",width/4,0,width);//histogram for chamber 4 (the bottom one)

			//loop over all data in chain
			Int_t nevent = chain->GetEntries();	//get the number of entries in the TChain
			for (Int_t i=0;i<nevent;i++) {
				chain->GetEntry(i);
				for(int x=0;x<width;x++){
					for(int y=0;y<height;y++){
						if(frame[x][y]>0){
							frameHisto->Fill(x,y,frame[x][y]);
							frameHistoAll->Fill(x,y,frame[x][y]);
							if(y>580 && y<610){
								chamber1->Fill(x,frame[x][y]);
								chamber1All->Fill(x,frame[x][y]);
							}
							else if(y>400 && y<440){
								chamber2->Fill(x,frame[x][y]);
								chamber2All->Fill(x,frame[x][y]);
							}
//.........这里部分代码省略.........

void makePUM0Table()
{
  const TString l1_input = "/mnt/hadoop/cms/store/user/luck/L1Emulator/minbiasHIanalyzer_withregions.root";
  TFile *lFile = TFile::Open(l1_input);
  TTree *l1Tree = (TTree*)lFile->Get("L1UpgradeAnalyzer/L1UpgradeTree");

  Int_t l1_event, l1_run;
  Int_t l1_num;
  Int_t region_hwPt[NREG], region_hwEta[NREG], region_hwPhi[NREG], region_tauVeto[NREG];

  l1Tree->SetBranchAddress("event",&l1_event);
  l1Tree->SetBranchAddress("run",&l1_run);
  l1Tree->SetBranchAddress("nRegions",&l1_num);
  l1Tree->SetBranchAddress("region_hwPt",region_hwPt);
  l1Tree->SetBranchAddress("region_hwEta",region_hwEta);
  l1Tree->SetBranchAddress("region_hwPhi",region_hwPhi);
  l1Tree->SetBranchAddress("region_tauVeto",region_tauVeto);

  const TString forest_input = "/mnt/hadoop/cms/store/user/luck/L1Emulator/minbiasForest_merged_v2/HiForest_PbPb_Data_minbias_fromSkim.root";
  TFile *fFile = TFile::Open(forest_input);
  TTree *fEvtTree = (TTree*)fFile->Get("hiEvtAnalyzer/HiTree");
  TTree *fSkimTree = (TTree*)fFile->Get("skimanalysis/HltTree");

  Int_t f_evt, f_run, f_lumi;
  Int_t hiBin;
  fEvtTree->SetBranchAddress("evt",&f_evt);
  fEvtTree->SetBranchAddress("run",&f_run);
  fEvtTree->SetBranchAddress("lumi",&f_lumi);
  fEvtTree->SetBranchAddress("hiBin",&hiBin);

  Int_t pcollisionEventSelection, pHBHENoiseFilter;
  fSkimTree->SetBranchAddress("pcollisionEventSelection",&pcollisionEventSelection);
  fSkimTree->SetBranchAddress("pHBHENoiseFilter",&pHBHENoiseFilter);

  TFile *outFile = new TFile("HI_PUM0_evtsel_out.root","RECREATE");

  std::map<Long64_t, Long64_t> kmap;

  // 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);

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

void Zlumi::Loop()
{
//   In a ROOT session, you can do:
//      Root > .L Zlumi.C
//      Root > Zlumi t
//      Root > t.GetEntry(12); // Fill t data members with entry number 12
//      Root > t.Show();       // Show values of entry 12
//      Root > t.Show(16);     // Read and show values of entry 16
//      Root > t.Loop();       // Loop on all entries
//

//     This is the loop skeleton where:
//    jentry is the global entry number in the chain
//    ientry is the entry number in the current Tree
//  Note that the argument to GetEntry must be:
//    jentry for TChain::GetEntry
//    ientry for TTree::GetEntry and TBranch::GetEntry
//
//       To read only selected branches, Insert statements like:
// METHOD1:
//    fChain->SetBranchStatus("*",0);  // disable all branches
//    fChain->SetBranchStatus("branchname",1);  // activate branchname
// METHOD2: replace line
//    fChain->GetEntry(jentry);       //read all branches
//by  b_branchname->GetEntry(ientry); //read only this branch
 
 gROOT->ForceStyle();
        tdrStyle();

  if (fChain == 0) return;

int minRun=0;
int maxRun=0;
int maxLS=0;
bool forminRun=true;


//TH2I * LumiSRun = new TH2I("LumiSRun", "LS vs Run", 3000, 0., 3000., 20000, 160000., 180000.);
//TH1I * test2 = new TH1I("test2","test2", 3000,0, 3000);


Long64_t nentries = fChain->GetEntriesFast();

Long64_t nbytes = 0, nb = 0;
for (Long64_t jentry=0; jentry<nentries;jentry++) {
	Long64_t ientry = LoadTree(jentry);
	if (ientry < 0) break;
	nb = fChain->GetEntry(jentry);   nbytes += nb;


	if(forminRun && (Run!=0)){minRun=Run; forminRun=false;}
	if((Run!=0) && (Run>maxRun)){maxRun=Run;}
	if((LS!=0) && (LS>maxLS)){maxLS=LS;}
	//printf("run %i ls %i \n",Run,LS);
}
cout << nentries << " nentries \n";

TH2I *LumiSRun = new TH2I("LumiSRun", "LS vs Run", maxLS, 0, maxLS, maxRun-minRun+2, minRun-1, maxRun+1);
TH1I *Runs = new TH1I("Runs","Run", maxRun-minRun+2, minRun-1, maxRun+1);
 Runs->Sumw2();
nbytes = 0; nb = 0;
for (Long64_t jentry=0; jentry<nentries;jentry++) {
	Long64_t ientry = LoadTree(jentry);
	if (ientry < 0) break;
	nb = fChain->GetEntry(jentry);   nbytes += nb;

	LumiSRun->Fill(LS,Run);
	//printf("run %i ls %i \n",Run,LS);
	Runs->Fill(Run);
	//test2->Fill(LS);
}
	printf("minRun %i maxRun %i \n",minRun,maxRun);
LumiSRun->Draw();
 for (int h=0;h<Runs->GetNbinsX();h++){
   Runs->SetBinError(h+1,sqrt(Runs->GetBinContent(h+1)) );
 }
Runs->Draw();
//test->Draw();
//test2->Draw();

TH1F *FileRuns = new TH1F("FileRuns","Run from Lumicalc", maxRun-minRun+2, minRun-1, maxRun+1);
TH1D *XsecDistro = new TH1D("XsecDistro","X sec distribution", 60, 0., 0.6);
//-------------
   const Int_t mpt = maxRun-minRun;
   int fileRun[mpt];
   double Lumi[mpt];

   int npt = 0;
   // read data file
   ifstream file;
   //file.open("./2011-run-lumi.txt");
   file.open("./LumiAeB-dav.txt");
   while (1) {

	   file >> fileRun[npt] >> Lumi[npt];
	   if ( ! file.good() ) break;
	   cout << "x = " << fileRun[npt] << " y = " << Lumi[npt] << endl;
	   
	   FileRuns->SetBinContent((fileRun[npt]-minRun+2),0.2);
	   npt++;
//.........这里部分代码省略.........