C++ TH1F::DrawCopy方法代码示例

//_____________________________________________________________________________
void ProofEventProc::Terminate()
{
   // The Terminate() function is the last function to be called during
   // a query. It always runs on the client, it can be used to present
   // the results graphically or save the results to file.

   // Check ranges
   CheckRanges();

   if (gROOT->IsBatch()) return;
   
   TCanvas* canvas = new TCanvas("event","event",800,10,700,780);
   canvas->Divide(2,2);
   TPad *pad1 = (TPad *) canvas->GetPad(1);
   TPad *pad2 = (TPad *) canvas->GetPad(2);
   TPad *pad3 = (TPad *) canvas->GetPad(3);
   TPad *pad4 = (TPad *) canvas->GetPad(4);

   // The number of tracks
   pad1->cd();
   pad1->SetLogy();
   TH1F *hi = dynamic_cast<TH1F*>(fOutput->FindObject("pz_dist"));
   if (hi) {
      hi->SetFillColor(30);
      hi->SetLineColor(9);
      hi->SetLineWidth(2);
      hi->DrawCopy();
   } else { Warning("Terminate", "no pz dist found"); }

   // The Pt distribution
   pad2->cd();
   pad2->SetLogy();
   TH1F *hf = dynamic_cast<TH1F*>(fOutput->FindObject("pt_dist"));
   if (hf) {
      hf->SetFillColor(30);
      hf->SetLineColor(9);
      hf->SetLineWidth(2);
      hf->DrawCopy();
   } else { Warning("Terminate", "no pt dist found"); }

   // The Px,Py distribution, color surface
   TH2F *h2f = dynamic_cast<TH2F*>(fOutput->FindObject("px_py"));
   if (h2f) {
      // Color surface
      pad3->cd();
      h2f->DrawCopy("SURF1 ");
      // Lego
      pad4->cd();
      h2f->DrawCopy("CONT2COL");
   } else {
      Warning("Terminate", "no px py found");
   }

   // Final update
   canvas->cd();
   canvas->Update();
}

void testIlcGlauberQuenching() {

  IlcFastGlauber g;
  g.Init(2);
  g.SetCentralityClass(0.00,0.10); // centrality (fraction of geometrical cross section)

  IlcQuenchingWeights afq;
  afq.InitMult();
  afq.SetQTransport(1.);
  afq.SetECMethod(0);
  afq.SetLengthMax(8.);
  afq.SampleEnergyLoss();

  TCanvas *c = new TCanvas("cELD","Energy Loss Distribution",0,0,800,500);
  c->Divide(2,1);

  for(Int_t itype=1;itype<=2;itype++){
    c->cd(itype);
    gPad->SetLogy();
    Char_t name[100];
    Char_t hname[100];
    if(itype==1){
      sprintf(name,"Energy Loss Distribution - Quarks;E_{loss} (GeV);#"); 
      sprintf(hname,"hQuarks"); 
    } else {
      sprintf(name,"Energy Loss Distribution - Gluons;E_{loss} (GeV);#"); 
      sprintf(hname,"hGluons"); 
    }

    TH1F *h = new TH1F(hname,name,100,0,200);

    for(Int_t i=0;i<10000;i++){
      if(i % 100 == 0) cout << "." << flush;
      Double_t ell;
      g.GetLength(ell);
      Double_t loss=afq.GetELossRandom(itype,ell);
      h->Fill(loss);
    }

    h->SetStats(kTRUE);
    if(itype==1){
      h->SetLineColor(1);
      h->DrawCopy();
    }
    else {
      h->SetLineColor(2);
      h->DrawCopy();
    }
    delete h;
  }

  c->Update();
}

// second example using ell distribution
void testFastIlcGlauberQuenching(Char_t *fname) {

  TFile f(fname);
  TH1F *hEll=(TH1F*)f.Get("hEll");
  if(!hEll) return;

  IlcQuenchingWeights afq;
  afq.InitMult();
  afq.SetQTransport(1.);
  afq.SetECMethod(0);
  afq.SetLengthMax(5.);
  afq.SampleEnergyLoss();

  TCanvas *c = new TCanvas("cELD","Energy Loss Distribution",0,0,800,500);
  c->Divide(2,1);

  for(Int_t itype=1;itype<=2;itype++){
    c->cd(itype);
    gPad->SetLogy();
    Char_t name[200];
    if(itype==1)
      sprintf(name,"Energy Loss Distribution - Quarks;E_{loss} (GeV);#"); 
    else 
      sprintf(name,"Energy Loss Distribution - Gluons;E_{loss} (GeV);#"); 
    Char_t hname[200];
    if(itype==1)
      sprintf(hname,"hQuarks"); 
    else 
      sprintf(hname,"hGluons"); 

    TH1F *h = afq.ComputeELossHisto(itype,1.,hEll);

    h->SetStats(kTRUE);
    if(itype==1){
      h->SetLineColor(1);
      h->DrawCopy();
    }
    else {
      h->SetLineColor(2);
      h->DrawCopy();
    }
    delete h;
  }

  c->Update();
}

//================================================
void toyMC_Cent()
{
  gStyle->SetOptStat(0);
  gStyle->SetOptFit(1);
  gStyle->SetStatY(0.9);                
  gStyle->SetStatX(0.9);  
  gStyle->SetStatW(0.2);                
  gStyle->SetStatH(0.2);

  myRandom = new TRandom3();
  TDatime *clock = new TDatime();
  myRandom->SetSeed(clock->GetTime());


  // Get true refMult vs Ncoll distribution
  TH2F *hRefMultVsNcoll = 0x0;

  // Get different efficiency curve vs. refMult
  const int nEff = 2;
  const int nMult = 2;
  TF1 *funcEff[nEff];
  funcEff[0] = new TF1("funcEff_0", "pol1",0,1000);
  funcEff[0]->SetParameters(0.9, -2e-4);
  funcEff[1] = new TF1("funcEff_1", "pol0",0,1000);
  funcEff[1]->SetParameter(0,0.6);
  TCanvas *c = new TCanvas("cEff","cEff",800,600);
  TH1F *hplot = new TH1F("hplot",";Mult;",100,0,1000);
  hplot->SetYTitle("Efficiency");
  hplot->DrawCopy();
  for(int i=0; i<2; i++)
    {
      funcEff[i]->SetLineColor(i+1);
      funcEff[i]->SetLineWidth(1.5);
      funcEff[i]->Draw("sames");
    }
  TPaveText *title = GetTitleText("TPC tracking efficiency");
  title->Draw();
  TLegend *leg = new TLegend(0.4,0.2,0.7,0.4);
  leg->SetBorderSize(0);
  leg->SetFillColor(0);
  leg->SetTextSize(0.04);
  leg->AddEntry(funcEff[0], "Efficiency 1", "L");
  leg->AddEntry(funcEff[1], "Efficiency 2", "L");
  leg->Draw();

  // Get ngTrack distribution to mimic efficiency distribution
  
  TH1F *hMcMult[nMult];
  
}

