C++ TH3D类代码示例

// 
// GetHistFile
//
// Sets pointers to the scaler and random-subtracted asymmetry 3D histograms
// for a certain input file (full/empty and perp/para).  Results are NOT
// dead-time corrected.
//
void GetHistFile( TFile* file, Bool_t cthflag = kFALSE)
{
	Double_t pa;
	TString prompt, random, scalers;

	TH3D* hP;
	TH3D* hR;

	pa = 0.0833;

	if ( cthflag == kFALSE ) {
//		prompt = "PhiCMCut2P_v_ThetaCMCut2P_v_TChanCut2P";
//		random = "PhiCMCut2R_v_ThetaCMCut2R_v_TChanCut2R";
		prompt = "PhiCMCut1P_v_ThetaCMCut1P_v_TChanCut1P";
		random = "PhiCMCut1R_v_ThetaCMCut1R_v_TChanCut1R";
	}
	else {
		prompt = "PhiCMCut1P_v_CosThetaCMCut1P_v_TChanCut1P";
		random = "PhiCMCut1R_v_CosThetaCMCut1R_v_TChanCut1R";
	}
	scalers = "SumScalers152to503";

	hP = (TH3D*) file->Get( prompt);
	hR = (TH3D*) file->Get( random);
	hsc = (TH1D*) file->Get( scalers);

	hS = (TH3D*) hP->Clone( "sub");
	hS->Sumw2();
	hS->Add( hR, -pa);
}

TH3D *plotCor(TTree *tSig,TString var,TString weight, TString etcut, TString htcut,TString njetcut, TString ph, TString bin)
{
  TCut etcut20=weight+" * ("+etcut+htcut+" && "+njetcut+")";

  TCut conv1=weight+" * ("+etcut+htcut+" &&"+njetcut+" && mHits<=1)";
  TCut conv2=weight+" * ("+etcut+htcut+" && "+njetcut+"&& mHits<=1 && mHits<2)";
  TCut conv3=weight+" * ("+etcut+htcut+" && "+njetcut+"&& mHits<=1 && mHits<2 && (dcot>0.02 || dcot<-0.02 || dist>0.02 || dist<-0.02))";
  
  tSig->Draw(var+">>hsqrti"+bin,conv1);
  TH3D *hsqrt = (TH3D*)gDirectory->Get("hsqrti");
  hsqrt->SetName(ph+var+htcut+"hsqrti");

  return hsqrt;
}

// ============================================================================
StatusCode Gaudi::Utils::Histos::fromXml
( TH3D& result , const std::string& input )
{
  //
  result.Reset() ;                                 // RESET old histogram
  //
  _Xml<TH3D> _xml ;
  std::auto_ptr<TH3D> histo = _xml ( input ) ;
  if ( 0 == histo.get() ) { return StatusCode::FAILURE ; }        // RETURN
  //
  result.Reset() ;
  histo->Copy ( result ) ;
  //
  return StatusCode::SUCCESS ;
}

void DoEffCor3D(TFile* a, int b, double c, double d, double e, double *f){
  char tmp[512];
  if(b == 0) sprintf(tmp,"eff_default");
  if(b == 1) sprintf(tmp,"eff_cowboy");
  if(b == 2) sprintf(tmp,"eff_sailor");
  TH3D *hEff = (TH3D*)a->Get(tmp);
  int sbin = 0;
  sbin = hEff->FindBin(c, d, e);
  f[0] = hEff->GetBinContent(sbin);
  f[1] = hEff->GetBinError(sbin);
  if(f[0] == 0) {
    cout<<"b : "<<b<<", c : "<<c<<", d : "<<d<<", e : "<<e<<endl;
    cout<<"sbin : "<<sbin<<", eff : "<<f[0]<<", err : "<<f[1]<<endl;
  }
  //return f;
}

int main ()
{
  TFile* testFile = new TFile("test.root", "RECREATE");
  
  TH3D* test = new TH3D("test","test", 50,-50,50, 50,-50,50, 50,-50,50);
  
  for (int i = -10; i <= 10; i++)
    {
      for (int j = -2; j <= 2; j++)
	{
	  for (int k = -2; k <= 2; k++)
	    {
	      test->Fill (i, j, k, 5);
	    }
	}
    }
  
  testFile->Write();

  return 0;
}

TH3D* readHist3D(TString nameHist,TString nameFile, int rebin)
{
 TFile* file = new TFile(nameFile);

 TH3D* hist = (TH3D*)file->Get(nameHist);
 hist->GetSumw2();
 // hist->SetLineWidth(2);
 // if(rebin>0) hist->RebinX(rebin); hist->RebinY(rebin);
 hist->GetXaxis()->SetTitleSize(.055);
 hist->GetYaxis()->SetTitleSize(.055);
 hist->GetXaxis()->SetLabelSize(.05);
 hist->GetYaxis()->SetLabelSize(.05);
 hist->SetStats(kFALSE);
 return hist;
}

