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

void DileptonAnalyzer::storeSignalEfficiencies( TFile & outputFile ) {
  // Store in a TProfile
  TProfile *sigEffOne = new TProfile("SignalEfficOneDecay","Signal efficiencies when one exotic decays in this channel",10,0.5,10.5,-9.9e9,9.9e9);
  TAxis *sigEffOneAxis = sigEffOne->GetXaxis();

  TProfile *sigEffTwo = new TProfile("SignalEfficTwoDecay","Signal efficiencies when two exotic decays in this channel",10,0.5,10.5,-9.9e9,9.9e9);
  TAxis *sigEffTwoAxis = sigEffTwo->GetXaxis();

  // Fill set of 3 bins in each TProfile for each exotic
  for ( int iExotic = 0; iExotic < nSignalParticles_; iExotic++ ) {

    TString axisTitle;

    sigEffOne->Fill( iExotic*nSignalParticles_+1, numEvents_oneSensitiveDecay_[iExotic] );
    axisTitle.Form("Number of events coming into macro for exotic %d",signalPdgId_[iExotic]);
    sigEffOneAxis->SetBinLabel( iExotic*nSignalParticles_+1, axisTitle.Data());

    sigEffOne->Fill( iExotic*nSignalParticles_+2, numExoticsRECO_oneSensitiveDecay_[iExotic] );
    axisTitle.Form("Number of candidates RECO'd for exotic %d",signalPdgId_[iExotic]);
    sigEffOneAxis->SetBinLabel( iExotic*nSignalParticles_+2, axisTitle.Data());

    sigEffOne->Fill( iExotic*nSignalParticles_+3, numExoticsCorrectRECO_oneSensitiveDecay_[iExotic] );
    axisTitle.Form("Number of exotics correctly RECO'd for exotic %d",signalPdgId_[iExotic]);
    sigEffOneAxis->SetBinLabel( iExotic*nSignalParticles_+3, axisTitle.Data());

    sigEffTwo->Fill( iExotic*nSignalParticles_+1, numEvents_twoSensitiveDecay_[iExotic] );
    axisTitle.Form("Number of events coming into macro for exotic %d",signalPdgId_[iExotic]);
    sigEffTwoAxis->SetBinLabel( iExotic*nSignalParticles_+1, axisTitle.Data());

    sigEffTwo->Fill( iExotic*nSignalParticles_+2, numExoticsRECO_twoSensitiveDecay_[iExotic] );
    axisTitle.Form("Number of candidates RECO'd for exotic %d",signalPdgId_[iExotic]);
    sigEffTwoAxis->SetBinLabel( iExotic*nSignalParticles_+2, axisTitle.Data());

    sigEffTwo->Fill( iExotic*nSignalParticles_+3, numExoticsCorrectRECO_twoSensitiveDecay_[iExotic] );
    axisTitle.Form("Number of exotics correctly RECO'd for exotic %d",signalPdgId_[iExotic]);
    sigEffTwoAxis->SetBinLabel( iExotic*nSignalParticles_+3, axisTitle.Data());
  }
  outputFile.cd();
  sigEffOne->Write();
  sigEffTwo->Write();
}

void profile()
{
  TCanvas* c1 = new TCanvas("canvas","nhitac",1200,600);
  c1->Divide(2,1);
  c1->cd(1);
  TFile *file1 = new TFile("test_nhitac_modified.root");
  //TH1F  *histo1 = new TH1F("nhitac1","nhitac1",60,0,300);
  TProfile *profile1 = new TProfile("profile","nhitac profile",100,0,100,0,300);

  TFile *file2 = new TFile("test_nhitac_unmodified.root");
  //TH1F  *histo2 = new TH1F("nhitac","nhitac",60,0,300);
  TProfile *profile2 = new TProfile("profile2","nhitac profile2",100,0,100,0,300);
  TProfile *ratio = new TProfile("ratio","nhitac ratio", 100,0,100,0,2);
  //histo2->SetMarkerStyle(31);
  TLegend *l1 = new TLegend(0.20, 0.60, 0.3, 0.7);
  l1->SetBorderSize(0);

  TTree* nhitac_modified = (TTree*)file1->Get("testnhitac");
  TTree* nhitac_unmodified = (TTree*)file2->Get("testnhitac");
  int testnhitac1, testnhitac2;
  nhitac_modified->SetBranchAddress("nhitac",&testnhitac1);
  nhitac_unmodified->SetBranchAddress("nhitac",&testnhitac2);
  for(int i= 0; i < nhitac_modified->GetEntries(); i++)
  {
    nhitac_modified->GetEntry(i);
    nhitac_unmodified->GetEntry(i);
    profile1->Fill(i,testnhitac1);
    profile2->Fill(i,testnhitac2);
    double temp = testnhitac2;
    ratio->Fill(i,testnhitac1/temp);
    //histo1->Fill(testnhitac1);
    //histo2->Fill(testnhitac2);
  }
  //histo1->SetMarkerColor(kRed);
  //histo1->SetLineColor(kRed);
  profile1->SetMarkerColor(kRed);
  profile1->SetMarkerStyle(5);
  profile2->SetMarkerColor(kBlue);
  l1->AddEntry(profile1,"modified od bad channel","P");
  l1->AddEntry(profile2,"old od bad channel","l");
  l1->SetTextSize(0.04);
  profile1->Draw("e1p");
  profile1->GetXaxis()->SetTitle("event #");
  profile1->GetYaxis()->SetTitle("nhitac");
  profile2->Draw("same");
  l1->Draw();
  c1->cd(2);
  ratio->GetXaxis()->SetTitle("event #");
  ratio->GetYaxis()->SetTitle("new/old");
  //ratio->SetMarkerStyle(5);
  //ratio->SetMarkerColor(kRed);
  ratio->Draw();
}

