C++ TProfile::GetBinContent方法代码示例

int main(int argc,  char * argv[]){
  if(argc < 3){cout<<" Usage: executable input_file.root outfile.txt"<<endl; return -4;}
  int sm_ = 18;
  TFile f(argv[1]);
  TProfile* pr = 0;
  pr = (TProfile*) f.Get("rel_timing_prof_conv_blu");
  if( !pr ){cout<<" timing profile not found in the root file"<<endl; return -1;}
  
  float time[1701]; for(int i=0;i<1701;i++){time[i]=-100;}
  float half_mean  = 0;

  for(int ch=1;ch<1701;ch++){ 
    if ( ch<101 || (ch-1)%20 > 9){half_mean=137.4/25.;}// half1
    else {half_mean=137.1/25.;}
    time[ch] = pr->GetBinContent(ch)+half_mean;
    if( fabs(pr->GetBinContent(ch))>0.2 || pr->GetBinEntries(ch)<500 ){
      cout<<"cry: "<<ch<<" entries: "<<pr->GetBinEntries(ch)<<" rel timing: "<<pr->GetBinContent(ch)<<endl;
    }
  }
  f.Close();
  
  ofstream txt_channels;
  txt_channels.open(argv[2],ios::out);
  for(int i=1;i<1701;i++){
    txt_channels <<sm_<<"   "<<setw(4)<<i<<" \t "<<setw(8)<<setprecision(7)<<time[i]<< endl;
  }
  txt_channels.close();
  
  return 0;
}

//------------------------------------------------------------------------------
double meanY( char* hs ) // for profiles only
{
  TObject* obj = gDirectory->Get(hs);
  if( !obj->InheritsFrom( "TProfile" ) ) {
    cout << hs << " is not profile plot" << endl;
    return 0;
  }
  TProfile* p = (TProfile*)obj;
  double sumw = 0;
  double sumy = 0;
  for( int i = 1; i <= p->GetNbinsX(); ++i ) { // bin 0 is underflow
    double w = p->GetBinEntries(i);
    sumw += w;
    sumy += w * p->GetBinContent(i);
  }
  if( sumw > 0.1 )
    return sumy/sumw;
  else
    return 0;
}