void plot1D( TH2F* hul ){

  cout << "PLOT1D" << endl;
  TFile* f = TFile::Open("stop_xsec.root");
  TH1F* hxsec = (TH1F*) f->Get("h_stop_xsec");

  TH1F* hulproj = (TH1F*) hul->ProjectionX("hulproj",1,1);

  TCanvas* cproj = new TCanvas();
  gPad->SetLogy();

  hulproj->Draw("hist");

  hxsec->SetLineColor(2);
  hxsec->DrawCopy("samehistl");

  hxsec->SetLineColor(4);
  hxsec->Scale(1.15);
  hxsec->DrawCopy("samehistl");

  hxsec->Scale(0.85/1.15);
  hxsec->DrawCopy("samehistl");

}

void display(){
  
  gROOT->ProcessLine(".x lhcbStyle.C");
  gStyle->SetPalette(1);
	
  string filename = "../analysed/SET9THL914.root";
	TFile* file = new TFile(filename.c_str(),"OPEN");

	int nBinsX=32; int nBinsY=32;
	int nPixelsX=64; int nPixelsY=64;
	
  std::cout<<"Making global intercept tot map"<<std::endl;
  // Global intercept map (not per pixel)
	TH2F* hToTMap = new TH2F("hToTMap","hToTMap",nBinsX,0,nPixelsX,nBinsY,0,nPixelsY);
	for(int x=0;x<nBinsX;x++){
		for(int y=0;y<nBinsY;y++){
			int id = x + y*nBinsX;
			// Get the ToT distribution for this location
			string histoname = "/tpanal/ClicpixAnalysis/hClusterTOTAssociated1pixMap"+makestring(id);
			TH1F* totMapHisto = 0;
			totMapHisto = (TH1F*)gDirectory->Get(histoname.c_str());
			if(totMapHisto == 0) continue;
			if(totMapHisto->GetEntries() < 20) continue;
	
			// Fit the tot distribution and enter the value in the ToT map
			totMapHisto->Fit("gaus","0q");
			TF1* fit = totMapHisto->GetFunction("gaus");
			hToTMap->Fill(x*floor(nPixelsX/nBinsX),y*floor(nPixelsY/nBinsY),fit->GetParameter(1));
			
			delete fit; delete totMapHisto;
		}
	}
	// Plot the map of ToT
	TCanvas* canv = new TCanvas();
	hToTMap->GetXaxis()->SetTitle("Column");
	hToTMap->GetYaxis()->SetTitle("Row");
  hToTMap->GetZaxis()->SetTitle("Mean charge (ToT)");
	hToTMap->SetMaximum(9);
	hToTMap->SetMinimum(5);
	hToTMap->DrawCopy("colz");

  
  std::cout<<"Making per pixel tot map"<<std::endl;
  // Per pixel map
  TH2F* hToTMapPerPixel = new TH2F("hToTMapPerPixel","hToTMapPerPixel",nPixelsX,0,nPixelsX,nPixelsY,0,nPixelsY);
  for(int x=0;x<nPixelsX;x++){
    for(int y=0;y<nPixelsY;y++){
      int id = x + y*nPixelsX;
      // Get the ToT distribution for this location
      string histoname = "/tpanal/ClicpixAnalysis/hClusterTOTAssociated1pixPixel"+makestring(id);
      TH1F* totMapHisto = 0;
      totMapHisto = (TH1F*)gDirectory->Get(histoname.c_str());
      if(totMapHisto == 0) continue;
      if(totMapHisto->GetEntries() < 10) continue;
      
      // Fit the tot distribution and enter the value in the ToT map
      totMapHisto->Fit("gaus","0q");
      TF1* fit = totMapHisto->GetFunction("gaus");
      hToTMapPerPixel->Fill(x,y,fit->GetParameter(1));
      
      delete fit; delete totMapHisto;
    }
  }
  // Plot the map of ToT
  TCanvas* canv2 = new TCanvas();
  hToTMapPerPixel->GetXaxis()->SetTitle("Column");
  hToTMapPerPixel->GetYaxis()->SetTitle("Row");
  hToTMapPerPixel->GetZaxis()->SetTitle("Mean charge (ToT)");
  hToTMapPerPixel->SetMaximum(10);
  hToTMapPerPixel->SetMinimum(0);
  hToTMapPerPixel->DrawCopy("colz");

  std::cout<<"Making pixel response tot maps"<<std::endl;
  // Pixel response maps
  int responseBins=100;
  double responseWidth=0.05;
  TH2F* hToTMapPixelResponseX = new TH2F("hToTMapPixelResponseX","hToTMapPixelResponseX",responseBins+1,-responseWidth,responseWidth,15,0,15);
  for(double rID=0;rID<responseBins;rID++){
    // Get the ToT distribution for this location
    string histoname = "/tpanal/ClicpixAnalysis/hClusterTOTAssociated1pixResponseX"+makestringD(rID);
    TH1F* totMapHisto = 0;
    totMapHisto = (TH1F*)gDirectory->Get(histoname.c_str());
    if(totMapHisto == 0) continue;
    if(totMapHisto->GetEntries() < 1) continue;
 		// loop over all tot bins
    for(int i=0;i<45;i++){
      hToTMapPixelResponseX->Fill(2.*responseWidth*(rID/100.)-responseWidth,i,totMapHisto->GetBinContent(totMapHisto->FindBin(i)));
    }
    delete totMapHisto;
  }
  // Plot the map of ToT
  TCanvas* canv3 = new TCanvas();
  hToTMapPixelResponseX->GetXaxis()->SetTitle("Track distance x (mm)");
  hToTMapPixelResponseX->GetYaxis()->SetTitle("ToT");
  hToTMapPixelResponseX->DrawCopy("colz");

  TH1F* hToTProfilePixelResponseY = new TH1F("hToTProfilePixelResponseY","hToTProfilePixelResponseY",responseBins+1,-responseWidth,responseWidth);
  TH2F* hToTMapPixelResponseY = new TH2F("hToTMapPixelResponseY","hToTMapPixelResponseY",responseBins+1,-responseWidth,responseWidth,16,0,16);
  for(double rID=0;rID<responseBins;rID++){
    // Get the ToT distribution for this location
//.........这里部分代码省略.........

void plotvn(){
    gStyle->SetOptStat(kFALSE);
    int icent = 0;
    int n = 3;
    int color[6] = {1,2,4,7,8,5};
    int style[12] = {20,21,24,25,26,27,29,30,31,32,33,34};
    TGraphErrors *gr[nsub][3][2];
    TGraphErrors *grraw[nsub][3][2];
    TString CNTEP, dire;
    for(int isub=0;isub<nsub;isub++){
    for(int idire=0;idire<3;idire++){
    for(int iCNTEP=0;iCNTEP<1;iCNTEP++){
        if(iCNTEP==0) CNTEP = "NoUseCNTEP";
        if(iCNTEP==1) CNTEP = "UseCNTEP";
        if(idire==0) dire = "";
        if(idire==1) dire = "_east";
        if(idire==2) dire = "_west";
        TString str = choosesub(isub);
        if(str=="ABORT") continue;
        gr[isub][idire][iCNTEP] = new TGraphErrors(Form("Result/%s/v%d_00_%d%s_%s.dat",CNTEP.Data(),n,icent,dire.Data(),str.Data()),"%lg %lg %lg");
        grraw[isub][idire][iCNTEP] = new TGraphErrors(Form("Result/%s/v%draw_00_%d%s_%s.dat",CNTEP.Data(),n,icent,dire.Data(),str.Data()),"%lg %lg %lg");
        SetStyle(*gr[isub][idire][iCNTEP], 1.2, color[idire+3*iCNTEP],style[isub]);
        SetStyle(*grraw[isub][idire][iCNTEP], 1.2, color[idire+3*iCNTEP],style[isub]);
    }
    }
    }
TH1F* h = new TH1F("h","",50,0,5);
h->GetXaxis()->SetRangeUser(0,3.2);
/*
TCanvas *c1 = new TCanvas("c1","c1",800,450);
iCNTEP = 0;
idire = 0;
c1->Divide(2);
c1->cd(1);
h->SetMinimum(0);
h->SetMaximum(0.3);
h->GetXaxis()->SetRangeUser(0,3.2);
//SetTitle(h,"","p_{T}","v_{2}^{raw}");
SetTitle(h,"","p_{T}","v_{2}");
h->DrawCopy();
TLegend *leg = new TLegend(0.2,0.7,0.5,0.9);
leg->SetFillColor(0);
leg->SetBorderSize(0);
for(int ilay = 0;ilay<4;ilay++){
leg->AddEntry(gr[ilay+8][idire][iCNTEP],Form("FVTX layer %d",ilay));
gr[ilay+8][idire][iCNTEP]->Draw("Psame");
}
gr[6][0][1]->Draw("Psame");
leg->AddEntry(gr[6][idire][iCNTEP],Form("FVTX -3.0<#eta<-1.0"));
leg->Draw("Psame");

c1->cd(2);
h->SetMinimum(0.8);
h->SetMaximum(1.2);
SetTitle(h,"","p_{T}","v_{2} ratio");
h->DrawCopy();
for(int ilay = 0;ilay<4;ilay++){
TGraphErrors *grr = (TGraphErrors*)DivideTwoGraphs(gr[ilay+8][idire][iCNTEP],gr[6][idire][iCNTEP]);
SetStyle(*grr,1.2,color[idire+3*iCNTEP],style[ilay+8]);
grr->Draw("Psame");
}
*/
TCanvas *c2 = new TCanvas("c2","c2",450,450);
iCNTEP = 0;
idire = 0;
//c2->Divide(2);
c2->cd(1);
if(n==1){
h->SetMinimum(-0.05);
h->SetMaximum(0.);
}
else if(n==2){
isub=5;
h->SetMinimum(0);
h->SetMaximum(0.3);
}
else if(n==3){
isub=6;
h->SetMinimum(-0.05);
h->SetMaximum(0.2);
}
//SetTitle(h,"","p_{T}","v_{2}^{raw}");
SetTitle(*h,"p_{T} (GeV/c)",Form("v_{%d}",n),"");
c2->SetLeftMargin(0.12);
h->GetXaxis()->CenterTitle();
h->GetYaxis()->CenterTitle();
h->GetYaxis()->SetTitleSize(0.06);
h->DrawCopy();
TLatex t;
t.SetNDC();
t.DrawLatex(0.6,0.82,"d+Au 62GeV");
TLegend *leg = new TLegend(0.2,0.7,0.5,0.88);
leg->SetFillColor(0);
leg->SetBorderSize(0);
leg->SetTextSize(0.04);
SetStyle(*gr[4][idire][iCNTEP], 1.2, 2,style[6]);
SetStyle(*gr[isub][idire][iCNTEP], 1.2, 4,style[6]);
gr[4][idire][iCNTEP]->Draw("Psame");
gr[isub][idire][iCNTEP]->Draw("Psame");
leg->AddEntry(gr[4][idire][iCNTEP],Form("BBCs"),"P");
//.........这里部分代码省略.........

//.........这里部分代码省略.........
          "galice.root");
    return;
  }
  runLoader->LoadHeader();
  runLoader->LoadKinematics();
  runLoader->LoadTrackRefs();

  AliLoader *dl = runLoader->GetDetectorLoader("ITS");

  //Trackf
  TTree *trackRefTree = 0x0; 
  TClonesArray *trackRef = new TClonesArray("AliTrackReference",1000);

  //  TH1F *hRef = new TH1F("","",100,0,100);
  TH1F *hR = new TH1F("","",100,0,100);
  if (hAllMC==0) hAllMC = new TH1F("","",100,0.1,2);
  Float_t ptmin = hAllMC->GetBinCenter(1)-hAllMC->GetBinWidth(1)/2;
  Float_t ptmax = hAllMC->GetBinCenter(hAllMC->GetNbinsX())+hAllMC->GetBinWidth(hAllMC->GetNbinsX())/2;
  //  Int_t nAllMC = 0;

  // Detector geometry
  TArrayD rmin(0);   TArrayD rmax(0); 
  GetDetectorRadii(&rmin,&rmax);
  TArrayI nLaySigs(rmin.GetSize());

  printf("Counting trackable MC tracks ...\n");
  
  for(Int_t iEv =0; iEv<runLoader->GetNumberOfEvents(); iEv++){
    Int_t nTrackableTracks = 0;
    runLoader->GetEvent(iEv);
    AliStack* stack = runLoader->Stack();  
    printf("+++ event %i (of %d) +++++++++++++++++++++++  # total MCtracks: %d \n",iEv,runLoader->GetNumberOfEvents()-1,stack->GetNtrack());

    trackRefTree=runLoader->TreeTR();
    TBranch *br = trackRefTree->GetBranch("TrackReferences");
    if(!br) {
      printf("no TR branch available , exiting \n");
      return;
    }
    br->SetAddress(&trackRef);

    // init the trackRef tree 
    trackRefTree=runLoader->TreeTR();
    trackRefTree->SetBranchAddress("TrackReferences",&trackRef);
 
    // Count trackable MC tracks
    for (Int_t iMC=0; iMC<stack->GetNtrack(); iMC++) {

      TParticle* particle = stack->Particle(iMC); 
      if (TMath::Abs(particle->Eta())>etaCut) continue;
      if (onlyPrims && !stack->IsPhysicalPrimary(iMC)) continue;
      if (onlyPion && TMath::Abs(particle->GetPdgCode())!=211) continue;


      Bool_t isTrackable = 0;
      nLaySigs.Reset(0);
 
      trackRefTree->GetEntry(stack->TreeKEntry(iMC));
      Int_t nref=trackRef->GetEntriesFast();
      for(Int_t iref =0; iref<nref; iref++){
	AliTrackReference *trR = (AliTrackReference*)trackRef->At(iref);
	if(!trR) continue;
	if(trR->DetectorId()!=AliTrackReference::kITS) continue;
	Float_t radPos = trR->R();
	hR->Fill(radPos);
	for (Int_t il=0; il<rmin.GetSize();il++) {
	  if (radPos>=rmin.At(il)-0.1 && radPos<=rmax.At(il)+0.1) {
	    //	    cout<<"  in Layer "<<il<<" "<<radPos;
	    nLaySigs.AddAt(1.,il);
	    //	    cout<<" "<<nLaySigs.At(il)<<endl;
	  }
	}
      }

      if (nLaySigs.GetSum()>=3) {
	isTrackable =1;
	//	cout<<nLaySigs.GetSum()<<endl;
      }
      
      if (isTrackable) {
	Double_t ptMC = particle->Pt();
	//	Double_t etaMC = particle->Eta();
	//	if (ptMC>ptmin&&ptMC<ptmax) {nTrackableTracks++;hAllMC->Fill(ptMC);}
	if (ptMC>ptmin) {nTrackableTracks++;hAllMC->Fill(ptMC);}

      }

      
    } // entries tracks MC
    printf(" -> trackable MC tracks: %d (%d)\n",nTrackableTracks,hAllMC->GetEntries());
  }//entries Events
  

  hR->DrawCopy();
  hAllMC->DrawCopy();
  runLoader->UnloadHeader();
  runLoader->UnloadKinematics();
  delete runLoader;

}