// Report cut flow findings to screen
void DileptonAnalyzer::reportCutFlow( TFile & outputFile ) {

  // Store in a TProfile
  TProfile *cutFlow = new TProfile("cutFlow","Cut flow for final set of analysis cuts",20,0.5,20.5,-9.9e9,9.9e9);
  TAxis *cutFlowAxis = cutFlow->GetXaxis();

  for ( unsigned int iCut = 0; iCut < cutNamesInOrder_.size(); iCut++ ) {
    cutFlow->Fill( iCut+1, cutFlowMap_[cutNamesInOrder_[iCut]]);
    cutFlowAxis->SetBinLabel( iCut+1, cutNamesInOrder_[iCut]);
  }

  outputFile.cd();
  cutFlow->Write();
}

TH1D * smartGausProfileXSQRTN (TH2F * strip, double width){
 TProfile * stripProfile = strip->ProfileX () ;
 
 // (from FitSlices of TH2.h)
 
 double xmin = stripProfile->GetXaxis ()->GetXmin () ;
 double xmax = stripProfile->GetXaxis ()->GetXmax () ;
 int profileBins = stripProfile->GetNbinsX () ;
 
 std::string name = strip->GetName () ;
 name += "_smartGaus_X" ; 
 TH1D * prof = new TH1D(name.c_str (),strip->GetTitle (),profileBins,xmin,xmax) ;
 
 int cut = 0 ; // minimum number of entries per fitted bin
 int nbins = strip->GetXaxis ()->GetNbins () ;
 int binmin = 1 ;
 int ngroup = 1 ; // bins per step
 int binmax = nbins ;
 
 // loop over the strip bins
 for (int bin=binmin ; bin<=binmax ; bin += ngroup) 
 {
  TH1D *hpy = strip->ProjectionY ("_temp",bin,bin+ngroup-1,"e") ;
  if (hpy == 0) continue ;
  int nentries = Int_t (hpy->GetEntries ()) ;
  if (nentries == 0 || nentries < cut) {delete hpy ; continue ;} 
  
  Int_t biny = bin + ngroup/2 ;
  
  TF1 * gaussian = new TF1 ("gaussian","gaus", hpy->GetMean () - width * hpy->GetRMS (), hpy->GetMean () + width * hpy->GetRMS ()) ; 
  gaussian->SetParameter (1,hpy->GetMean ()) ;
  gaussian->SetParameter (2,hpy->GetRMS ()) ;
  hpy->Fit ("gaussian","RQL") ;           
  
  //       hpy->GetXaxis ()->SetRangeUser ( hpy->GetMean () - width * hpy->GetRMS (), hpy->GetMean () + width * hpy->GetRMS ()) ;         
  prof->Fill (strip->GetXaxis ()->GetBinCenter (biny), gaussian->GetParameter (1)) ;       
  prof->SetBinError (biny,gaussian->GetParameter (2) / sqrt(hpy->GetEntries())) ;
//   prof->SetBinError (biny,gaussian->GetParError (1)) ;
  
  delete gaussian ;
  delete hpy ;
 } // loop over the bins
 
 delete stripProfile ;
 return prof ;
}

//------------------------------------------------------------//
// Get TProfile
//------------------------------------------------------------//
TProfile* getProfile(TFile* file, TString pname, TString xtitle,
		     TString ytitle, int color, int marker)
{
  
  TProfile* prof = (TProfile*) (file->Get(pname.Data())->Clone(Form("%s_%i",pname.Data(),color)));
  prof->GetXaxis()->SetTitle(xtitle.Data());
  prof->GetYaxis()->SetTitle(ytitle.Data());
  prof->SetMarkerStyle(marker);
  prof->SetMarkerColor(color);
  //prof->SetMarkerSize(0.5);
  prof->SetLineColor(color);
  prof->SetTitle("");
  prof->SetStats(0);
  prof->GetYaxis()->SetTitleOffset(1.5);
  
  return prof;

}

TProfile* plotIsoPerformance( TFile* ftt, 
			      const char* signal,     // histogram name
			      const char* background, // histogram name
			      const char* name,       // unique name
			      bool reverse = false,   // normally signal near zero bin, reverse means signal is around max bin
			      double bkg_eff_min = 0,
			      double bkg_eff_max = 1,
			      double sig_eff_min = 0,
			      double sig_eff_max = 0
			      )
{ 
  TH1F* S = dynamic_cast<TH1F*>(ftt->Get(signal));
  if ( ! S ) {
    std::cout << "Error: histogram not found " << signal << std::endl;
    return 0;
  }
  TH1F* B = dynamic_cast<TH1F*>(ftt->Get(background));
  if ( ! B ) {
    std::cout << "Error: histogram not found " << background << std::endl;
    return 0;
  }
  char buf[1024];
  sprintf(buf,"c_%s",name);
  // TCanvas* c = new TCanvas(buf,buf,500,500);
  sprintf(buf,"p_%s",name);
  TProfile* p = new TProfile(buf,buf,50,bkg_eff_min,bkg_eff_max,sig_eff_min,sig_eff_max);
  p->SetLineColor(kBlue);
  p->SetLineWidth(2);
  p->SetMarkerStyle(20);
  p->SetMarkerSize(1);
  p->GetXaxis()->SetTitle("Background Efficiency");
  p->GetYaxis()->SetTitle("Signal Efficiency");
  p->SetStats(kFALSE);
  for( int i=0;i<=S->GetNbinsX()+1; ++i )
    if ( reverse ) 
      p->Fill(B->Integral(i,B->GetNbinsX()+1)/B->Integral(0,B->GetNbinsX()+1),
	      S->Integral(i,S->GetNbinsX()+1)/S->Integral(0,S->GetNbinsX()+1));
    else
      p->Fill(B->Integral(0,i)/B->Integral(0,B->GetNbinsX()+1),
	      S->Integral(0,i)/S->Integral(0,S->GetNbinsX()+1));
  // p->Draw();
  return p;
}