void plotAnaMult4(char *infname="ana.root")
{    
   TFile *inf = new TFile(infname);
   
   TTree *tData = (TTree*) inf->Get("Data_tree");
   TTree *tMC   = (TTree*) inf->Get("MC_tree");
   TTree *tPP   = (TTree*) inf->Get("PP_tree");

   TCanvas *c1 = new TCanvas("c","",1200,700);
//   c->Divide(4,1);
   TCut recoCut = "leadingJetPt>120&&subleadingJetPt>50&&dphi>5*3.14159265358979/6.&&abs(leadingJetEta)<1.6&&abs(subleadingJetEta)<1.6";
   TCut genCut = "genleadingJetPt>120&&gensubleadingJetPt>50&&genDphi>5*3.14159265358979/6.&&abs(genleadingJetEta)<1.6&&abs(gensubleadingJetEta)<1.6";

   
   const int nPtBin=5;
   Float_t PtBins[nPtBin+1] = {0.001,0.5,1,1.5,2,3.2};
   
   Int_t centBin[5] = {200,100,60,20,0};
   makeMultiPanelCanvas(c1,4,2,0.0,0.0,0.2,0.2,0.02);
   
   for (int i=0;i<4;i++) {
      c1->cd(i+1);
      TH1D * empty=new TH1D("empty","",100,0,3.19);
//      empty->Fill(0.5,1000); 
      empty->SetMaximum(39.99); 
      empty->SetMinimum(0.001); 
      empty->SetNdivisions(105,"X");
      empty->GetXaxis()->SetTitleSize(28);
      empty->GetXaxis()->SetTitleFont(43); 
      empty->GetXaxis()->SetTitleOffset(1.8);
      empty->GetXaxis()->SetLabelSize(22);
      empty->GetXaxis()->SetLabelFont(43);
      empty->GetYaxis()->SetTitleSize(28);
      empty->GetYaxis()->SetTitleFont(43); 
      empty->GetYaxis()->SetTitleOffset(1.8);
      empty->GetYaxis()->SetLabelSize(22);
      empty->GetYaxis()->SetLabelFont(43);
      empty->GetXaxis()->CenterTitle();
      empty->GetYaxis()->CenterTitle();

      empty->SetXTitle("#Delta#eta_{1,2}");
      empty->SetYTitle("Multiplicity Difference");
   
      TProfile *pData = new TProfile(Form("pData%d",i),"",nPtBin,PtBins);
      TProfile *pMC = new TProfile(Form("pMC%d",i),"",nPtBin,PtBins);
      TProfile *pPP = new TProfile(Form("pPP%d",i),"",nPtBin,PtBins);
      TProfile *pGen = new TProfile(Form("pGen%d",i),"",nPtBin,PtBins);
      TCut centCut = Form("hiBin>=%d&&hiBin<%d",centBin[i+1],centBin[i]);
      tData->Draw(Form("-multDiff:abs(leadingJetEta-subleadingJetEta)>>pData%d",i),recoCut&&centCut);
      tPP->Draw(Form("-multDiff:abs(leadingJetEta-subleadingJetEta)>>pPP%d",i),recoCut);
      tMC->Draw(Form("-multDiff:abs(leadingJetEta-subleadingJetEta)>>pMC%d",i),recoCut&&centCut);  
      tMC->Draw(Form("-genMultDiff:abs(genleadingJetEta-gensubleadingJetEta)>>pGen%d",i),genCut&&centCut);  

      pMC->SetLineColor(2);
      pMC->SetMarkerColor(2);
      pPP->SetLineColor(4);
      pPP->SetMarkerColor(4);
      
//      pData->SetAxisRange(0,50,"Y");
//      pData->SetAxisRange(0,0.49,"X");
      empty->Draw();

      double diff=0;
      if (i==0) diff=0.1;
      drawText(Form("%d-%d %%",(int)(0.5*centBin[i+1]),(int)(0.5*centBin[i])),0.22+diff,0.65);
      if (i==0) drawText("PbPb #sqrt{s_{NN}}=2.76 TeV 150/#mub",0.22+diff,0.85);
      if (i==0) drawText("pp      #sqrt{s_{NN}}=2.76 TeV 5.3/pb",0.22+diff,0.75);
      if (i==3) drawText("CMS Preliminary",0.3+diff,0.85);
            
      Float_t sys[4]={1,1,2.5,3};
      

      for (int j=1;j<=pData->GetNbinsX();j++)
      {
         TBox *b = new TBox(pData->GetBinLowEdge(j),pData->GetBinContent(j)-sys[i],pData->GetBinLowEdge(j+1),pData->GetBinContent(j)+sys[i]);
	 //b->SetFillColor(kGray);
	 b->SetFillStyle(0);
	 b->SetLineColor(1);
	 b->Draw();
         TBox *b2 = new TBox(pPP->GetBinLowEdge(j),pPP->GetBinContent(j)-1,pPP->GetBinLowEdge(j+1),pPP->GetBinContent(j)+1);
	 //b2->SetFillColor(65);
	 b2->SetFillStyle(0);
	 b2->SetLineColor(4);
	 b2->Draw();
      }
      pData->Draw("same");
      pMC->Draw("same");
      pPP->Draw("same");
      pGen->SetMarkerColor(4);
      pGen->SetMarkerStyle(24);
      pData->Draw("same");
      
//      pGen->Draw(" same");

      c1->cd(5+i);
      TH1D *empty2 = (TH1D*)empty->Clone("empty2");
      empty2->SetYTitle("PbPb - pp");
      empty2->SetMinimum(-5);
      empty2->SetMaximum(+29.99);
      empty2->Draw();
//.........这里部分代码省略.........

void process()
{
   Int_t raw[512];   // buffer for input signal and bkg trees
   // buffers for output trees
   Int_t sig[512];
   Int_t cmsig[512];
   Int_t cm[16];

   // pedestal
   const char* fbkg_name = "Raw_Data_FZ320P_05_MSSD_2_250V_K237_Pedestal.dat-events.root";
   TFile* fbkg = TFile::Open(fbkg_name);
   if (!fbkg) cout<< "File not found: " << fbkg <<endl<<exitl;

   TTree* tree = (TTree*) fbkg->Get("etree");
   tree->SetBranchAddress("raw", &raw);

   TH2* h2d = (TH2*) fbkg->Get("h2d");
   new TCanvas;
   h2d->Draw();

   TProfile* profile = (TProfile*) fbkg->Get("profile"); 
   //new TCanvas;
   //profile->Draw();

   Int_t pedestal[512];
   for (int i=0; i<512; i++) {
      // pedestal[i] = profile->GetBinContent(i+1) - 0.5;
      // pedestal[i] = profile->GetBinContent(i+1) + 0.5;
      pedestal[i] = profile->GetBinContent(i+1);
   }

   // cout << "\nPedestals for every channel\n" << endl;
   // for (int i=0; i<512; i++) {
   //    cout << pedestal[i] << " ";
   //    if (i>0 && (i+1)%128==0)
   //       cout << endl;
   // }

   cout<< "processing bkg" <<endl;

   // output file with tree
   const char* obfname = "FZ320P_05_MSSD_2-bkg.root";
   TFile* obfile = TFile::Open(obfname, "recreate");

   TTree* btree = new TTree("btree", "btree");
   btree->Branch("sig", &sig, "sig[512]/I");
   btree->Branch("cmsig", &cmsig, "cmsig[512]/I");
   btree->Branch("cm", &cm, "cm[16]/I");
   btree->SetMarkerStyle(6);
   btree->SetMarkerColor(2);

   for (int jentry=0; jentry<tree->GetEntries(); ++jentry)
   {
      tree->GetEvent(jentry);
      // sig
      for (int i=0; i<512; ++i) {
         sig[i] = raw[i] - pedestal[i];
      }
      // calc common mode
      Int_t group32[32];
      Int_t index32[32];
      for (int igroup=0; igroup<16; ++igroup) {
         for (int istrip=0; istrip<32; ++istrip)   // istrip is number inside group of 32
         {
            group32[istrip] = sig[igroup*32 + istrip];
         }
         // sort array group32 in ascending order
         TMath::Sort(32, group32, index32, kFALSE);
         Int_t median = group32[index32[14]];
         cm[igroup] = median;
      }
      // subtract common mode
      for (int istrip=0; istrip<512; ++istrip) {
         Int_t igroup = istrip/32;
         cmsig[istrip] = sig[istrip] - cm[igroup];
      }
      // Fill sig, cmsig, cm
      btree->Fill();
   }
   obfile->Write();

   /////////////////////////////////////////////////
   //
   // signal tree
   //
   /////////////////////////////////////////////////

   cout<< "processing signal" <<endl;
   
   TChain* chain = new TChain("etree");
   chain->Add("Raw_Data_FZ320P_05_MSSD_250V_K237_Position_1.dat-events.root");
   chain->Add("Raw_Data_FZ320P_05_MSSD_250V_K237_Position_2.dat-events.root");
   chain->Add("Raw_Data_FZ320P_05_MSSD_2_250V_K237_Position_3.dat-events.root");
   chain->Add("Raw_Data_FZ320P_05_MSSD_2_250V_K237_Position_4.dat-events.root");
   chain->SetBranchAddress("raw", &raw);

   // output file with tree
   const char* osfname = "FZ320P_05_MSSD_2-signal.root";
   TFile* osfile = TFile::Open(osfname, "recreate");

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

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

  if(ppPbPb){
    for(int i=0;i<tB->GetEntries();i++){
      tB->GetEntry(i);
      if(fabs(jtetaB)<2 && binB>=cbinlo && binB<cbinhi && abs(refparton_flavorForBB)==5)
	hB->Fill(refptB,jtptB/refptB,weightB);
    }
    for(int i=0;i<tC->GetEntries();i++){
      tC->GetEntry(i);
      if(fabs(jtetaC)<2 && binC>=cbinlo && binC<cbinhi && abs(refparton_flavorForBC)==4)
	hC->Fill(refptC,jtptC/refptC,weightC);
    }
  }

  
 
  hL->SetMinimum(0.);
  
  hL->SetLineColor(kBlue);
  hB->SetLineColor(kRed);
  hC->SetLineColor(kGreen);

  hL->SetMarkerColor(kBlue);
  hB->SetMarkerColor(kRed);
  hC->SetMarkerColor(kGreen);
  
  //hL->SetMarkerStyle(4);
  //hB->SetMarkerStyle(4);
  //hC->SetMarkerStyle(4);
  
  hL->SetXTitle("genJet p_{T} (GeV/c)");
  hL->SetYTitle("<reco p_{T} / gen p_{T} >");

  hL->GetXaxis()->SetRangeUser(50.,199.999);
  hL->GetYaxis()->SetRangeUser(0.5,1.05);
  
  TCanvas *c1=new TCanvas("c1","c1",800,600);
  c1->SetGridx(1);
  c1->SetGridy(1);

  hL->Draw("e1");
  hB->Draw("e1,same");
  hC->Draw("e1,same");

  TLegend *leg=new TLegend(0.4,0.15,0.9,0.45);
  leg->SetBorderSize(0);
  leg->SetFillStyle(0);
  if(ppPbPb&&cbinlo==0&&cbinhi==40)leg->SetHeader("Pythia+Hydjet, 0-100%");
  leg->AddEntry(hL,"Inclusive jets","pl");
  leg->AddEntry(hC,"c-jets","pl");
  leg->AddEntry(hB,"b-jets","pl");
  leg->Draw();

  TCanvas *c2=new TCanvas("c2","c2",1);
  /*
  TH1F *hL2 = (TH1F*)hL->Clone("hL2");
  TH1F *hB2 = (TH1F*)hB->Clone("hB2");
  hL2->Add(hB2,-1);
  hL2->Draw();
  */

  TH1F  *hcorr = new TH1F("hcorr","hcorr",250,50,300);
  hcorr->Sumw2();

  for(int i=0;i<hL->GetNbinsX();i++){
    cout<<" b resp "<<hB->GetBinContent(i+1)<<endl;
    cout<<" l resp "<<hL->GetBinContent(i+1)<<endl;
    cout<<" l offset "<<1.-hL->GetBinContent(i+1)<<endl;
    cout<<" corrected b resp "<<hB->GetBinContent(i+1)+1.-hL->GetBinContent(i+1)<<endl;
    float jesOffset = 1.-hL->GetBinContent(i+1);

    hcorr->SetBinContent(i+1,hB->GetBinContent(i+1)+jesOffset);

    hcorr->SetBinError(i+1,sqrt(hB->GetBinError(i+1)*hB->GetBinError(i+1)+hL->GetBinError(i+1)*hL->GetBinError(i+1)));


  }

  hcorr->SetMinimum(0.5);
  hcorr->SetMaximum(1.1);
      
  hcorr->SetLineColor(kRed);
  hcorr->SetMarkerColor(kRed);
  hcorr->SetMarkerStyle(4);
  hcorr->Draw();

  TF1 *fCorr = new TF1("fCorr","[0]+[1]*log(x)+[2]*log(x)*log(x)",50,300);
  fCorr->SetLineWidth(1);
  fCorr->SetLineColor(kBlue);
  hcorr->Fit(fCorr);

  TFile *fout;
  if(ppPbPb) fout =new TFile(Form("bJEShistos/bJetScale_PbPb_Cent_fineBin_%d_%d.root",cbinlo,cbinhi),"recreate");
  else fout =new TFile("bJEShistos/bJetScale_PP_fineBin.root","recreate");
  hcorr->Write();
  fCorr->Write();
  fout->Close();

}

void read_tpvar(const char* file, int YS, 
      const char* inputtable_pt, const char* outputtable_pt, 
      const char* inputtable_rap, const char* outputtable_rap, 
      const char* inputtable_cent="nocentrality", const char* outputtable_cent="nocentrality")
{
   // YS decides on the binning.
   // YS=1 -> fine binning
   // YS=2 -> coarse binning
   // YS=0 -> both
   // YS=4 -> both, except centrality where only fine (for 1S in pbpb)
   TFile *f = new TFile("tpvar.root","RECREATE");
   TTree *tr = new TTree("tree","tree");
   tr->ReadFile(file,"name/C:binlow/F:binhigh/F:eff/F:stat/F:syst/F");

   unsigned int NPTNS = YS==1 ? NPT1S : NPT2S;
   unsigned int NRAPNS = YS==1 ? NRAP1S : NRAP2S;
   unsigned int NCENTNS = YS==1 ? NCENT1S : NCENT2S;

   double *ptbins_NS = YS==1 ? ptbins_1S : ptbins_2S;
   double *rapbins_NS = YS==1 ? rapbins_1S : rapbins_2S;
   int *centbins_NS = YS==1 ? centbins_1S : centbins_2S;

   bool docentrality=strcmp(inputtable_cent,"nocentrality");
   ifstream intable_pt(inputtable_pt);
   ifstream intable_rap(inputtable_rap);
   ifstream intable_cent; if (docentrality) intable_cent.open(inputtable_cent);
   ofstream outtable_pt(outputtable_pt);
   ofstream outtable_rap(outputtable_rap);
   ofstream outtable_cent; if (docentrality) outtable_cent.open(outputtable_cent);
   double binmin, binmax, eff0, efferr, effmean;

   TString var("pt_SF");
   tr->Draw("eff:1>>htemp(1,0,2)",Form("name==\"%s\"&&abs(binlow-%f)<.1&&abs(binhigh-%f)<.1",var.Data(),ptbins_NS[0],ptbins_NS[NPTNS]),"PROFs");
   // tr->Draw("eff:1>>htemp(1,0,2)",Form("name==\"%s\"&&abs(binlow-%f)<.1&&abs(binhigh-%f)<.1",var.Data(),ptbins_NS[0],40.),"PROFs"); // ugly fix
   TProfile *htemp = (TProfile*) gDirectory->Get("htemp");
   binmin=ptbins_NS[0];
   binmax=ptbins_NS[NPTNS]; //40.;//(ugly fix)
   eff0=firsteff(tr,var.Data(),binmin,binmax);
   efferr=htemp->GetBinError(1);
   effmean=htemp->GetBinContent(1);
   cout << setprecision(3)<<fixed << var << " " << binmin << " " << binmax << " " << eff0 << " " << efferr << " " << effmean << endl;
   if (dotables) print(intable_pt,outtable_pt,binmin,binmax,eff0,efferr);
   delete htemp;

   skip_lines(intable_pt,1);// \hline

   if (YS==0||YS==4)
   {
      NPTNS=NPT2S; ptbins_NS=ptbins_2S;
   }
   for (unsigned int i=0; i<NPTNS; i++)
   {
      // cout << setprecision(3)<<fixed << Form("name==\"%s\"&&binlow==%f&&binhigh==%f",var.Data(),ptbins_1S[i],ptbins_1S[i+1]) << endl;
      tr->Draw("eff:1>>htemp(1,0,2)",Form("name==\"%s\"&&abs(binlow-%f)<.1&&abs(binhigh-%f)<.1",var.Data(),ptbins_1S[i],ptbins_1S[i+1]),"PROFs");
      htemp = (TProfile*) gDirectory->Get("htemp");
      binmin=ptbins_NS[i];
      binmax=ptbins_NS[i+1];
      eff0=firsteff(tr,var.Data(),ptbins_NS[i],ptbins_NS[i+1]);
      efferr=htemp->GetBinError(1);
      effmean=htemp->GetBinContent(1);
      cout << setprecision(3)<<fixed << var << " " << binmin << " " << binmax << " " << eff0 << " " << efferr << " " << effmean << endl;
      if (dotables) print(intable_pt,outtable_pt,binmin,binmax,eff0,efferr);
      delete htemp;
   }

   if (YS==0||YS==4)
   {
      NPTNS=NPT1S; ptbins_NS=ptbins_1S;
      outtable_pt << "\\hline" << endl;
      skip_lines(intable_pt,1);// \hline

      for (unsigned int i=0; i<NPT1S; i++)
      {
         // cout << setprecision(3)<<fixed << Form("name==\"%s\"&&binlow==%f&&binhigh==%f",var.Data(),ptbins_1S[i],ptbins_1S[i+1]) << endl;
         tr->Draw("eff:1>>htemp(1,0,2)",Form("name==\"%s\"&&abs(binlow-%f)<.1&&abs(binhigh-%f)<.1",var.Data(),ptbins_1S[i],ptbins_1S[i+1]),"PROFs");
         htemp = (TProfile*) gDirectory->Get("htemp");
         binmin=ptbins_NS[i];
         binmax=ptbins_NS[i+1];
         eff0=firsteff(tr,var.Data(),ptbins_NS[i],ptbins_NS[i+1]);
         efferr=htemp->GetBinError(1);
         effmean=htemp->GetBinContent(1);
         cout << setprecision(3)<<fixed << var << " " << binmin << " " << binmax << " " << eff0 << " " << efferr << " " << effmean << endl;
         if (dotables) print(intable_pt,outtable_pt,binmin,binmax,eff0,efferr);
         delete htemp;
      }
   }

   var = TString("rapidity_SF");
   tr->Draw("eff:1>>htemp(1,0,2)",Form("name==\"%s\"&&abs(binlow-%f)<.1&&abs(binhigh-%f)<.1",var.Data(),rapbins_NS[0],rapbins_NS[NRAPNS]),"PROFs");
   htemp = (TProfile*) gDirectory->Get("htemp");
   binmin=rapbins_NS[0];
   binmax=rapbins_NS[NRAPNS];
   eff0=firsteff(tr,var.Data(),rapbins_NS[0],rapbins_NS[NRAPNS]);
   efferr=htemp->GetBinError(1);
   effmean=htemp->GetBinContent(1);
   cout << setprecision(3)<<fixed << var << " " << binmin << " " << binmax << " " << eff0 << " " << efferr << " " << effmean << endl;
   delete htemp;

   if (YS==0||YS==4)
   {
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  Float_t z;
  
  //T->Branch("nevento",&nevento,"nevento/I");
  //T->SetBranchAddress("ntracks",&ntracks);
  T->SetBranchAddress("ncluster",&ncluster);
  T->SetBranchAddress("tempo",tempo);
  T->SetBranchAddress("DeltaX",DeltaX);
  T->SetBranchAddress("DeltaZ",DeltaZ);
  T->SetBranchAddress("ChannelTOF",ChannelTOF);
  T->SetBranchAddress("impulso_trasv",&impulso_trasv);
  T->SetBranchAddress("exp_time_pi",exp_time_pi);
  T->SetBranchAddress("L",L);
  T->SetBranchAddress("TOT",TOT);
  T->SetBranchAddress("res",res);
  T->SetBranchAddress("charge",&charge);
  T->SetBranchAddress("phi",&phi);
  T->SetBranchAddress("eta",&eta);
  T->SetBranchAddress("secPhi",&secAngle);
  T->SetBranchAddress("cval",cval);
  T->SetBranchAddress("thetay",&thetay);
  T->SetBranchAddress("StartTime",&StartTime);
  T->SetBranchAddress("StartTimeRes",&StartTimeRes);
  
  
  Int_t nentries = (Int_t)T->GetEntries();
  
  for(Int_t i=0;i<nentries;i++) {
    T->GetEntry(i);
    
    for(Int_t ip=0;ip<ncluster;ip++){
      tempo[ip] -= StartTime;
      Int_t strip=ChannelTOF[0]/96;
      if(kCal){
	DeltaX[ip] -= hCalX->GetBinContent(strip+1);
	DeltaZ[ip] -= hCalZ->GetBinContent(strip+1);

      }
    }
    
    if(ncluster == 1){
      
      if(impulso_trasv>1.3){
	hx->Fill(int(ChannelTOF[0]/96),DeltaX[0]);
	hz->Fill(int(ChannelTOF[0]/96),DeltaZ[0]);

	if(!kCal){
	  hCalX->Fill(int(ChannelTOF[0]/96),DeltaX[0]);
	  hCalZ->Fill(int(ChannelTOF[0]/96),DeltaZ[0]);	  
	}

	h2resxz1->Fill(DeltaX[0],DeltaZ[0]);
	h2resxz1->GetXaxis()->SetTitle("Dx1 (cm)");
	h2resxz1->GetYaxis()->SetTitle("Dz1 (cm)");
	
	hresx1->Fill(DeltaX[0]);
	hresx1->GetXaxis()->SetTitle("Dx1 (cm)");
	
	hresz1->Fill(DeltaZ[0]);
	hresz1->GetXaxis()->SetTitle("Dz1 (cm)");
	
	hresdist1->Fill(sqrt(DeltaX[0]*DeltaX[0]+DeltaZ[0]*DeltaZ[0]));
	hresdist1->GetXaxis()->SetTitle("sqrt(Dx^2+Dz^2) (cm)");
      }

      if(impulso_trasv>0.8 && impulso_trasv<1.){ // serve per gli exp time
	//if( TMath::Abs(DeltaZ[0])<1.75){

void extractFlowVZERO(Int_t icentr,const char *type,Int_t arm,Float_t pTh,Bool_t isMC,Int_t addbin){
  LoadLib();

  char name[100];
  snprintf(name,100,"AnalysisResults%s.root",type);
  if(!fo) fo = new TFile(name);

  new TCanvas();

  Int_t cMin[] = {0,5,10,20,30,40,50,60,70};
  Int_t cMax[] = {5,10,20,30,40,50,60,70,80};

  if(kNUOcorr){ // Compute correction for NUO in TOF
    compareTPCTOF(icentr,0,arm,pTh,addbin);
//     compareTPCTOF(icentr,1,arm,pTh,addbin);
//     compareTPCTOF(icentr,2,arm,pTh,addbin);
//     compareTPCTOF(icentr,3,arm,pTh,addbin);
//     compareTPCTOF(icentr,4,arm,pTh,addbin);
  }

  TProfile *pAll;
  pAll=extractFlowVZEROsingle(icentr,0,arm,0,pTh,addbin,"all",0,1);
  pAll->SetMarkerStyle(24);
  TProfile *pPiTOF,*pPiTPC,*pPiTPC2;
  pPiTOF=extractFlowVZEROsingle(icentr,1,arm,0,pTh,addbin,"piTOF",1,1);
  pPiTPC=extractFlowVZEROsingle(icentr,1,arm,0,pTh,addbin,"piTPC",0,0);
  pPiTPC2=extractFlowVZEROsingle(icentr,1,arm,0,pTh,addbin,"piTPC2",2,2);
  pPiTPC->Add(pPiTPC2);

  TH1D *hPi = pPiTOF->ProjectionX("hPi");
  hPi->SetLineColor(4);
  hPi->SetMarkerColor(4);
  hPi->SetMarkerStyle(20);
  for(Int_t i=1;i <=hPi->GetNbinsX();i++){
    Float_t x = hPi->GetBinCenter(i);
    if(x < 0.2){
      hPi->SetBinContent(i,0);
      hPi->SetBinError(i,0);
    }
    else if(x < 0.5){
      hPi->SetBinContent(i,pPiTPC->GetBinContent(i));
      hPi->SetBinError(i,pPiTPC->GetBinError(i));      
    }
    else{
      if(kNUOcorr){
	hPi->SetBinContent(i,pPiTOF->GetBinContent(i) + hNUO[0]->GetBinContent(i));
	hPi->SetBinError(i,pPiTOF->GetBinError(i));      
      }
      else{
	hPi->SetBinContent(i,pPiTOF->GetBinContent(i));
	hPi->SetBinError(i,pPiTOF->GetBinError(i));      
      }
    }
  }

  TProfile *pElTOF,*pElTPC,*pElTPC2;
  pElTOF=extractFlowVZEROsingle(icentr,4,arm,0,pTh,addbin,"piTOF",1,1);
  pElTPC=extractFlowVZEROsingle(icentr,4,arm,0,pTh,addbin,"piTPC",0,0);
  pElTPC2=extractFlowVZEROsingle(icentr,4,arm,0,pTh,addbin,"piTPC2",2,2);
  pElTPC->Add(pElTPC2);

  TH1D *hEl = pElTOF->ProjectionX("hEl");
  hEl->SetLineColor(6);
  hEl->SetMarkerColor(6);
  hEl->SetMarkerStyle(20);
  for(Int_t i=1;i <=hEl->GetNbinsX();i++){
    Float_t x = hEl->GetBinCenter(i);
    if(x < 0.2){
      hEl->SetBinContent(i,0);
      hEl->SetBinError(i,0);
    }
    else if(x < 0.3){
      hEl->SetBinContent(i,pElTPC->GetBinContent(i));
      hEl->SetBinError(i,pElTPC->GetBinError(i));      
    }
    else{
      if(kNUOcorr){
	hEl->SetBinContent(i,pElTOF->GetBinContent(i) + hNUO[0]->GetBinContent(i));
	hEl->SetBinError(i,pElTOF->GetBinError(i));      
      }
      else{
	hEl->SetBinContent(i,pElTOF->GetBinContent(i));
	hEl->SetBinError(i,pElTOF->GetBinError(i));      
      }
    }
  }

  TProfile *pKTOF,*pKTPC,*pKTPC2;
  pKTOF=extractFlowVZEROsingle(icentr,2,arm,0,pTh,addbin,"kaTOF",1,1);
  pKTPC=extractFlowVZEROsingle(icentr,2,arm,0,pTh,addbin,"kaTPC",0,0);
  pKTPC2=extractFlowVZEROsingle(icentr,2,arm,0,pTh,addbin,"kaTPC2",2,2);
  pKTPC->Add(pKTPC2);

  TH1D *hK = pKTOF->ProjectionX("hKa");
  hK->SetLineColor(1);
  hK->SetMarkerColor(1);
  hK->SetMarkerStyle(21);
  for(Int_t i=1;i <=hK->GetNbinsX();i++){
    Float_t x = hK->GetBinCenter(i);
    if(x < 0.25){
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    hDen1D->SetName(hName); 
    
    if(opt.Contains("comov"))
      hDen1D->GetXaxis()->SetTitle("k_{p}#zeta");
    else
      hDen1D->GetXaxis()->SetTitle("k_{p}z");
  
    hDen1D->GetYaxis()->SetTitle("n_{b}/n_{0}");
  }

  // On-axis beam density vs \zeta vs time! _________________________________
  TH2F *hDen1DvsTime = NULL; 
  if(hDen1D) {
    char hName[24];
    sprintf(hName,"hDen1DvsTime");
    TH2F *hDen1DvsTimeOld = (TH2F*) ifile->Get(hName);

    Int_t nBins   = 1;
    Float_t edge0 = Time-0.5;
    Float_t edge1 = Time+0.5;
    if(hDen1DvsTimeOld!=NULL) {
      nBins = hDen1DvsTimeOld->GetNbinsX()+1;
      Float_t binwidth =  (Time - hDen1DvsTimeOld->GetXaxis()->GetBinCenter(1))/(nBins-1);
      edge0 = hDen1DvsTimeOld->GetXaxis()->GetBinCenter(1) - binwidth/2.;
      edge1 = Time + binwidth/2.;
    }
    hDen1DvsTime = new TH2F("temp","",nBins,edge0,edge1,
		       	hDen1D->GetNbinsX(),
			hDen1D->GetBinLowEdge(1),
			hDen1D->GetBinLowEdge(hDen1D->GetNbinsX()+1));
    
    for(Int_t ix=1;ix<hDen1DvsTime->GetNbinsX();ix++) {
      for(Int_t iy=1;iy<hDen1DvsTime->GetNbinsY();iy++) {
	hDen1DvsTime->SetBinContent(ix,iy,hDen1DvsTimeOld->GetBinContent(ix,iy));
      }
    }  
    delete hDen1DvsTimeOld;
  
    // Fill last bin with the newest values.
    for(Int_t iy=1;iy<=hDen1D->GetNbinsX();iy++) {
      hDen1DvsTime->SetBinContent(nBins,iy,hDen1D->GetBinContent(iy));
    }   

    hDen1DvsTime->GetZaxis()->SetTitle("n_{b}/n_{0}");
    hDen1DvsTime->GetYaxis()->SetTitle("k_{p}#zeta");
    hDen1DvsTime->GetXaxis()->SetTitle("k_{p}z");
    hDen1DvsTime->GetZaxis()->CenterTitle();
    hDen1DvsTime->GetYaxis()->CenterTitle();
    hDen1DvsTime->GetXaxis()->CenterTitle();
    hDen1DvsTime->SetName(hName);

    // Change the range of z axis 
    Float_t Denmax = hDen1DvsTime->GetMaximum();
    hDen1DvsTime->GetZaxis()->SetRangeUser(0,Denmax); 
    hDen1DvsTime->Write(hName,TObject::kOverwrite);

  }

  // RMS (vs z) of the beam's charge distribution: 
  TProfile *hDen2Dprof = NULL;
  TH1F *hRms = NULL;
  Double_t axisPos = x2Mid;
  if(hDen2D) {
    TString pname = hDen2D->GetName();
    pname += "_pfx";
    

//.........这里部分代码省略.........
	      // (response difference between pT1 and vecpT2~10%)
	      // => 0.25*0.10 = 2.5%
	      // For data/MC, estimate uncertainty as half of this
	      // => 1.25%
	      ga->SetPoint(n, 2.5, 0.);
	      if (itype!=2) ga->SetPointError(n, 0, 0.025);
	      if (itype==2) ga->SetPointError(n, 0, 0.0125);
	    } // MPF
	  } // use priors

	  if (ga->GetN()>2) {

	    f1->SetRange(minalpha, 3.);
	    ga->Fit(f1,"QRN");

	    if (f1->GetNDF()>=0) {
	      f1->DrawClone("SAME");
	      f1->SetRange(0,0.4);
	      f1->SetLineStyle(kDashed);
	      f1->DrawClone("SAME");

	      // Store results
	      TGraphErrors *gk = gkmap[ct][cm][cs];
	      if (!gk) {
		gk = new TGraphErrors(0);
		gk->SetMarkerStyle(ga->GetMarkerStyle());
		gk->SetMarkerColor(ga->GetMarkerColor());
		gk->SetLineColor(ga->GetLineColor());
		gkmap[ct][cm][cs] = gk;
	      }
	      int n = gk->GetN();
	      TProfile *ppt = (isample==0 ? ppt2 : ppt1);
	      if (isample==3) { ppt = ppt4; } // pas-v6
	      double pt = ppt->GetBinContent(ppt->FindBin(ipt));
	      gk->SetPoint(n, pt, f1->GetParameter(1));
	      gk->SetPointError(n, 0, f1->GetParError(1));
	    } // f1->GetNDF()>=0
	  } // ga->GetN()>2
	} // for itpt
	
      } // for isample
      
      c2->SaveAs(Form("pdf/softrad_3x3_%s_%s_vsalpha.pdf",ct,cm));
      
    }
  }


  cout << "Drawing plots of kFSR vs pT" << endl;

  // 2x6 plots
  for (int itype = 0; itype != ntypes; ++itype) {
    for (int  imethod = 0; imethod != nmethods; ++imethod) {

      const char *ct = types[itype];
      const char *cm = methods[imethod];

      TMultiGraph *mgk = new TMultiGraph();

      int ipad = ntypes*imethod + itype + 1; assert(ipad<=6);
      c3->cd(ipad);
      gPad->SetLogx();
      h3->SetMaximum(imethod==0 ? 0.05 : (itype!=2 ? 0.1 : 0.25));
      h3->SetMinimum(imethod==0 ? -0.05 : (itype!=2 ? -0.4 : -0.25));
      h3->SetYTitle(Form("k_{FSR} = dR/d#alpha (%s)",ct));
      h3->DrawClone("AXIS");

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

          fmc->Close();
        } // dotrigeff


        // Scale data to account for time dependence
        bool dotimedep = ((string(obj2->GetName())=="hpt") || isjk || isjet);
        TH1D *htimedep = (TH1D*)outdir->FindObject("htimedep");
        TH1D *htimefit = (TH1D*)outdir->FindObject("htimefit");
        TH1D *hpt_notimedep = 0, *hpt_withtimedep = 0;
        double ktime = 1.;

        if (!htimedep) {

          TH1D *h = (TH1D*)dir->Get("hpt");
          TH1D *hsel = (TH1D*)dir->Get("hselpt");
          TH1D *hpre = (TH1D*)dir->Get("hpt_pre");
          //TH1D *hlumi0 = (TH1D*)dir->Get("../jt450/hlumi");

          // Fix luminosity for unprescaled trigger
          //string strg = dir->GetName();
          //double lum0 = triglumi["jt450"];
          //for (int i = 1; i != hlumi0->GetNbinsX()+1; ++i) {
          //hlumi0->SetBinContent(i, lum0);
          //}

          outdir->cd();
          if (h) hpt_notimedep = (TH1D*)h->Clone("hpt_notimedep");
          if (hpre && h) htimedep = (TH1D*)hpre->Clone("htimedep");
          if (hpre && h) htimedep->Divide(hpre,h);//,1,1,"B");

          // Figure out trigger luminosities
          double lumi = 0;
          if (hlumi) lumi = hlumi->GetBinContent(1);
          double lumi0 = 0;
          if (hlumi0) lumi0 = hlumi0->GetBinContent(1);
          if (htimedep && lumi && lumi0) {
            htimedep->Scale(lumi / lumi0);
          }

          // Find proper pT range and fit
          double minpt = 0.;
          double maxpt = 6500.;
          if (hsel) {
            for (int i = 1; i != hsel->GetNbinsX()+1; ++i) {
              if (hsel->GetBinContent(i)!=0 &&
                  hsel->GetBinLowEdge(i)>=_jp_xmin57) {
                if (minpt<20) minpt = hsel->GetBinLowEdge(i);
                maxpt = hsel->GetBinLowEdge(i+1);
              }
            }
          }
          TF1 *ftmp = new TF1("ftmp","[0]",minpt,maxpt);
          ftmp->SetParameter(0,1);
          if (htimedep && htimedep->Integral()>0) htimedep->Fit(ftmp,"QRN");

          if (htimedep && ftmp->GetParameter(0)>0)
            ktime = 1./ftmp->GetParameter(0);

          if (htimedep) {
            outdir->cd();
            htimefit = (TH1D*)hsel->Clone("htimefit");
            hpt_withtimedep = (TH1D*)h->Clone("hpt_withtimedep");

            for (int i = 1; i != htimefit->GetNbinsX()+1; ++i) {