void plotMerged(Bool_t onlyPlot=0) {

  gStyle->SetPalette(1);
 
  TFile f("histoSum.root","UPDATE");

  TH1F* hAllMC = f.Get("allMC");
  TH1F* hAllFound= f.Get("allFound");
  TH1F* hImperfect= f.Get("imperfect");
  TH1F* hPerfect= f.Get("perfect");
  TH1F* hNoMCTrack= f.Get("noMCtrack");
  
  
  // have to be recalculated
  TH1F* hPurity = f.Get("purity");
  TH1F* hEff= f.Get("efficiency");
  TH1F* hFake= f.Get("fake");
  TH1F* hAnna= f.Get("annaEff");

  TH2D* h2Ddca= f.Get("dca2D");
  TGraphErrors *d0= f.Get("dca");

  TH2D* h2Dpt= f.Get("dPt2D");
  TGraphErrors *gPt= f.Get("dPt");


  if (!onlyPlot) {
    /*    // Get Errors right
    hAllMC->Sumw2();   
    hAllFound->Sumw2();
    hPerfect->Sumw2(); 
    hImperfect->Sumw2(); 
    h2Dpt->Sumw2();
    h2Ddca->Sumw2();
    */

    // Efficiencies - and normalize to 100%
    
    TF1 f1("f1","100+x*0",0.,1.e3);
    
    hPurity->Divide(hPerfect,hAllFound,1,1,"b"); 
    hPurity->Multiply(&f1);
    hPurity->SetMarkerColor(kGreen);
    hPurity->SetMarkerStyle(21);
    hPurity->GetXaxis()->SetTitle("transverse momentum p_{t} (GeV)");
    hPurity->SetStats(0);
    
    hPurity->GetYaxis()->SetRangeUser(0,100);
    hPurity->SetTitle("Efficiency & Purity");
    
    hEff->Divide(hPerfect,hAllMC,1,1,"b");
    hEff->Multiply(&f1);
    hEff->GetXaxis()->SetTitle("transverse momentum p_{t} (GeV)");
    hEff->SetMarkerColor(kBlue);
    hEff->SetMarkerStyle(21);
    hEff->SetStats(0);
    
    hFake->Divide(hImperfect,hAllMC,1,1,"b");
    hFake->Multiply(&f1);
    hFake->GetXaxis()->SetTitle("transverse momentum p_{t} (GeV)");
    hFake->SetMarkerColor(kRed);
    hFake->SetMarkerStyle(21);
    hFake->SetStats(0);
    
    hAnna->Divide(hAllFound,hAllMC,1,1,"b");
    hAnna->Multiply(&f1);
    hAnna->GetXaxis()->SetTitle("transverse momentum p_{t} (GeV)");
    hAnna->SetMarkerColor(kBlack);
    hAnna->SetMarkerStyle(21);
    hAnna->SetStats(0);
    
    
    // Impact parameter resolution ---------------
    TCanvas *c3 = new TCanvas("c3","impact");//,200,10,900,900);
    c3->Divide(2,1); c3->cd(1);
    h2Ddca->DrawCopy("colz");
    h2Ddca->FitSlicesY() ;
    TH2D *dcaM = (TH2D*)gDirectory->Get("dca2D_1"); dcaM->Draw("same");
    TH2D *dcaRMS = (TH2D*)gDirectory->Get("dca2D_2"); //dcaRMS->Draw();
    TGraphErrors *d0 = new TGraphErrors(); 
    for (Int_t ibin =1; ibin<=dcaRMS->GetXaxis()->GetNbins(); ibin++) {
      d0->SetPoint(     ibin-1,dcaRMS->GetBinCenter(ibin),dcaRMS->GetBinContent(ibin)*1e4); // microns
      d0->SetPointError(ibin-1,0,dcaRMS->GetBinError(ibin)*1e4); // microns
    }
    d0->SetMarkerStyle(21);
    d0->SetMaximum(200);  d0->SetMinimum(0);
    d0->GetXaxis()->SetTitle("transverse momentum p_{t} (GeV)");
    d0->GetYaxis()->SetTitle("R-#phi Pointing Resolution (#mum)");
    d0->SetName("dca");  d0->SetTitle("DCAvsPt");
    //  c3->cd(1);  h2Ddca->Draw("surf2");
    c3->cd(2);  d0->Draw("APE");
    
    // PT RESOLUTION ------------
    TCanvas *c4 = new TCanvas("c4","pt resolution");//,200,10,900,900);  
    c4->Divide(2,1); c4->cd(1);
    h2Dpt->DrawCopy("colz");
    h2Dpt->FitSlicesY() ;
    TH2D *dPtM = (TH2D*)gDirectory->Get("dPt2D_1"); dPtM->Draw("same");
    TH2D *dPtRMS = (TH2D*)gDirectory->Get("dPt2D_2"); // dPtRMS->Draw("");
    TGraphErrors *gPt = new TGraphErrors(); 
//.........这里部分代码省略.........

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

	Double_t x, y, chtmp, x1, x2, y1, y2;

	/*TList *grCellCalibList = OpenCalibFile("CalibrationData1000events.root");

	TGraphErrors *grCellCalib;
	TGraphErrors *grDataSubCalib = new TGraphErrors(DOMINO_NCELL);
	grDataSubCalib->SetTitle("Data after calibration correction");
	grDataSub = (TGraphErrors *) grDataSubList->At(anaChannel);


	for(ch = 0; ch < DOMINO_NCELL; ch++) {
		grCellCalib = ((TGraphErrors *) grCellCalibList->At(ch));
		grCellCalib->Fit("pol3", "Q");
		TF1 *pol3fit = ((TF1 *) grCellCalib->GetFunction("pol3"));
		grDataSub->GetPoint(ch, x, y);
		chtmp = y - (Double_t)(pol3fit->Eval(y/3.25));
		grDataSubCalib->SetPoint(ch, x, chtmp);
	}

	TCanvas *cGrTest = new TCanvas("grTest", "test per vedere i dati", 1000,1000);

	grDataSubCalib->Draw("APEL");*/


	TString Title = "Charge Distribution per channel";
	gStyle->SetOptFit(111);
	TCanvas *cdistch = new TCanvas("cdistch", Title, 1000, 1000);
	cdistch->Divide(3, 3);
	for (int i = 0; i < DOMINO_NCH; i++) {
		cdistch->cd(i + 1);
		TH1 *dhist = (TH1 *) DistChList->At(i);
		dhist->DrawCopy();
		dhist->SetLineWidth(1);
		dhist->Fit("gaus", "Q");
		dhist->GetFunction("gaus")->SetLineColor(4);
		dhist->GetFunction("gaus")->SetLineWidth(2);
	}

	TString Title = "Charge Distribution Pedestals Subtracted per channel";
	TCanvas *cdistchsub = new TCanvas("cdistchsub", Title, 1000, 1000);
	cdistchsub->Divide(3, 3);
	for (int i = 0; i < DOMINO_NCH; i++) {
		cdistchsub->cd(i + 1);
		TH1 *dsubhist = (TH1 *) DistChSubList->At(i);
		dsubhist->DrawCopy();
		dsubhist->SetLineWidth(1);
		dsubhist->Fit("gaus", "Q");
		dsubhist->GetFunction("gaus")->SetLineColor(4);
		dsubhist->GetFunction("gaus")->SetLineWidth(2);
	}

	TString Title = "Charge Distribution Pedestals and Ch0 Subtracted per channel";
	TCanvas *cdistch0sub = new TCanvas("cdistch0sub", Title, 1000, 1000);
	cdistch0sub->Divide(3, 3);
	for (int i = 0; i < DOMINO_NCH; i++) {
		cdistch0sub->cd(i + 1);
		TH1 *dch0subhist = (TH1 *) DistCh0SubList->At(i);
		dch0subhist->DrawCopy();
		dch0subhist->SetLineWidth(1);
		dch0subhist->Fit("gaus", "Q");
		dch0subhist->GetFunction("gaus")->SetLineColor(4);
		dch0subhist->GetFunction("gaus")->SetLineWidth(2);
	}

//.........这里部分代码省略.........
  gStyle->SetOptStat(0);
  gStyle->SetOptTitle(0);

  c1->Range(-38.59656,-12.02472,644.9871,2.486126);
  c1->SetFillColor(0);
  c1->SetBorderMode(0);
  c1->SetBorderSize(2);
  c1->SetTickx(1);
  c1->SetTicky(1);
  c1->SetLeftMargin(0.1552174);
  c1->SetRightMargin(0.02034783);
  c1->SetTopMargin(0.08166969);
  c1->SetBottomMargin(0.1288566);
  c1->SetFrameFillStyle(0);
  c1->SetFrameBorderMode(0);

  c1->cd();
  c1->cd()->SetLogy(); 
 

  TH1F *hDum1 =  new TH1F("hdum1","hdum1",13,StartPoint,EndPoint);
  hDum1->SetMinimum(1.E3);
  hDum1->SetMaximum(1.E8);
  hDum1->GetXaxis()->SetTitle("p_{T} [GeV/c]");
  hDum1->GetXaxis()->SetTitleSize(0.05);
  hDum1->GetXaxis()->SetLabelSize(0.045);
  hDum1->GetXaxis()->SetTitleOffset(1.18);
  hDum1->GetXaxis()->CenterTitle();
  hDum1->GetYaxis()->SetTitle("dN_{jet}/dp_{T} [(GeV/c)^{-1}]");
  hDum1->GetYaxis()->SetTitleSize(0.05);
  hDum1->GetYaxis()->SetLabelSize(0.045);
  hDum1->GetYaxis()->SetTitleOffset(1.45);
  hDum1->GetYaxis()->CenterTitle();
  hDum1->DrawCopy();

  hjet_Rebin_Jet20->GetXaxis()->SetRangeUser(StartPoint+1,EndPoint-2);
  hjet_Rebin_Jet20->SetMarkerStyle(20);
  hjet_Rebin_Jet20->SetMarkerColor(1);
  hjet_Rebin_Jet20->SetMarkerSize(1.0);
  hjet_Rebin_Jet20->SetLineColor(1);
  hjet_Rebin_Jet20->DrawCopy("same");

  hjet_Rebin_Jet40->GetXaxis()->SetRangeUser(StartPoint+1,EndPoint-2);
  hjet_Rebin_Jet40->SetMarkerStyle(20);
  hjet_Rebin_Jet40->SetMarkerColor(2);
  hjet_Rebin_Jet40->SetMarkerSize(1.0);
  hjet_Rebin_Jet40->SetLineColor(2);
  hjet_Rebin_Jet40->DrawCopy("same");

  for(int itrig=2; itrig<5; itrig++)
  {	 	  
    hjet[itrig]->Rebin(2);
    normalizeByBinWidth(hjet[itrig]);
    hjet[itrig]->GetXaxis()->SetRangeUser(StartPoint+1,EndPoint-2);	   
    hjet[itrig]->SetMarkerStyle(20);
    hjet[itrig]->SetMarkerColor(1+itrig);
    hjet[itrig]->SetMarkerSize(1.0);
    hjet[itrig]->SetLineColor(1+itrig);
    hjet[itrig]->DrawCopy("same");
  }
  hcombine->Rebin(2);
  normalizeByBinWidth(hcombine);
  hcombine->GetXaxis()->SetRangeUser(StartPoint+1,EndPoint-2);
  hcombine->SetMarkerStyle(24);
  hcombine->SetMarkerColor(1);
  hcombine->SetMarkerSize(1.5);

void plotFrac(TList* HistList, TH1F* compT, std::string name, bool fracAll){

  gROOT->SetBatch();
  system("mkdir -p plots/ada/fracs");
  system("mkdir -p plots/grad/fracs");

  std::string bdt;
  TString str = compT->GetName();
  if (str.Contains("ada")) bdt="ada";
  if (str.Contains("grad")) bdt="grad";

  int nHists=HistList->GetEntries();
  TH1F *comp = linearBin(compT);
  TH1F *uandd[nHists];
  
  gROOT->SetStyle("Plain");
  gROOT->ForceStyle();
  gStyle->SetOptStat(0);
  TCanvas *canv = new TCanvas("","",700,700);
  
  canv->Divide(1,2,0.1,0.1);
  canv->cd(1);
  gPad->SetPad(0.01,0.3,0.99,0.99);
  gPad->SetBottomMargin(0.0);
  gPad->SetLeftMargin(0.1);
  comp->SetLineColor(1);
  comp->SetFillColor(kGray);
  comp->Draw("hist");
  comp->SetTitle(("Up, down and interpolated templates for "+bdt+" "+name).c_str());
  comp->GetYaxis()->SetRangeUser(0.0,((TH1F*)HistList->At(0))->GetMaximum()+0.5);
  comp->GetYaxis()->SetTitle("Events / bin");
  comp->GetXaxis()->SetTitle("Category");
  int mass;
  if (name=="syst120") mass=120;
  if (name=="syst135") mass=135;
  TLegend *leg = new TLegend(0.6,0.5,0.88,0.88);
  leg->SetLineColor(0);
  leg->SetFillColor(0);
  if (name=="syst120" || name=="syst135"){
    leg->AddEntry(comp,Form("True %d signal binned at %d",mass,mass),"f");
  }
  
  canv->cd(2);
  gPad->SetPad(0.01,0.01,0.99,0.3);
  gPad->SetTopMargin(0.0);
  gPad->SetBottomMargin(0.2);
  gPad->SetLeftMargin(0.1);
  TF1 *line = new TF1("line","0.0",0.,comp->GetBinLowEdge(comp->GetNbinsX()+1));
  line->SetLineColor(kBlack);
  
  for (int i=0; i<nHists; i++){
    TH1F *temp = linearBin((TH1F*)HistList->At(i));
    temp->SetLineColor((i+1)*2);
    if (!fracAll && i==nHists-1) temp->SetLineColor(kBlack); 
    canv->cd(1);
    temp->DrawCopy("same e");
    uandd[i]= (TH1F*)temp->Clone();
    uandd[i]->Add(comp,-1);
    uandd[i]->Divide(temp);
    uandd[i]->Scale(-100.);
    //uandd[i]->Divide(comp);
    if (name=="syst120" || name=="syst135") {
      if (i==0) leg->AddEntry(uandd[i],Form("True %d signal binned at %d",mass+5,mass),"lep");
      if (i==1) leg->AddEntry(uandd[i],Form("True %d signal binned at %d",mass-5,mass),"lep");
      if (i==nHists-1) leg->AddEntry(uandd[i],Form("Interpolated %d signal",mass),"lep");
    }
    //uandd[i]->Scale(100.0);
    canv->cd(2);
    if (fracAll){
      if (i==0) uandd[i]->Draw("e");
      else {
        uandd[i]->Draw("same e");
        line->Draw("same");
      }
    }
    else {
      if (i==nHists-1) {
        uandd[i]->Draw("e");
        line->Draw("same");
      }
    }
    uandd[i]->SetTitle("");
    uandd[i]->GetYaxis()->SetLabelSize(0.08);
    uandd[i]->GetYaxis()->SetRangeUser(-50.,50.);
    uandd[i]->GetYaxis()->SetTitle("#frac{#Delta_{int}}{int} %");
    uandd[i]->GetYaxis()->SetTitleOffset(0.4);
    uandd[i]->GetYaxis()->SetTitleSize(0.08);
    uandd[i]->GetXaxis()->SetLabelSize(0.08);
    uandd[i]->GetXaxis()->SetTitle("BDT output");
    uandd[i]->GetXaxis()->SetTitleSize(0.08);
    canv->cd(1);
    if (name=="syst120" || name=="syst135") leg->Draw("same");
  }
  

  canv->Print(("plots/"+bdt+"/fracs/"+name+".png").c_str(),"png");
  delete canv;
  
  fracCalls++;
  
//.........这里部分代码省略.........

//_____________________________________________________________________________
void ProofSimple::Terminate()
{
   // The Terminate() function is the last function to be called during
   // a query. It always runs on the client, it can be used to present
   // the results graphically or save the results to file.

   //
   // Create a canvas, with 100 pads
   //
   TCanvas *c1 = (TCanvas *) gDirectory->FindObject("c1");
   if (c1) {
      gDirectory->Remove(c1);
      delete c1;
   }
   c1 = new TCanvas("c1","Proof ProofSimple canvas",200,10,700,700);
   Int_t nside = (Int_t)TMath::Sqrt((Float_t)fNhist);
   nside = (nside*nside < fNhist) ? nside+1 : nside;
   c1->Divide(nside,nside,0,0);

   Bool_t tryfc = kFALSE;
   TH1F *h = 0;
   for (Int_t i=0; i < fNhist; i++) {
      if (!(h = dynamic_cast<TH1F *>(TProof::GetOutput(Form("h%d",i), fOutput)))) {
         // Not found: try TFileCollection
         tryfc = kTRUE;
         break;
      }
      c1->cd(i+1);
      h->DrawCopy();
   }

   // If the histograms are not found they may be in files: is there a file collection?
   if (tryfc && GetHistosFromFC(c1) != 0) {
      Warning("Terminate", "histograms not found");
   } else {
      // Final update
      c1->cd();
      c1->Update();
   }

   // Analyse hlab, if there
   if (fHLab && !gROOT->IsBatch()) {
      // Printout
      Int_t nb = fHLab->GetNbinsX();
      if (nb > 0) {
         Double_t entb = fHLab->GetEntries() / nb;
         if (entb) {
            for (Int_t i = 0; i < nb; i++) {
               TString lab = TString::Format("hl%d", i);
               Int_t ib = fHLab->GetXaxis()->FindBin(lab);
               Info("Terminate","  %s [%d]:\t%f", lab.Data(), ib, fHLab->GetBinContent(ib)/entb); 
            }
         } else
            Warning("Terminate", "no entries in the hlab histogram!");
      }
   }

   // Process the ntuple, if required
   if (fHasNtuple != 1 || !fPlotNtuple) return;

   if (!(fNtp = dynamic_cast<TNtuple *>(TProof::GetOutput("ntuple", fOutput)))) {
      // Get the ntuple from the file
      if ((fProofFile =
            dynamic_cast<TProofOutputFile*>(fOutput->FindObject("SimpleNtuple.root")))) {

         TString outputFile(fProofFile->GetOutputFileName());
         TString outputName(fProofFile->GetName());
         outputName += ".root";
         Printf("outputFile: %s", outputFile.Data());

         // Read the ntuple from the file
         fFile = TFile::Open(outputFile);
         if (fFile) {
            Printf("Managed to open file: %s", outputFile.Data());
            fNtp = (TNtuple *) fFile->Get("ntuple");
         } else {
            Error("Terminate", "could not open file: %s", outputFile.Data());
         }
         if (!fFile) return; 

      } else {
         Error("Terminate", "TProofOutputFile not found");
         return;
      }
   }
   // Plot ntuples
   if (fNtp) PlotNtuple(fNtp, "proof ntuple");
}

void plotfigure4(bool nodata = false, bool prelim = false, bool AMPT=true, bool SONIC=true, 
                 bool superSONIC=true, bool IPGLASMA=true, bool BOZEK=true) {

  gROOT->SetStyle("Plain");
  gStyle->SetOptStat(0);
  gStyle->SetOptTitle(0);

  TCanvas *c1 = new TCanvas("c1","c1",10,10,700,700);

  // create template histogram
  TH1F *htemplate = new TH1F("htemplate","htemplate",10,0.0,3.5);
  htemplate->SetMaximum(0.30); //25);
  htemplate->SetXTitle("p_{T} [GeV/c]");
  //  htemplate->SetYTitle("v_{n} Coefficient");
  htemplate->SetYTitle("v_{n}");
  htemplate->GetYaxis()->SetTitleOffset(1.2);
  htemplate->DrawCopy();

  // grab AMPT points (http://arxiv.org/abs/1501.06880)
  TFile *fin_ampt = new TFile("AMPT_3HeAu_central_vn_new.root");
  // copy into new TGraphs for easier manipulcation
  //TGraphErrors *tamptv2 = gv2->Clone();
  //TGraphErrors *tamptv3 = gv3->Clone();

  TGraphErrors *tamptv22 = gv2->Clone();
  TGraphErrors *tamptv32 = gv3->Clone();
  
  TGraph *tamptv2 = new TGraph();
  TGraph *tamptv3 = new TGraph();
  for (int i=0;i<gv2->GetN();i++) {
    double x; double y;
    gv2->GetPoint(i,x,y);
    tamptv2->SetPoint(i,x,y);
  }
  for (int i=0;i<gv3->GetN();i++) {
    double x; double y;
    gv3->GetPoint(i,x,y);
    tamptv3->SetPoint(i,x,y);
  }
  
  

  
  
  tamptv22->SetFillColor(42);
  tamptv32->SetFillColor(42);

  tamptv22->SetLineColor(46);
  tamptv22->SetLineWidth(3);
  tamptv22->SetLineStyle(1);
  if (AMPT) tamptv22->Draw("3,same");
  tamptv32->SetLineColor(46);
  tamptv32->SetLineWidth(3);
  tamptv32->SetLineStyle(1);
  if (AMPT) tamptv32->Draw("3,same");

  tamptv2->SetLineColor(46);
  tamptv2->SetLineWidth(3);
  tamptv2->SetLineStyle(1);
  if (AMPT) tamptv2->Draw("l,same");
  tamptv3->SetLineColor(46);
  tamptv3->SetLineWidth(3);
  tamptv3->SetLineStyle(1);
  if (AMPT) tamptv3->Draw("l,same");
  

  // IPGlasma curves (http://arxiv.org/abs/1407.7557)
  TGraph *tipglasmav2 = new TGraph("ipglasmaflowv2.dat","%lg, %lg");
  tipglasmav2->SetLineColor(6);
  tipglasmav2->SetLineWidth(3);
  tipglasmav2->SetLineStyle(1);
  if (IPGLASMA) tipglasmav2->Draw("l,same");
  TGraph *tipglasmav3 = new TGraph("ipglasmaflowv3.dat","%lg, %lg");
  tipglasmav3->SetLineColor(6);
  tipglasmav3->SetLineWidth(3);
  tipglasmav3->SetLineStyle(1);
  if (IPGLASMA) tipglasmav2->Draw("l,same");
  if (IPGLASMA) tipglasmav3->Draw("l,same"); 

  // Now Paul Romatschke result without and with preflow
  // Refs:  http://arxiv.org/abs/1502.04745

  double sonic_beta2[100][13];
  double sonic_beta3[100][13];
  double supersonic_beta2[100][13];
  double supersonic_beta3[100][13];
  // pt:yield:yields:v1:v1s:v2:v2s:v3:v3s:v4:v4s:v5:v5s 

  ifstream myfile;
  myfile.open("res_he3config-42mb-s63-sigma0.4-scal2.00-beta2-nopf-av_unid_vn.dat");
  int j=0;
  string dummyLine;
  getline(myfile, dummyLine); // skip first comment line
  while (myfile) {
    myfile >> sonic_beta2[j][0] >> sonic_beta2[j][1] >> sonic_beta2[j][2] >>
		   sonic_beta2[j][3] >> sonic_beta2[j][4] >> sonic_beta2[j][5] >>
		   sonic_beta2[j][6] >> sonic_beta2[j][7] >> sonic_beta2[j][8] >>
		   sonic_beta2[j][9] >> sonic_beta2[j][10] >> sonic_beta2[j][11] >>
		   sonic_beta2[j][12];
    cout << "pt = " << sonic_beta2[j][0] << " v2 = " << sonic_beta2[j][5] << endl;
//.........这里部分代码省略.........

void scan2() {

    /*  TFile f1("hist.root");
      TFile f2("SysTot.root");

      TH1D *datahist = (TH1D*)f1.Get("datahist");
      TH1D *fithist = (TH1D*)f1.Get("fithist");

      TH1D *totlow = (TH1D*)f2.Get("totlow");
      TH1D *tothigh = (TH1D*)f2.Get("tothigh");
    */

    TF1 *fun1 = new TF1("fun",fun,0,1000,3);

    TFile f("Output.root");
    TH1F *datahist = dynamic_cast<TH1F*> (f.Get("InvMass"));
    TH1F *fithist = dynamic_cast<TH1F*> (f.Get("hist_err"));
    TH1F *tothigh = dynamic_cast<TH1F*> (f.Get("SystPlus"));
    TH1F *totlow = dynamic_cast<TH1F*> (f.Get("SystMinus"));
    new TCanvas();
    datahist->DrawCopy();
    fithist->DrawCopy("same");
    gPad->SetEditable(0);

    /*	if (datahist->IsZombie()) { std::cout << "Err! datahist not found " << std::endl; return; }
    	if (fithist->IsZombie()) { std::cout << "Err! err_hist not found " << std::endl; return; }
    	if (tothigh->IsZombie()) { std::cout << "Err! syst_plus hist not found " << std::endl; return; }
    	if (totlow->IsZombie()) { std::cout << "Err! syst_minus hist not found " << std::endl; return; }
    */

    TH1D *p = new TH1D("p","Probability",100,150,900);

    for (int i = 1; i<=fithist->GetNbinsX(); i++)
    {
        double mass = fithist->GetBinCenter(i);

        if (mass>440 && mass<450)
        {
            double sigma = sqrt(2.0)*sqrt(pow(0.135*sqrt(mass/2.0),2)+
                                          pow(0.02*(mass/2.0),2));
            cout<< " ===================== mass +/- sigma = " << mass<<"+/-"<<sigma<<endl;
            int bin1 = fithist->FindBin(mass-sigma/2);
            int bin2 = fithist->FindBin(mass+sigma/2);
            cout<<mass<<" "<<bin1<<" "<<bin2<<endl;
            double data = 0;
            double bg = 0;
            double err = 0;

            for (int j = bin1; j<=bin2; j++)
            {
                data+=datahist->GetBinContent(j);
                bg+=fithist->GetBinContent(j);
                double err1 = -totlow->GetBinContent(j);
                double err2 = tothigh->GetBinContent(j);
                err+=TMath::Max(err1,err2)*bg; //why multiply by bg???
            }
            cout << "Total Data/Bg+/-err in mass window[" << mass << "] = "<< data <<"/ "<< bg << "+/-" << err <<endl;
            double prob = 0;
            fun1->SetParameter(0,bg);
            fun1->SetParameter(1,err);

            for (int j = int(data+0.001); j<100; j++) {
                fun1->SetParameter(2,j);
                //fun1->Print();
                //cout << "Evaluating Intrgral for j = " << j << " from x0= " << TMath::Max(0.0,bg-10*err) << " to x1 = " << bg+10*err << endl;
                double val = fun1->Integral(TMath::Max(0.0,bg-10*err),bg+10*err);
                //double val = fun1->Integral(TMath::Max(0.0,bg-2*err),bg+2*err);
                /*for (int z=TMath::Max(0.0,bg-2*err); z < bg+2*err; ++z)
                {
                	if (c<4)
                	{
                	std::cout << "func at [" << z << "]=" << 	fun1->Eval(z) << std::endl;
                	}
                }
                */
                prob += val;
            }
            cout<< "Prob for mass[" << mass<<"]="<< prob <<endl;
            p->SetBinContent(p->FindBin(mass),prob);
        }
    }
    /*
    delete gRandom;
    gRandom = (TRandom*) new TRandom3;
    gRandom->SetSeed(3);

    TH1D *minp = new TH1D("minp","Minimum Probability of Each PseudoExpt",100,0,0.2);

    //int nexp = 50000;
    int nexp = 10;

    TH1D *htemp = (TH1D*)datahist->Clone("htemp");

    for (int iexp = 0; iexp<nexp; iexp++){
      //if (iexp%10==0) cout<<iexp<<endl;
      //generate pseudo-experiments
      htemp->Reset();
      for (int i = 1; i<=htemp->GetNbinsX(); i++){
        double mass = htemp->GetBinCenter(i);
        if (mass>150&&mass<650){
//.........这里部分代码省略.........