TCanvas *Plot2D_profileX_my(TChain *data, TString branchname, TString binning,TString selection,TString opt,TString xLabel, TString yLabel){
  //type == 0: data only
  //type == 1: MC only
  //type == 2: data/MC
   
  TCanvas *c = new TCanvas("c","");
  data->Draw(branchname+">>data_hist"+binning,selection,opt);
  TH2F *d = (TH2F *) gROOT->FindObject("data_hist");

  TCanvas *c1 = new TCanvas("c1","");
  TProfile *prof = d->ProfileX("prof",1,-1,"s");
  prof->SetMarkerStyle(20);
  prof->SetMarkerSize(1);
  prof->Draw();
  prof->GetYaxis()->SetTitle(yLabel);
  prof->GetXaxis()->SetTitle(xLabel);

  return c1;
}

void Fake100PeVShower(int opt)
{

  // opt == 0 implies save to root file
  // opt == 1 implies plot and save to canvas
  TFile* file = NULL;
  TCanvas* c  = NULL;

  if(opt == 0)
    file = new TFile("fake100PeVShower.root","recreate");
  if(opt == 1)
    c = makeCanvas("c");

  TString gaus = "1.5e7*TMath::Exp(-(pow(x-10,2)/9))";
  TString bump = "5e6*TMath::Exp(-(pow(x-17,2)/20))";
  TF1* f = new TF1("f",(gaus+"+"+bump).Data(),0,30);

  int nbins = 3000;
  float xmin  = 0;
  float xmax  = 30;
  TProfile* prof = new TProfile("prof","",nbins,xmin,xmax);
  prof->SetStats(0);
  prof->SetTitle("");
  prof->GetYaxis()->SetTitle("Charge excess / 1e7");
  prof->GetXaxis()->SetTitle("z [m]");

  float step = xmax / nbins;
  for(int i =0; i<nbins; ++i){
    float x = step * i;
    if(opt == 1 ) prof->Fill(x, f->Eval(x)/1e7);
    else          prof->Fill(x, f->Eval(x));
  }
    
  if(opt == 1){
    prof->Draw();
    c->SaveAs("JaimeCheck/FakeProfile100PeV.png");
  }
  if(opt == 0){
    file->Write();
    file->Close();
  }
}

void plotQApid(Int_t ic,Float_t pt,Int_t addbin){
  char name[100];
  char stringa[200];
  LoadLib();

  snprintf(name,100,"AnalysisResults.root");
  if(!fo) fo = new TFile(name);
  snprintf(name,100,"contVZEROv2");
  TList *cont = (TList *) fo->Get(name);
  AliFlowVZEROResults *pidqa = cont->FindObject("qaPID");
  Float_t xval[2] = {2.,pt+0.00001};
  Float_t xval2[2] = {2.+addbin,pt+0.00001};

  TProfile *proTPCpi = pidqa->GetV2(0,xval,xval2);
  TProfile *proTOFpi = pidqa->GetV2(1,xval,xval2);
  TProfile *proTPCka = pidqa->GetV2(2,xval,xval2);
  TProfile *proTOFka = pidqa->GetV2(3,xval,xval2);
  TProfile *proTPCpr = pidqa->GetV2(4,xval,xval2);
  TProfile *proTOFpr = pidqa->GetV2(5,xval,xval2);

  proTPCpi->SetName("hPiTPC");
  proTOFpi->SetName("hPiTOF");
  proTPCka->SetName("hKaTPC");
  proTOFka->SetName("hKaTOF");
  proTPCpr->SetName("hPrTPC");
  proTOFpr->SetName("hPrTOF");
  
  proTPCpi->GetXaxis()->SetTitle("dE/dx - dE/dx_{calc}^{#pi} (a.u)");
  proTOFpi->GetXaxis()->SetTitle("t_{tof} - t_{calc}^{#pi} (ps)");
  proTPCka->GetXaxis()->SetTitle("dE/dx - dE/dx_{calc}^{K} (a.u)");
  proTOFka->GetXaxis()->SetTitle("t_{tof} - t_{calc}^{K} (ps)");
  proTPCpr->GetXaxis()->SetTitle("dE/dx - dE/dx_{calc}^{p} (a.u)");
  proTOFpr->GetXaxis()->SetTitle("t_{tof} - t_{calc}^{p} (ps)");

  TCanvas *c = new TCanvas("cVproj","cVproj");
  c->Divide(2,3);
  c->cd(1);
  proTPCpi->Draw();
  c->cd(2);
  proTOFpi->Draw();
  c->cd(3);
  proTPCka->Draw();
  c->cd(4);
  proTOFka->Draw();
  c->cd(5);
  proTPCpr->Draw();
  c->cd(6);
  proTOFpr->Draw();
}