void projectPbPbMBMC()
{
  gStyle->SetTextSize(0.05);
  gStyle->SetTextFont(42);
  gStyle->SetPadRightMargin(0.043);
  gStyle->SetPadLeftMargin(0.18);
  gStyle->SetPadTopMargin(0.1);
  gStyle->SetPadBottomMargin(0.145);
  gStyle->SetTitleX(.0f);
  gStyle->SetOptFit(1111);
  gStyle->SetOptStat(0);
  gStyle->SetOptTitle(0);

  TFile* infMCP = new TFile("/data/HeavyFlavourRun2/MC2015/Dntuple/PbPb/ntD_EvtBase_20160513_DfinderMC_PbPb_20160502_dPt1tkPt0p5_D0_prompt_Dpt2Dy1p1tkPt0p7tkEta2Decay2p9Dalpha0p14Skim_pthatweight.root");
  TFile* infMCNP = new TFile("/data/HeavyFlavourRun2/MC2015/Dntuple/PbPb/ntD_EvtBase_20160513_DfinderMC_PbPb_20160502_dPt1tkPt0p5_D0_nonprompt_Dpt2Dy1p1tkPt0p7tkEta2Decay2p9Dalpha0p14Skim_pthatweight.root");
  
  TTree* tMCP = (TTree*)infMCP->Get("ntDkpi");
  TTree* tMCNP = (TTree*)infMCNP->Get("ntDkpi");
  TTree* tMCPHI = (TTree*)infMCP->Get("ntHi");
  TTree* tMCNPHI = (TTree*)infMCNP->Get("ntHi");
  TTree* tMCPSkim = (TTree*)infMCP->Get("ntSkim");
  TTree* tMCNPSkim = (TTree*)infMCNP->Get("ntSkim");
  tMCP->AddFriend(tMCPHI);
  tMCNP->AddFriend(tMCNPHI);
  tMCP->AddFriend(tMCPSkim);
  tMCNP->AddFriend(tMCNPSkim);
  
  TFile* outf = new TFile("bFeedDownPbPbMBMC.hist.root","recreate");
  
  TCut cutPbPbMB = "pclusterCompatibilityFilter&&pprimaryVertexFilter&&phfCoincFilter3&&Dy>-1.&&Dy<1.&&Dtrk1highPurity&&Dtrk2highPurity&&Dtrk1Pt>1.0&&Dtrk2Pt>1.0&&Dtrk1PtErr/Dtrk1Pt<0.3&&Dtrk2PtErr/Dtrk2Pt<0.3&&abs(Dtrk1Eta)<1.5&&abs(Dtrk2Eta)<1.5&&((DlxyBS/DlxyBSErr)>1.5&&Dalpha<0.12&&((Dpt>1&&Dpt<2&&(DsvpvDistance/DsvpvDisErr)>6.0&&Dchi2cl>0.25)||(Dpt>2&&Dpt<4&&(DsvpvDistance/DsvpvDisErr)>5.86&&Dchi2cl>0.224)||(Dpt>4&&Dpt<5&&(DsvpvDistance/DsvpvDisErr)>5.46&&Dchi2cl>0.196)||(Dpt>5&&Dpt<6&&(DsvpvDistance/DsvpvDisErr)>4.86&&Dchi2cl>0.170)||(Dpt>6&&Dpt<8&&(DsvpvDistance/DsvpvDisErr)>4.54&&Dchi2cl>0.125)||(Dpt>8&&Dpt<10&&(DsvpvDistance/DsvpvDisErr)>4.42&&Dchi2cl>0.091)||(Dpt>10&&Dpt<15&&(DsvpvDistance/DsvpvDisErr)>4.06&&Dchi2cl>0.069)||(Dpt>15&&Dpt<20&&(DsvpvDistance/DsvpvDisErr)>3.71&&Dchi2cl>0.056)||(Dpt>20&&Dpt<25&&(DsvpvDistance/DsvpvDisErr)>3.25&&Dchi2cl>0.054)||(Dpt>25&&(DsvpvDistance/DsvpvDisErr)>2.97&&Dchi2cl>0.050)))";

  TCut cutmc = "(Dgen==23333||Dgen==23344)";
  TCut cutmcSignal = "Dgen==23333";
  TCut cutmcSwapped = "Dgen==23344";
  TCut cutpt = "Dpt<20";
  TCut cutSignal = "abs(Dmass-1.8649)<0.025";
  TCut cutSideband = "abs(Dmass-1.8649)>0.075&&abs(Dmass-1.8649)<0.1";
  TCut cutPrompt = "DgenBAncestorpt<=0";
  TCut cutNonPrompt ="DgenBAncestorpt>0";
  TCut weightfunctionreco = "pthatweight";

  const int nBinX = 14;
  Float_t binsX[nBinX+1] = {2.,3.,4.,5.,6.,8.,10.,12.5,15.0,20.,25.,30.,40.,60.,100};

  const int nBinY = 20;
  Float_t binsY[nBinY+1];
  float firstBinYWidth = 0.001;
  float binYWidthRatio = 1.27;
  binsY[0]=0;
  for(int i=1; i<=nBinY; i++)
    binsY[i] = binsY[i-1]+firstBinYWidth*pow(binYWidthRatio,i-1);
  cout<<"last y bin: "<<binsY[nBinY]<<endl;

  const int nBinZ = 20;
  Float_t binsZ[nBinZ+1];
  float firstBinZ = 3.5;
  float binWidthZ = 5;
  for(int i=0; i<=nBinZ; i++)
    binsZ[i] = firstBinZ+binWidthZ*i;
  
  const int nBinM = 60;
  Float_t binsM[nBinM];
  float minMassBin = 1.7;
  float massBinWidth = 0.005;
  for(int i=0; i<=nBinM; i++)
    binsM[i] = minMassBin + massBinWidth*i;

  TH3D* hMCPSignal = new TH3D("hMCPSignal",";p_{T} (GeV/c);D0 DCA (cm);flight distance significance",nBinX,binsX,nBinY,binsY,nBinZ,binsZ);
  TH3D* hMCNPSignal = new TH3D("hMCNPSignal",";p_{T} (GeV/c);D0 DCA (cm);flight distance significance",nBinX,binsX,nBinY,binsY,nBinZ,binsZ);
  TH3D* hPtMD0DcaMCPSignal = new TH3D("hPtMD0DcaMCPSignal",";p_{T} (GeV/c);m (GeV/c^{2});D^{0} DCA (cm)",nBinX,binsX,nBinM,binsM,nBinY,binsY);
  TH3D* hPtMD0DcaMCPSwapped = new TH3D("hPtMD0DcaMCPSwapped",";p_{T} (GeV/c);m (GeV/c^{2});D^{0} DCA (cm)",nBinX,binsX,nBinM,binsM,nBinY,binsY);
  
  hMCPSignal->Sumw2();
  hMCNPSignal->Sumw2();
  hPtMD0DcaMCPSignal->Sumw2();
  hPtMD0DcaMCPSwapped->Sumw2();
  
  tMCP->Draw("DsvpvDistance/DsvpvDisErr:DsvpvDistance*sin(Dalpha):Dpt>>hMCPSignal",weightfunctionreco*(cutpt&&cutPbPbMB&&cutSignal&&cutmc&&cutPrompt));
  tMCNP->Draw("DsvpvDistance/DsvpvDisErr:DsvpvDistance*sin(Dalpha):Dpt>>hMCNPSignal",weightfunctionreco*(cutpt&&cutPbPbMB&&cutSignal&&cutmc&&cutNonPrompt));
  tMCP->Draw("DsvpvDistance*sin(Dalpha):Dmass:Dpt>>hPtMD0DcaMCPSignal",weightfunctionreco*(cutpt&&cutPbPbMB&&cutmcSignal&&cutPrompt));
  tMCP->Draw("DsvpvDistance*sin(Dalpha):Dmass:Dpt>>hPtMD0DcaMCPSwapped",weightfunctionreco*(cutpt&&cutPbPbMB&&cutmcSwapped&&cutPrompt));
  
  outf->Write();
}

