/* Float_t quadFun(Double_t *x, Double_t *par){ if (x[0] > par[3] && x[0] < par[4]){ TF1::RejectPoint(); return 0; } return par[0] + par[1]*x[0] + par[2]*x[0]*x[0] ; //a+bx+cx**2 } */ void processResultsAllTrigPt(Int_t specie = 0, Int_t centMin = 0, Int_t centMax = 10, Bool_t bAnalyticCorr= 1 /*kTRUE*/, Bool_t bMix = 0, //Bool_t bMerged = 1, Bool_t bAll = 1){ //Note: not merged for the centrality bin of 10-40% for some cases gROOT->LoadMacro("fitHisto.C"); if(!bAnalyticCorr){ if (!bMix) { cout<<"Can not proceed further if the option is non-mixed!"<SetParameters(1., 0.625); TString specieName[] = {"K0s", "Lambda", "AntiLambda"}; // defining bins for Z vertex const Int_t nVtx = 8; Double_t vertexBins[] = {-8., -6., -4., -2., 0., 2., 4., 6., 8.}; // pt bins of associated particles for the analysis const Int_t nPtAsso = 4; const Double_t PtBins[5] = {1.0, 2.0, 3.0, 4.0, 8.0}; // const Int_t nPtAsso = 1; //const Double_t PtBins[2] = {1.0, 2.0}; // pt bins of trigger particles for the analysis const Int_t nPtTrig = 5;//4 const Double_t PtBinsV0Xi[5] = {1.0, 2.0, 4.0, 8.0, 15.0}; // Int_t rebin = 2; const Int_t nPhi = 20; const Int_t nEta = 32; // defining bins for centrality //const Int_t nCent = 2; //Double_t centBins[] = {0.,10.,20.,30.,40.,50.,60., 90.}; Double_t massMean = 0.493677; Double_t massRes = 0.005; //if(specie == 1 || specie == 2){ // massMean = 1.115683; // massRes = 0.0015; //} Double_t fitRangeLow = massMean-15.*massRes; Double_t fitRangeHigh = massMean+15.*massRes; //TH1F * hNtrig; TH3F * h3EtaPhiSB[nVtx][nPtTrig][nPtAsso]; TH2F * h2EtaPhiSB[nVtx][nPtTrig][nPtAsso]; TH2F * h2EtaPhiSBMix[nVtx][nPtTrig][nPtAsso]; TH3F * h3EtaPhiSBCor[nVtx][nPtTrig][nPtAsso]; //corrected by acceptance TH3F * h3EtaPhiSBCorAll[nPtTrig][nPtAsso] = {{NULL}}; // sum up vtx bins TH2F * h2EtaPhiSBPeakCor[nVtx][nPtTrig][nPtAsso]; TH2F * h2EtaPhiSBPeakCorAll[nPtTrig][nPtAsso] = {{NULL}}; TH2F * h2EtaPhiS[nVtx][nPtTrig][nPtAsso]; //after background subtraction TH2F * h2EtaPhiSPerTrig[nPtTrig][nPtAsso] = {{NULL}}; // per trigger phi correlation TH1F * hMass[nVtx][nPtTrig][nPtAsso]; THnSparseF *fHistTrigSibAllK0s; THnSparseF *fHistdPhidEtaSibK0s[nVtx]; THnSparseF *fHistdPhidEtaMixK0s[nVtx]; TFile * file = new TFile("ResultData18JUL.root", "READ"); /* TFile * file; if(bMerged){ if(bAll) file = new TFile("ResultData18JUL.root", "READ"); else file = new TFile(Form("MergedResultsCent%d_%d.root", centMin, centMax), "READ"); */ //gDirectory->cd(Form("Cent%d_%d", centMin, centMax)); TIter nextkey(gDirectory->GetListOfKeys()); TKey* key; while (Key = (TKey*)nextkey()){ if(!strcmp( key->GetName(), Form("OutputCent%d_%d", centMin, centMax))){ TList * V0Xi = (TList*)key->ReadObj(); fHistTrigSibAllK0s = (THnSparseF *)V0Xi->FindObject("K0s")->FindObject("fHistTrigSibAllK0s"); if(!fHistTrigSibAllK0s) cout<<"Does not work!"<FindObject("K0s")->FindObject("fHistdPhidEtaSibK0s"); if(!fHistdPhidEtaSibK0s) cout<<"Does not work!"<FindObject("K0s")->FindObject("fHistdPhidEtaMixK0s"); if(!fHistdPhidEtaMixK0s) cout<<"Does not work!"<FindObject("K0s"))->FindObject("fHistTrigSibAllK0s"); if(!fHistTrigSibAllK0s) cout<<"Does not work!"<FindObject(Form("K0s")))->FindObject(Form("fHistdPhidEtaSibK0s")); // // fHistdPhidEtaSibK0s[i]== (THnSparseF*)(V0Xi->FindObject("AnalysisK0s"))->FindObject(Form("fHistdPhidEtaSibK0sVtx%d",i)); if(!fHistdPhidEtaSibK0s[i]) cout<<"Sib Does not work!"<FindObject(Form("K0s"))->FindObject(Form("fHistdPhidEtaMixK0s")); // fHistdPhidEtaMixK0s[i]== (THnSparseF*)V0Xi->FindObject("AnalysisK0s")->FindObject(Form("fHistdPhidEtaMixK0sVtx%d",i)); if(!fHistdPhidEtaMixK0s[i]) cout<<"Mix Does not work!"<>fileName){ file = new TFile(fileName, "READ"); TIter nextkey(gDirectory->GetListOfKeys()); TKey *key; while (key = (TKey*)nextkey()){ if( !strcmp( key->GetName(), Form("OutputCent%d_%d", centMin, centMax))){ TList * V0Xi = (TList*)key->ReadObj(); if(!fHistTrigSibAllK0s) fHistTrigSibAllK0sSB *= (TH2F*)V0Xi->FindObject("fHistTrigSibAllK0s"); else{ TH2F * tmp = (TH2F*)V0Xi->FindObject("fHistTrigSibAllK0s"); fHistTrigSibAllK0s->Add(tmp); } if(!fHistTrigSibAllK0s) cout<<"Does not work!"<FindObject(Form("fHistdPhidEtaSib%sVtx%d", specieName[specie].Data(), i)); else{ THnSparseF * tmpHn = (THnSparseF*)V0Xi->FindObject(Form("fHistdPhidEtaSib%sVtx%d", specieName[specie].Data(), i)); fHistdPhidEtaSibK0s[i]->Add(tmpHn); } if(!fHistdPhidEtaSibK0s[i]) cout<<"Does not work!"<FindObject(Form("fHistdPhidEtaMix%sVtx%d", specieName[specie].Data(), i)); else{ THnSparseF * tmpHnMix = (THnSparseF*)V0Xi->FindObject(Form("fHistdPhidEtaMix%sVtx%d", specieName[specie].Data(), i)); fHistdPhidEtaMixK0s[i]->Add(tmpHnMix); } if(!fHistdPhidEtaMixK0s[i]) cout<<"Does not work!"<ProjectionX(); hNtrig->SetName(Form("hPtTrigCent%d_%d", centMin, centMax)); hNtrig->SetTitle(Form("Cent: [%d, %d]; p_{T}^{Trig} [GeV/c]", centMin, centMax)); */ for(Int_t iVtx = 0; iVtx < nVtx; iVtx++) for(Int_t iVtx = 0; iVtx < nVtx; iVtx++) for(Int_t iPtT = 0; iPtT < nPtTrig; iPtT++) for(Int_t iPtA = 0; iPtA < nPtAsso; iPtA++){ if(PtBins[iPtA] >= PtBinsV0Xi[iPtT+1]) break; Int_t ptTrigBinLow = fHistTrigSibAllK0s[iVtx]->GetAxis(0)->FindBin(PtBinsV0Xi[iPtT]+0.01); Int_t ptTrigBinHigh = fHistTrigSibAllK0s[iVtx]->GetAxis(0)->FindBin(PtBinsV0Xi[iPtT+1]-0.01); //changwan change /* Int_t ptAssoBinLow = fHistdPhidEtaSibK0s[iVtx]->GetAxis(2)->FindBin(PtBins[iPtA]+0.01); Int_t ptAssoBinHigh = fHistdPhidEtaSibK0s[iVtx]->GetAxis(2)->FindBin(PtBins[iPtA+1]-0.01); */ fHistTrigSibAllK0s[iVtx]->GetAxis(0)->SetRange(ptTrigBinLow, ptTrigBinHigh); // fHistdPhidEtaSibK0s[iVtx]->GetAxis(2)->SetRange(ptAssoBinLow, ptAssoBinHigh); hMass[iVtx][iPtT][iPtA] = (TH1F*)fHistTrigSibAllK0s[iVtx]->Projection(0, "OE"); hMass[iVtx][iPtT][iPtA] ->SetName(Form("MassCent%d_%dVtx%.1f_%.1fPtTr%.1f_%.1fAs%.1f", centMin, centMax, vertexBins[iVtx], vertexBins[iVtx+1], PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1], PtBins[iPtA])); hMass[iVtx][iPtT][iPtA] ->SetTitle(Form("Mass, Cent:[%d, %d], Vtx:[%.1f, %.1f], p_{T}^{Tr}:[%.1f, %.1f], p_{T}^{As}:[%.1f, %.1f]; mass [GeV/c^{2}]", centMin, centMax, vertexBins[iVtx], vertexBins[iVtx+1], PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1], PtBins[iPtA], PtBins[iPtA+1])); } for(Int_t iVtx = 0; iVtx < nVtx; iVtx++) for(Int_t iPtT = 0; iPtT < nPtTrig; iPtT++) for(Int_t iPtA = 0; iPtA < nPtAsso; iPtA++){ if(PtBins[iPtA] >= PtBinsV0Xi[iPtT+1]) break; Int_t ptTrigBinLow = fHistdPhidEtaSibK0s[iVtx]->GetAxis(2)->FindBin(PtBinsV0Xi[iPtT]+0.01); Int_t ptTrigBinHigh = fHistdPhidEtaSibK0s[iVtx]->GetAxis(2)->FindBin(PtBinsV0Xi[iPtT+1]-0.01); //changwan change Int_t ptAssoBinLow = fHistdPhidEtaSibK0s[iVtx]->GetAxis(2)->FindBin(PtBins[iPtA]+0.01); Int_t ptAssoBinHigh = fHistdPhidEtaSibK0s[iVtx]->GetAxis(2)->FindBin(PtBins[iPtA+1]-0.01); fHistdPhidEtaSibK0s[iVtx]->GetAxis(2)->SetRange(ptTrigBinLow, ptTrigBinHigh);//changwan change // fHistdPhidEtaSibK0s[iVtx]->GetAxis(2)->SetRange(ptAssoBinLow, ptAssoBinHigh); // Int_t binLow = fHistdPhidEtaSibK0s[iVtx]->GetAxis(4)->FindBin(massMean-3.*massRes); //Int_t binHigh = fHistdPhidEtaSibK0s[iVtx]->GetAxis(4)->FindBin(massMean+3.*massRes); //fHistdPhidEtaSibK0s[iVtx]->GetAxis(4)->SetRange(binLow, binHigh); h3EtaPhiSB[iVtx][iPtT][iPtA] = (TH3F*)fHistdPhidEtaSibK0s[iVtx]->Projection(5, 2, 3, "OE"); //changwan change h3EtaPhiSB[iVtx][iPtT][iPtA] ->SetName(Form("EtaPhiSBCent%d_%dVtx%.1f_%.1fPtTr%.1f_%.1fAs%.1f", centMin, centMax, vertexBins[iVtx], vertexBins[iVtx+1], PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1], PtBins[iPtA])); h3EtaPhiSB[iVtx][iPtT][iPtA] ->SetTitle(Form("S. R., Cent:[%d, %d], Vtx:[%.1f, %.1f], p_{T}^{Tr}:[%.1f, %.1f], p_{T}^{As}:[%.1f, %.1f];#Delta #phi; #Delta#eta; mass", centMin, centMax, vertexBins[iVtx], vertexBins[iVtx+1], PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1], PtBins[iPtA], PtBins[iPtA+1])); //h3EtaPhiSB[iVtx][iPtT][iPtA]->RebinX(rebin); } for(Int_t iVtx = 0; iVtx < nVtx; iVtx++) for(Int_t iPtT = 0; iPtT < nPtTrig; iPtT++) for(Int_t iPtA = 0; iPtA < nPtAsso; iPtA++){ if(PtBins[iPtA] >= PtBinsV0Xi[iPtT+1]) break; Int_t binLow = h3EtaPhiSB[iVtx][iPtT][iPtA]->GetZaxis()->FindBin(massMean-3.*massRes); Int_t binHigh = h3EtaPhiSB[iVtx][iPtT][iPtA]->GetZaxis()->FindBin(massMean+3.*massRes); h3EtaPhiSB[iVtx][iPtT][iPtA]->GetZaxis()->SetRange(binLow, binHigh); h2EtaPhiSB[iVtx][iPtT][iPtA] = (TH2F*)h3EtaPhiSB[iVtx][iPtT][iPtA]->Project3D("xy"); h2EtaPhiSB[iVtx][iPtT][iPtA] ->SetName(Form("h2EtaPhiSBCent%d_%dVtx%.1f_%.1fPtTr%.1f_%.1fAs%.1f", centMin, centMax, vertexBins[iVtx], vertexBins[iVtx+1], PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1], PtBins[iPtA])); h2EtaPhiSB[iVtx][iPtT][iPtA] ->SetTitle(Form("S. R., Cent:[%d, %d], Vtx:[%.1f, %.1f], p_{T}^{Tr}:[%.1f, %.1f], p_{T}^{As}:[%.1f, %.1f];#Delta#eta; #Delta#phi", centMin, centMax, vertexBins[iVtx], vertexBins[iVtx+1], PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1], PtBins[iPtA], PtBins[iPtA+1])); //} if(bMix) for(Int_t iVtx = 0; iVtx < nVtx; iVtx++) for(Int_t iPtT = 0; iPtT < nPtTrig; iPtT++) for(Int_t iPtA = 0; iPtA < nPtAsso; iPtA++){ if(PtBins[iPtA] >= PtBinsV0Xi[iPtT+1]) break; Int_t ptTrigBinLow = fHistdPhidEtaMixK0s[iVtx]->GetAxis(2)->FindBin(PtBinsV0Xi[iPtT]+0.01); Int_t ptTrigBinHigh = fHistdPhidEtaMixK0s[iVtx]->GetAxis(2)->FindBin(PtBinsV0Xi[iPtT+1]-0.01); //changwan change //Int_t ptAssoBinLow = fHistdPhidEtaSibK0s[iVtx]->GetAxis(2)->FindBin(PtBins[iPtA]+0.01); // Int_t ptAssoBinHigh = fHistdPhidEtaSibK0s[iVtx]->GetAxis(2)->FindBin(PtBins[iPtA+1]-0.01); fHistdPhidEtaMixK0s[iVtx]->GetAxis(2)->SetRange(ptTrigBinLow, ptTrigBinHigh);//changwan change // fHistdPhidEtaMixK0s[iVtx]->GetAxis(2)->SetRange(ptAssoBinLow, ptAssoBinHigh); Int_t binLow = fHistdPhidEtaSibK0s[iVtx]->GetAxis(2)->FindBin(massMean-3.*massRes); Int_t binHigh = fHistdPhidEtaSibK0s[iVtx]->GetAxis(2)->FindBin(massMean+3.*massRes); fHistdPhidEtaMixK0s[iVtx]->GetAxis(2)->SetRange(binLow, binHigh); h2EtaPhiSBMix[iVtx][iPtT][iPtA] = (TH2F*)fHistdPhidEtaMixK0s[iVtx]->Projection(5, 2, 3, "OE"); h2EtaPhiSBMix[iVtx][iPtT][iPtA] ->SetName(Form("EtaPhiSBMixCent%d_%dVtx%.1f_%.1fPtTr%.1f_%.1fAs%.1f", centMin, centMax, vertexBins[iVtx], vertexBins[iVtx+1], PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1], PtBins[iPtA])); h2EtaPhiSBMix[iVtx][iPtT][iPtA] ->SetTitle(Form("Cent:[%d, %d], Vtx:[%.1f, %.1f], p_{T}^{Tr}:[%.1f, %.1f], p_{T}^{As}:[%.1f, %.1f];#Delta#eta; #Delta#phi", centMin, centMax, vertexBins[iVtx], vertexBins[iVtx+1], PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1], PtBins[iPtA], PtBins[iPtA+1])); //h2EtaPhiSBMix[iVtx][iPtT][iPtA]->RebinY(rebin); } //corrected for the acceptance effect for(Int_t iVtx = 0; iVtx < nVtx; iVtx++) for(Int_t iPtT = 0; iPtT < nPtTrig; iPtT++) for(Int_t iPtA = 0; iPtA < nPtAsso; iPtA++){ if(PtBins[iPtA] >= PtBinsV0Xi[iPtT+1]) break; h3EtaPhiSBCor[iVtx][iPtT][iPtA] = (TH3F*) h3EtaPhiSB[iVtx][iPtT][iPtA]->Clone(); h3EtaPhiSBCor[iVtx][iPtT][iPtA] ->SetName(Form("EtaPhiSBCorCent%d_%dVtx%.1f_%.1fPtTr%.1f_%.1fAs%.1f", centMin, centMax, vertexBins[iVtx], vertexBins[iVtx+1], PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1], PtBins[iPtA])); h3EtaPhiSBCor[iVtx][iPtT][iPtA] ->SetTitle(Form("Cent:[%d, %d], Vtx:[%.1f, %.1f], p_{T}^{Tr}:[%.1f, %.1f],p_{T}^{As}:[%.1f, %.1f];#Delta#phi; #Delta#eta; mass", centMin, centMax, vertexBins[iVtx], vertexBins[iVtx+1], PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1], PtBins[iPtA], PtBins[iPtA+1])); if(h3EtaPhiSBCor[iVtx][iPtT][iPtA]){ for(Int_t iPhi = 0; iPhi < h3EtaPhiSBCor[iVtx][iPtT][iPtA]->GetXaxis()->GetNbins(); iPhi++) for(Int_t iEta = 0; iEta < h3EtaPhiSBCor[iVtx][iPtT][iPtA]->GetYaxis()->GetNbins(); iEta++){ //note 2D projection x,y axes are in different order than 3d Double_t corr = 0; if(!bAnalyticCorr && bMix) corr = h2EtaPhiSBMix[iVtx][iPtT][iPtA]->GetBinContent(iEta+1, iPhi+1); else corr = fAccCorr->Eval(h3EtaPhiSBCor[iVtx][iPtT][iPtA]->GetYaxis()->GetBinCenter(iEta+1)); if(corr == 0) continue; for(Int_t iMass = 0; iMass < h3EtaPhiSBCor[iVtx][iPtT][iPtA]->GetZaxis()->GetNbins(); iMass++){ h3EtaPhiSBCor[iVtx][iPtT][iPtA]->SetBinContent(iPhi+1, iEta+1, iMass+1, h3EtaPhiSBCor[iVtx][iPtT][iPtA]->GetBinContent(iPhi+1, iEta+1, iMass+1) /corr); h3EtaPhiSBCor[iVtx][iPtT][iPtA]->SetBinError(iPhi+1, iEta+1, iMass+1, h3EtaPhiSBCor[iVtx][iPtT][iPtA]->GetBinError(iPhi+1, iEta+1, iMass+1) /corr); } } if(!bAnalyticCorr && bMix) h3EtaPhiSBCor[iVtx][iPtT][iPtA] ->Scale(h2EtaPhiSBMix[iVtx][iPtT][iPtA]->GetMaximum()); } h2EtaPhiSBPeakCor[iVtx][iPtT][iPtA] = (TH2F*) h2EtaPhiSB[iVtx][iPtT][iPtA]->Clone(); h2EtaPhiSBPeakCor[iVtx][iPtT][iPtA] ->SetName(Form("h2EtaPhiSBPeakCorCent%d_%dVtx%.1f_%.1fPtTr%.1f_%.1fAs%.1f", centMin, centMax, vertexBins[iVtx], vertexBins[iVtx+1], PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1], PtBins[iPtA])); h2EtaPhiSBPeakCor[iVtx][iPtT][iPtA] ->SetTitle(Form("Cent:[%d, %d], Vtx:[%.1f, %.1f], p_{T}^{Tr}:[%.1f, %.1f],p_{T}^{As}:[%.1f, %.1f];#Delta#eta; #Delta#phi", centMin, centMax, vertexBins[iVtx], vertexBins[iVtx+1], PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1], PtBins[iPtA], PtBins[iPtA+1])); if(h2EtaPhiSBPeakCor[iVtx][iPtT][iPtA]){ for(Int_t iPhi = 0; iPhi < h3EtaPhiSBCor[iVtx][iPtT][iPtA]->GetXaxis()->GetNbins(); iPhi++) for(Int_t iEta = 0; iEta < h3EtaPhiSBCor[iVtx][iPtT][iPtA]->GetYaxis()->GetNbins(); iEta++){ Double_t corr = 0; if(!bAnalyticCorr && bMix) corr = h2EtaPhiSBMix[iVtx][iPtT][iPtA]->GetBinContent(iEta+1, iPhi+1); else corr = fAccCorr->Eval(h3EtaPhiSBCor[iVtx][iPtT][iPtA]->GetYaxis()->GetBinCenter(iEta+1)); if(corr == 0) continue; h2EtaPhiSBPeakCor[iVtx][iPtT][iPtA]->SetBinContent(iEta+1, iPhi+1, h2EtaPhiSBPeakCor[iVtx][iPtT][iPtA]->GetBinContent(iEta+1, iPhi+1) /corr); h2EtaPhiSBPeakCor[iVtx][iPtT][iPtA]->SetBinError(iEta+1, iPhi+1, h2EtaPhiSBPeakCor[iVtx][iPtT][iPtA]->GetBinError(iEta+1, iPhi+1) /corr); } if(!bAnalyticCorr && bMix) h2EtaPhiSBPeakCor[iVtx][iPtT][iPtA] ->Scale(h2EtaPhiSBMix[iVtx][iPtT][iPtA]->GetMaximum()); } } //sum up vtx bins for(Int_t iPtT = 0; iPtT < nPtTrig; iPtT++) for(Int_t iPtA = 0; iPtA < nPtAsso; iPtA++){ if(PtBins[iPtA] >= PtBinsV0Xi[iPtT+1]) break; for(Int_t iVtx = 0; iVtx < nVtx; iVtx++){ if(iVtx == 0 && !h3EtaPhiSBCorAll[iPtT][iPtA]){ h3EtaPhiSBCorAll[iPtT][iPtA] = (TH3F*)h3EtaPhiSBCor[iVtx][iPtT][iPtA]->Clone(); h3EtaPhiSBCorAll[iPtT][iPtA]->SetName(Form("EtaPhiSBCorAllCent%d_%dPtTr%.1f_%.1fAs%.1f", centMin, centMax, PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1], PtBins[iPtA])); h3EtaPhiSBCorAll[iPtT][iPtA] ->SetTitle(Form("Cent:[%d, %d], p_{T}^{Tr}:[%.1f, %.1f],p_{T}^{As}:[%.1f, %.1f];#Delta#phi; #Delta#eta; mass", centMin, centMax, PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1], PtBins[iPtA], PtBins[iPtA+1])); h2EtaPhiSBPeakCorAll[iPtT][iPtA] = (TH2F*)h2EtaPhiSBPeakCor[iVtx][iPtT][iPtA]->Clone(); h2EtaPhiSBPeakCorAll[iPtT][iPtA]->SetName(Form("h2EtaPhiSBPeakCorAllCent%d_%dPtTr%.1f_%.1fAs%.1f", centMin, centMax, PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1], PtBins[iPtA])); h2EtaPhiSBPeakCorAll[iPtT][iPtA] ->SetTitle(Form("Cent:[%d, %d], p_{T}^{Tr}:[%.1f, %.1f],p_{T}^{As}:[%.1f, %.1f];#Delta#eta; #Delta#phi", centMin, centMax, PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1], PtBins[iPtA], PtBins[iPtA+1])); } else{ //cout <<"iVtx:iPtT:iPtA = "<Add(h3EtaPhiSBCor[iVtx][iPtT][iPtA], 1.); h2EtaPhiSBPeakCorAll[iPtT][iPtA]->Add(h2EtaPhiSBPeakCor[iVtx][iPtT][iPtA], 1.); } } } //projection of mass to extract the yield in dPhi&dEta bins TH1F * h1PhiMass[nPtTrig][nPtAsso][nPhi]; //integrate eta from -0.7, .7 for peak region TH1F * h1PhiBulkMass[nPtTrig][nPtAsso][nPhi]; //integrate eta from [-1.2, -0.7]U[0.7, 1.2] for bulk TH1F * h1EtaMass[nPtTrig][nPtAsso][nEta]; //integrate phi from [-.5pi, .5pi] for peak region TH1F * h1Bulk1Mass[nPtTrig][nPtAsso]; // integrate eta from [0.7, 1.2]U[-1.2, -0.7] && phi from [-1., 1] TH1F * h1Bulk2Mass[nPtTrig][nPtAsso]; // integrate eta from [-1, 1] && phi from [1.05, 2.09] TH1F * h1PeakNSMass[nPtTrig][nPtAsso]; // integrate eta from [-0.7, 0.7] && phi from [-1., 1] TH1F * h1PeakASMass[nPtTrig][nPtAsso]; // integrate eta from [-1.5, 1.5] && phi from [0.5pi, 1.5pi] TH2F * h2EtaPhiSBCorAll[nPtTrig][nPtAsso]; TH1F * h1Phi[nPtTrig][nPtAsso]; TH1F * h1PhiPerTrig[nPtTrig][nPtAsso]; TH1F * h1PhiBulk[nPtTrig][nPtAsso]; TH1F * h1PhiBulkPerTrig[nPtTrig][nPtAsso]; TH1F * h1Eta[nPtTrig][nPtAsso]; TH1F * h1EtaPerTrig[nPtTrig][nPtAsso]; TH1F * h1Bulk1[nPtTrig]; TH1F * h1Bulk2[nPtTrig]; //to get NS ans AS yields without fitting phi or eta TH1F * h1PeakNS[nPtTrig]; TH1F * h1PeakAS[nPtTrig]; TH1F * h1JetNS[nPtTrig]; //after subtraction of bulk contribution TH1F * h1JetAS[nPtTrig]; for(Int_t i= 0; i < nPtTrig; i++){ h1Bulk1[i] = new TH1F(Form("h1Bulk1%d", i), "; p_{T} (GeV/c); dN/dp_{T}", nPtAsso, PtBins); h1Bulk1[i]->Sumw2(); h1Bulk2[i] = new TH1F(Form("h1Bulk2%d", i), "; p_{T} (GeV/c); dN/dp_{T}", nPtAsso, PtBins); h1Bulk2[i]->Sumw2(); h1PeakNS[i] = new TH1F(Form("h1PeakNS%d", i), "; p_{T} (GeV/c); dN/dp_{T}", nPtAsso, PtBins); h1PeakNS[i]->Sumw2(); h1PeakAS[i] = new TH1F(Form("h1PeakAS%d", i), "; p_{T} (GeV/c); dN/dp_{T}", nPtAsso, PtBins); h1PeakAS[i]->Sumw2(); h1JetNS[i] = new TH1F(Form("h1JetNS%d", i), "; p_{T} (GeV/c); dN/dp_{T}", nPtAsso, PtBins); h1JetNS[i]->Sumw2(); h1JetAS[i] = new TH1F(Form("h1JetAS%d", i), "; p_{T} (GeV/c); dN/dp_{T}", nPtAsso, PtBins); h1JetAS[i]->Sumw2(); } TH1F * h1Bulk1PerTrig[nPtTrig]; TH1F * h1Bulk2PerTrig[nPtTrig]; //to store per trigger yields for NS and AS jets TH1F * h1JetNSPerTrig[nPtTrig]; TH1F * h1JetASPerTrig[nPtTrig]; TH1F * h1JetEtaPerTrig[nPtTrig]; //yield in jet from eta distribution TH1F * h1JetPhiPerTrig[nPtTrig]; //yield in jet from phi distribution for(Int_t i= 0; i < nPtTrig; i++){ h1JetEtaPerTrig[i] = new TH1F(Form("h1JetEtaPerTrig%d", i), "; p_{T} (GeV/c); #frac{1}{N_{Trig}}#frac{dN}{dp_{T}}", nPtAsso, PtBins); h1JetPhiPerTrig[i] = new TH1F(Form("h1JetEtaPerTrig%d", i), "; p_{T} (GeV/c); #frac{1}{N_{Trig}}#frac{dN}{dp_{T}}", nPtAsso, PtBins); } TCanvas * cPhiMass[nPtTrig][nPtAsso]; TCanvas * cPhiBulkMass[nPtTrig][nPtAsso]; TCanvas * cEtaMass[nPtTrig][nPtAsso]; TCanvas * cBulk1Mass = new TCanvas("bulk1Mass", "", 800, 600); TCanvas * cBulk2Mass = new TCanvas("bulk2Mass", "", 800, 600);; cBulk1Mass->Divide(nPtAsso, nPtTrig); cBulk2Mass->Divide(nPtAsso, nPtTrig); TCanvas * cPeakNSMass = new TCanvas("PeakNSMass", "", 800, 600); TCanvas * cPeakASMass = new TCanvas("PeakASMass", "", 800, 600); cPeakNSMass->Divide(nPtAsso, nPtTrig); cPeakASMass->Divide(nPtAsso, nPtTrig); Char_t name[256]; Int_t iPad = 0; for(Int_t iPtT = 0; iPtT < nPtTrig; iPtT++){ for(Int_t iPtA = 0; iPtA < nPtAsso; iPtA++){ if(PtBins[iPtA] >= PtBinsV0Xi[iPtT+1]) break; cPhiMass[iPtT][iPtA] = new TCanvas(Form("cPhiMass%d%d", iPtT, iPtA), "", 800, 600); cPhiMass[iPtT][iPtA]->Divide(9, 4); //cPhiMass[iPtT][iPtA]->Divide(6, 3); cPhiBulkMass[iPtT][iPtA] = new TCanvas(Form("cPhiBulkMass%d%d", iPtT, iPtA), "", 800, 600); cPhiBulkMass[iPtT][iPtA]->Divide(9, 4); //cPhiBulkMass[iPtT][iPtA]->Divide(6, 3); cEtaMass[iPtT][iPtA] = new TCanvas(Form("cEtaMass%d%d", iPtT, iPtA), "", 800, 600); cEtaMass[iPtT][iPtA]->Divide(8, 4); h2EtaPhiSBCorAll[iPtT][iPtA] = (TH2F*)h3EtaPhiSBCorAll[iPtT][iPtA]->Project3D("xy"); h2EtaPhiSBCorAll[iPtT][iPtA]->SetName(Form("h2EtaPhiSBCorAllCent%d_%dPtTr%.1f_%.1fAs%.1f", centMin, centMax, PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1], PtBins[iPtA])); h2EtaPhiSBCorAll[iPtT][iPtA]->SetTitle(Form("Cent:[%d, %d], p_{T}^{Tr}:[%.1f, %.1f],p_{T}^{As}:[%.1f, %.1f];#Delta#eta; #Delta#phi", centMin, centMax, PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1], PtBins[iPtA], PtBins[iPtA+1])); //setup an empty histogram to store final results after background subtraction h1Phi[iPtT][iPtA] = (TH1F*)h3EtaPhiSBCorAll[iPtT][iPtA]->ProjectionX(); h1Phi[iPtT][iPtA]->SetName(Form("h1PhiCent%d_%dPtTr%.1f_%.1fAs%.1f", centMin, centMax, PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1], PtBins[iPtA])); h1Phi[iPtT][iPtA]->SetTitle(Form("Cent:[%d, %d], p_{T}^{Tr}:[%.1f, %.1f],p_{T}^{As}:[%.1f, %.1f]; #Delta#phi; dN/d#Delta#phi", centMin, centMax, PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1], PtBins[iPtA], PtBins[iPtA+1])); h1Phi[iPtT][iPtA]->Reset(); h1Phi[iPtT][iPtA]->Sumw2(); h1PhiBulk[iPtT][iPtA] = (TH1F*)h1Phi[iPtT][iPtA]->Clone();//h3EtaPhiSBCorAll[iPtT][iPtA]->ProjectionX(); h1PhiBulk[iPtT][iPtA]->SetName(Form("h1PhiBulkCent%d_%dPtTr%.1f_%.1fAs%.1f", centMin, centMax, PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1], PtBins[iPtA])); h1PhiBulk[iPtT][iPtA]->SetTitle(Form("Cent:[%d, %d], p_{T}^{Tr}:[%.1f, %.1f],p_{T}^{As}:[%.1f, %.1f]; #Delta#phi; dN/d#Delta#phi", centMin, centMax, PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1], PtBins[iPtA], PtBins[iPtA+1])); h1PhiBulk[iPtT][iPtA]->Reset(); h1PhiBulk[iPtT][iPtA]->Sumw2(); h1Eta[iPtT][iPtA] = (TH1F*)h3EtaPhiSBCorAll[iPtT][iPtA]->ProjectionY(); h1Eta[iPtT][iPtA]->SetName(Form("h1EtaCent%d_%dPtTr%.1f_%.1fAs%.1f", centMin, centMax, PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1], PtBins[iPtA])); h1Eta[iPtT][iPtA]->SetTitle(Form("Cent:[%d, %d], p_{T}^{Tr}:[%.1f, %.1f],p_{T}^{As}:[%.1f, %.1f]; #Delta#eta; dN/d#Delta#eta", centMin, centMax, PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1], PtBins[iPtA], PtBins[iPtA+1])); h1Eta[iPtT][iPtA]->Reset(); h1Eta[iPtT][iPtA]->Sumw2(); //projection in phi for peak region with S+B Int_t etaLow = h3EtaPhiSBCorAll[iPtT][iPtA]->GetYaxis()->FindBin(-.7+0.005); Int_t etaHigh = h3EtaPhiSBCorAll[iPtT][iPtA]->GetYaxis()->FindBin(.7-0.005); Double_t etaEdgeLow = h3EtaPhiSBCorAll[iPtT][iPtA]->GetYaxis()->GetBinLowEdge(etaLow); Double_t etaEdgeHigh = h3EtaPhiSBCorAll[iPtT][iPtA]->GetYaxis()->GetBinLowEdge(etaHigh+1); Double_t peakEtaRange = etaEdgeHigh-etaEdgeLow; for(Int_t iPhi = 0; iPhi < h3EtaPhiSBCorAll[iPtT][iPtA]->GetXaxis()->GetNbins(); iPhi++){ sprintf(name, "h1PhiMassCent%d_%dPtTr%.1f_%.1fAs%.1fPhi%d",centMin, centMax, PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1], PtBins[iPtA], iPhi); h1PhiMass[iPtT][iPtA][iPhi] = (TH1F*)h3EtaPhiSBCorAll[iPtT][iPtA]->ProjectionZ(name, iPhi+1, iPhi+1, etaLow, etaHigh); h1PhiMass[iPtT][iPtA][iPhi]->SetTitle(Form("Cent:[%d, %d], p_{T}^{Tr}:[%.1f, %.1f],p_{T}^{As}:[%.1f, %.1f], #Delta#phi:bin %d; mass (GeV/c^{2})", centMin, centMax, PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1], PtBins[iPtA], PtBins[iPtA+1], iPhi+1)); cPhiMass[iPtT][iPtA]->cd(iPhi+1); h1PhiMass[iPtT][iPtA][iPhi]->Draw(); Double_t sig, sigErr; fitHist(h1PhiMass[iPtT][iPtA][iPhi], massMean, massRes, sig, sigErr); h1Phi[iPtT][iPtA]->SetBinContent(iPhi+1, sig); h1Phi[iPtT][iPtA]->SetBinError(iPhi+1, sigErr); } //normalization, but not for the eta range to extract the yield //h1Phi[iPtT][iPtA]->Scale(1./(etaEdgeHigh-etaEdgeLow)); //h1Phi[iPtT][iPtA]->Scale(1./h1Phi[iPtT][iPtA]->GetBinWidth(1)); //projection in phi for bulk region in order to extract background in the peak region etaLow = h3EtaPhiSBCorAll[iPtT][iPtA]->GetYaxis()->FindBin(-1.2+0.005); etaHigh = h3EtaPhiSBCorAll[iPtT][iPtA]->GetYaxis()->FindBin(-.7-0.005); Int_t etaLow1 = h3EtaPhiSBCorAll[iPtT][iPtA]->GetYaxis()->FindBin(.7+0.005); Int_t etaHigh1 = h3EtaPhiSBCorAll[iPtT][iPtA]->GetYaxis()->FindBin(1.2-0.005); etaEdgeLow = h3EtaPhiSBCorAll[iPtT][iPtA]->GetYaxis()->GetBinLowEdge(etaLow); etaEdgeHigh = h3EtaPhiSBCorAll[iPtT][iPtA]->GetYaxis()->GetBinLowEdge(etaHigh+1); Double_t bkgEtaRange = 2.*(etaEdgeHigh-etaEdgeLow); for(Int_t iPhi = 0; iPhi < h3EtaPhiSBCorAll[iPtT][iPtA]->GetXaxis()->GetNbins(); iPhi++){ sprintf(name, "h1PhiBulkMassCent%d_%dPtTr%.1f_%.1fAs%.1fPhi%d",centMin, centMax, PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1], PtBins[iPtA], iPhi); h1PhiBulkMass[iPtT][iPtA][iPhi] = (TH1F*)h3EtaPhiSBCorAll[iPtT][iPtA]->ProjectionZ(name, iPhi+1, iPhi+1, etaLow, etaHigh); h1PhiBulkMass[iPtT][iPtA][iPhi] ->SetTitle(Form("Cent:[%d, %d], p_{T}^{Tr}:[%.1f, %.1f], p_{T}^{As}:[%.1f, %.1f], #Delta#phi:bin %d; mass (GeV/c^{2})", centMin, centMax, PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1], PtBins[iPtA], PtBins[iPtA+1], iPhi+1)); TH1F * htemp = (TH1F*)h3EtaPhiSBCorAll[iPtT][iPtA]->ProjectionZ(Form("htemp%d%d%d", iPtT, iPtA, iPhi), iPhi+1, iPhi+1, etaLow1, etaHigh1); h1PhiBulkMass[iPtT][iPtA][iPhi]->Add(htemp, 1.); cPhiBulkMass[iPtT][iPtA]->cd(iPhi+1); h1PhiBulkMass[iPtT][iPtA][iPhi]->Draw(); Double_t sig, sigErr; fitHist(h1PhiBulkMass[iPtT][iPtA][iPhi], massMean, massRes, sig, sigErr); h1PhiBulk[iPtT][iPtA]->SetBinContent(iPhi+1, sig); h1PhiBulk[iPtT][iPtA]->SetBinError(iPhi+1, sigErr); } //normalization h1PhiBulk[iPtT][iPtA]->Scale(peakEtaRange/bkgEtaRange); //scale to the same eta range as peak region for comparison // h1PhiBulk[iPtT][iPtA]->Scale(1./h1PhiBulk[iPtT][iPtA]->GetBinWidth(1)); //projecton in eta Int_t phiLow = h3EtaPhiSBCorAll[iPtT][iPtA]->GetXaxis()->FindBin(-0.5*TMath::Pi()+0.005); Int_t phiHigh = h3EtaPhiSBCorAll[iPtT][iPtA]->GetYaxis()->FindBin(0.5*TMath::Pi()-0.005); Double_t phiEdgeLow = h3EtaPhiSBCorAll[iPtT][iPtA]->GetXaxis()->GetBinLowEdge(phiLow); Double_t phiEdgeHigh = h3EtaPhiSBCorAll[iPtT][iPtA]->GetXaxis()->GetBinLowEdge(phiHigh+1); for(Int_t iEta = 0; iEta < h3EtaPhiSBCorAll[iPtT][iPtA]->GetYaxis()->GetNbins(); iEta++){ sprintf(name, "h1EtaMassCent%d_%dPtTr%.1f_%.1fAs%.1fEta%d",centMin, centMax, PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1], PtBins[iPtA], iEta); h1EtaMass[iPtT][iPtA][iEta] = (TH1F*)h3EtaPhiSBCorAll[iPtT][iPtA]->ProjectionZ(name, phiLow, phiHigh, iEta+1, iEta+1); h1EtaMass[iPtT][iPtA][iEta]->SetTitle(Form("Cent:[%d, %d], p_{T}^{Tr}:[%.1f, %.1f],p_{T}^{As}:[%.1f, %.1f], #Delta#eta:bin %d; mass (GeV/c^{2})", centMin, centMax, PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1], PtBins[iPtA], PtBins[iPtA+1], iEta+1)); cEtaMass[iPtT][iPtA]->cd(iEta+1); h1EtaMass[iPtT][iPtA][iEta]->Draw(); Double_t sig, sigErr; fitHist(h1EtaMass[iPtT][iPtA][iEta], massMean, massRes, sig, sigErr); h1Eta[iPtT][iPtA]->SetBinContent(iEta+1, sig); h1Eta[iPtT][iPtA]->SetBinError(iEta+1, sigErr);} //normalization //h1Eta[iPtT][iPtA]->Scale(1./(phiEdgeHigh-phiEdgeLow)); //not need for the yield extraction //h1Eta[iPtT][iPtA]->Scale(1./h1Eta[iPtT][iPtA]->GetBinWidth(1)); //projection in bulk1 etaLow = h3EtaPhiSBCorAll[iPtT][iPtA]->GetYaxis()->FindBin(-1.2+0.005); etaHigh = h3EtaPhiSBCorAll[iPtT][iPtA]->GetYaxis()->FindBin(-.7-0.005); phiLow = h3EtaPhiSBCorAll[iPtT][iPtA]->GetXaxis()->FindBin(-1.+0.005); phiHigh = h3EtaPhiSBCorAll[iPtT][iPtA]->GetXaxis()->FindBin(1.-0.005); //second region etaLow1 = h3EtaPhiSBCorAll[iPtT][iPtA]->GetYaxis()->FindBin(.7+0.005); etaHigh1 = h3EtaPhiSBCorAll[iPtT][iPtA]->GetYaxis()->FindBin(1.2-0.005); etaEdgeLow = h3EtaPhiSBCorAll[iPtT][iPtA]->GetYaxis()->GetBinLowEdge(etaLow); etaEdgeHigh = h3EtaPhiSBCorAll[iPtT][iPtA]->GetYaxis()->GetBinLowEdge(etaHigh+1); phiEdgeLow = h3EtaPhiSBCorAll[iPtT][iPtA]->GetXaxis()->GetBinLowEdge(phiLow); phiEdgeHigh = h3EtaPhiSBCorAll[iPtT][iPtA]->GetXaxis()->GetBinLowEdge(phiHigh+1); Double_t norm1 = 2.*(etaEdgeHigh-etaEdgeLow)*(phiEdgeHigh-phiEdgeLow); sprintf(name, "h1Bulk1MassCent%d_%dPtTr%.1f_%.1fAs%.1f",centMin, centMax, PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1], PtBins[iPtA]); h1Bulk1Mass[iPtT][iPtA] = (TH1F*) h3EtaPhiSBCorAll[iPtT][iPtA]->ProjectionZ(name, phiLow, phiHigh, etaLow, etaHigh); h1Bulk1Mass[iPtT][iPtA] ->SetTitle(Form("Cent:[%d, %d], p_{T}^{Tr}:[%.1f, %.1f],p_{T}^{As}:[%.1f, %.1f]; mass (GeV/c^{2})", centMin, centMax, PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1], PtBins[iPtA], PtBins[iPtA+1])); TH1F * h1temp = (TH1F*) h3EtaPhiSBCorAll[iPtT][iPtA]->ProjectionZ(Form("h1temp%d%d", iPtT, iPtA), phiLow, phiHigh, etaLow1, etaHigh1); h1Bulk1Mass[iPtT][iPtA]->Add(h1temp, 1.); cBulk1Mass->cd(iPad+1); h1Bulk1Mass[iPtT][iPtA]->Draw(); Double_t sig1, sigErr1; fitHist(h1Bulk1Mass[iPtT][iPtA], massMean, massRes, sig1, sigErr1); Double_t sigBulk1 = sig1; Double_t sigErrBulk1 = sigErr1; h1Bulk1[iPtT]->SetBinContent(iPtA+1, sig1/h1Bulk1[iPtT]->GetBinWidth(iPtA+1)); h1Bulk1[iPtT]->SetBinError(iPtA+1, sigErr1/h1Bulk1[iPtT]->GetBinWidth(iPtA+1)); //projection in bulk2 etaLow = h3EtaPhiSBCorAll[iPtT][iPtA]->GetYaxis()->FindBin(-1.+0.005); etaHigh = h3EtaPhiSBCorAll[iPtT][iPtA]->GetYaxis()->FindBin(1.-0.005); phiLow = h3EtaPhiSBCorAll[iPtT][iPtA]->GetXaxis()->FindBin(1.05+0.005); phiHigh = h3EtaPhiSBCorAll[iPtT][iPtA]->GetXaxis()->FindBin(2.09-0.005); etaEdgeLow = h3EtaPhiSBCorAll[iPtT][iPtA]->GetYaxis()->GetBinLowEdge(etaLow); etaEdgeHigh = h3EtaPhiSBCorAll[iPtT][iPtA]->GetYaxis()->GetBinLowEdge(etaHigh+1); phiEdgeLow = h3EtaPhiSBCorAll[iPtT][iPtA]->GetXaxis()->GetBinLowEdge(phiLow); phiEdgeHigh = h3EtaPhiSBCorAll[iPtT][iPtA]->GetXaxis()->GetBinLowEdge(phiHigh+1); Double_t norm2 = (etaEdgeHigh-etaEdgeLow)*(phiEdgeHigh-phiEdgeLow); sprintf(name, "h1Bulk2MassCent%d_%dPtTr%.1f_%.1fAs%.1f",centMin, centMax, PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1], PtBins[iPtA]); h1Bulk2Mass[iPtT][iPtA] = (TH1F*) h3EtaPhiSBCorAll[iPtT][iPtA]->ProjectionZ(name, phiLow, phiHigh, etaLow, etaHigh); h1Bulk2Mass[iPtT][iPtA] ->SetTitle(Form("Cent:[%d, %d], p_{T}^{Tr}:[%.1f, %.1f],p_{T}^{As}:[%.1f, %.1f]; mass (GeV/c^{2})", centMin, centMax, PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1], PtBins[iPtA], PtBins[iPtA+1])); cBulk2Mass->cd(iPad+1); h1Bulk2Mass[iPtT][iPtA]->Draw(); fitHist(h1Bulk2Mass[iPtT][iPtA], massMean, massRes, sig1, sigErr1); h1Bulk2[iPtT]->SetBinContent(iPtA+1, sig1/h1Bulk2[iPtT]->GetBinWidth(iPtA+1)); h1Bulk2[iPtT]->SetBinError(iPtA+1, sigErr1/h1Bulk2[iPtT]->GetBinWidth(iPtA+1)); //projection in NS peak etaLow = h3EtaPhiSBCorAll[iPtT][iPtA]->GetYaxis()->FindBin(-.7+0.005); etaHigh = h3EtaPhiSBCorAll[iPtT][iPtA]->GetYaxis()->FindBin(.7-0.005); phiLow = h3EtaPhiSBCorAll[iPtT][iPtA]->GetXaxis()->FindBin(-1.+0.005); phiHigh = h3EtaPhiSBCorAll[iPtT][iPtA]->GetXaxis()->FindBin(1.-0.005); etaEdgeLow = h3EtaPhiSBCorAll[iPtT][iPtA]->GetYaxis()->GetBinLowEdge(etaLow); etaEdgeHigh = h3EtaPhiSBCorAll[iPtT][iPtA]->GetYaxis()->GetBinLowEdge(etaHigh+1); phiEdgeLow = h3EtaPhiSBCorAll[iPtT][iPtA]->GetXaxis()->GetBinLowEdge(phiLow); phiEdgeHigh = h3EtaPhiSBCorAll[iPtT][iPtA]->GetXaxis()->GetBinLowEdge(phiHigh+1); Double_t norm3 = (etaEdgeHigh-etaEdgeLow)*(phiEdgeHigh-phiEdgeLow); sprintf(name, "h1PeakNSMassCent%d_%dPtTr%.1f_%.1fAs%.1f",centMin, centMax, PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1], PtBins[iPtA]); h1PeakNSMass[iPtT][iPtA] = (TH1F*) h3EtaPhiSBCorAll[iPtT][iPtA]->ProjectionZ(name, phiLow, phiHigh, etaLow, etaHigh); h1PeakNSMass[iPtT][iPtA] ->SetTitle(Form("Cent:[%d, %d], p_{T}^{Tr}:[%.1f, %.1f],p_{T}^{As}:[%.1f, %.1f]; mass (GeV/c^{2})", centMin, centMax, PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1], PtBins[iPtA], PtBins[iPtA+1])); cPeakNSMass->cd(iPad+1); h1PeakNSMass[iPtT][iPtA]->Draw(); fitHist(h1PeakNSMass[iPtT][iPtA], massMean, massRes, sig1, sigErr1); h1PeakNS[iPtT]->SetBinContent(iPtA+1, sig1/h1PeakNS[iPtT]->GetBinWidth(iPtA+1)); h1PeakNS[iPtT]->SetBinError(iPtA+1, sigErr1/h1PeakNS[iPtT]->GetBinWidth(iPtA+1)); Double_t psRatio = norm3/norm1; //phase space ratio Double_t sigJetNS = sig1-sigBulk1*psRatio; Double_t sigErrJetNS = TMath::Sqrt(sigErr1*sigErr1+sigErrBulk1*sigErrBulk1*psRatio*psRatio); h1JetNS[iPtT]->SetBinContent(iPtA+1, sigJetNS/h1PeakNS[iPtT]->GetBinWidth(iPtA+1)); h1JetNS[iPtT]->SetBinError(iPtA+1, sigErrJetNS/h1PeakNS[iPtT]->GetBinWidth(iPtA+1)); //projection in AS peak etaLow = h3EtaPhiSBCorAll[iPtT][iPtA]->GetYaxis()->FindBin(-1.5+0.005); etaHigh = h3EtaPhiSBCorAll[iPtT][iPtA]->GetYaxis()->FindBin(1.5-0.005); phiLow = h3EtaPhiSBCorAll[iPtT][iPtA]->GetXaxis()->FindBin(-.5*TMath::Pi()+0.005); phiHigh = h3EtaPhiSBCorAll[iPtT][iPtA]->GetXaxis()->FindBin(1.5*TMath::Pi()-0.005); etaEdgeLow = h3EtaPhiSBCorAll[iPtT][iPtA]->GetYaxis()->GetBinLowEdge(etaLow); etaEdgeHigh = h3EtaPhiSBCorAll[iPtT][iPtA]->GetYaxis()->GetBinLowEdge(etaHigh+1); phiEdgeLow = h3EtaPhiSBCorAll[iPtT][iPtA]->GetXaxis()->GetBinLowEdge(phiLow); phiEdgeHigh = h3EtaPhiSBCorAll[iPtT][iPtA]->GetXaxis()->GetBinLowEdge(phiHigh+1); Double_t norm4 = (etaEdgeHigh-etaEdgeLow)*(phiEdgeHigh-phiEdgeLow); sprintf(name, "h1PeakASMassCent%d_%dPtTr%.1f_%.1fAs%.1f",centMin, centMax, PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1], PtBins[iPtA]); h1PeakASMass[iPtT][iPtA] = (TH1F*) h3EtaPhiSBCorAll[iPtT][iPtA]->ProjectionZ(name, phiLow, phiHigh, etaLow, etaHigh); h1PeakASMass[iPtT][iPtA] ->SetTitle(Form("Cent:[%d, %d], p_{T}^{Tr}:[%.1f, %.1f],p_{T}^{As}:[%.1f, %.1f]; mass (GeV/c^{2})", centMin, centMax, PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1], PtBins[iPtA], PtBins[iPtA+1])); cPeakASMass->cd(iPad+1); h1PeakASMass[iPtT][iPtA]->Draw(); fitHist(h1PeakASMass[iPtT][iPtA], massMean, massRes, sig1, sigErr1); h1PeakAS[iPtT]->SetBinContent(iPtA+1, sig1/h1PeakAS[iPtT]->GetBinWidth(iPtA+1)); h1PeakAS[iPtT]->SetBinError(iPtA+1, sigErr1/h1PeakAS[iPtT]->GetBinWidth(iPtA+1)); psRatio = norm4/norm1; //phase space ratio Double_t sigJetAS = sig1-sigBulk1*psRatio; Double_t sigErrJetAS = TMath::Sqrt(sigErr1*sigErr1+sigErrBulk1*sigErrBulk1*psRatio*psRatio); h1JetAS[iPtT]->SetBinContent(iPtA+1, sigJetAS/h1PeakNS[iPtT]->GetBinWidth(iPtA+1)); h1JetAS[iPtT]->SetBinError(iPtA+1, sigErrJetAS/h1PeakNS[iPtT]->GetBinWidth(iPtA+1)); iPad++; } //normalization h1Bulk1[iPtT]->Scale(1./norm1); h1Bulk2[iPtT]->Scale(1./norm2); //Not should be divided by the p_T bin width if not equal to 1. } //calculate the per trigger yields for(Int_t iPtT = 0; iPtT < nPtTrig; iPtT++){ Int_t ptTrigBinLow = hNtrig->GetXaxis()->FindBin(PtBinsV0Xi[iPtT]+0.01); Int_t ptTrigBinHigh = hNtrig->GetXaxis()->FindBin(PtBinsV0Xi[iPtT+1]-0.01); Double_t ntrig = hNtrig->Integral(ptTrigBinLow, ptTrigBinHigh); h1Bulk1PerTrig[iPtT] = (TH1F*)h1Bulk1[iPtT]->Clone(Form("h1Bulk1PerTrig%d", iPtT)); h1Bulk1PerTrig[iPtT]->SetTitle("; p_{T} (GeV/c); #frac{1}{N_{Trig}}#frac{dN}{d#Delta#phi d#Delta#eta dp_{T}}"); h1Bulk2PerTrig[iPtT] = (TH1F*)h1Bulk2[iPtT]->Clone(Form("h1Bulk2PerTrig%d", iPtT)); h1Bulk2PerTrig[iPtT]->SetTitle("; p_{T} (GeV/c); #frac{1}{N_{Trig}}#frac{dN}{d#Delta#phi d#Delta#eta dp_{T}}"); h1JetNSPerTrig[iPtT] = (TH1F*)h1JetNS[iPtT]->Clone(Form("h1JetNSPerTrig%d", iPtT)); h1JetNSPerTrig[iPtT]->SetTitle("; p_{T} (GeV/c); #frac{1}{N_{Trig}}#frac{dN}{dp_{T}}"); h1JetASPerTrig[iPtT] = (TH1F*)h1JetAS[iPtT]->Clone(Form("h1JetASPerTrig%d", iPtT)); h1JetASPerTrig[iPtT]->SetTitle("; p_{T} (GeV/c); #frac{1}{N_{Trig}}#frac{dN}{dp_{T}}"); if(ntrig!=0){ h1Bulk1PerTrig[iPtT]->Scale(1./ntrig); //shall scaled by 2pi? No! h1Bulk2PerTrig[iPtT]->Scale(1./ntrig); //shall scaled by 2pi? h1JetNSPerTrig[iPtT]->Scale(1./ntrig); h1JetASPerTrig[iPtT]->Scale(1./ntrig); } for(Int_t iPtA = 0; iPtA < nPtAsso; iPtA++){ if(PtBins[iPtA] >= PtBinsV0Xi[iPtT+1]) break; h1PhiPerTrig[iPtT][iPtA] = (TH1F*)h1Phi[iPtT][iPtA]->Clone(Form("h1PhiPerTrigCent%d_%dPtTr%.1f_%.1fAs%.1f", centMin, centMax, PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1], PtBins[iPtA])); h1PhiPerTrig[iPtT][iPtA] ->SetTitle(Form("Cent:[%d, %d], p_{T}^{Tr}:[%.1f, %.1f],p_{T}^{As}:[%.1f, %.1f]; d#Delta#phi; #frac{1}{N_{Trig}}#frac{dN}{d#Delta#phi}", centMin, centMax, PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1], PtBins[iPtA], PtBins[iPtA+1])); h1PhiBulkPerTrig[iPtT][iPtA] = (TH1F*)h1PhiBulk[iPtT][iPtA]->Clone(Form("h1PhiBulkPerTrigCent%d_%dPtTr%.1f_%.1fAs%.1f", centMin, centMax, PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1], PtBins[iPtA])); h1PhiBulkPerTrig[iPtT][iPtA] ->SetTitle(Form("Cent:[%d, %d], p_{T}^{Tr}:[%.1f, %.1f],p_{T}^{As}:[%.1f, %.1f]; d#Delta#phi; #frac{1}{N_{Trig}}#frac{dN}{d#Delta#phi}", centMin, centMax, PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1], PtBins[iPtA], PtBins[iPtA+1])); h1EtaPerTrig[iPtT][iPtA] = (TH1F*)h1Eta[iPtT][iPtA]->Clone(Form("h1EtaPerTrigCent%d_%dPtTr%.1f_%.1fAs%.1f", centMin, centMax, PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1], PtBins[iPtA])); h1EtaPerTrig[iPtT][iPtA] ->SetTitle(Form("Cent:[%d, %d], p_{T}^{Tr}:[%.1f, %.1f],p_{T}^{As}:[%.1f, %.1f]; d#Delta#eta; #frac{1}{N_{Trig}}#frac{dN}{d#Delta#eta}", centMin, centMax, PtBinsV0Xi[iPtT], PtBinsV0Xi[iPtT+1], PtBins[iPtA], PtBins[iPtA+1])); if(ntrig != 0){ //also divided by the bin width, not useful for the yield extraction h1PhiPerTrig[iPtT][iPtA]->Scale(1./ntrig/h1PhiPerTrig[iPtT][iPtA]->GetBinWidth(1)); h1PhiBulkPerTrig[iPtT][iPtA]->Scale(1./ntrig/h1PhiBulkPerTrig[iPtT][iPtA]->GetBinWidth(1)); h1EtaPerTrig[iPtT][iPtA]->Scale(1./ntrig/h1EtaPerTrig[iPtT][iPtA]->GetBinWidth(1)); } } } //extract yield in jet TCanvas *cJetEta = new TCanvas("cJetEta", "", 800, 600); cJetEta->Divide(nPtAsso, nPtTrig); iPad = 0; for(Int_t iPtT = 0; iPtT < nPtTrig; iPtT++){ for(Int_t iPtA = 0; iPtA < nPtAsso; iPtA++){ if(PtBins[iPtA] >= PtBinsV0Xi[iPtT+1]) break; cJetEta->cd(iPad+1); h1EtaPerTrig[iPtT][iPtA]->Draw(); Double_t sig, sigErr; //already normalized fitHist( h1EtaPerTrig[iPtT][iPtA], 0, 0.15, sig, sigErr, kTRUE); h1JetEtaPerTrig[iPtT]->SetBinContent(iPtA+1, sig); h1JetEtaPerTrig[iPtT]->SetBinError(iPtA+1, sigErr); iPad++; } } TFile * outFile = 0x0; if(!bAnalyticCorr) outFile = new TFile(Form("AnalysisResultsSpecie%dCent%d_%dAll.root", centMin, centMax ), //, specie, centMin, centMax "RECREATE"); else outFile = new TFile(Form("AnalysisResultsSpecie%dCent%d_%dAllAnalytic.root", centMin, centMax), "RECREATE"); outFile->cd(); hNtrig->Write(); for(Int_t iVtx = 0; iVtx < nVtx; iVtx++){ for(Int_t iPtT = 0; iPtT < nPtTrig; iPtT++) for(Int_t iPtA = 0; iPtA < nPtAsso; iPtA++){ if(PtBins[iPtA] >= PtBinsV0Xi[iPtT+1]) break; hMass[iVtx][iPtT][iPtA]->Write(); h3EtaPhiSB[iVtx][iPtT][iPtA]->Write(); h2EtaPhiSB[iVtx][iPtT][iPtA]->Write(); if(bMix) h2EtaPhiSBMix[iVtx][iPtT][iPtA]->Write(); h3EtaPhiSBCor[iVtx][iPtT][iPtA]->Write(); h2EtaPhiSBPeakCor[iVtx][iPtT][iPtA]->Write(); } } for(Int_t iPtT = 0; iPtT < nPtTrig; iPtT++) for(Int_t iPtA = 0; iPtA < nPtAsso; iPtA++){ if(PtBins[iPtA] >= PtBinsV0Xi[iPtT+1]) break; h3EtaPhiSBCorAll[iPtT][iPtA]->Write(); h2EtaPhiSBPeakCorAll[iPtT][iPtA]->Write(); h2EtaPhiSBCorAll[iPtT][iPtA]->Write(); } for(Int_t iPtT = 0; iPtT < nPtTrig; iPtT++) for(Int_t iPtA = 0; iPtA < nPtAsso; iPtA++){ if(PtBins[iPtA] >= PtBinsV0Xi[iPtT+1]) break; for(Int_t iPhi = 0; iPhi < h3EtaPhiSBCorAll[iPtT][iPtA]->GetXaxis()->GetNbins(); iPhi++){ h1PhiMass[iPtT][iPtA][iPhi]->Write(); h1PhiBulkMass[iPtT][iPtA][iPhi]->Write(); } for(Int_t iEta = 0; iEta < h3EtaPhiSBCorAll[iPtT][iPtA]->GetYaxis()->GetNbins(); iEta++){ h1EtaMass[iPtT][iPtA][iEta]->Write(); } h1Bulk1Mass[iPtT][iPtA]->Write(); h1Bulk2Mass[iPtT][iPtA]->Write(); h1PeakNSMass[iPtT][iPtA]->Write(); h1PeakASMass[iPtT][iPtA]->Write(); } for(Int_t iPtT = 0; iPtT < nPtTrig; iPtT++){ for(Int_t iPtA = 0; iPtA < nPtAsso; iPtA++){ if(PtBins[iPtA] >= PtBinsV0Xi[iPtT+1]) break; h1Phi[iPtT][iPtA]->Write(); h1PhiBulk[iPtT][iPtA]->Write(); h1Eta[iPtT][iPtA]->Write(); h1PhiPerTrig[iPtT][iPtA]->Write(); h1PhiBulkPerTrig[iPtT][iPtA]->Write(); h1EtaPerTrig[iPtT][iPtA]->Write(); } h1Bulk1[iPtT]->Write(); h1Bulk2[iPtT]->Write(); h1Bulk1PerTrig[iPtT]->Write(); h1Bulk2PerTrig[iPtT]->Write(); h1PeakNS[iPtT]->Write(); h1PeakAS[iPtT]->Write(); h1JetNS[iPtT]->Write(); h1JetAS[iPtT]->Write(); h1JetNSPerTrig[iPtT]->Write(); h1JetASPerTrig[iPtT]->Write(); h1JetEtaPerTrig[iPtT]->Write(); } }