void Dcurvature(const char *chargechoice = "plus", double ptmin=20){

  TString sign=chargechoice;

  gStyle->SetPalette(1);
  gStyle->SetOptStat("e");

//  double ptmin=20;
  double ptmax=200;
  double etamax=2.1; 
  int ptbins=(ptmax-ptmin)/10;
  double ptbinwidth=(ptmax-ptmin)/ptbins;
  int etabins=21;
  int phibins=21;


  TProfile3D *khistptetaphi = new TProfile3D("DeltaCurv(pt,eta,phi)","DeltaCurv(pt,eta,phi) "+sign,ptbins,ptmin,ptmax,etabins,-etamax,etamax,phibins,-3.14,3.14);
  khistptetaphi->GetXaxis()->SetTitle("Pt");
  khistptetaphi->GetYaxis()->SetTitle("Eta");
  khistptetaphi->GetZaxis()->SetTitle("Phi");

  TProfile2D *khistpteta = new TProfile2D("DeltaCurv(pt,eta)","DeltaCurv(pt,eta) "+sign,ptbins,ptmin,ptmax,etabins,-etamax,etamax);
  khistpteta->GetXaxis()->SetTitle("Pt");
  khistpteta->GetYaxis()->SetTitle("Eta");
  TProfile2D *khistptphi = new TProfile2D("DeltaCurv(pt,phi)","DeltaCurv(pt,phi) "+sign,ptbins,ptmin,ptmax,phibins,-3.14,3.14);
  khistptphi->GetXaxis()->SetTitle("Pt");
  khistptphi->GetYaxis()->SetTitle("Phi");

  TProfile *khistpt = new TProfile("DeltaCurv(pt)","DeltaCurv(pt) "+sign,ptbins,ptmin,ptmax);
  khistpt->GetXaxis()->SetTitle("Pt");
  khistpt->SetAxisRange(-0.001,0.001,"Y");
  TProfile *khisteta = new TProfile("DeltaCurv(eta)","DeltaCurv(eta) "+sign,etabins,-etamax,etamax);
  khistpt->GetXaxis()->SetTitle("Eta");
  TProfile *khistphi = new TProfile("DeltaCurv(phi)","DeltaCurv(phi) "+sign,phibins,-3.14,3.14);
  khistpt->GetXaxis()->SetTitle("Phi");

  TObjArray *khistptbins= new TObjArray();
  for (int i=0; i<ptbins; i++) {
    TString name="DeltaCurv(eta,phi), pt bin ";
    name+=int(ptmin+i*(ptmax-ptmin)/ptbins);
    name+=TString(", charge ")+sign;
    TProfile2D *ist = new TProfile2D(name.Data(),name.Data(),phibins,-3.14,3.14,etabins,-etamax,etamax);
    ist->SetAxisRange(-0.002,0.002,"Z");
    khistptbins->Add(ist);
  }

  TFile *f = new TFile("RecoRoutines_Z-selection_ZJets_TuneZ2_7TeV_alpgen_tauola.rew8.corr1.root","read");

  TTree *tree;

  f->GetObject("SingleMuPtScale/"+sign+"muonstree",tree);

  Double_t MCPt;
  Double_t MCEta;
  Double_t MCPhi;
  Double_t RecoPt;
  Double_t RecoEta;
  Double_t RecoPhi;
  Double_t EvWeight;

  tree->SetBranchAddress("RecoPt",&RecoPt);
  tree->SetBranchAddress("RecoEta",&RecoEta);
  tree->SetBranchAddress("RecoPhi",&RecoPhi);
  tree->SetBranchAddress("MCPt",&MCPt);
  tree->SetBranchAddress("MCEta",&MCEta);
  tree->SetBranchAddress("MCPhi",&MCPhi);
  tree->SetBranchAddress("EvWeight",&EvWeight);

  long nentries = tree->GetEntriesFast();

  for (int i=0; i<nentries; i++){
    tree->GetEntry(i);
    if (RecoPt<ptmin || RecoPt>ptmax || RecoEta<-etamax || RecoEta>etamax) continue;
    double quantity=(MCPt-RecoPt)/MCPt/RecoPt;
    khistptetaphi->Fill(RecoPt,RecoEta,RecoPhi,quantity,EvWeight);
    khistpteta->Fill(RecoPt,RecoEta,quantity,EvWeight);
    khistptphi->Fill(RecoPt,RecoPhi,quantity,EvWeight);
    ((TProfile2D*)(khistptbins->At(int((RecoPt-ptmin)/ptbinwidth))))->Fill(RecoPhi,RecoEta,quantity,EvWeight);
    khistpt->Fill(RecoPt,quantity,EvWeight);
    khisteta->Fill(RecoEta,quantity,EvWeight);
    khistphi->Fill(RecoPhi,quantity,EvWeight);
  }



  TCanvas *c1 = new TCanvas();
  khistptetaphi->Draw("BOX");

  TCanvas *c2 = new TCanvas();
  c2->Divide(2,1);
  c2->cd(1);
  khistpteta->Draw("BOX");
  c2->cd(2);
  khistptphi->Draw("BOX");

  TCanvas *c2b = new TCanvas();
  c2b->Divide(3,1);
  c2b->cd(1);
  khistpt->Draw();
  c2b->cd(2);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
            if (mctrack->GetPdgId() == 2212  )  
            {
                v.SetPtEtaPhiM (track->GetPt(0), track->GetEta(0), track->GetPhi(0), 0.9386);
                hEffP->Fill ( p, ret[0] > 0.9 );
                hEffPSigma->Fill ( p,  sigmas[0] < 3&& sigmas[1] > 2 && sigmas[2] > 2 );
                hEffPtYP -> Fill( v.Rapidity() - 1.52, v.Pt(),  ret[0] > 0.9 );
            }
            
            if (mctrack->GetPdgId() == 321  )  
            {
                v.SetPtEtaPhiM (track->GetPt(0), track->GetEta(0), track->GetPhi(0), 0.5);
                hEffK->Fill ( p, ret[1] > 0.9 );
                hEffKSigma->Fill ( p,  sigmas[1] < 3&& sigmas[2] > 2 && sigmas[0] > 2 ); 
                hEffPtYK -> Fill(  v.Rapidity() - 1.52, v.Pt(),  ret[1] > 0.9 );
            }
            
            if (mctrack->GetPdgId() == 211  )  
            {
                v.SetPtEtaPhiM (track->GetPt(0), track->GetEta(0), track->GetPhi(0), 0.14);
                hEffPi->Fill ( p, ret[2] > 0.9 );
                hEffPiSigma->Fill ( p, sigmas[2] < 3&& sigmas[0] > 2 && sigmas[1] > 2 );  
                hEffPtYPi -> Fill( v.Rapidity() - 1.52, v.Pt(),  ret[2] > 0.9  );
            }
        }
    }

    hEffP       -> SetMarkerStyle(21);      hEffP       -> SetMarkerColor(kRed);     hEffP       -> SetLineColor(kRed); 
    hEffPi      -> SetMarkerStyle(21);      hEffPi      -> SetMarkerColor(kRed);     hEffPi      -> SetLineColor(kRed); 
    hEffK       -> SetMarkerStyle(21);      hEffK       -> SetMarkerColor(kRed);     hEffK       -> SetLineColor(kRed); 
    hEffPSigma  -> SetMarkerStyle(22);      hEffPSigma  -> SetMarkerColor(kBlue);    hEffPSigma  -> SetLineColor(kBlue);
    hEffPiSigma -> SetMarkerStyle(22);      hEffPiSigma -> SetMarkerColor(kBlue);    hEffPiSigma -> SetLineColor(kBlue);
    hEffKSigma  -> SetMarkerStyle(22);      hEffKSigma  -> SetMarkerColor(kBlue);    hEffKSigma  -> SetLineColor(kBlue);

    hEffP->GetXaxis()->SetTitle( "p, GeV/c" );
    hEffP->GetYaxis()->SetTitle(  "Efficiency"  );    
    hEffP->GetYaxis()->SetRangeUser(0.0, 1.);

    hEffK->GetXaxis()->SetTitle( "p, GeV/c" );
    hEffK->GetYaxis()->SetTitle(  "Efficiency"  );    
    hEffK->GetYaxis()->SetRangeUser(0.0, 1.);

    hEffPi->GetXaxis()->SetTitle( "p, GeV/c" );
    hEffPi->GetYaxis()->SetTitle(  "Efficiency"  );    
    hEffPi->GetYaxis()->SetRangeUser(0.0, 1.);
    
    hEffPtYP->GetYaxis()->SetTitle( "p_{T}, GeV/c" );
    hEffPtYP->GetXaxis()->SetTitle(  "Rapidity"  );
    hEffPtYK->GetYaxis()->SetTitle( "p_{T}, GeV/c" );
    hEffPtYK->GetXaxis()->SetTitle(  "Rapidity"  );
    hEffPtYPi->GetYaxis()->SetTitle( "p_{T}, GeV/c" );
    hEffPtYPi->GetXaxis()->SetTitle(  "Rapidity"  );
    
    TCanvas *c1 = new TCanvas ("c1", "c1", 1400, 700);
    c1->Divide(3,1);

    c1->cd(1);
    hEffP->Draw();
    hEffPSigma->Draw("same");
    c1->cd(2);
    hEffK->Draw();
    hEffKSigma->Draw("same");
    c1->cd(3);
    hEffPi->Draw();
    hEffPiSigma->Draw("same");