void MomentumResSmear(){

	TFile* file1 = new TFile("/Users/kongkong/2014Research/ROOT_file/newEPOSsample/EPOS_TH3D_Nov21_vtx_2014.root");
	TFile* file2 = new TFile("/Users/kongkong/2014Research/ROOT_file/newEPOSsample/EPOS_TH3D_Nov21_vtx_pTsmear_2014.root");

	TH3D* genKS = (TH3D*)file1->Get("ana/genKS_underlying");
	TH3D* genLA = (TH3D*)file1->Get("ana/genLA_underlying");

	TH3D* genKS_smear = (TH3D*)file2->Get("ana/genKS_underlying");
	TH3D* genLA_smear = (TH3D*)file2->Get("ana/genLA_underlying");

	//cout << "number of entries: " << genKS->GetEntries() << endl;
	//cout << "number of entries: " << genKS_smear->GetEntries() << endl;


	TH1D* genks_mass[5][28];
	TH1D* genla_mass[5][20];

	TH1D* genks_smear_mass[5][28];
	TH1D* genla_smear_mass[5][28];



	double ks_pTbinsBound[29] = {0,1,2,3,4,5,6,7,8,9,10,12,14,16,18,20,22,24,26,28,30,34,38,42,46,50,56,66,90};
    double ks_ptbins[29] = {0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.2,1.4,1.6,1.8,2.0,2.2,2.4,2.6,2.8,3.0,3.4,3.8,4.2,4.6,5.0,5.6,6.6,9.0};
    double ks_binwidth[28] = {0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.4,0.4,0.4,0.4,0.4,0.6,1.0,2.4};
    double ks_ptbincenter[28] = {0.05,0.15,0.25,0.35,0.45,0.55,0.65,0.75,0.85,0.95,1.1,1.3,1.5,1.7,1.9,2.1,2.3,2.5,2.7,2.9,3.2,3.6,4.0,4.4,4.8,5.3,6.1,7.8};

    double la_pTbinsBound[21] = {6,8,10,12,14,16,18,20,22,24,26,28,30,34,38,42,46,50,56,66,90};
    double la_ptbins[21] = {0.6,0.8,1.0,1.2,1.4,1.6,1.8,2.0,2.2,2.4,2.6,2.8,3.0,3.4,3.8,4.2,4.6,5.0,5.6,6.6,9.0};
    double la_ptbincenter[20] = {0.7,0.9,1.1,1.3,1.5,1.7,1.9,2.1,2.3,2.5,2.7,2.9,3.2,3.6,4.0,4.4,4.8,5.3,6.1,7.8};
    double la_binwidth[20] = {0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.4,0.4,0.4,0.4,0.4,0.6,1.0,2.4};

    double rpybins[6] = {6,16,26,35,44,55};
    double rpybinwidth[5] = {1.0,0.97,0.9,0.93,1.0};

	stringstream ksName;
	stringstream laName;

	for(int rpy = 0; rpy < 5; rpy++){

		for(int pt = 0; pt < 28; pt++){

			ksName.str("");
			ksName << "genks_mass";
			ksName << "_";
			ksName << rpy+1;
			ksName << "_";
			ksName << pt+1;

			genks_mass[rpy][pt] = genKS->ProjectionZ(ksName.str().c_str(),rpybins[rpy]+1,rpybins[rpy+1],ks_pTbinsBound[pt]+1,ks_pTbinsBound[pt+1]);

			ksName.str("");
			ksName << "genks_smear_mass";
			ksName << "_";
			ksName << rpy+1;
			ksName << "_";
			ksName << pt+1;

			genks_smear_mass[rpy][pt] = genKS_smear->ProjectionZ(ksName.str().c_str(),rpybins[rpy]+1,rpybins[rpy+1],ks_pTbinsBound[pt]+1,ks_pTbinsBound[pt+1]);
		}

		for(pt = 0; pt < 20; pt++){

			laName.str("");
			laName << "genla_mass";
			laName << "_";
			laName << rpy+1;
			laName << "_";
			laName << pt+1;

			genla_mass[rpy][pt] = genLA->ProjectionZ(laName.str().c_str(),rpybins[rpy]+1,rpybins[rpy+1],la_pTbinsBound[pt]+1,la_pTbinsBound[pt+1]);

			laName.str("");
			laName << "genla_smear_mass";
			laName << "_";
			laName << rpy+1;
			laName << "_";
			laName << pt+1;

			genla_smear_mass[rpy][pt] = genLA_smear->ProjectionZ(laName.str().c_str(),rpybins[rpy]+1,rpybins[rpy+1],la_pTbinsBound[pt]+1,la_pTbinsBound[pt+1]);


		}
	}

   	int genksYield[5][28];
   	int genlaYield[5][20];

  	int genksSmearYield[5][28];
  	int genlaSmearYield[5][20];

   	for(rpy = 0; rpy < 5; rpy++){

   		for( pt = 0; pt < 28; pt++){

   			genksYield[rpy][pt] = genks_mass[rpy][pt]->GetEntries();
   			genksSmearYield[rpy][pt] = genks_smear_mass[rpy][pt]->GetEntries();
   		}
	   	
//.........这里部分代码省略.........

void correlate(
    const char* infname = "/data_CMS/cms/yilmaz/HiForest_HYDJET_Track8_Jet21_STARTHI53_LV1_merged_forest_0.root",
    const char* outname = "histograms_01.root",
    bool MC = 1,
    bool PbPb = 1,
    double jetEtaMax = 1.6,
    int centIndex = 0, int etaBin = 0, int leadJetPtBin = 0, int trackPtBin = 0
) {

    cout<<"Begin"<<endl;

    bool mini = 1;

    int Nevents = 1000;

    bool usePF = 1;
    int R = 3;

    double etaCOM = -0.465;

    bool doFlow = 0;
    bool doInclusiveJets = 1;
    bool doTracks = 1;
    bool doGenParticles = 1;

    bool fillTracks = 0;

    if(mini) {
        fillTracks = 0;
        doTracks = 0;
    }

    double pi = TMath::Pi();

    double trkMin = 0.5;

    double leadPtMin = 120;
    double subleadPtMin = 50;
    double dphiMin = 2.*pi/3.;

    int frame = 1; // Dijet frame for z

    int analysisId = 0;

    double ptSubLeadMin = 30;
    double ptLeadMin = 120;

    double etLeadMin[10] = {100, 100,120,150,180,200  };
    double etLeadMax[10] = {1000,120,150,180,200,1000 };

    double tkMin[10] = {4.,  4.,5.,7.,10., 20.};
    double tkMax[10] = {100.,5.,7.,10.,20.,1000.};

    TF1* gaus = new TF1("gaus","gaus",-5,5);
    gaus->SetParameter(0,1);
    gaus->SetParameter(1,0);
    gaus->SetParameter(2,1);

    TRandom* engin = new TRandom();

    //  double fitMin[10] = {100,80,60,50,100};
    //  double fitMax[10] = {250,200,150,100,250};

    double fitMin[10] = {100,90,80,70,100};
    double fitMax[10] = {300,300,300,300,300};

    TH1::SetDefaultSumw2();
    TH2::SetDefaultSumw2();
    cout<<"x"<<endl;

    string name[10] = {"c0to10","c10to20","c20to30","c30to50","c50to100","c0to30","c30to100","c0to100"};
    cout<<"x"<<endl;

    TFile* outf = new TFile(outname,"recreate");

    TNtuple *nt;
    nt = new TNtuple("nt","nt","x");

    TH3D* hAxisLead = new TH3D("hAxisLead","",500,0,500,160,-1.6,1.6,200,-pi,pi);
    TH3D* hAxisSubLead = new TH3D("hAxisSubLead","",500,0,500,160,-1.6,1.6,200,-pi,pi);

    TH3D* hCorrLead = new TH3D("hCorrLead","",500,0,500,160,-1.6,1.6,200,-pi,pi);
    TH3D* hCorrSubLead = new TH3D("hCorrSubLead","",500,0,500,160,-1.6,1.6,200,-pi,pi);

    TH3D* hGenParticleLead = new TH3D("hGenParticleLead","",500,0,500,160,-1.6,1.6,200,-pi,pi);
    TH3D* hGenParticleSubLead = new TH3D("hGenParticleSubLead","",500,0,500,160,-1.6,1.6,200,-pi,pi);

    TH1D* hPtLead = new TH1D("hPtLead","",500,0,500);
    TH1D* hPtSubLead = new TH1D("hPtSubLead","",500,0,500);

    bool pp = 0;

    double ajMin[10] = {0,    0.11, 0.22, 0.33, 0. , 0, 0};
    double ajMax[10] = {0.11, 0.22, 0.33, 1.,   1. , 0 ,0};

    int nEta = 4;
    double etaMin[10] = {0,  0,    0.5,  1.};
    double etaMax[10] = {5., 0.5,  1.,   2.};

    outf->cd();
//.........这里部分代码省略.........

void pileup()
{
//  TH3D* h = (TH3D*)GetHist("/net/hisrv0001/home/davidlw/scratch1/DiHadronCorrelations/outputs_312/HIData_Minbias_2760GeV/merged/HIData_Minbias_2760GeV_PPRereco_INCLMULT_nmin-1_nmax-1_etatrg-2.4-2.4_etaass-2.4-2.4_centmin-1_centmax-1.root","multmaxvssecvszvtxsep");
  TH3D* h = (TH3D*)GetHist("/net/hisrv0001/home/davidlw/scratch1/DiHadronCorrelations/outputs_312/HydjetReco_Minbias_2760GeV/merged/HydjetReco_Minbias_2760GeV_PeripheralNew_INCLMULT_nmin-1_nmax-1_etatrg-2.4-2.4_etaass-2.4-2.4_centmin-1_centmax-1.root","multmaxvssecvszvtxsep");
//  TH3D* h = (TH3D*)GetHist("/net/hisrv0001/home/davidlw/scratch1/DiHadronCorrelations/outputs_312/PAData_Minbias_5TeV/merged/PAData_Minbias_5TeV_HM_PromptReco_INCLMULTAllPURej3Vtx_nmin-1_nmax-1_etatrg-2.4-2.4_etaass-2.4-2.4_centmin-1_centmax-1.root","multmaxvssecvszvtxsep");

  TF1* func[40];
  TGraph* xx_yy[40];
  double par0_arr[40];
  double par1_arr[40];
  double par2_arr[40];
  double dz_arr[40];
//  double binwidth=0.2;
  double binedges[16]={0.0,0.2,0.3,0.4,0.5,0.6,0.8,1.0,1.4,1.8,2.4,3.0,4.2,5.0,7.0,10.0};

  TCanvas* c16 = new TCanvas("c16","c16",800,600);
  c16->Divide(4,4);

  for(int i=0;i<15;i++)
  {
//    h->GetZaxis()->SetRange(h->GetZaxis()->FindBin(binwidth*i),h->GetZaxis()->FindBin(binwidth*i+binwidth));
    h->GetZaxis()->SetRange(h->GetZaxis()->FindBin(binedges[i]),h->GetZaxis()->FindBin(binedges[i+1]));
    h->GetYaxis()->SetRange(h->GetYaxis()->FindBin(2),h->GetYaxis()->FindBin(50));
    TH2D* h2D = (TH2D*)h->Project3DProfile(Form("yx%d",i));
    h2D->Rebin2D(10,1);
    TProfile* hprof_x = (TProfile*)h2D->ProfileX("hprof_x",-1,-1,"s");
    double xx[500];
    double yy[500];
    for(int ibin=1;ibin<=hprof_x->GetNbinsX();ibin++)
    {
      yy[ibin-1] = hprof_x->GetBinContent(ibin)+2*hprof_x->GetBinError(ibin);
      xx[ibin-1] = hprof_x->GetBinCenter(ibin);
    }
  
// Hydjet
//    TF1* func = new TF1("func","[0]+[1]*(x-[2])*(x-[2])",0,500);
//    func->SetParameters(20,-0.1,150);
  func[i] = new TF1(Form("func_%d",i),"[0]*x**[1]",0,500);
  func[i]->SetParameters(1,0.5);
  func[i]->SetParLimits(0,0,10000);
  func[i]->SetParLimits(1,0,10000);
  if(i==0) func[i]->FixParameter(0,8); //pPb
//  else if(i>11) func[i]->FixParameter(0,20); //pPb
  else func[i]->FixParameter(0,2.); //pPb

//  func->FixParameter(0,5); //PbPb

//    TF1* func = new TF1("func","[0]+[1]*x",40,500);
//    double par1 = (hprof_x->GetBinContent(6)+hprof_x->GetBinError(6)-hprof_x->GetBinContent(1)-hprof_x->GetBinError(1))/(hprof_x->GetBinCenter(6)-hprof_x->GetBinCenter(1));
//    double par0 = hprof_x->GetBinContent(6)+hprof_x->GetBinError(6)-par1*hprof_x->GetBinCenter(6);
//    func->SetParameters(par0,par1);

    xx_yy[i] = new TGraph(hprof_x->GetNbinsX(),xx,yy);
    xx_yy[i]->Fit(Form("func_%d",i),"NO","",0,300);

    c16->cd(i+1);
    xx_yy[i]->Draw("AP"); 
    func[i]->SetLineColor(kRed);
    func[i]->Draw("Lsame");   

    par2_arr[i] = func[i]->GetParameter(2);    
    par1_arr[i] = func[i]->GetParameter(1);
    par0_arr[i] = func[i]->GetParameter(0);
    dz_arr[i] = (binedges[i]+binedges[i+1])/2;
  }

  TGraph* gr_par0 = new TGraph(15,dz_arr,par0_arr);
  TGraph* gr_par1 = new TGraph(15,dz_arr,par1_arr);
  TGraph* gr_par2 = new TGraph(15,dz_arr,par2_arr);

  TCanvas* cc = new TCanvas("cc","",950,400);  
  cc->Divide(3,1);
  cc->cd(1);
  gr_par0->Draw("AP");
  TF1* func_par0 = new TF1("func_par0","[2]*exp(-x**2/[0])*x**[3]+[1]",0.0,10.0);
  func_par0->SetParameters(0.2,0,5,-1);
  func_par0->FixParameter(2,0);
  gr_par0->Fit("func_par0","NO","",0.3,4.0);
  func_par0->Draw("Lsame");
  cc->cd(2);
  gr_par1->Draw("AP");
//  TF1* func_par1 = new TF1("func_par1","[2]*exp(-x**2/[0])*x**[3]+[1]",0.,10.0);
//  func_par1->SetParameters(0.2,0,-0.0005,0);
//  TF1* func_par1 = new TF1("func_par1","[2]*exp(-x/[0])*x**[3]+[1]",0.,10.0);
  TF1* func_par1 = new TF1("func_par1","[1]*exp(-x**[2]/[0])*x**[3]",0.,10.0);
  func_par1->SetParameters(0.1,0.6,2,1);
//  func_par1->FixParameter(2,1);
//  func_par1->FixParameter(3,1);

  gr_par1->Fit("func_par1","NO","",0.3,4.0);
  func_par1->Draw("Lsame");
  cc->cd(3);
  gr_par2->Draw("AP"); 
  TF1* func_par2 = new TF1("func_par2","[2]*exp(-x**2/[0])*x**[3]+[1]",0.,10.0);
  func_par2->SetParameters(1.2,170,120,1);
  gr_par2->Fit("func_par2","NO","",0.2,2.0);
  func_par2->Draw("Lsame");

cout<<Form("par0=%.2f*exp(-x**2/%.2f)*(x**%.2f)+%.2f",func_par0->GetParameter(2),func_par0->GetParameter(0),func_par0->GetParameter(3),func_par0->GetParameter(1))<<endl;
cout<<Form("par1=%.6f*exp(-x/%.6f)*(x**%.6f)+%.2f",func_par1->GetParameter(2),func_par1->GetParameter(0),func_par1->GetParameter(3),func_par1->GetParameter(1))<<endl;
//.........这里部分代码省略.........

void pileup1(double dzmin, double dzmax)
{
  TH1D* hdz = (TH1D*)GetHist("/net/hisrv0001/home/davidlw/scratch1/DiHadronCorrelations/outputs_312/PAData_Minbias_5TeV/merged/PAData_Minbias_5TeV_MBNEW_INCLMULT_nmin-1_nmax-1_etatrg-2.4-2.4_etaass-2.4-2.4_centmin-1_centmax-1.root","xyzvtxsepmin");
  TH3D* h1 = (TH3D*)GetHist("/net/hisrv0001/home/davidlw/scratch1/DiHadronCorrelations/outputs_312/PAData_Minbias_5TeV/merged/PAData_Minbias_5TeV_MBNEW_INCLMULT_nmin-1_nmax-1_etatrg-2.4-2.4_etaass-2.4-2.4_centmin-1_centmax-1.root","multmaxvssecvszvtxsep"); 
  TH2D* h21 = (TH2D*)GetHist("/net/hisrv0001/home/davidlw/scratch1/DiHadronCorrelations/outputs_312/PAData_Minbias_5TeV/merged/PAData_Minbias_5TeV_MBNEW_INCLMULT_nmin-1_nmax-1_etatrg-2.4-2.4_etaass-2.4-2.4_centmin-1_centmax-1.root","multmaxvssec");
  TH2D* hNvtxVsMult = (TH2D*)GetHist("/net/hisrv0001/home/davidlw/scratch1/DiHadronCorrelations/outputs_312/PAData_Minbias_5TeV/merged/PAData_Minbias_5TeV_MBNEW_INCLMULT_nmin-1_nmax-1_etatrg-2.4-2.4_etaass-2.4-2.4_centmin-1_centmax-1.root","nvtxvsnmult");

  TH1D* htrk_max1 = h1->ProjectionX();
  h1->GetZaxis()->SetRange(h1->GetZaxis()->FindBin(dzmin),h1->GetZaxis()->FindBin(dzmax));
  h1->GetYaxis()->SetRange(h1->GetYaxis()->FindBin(2),h1->GetYaxis()->FindBin(50));
  TH2D* h2D1 = (TH2D*)h1->Project3DProfile("yx100");
  h2D1->Rebin2D(5,1);
  TProfile* hprof1_x = (TProfile*)h2D1->ProfileX("hprof1_x",-1,-1,"s");
  TH1D* h1D1 = (TH1D*)h2D1->ProjectionY("h1D1",13,15,"e");
/*
  TH1D* hdz = (TH1D*)GetHist("/net/hisrv0001/home/davidlw/scratch1/DiHadronCorrelations/outputs_312/HIData_Minbias_2760GeV/merged/HIData_Minbias_2760GeV_PPRereco_INCLMULT_nmin-1_nmax-1_etatrg-2.4-2.4_etaass-2.4-2.4_centmin-1_centmax-1.root","xyzvtxsepmin");
  TH3D* h = (TH3D*)GetHist("/net/hisrv0001/home/davidlw/scratch1/DiHadronCorrelations/outputs_312/HIData_Minbias_2760GeV/merged/HIData_Minbias_2760GeV_PPRereco_INCLMULT_nmin-1_nmax-1_etatrg-2.4-2.4_etaass-2.4-2.4_centmin-1_centmax-1.root","multmaxvssecvszvtxsep");
  TH2D* h2 = (TH2D*)GetHist("/net/hisrv0001/home/davidlw/scratch1/DiHadronCorrelations/outputs_312/HIData_Minbias_2760GeV/merged/HIData_Minbias_2760GeV_PPRereco_INCLMULT_nmin-1_nmax-1_etatrg-2.4-2.4_etaass-2.4-2.4_centmin-1_centmax-1.root","multmaxvssec");
  TH2D* hNvtxVsMult = (TH2D*)GetHist("/net/hisrv0001/home/davidlw/scratch1/DiHadronCorrelations/outputs_312/HIData_Minbias_2760GeV/merged/HIData_Minbias_2760GeV_PPRereco_INCLMULT_nmin-1_nmax-1_etatrg-2.4-2.4_etaass-2.4-2.4_centmin-1_centmax-1.root","nvtxvsnmult");
*/

  TH1D* hdz = (TH1D*)GetHist("/net/hisrv0001/home/davidlw/scratch1/DiHadronCorrelations/outputs_312/HydjetReco_Minbias_2760GeV/merged/HydjetReco_Minbias_2760GeV_PeripheralNew_INCLMULT_nmin-1_nmax-1_etatrg-2.4-2.4_etaass-2.4-2.4_centmin-1_centmax-1.root","xyzvtxsepmin");
  TH3D* h = (TH3D*)GetHist("/net/hisrv0001/home/davidlw/scratch1/DiHadronCorrelations/outputs_312/HydjetReco_Minbias_2760GeV/merged/HydjetReco_Minbias_2760GeV_PeripheralNew_INCLMULT_nmin-1_nmax-1_etatrg-2.4-2.4_etaass-2.4-2.4_centmin-1_centmax-1.root","multmaxvssecvszvtxsep");
  TH2D* h2 = (TH2D*)GetHist("/net/hisrv0001/home/davidlw/scratch1/DiHadronCorrelations/outputs_312/HydjetReco_Minbias_2760GeV/merged/HydjetReco_Minbias_2760GeV_PeripheralNew_INCLMULT_nmin-1_nmax-1_etatrg-2.4-2.4_etaass-2.4-2.4_centmin-1_centmax-1.root","multmaxvssec");
  TH2D* hNvtxVsMult = (TH2D*)GetHist("/net/hisrv0001/home/davidlw/scratch1/DiHadronCorrelations/outputs_312/HydjetReco_Minbias_2760GeV/merged/HydjetReco_Minbias_2760GeV_PeripheralNew_INCLMULT_nmin-1_nmax-1_etatrg-2.4-2.4_etaass-2.4-2.4_centmin-1_centmax-1.root","nvtxvsnmult");

/*
  TH1D* hdz = (TH1D*)GetHist("/net/hisrv0001/home/davidlw/scratch1/DiHadronCorrelations/outputs_312/HijingReco_Minbias_2760GeV/merged/HijingReco_Minbias_2760GeV_53x_INCLMULT_nmin-1_nmax-1_etatrg-2.4-2.4_etaass-2.4-2.4_centmin-1_centmax-1.root","xyzvtxsepmin");
  TH3D* h = (TH3D*)GetHist("/net/hisrv0001/home/davidlw/scratch1/DiHadronCorrelations/outputs_312/HijingReco_Minbias_2760GeV/merged/HijingReco_Minbias_2760GeV_53x_INCLMULT_nmin-1_nmax-1_etatrg-2.4-2.4_etaass-2.4-2.4_centmin-1_centmax-1.root","multmaxvssecvszvtxsep");
  TH2D* hNvtxVsMult = (TH2D*)GetHist("/net/hisrv0001/home/davidlw/scratch1/DiHadronCorrelations/outputs_312/HijingReco_Minbias_2760GeV/merged/HijingReco_Minbias_2760GeV_53x_INCLMULT_nmin-1_nmax-1_etatrg-2.4-2.4_etaass-2.4-2.4_centmin-1_centmax-1.root","nvtxvsnmult");
*/
/*
  TH1D* hdz = (TH1D*)GetHist("/net/hisrv0001/home/davidlw/scratch1/DiHadronCorrelations/outputs_312/PAData_Minbias_5TeV/merged/PAData_Minbias_5TeV_MBNEW_INCLMULT_nmin-1_nmax-1_etatrg-2.4-2.4_etaass-2.4-2.4_centmin-1_centmax-1.root","xyzvtxsepmin");
  TH3D* h = (TH3D*)GetHist("/net/hisrv0001/home/davidlw/scratch1/DiHadronCorrelations/outputs_312/PAData_Minbias_5TeV/merged/PAData_Minbias_5TeV_MBNEW_INCLMULT_nmin-1_nmax-1_etatrg-2.4-2.4_etaass-2.4-2.4_centmin-1_centmax-1.root","multmaxvssecvszvtxsep"); 
  TH2D* h2 = (TH2D*)GetHist("/net/hisrv0001/home/davidlw/scratch1/DiHadronCorrelations/outputs_312/PAData_Minbias_5TeV/merged/PAData_Minbias_5TeV_MBNEW_INCLMULT_nmin-1_nmax-1_etatrg-2.4-2.4_etaass-2.4-2.4_centmin-1_centmax-1.root","multmaxvssec");
  TH2D* hNvtxVsMult = (TH2D*)GetHist("/net/hisrv0001/home/davidlw/scratch1/DiHadronCorrelations/outputs_312/PAData_Minbias_5TeV/merged/PAData_Minbias_5TeV_MBNEW_INCLMULT_nmin-1_nmax-1_etatrg-2.4-2.4_etaass-2.4-2.4_centmin-1_centmax-1.root","nvtxvsnmult");
*/
/*
  TH1D* hdz = (TH1D*)GetHist("/net/hisrv0001/home/davidlw/scratch1/DiHadronCorrelations/outputs_312/PAData_Minbias_5TeV/merged/PAData_Minbias_5TeV_HM_PromptSkim_INCLMULT220_nmin-1_nmax-1_etatrg-2.4-2.4_etaass-2.4-2.4_centmin-1_centmax-1.root","xyzvtxsepmin");
  TH3D* h = (TH3D*)GetHist("/net/hisrv0001/home/davidlw/scratch1/DiHadronCorrelations/outputs_312/PAData_Minbias_5TeV/merged/PAData_Minbias_5TeV_HM_PromptSkim_INCLMULT220_nmin-1_nmax-1_etatrg-2.4-2.4_etaass-2.4-2.4_centmin-1_centmax-1.root","multmaxvssecvszvtxsep"); 
  TH2D* h2 = (TH2D*)GetHist("/net/hisrv0001/home/davidlw/scratch1/DiHadronCorrelations/outputs_312/PAData_Minbias_5TeV/merged/PAData_Minbias_5TeV_HM_PromptSkim_INCLMULT220_nmin-1_nmax-1_etatrg-2.4-2.4_etaass-2.4-2.4_centmin-1_centmax-1.root","multmaxvssec"); 
  TH2D* hNvtxVsMult = (TH2D*)GetHist("/net/hisrv0001/home/davidlw/scratch1/DiHadronCorrelations/outputs_312/PAData_Minbias_5TeV/merged/PAData_Minbias_5TeV_HM_PromptSkim_INCLMULT220_nmin-1_nmax-1_etatrg-2.4-2.4_etaass-2.4-2.4_centmin-1_centmax-1.root","nvtxvsnmult");
*/
/*
  TH1D* hdz = (TH1D*)GetHist("/net/hisrv0001/home/davidlw/scratch1/DiHadronCorrelations/outputs_312/PAData_Minbias_5TeV/merged/PAData_Minbias_5TeV_HM_PromptRecoNew_INCLMULT220PURej2sigmaMode1_nmin-1_nmax-1_etatrg-2.4-2.4_etaass-2.4-2.4_centmin-1_centmax-1.root","xyzvtxsepmin");
  TH3D* h = (TH3D*)GetHist("/net/hisrv0001/home/davidlw/scratch1/DiHadronCorrelations/outputs_312/PAData_Minbias_5TeV/merged/PAData_Minbias_5TeV_HM_PromptRecoNew_INCLMULT220PURej2sigmaMode1_nmin-1_nmax-1_etatrg-2.4-2.4_etaass-2.4-2.4_centmin-1_centmax-1.root","multmaxvssecvszvtxsep"); 
  TH2D* h2 = (TH2D*)GetHist("/net/hisrv0001/home/davidlw/scratch1/DiHadronCorrelations/outputs_312/PAData_Minbias_5TeV/merged/PAData_Minbias_5TeV_HM_PromptRecoNew_INCLMULT220PURej2sigmaMode1_nmin-1_nmax-1_etatrg-2.4-2.4_etaass-2.4-2.4_centmin-1_centmax-1.root","multmaxvssec"); 
  TH2D* hNvtxVsMult = (TH2D*)GetHist("/net/hisrv0001/home/davidlw/scratch1/DiHadronCorrelations/outputs_312/PAData_Minbias_5TeV/merged/PAData_Minbias_5TeV_HM_PromptRecoNew_INCLMULT220PURej2sigmaMode1_nmin-1_nmax-1_etatrg-2.4-2.4_etaass-2.4-2.4_centmin-1_centmax-1.root","nvtxvsnmult");
*/
  TH1D* htrk_max = h->ProjectionX();
  h->GetZaxis()->SetRange(h->GetZaxis()->FindBin(dzmin),h->GetZaxis()->FindBin(dzmax));
  h->GetYaxis()->SetRange(h->GetYaxis()->FindBin(2),h->GetYaxis()->FindBin(241));
//  h->Rebin3D(10,1);
  TH2D* h2D = (TH2D*)h->Project3DProfile("yx");
  h2D->Rebin2D(5,1);
  TProfile* hprof_x = (TProfile*)h2D->ProfileX("hprof_x",-1,-1,"s");
//  TH1D* h1D = (TH1D*)h2D->ProjectionY("h1D",h2D->GetXaxis()->FindBin(200),h2D->GetXaxis()->FindBin(220),"e");
  TH1D* h1D = (TH1D*)h2->ProjectionY("h1D",h2->GetXaxis()->FindBin(200),h2->GetXaxis()->FindBin(220),"e");
   
  double xx[1000];
  double yy[1000];
  for(int i=1;i<=hprof_x->GetNbinsX();i++)
  {
    yy[i-1] = hprof_x->GetBinContent(i)+2*hprof_x->GetBinError(i);
    xx[i-1] = hprof_x->GetBinCenter(i);
  }

  TGraph* xx_yy = new TGraph(hprof_x->GetNbinsX(),xx,yy);

  int ibin = hprof_x->GetXaxis()->FindBin(200);
/*
  TF1* func = new TF1("func","[0]+[1]*x",0,500);
  double par1 = (hprof_x->GetBinContent(ibin)+hprof_x->GetBinError(ibin)-hprof_x->GetBinContent(1)-hprof_x->GetBinError(1))/(hprof_x->GetBinCenter(ibin)-hprof_x->GetBinCenter(1));
  double par0 = hprof_x->GetBinContent(ibin)+hprof_x->GetBinError(ibin)-par1*hprof_x->GetBinCenter(ibin);
  func->SetParameters(par0,par1);
*/
//  TF1* func = new TF1("func","[0]+[1]*(x-[2])*(x-[2])",0,500);
//  func->SetParameters(30,-0.1,150);
  TF1* func = new TF1("func","[0]*x**[1]",0,500);
  func->SetParameters(0.1,0.5);
  func->SetParLimits(0,0,10000);
  func->SetParLimits(1,0,10000);
//  func->FixParameter(0,5); small dz
  func->FixParameter(0,2);
  xx_yy->Fit("func","NO","",10,400);

  TCanvas* cc6 = new TCanvas("cc6","cc6",550,500);
  TH2D* htmp = new TH2D("htmp",";N_{trk}^{primary};N_{trk}^{secondary}",500,0,500,50,0,50);
  htmp->Draw("");
  xx_yy->Draw("Psame");
  func->SetLineColor(2);
  func->Draw("Lsame");
  TLatex* latex = new TLatex();
  latex->SetNDC();
  latex->SetTextSize(latex->GetTextSize()*0.75);
  latex->DrawLatex(0.23,0.87,"PbPb HYDJET, 1<|dz|<2 cm");
//  SaveCanvas(cc6,"pPb/corr","NpriVsNsecProfile_2Sigma_hydjet_dz1020_ForAN");

  TCanvas* cc4 = new TCanvas("cc4","cc4",550,500);
  cc4->SetLogz();
//.........这里部分代码省略.........

void GetYield()
{
	Int_t xmin, xmax, ymin, ymax, zmin, zmax;
	Double_t count;
	Char_t name[64];


	TH3D *h = (TH3D*)full.Get( "THits_TGG_ThetaCM_P");
	xmin = h->GetXaxis()->FindBin( 0);
	xmax = h->GetXaxis()->FindBin( 180);
	ymin = h->GetYaxis()->FindBin( 0);
	ymax = h->GetYaxis()->FindBin( 180);
	zmin = h->GetZaxis()->FindBin( 0);
	zmax = h->GetZaxis()->FindBin( 352);
	count = h->Integral(xmin, xmax, ymin, ymax, zmin, zmax);
	cout << count << endl;

	count = 0;
	TH3D *h1 = (TH3D*)full.Get( "THits_TGG_Phi_1perp_P");
	xmin = h1->GetXaxis()->FindBin( -180);
	xmax = h1->GetXaxis()->FindBin( 180);
	ymin = h1->GetYaxis()->FindBin( 0);
	ymax = h1->GetYaxis()->FindBin( 180);
	zmin = h1->GetZaxis()->FindBin( 0);
	zmax = h1->GetZaxis()->FindBin( 352);
	count += h1->Integral(xmin, xmax, ymin, ymax, zmin, zmax);
	cout << count << endl;
	TH3D *h2 = (TH3D*)full.Get( "THits_TGG_Phi_2perp_P");
	count += h2->Integral(xmin, xmax, ymin, ymax, zmin, zmax);
	cout << count << endl;
	TH3D *h3 = (TH3D*)full.Get( "THits_TGG_Phi_3perp_P");
	count += h3->Integral(xmin, xmax, ymin, ymax, zmin, zmax);
	cout << count << endl;
	TH3D *h4 = (TH3D*)full.Get( "THits_TGG_Phi_4perp_P");
	count += h4->Integral(xmin, xmax, ymin, ymax, zmin, zmax);
	cout << count << endl;
	TH3D *h5 = (TH3D*)full.Get( "THits_TGG_Phi_5perp_P");
	count += h5->Integral(xmin, xmax, ymin, ymax, zmin, zmax);
	cout << count << endl;
	TH3D *h6 = (TH3D*)full.Get( "THits_TGG_Phi_6perp_P");
	count += h6->Integral(xmin, xmax, ymin, ymax, zmin, zmax);
	cout << count << endl;
	TH3D *h7 = (TH3D*)full.Get( "THits_TGG_Phi_7perp_P");
	count += h7->Integral(xmin, xmax, ymin, ymax, zmin, zmax);
	cout << count << endl;
	TH3D *h8 = (TH3D*)full.Get( "THits_TGG_Phi_8perp_P");
	count += h8->Integral(xmin, xmax, ymin, ymax, zmin, zmax);
	cout << count << endl;
	TH3D *h9 = (TH3D*)full.Get( "THits_TGG_Phi_9perp_P");
	count += h9->Integral(xmin, xmax, ymin, ymax, zmin, zmax);
	cout << count << endl;

	TH3D *h1 = (TH3D*)full.Get( "THits_TGG_Phi_1para_P");
	count += h1->Integral(xmin, xmax, ymin, ymax, zmin, zmax);
	cout << count << endl;
	TH3D *h2 = (TH3D*)full.Get( "THits_TGG_Phi_2para_P");
	count += h2->Integral(xmin, xmax, ymin, ymax, zmin, zmax);
	cout << count << endl;
	TH3D *h3 = (TH3D*)full.Get( "THits_TGG_Phi_3para_P");
	count += h3->Integral(xmin, xmax, ymin, ymax, zmin, zmax);
	cout << count << endl;
	TH3D *h4 = (TH3D*)full.Get( "THits_TGG_Phi_4para_P");
	count += h4->Integral(xmin, xmax, ymin, ymax, zmin, zmax);
	cout << count << endl;
	TH3D *h5 = (TH3D*)full.Get( "THits_TGG_Phi_5para_P");
	count += h5->Integral(xmin, xmax, ymin, ymax, zmin, zmax);
	cout << count << endl;
	TH3D *h6 = (TH3D*)full.Get( "THits_TGG_Phi_6para_P");
	count += h6->Integral(xmin, xmax, ymin, ymax, zmin, zmax);
	cout << count << endl;
	TH3D *h7 = (TH3D*)full.Get( "THits_TGG_Phi_7para_P");
	count += h7->Integral(xmin, xmax, ymin, ymax, zmin, zmax);
	cout << count << endl;
	TH3D *h8 = (TH3D*)full.Get( "THits_TGG_Phi_8para_P");
	count += h8->Integral(xmin, xmax, ymin, ymax, zmin, zmax);
	cout << count << endl;
	TH3D *h9 = (TH3D*)full.Get( "THits_TGG_Phi_9para_P");
	count += h9->Integral(xmin, xmax, ymin, ymax, zmin, zmax);
	cout << count << endl;

}

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



  ////////////////////////
  //                    //
  // Declare Histograms //
  //                    //
  ////////////////////////
  //
  //
  //  For Using DownStream Scripts, please adhere to the conventions:
  //
  //
  //    histogram_name = "your_name_here"+"__"+regionList[i]+"__genClassy_"+genClassyObject.label+"__systematic_"+sysInfoObject.label;
  //
  //
  //  Where regionList is the list of "SR", "CR0b", "CR0b_tightBTagHighMlb" or "CR2l"
  //
  //    And systematicList[0] is the nominal selection
  //
  //    And if there is andditional selection involved in this histogram, please refer to it in "you name here"
  //

  cout << "    Preparing histograms" << endl << endl;


	//
	// Declare yield histograms
	//
  f_output->cd(); // All yield histos will belong to the output file
	for( analyzer* thisAnalyzer : analyzers ) {

		TH1D* h_template = thisAnalyzer->GetYieldTemplate();
		TH3D* h_template_signal = thisAnalyzer->GetYieldTemplateSignal();

		for( int iClassy=0; iClassy<nGenClassy; iClassy++ ) {
			for( int iSys=0; iSys<nSystematics; iSys++ ) {

				int histIndex = iClassy*nSystematics + iSys;

				// Gen and Systematic String
				TString reg_gen_sys_name = "__";
				reg_gen_sys_name += thisAnalyzer->GetLabel();
				reg_gen_sys_name += "__genClassy_";
				reg_gen_sys_name += genClassyList[iClassy].label;
				reg_gen_sys_name += "__systematic_";
				reg_gen_sys_name += systematicList[iSys].label;

				TH1* h_tmp = 0;

				TString yieldname = "h_yields";
				yieldname += reg_gen_sys_name;
				if( sample.isSignalScan ) h_tmp = (TH3D*)h_template_signal->Clone( yieldname );
				else                      h_tmp = (TH1D*)h_template->Clone( yieldname );
				thisAnalyzer->SetYieldHistogram( histIndex, h_tmp );

			} // End loop over systematics
		} // End loop over genClassys
	} // End loop over analyzers

  
  
  

  //
  // Declare non-yield histograms

void tight_loose_ratioplot()
{
  DEBUG("tight_loose_ratioplot() start");
  MakeCanvas(1,1);

  DEBUG("reading histograms");
  // "plot" histograms, i.e. read into memory
  cd(1);
  top();
  plot3("TightMuons");
  cd(2);
  top();
  plot3("LooseMuons");
  cd(1);

  TLegend * leg = new TLegend(0.53, 0.52, 0.77, 0.84);
  setopt(leg);
  Bool_t first = kTRUE;

  // get summed histograms
  vector<Int_t> entries;
  TH3D * hTightSum = 0;
  TH3D * hLooseSum = 0;
  for (Int_t i = 0; i < gMaxProcess; i++) {
    Int_t process = gOrder[0][i];
    // not existing
    TH3D * hTight3 = gHisto3[0][process];
    TH3D * hLoose3 = gHisto3[1][process];
    if (hTight3 == 0 || hLoose3 == 0) {
      ERROR("could not get histogram # " << i << " from memory");
      continue;
    }
    // if joined, it was already considered in previous iteration
    if (gProcess[process].join)
      continue;
    if (strncmp(gProcess[process].fname, "qcd", 3) &&
	strncmp(gProcess[process].fname, "dyll", 4) &&
	strncmp(gProcess[process].fname, "ttjets", 6) &&
	strncmp(gProcess[process].fname, "wjetstolnu", 10) &&
	strncmp(gProcess[process].fname, "data", 3))
      continue;
    DEBUG("Creating histos for process " << gProcess[process].fname);
    hTightSum = new TH3D(*hTight3);
    hLooseSum = new TH3D(*hLoose3);
    // check if histograms are joined -> we need to add statistics
    for (Int_t j = i+1; j < gMaxProcess; j++) {
      Int_t proc = gOrder[0][j];
      if (gOrder[1][j] != proc) {
	ERROR("wrong order - need to restart!");
	return;
      }
      // only add joined histograms
      if (!gProcess[proc].join)
	break;
      // not existing
      if (gHisto3[0][proc] == 0 || gHisto3[1][proc] == 0)
	continue;
      DEBUG("Adding histos for process " << gProcess[proc].fname);
      hTightSum->Add(gHisto3[0][proc], 1.);
      hLooseSum->Add(gHisto3[1][proc], 1.);
    }
    DEBUG("Getting ratio");
    TH1D * histo = get_1d_ratio(hTightSum, hLooseSum);
    if (histo == 0) {
      ERROR("Division failed");
      return;
    }
    check_1d_ratio_histogram(histo);
    setopt(histo);
    histo->SetLineColor(gProcess[process].lcolor);
    histo->SetLineStyle(gProcess[process].lstyle);
    DEBUG("Adding to legend");
    if (process != gMaxProcess-1) {
      leg->AddEntry(histo, gProcess[process].tname, "l");
    }
    else {
      leg->AddEntry(histo, gProcess[process].tname);
    }
    DEBUG("Draw");
    if (first) {
      histo->SetMaximum(1.);
      histo->SetMinimum(-0.1);
      histo->SetXTitle("p_{T}(#mu) [GeV]");
      histo->SetYTitle("T/L ratio");
      histo->Draw("ehisto");
      first = kFALSE;
    }
    else {
      if (process == gMaxProcess-1) {
	histo->SetMarkerStyle(gProcess[gMaxProcess-1].marker);
	histo->Draw("epsame");
      }
      else {
	histo->Draw("ehistosame");
      }
    }
  }
  // add subtracted data histogram
  TH1D * hdata_subtracted = get_subtracted_tight_loose_ratio(false, false);
  if (hdata_subtracted == 0)
//.........这里部分代码省略.........

void analyze(Int_t step){

  // TOF propagation factors (TOF efficiencies)
  fEfficiencyPiTOF = new TF1("fEfficiencyPiTOF","(x > 0.3)*0.7",0,10);
  fEfficiencyKaTOF = new TF1("fEfficiencyKaTOF","(x > 0.3)*(x-0.3)*(x<1) + (x>1)*0.7",0,10);
  fEfficiencyPrTOF = new TF1("fEfficiencyPrTOF","(x > 0.3)*(x-0.3)*(x<1) + (x>1)*0.7",0,10);

  // teoretical separation (perfect if equal to the one simualted in sim.C)
  fseparation = new TF1("f","[0]+[1]/x",0,100);
  fseparation->SetParameter(0,0.);
  fseparation->SetParameter(1,7.);
 
  fseparationPiKa = new TF1("fPiKa","[0]+[1]/TMath::Power(x,2.5)",0,100);
  fseparationPiKa->SetParameter(0,2.34);
  fseparationPiKa->SetParameter(1,10);

  fseparationKaPr = new TF1("fKaPr","[0]+[1]/TMath::Power(x,2.5)",0,100);
  fseparationKaPr->SetParameter(0,1);
  fseparationKaPr->SetParameter(1,56);

  // x=p, y=pt/p (normalized at the number of sigma assuming 80 ps resolution)
  fTOFpi = new TF2("fTOFpi","3.7/y*(sqrt(x*x+0.0193210)/x-1)*37.47405725",0.3,10,0.5,1);
  fTOFka = new TF2("fTOFka","3.7/y*(sqrt(x*x+0.243049)/x-1)*37.47405725",0.3,10,0.5,1);
  fTOFpr = new TF2("fTOFpr","3.7/y*(sqrt(x*x+0.879844)/x-1)*37.47405725",0.3,10,0.5,1);


  // x=p, already normalized in number of sigma (sigma assumed 3.5=7% of the MIP)
  fTPCpi = new TF1("fTPCpi",BetheBlochAleph,0,10,6);
  fTPCpi->SetParameter(0,fKp1);
  fTPCpi->SetParameter(1,fKp2);
  fTPCpi->SetParameter(2,fKp3);
  fTPCpi->SetParameter(3,fKp4);
  fTPCpi->SetParameter(4,fKp5);
  fTPCpi->SetParameter(5,0.139);
  fTPCka = new TF1("fTPCka",BetheBlochAleph,0,10,6);
  fTPCka->SetParameter(0,fKp1);
  fTPCka->SetParameter(1,fKp2);
  fTPCka->SetParameter(2,fKp3);
  fTPCka->SetParameter(3,fKp4);
  fTPCka->SetParameter(4,fKp5);
  fTPCka->SetParameter(5,0.493);
  fTPCpr = new TF1("fTPCpr",BetheBlochAleph,0,10,6);
  fTPCpr->SetParameter(0,fKp1);
  fTPCpr->SetParameter(1,fKp2);
  fTPCpr->SetParameter(2,fKp3);
  fTPCpr->SetParameter(3,fKp4);
  fTPCpr->SetParameter(4,fKp5);
  fTPCpr->SetParameter(5,0.938);

  Float_t width = 1.0;
  addshift =0;

  invwidth = 1./width;

  Float_t widthTOF = 1.0;
  addshiftTOF =0;

  invwidthTOF = 1./widthTOF;

  TH1D *priorsPt[6];
  TH1D *newpriorsPt[6];
  TH1D *truePt[6];
  TH1D *allPtPos = new TH1D("allPtP","All positive;p_{T} (GeV/#it{c});N",100,0,10);
  TH1D *allPtNeg = new TH1D("allPtN","All negative;p_{T} (GeV/#it{c});N",100,0,10);

  TH3D *priorsKs[3][3];
  TH3D *newpriorsKs[3][3];
  TH3D *truePidKs[3][3];
  TH3D *trueKs;

  TH2D *priorsPhi[3][3];
  TH2D *newpriorsPhi[3][3];
  TH2D *truePidPhi[3][3];
  TH2D *truePhi;

  TH3D *priorsLc[3][3][3];
  TH3D *newpriorsLc[3][3][3];
  TH3D *truePidLc[3][3][3];
  TH3D *trueLc,*mypidLc;

  TH3D *priorsLcbar[3][3][3];
  TH3D *newpriorsLcbar[3][3][3];
  TH3D *truePidLcbar[3][3][3];
  TH3D *trueLcbar,*mypidLcbar;

  Int_t nbinPtFrKa = 8;
  Int_t nbinPtFrPi = 8;
  Int_t nbinY = 1;
  Int_t nbinpol=nbinPtFrKa*nbinPtFrPi*nbinY;
  Double_t normbin = 1./nbinpol;
  Int_t nbinmlc = 100;
  Int_t nbinptlc = 10;

  const char *spec[3] = {"Pi","Ka","Pr"};

  if(step==0){
    priorsPt[0] = new TH1D("oldpriorsPtPiP","Pion (+) priors;p_{T} (GeV/#it{c});N",100,0,10);
    for(Int_t i=1;i<=100;i++)
      priorsPt[0]->SetBinContent(i,1);
    priorsPt[1] = new TH1D("oldpriorsPtKaP","Kaon (+) priors;p_{T} (GeV/#it{c});N",100,0,10);
//.........这里部分代码省略.........