//     c1->cd(4);

//.........这里部分代码省略.........
      ptJetResWithEta->Fill(genJetTau1eta, (genJetTau1pt-jetTau1pt)/genJetTau1pt);
      ptJetResCorrWithEta->Fill(genJetTau1eta, (genJetTau1pt-corrJetTau1pt)/genJetTau1pt);

    }
    else if ( (genTau1pt != 999) && (jetTau1pt == 999) ) {
      tauTagEffWithPT->Fill(genTau1pt, 0);
      tauTagEffWithEta->Fill(genTau1eta, 0);
    }
    if ( (genTau2pt != 999) && (jetTau2pt != 999) ) {
      tauTagEffWithPT->Fill(genTau2pt, 1);
      //ptScaleWithPT->Fill(genTau2pt, jetTau2pt/genTau2pt);
      //ptJetScaleWithPT->Fill(genJetTau2pt, jetTau2pt/genJetTau2pt);
      //ptResWithPT->Fill(genTau2pt, (genTau2pt-jetTau2pt)/genTau2pt);
      ptJetResWithPT->Fill(genJetTau2pt, (genJetTau2pt-jetTau2pt)/genJetTau2pt);
      ptJetResCorrWithPT->Fill(genJetTau2pt, (genJetTau2pt-corrJetTau2pt)/genJetTau2pt);

      tauTagEffWithEta->Fill(genTau2eta, 1);
      //ptScaleWithEta->Fill(genTau2eta, jetTau2pt/genTau2pt);
      //ptJetScaleWithEta->Fill(genJetTau2eta, jetTau2pt/genJetTau2pt);
      //ptResWithEta->Fill(genTau2eta, (genTau2pt-jetTau2pt)/genTau2pt);
      ptJetResWithEta->Fill(genJetTau2eta, (genJetTau2pt-jetTau2pt)/genJetTau2pt);
      ptJetResCorrWithEta->Fill(genJetTau2eta, (genJetTau2pt-corrJetTau2pt)/genJetTau2pt);

    }
    else if ( (genTau2pt != 999) && (jetTau2pt == 999) ) {
      tauTagEffWithPT->Fill(genTau2pt, 0);
      tauTagEffWithEta->Fill(genTau2eta, 0);
    }

  }

  TCanvas *c1 = MakeCanvas("c1", "Tau Tagging Efficiency", 800, 600);
  tauTagEffWithPT->SetMarkerStyle(1);
  tauTagEffWithPT->GetXaxis()->SetTitle("P_{T} of Generator Tau");
  tauTagEffWithPT->GetYaxis()->SetTitle("Efficiency");
  tauTagEffWithPT->Draw();
  
  c1->SaveAs("tauTagEffWithPT.png");

  TCanvas *c2 = MakeCanvas("c2", "Tau Tagging Efficiency", 800, 600);
  tauTagEffWithEta->SetMarkerStyle(1);
  tauTagEffWithEta->GetXaxis()->SetTitle("Eta of Generator Tau");
  tauTagEffWithEta->GetYaxis()->SetTitle("Efficiency");
  tauTagEffWithEta->Draw();

  c2->SaveAs("tauTagEffWithEta.png");
  /*  
  TCanvas *c3 = MakeCanvas("c3", "Tau P_{T} Scale (reco Tau Jet P_{T}/gen Tau P_{T})", 800, 600);
  ptScaleWithPT->SetMarkerStyle(1);
  ptScaleWithPT->GetXaxis()->SetTitle("P_{T} of Generator Tau");
  ptScaleWithPT->GetYaxis()->SetTitle("P_{T} Scale");
  ptScaleWithPT->Draw();

  c3->SaveAs("ptScaleWithPT.png");

  TCanvas *c4 = MakeCanvas("c4", "Tau P_{T} Scale (reco Tau Jet P_{T}/gen Tau P_{T})", 800, 600);
  ptScaleWithEta->SetMarkerStyle(1);
  ptScaleWithEta->GetXaxis()->SetTitle("Eta of Generator Tau");
  ptScaleWithEta->GetYaxis()->SetTitle("P_{T} Scale");
  ptScaleWithEta->Draw();

  c4->SaveAs("ptScaleWithEta.png");

  TCanvas *c5 = MakeCanvas("c5", "Tau P_{T} Scale (reco Tau Jet P_{T}/gen Tau Jet P_{T})", 800, 600);
  ptJetScaleWithPT->SetMarkerStyle(1);
  ptJetScaleWithPT->GetXaxis()->SetTitle("P_{T} of Generator Tau Jet");

void
Analyze(const TString mode="CLOSED",
        UShort_t       maxH=6,
        Bool_t         /*doLoops*/=false,
        Int_t ifile=0   )
{
#ifdef __CINT__
  gROOT->LoadMacro("correlations/Types.hh++");
  gROOT->LoadMacro("correlations/Result.hh++");
  gROOT->LoadMacro("correlations/QVector.hh++");
  gROOT->LoadMacro("correlations/recursive/FromQVector.hh++");
  gROOT->LoadMacro("correlations/recurrence/FromQVector.hh++");
  gROOT->LoadMacro("correlations/closed/FromQVector.hh++");
  gROOT->LoadMacro("correlations/test/ReadData.hh++");
#endif
  // --- Setup of harmonics, etc -------------------------------------
  gRandom->SetSeed(54321);
  UShort_t emode = 0;
  if      (mode.EqualTo("closed",     TString::kIgnoreCase)) emode = 0;
  else if (mode.EqualTo("recurrence", TString::kIgnoreCase)) emode = 1;
  else if (mode.EqualTo("recursive",  TString::kIgnoreCase)) emode = 2;
  else
    Warning("Analyze", "Mode %s unknown, assuming CLOSED", mode.Data());

  correlations::QVector        q[nbin];
  correlations::FromQVector*   c[nbin];
  correlations::HarmonicVector h(maxH);
  for (UShort_t i = 0; i < maxH; i++) {
    // Generate random harmonicx
   // h[i] = -6 + gRandom->Integer(12);
   h[0] = 2;
   h[1] = -2;
   h[2] = 2;
   h[3] = -2;
   h[4] = 2;
   h[5] = -2;
   h[6] = -2;
   h[7] = 2;
   // Printf("h_%d:\t%d", i, h[i]);
  }

  // Resize the Q-vector to fit the harmonics
    for(int ibin=0;ibin<nbin;ibin++){
    q[ibin] = correlations::QVector(0,0,false);
    q[ibin].resize(h);
  switch (emode) {
  case 0: c[ibin] = new correlations::closed::FromQVector(q[ibin]); break;
  case 1: c[ibin] = new correlations::recurrence::FromQVector(q[ibin]); break;
  case 2: c[ibin] = new correlations::recursive::FromQVector(q[ibin]); break;
    }
  }
  //Printf("Correlator: %s", c->name());


  // --- Some histograms ---------------------------------------------
  TH1* sumreals[nbin];
  TH1* sumimags[nbin];
  TH1* weights[nbin];
  TVectorD tottrk;
  TVectorD Nevent;
  tottrk.ResizeTo(nbin);
  tottrk.Zero();
  Nevent.ResizeTo(nbin);
  Nevent.Zero();

  //TH1*      reals  = new TH1D("reals", "Re(C{n})", maxH-2+1, 2+.5, maxH+1+.5);
  //TH1*      imags  = static_cast<TH1*>(reals->Clone("imags"));
    for(int ibin=0;ibin<nbin;ibin++){
  sumreals[ibin]  = new TH1D(Form("sumreals_%d",ibin), "Re(C{n})", maxH-2+1, 2+.5, maxH+1+.5);
  sumimags[ibin]  = static_cast<TH1*>(sumreals[ibin]->Clone(Form("sumimags_%d",ibin)));
  weights[ibin]  = static_cast<TH1*>(sumreals[ibin]->Clone(Form("weights_%d",ibin)));
    }
  TProfile* timing = new TProfile("timing", "Timing", maxH-2+1, 2+.5,maxH+1+.5);
  TH1*      hs     = new TH1I("harmonics", "Harmonics", maxH, 1.5, maxH+1.5);
  /*reals->SetFillColor(kGreen+1);
  reals->SetFillStyle(3001);
  reals->SetStats(0);
  imags->SetTitle("Im(C{n})");
  imags->SetFillColor(kBlue+1);
  imags->SetFillStyle(3001);
  imags->SetStats(0);
  timing->SetFillColor(kRed+1);
  timing->SetFillStyle(3001);
  timing->SetStats(0);
  hs->SetFillColor(kMagenta+1);
  hs->SetFillStyle(3001);
  hs->SetStats(0);*/
  for (UShort_t i = 0; i < maxH-1; i++) {
    TString label = TString::Format("C{%d}", i+2);
//    reals->GetXaxis()->SetBinLabel(i+1, label);
//    imags->GetXaxis()->SetBinLabel(i+1, label);
    timing->GetXaxis()->SetBinLabel(i+1, label);
    hs->GetXaxis()->SetBinLabel(i+1,Form("h_{%d}", i+1));
    hs->SetBinContent(i+1, h[i]);
  }
  hs->GetXaxis()->SetBinLabel(maxH,Form("h_{%d}", maxH));
  hs->SetBinContent(maxH, h[maxH-1]);

  TStopwatch timer;
 
//.........这里部分代码省略.........

void balanceMetVsAj(TString infname = "dj_HCPR-J50U-hiGoodMergedTracks_OfficialSelv2_Final0_120_50.root",TCut myCut = "cent<30", char *title = "",bool drawLegend = false)
{
   TFile *inf = new TFile(infname);
   TTree *t = (TTree*)inf->Get("ntjt");

   const int nBin = 3;
   double bins[nBin+1] = {0.5,4,8,1000};  
   double colors[nBin] = {kOrange-8,kBlue-3,kRed};

   // Selection cut
   TCut evtCut = "nljet>100&&abs(nljetacorr)<2&&aljet>50&&abs(aljetacorr)<2&&jdphi>2./3*TMath::Pi()&&!maskEvt";

   cout << "Sel evt: " << t->GetEntries(evtCut&&myCut) << endl;

   TProfile *p[nBin];

   for (int i=0;i<nBin;i++)
   {
      p[i] = new TProfile(Form("p%d",i),"",5,0.0001,0.49999);
      t->Project(Form("p%d",i),Form("-1*metxMerged%d:Aj",i), "1"*(evtCut&&myCut));
      p[i]->SetLineColor(colors[i]);     
      p[i]->SetMarkerColor(colors[i]);     
   }

   TProfile *pall = new TProfile("pall","",5,0.0001,0.49999);
   pall->SetXTitle("A_{J}");
   pall->SetYTitle("<#slash{p}_{T}^{#parallel}> (GeV/c)");
   pall->GetXaxis()->CenterTitle();
   pall->GetYaxis()->CenterTitle();
   pall->SetNdivisions(505);
   t->Project("pall","-1*metx:Aj", "1"*(evtCut&&myCut));
   pall->SetAxisRange(-50,50,"Y");
   pall->SetMarkerSize(1);
   pall->Draw("E");
   for (int i=0;i<nBin;++i) {
      p[i]->SetLineWidth(3);
      p[i]->Draw("e hist same");
   }
   pall->Draw("E same");

   // Legend
   TLegend *leg = new TLegend(0.10,0.68,0.70,0.96);
   leg->SetFillStyle(0);
   leg->SetBorderSize(0);
   leg->SetTextFont(63);
   leg->SetTextSize(16);
   leg->AddEntry(pall,Form("> %.1f GeV/c",bins[0]),"pl");
   for (int i=0;i<nBin;++i) {
      if (i!=nBin-1){
         leg->AddEntry(p[i],Form("%.1f - %.1f GeV/c",bins[i],bins[i+1]),"l");
      } else {
         leg->AddEntry(p[i],Form("> %.1f GeV/c",bins[i]),"l");
      }
   }

   if (drawLegend) leg->Draw();

   TLine * l0 = new TLine(0,0,0.5,0);
   l0->SetLineStyle(2);
   l0->Draw();

   TText *titleText = new TText(0.3,30,title);
   titleText->Draw();
}

void draw_clouds_profiles()
{
    //        gROOT->ProcessLine(".L ../utils.C");
    //        gROOT->ProcessLine(".L AliLRCFit.cxx");
    TFile *f[8];

    //    f[0] = new TFile( "output_classesByV0M_LHC10h.root" );
    //    f[0] = new TFile( "output_classesByV0M_LHC10h_c10_5_1.root" );
        f[0] = new TFile( "output_classesByV0M_LHC10h_c10_5_25_1_05.root" );
//    f[0] = new TFile( "output_classesByV0M_LHC11h_FemtoPlus_c10_5_CUT_OUTLIERS.root" );

//    f[0] = new TFile( "output_classesByV0M_LHC15o_fieldMM_c10_5_CUT_OUTLIERS.root" );
//    f[0] = new TFile( "output_classesByV0M_LHC15o_fieldPP_c10_5_CUT_OUTLIERS.root" );

//    const int nCW = 2; //nCentrWidths
//    const double cWidths[nCW] = { 10, 5 }; //width of the centrality bins
//    const double cStep[nCW] = { 5, 2.5 }; //centrality bins step
//    const int nCentrBins[nCW] = { 17, 35 }; //n centrality bins

    //    const int nCW = 3; //nCentrWidths
    //    const double cWidths[nCW] = { 10, 5, 1.0 }; //width of the centrality bins
    //    const double cStep[nCW] = { 5, 2.5, 1.0 }; //centrality bins step
    //    const int nCentrBins[nCW] = { 17, 35, 90 }; //n centrality bins

    //    const int nCW = 4; //nCentrWidths
    //    const double cWidths[nCW] = { 10, 5, 1.0, 0.5 }; //width of the centrality bins
    //    const double cStep[nCW] = { 5, 2.5, 1.0, 1.0 }; //centrality bins step
    //    const int nCentrBins[nCW] = { 17, 35, 90, 90 }; //n centrality bins

        const int nCW = 5; //nCentrWidths
        const double cWidths[nCW] = { 10, 5, 2.5, 1.0, 0.5 }; //width of the centrality bins
        const double cStep[nCW] = { 5, 2.5,  2.5, 1.0, 1.0 }; //centrality bins step
        const int nCentrBins[nCW] = { 17, 35, 36,  90, 90 }; //n centrality bins


    TH2D *hist2D;//[200][3];
    TProfile *profile;//[200][3];
    int cW = 2;
    int etaW = 1;
    int phiW = 0;

    const int kCorrType = 1; //0-NN, 1-PtPt, 2-PtN

    TCanvas *canv_tmp_for_fit = new TCanvas("canv_tmp_for_fit","canv_tmp_for_fit",50,50,300,300 );
    TCanvas *canv_2D_clouds = new TCanvas("canv_2D_clouds","canv_2D_clouds",150,250,1400,600 );
    tuneCanvas(canv_2D_clouds);
    canv_2D_clouds->Divide(2,1);

    gStyle->SetOptStat( kFALSE );

    TGraph *grFromFit2D = new TGraph;

    bool firstDraw = true;
    //    for ( int cBin = 0; cBin < nCentrBins[cW]; cBin++ )
    for ( int cBin = nCentrBins[cW]-1; cBin >= 0; cBin-- )
    {
                if (cBin%2!=0)
                    continue;

        //        cout << "cBin=" << cBin << endl;

        float cBinMin = cStep[cW] * cBin;
        float cBinMax = cWidths[cW] + cStep[cW] * cBin;

        // ##### pad 1 - clouds
        tunePad( canv_2D_clouds->cd(1) );
        if ( kCorrType == 0 )
        {
            hist2D = (TH2D*)f[0]->Get( Form("hist2D_NN_c%.1f-%.1f_etaW_%d_phiW_%d", cBinMin, cBinMax, etaW, phiW) );
            hist2D->SetTitle( "");
            hist2D->GetXaxis()->SetTitle( "N_{ch} Forward");
            hist2D->GetYaxis()->SetTitle( "N_{ch} Backward");
            hist2D->GetXaxis()->SetRangeUser(0,650);
            hist2D->GetYaxis()->SetRangeUser(0,650);
        }
        else if ( kCorrType == 1 )
        {
            hist2D = (TH2D*)f[0]->Get( Form("hist2D_PtPt_c%.1f-%.1f_etaW_%d_phiW_%d", cBinMin, cBinMax, etaW, phiW) );
            hist2D->SetTitle( "");
            hist2D->GetXaxis()->SetTitle( "#LTp_{T}#GT Forward");
            hist2D->GetYaxis()->SetTitle( "#LTp_{T}#GT Backward");
        }     hist2D->SetMarkerColor(kOrange-9+cBin);
        tuneHist2D_onPad(hist2D);
        hist2D->GetXaxis()->CenterTitle();
        hist2D->GetYaxis()->CenterTitle();

        //        removeBinsWithFewEntries(hist2D);

        hist2D->DrawCopy( firstDraw ? "" : "same" );

        // ##### pad 2 - profiles
        tunePad( canv_2D_clouds->cd(2) );
        profile = hist2D->ProfileX();    //(TProfile*)f[0]->Get( Form("hist2D_c%.1f-%.1f_etaW_%d_phiW_%d_pfx", cBinMin, cBinMax, etaW, phiW) );
        if ( kCorrType == 0 )
        {
            profile->SetTitle( "");
            profile->GetYaxis()->SetTitle( "#LTN_{ch}#GT Backward");
            profile->GetXaxis()->SetRangeUser(0,650);
            profile->GetYaxis()->SetRangeUser(0,650);
        }
//.........这里部分代码省略.........