//________________________________________________________________________________________________________________________________________________
AliAnaPi0::AliAnaPi0() : AliAnaPartCorrBaseClass(),
-fDoOwnMix(kFALSE),fNCentrBin(0),//fNZvertBin(0),fNrpBin(0),
-fNmaxMixEv(0), fCalorimeter(""),
-fNModules(12), fUseAngleCut(kFALSE), fUseAngleEDepCut(kFALSE),fAngleCut(0), fAngleMaxCut(7.),fEventsList(0x0), fMultiCutAna(kFALSE), fMultiCutAnaSim(kFALSE),
-fNPtCuts(0),fNAsymCuts(0), fNCellNCuts(0),fNPIDBits(0), fMakeInvPtPlots(kFALSE), fSameSM(kFALSE),
+fDoOwnMix(kFALSE), fCalorimeter(""),
+fNModules(12), fUseAngleCut(kFALSE), fUseAngleEDepCut(kFALSE),fAngleCut(0), fAngleMaxCut(7.),fEventsList(0x0),
+fMultiCutAna(kFALSE), fMultiCutAnaSim(kFALSE),
+fNPtCuts(0),fNAsymCuts(0), fNCellNCuts(0),fNPIDBits(0),
+fMakeInvPtPlots(kFALSE), fSameSM(kFALSE), fFillSMCombinations(kFALSE), fCheckConversion(kFALSE),
fUseTrackMultBins(kFALSE),fUsePhotonMultBins(kFALSE),fUseAverCellEBins(kFALSE), fUseAverClusterEBins(kFALSE),
fUseAverClusterEDenBins(0), //fUseAverClusterPairRBins(0), fUseAverClusterPairRWeightBins(0), fUseEMaxBins(0),
fFillBadDistHisto(kFALSE),
fhReInvPt1(0x0), fhMiInvPt1(0x0), fhReInvPt2(0x0), fhMiInvPt2(0x0), fhReInvPt3(0x0), fhMiInvPt3(0x0),
fhRePtNCellAsymCuts(0x0), fhRePtNCellAsymCutsSM0(0x0), fhRePtNCellAsymCutsSM1(0x0), fhRePtNCellAsymCutsSM2(0x0), fhRePtNCellAsymCutsSM3(0x0), fhMiPtNCellAsymCuts(0x0),
fhRePIDBits(0x0),fhRePtMult(0x0), fhRePtAsym(0x0), fhRePtAsymPi0(0x0),fhRePtAsymEta(0x0),
-fhEvents(0x0), fhCentrality(0x0),
+fhEvents(0x0), fhCentrality(0x0),fhCentralityNoPair(0x0),
+fhEventPlaneAngle(0x0), fhEventPlaneResolution(0x0),
fhRealOpeningAngle(0x0),fhRealCosOpeningAngle(0x0), fhMixedOpeningAngle(0x0),fhMixedCosOpeningAngle(0x0),
fhPrimPi0Pt(0x0), fhPrimPi0AccPt(0x0), fhPrimPi0Y(0x0), fhPrimPi0AccY(0x0), fhPrimPi0Phi(0x0), fhPrimPi0AccPhi(0x0),
fhPrimPi0OpeningAngle(0x0), fhPrimPi0CosOpeningAngle(0x0),
// Remove event containers
if(fDoOwnMix && fEventsList){
- for(Int_t ic=0; ic<fNCentrBin; ic++){
+ for(Int_t ic=0; ic<GetNCentrBin(); ic++){
for(Int_t iz=0; iz<GetNZvertBin(); iz++){
for(Int_t irp=0; irp<GetNRPBin(); irp++){
fEventsList[ic*GetNZvertBin()*GetNRPBin()+iz*GetNRPBin()+irp]->Delete() ;
AddToHistogramsName("AnaPi0_");
fNModules = 12; // set maximum to maximum number of EMCAL modules
- fNCentrBin = 1;
-// fNZvertBin = 1;
-// fNrpBin = 1;
- fNmaxMixEv = 10;
fCalorimeter = "PHOS";
fUseAngleCut = kFALSE;
fMultiCutAna = kFALSE;
- fNPtCuts = 3;
+ fNPtCuts = 1;
fPtCuts[0] = 0.; fPtCuts[1] = 0.3; fPtCuts[2] = 0.5;
for(Int_t i = fNPtCuts; i < 10; i++)fPtCuts[i] = 0.;
- fNAsymCuts = 4;
- fAsymCuts[0] = 1.; fAsymCuts[1] = 0.8; fAsymCuts[2] = 0.6; fAsymCuts[3] = 0.1;
+ fNAsymCuts = 2;
+ fAsymCuts[0] = 1.; fAsymCuts[1] = 0.7; //fAsymCuts[2] = 0.6; // fAsymCuts[3] = 0.1;
for(Int_t i = fNAsymCuts; i < 10; i++)fAsymCuts[i] = 0.;
- fNCellNCuts = 3;
+ fNCellNCuts = 1;
fCellNCuts[0] = 0; fCellNCuts[1] = 1; fCellNCuts[2] = 2;
for(Int_t i = fNCellNCuts; i < 10; i++)fCellNCuts[i] = 0;
- fNPIDBits = 2;
+ fNPIDBits = 1;
fPIDBits[0] = 0; fPIDBits[1] = 2; // fPIDBits[2] = 4; fPIDBits[3] = 6;// check, no cut, dispersion, neutral, dispersion&&neutral
for(Int_t i = fNPIDBits; i < 10; i++)fPIDBits[i] = 0;
char onePar[buffersize] ;
snprintf(onePar,buffersize,"--- AliAnaPi0 ---\n") ;
parList+=onePar ;
- snprintf(onePar,buffersize,"Number of bins in Centrality: %d \n",fNCentrBin) ;
+ snprintf(onePar,buffersize,"Number of bins in Centrality: %d \n",GetNCentrBin()) ;
parList+=onePar ;
snprintf(onePar,buffersize,"Number of bins in Z vert. pos: %d \n",GetNZvertBin()) ;
parList+=onePar ;
snprintf(onePar,buffersize,"Number of bins in Reac. Plain: %d \n",GetNRPBin()) ;
parList+=onePar ;
- snprintf(onePar,buffersize,"Depth of event buffer: %d \n",fNmaxMixEv) ;
+ snprintf(onePar,buffersize,"Depth of event buffer: %d \n",GetNMaxEvMix()) ;
parList+=onePar ;
- snprintf(onePar,buffersize,"Pair in same Module: %d ; TrackMult as centrality: %d; PhotonMult as centrality: %d; cluster E as centrality: %d; cell as centrality: %d; Fill InvPt histos %d\n",
- fSameSM, fUseTrackMultBins, fUsePhotonMultBins, fUseAverClusterEBins, fUseAverCellEBins, fMakeInvPtPlots) ;
+ snprintf(onePar,buffersize,"Pair in same Module: %d ; Fill Different SM histos %d; CheckConversions %d; TrackMult as centrality: %d; PhotonMult as centrality: %d; cluster E as centrality: %d; cell as centrality: %d; Fill InvPt histos %d\n",
+ fSameSM, fFillSMCombinations, fCheckConversion, fUseTrackMultBins, fUsePhotonMultBins, fUseAverClusterEBins, fUseAverCellEBins, fMakeInvPtPlots) ;
parList+=onePar ;
snprintf(onePar,buffersize,"Select pairs with their angle: %d, edep %d, min angle %2.3f, max angle %2.3f,\n",fUseAngleCut, fUseAngleEDepCut,fAngleCut,fAngleMaxCut) ;
parList+=onePar ;
// store them in fOutputContainer
//create event containers
- fEventsList = new TList*[fNCentrBin*GetNZvertBin()*GetNRPBin()] ;
+ fEventsList = new TList*[GetNCentrBin()*GetNZvertBin()*GetNRPBin()] ;
- for(Int_t ic=0; ic<fNCentrBin; ic++){
+ for(Int_t ic=0; ic<GetNCentrBin(); ic++){
for(Int_t iz=0; iz<GetNZvertBin(); iz++){
for(Int_t irp=0; irp<GetNRPBin(); irp++){
fEventsList[ic*GetNZvertBin()*GetNRPBin()+iz*GetNRPBin()+irp] = new TList() ;
}
- fhRe1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
- fhMi1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
+ fhRe1 = new TH2D*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
+ fhMi1 = new TH2D*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
if(fFillBadDistHisto){
- fhRe2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
- fhRe3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
- fhMi2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
- fhMi3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
+ fhRe2 = new TH2D*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
+ fhRe3 = new TH2D*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
+ fhMi2 = new TH2D*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
+ fhMi3 = new TH2D*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
}
if(fMakeInvPtPlots) {
- fhReInvPt1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
- fhMiInvPt1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
+ fhReInvPt1 = new TH2D*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
+ fhMiInvPt1 = new TH2D*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
if(fFillBadDistHisto){
- fhReInvPt2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
- fhReInvPt3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
- fhMiInvPt2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
- fhMiInvPt3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
+ fhReInvPt2 = new TH2D*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
+ fhReInvPt3 = new TH2D*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
+ fhMiInvPt2 = new TH2D*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
+ fhMiInvPt3 = new TH2D*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
}
}
Int_t ntrmmax = GetHistoTrackMultiplicityMax();
Int_t ntrmmin = GetHistoTrackMultiplicityMin();
- fhAverTotECluster = new TH1F("hAverTotECluster","hAverTotECluster",200,0,50) ;
- fhAverTotECluster->SetXTitle("E_{cluster, aver. SM} (GeV)");
- outputContainer->Add(fhAverTotECluster) ;
-
- fhAverTotECell = new TH1F("hAverTotECell","hAverTotECell",200,0,50) ;
- fhAverTotECell->SetXTitle("E_{cell, aver. SM} (GeV)");
- outputContainer->Add(fhAverTotECell) ;
-
- fhAverTotECellvsCluster = new TH2F("hAverTotECellvsCluster","hAverTotECellvsCluster",200,0,50,200,0,50) ;
- fhAverTotECellvsCluster->SetYTitle("E_{cell, aver. SM} (GeV)");
- fhAverTotECellvsCluster->SetXTitle("E_{cluster, aver. SM} (GeV)");
- outputContainer->Add(fhAverTotECellvsCluster) ;
-
- fhEDensityCluster = new TH1F("hEDensityCluster","hEDensityCluster",200,0,50) ;
- fhEDensityCluster->SetXTitle("#Sigma E_{cluster} / N_{cluster} (GeV)");
- outputContainer->Add(fhEDensityCluster) ;
-
- fhEDensityCell = new TH1F("hEDensityCell","hEDensityCell",200,0,50) ;
- fhEDensityCell->SetXTitle("#Sigma E_{cell} / N_{cell} (GeV)");
- outputContainer->Add(fhEDensityCell) ;
-
- fhEDensityCellvsCluster = new TH2F("hEDensityCellvsCluster","hEDensityCellvsCluster",200,0,50,200,0,50) ;
- fhEDensityCellvsCluster->SetYTitle("#Sigma E_{cell} / N_{cell} (GeV)");
- fhEDensityCellvsCluster->SetXTitle("#Sigma E_{cluster} / N_{cluster} (GeV)");
- outputContainer->Add(fhEDensityCellvsCluster) ;
-
-// fhClusterPairDist = new TH1F("hClusterPairDist","Distance between clusters",250,0,750) ;
-// fhClusterPairDist->SetXTitle("#sqrt{(x_{1}-x_{2})^2+(z_{1}-z_{2})^2} (cm)");
-// outputContainer->Add(fhClusterPairDist) ;
-//
-// fhClusterPairDistWeight = new TH1F("hClusterPairDistWeighted","Distance between clusters, weighted by pair energy",200,0,400) ;
-// fhClusterPairDistWeight->SetXTitle("#sqrt{(x_{1}E_{1}-x_{2}E_{2})^{2}+(z_{1}E_{1}-z_{2}E_{2})^{2}}/ (E_{1}+E_{2}) (cm)");
-// outputContainer->Add(fhClusterPairDistWeight) ;
-//
-// fhAverClusterPairDist = new TH1F("hAverClusterPairDist","Average distance between clusters",250,0,750) ;
-// fhAverClusterPairDist->SetXTitle("#Sigma (#sqrt{(x_{1}-x_{2})^{2}+(z_{1}-z_{2})^{2}}) / N_{pairs} (cm)");
-// outputContainer->Add(fhAverClusterPairDist) ;
-//
-// fhAverClusterPairDistWeight = new TH1F("hAverClusterPairDistWeighted","Average distance between clusters, weighted by pair energy",200,0,400) ;
-// fhAverClusterPairDistWeight->SetXTitle("#Sigma (#sqrt{(x_{1}E_{1}-x_{2}E_{2})^{2}+(z_{1}E_{1}-z_{2}E_{2})^{2}}/ (E_{1}+E_{2})) / N_{pairs} (cm)");
-// outputContainer->Add(fhAverClusterPairDistWeight) ;
-//
-// fhAverClusterPairDistvsAverE = new TH2F("hAverClusterPairDistvsAverE","Average distance between clusters",250,0,750,200,0,50) ;
-// fhAverClusterPairDistvsAverE->SetXTitle("#Sigma (#sqrt{(x_{1}-x_{2})^{2}+(z_{1}-z_{2})^{2}}) / N_{pairs} (cm)");
-// fhAverClusterPairDistvsAverE->SetYTitle("#Sigma E_{cluster} / N_{cluster} (GeV)");
-// outputContainer->Add(fhAverClusterPairDistvsAverE) ;
-//
-// fhAverClusterPairDistWeightvsAverE = new TH2F("hAverClusterPairDistWeightedvsAverE","Average distance between clusters, weighted by pair energy",200,0,400,200,0,50) ;
-// fhAverClusterPairDistWeightvsAverE->SetXTitle("#Sigma (#sqrt{(x_{1}E_{1}-x_{2}E_{2})^2+(z_{1}E_{1}-z_{2}E_{2})^2}/ (E_{1}+E_{2})) / N_{pairs} (cm/GeV)");
-// fhAverClusterPairDistWeightvsAverE->SetYTitle("#Sigma E_{cluster} / N_{cluster} (GeV)");
-// outputContainer->Add(fhAverClusterPairDistWeightvsAverE) ;
-
-// fhAverClusterPairDistvsN = new TH2F("hAverClusterPairDistvsN","Average distance between clusters",250,0,750,200,0,50) ;
-// fhAverClusterPairDistvsN->SetXTitle("#Sigma (#sqrt{(x_{1}-x_{2})^{2}+(z_{1}-z_{2})^{2}}) / N_{pairs} (cm)");
-// fhAverClusterPairDistvsN->SetYTitle("N_{cluster}");
-// outputContainer->Add(fhAverClusterPairDistvsN) ;
-//
-// fhAverClusterPairDistWeightvsN = new TH2F("hAverClusterPairDistWeightedvsN","Average distance between clusters, weighted by pair energy",200,0,400,200,0,50) ;
-// fhAverClusterPairDistWeightvsN->SetXTitle("#Sigma (#sqrt{(x_{1}E_{1}-x_{2}E_{2})^{2}+(z_{1}E_{1}-z_{2}E_{2})^{2}}/ (E_{1}+E_{2})) / N_{pairs} (cm)");
-// fhAverClusterPairDistWeightvsN->SetYTitle("N_{cluster}");
-// outputContainer->Add(fhAverClusterPairDistWeightvsN) ;
-
-// fhMaxEvsClustMult = new TH2F("hMaxEvsClustMult","",nptbins,ptmin,ptmax,50,0,50) ;
-// fhMaxEvsClustMult->SetXTitle("E_{max}");
-// fhMaxEvsClustMult->SetYTitle("N_{cluster}");
-// outputContainer->Add(fhMaxEvsClustMult) ;
-//
-// fhMaxEvsClustEDen = new TH2F("hMaxEvsClustEDen","",nptbins,ptmin,ptmax,200,0,50) ;
-// fhMaxEvsClustEDen->SetXTitle("E_{max}");
-// fhMaxEvsClustEDen->SetYTitle("#Sigma E_{cluster} / N_{cluster} (GeV)");
-// outputContainer->Add(fhMaxEvsClustEDen) ;
-
- fhReConv = new TH2D("hReConv","Real Pair with one recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
- fhReConv->SetXTitle("p_{T} (GeV/c)");
- fhReConv->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
- outputContainer->Add(fhReConv) ;
-
- fhReConv2 = new TH2D("hReConv2","Real Pair with 2 recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
- fhReConv2->SetXTitle("p_{T} (GeV/c)");
- fhReConv2->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
- outputContainer->Add(fhReConv2) ;
-
- if(fDoOwnMix){
- fhMiConv = new TH2D("hMiConv","Mixed Pair with one recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
- fhMiConv->SetXTitle("p_{T} (GeV/c)");
- fhMiConv->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
- outputContainer->Add(fhMiConv) ;
-
- fhMiConv2 = new TH2D("hMiConv2","Mixed Pair with 2 recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
- fhMiConv2->SetXTitle("p_{T} (GeV/c)");
- fhMiConv2->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
- outputContainer->Add(fhMiConv2) ;
+ if(GetNCentrBin() > 1 && (fUseAverCellEBins||fUseAverClusterEBins||fUseAverClusterEDenBins)){
+
+ fhAverTotECluster = new TH1F("hAverTotECluster","hAverTotECluster",200,0,50) ;
+ fhAverTotECluster->SetXTitle("E_{cluster, aver. SM} (GeV)");
+ outputContainer->Add(fhAverTotECluster) ;
+
+ fhAverTotECell = new TH1F("hAverTotECell","hAverTotECell",200,0,50) ;
+ fhAverTotECell->SetXTitle("E_{cell, aver. SM} (GeV)");
+ outputContainer->Add(fhAverTotECell) ;
+
+ fhAverTotECellvsCluster = new TH2F("hAverTotECellvsCluster","hAverTotECellvsCluster",200,0,50,200,0,50) ;
+ fhAverTotECellvsCluster->SetYTitle("E_{cell, aver. SM} (GeV)");
+ fhAverTotECellvsCluster->SetXTitle("E_{cluster, aver. SM} (GeV)");
+ outputContainer->Add(fhAverTotECellvsCluster) ;
+
+ fhEDensityCluster = new TH1F("hEDensityCluster","hEDensityCluster",200,0,50) ;
+ fhEDensityCluster->SetXTitle("#Sigma E_{cluster} / N_{cluster} (GeV)");
+ outputContainer->Add(fhEDensityCluster) ;
+
+ fhEDensityCell = new TH1F("hEDensityCell","hEDensityCell",200,0,50) ;
+ fhEDensityCell->SetXTitle("#Sigma E_{cell} / N_{cell} (GeV)");
+ outputContainer->Add(fhEDensityCell) ;
+
+ fhEDensityCellvsCluster = new TH2F("hEDensityCellvsCluster","hEDensityCellvsCluster",200,0,50,200,0,50) ;
+ fhEDensityCellvsCluster->SetYTitle("#Sigma E_{cell} / N_{cell} (GeV)");
+ fhEDensityCellvsCluster->SetXTitle("#Sigma E_{cluster} / N_{cluster} (GeV)");
+ outputContainer->Add(fhEDensityCellvsCluster) ;
+
+ // fhClusterPairDist = new TH1F("hClusterPairDist","Distance between clusters",250,0,750) ;
+ // fhClusterPairDist->SetXTitle("#sqrt{(x_{1}-x_{2})^2+(z_{1}-z_{2})^2} (cm)");
+ // outputContainer->Add(fhClusterPairDist) ;
+ //
+ // fhClusterPairDistWeight = new TH1F("hClusterPairDistWeighted","Distance between clusters, weighted by pair energy",200,0,400) ;
+ // fhClusterPairDistWeight->SetXTitle("#sqrt{(x_{1}E_{1}-x_{2}E_{2})^{2}+(z_{1}E_{1}-z_{2}E_{2})^{2}}/ (E_{1}+E_{2}) (cm)");
+ // outputContainer->Add(fhClusterPairDistWeight) ;
+ //
+ // fhAverClusterPairDist = new TH1F("hAverClusterPairDist","Average distance between clusters",250,0,750) ;
+ // fhAverClusterPairDist->SetXTitle("#Sigma (#sqrt{(x_{1}-x_{2})^{2}+(z_{1}-z_{2})^{2}}) / N_{pairs} (cm)");
+ // outputContainer->Add(fhAverClusterPairDist) ;
+ //
+ // fhAverClusterPairDistWeight = new TH1F("hAverClusterPairDistWeighted","Average distance between clusters, weighted by pair energy",200,0,400) ;
+ // fhAverClusterPairDistWeight->SetXTitle("#Sigma (#sqrt{(x_{1}E_{1}-x_{2}E_{2})^{2}+(z_{1}E_{1}-z_{2}E_{2})^{2}}/ (E_{1}+E_{2})) / N_{pairs} (cm)");
+ // outputContainer->Add(fhAverClusterPairDistWeight) ;
+ //
+ // fhAverClusterPairDistvsAverE = new TH2F("hAverClusterPairDistvsAverE","Average distance between clusters",250,0,750,200,0,50) ;
+ // fhAverClusterPairDistvsAverE->SetXTitle("#Sigma (#sqrt{(x_{1}-x_{2})^{2}+(z_{1}-z_{2})^{2}}) / N_{pairs} (cm)");
+ // fhAverClusterPairDistvsAverE->SetYTitle("#Sigma E_{cluster} / N_{cluster} (GeV)");
+ // outputContainer->Add(fhAverClusterPairDistvsAverE) ;
+ //
+ // fhAverClusterPairDistWeightvsAverE = new TH2F("hAverClusterPairDistWeightedvsAverE","Average distance between clusters, weighted by pair energy",200,0,400,200,0,50) ;
+ // fhAverClusterPairDistWeightvsAverE->SetXTitle("#Sigma (#sqrt{(x_{1}E_{1}-x_{2}E_{2})^2+(z_{1}E_{1}-z_{2}E_{2})^2}/ (E_{1}+E_{2})) / N_{pairs} (cm/GeV)");
+ // fhAverClusterPairDistWeightvsAverE->SetYTitle("#Sigma E_{cluster} / N_{cluster} (GeV)");
+ // outputContainer->Add(fhAverClusterPairDistWeightvsAverE) ;
+
+ // fhAverClusterPairDistvsN = new TH2F("hAverClusterPairDistvsN","Average distance between clusters",250,0,750,200,0,50) ;
+ // fhAverClusterPairDistvsN->SetXTitle("#Sigma (#sqrt{(x_{1}-x_{2})^{2}+(z_{1}-z_{2})^{2}}) / N_{pairs} (cm)");
+ // fhAverClusterPairDistvsN->SetYTitle("N_{cluster}");
+ // outputContainer->Add(fhAverClusterPairDistvsN) ;
+ //
+ // fhAverClusterPairDistWeightvsN = new TH2F("hAverClusterPairDistWeightedvsN","Average distance between clusters, weighted by pair energy",200,0,400,200,0,50) ;
+ // fhAverClusterPairDistWeightvsN->SetXTitle("#Sigma (#sqrt{(x_{1}E_{1}-x_{2}E_{2})^{2}+(z_{1}E_{1}-z_{2}E_{2})^{2}}/ (E_{1}+E_{2})) / N_{pairs} (cm)");
+ // fhAverClusterPairDistWeightvsN->SetYTitle("N_{cluster}");
+ // outputContainer->Add(fhAverClusterPairDistWeightvsN) ;
+
+ // fhMaxEvsClustMult = new TH2F("hMaxEvsClustMult","",nptbins,ptmin,ptmax,50,0,50) ;
+ // fhMaxEvsClustMult->SetXTitle("E_{max}");
+ // fhMaxEvsClustMult->SetYTitle("N_{cluster}");
+ // outputContainer->Add(fhMaxEvsClustMult) ;
+ //
+ // fhMaxEvsClustEDen = new TH2F("hMaxEvsClustEDen","",nptbins,ptmin,ptmax,200,0,50) ;
+ // fhMaxEvsClustEDen->SetXTitle("E_{max}");
+ // fhMaxEvsClustEDen->SetYTitle("#Sigma E_{cluster} / N_{cluster} (GeV)");
+ // outputContainer->Add(fhMaxEvsClustEDen) ;
+ }//counting and average histograms
+
+ if(fCheckConversion){
+ fhReConv = new TH2D("hReConv","Real Pair with one recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+ fhReConv->SetXTitle("p_{T} (GeV/c)");
+ fhReConv->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+ outputContainer->Add(fhReConv) ;
+
+ fhReConv2 = new TH2D("hReConv2","Real Pair with 2 recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+ fhReConv2->SetXTitle("p_{T} (GeV/c)");
+ fhReConv2->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+ outputContainer->Add(fhReConv2) ;
+
+ if(fDoOwnMix){
+ fhMiConv = new TH2D("hMiConv","Mixed Pair with one recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+ fhMiConv->SetXTitle("p_{T} (GeV/c)");
+ fhMiConv->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+ outputContainer->Add(fhMiConv) ;
+
+ fhMiConv2 = new TH2D("hMiConv2","Mixed Pair with 2 recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+ fhMiConv2->SetXTitle("p_{T} (GeV/c)");
+ fhMiConv2->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+ outputContainer->Add(fhMiConv2) ;
+ }
}
- for(Int_t ic=0; ic<fNCentrBin; ic++){
+ for(Int_t ic=0; ic<GetNCentrBin(); ic++){
for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
Int_t index = ((ic*fNPIDBits)+ipid)*fNAsymCuts + iasym;
}// pid bit loop
fhRePtNCellAsymCuts = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
- fhRePtNCellAsymCutsSM0 = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
- fhRePtNCellAsymCutsSM1 = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
- fhRePtNCellAsymCutsSM2 = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
- fhRePtNCellAsymCutsSM3 = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
fhMiPtNCellAsymCuts = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
+
+ if(fFillSMCombinations){
+ fhRePtNCellAsymCutsSM0 = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
+ fhRePtNCellAsymCutsSM1 = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
+ fhRePtNCellAsymCutsSM2 = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
+ fhRePtNCellAsymCutsSM3 = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
+ }
+
for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
for(Int_t icell=0; icell<fNCellNCuts; icell++){
for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
fhRePtNCellAsymCuts[index]->SetXTitle("p_{T} (GeV/c)");
fhRePtNCellAsymCuts[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
outputContainer->Add(fhRePtNCellAsymCuts[index]) ;
-
- snprintf(key, buffersize,"hRe_pt%d_cell%d_asym%d_SM0",ipt,icell,iasym) ;
- snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f, SM 0 ",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
- fhRePtNCellAsymCutsSM0[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
- fhRePtNCellAsymCutsSM0[index]->SetXTitle("p_{T} (GeV/c)");
- fhRePtNCellAsymCutsSM0[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
- outputContainer->Add(fhRePtNCellAsymCutsSM0[index]) ;
-
- snprintf(key, buffersize,"hRe_pt%d_cell%d_asym%d_SM1",ipt,icell,iasym) ;
- snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f, SM 1 ",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
- fhRePtNCellAsymCutsSM1[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
- fhRePtNCellAsymCutsSM1[index]->SetXTitle("p_{T} (GeV/c)");
- fhRePtNCellAsymCutsSM1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
- outputContainer->Add(fhRePtNCellAsymCutsSM1[index]) ;
-
- snprintf(key, buffersize,"hRe_pt%d_cell%d_asym%d_SM2",ipt,icell,iasym) ;
- snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f, SM 2 ",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
- fhRePtNCellAsymCutsSM2[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
- fhRePtNCellAsymCutsSM2[index]->SetXTitle("p_{T} (GeV/c)");
- fhRePtNCellAsymCutsSM2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
- outputContainer->Add(fhRePtNCellAsymCutsSM2[index]) ;
-
- snprintf(key, buffersize,"hRe_pt%d_cell%d_asym%d_SM3",ipt,icell,iasym) ;
- snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f, SM 3 ",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
- fhRePtNCellAsymCutsSM3[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
- fhRePtNCellAsymCutsSM3[index]->SetXTitle("p_{T} (GeV/c)");
- fhRePtNCellAsymCutsSM3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
- outputContainer->Add(fhRePtNCellAsymCutsSM3[index]) ;
snprintf(key, buffersize,"hMi_pt%d_cell%d_asym%d",ipt,icell,iasym) ;
snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
fhMiPtNCellAsymCuts[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
fhMiPtNCellAsymCuts[index]->SetXTitle("p_{T} (GeV/c)");
fhMiPtNCellAsymCuts[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
- outputContainer->Add(fhMiPtNCellAsymCuts[index]) ;
+ outputContainer->Add(fhMiPtNCellAsymCuts[index]) ;
+
+ if(fFillSMCombinations){
+ snprintf(key, buffersize,"hRe_pt%d_cell%d_asym%d_SM0",ipt,icell,iasym) ;
+ snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f, SM 0 ",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
+ fhRePtNCellAsymCutsSM0[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+ fhRePtNCellAsymCutsSM0[index]->SetXTitle("p_{T} (GeV/c)");
+ fhRePtNCellAsymCutsSM0[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+ outputContainer->Add(fhRePtNCellAsymCutsSM0[index]) ;
+
+ snprintf(key, buffersize,"hRe_pt%d_cell%d_asym%d_SM1",ipt,icell,iasym) ;
+ snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f, SM 1 ",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
+ fhRePtNCellAsymCutsSM1[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+ fhRePtNCellAsymCutsSM1[index]->SetXTitle("p_{T} (GeV/c)");
+ fhRePtNCellAsymCutsSM1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+ outputContainer->Add(fhRePtNCellAsymCutsSM1[index]) ;
+
+ snprintf(key, buffersize,"hRe_pt%d_cell%d_asym%d_SM2",ipt,icell,iasym) ;
+ snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f, SM 2 ",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
+ fhRePtNCellAsymCutsSM2[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+ fhRePtNCellAsymCutsSM2[index]->SetXTitle("p_{T} (GeV/c)");
+ fhRePtNCellAsymCutsSM2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+ outputContainer->Add(fhRePtNCellAsymCutsSM2[index]) ;
+
+ snprintf(key, buffersize,"hRe_pt%d_cell%d_asym%d_SM3",ipt,icell,iasym) ;
+ snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f, SM 3 ",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
+ fhRePtNCellAsymCutsSM3[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+ fhRePtNCellAsymCutsSM3[index]->SetXTitle("p_{T} (GeV/c)");
+ fhRePtNCellAsymCutsSM3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+ outputContainer->Add(fhRePtNCellAsymCutsSM3[index]) ;
+ }
}
}
}
- fhRePtMult = new TH3D*[fNAsymCuts] ;
- for(Int_t iasym = 0; iasym<fNAsymCuts; iasym++){
- fhRePtMult[iasym] = new TH3D(Form("hRePtMult_asym%d",iasym),Form("(p_{T},C,M)_{#gamma#gamma}, A<%1.2f",fAsymCuts[iasym]),
- nptbins,ptmin,ptmax,ntrmbins,ntrmmin,ntrmmax,nmassbins,massmin,massmax);
- fhRePtMult[iasym]->SetXTitle("p_{T} (GeV/c)");
- fhRePtMult[iasym]->SetYTitle("Track multiplicity");
- fhRePtMult[iasym]->SetZTitle("m_{#gamma,#gamma} (GeV/c^{2})");
- outputContainer->Add(fhRePtMult[iasym]) ;
+ if(ntrmbins!=0){
+ fhRePtMult = new TH3D*[fNAsymCuts] ;
+ for(Int_t iasym = 0; iasym<fNAsymCuts; iasym++){
+ fhRePtMult[iasym] = new TH3D(Form("hRePtMult_asym%d",iasym),Form("(p_{T},C,M)_{#gamma#gamma}, A<%1.2f",fAsymCuts[iasym]),
+ nptbins,ptmin,ptmax,ntrmbins,ntrmmin,ntrmmax,nmassbins,massmin,massmax);
+ fhRePtMult[iasym]->SetXTitle("p_{T} (GeV/c)");
+ fhRePtMult[iasym]->SetYTitle("Track multiplicity");
+ fhRePtMult[iasym]->SetZTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+ outputContainer->Add(fhRePtMult[iasym]) ;
+ }
}
-
}// multi cuts analysis
- fhEvents=new TH3D("hEvents","Number of events",fNCentrBin,0.,1.*fNCentrBin,
+ fhEvents=new TH3D("hEvents","Number of events",GetNCentrBin(),0.,1.*GetNCentrBin(),
GetNZvertBin(),0.,1.*GetNZvertBin(),GetNRPBin(),0.,1.*GetNRPBin()) ;
fhEvents->SetXTitle("Centrality bin");
fhEvents->SetZTitle("RP bin");
outputContainer->Add(fhEvents) ;
- fhCentrality=new TH1D("hCentralityBin","Number of events in centrality bin",fNCentrBin*10,0.,1.*fNCentrBin) ;
- fhCentrality->SetXTitle("Centrality bin");
- outputContainer->Add(fhCentrality) ;
-
+ if(GetNCentrBin()>1){
+ fhCentrality=new TH1D("hCentralityBin","Number of events in centrality bin",GetNCentrBin(),0.,1.*GetNCentrBin()) ;
+ fhCentrality->SetXTitle("Centrality bin");
+ outputContainer->Add(fhCentrality) ;
+
+ fhCentralityNoPair=new TH1D("hCentralityBinNoPair","Number of events in centrality bin, with no cluster pairs",GetNCentrBin(),0.,1.*GetNCentrBin()) ;
+ fhCentralityNoPair->SetXTitle("Centrality bin");
+ outputContainer->Add(fhCentralityNoPair) ;
+ }
+
+ if(GetNRPBin() > 1 ){
+
+ fhEventPlaneAngle=new TH1D("hEventPlaneAngleBin","Number of events in centrality bin",100,0.,TMath::TwoPi()) ;
+ fhEventPlaneAngle->SetXTitle("EP angle (rad)");
+ outputContainer->Add(fhEventPlaneAngle) ;
+
+ if(GetNCentrBin()>1){
+ fhEventPlaneResolution=new TH2D("hEventPlaneResolution","Event plane resolution",GetNCentrBin(),0,GetNCentrBin(),100,0.,TMath::TwoPi()) ;
+ fhEventPlaneResolution->SetYTitle("Resolution");
+ fhEventPlaneResolution->SetXTitle("Centrality Bin");
+ outputContainer->Add(fhEventPlaneResolution) ;
+ }
+ }
+
fhRealOpeningAngle = new TH2D
("hRealOpeningAngle","Angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,300,0,TMath::Pi());
fhRealOpeningAngle->SetYTitle("#theta(rad)");
//Histograms filled only if MC data is requested
if(IsDataMC()){
//Pi0
- fhPrimPi0Pt = new TH1D("hPrimPi0Pt","Primary pi0 pt",nptbins,ptmin,ptmax) ;
+ fhPrimPi0Pt = new TH1D("hPrimPi0Pt","Primary pi0 pt, Y<1",nptbins,ptmin,ptmax) ;
fhPrimPi0AccPt = new TH1D("hPrimPi0AccPt","Primary pi0 pt with both photons in acceptance",nptbins,ptmin,ptmax) ;
fhPrimPi0Pt ->SetXTitle("p_{T} (GeV/c)");
fhPrimPi0AccPt->SetXTitle("p_{T} (GeV/c)");
outputContainer->Add(fhPrimPi0Pt) ;
outputContainer->Add(fhPrimPi0AccPt) ;
- fhPrimPi0Y = new TH2D("hPrimPi0Rapidity","Rapidity of primary pi0",nptbins,ptmin,ptmax, netabins,etamin,etamax) ;
+ Int_t netabinsopen = TMath::Nint(netabins*4/(etamax-etamin));
+ fhPrimPi0Y = new TH2D("hPrimPi0Rapidity","Rapidity of primary pi0",nptbins,ptmin,ptmax,netabinsopen,-2, 2) ;
fhPrimPi0Y ->SetYTitle("Rapidity");
fhPrimPi0Y ->SetXTitle("p_{T} (GeV/c)");
outputContainer->Add(fhPrimPi0Y) ;
- fhPrimPi0AccY = new TH2D("hPrimPi0AccRapidity","Rapidity of primary pi0",nptbins,ptmin,ptmax, netabins,etamin,etamax) ;
+ fhPrimPi0AccY = new TH2D("hPrimPi0AccRapidity","Rapidity of primary pi0",nptbins,ptmin,ptmax,netabins,etamin,etamax) ;
fhPrimPi0AccY->SetYTitle("Rapidity");
fhPrimPi0AccY->SetXTitle("p_{T} (GeV/c)");
outputContainer->Add(fhPrimPi0AccY) ;
- fhPrimPi0Phi = new TH2D("hPrimPi0Phi","Azimuthal of primary pi0",nptbins,ptmin,ptmax, nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
+ Int_t nphibinsopen = TMath::Nint(nphibins*TMath::TwoPi()/(phimax-phimin));
+ fhPrimPi0Phi = new TH2D("hPrimPi0Phi","Azimuthal of primary pi0, Y<1",nptbins,ptmin,ptmax,nphibinsopen,0,360) ;
fhPrimPi0Phi->SetYTitle("#phi (deg)");
fhPrimPi0Phi->SetXTitle("p_{T} (GeV/c)");
outputContainer->Add(fhPrimPi0Phi) ;
- fhPrimPi0AccPhi = new TH2D("hPrimPi0AccPhi","Azimuthal of primary pi0 with accepted daughters",nptbins,ptmin,ptmax, nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
+ fhPrimPi0AccPhi = new TH2D("hPrimPi0AccPhi","Azimuthal of primary pi0 with accepted daughters",nptbins,ptmin,ptmax,
+ nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
fhPrimPi0AccPhi->SetYTitle("#phi (deg)");
fhPrimPi0AccPhi->SetXTitle("p_{T} (GeV/c)");
outputContainer->Add(fhPrimPi0AccPhi) ;
}
}
- TString pairnamePHOS[] = {"(0-1)","(0-2)","(1-2)","(0-3)","(0-4)","(1-3)","(1-4)","(2-3)","(2-4)","(3-4)"};
- for(Int_t imod=0; imod<fNModules; imod++){
- //Module dependent invariant mass
- snprintf(key, buffersize,"hReMod_%d",imod) ;
- snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for Module %d",imod) ;
- fhReMod[imod] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
- fhReMod[imod]->SetXTitle("p_{T} (GeV/c)");
- fhReMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
- outputContainer->Add(fhReMod[imod]) ;
- if(fCalorimeter=="PHOS"){
- snprintf(key, buffersize,"hReDiffPHOSMod_%d",imod) ;
- snprintf(title, buffersize,"Real pairs PHOS, clusters in different Modules: %s",(pairnamePHOS[imod]).Data()) ;
- fhReDiffPHOSMod[imod] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
- fhReDiffPHOSMod[imod]->SetXTitle("p_{T} (GeV/c)");
- fhReDiffPHOSMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
- outputContainer->Add(fhReDiffPHOSMod[imod]) ;
- }
- else{//EMCAL
- if(imod<fNModules/2){
- snprintf(key, buffersize,"hReSameSectorEMCAL_%d",imod) ;
- snprintf(title, buffersize,"Real pairs EMCAL, clusters in same sector, SM(%d,%d)",imod*2,imod*2+1) ;
- fhReSameSectorEMCALMod[imod] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
- fhReSameSectorEMCALMod[imod]->SetXTitle("p_{T} (GeV/c)");
- fhReSameSectorEMCALMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
- outputContainer->Add(fhReSameSectorEMCALMod[imod]) ;
- }
- if(imod<fNModules-2){
- snprintf(key, buffersize,"hReSameSideEMCAL_%d",imod) ;
- snprintf(title, buffersize,"Real pairs EMCAL, clusters in same side SM(%d,%d)",imod, imod+2) ;
- fhReSameSideEMCALMod[imod] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
- fhReSameSideEMCALMod[imod]->SetXTitle("p_{T} (GeV/c)");
- fhReSameSideEMCALMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
- outputContainer->Add(fhReSameSideEMCALMod[imod]) ;
- }
- }//EMCAL
-
- if(fDoOwnMix){
- snprintf(key, buffersize,"hMiMod_%d",imod) ;
- snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for Module %d",imod) ;
- fhMiMod[imod] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
- fhMiMod[imod]->SetXTitle("p_{T} (GeV/c)");
- fhMiMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
- outputContainer->Add(fhMiMod[imod]) ;
-
+ if(fFillSMCombinations){
+ TString pairnamePHOS[] = {"(0-1)","(0-2)","(1-2)","(0-3)","(0-4)","(1-3)","(1-4)","(2-3)","(2-4)","(3-4)"};
+ for(Int_t imod=0; imod<fNModules; imod++){
+ //Module dependent invariant mass
+ snprintf(key, buffersize,"hReMod_%d",imod) ;
+ snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for Module %d",imod) ;
+ fhReMod[imod] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+ fhReMod[imod]->SetXTitle("p_{T} (GeV/c)");
+ fhReMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+ outputContainer->Add(fhReMod[imod]) ;
if(fCalorimeter=="PHOS"){
- snprintf(key, buffersize,"hMiDiffPHOSMod_%d",imod) ;
- snprintf(title, buffersize,"Mixed pairs PHOS, clusters in different Modules: %s",(pairnamePHOS[imod]).Data()) ;
- fhMiDiffPHOSMod[imod] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
- fhMiDiffPHOSMod[imod]->SetXTitle("p_{T} (GeV/c)");
- fhMiDiffPHOSMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
- outputContainer->Add(fhMiDiffPHOSMod[imod]) ;
- }//PHOS
+ snprintf(key, buffersize,"hReDiffPHOSMod_%d",imod) ;
+ snprintf(title, buffersize,"Real pairs PHOS, clusters in different Modules: %s",(pairnamePHOS[imod]).Data()) ;
+ fhReDiffPHOSMod[imod] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+ fhReDiffPHOSMod[imod]->SetXTitle("p_{T} (GeV/c)");
+ fhReDiffPHOSMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+ outputContainer->Add(fhReDiffPHOSMod[imod]) ;
+ }
else{//EMCAL
if(imod<fNModules/2){
- snprintf(key, buffersize,"hMiSameSectorEMCALMod_%d",imod) ;
- snprintf(title, buffersize,"Mixed pairs EMCAL, clusters in same sector, SM(%d,%d)",imod*2,imod*2+1) ;
- fhMiSameSectorEMCALMod[imod] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
- fhMiSameSectorEMCALMod[imod]->SetXTitle("p_{T} (GeV/c)");
- fhMiSameSectorEMCALMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
- outputContainer->Add(fhMiSameSectorEMCALMod[imod]) ;
+ snprintf(key, buffersize,"hReSameSectorEMCAL_%d",imod) ;
+ snprintf(title, buffersize,"Real pairs EMCAL, clusters in same sector, SM(%d,%d)",imod*2,imod*2+1) ;
+ fhReSameSectorEMCALMod[imod] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+ fhReSameSectorEMCALMod[imod]->SetXTitle("p_{T} (GeV/c)");
+ fhReSameSectorEMCALMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+ outputContainer->Add(fhReSameSectorEMCALMod[imod]) ;
}
if(imod<fNModules-2){
- snprintf(key, buffersize,"hMiSameSideEMCALMod_%d",imod) ;
- snprintf(title, buffersize,"Mixed pairs EMCAL, clusters in same side SM(%d,%d)",imod, imod+2) ;
- fhMiSameSideEMCALMod[imod] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
- fhMiSameSideEMCALMod[imod]->SetXTitle("p_{T} (GeV/c)");
- fhMiSameSideEMCALMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
- outputContainer->Add(fhMiSameSideEMCALMod[imod]) ;
+ snprintf(key, buffersize,"hReSameSideEMCAL_%d",imod) ;
+ snprintf(title, buffersize,"Real pairs EMCAL, clusters in same side SM(%d,%d)",imod, imod+2) ;
+ fhReSameSideEMCALMod[imod] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+ fhReSameSideEMCALMod[imod]->SetXTitle("p_{T} (GeV/c)");
+ fhReSameSideEMCALMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+ outputContainer->Add(fhReSameSideEMCALMod[imod]) ;
}
- }//EMCAL
- }
- }
+ }//EMCAL
+
+ if(fDoOwnMix){
+ snprintf(key, buffersize,"hMiMod_%d",imod) ;
+ snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for Module %d",imod) ;
+ fhMiMod[imod] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+ fhMiMod[imod]->SetXTitle("p_{T} (GeV/c)");
+ fhMiMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+ outputContainer->Add(fhMiMod[imod]) ;
+
+ if(fCalorimeter=="PHOS"){
+ snprintf(key, buffersize,"hMiDiffPHOSMod_%d",imod) ;
+ snprintf(title, buffersize,"Mixed pairs PHOS, clusters in different Modules: %s",(pairnamePHOS[imod]).Data()) ;
+ fhMiDiffPHOSMod[imod] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+ fhMiDiffPHOSMod[imod]->SetXTitle("p_{T} (GeV/c)");
+ fhMiDiffPHOSMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+ outputContainer->Add(fhMiDiffPHOSMod[imod]) ;
+ }//PHOS
+ else{//EMCAL
+ if(imod<fNModules/2){
+ snprintf(key, buffersize,"hMiSameSectorEMCALMod_%d",imod) ;
+ snprintf(title, buffersize,"Mixed pairs EMCAL, clusters in same sector, SM(%d,%d)",imod*2,imod*2+1) ;
+ fhMiSameSectorEMCALMod[imod] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+ fhMiSameSectorEMCALMod[imod]->SetXTitle("p_{T} (GeV/c)");
+ fhMiSameSectorEMCALMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+ outputContainer->Add(fhMiSameSectorEMCALMod[imod]) ;
+ }
+ if(imod<fNModules-2){
+ snprintf(key, buffersize,"hMiSameSideEMCALMod_%d",imod) ;
+ snprintf(title, buffersize,"Mixed pairs EMCAL, clusters in same side SM(%d,%d)",imod, imod+2) ;
+ fhMiSameSideEMCALMod[imod] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+ fhMiSameSideEMCALMod[imod]->SetXTitle("p_{T} (GeV/c)");
+ fhMiSameSideEMCALMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+ outputContainer->Add(fhMiSameSideEMCALMod[imod]) ;
+ }
+ }//EMCAL
+ }// own mix
+ }//loop combinations
+ } // SM combinations
// for(Int_t i = 0; i < outputContainer->GetEntries() ; i++){
//
printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
AliAnaPartCorrBaseClass::Print(" ");
- printf("Number of bins in Centrality: %d \n",fNCentrBin) ;
+ printf("Number of bins in Centrality: %d \n",GetNCentrBin()) ;
printf("Number of bins in Z vert. pos: %d \n",GetNZvertBin()) ;
printf("Number of bins in Reac. Plain: %d \n",GetNRPBin()) ;
- printf("Depth of event buffer: %d \n",fNmaxMixEv) ;
+ printf("Depth of event buffer: %d \n",GetNMaxEvMix()) ;
printf("Pair in same Module: %d \n",fSameSM) ;
printf("Cuts: \n") ;
// printf("Z vertex position: -%2.3f < z < %2.3f \n",GetZvertexCut(),GetZvertexCut()) ; //It crashes here, why?
if(GetReader()->ReadStack()){
AliStack * stack = GetMCStack();
if(stack){
- for(Int_t i=0 ; i<stack->GetNprimary(); i++){
+ for(Int_t i=0 ; i<stack->GetNtrack(); i++){
TParticle * prim = stack->Particle(i) ;
Int_t pdg = prim->GetPdgCode();
+ //printf("i %d, %s %d %s %d \n",i, stack->Particle(i)->GetName(), stack->Particle(i)->GetPdgCode(),
+ // prim->GetName(), prim->GetPdgCode());
+
if( pdg == 111 || pdg == 221){
Double_t pi0Pt = prim->Pt() ;
- //printf("pi0, pt %2.2f\n",pi0Pt);
if(prim->Energy() == TMath::Abs(prim->Pz())) continue ; //Protection against floating point exception
Double_t pi0Y = 0.5*TMath::Log((prim->Energy()-prim->Pz())/(prim->Energy()+prim->Pz())) ;
Double_t phi = TMath::RadToDeg()*prim->Phi() ;
if(pdg == 111){
if(TMath::Abs(pi0Y) < 1.0){
- fhPrimPi0Pt->Fill(pi0Pt) ;
+ fhPrimPi0Pt ->Fill(pi0Pt) ;
+ fhPrimPi0Phi->Fill(pi0Pt, phi) ;
}
fhPrimPi0Y ->Fill(pi0Pt, pi0Y) ;
- fhPrimPi0Phi->Fill(pi0Pt, phi) ;
}
else if(pdg == 221){
if(TMath::Abs(pi0Y) < 1.0){
- fhPrimEtaPt->Fill(pi0Pt) ;
+ fhPrimEtaPt ->Fill(pi0Pt) ;
+ fhPrimEtaPhi->Fill(pi0Pt, phi) ;
}
fhPrimEtaY ->Fill(pi0Pt, pi0Y) ;
- fhPrimEtaPhi->Fill(pi0Pt, phi) ;
}
//Origin of meson
Int_t momindex = prim->GetFirstMother();
- if(momindex < 0) continue;
- TParticle* mother = stack->Particle(momindex);
- Int_t mompdg = TMath::Abs(mother->GetPdgCode());
- Int_t momstatus = mother->GetStatusCode();
- if(pdg == 111){
- if (momstatus == 21)fhPrimPi0PtOrigin->Fill(pi0Pt,0.5);//parton
- else if(mompdg < 22 ) fhPrimPi0PtOrigin->Fill(pi0Pt,1.5);//quark
- else if(mompdg > 2100 && mompdg < 2210) fhPrimPi0PtOrigin->Fill(pi0Pt,2.5);// resonances
- else if(mompdg == 221) fhPrimPi0PtOrigin->Fill(pi0Pt,8.5);//eta
- else if(mompdg == 331) fhPrimPi0PtOrigin->Fill(pi0Pt,9.5);//eta prime
- else if(mompdg == 213) fhPrimPi0PtOrigin->Fill(pi0Pt,4.5);//rho
- else if(mompdg == 223) fhPrimPi0PtOrigin->Fill(pi0Pt,5.5);//omega
- else if(mompdg >= 310 && mompdg <= 323) fhPrimPi0PtOrigin->Fill(pi0Pt,6.5);//k0S, k+-,k*
- else if(mompdg == 130) fhPrimPi0PtOrigin->Fill(pi0Pt,6.5);//k0L
- else if(momstatus == 11 || momstatus == 12 ) fhPrimPi0PtOrigin->Fill(pi0Pt,3.5);//resonances
- else fhPrimPi0PtOrigin->Fill(pi0Pt,7.5);//other?
- }//pi0
- else {
- if (momstatus == 21 ) fhPrimEtaPtOrigin->Fill(pi0Pt,0.5);//parton
- else if(mompdg < 22 ) fhPrimEtaPtOrigin->Fill(pi0Pt,1.5);//quark
- else if(mompdg > 2100 && mompdg < 2210) fhPrimEtaPtOrigin->Fill(pi0Pt,2.5);//qq resonances
- else if(mompdg == 331) fhPrimEtaPtOrigin->Fill(pi0Pt,5.5);//eta prime
- else if(momstatus == 11 || momstatus == 12 ) fhPrimEtaPtOrigin->Fill(pi0Pt,3.5);//resonances
- else fhPrimEtaPtOrigin->Fill(pi0Pt,4.5);//stable, conversions?
- //printf("Other Meson pdg %d, Mother %s, pdg %d, status %d\n",pdg, TDatabasePDG::Instance()->GetParticle(mompdg)->GetName(),mompdg, momstatus );
- }
-
+ if(momindex >= 0) {
+ TParticle* mother = stack->Particle(momindex);
+ Int_t mompdg = TMath::Abs(mother->GetPdgCode());
+ Int_t momstatus = mother->GetStatusCode();
+ if(pdg == 111){
+ if (momstatus == 21)fhPrimPi0PtOrigin->Fill(pi0Pt,0.5);//parton
+ else if(mompdg < 22 ) fhPrimPi0PtOrigin->Fill(pi0Pt,1.5);//quark
+ else if(mompdg > 2100 && mompdg < 2210) fhPrimPi0PtOrigin->Fill(pi0Pt,2.5);// resonances
+ else if(mompdg == 221) fhPrimPi0PtOrigin->Fill(pi0Pt,8.5);//eta
+ else if(mompdg == 331) fhPrimPi0PtOrigin->Fill(pi0Pt,9.5);//eta prime
+ else if(mompdg == 213) fhPrimPi0PtOrigin->Fill(pi0Pt,4.5);//rho
+ else if(mompdg == 223) fhPrimPi0PtOrigin->Fill(pi0Pt,5.5);//omega
+ else if(mompdg >= 310 && mompdg <= 323) fhPrimPi0PtOrigin->Fill(pi0Pt,6.5);//k0S, k+-,k*
+ else if(mompdg == 130) fhPrimPi0PtOrigin->Fill(pi0Pt,6.5);//k0L
+ else if(momstatus == 11 || momstatus == 12 ) fhPrimPi0PtOrigin->Fill(pi0Pt,3.5);//resonances
+ else fhPrimPi0PtOrigin->Fill(pi0Pt,7.5);//other?
+ }//pi0
+ else {
+ if (momstatus == 21 ) fhPrimEtaPtOrigin->Fill(pi0Pt,0.5);//parton
+ else if(mompdg < 22 ) fhPrimEtaPtOrigin->Fill(pi0Pt,1.5);//quark
+ else if(mompdg > 2100 && mompdg < 2210) fhPrimEtaPtOrigin->Fill(pi0Pt,2.5);//qq resonances
+ else if(mompdg == 331) fhPrimEtaPtOrigin->Fill(pi0Pt,5.5);//eta prime
+ else if(momstatus == 11 || momstatus == 12 ) fhPrimEtaPtOrigin->Fill(pi0Pt,3.5);//resonances
+ else fhPrimEtaPtOrigin->Fill(pi0Pt,4.5);//stable, conversions?
+ //printf("Other Meson pdg %d, Mother %s, pdg %d, status %d\n",pdg, TDatabasePDG::Instance()->GetParticle(mompdg)->GetName(),mompdg, momstatus );
+ }
+ } // pi0 has mother
//Check if both photons hit Calorimeter
- if(prim->GetNDaughters()!=2) return; //Only interested in 2 gamma decay
+ if(prim->GetNDaughters()!=2) continue; //Only interested in 2 gamma decay
Int_t iphot1=prim->GetFirstDaughter() ;
Int_t iphot2=prim->GetLastDaughter() ;
if(iphot1>-1 && iphot1<stack->GetNtrack() && iphot2>-1 && iphot2<stack->GetNtrack()){
TLorentzVector lv1, lv2;
phot1->Momentum(lv1);
phot2->Momentum(lv2);
-
Bool_t inacceptance = kFALSE;
if(fCalorimeter == "PHOS"){
if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet()){
}//stack exists and data is MC
}//read stack
else if(GetReader()->ReadAODMCParticles()){
-
TClonesArray * mcparticles = GetReader()->GetAODMCParticles(0);
if(mcparticles){
Int_t nprim = mcparticles->GetEntriesFast();
+
for(Int_t i=0; i < nprim; i++)
{
AliAODMCParticle * prim = (AliAODMCParticle *) mcparticles->At(i);
- // Only generator particles
- if( prim->GetStatus() == 0) break;
+
+ // Only generator particles, when they come from PYTHIA, PHOJET, HERWIG ...
+ //if( prim->GetStatus() == 0 && (GetMCAnalysisUtils()->GetMCGenerator()).Length()!=0) break;
Int_t pdg = prim->GetPdgCode();
if( pdg == 111 || pdg == 221){
Double_t pi0Pt = prim->Pt() ;
- //printf("pi0, pt %2.2f\n",pi0Pt);
- if(prim->E() == TMath::Abs(prim->Pz())) continue ; //Protection against floating point exception
+ //printf("pi0, pt %2.2f, eta %f, phi %f\n",pi0Pt, prim->Eta(), prim->Phi());
+ if(prim->E() == TMath::Abs(prim->Pz())) continue ; //Protection against floating point exception
+
Double_t pi0Y = 0.5*TMath::Log((prim->E()-prim->Pz())/(prim->E()+prim->Pz())) ;
Double_t phi = TMath::RadToDeg()*prim->Phi() ;
if(pdg == 111){
- if(TMath::Abs(pi0Y) < 0.5){
- fhPrimPi0Pt->Fill(pi0Pt) ;
+ if(TMath::Abs(pi0Y) < 1){
+ fhPrimPi0Pt->Fill(pi0Pt) ;
+ fhPrimPi0Phi->Fill(pi0Pt, phi) ;
}
fhPrimPi0Y ->Fill(pi0Pt, pi0Y) ;
- fhPrimPi0Phi->Fill(pi0Pt, phi) ;
}
else if(pdg == 221){
- if(TMath::Abs(pi0Y) < 0.5){
- fhPrimEtaPt->Fill(pi0Pt) ;
+ if(TMath::Abs(pi0Y) < 1){
+ fhPrimEtaPt->Fill(pi0Pt) ;
+ fhPrimEtaPhi->Fill(pi0Pt, phi) ;
}
fhPrimEtaY ->Fill(pi0Pt, pi0Y) ;
- fhPrimEtaPhi->Fill(pi0Pt, phi) ;
}
//Origin of meson
Int_t momindex = prim->GetMother();
- if(momindex < 0) continue;
- AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
- Int_t mompdg = TMath::Abs(mother->GetPdgCode());
- Int_t momstatus = mother->GetStatus();
+ if(momindex >= 0) {
+ AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
+ Int_t mompdg = TMath::Abs(mother->GetPdgCode());
+ Int_t momstatus = mother->GetStatus();
if(pdg == 111){
if (momstatus == 21) fhPrimPi0PtOrigin->Fill(pi0Pt,0.5);//parton
else if(mompdg < 22 ) fhPrimPi0PtOrigin->Fill(pi0Pt,1.5);//quark
else fhPrimEtaPtOrigin->Fill(pi0Pt,4.5);//stable, conversions?
//printf("Other Meson pdg %d, Mother %s, pdg %d, status %d\n",pdg, TDatabasePDG::Instance()->GetParticle(mompdg)->GetName(),mompdg, momstatus );
}
+ }//pi0 has mother
//Check if both photons hit Calorimeter
- if(prim->GetNDaughters()!=2) return; //Only interested in 2 gamma decay
+ if(prim->GetNDaughters()!=2) continue; //Only interested in 2 gamma decay
Int_t iphot1=prim->GetDaughter(0) ;
Int_t iphot2=prim->GetDaughter(1) ;
if(iphot1>-1 && iphot1<nprim && iphot2>-1 && iphot2<nprim){
Int_t ancPDG = 0;
Int_t ancStatus = 0;
TLorentzVector ancMomentum;
+ TVector3 prodVertex;
Int_t ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(index1, index2,
- GetReader(), ancPDG, ancStatus,ancMomentum);
+ GetReader(), ancPDG, ancStatus,ancMomentum, prodVertex);
Int_t momindex = -1;
Int_t mompdg = -1;
fhMCOrgDeltaEta[10]->Fill(pt,deta);
fhMCOrgDeltaPhi[10]->Fill(pt,dphi);
}
- else {
- printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Check ** Common ancestor label %d, pdg %d, name %s, status %d; \n",
- ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
- }
+ // else {
+ // printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Check ** Common ancestor label %d, pdg %d, name %s, status %d; \n",
+ // ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
+ // }
}//status 21
- else {
- printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Check *** Common ancestor label %d, pdg %d, name %s, status %d; \n",
- ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
- }
+ //else {
+ // printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Check *** Common ancestor label %d, pdg %d, name %s, status %d; \n",
+ // ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
+ // }
}////Partons, colliding protons, strings, intermediate corrections
}//ancLabel > -1
else { //ancLabel <= -1
}
}
+//____________________________________________________________________________________________________________________________________________________
+void AliAnaPi0::CountAndGetAverages(Int_t &nClus,Int_t &nCell, Float_t &eClusTot,Float_t &eCellTot, Float_t &eDenClus,Float_t &eDenCell) {
+// Count the number of clusters and cells, deposited energy, and do some averages in case multiplicity bins dependent on such numbers
+// are requested
+ if(fCalorimeter=="EMCAL"){
+ nClus = GetEMCALClusters() ->GetEntriesFast();
+ nCell = GetEMCALCells()->GetNumberOfCells();
+ for(Int_t icl=0; icl < nClus; icl++) {
+ Float_t e1 = ((AliVCluster*)GetEMCALClusters()->At(icl))->E();
+ eClusTot += e1;
+ // if(e1 > emax) emax = e1;
+ // ((AliVCluster*)GetEMCALClusters()->At(icl))->GetPosition(pos1);
+ // for(Int_t icl2=icl+1; icl2 < nClus; icl2++) {
+ // Float_t e2 = ((AliVCluster*)GetEMCALClusters()->At(icl2))->E();
+ // ((AliVCluster*)GetEMCALClusters()->At(icl2))->GetPosition(pos2);
+ // rtmp = TMath::Sqrt((pos1[0]-pos2[0])*(pos1[0]-pos2[0]) + (pos1[2]-pos2[2])*(pos1[2]-pos2[2]));
+ // rtmpw = TMath::Sqrt((pos1[0]*e1-pos2[0]*e2)*(pos1[0]*e1-pos2[0]*e2) + (pos1[2]*e1-pos2[2]*e2)*(pos1[2]*e1-pos2[2]*e2))/(e1+e2);
+ // rxz += rtmp;
+ // rxzw += rtmpw;
+ // ncomb++;
+ // fhClusterPairDist ->Fill(rtmp);
+ // fhClusterPairDistWeight->Fill(rtmpw);
+ // //printf("Distance: %f; weighted %f\n ",rtmp,rtmp/(e1+((AliVCluster*)GetEMCALClusters()->At(icl2))->E()));
+ //
+ // }// second cluster loop
+ }// first cluster
+
+ for(Int_t jce=0; jce < nCell; jce++) eCellTot += GetEMCALCells()->GetAmplitude(jce);
+ }
+ else {
+ nClus = GetPHOSClusters()->GetEntriesFast();
+ nCell = GetPHOSCells() ->GetNumberOfCells();
+ for(Int_t icl=0; icl < nClus; icl++) {
+ Float_t e1 = ((AliVCluster*)GetPHOSClusters()->At(icl))->E();
+ eClusTot += e1;
+ // ((AliVCluster*)GetPHOSClusters()->At(icl))->GetPosition(pos1);
+ // for(Int_t icl2=icl+1; icl2 < nClus; icl2++) {
+ // Float_t e2 = ((AliVCluster*)GetPHOSClusters()->At(icl2))->E();
+ // ((AliVCluster*)GetPHOSClusters()->At(icl2))->GetPosition(pos2);
+ // rtmp = TMath::Sqrt((pos1[0]-pos2[0])*(pos1[0]-pos2[0]) + (pos1[2]-pos2[2])*(pos1[2]-pos2[2]));
+ // rtmpw = TMath::Sqrt((pos1[0]*e1-pos2[0]*e2)*(pos1[0]*e1-pos2[0]*e2) + (pos1[2]*e1-pos2[2]*e2)*(pos1[2]*e1-pos2[2]*e2))/(e1+e2);
+ // rxz += rtmp;
+ // rxzw += rtmpw;
+ // ncomb++;
+ // fhClusterPairDist ->Fill(rtmp);
+ // fhClusterPairDistWeight->Fill(rtmpw);
+ // }// second cluster loop
+ }// first cluster
+ for(Int_t jce=0; jce < nCell; jce++) eCellTot += GetPHOSCells()->GetAmplitude(jce);
+ }
+ if(GetDebug() > 1)
+ printf("AliAnaPi0::MakeAnalysisFillHistograms() - # Clusters %d, sum cluster E per SM %f,# Cells %d, sum cell E per SM %f\n", nClus,eClusTot,nCell,eCellTot);
+
+ //Fill histograms with "energy density", ncell and nclust will be > 0 since there are at least 2 "photons"
+ eDenClus = eClusTot/nClus;
+ eDenCell = eCellTot/nCell;
+ fhEDensityCluster ->Fill(eDenClus);
+ fhEDensityCell ->Fill(eDenCell);
+ fhEDensityCellvsCluster->Fill(eDenClus, eDenCell);
+ //Fill the average number of cells or clusters per SM
+ eClusTot /=fNModules;
+ eCellTot /=fNModules;
+ fhAverTotECluster ->Fill(eClusTot);
+ fhAverTotECell ->Fill(eCellTot);
+ fhAverTotECellvsCluster->Fill(eClusTot, eCellTot);
+ //printf("Average Cluster: E %f, density %f; Average Cell E %f, density %f\n ",eClusTot,eDenClus,eCellTot,eDenCell);
+
+ // //Average weighted pair distance
+ // rxz /= ncomb;
+ // rxzw /= ncomb;
+ //
+ // fhAverClusterPairDist ->Fill(rxz );
+ // fhAverClusterPairDistWeight ->Fill(rxzw);
+ // fhAverClusterPairDistvsAverE ->Fill(rxz ,eDenClus);
+ // fhAverClusterPairDistWeightvsAverE->Fill(rxzw,eDenClus);
+ // fhAverClusterPairDistvsN ->Fill(rxz ,nClus);
+ // fhAverClusterPairDistWeightvsN ->Fill(rxzw,nClus);
+ //
+ // //emax
+ // fhMaxEvsClustEDen->Fill(emax,eDenClus);
+ // fhMaxEvsClustMult->Fill(emax,nPhot);
+
+ //printf("Average Distance: %f; weighted %f\n ",rxz,rxzw);
+
+}
+
//____________________________________________________________________________________________________________________________________________________
void AliAnaPi0::MakeAnalysisFillHistograms()
{
// Float_t pos2[3];
// Float_t emax = 0;
+ if(GetNCentrBin() > 1 && (fUseAverCellEBins||fUseAverClusterEBins||fUseAverClusterEDenBins))
+ CountAndGetAverages(nClus,nCell,eClusTot,eCellTot,eDenClus,eDenCell);
+
+
if(GetDebug() > 1)
printf("AliAnaPi0::MakeAnalysisFillHistograms() - Photon entries %d\n", nPhot);
//If less than photon 2 entries in the list, skip this event
- if(nPhot < 2 ) return ;
-
- // Count the number of clusters and cells, in case multiplicity bins dependent on such numbers
- // are requested
- if(fCalorimeter=="EMCAL"){
- nClus = GetEMCALClusters() ->GetEntriesFast();
- nCell = GetEMCALCells()->GetNumberOfCells();
- for(Int_t icl=0; icl < nClus; icl++) {
- Float_t e1 = ((AliVCluster*)GetEMCALClusters()->At(icl))->E();
- eClusTot += e1;
-// if(e1 > emax) emax = e1;
-// ((AliVCluster*)GetEMCALClusters()->At(icl))->GetPosition(pos1);
-// for(Int_t icl2=icl+1; icl2 < nClus; icl2++) {
-// Float_t e2 = ((AliVCluster*)GetEMCALClusters()->At(icl2))->E();
-// ((AliVCluster*)GetEMCALClusters()->At(icl2))->GetPosition(pos2);
-// rtmp = TMath::Sqrt((pos1[0]-pos2[0])*(pos1[0]-pos2[0]) + (pos1[2]-pos2[2])*(pos1[2]-pos2[2]));
-// rtmpw = TMath::Sqrt((pos1[0]*e1-pos2[0]*e2)*(pos1[0]*e1-pos2[0]*e2) + (pos1[2]*e1-pos2[2]*e2)*(pos1[2]*e1-pos2[2]*e2))/(e1+e2);
-// rxz += rtmp;
-// rxzw += rtmpw;
-// ncomb++;
-// fhClusterPairDist ->Fill(rtmp);
-// fhClusterPairDistWeight->Fill(rtmpw);
-// //printf("Distance: %f; weighted %f\n ",rtmp,rtmp/(e1+((AliVCluster*)GetEMCALClusters()->At(icl2))->E()));
-//
-// }// second cluster loop
- }// first cluster
+ if(nPhot < 2 ) {
- for(Int_t jce=0; jce < nCell; jce++) eCellTot += GetEMCALCells()->GetAmplitude(jce);
- }
- else {
- nClus = GetPHOSClusters()->GetEntriesFast();
- nCell = GetPHOSCells() ->GetNumberOfCells();
- for(Int_t icl=0; icl < nClus; icl++) {
- Float_t e1 = ((AliVCluster*)GetPHOSClusters()->At(icl))->E();
- eClusTot += e1;
-// ((AliVCluster*)GetPHOSClusters()->At(icl))->GetPosition(pos1);
-// for(Int_t icl2=icl+1; icl2 < nClus; icl2++) {
-// Float_t e2 = ((AliVCluster*)GetPHOSClusters()->At(icl2))->E();
-// ((AliVCluster*)GetPHOSClusters()->At(icl2))->GetPosition(pos2);
-// rtmp = TMath::Sqrt((pos1[0]-pos2[0])*(pos1[0]-pos2[0]) + (pos1[2]-pos2[2])*(pos1[2]-pos2[2]));
-// rtmpw = TMath::Sqrt((pos1[0]*e1-pos2[0]*e2)*(pos1[0]*e1-pos2[0]*e2) + (pos1[2]*e1-pos2[2]*e2)*(pos1[2]*e1-pos2[2]*e2))/(e1+e2);
-// rxz += rtmp;
-// rxzw += rtmpw;
-// ncomb++;
-// fhClusterPairDist ->Fill(rtmp);
-// fhClusterPairDistWeight->Fill(rtmpw);
-// }// second cluster loop
- }// first cluster
- for(Int_t jce=0; jce < nCell; jce++) eCellTot += GetPHOSCells()->GetAmplitude(jce);
+ if(GetDebug() > 2)
+ printf("AliAnaPi0::MakeAnalysisFillHistograms() - nPhotons %d, cent bin %d continue to next event\n",nPhot, GetEventCentrality());
+
+ if(GetNCentrBin() > 1) fhCentralityNoPair->Fill(GetEventCentrality() * GetNCentrBin() / GetReader()->GetCentralityOpt());
+
+ return ;
}
- if(GetDebug() > 1)
- printf("AliAnaPi0::MakeAnalysisFillHistograms() - # Clusters %d, sum cluster E per SM %f,# Cells %d, sum cell E per SM %f\n", nClus,eClusTot,nCell,eCellTot);
-
- //Fill histograms with "energy density", ncell and nclust will be > 0 since there are at least 2 "photons"
- eDenClus = eClusTot/nClus;
- eDenCell = eCellTot/nCell;
- fhEDensityCluster ->Fill(eDenClus);
- fhEDensityCell ->Fill(eDenCell);
- fhEDensityCellvsCluster->Fill(eDenClus, eDenCell);
- //Fill the average number of cells or clusters per SM
- eClusTot /=fNModules;
- eCellTot /=fNModules;
- fhAverTotECluster ->Fill(eClusTot);
- fhAverTotECell ->Fill(eCellTot);
- fhAverTotECellvsCluster->Fill(eClusTot, eCellTot);
- //printf("Average Cluster: E %f, density %f; Average Cell E %f, density %f\n ",eClusTot,eDenClus,eCellTot,eDenCell);
-
-// //Average weighted pair distance
-// rxz /= ncomb;
-// rxzw /= ncomb;
-//
-// fhAverClusterPairDist ->Fill(rxz );
-// fhAverClusterPairDistWeight ->Fill(rxzw);
-// fhAverClusterPairDistvsAverE ->Fill(rxz ,eDenClus);
-// fhAverClusterPairDistWeightvsAverE->Fill(rxzw,eDenClus);
-// fhAverClusterPairDistvsN ->Fill(rxz ,nClus);
-// fhAverClusterPairDistWeightvsN ->Fill(rxzw,nClus);
-//
-// //emax
-// fhMaxEvsClustEDen->Fill(emax,eDenClus);
-// fhMaxEvsClustMult->Fill(emax,nPhot);
-
- //printf("Average Distance: %f; weighted %f\n ",rxz,rxzw);
-
//Init variables
Int_t module1 = -1;
Int_t curRPBin = 0 ;
Int_t curZvertBin = 0 ;
+ //Get shower shape information of clusters
+ TObjArray *clusters = 0;
+ if (fCalorimeter="EMCAL") clusters = GetEMCALClusters();
+ else if(fCalorimeter="PHOS" ) clusters = GetPHOSClusters() ;
+
//---------------------------------
//First loop on photons/clusters
//---------------------------------
if(fUseTrackMultBins){ // Track multiplicity bins
//printf("track mult %d\n",GetTrackMultiplicity());
curCentrBin = (GetTrackMultiplicity()-1)/5;
- if(curCentrBin > fNCentrBin-1) curCentrBin=fNCentrBin-1;
+ if(curCentrBin > GetNCentrBin()-1) curCentrBin=GetNCentrBin()-1;
//printf("track mult bin %d\n",curCentrBin);
}
else if(fUsePhotonMultBins){ // Photon multiplicity bins
//printf("photon mult %d cluster mult %d\n",nPhot, nClus);
- curRPBin = nPhot-2;
- if(curRPBin > GetNRPBin() -1) curRPBin=GetNRPBin()-1;
+ curCentrBin = nPhot-2;
+ if(curCentrBin > GetNCentrBin() -1) curCentrBin=GetNCentrBin()-1;
//printf("photon mult bin %d\n",curRPBin);
}
else if(fUseAverClusterEBins){ // Cluster average energy bins
//Bins for pp, if needed can be done in a more general way
- curCentrBin = (Int_t) eClusTot/10 * fNCentrBin;
- if(curCentrBin > fNCentrBin-1) curCentrBin=fNCentrBin-1;
+ curCentrBin = (Int_t) eClusTot/10 * GetNCentrBin();
+ if(curCentrBin > GetNCentrBin()-1) curCentrBin=GetNCentrBin()-1;
//printf("cluster E average %f, bin %d \n",eClusTot,curCentrBin);
}
else if(fUseAverCellEBins){ // Cell average energy bins
//Bins for pp, if needed can be done in a more general way
- curCentrBin = (Int_t) eCellTot/10*fNCentrBin;
- if(curCentrBin > fNCentrBin-1) curCentrBin=fNCentrBin-1;
+ curCentrBin = (Int_t) eCellTot/10*GetNCentrBin();
+ if(curCentrBin > GetNCentrBin()-1) curCentrBin=GetNCentrBin()-1;
//printf("cell E average %f, bin %d \n",eCellTot,curCentrBin);
}
else if(fUseAverClusterEDenBins){ // Energy density bins
//Bins for pp, if needed can be done in a more general way
- curCentrBin = (Int_t) eDenClus/10*fNCentrBin;
- if(curCentrBin > fNCentrBin-1) curCentrBin=fNCentrBin-1;
+ curCentrBin = (Int_t) eDenClus/10*GetNCentrBin();
+ if(curCentrBin > GetNCentrBin()-1) curCentrBin=GetNCentrBin()-1;
//printf("cluster Eden average %f, bin %d \n",eDenClus,curCentrBin);
}
// else if(fUseAverClusterPairRBins){ // Cluster average distance bins
// //Bins for pp, if needed can be done in a more general way
-// curCentrBin = rxz/650*fNCentrBin;
-// if(curCentrBin > fNCentrBin-1) curCentrBin=fNCentrBin-1;
+// curCentrBin = rxz/650*GetNCentrBin();
+// if(curCentrBin > GetNCentrBin()-1) curCentrBin=GetNCentrBin()-1;
// //printf("cluster pair R average %f, bin %d \n",rxz,curCentrBin);
// }
// else if(fUseAverClusterPairRWeightBins){ // Cluster average distance bins
// //Bins for pp, if needed can be done in a more general way
-// curCentrBin = rxzw/350*fNCentrBin;
-// if(curCentrBin > fNCentrBin-1) curCentrBin=fNCentrBin-1;
+// curCentrBin = rxzw/350*GetNCentrBin();
+// if(curCentrBin > GetNCentrBin()-1) curCentrBin=GetNCentrBin()-1;
// //printf("cluster pair rW average %f, bin %d \n",rxzw,curCentrBin);
// }
// else if(fUseEMaxBins){ // Cluster average distance bins
// //Bins for pp, if needed can be done in a more general way
-// curCentrBin = emax/20*fNCentrBin;
-// if(curCentrBin > fNCentrBin-1) curCentrBin=fNCentrBin-1;
+// curCentrBin = emax/20*GetNCentrBin();
+// if(curCentrBin > GetNCentrBin()-1) curCentrBin=GetNCentrBin()-1;
// //printf("cluster pair rW average %f, bin %d \n",rxzw,curCentrBin);
// }
else { //Event centrality
- curCentrBin = GetEventCentrality();
- //printf("curCentrBin %d\n",curCentrBin);
+ // Centrality task returns at maximum 10, 20 or 100, depending on option chosen and
+ // number of bins, the bin has to be corrected
+ curCentrBin = GetEventCentrality() * GetNCentrBin() / GetReader()->GetCentralityOpt();
+ if(GetDebug() > 0 )printf("AliAnaPi0::MakeAnalysisFillHistograms() - curCentrBin %d, centrality %d, n bins %d, max bin from centrality %d\n",
+ curCentrBin, GetEventCentrality(), GetNCentrBin(), GetReader()->GetCentralityOpt());
}
- if (curCentrBin < 0 || curCentrBin >= fNCentrBin){
+ if (curCentrBin < 0 || curCentrBin >= GetNCentrBin()){
if(GetDebug() > 0)
- printf("AliAnaPi0::MakeAnalysisFillHistograms() - Centrality bin <%d> not expected, n bins <%d> , return\n",curCentrBin,fNCentrBin);
+ printf("AliAnaPi0::MakeAnalysisFillHistograms() - Centrality bin <%d> not expected, n bins <%d> , return\n",curCentrBin,GetNCentrBin());
return;
}
//Reaction plane bin
curRPBin = 0 ;
-
+ if(GetNRPBin()>1 && GetEventPlane()){
+ Float_t epAngle = GetEventPlane()->GetEventplane(GetEventPlaneMethod());
+ fhEventPlaneAngle->Fill(epAngle);
+ curRPBin = TMath::Nint(epAngle*(GetNRPBin()-1)/TMath::Pi());
+ if(curRPBin >= GetNRPBin()) printf("RP Bin %d out of range %d\n",curRPBin,GetNRPBin());
+ //printf("RP: %d, %f, angle %f, n bin %d\n", curRPBin,epAngle*(GetNRPBin()-1)/TMath::Pi(),epAngle,GetNRPBin());
+ }
+
//Get vertex z bin
curZvertBin = (Int_t)(0.5*GetNZvertBin()*(vert[2]+GetZvertexCut())/GetZvertexCut()) ;
//Fill event bin info
fhEvents ->Fill(curCentrBin+0.5,curZvertBin+0.5,curRPBin+0.5) ;
- fhCentrality->Fill(curCentrBin);
+ if(GetNCentrBin() > 1) {
+ fhCentrality->Fill(curCentrBin);
+ if(GetNRPBin() > 1 && GetEventPlane()) fhEventPlaneResolution->Fill(curCentrBin,TMath::Cos(2.*GetEventPlane()->GetQsubRes()));
+ }
currentEvtIndex = evtIndex1 ;
if(GetDebug() > 1)
printf("AliAnaPi0::MakeAnalysisFillHistograms() - Centrality %d, Vertex Bin %d, RP bin %d \n",curCentrBin,curRPBin,curZvertBin);
//Get (Super)Module number of this cluster
module1 = GetModuleNumber(p1);
+
+ //Get original cluster time,
+ AliVCluster *cluster1 = 0;
+ Bool_t bFound1 = kFALSE;
+ Int_t caloLabel1 = p1->GetCaloLabel(0);
+ Bool_t iclus1 =-1;
+ if(clusters){
+ for(Int_t iclus = 0; iclus < clusters->GetEntriesFast(); iclus++){
+ AliVCluster *cluster= dynamic_cast<AliVCluster*> (clusters->At(iclus));
+ if(cluster){
+ if (cluster->GetID()==caloLabel1) {
+ bFound1 = kTRUE ;
+ cluster1 = cluster;
+ iclus1 = iclus;
+ }
+ }
+ if(bFound1) break;
+ }
+ }// calorimeter clusters loop
+
//---------------------------------
//Second loop on photons/clusters
//---------------------------------
for(Int_t i2=i1+1; i2<nPhot; i2++){
AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i2)) ;
-
+
//In case of mixing frame, check we are not in the same event as the first cluster
Int_t evtIndex2 = GetEventIndex(p2, vert) ;
if ( evtIndex2 == -1 )
continue ;
if (GetMixedEvent() && (evtIndex1 == evtIndex2))
continue ;
+
+ //------------------------------------------
+ //Check if time of clusters is similar
+ AliVCluster *cluster2 = 0;
+ Bool_t bFound2 = kFALSE;
+ Int_t caloLabel2 = p2->GetCaloLabel(0);
+ if(clusters){
+ for(Int_t iclus = iclus1+1; iclus < clusters->GetEntriesFast(); iclus++){
+ AliVCluster *cluster= dynamic_cast<AliVCluster*> (clusters->At(iclus));
+ if(cluster){
+ if(cluster->GetID()==caloLabel2) {
+ bFound2 = kTRUE ;
+ cluster2 = cluster;
+ }
+ }
+ if(bFound2) break;
+ }// calorimeter clusters loop
+ }
+ Float_t tof1 = -1;
+ if(cluster1 && bFound1){
+ tof1 = cluster1->GetTOF()*1e9;
+ }
+// else printf("cluster1 not available: calo label %d / %d, cluster ID %d\n",
+// p1->GetCaloLabel(0),(GetReader()->GetInputEvent())->GetNumberOfCaloClusters()-1,cluster1->GetID());
+
+ Float_t tof2 = -1;
+ if(cluster2 && bFound2){
+ tof2 = cluster2->GetTOF()*1e9;
+ }
+// else printf("cluster2 not available: calo label %d / %d, cluster ID %d\n",
+// p2->GetCaloLabel(0),(GetReader()->GetInputEvent())->GetNumberOfCaloClusters()-1,cluster2->GetID());
+
+ if(clusters){
+ Double_t t12diff = tof1-tof2;
+ if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
+ }
+ //------------------------------------------
+
//printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 2 Evt %d Vertex : %f,%f,%f\n",evtIndex2, GetVertex(evtIndex2)[0] ,GetVertex(evtIndex2)[1],GetVertex(evtIndex2)[2]);
-
+
//Get the momentum of this cluster
TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
//Get module number
//-------------------------------------------------------------------------------------------------
//Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
//-------------------------------------------------------------------------------------------------
- if(a < fAsymCuts[0]){
+ if(a < fAsymCuts[0] && fFillSMCombinations){
if(module1==module2 && module1 >=0 && module1<fNModules)
fhReMod[module1]->Fill(pt,m) ;
if(ok){
//Check if one of the clusters comes from a conversion
- if (p1->IsTagged() && p2->IsTagged()) fhReConv2->Fill(pt,m);
- else if(p1->IsTagged() || p2->IsTagged()) fhReConv ->Fill(pt,m);
+ if(fCheckConversion){
+ if (p1->IsTagged() && p2->IsTagged()) fhReConv2->Fill(pt,m);
+ else if(p1->IsTagged() || p2->IsTagged()) fhReConv ->Fill(pt,m);
+ }
//Fill histograms for different bad channel distance, centrality, assymmetry cut and pid bit
for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell]){
fhRePtNCellAsymCuts[index]->Fill(pt,m) ;
//printf("ipt %d, icell%d, iasym %d, name %s\n",ipt, icell, iasym, fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->GetName());
- if(module1==module2){
+ if(fFillSMCombinations && module1==module2){
if (module1==0) fhRePtNCellAsymCutsSM0[index]->Fill(pt,m) ;
else if(module1==1) fhRePtNCellAsymCutsSM1[index]->Fill(pt,m) ;
else if(module1==2) fhRePtNCellAsymCutsSM2[index]->Fill(pt,m) ;
else if(module1==3) fhRePtNCellAsymCutsSM3[index]->Fill(pt,m) ;
- else printf("AliAnaPi0::FillHistograms() - WRONG SM NUMBER\n");
+ //else printf("AliAnaPi0::FillHistograms() - WRONG SM NUMBER\n");
}
}
}// pid bit cut loop
}// icell loop
}// pt cut loop
- for(Int_t iasym = 0; iasym < fNAsymCuts; iasym++){
- if(a < fAsymCuts[iasym])fhRePtMult[iasym]->Fill(pt,GetTrackMultiplicity(),m) ;
+ if(GetHistoTrackMultiplicityBins()){
+ for(Int_t iasym = 0; iasym < fNAsymCuts; iasym++){
+ if(a < fAsymCuts[iasym])fhRePtMult[iasym]->Fill(pt,GetTrackMultiplicity(),m) ;
+ }
}
}// multiple cuts analysis
}// ok if same sm
//-------------------------------------------------------------------------------------------------
//Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
//-------------------------------------------------------------------------------------------------
- if(a < fAsymCuts[0]){
+ if(a < fAsymCuts[0] && fFillSMCombinations){
if(module1==module2 && module1 >=0 && module1<fNModules)
fhMiMod[module1]->Fill(pt,m) ;
if(ok){
//Check if one of the clusters comes from a conversion
- if (p1->IsTagged() && p2->IsTagged()) fhMiConv2->Fill(pt,m);
- else if(p1->IsTagged() || p2->IsTagged()) fhMiConv ->Fill(pt,m);
-
+ if(fCheckConversion){
+ if (p1->IsTagged() && p2->IsTagged()) fhMiConv2->Fill(pt,m);
+ else if(p1->IsTagged() || p2->IsTagged()) fhMiConv ->Fill(pt,m);
+ }
//Fill histograms for different bad channel distance, centrality, assymmetry cut and pid bit
for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton))){
if(currentEvent->GetEntriesFast()>0){
evMixList->AddFirst(currentEvent) ;
currentEvent=0 ; //Now list of particles belongs to buffer and it will be deleted with buffer
- if(evMixList->GetSize()>=fNmaxMixEv)
+ if(evMixList->GetSize() >= GetNMaxEvMix())
{
TClonesArray * tmp = (TClonesArray*) (evMixList->Last()) ;
evMixList->RemoveLast() ;
// the first one and then add the next.
Int_t index = outputList->IndexOf(outputList->FindObject(GetAddedHistogramsStringToName()+"hRe_cen0_pid0_dist1"));
- if(!fhRe1) fhRe1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
- if(!fhRe2) fhRe2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
- if(!fhRe3) fhRe3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
- if(!fhMi1) fhMi1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
- if(!fhMi2) fhMi2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
- if(!fhMi3) fhMi3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
- if(!fhReInvPt1) fhReInvPt1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
- if(!fhReInvPt2) fhReInvPt2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
- if(!fhReInvPt3) fhReInvPt3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
- if(!fhMiInvPt1) fhMiInvPt1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
- if(!fhMiInvPt2) fhMiInvPt2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
- if(!fhMiInvPt3) fhMiInvPt3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
- if(!fhReMod) fhReMod = new TH2D*[fNModules] ;
- if(!fhReDiffPHOSMod) fhReDiffPHOSMod = new TH2D*[fNModules] ;
- if(!fhReSameSectorEMCALMod)fhReSameSectorEMCALMod = new TH2D*[fNModules/2] ;
- if(!fhReSameSideEMCALMod) fhReSameSideEMCALMod = new TH2D*[fNModules-2] ;
- if(!fhMiMod) fhMiMod = new TH2D*[fNModules] ;
- if(!fhMiDiffPHOSMod) fhMiDiffPHOSMod = new TH2D*[fNModules] ;
- if(!fhMiSameSectorEMCALMod)fhMiSameSectorEMCALMod = new TH2D*[fNModules/2] ;
- if(!fhMiSameSideEMCALMod) fhMiSameSideEMCALMod = new TH2D*[fNModules-2] ;
-
- fhReConv = (TH2D*) outputList->At(index++);
- fhMiConv = (TH2D*) outputList->At(index++);
- fhReConv2 = (TH2D*) outputList->At(index++);
- fhMiConv2 = (TH2D*) outputList->At(index++);
-
- for(Int_t ic=0; ic<fNCentrBin; ic++){
+ if(!fhRe1) fhRe1 = new TH2D*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
+ if(!fhRe2) fhRe2 = new TH2D*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
+ if(!fhRe3) fhRe3 = new TH2D*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
+ if(!fhMi1) fhMi1 = new TH2D*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
+ if(!fhMi2) fhMi2 = new TH2D*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
+ if(!fhMi3) fhMi3 = new TH2D*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
+ if(!fhReInvPt1) fhReInvPt1 = new TH2D*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
+ if(!fhReInvPt2) fhReInvPt2 = new TH2D*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
+ if(!fhReInvPt3) fhReInvPt3 = new TH2D*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
+ if(!fhMiInvPt1) fhMiInvPt1 = new TH2D*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
+ if(!fhMiInvPt2) fhMiInvPt2 = new TH2D*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
+ if(!fhMiInvPt3) fhMiInvPt3 = new TH2D*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
+
+ if(fFillSMCombinations){
+ if(!fhReMod) fhReMod = new TH2D*[fNModules] ;
+ if(!fhReDiffPHOSMod) fhReDiffPHOSMod = new TH2D*[fNModules] ;
+ if(!fhReSameSectorEMCALMod)fhReSameSectorEMCALMod = new TH2D*[fNModules/2] ;
+ if(!fhReSameSideEMCALMod) fhReSameSideEMCALMod = new TH2D*[fNModules-2] ;
+ if(!fhMiMod) fhMiMod = new TH2D*[fNModules] ;
+ if(!fhMiDiffPHOSMod) fhMiDiffPHOSMod = new TH2D*[fNModules] ;
+ if(!fhMiSameSectorEMCALMod)fhMiSameSectorEMCALMod = new TH2D*[fNModules/2] ;
+ if(!fhMiSameSideEMCALMod) fhMiSameSideEMCALMod = new TH2D*[fNModules-2] ;
+ }
+
+ if(fCheckConversion){
+ fhReConv = (TH2D*) outputList->At(index++);
+ fhMiConv = (TH2D*) outputList->At(index++);
+ fhReConv2 = (TH2D*) outputList->At(index++);
+ fhMiConv2 = (TH2D*) outputList->At(index++);
+ }
+
+ for(Int_t ic=0; ic<GetNCentrBin(); ic++){
for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
Int_t ihisto = ((ic*fNPIDBits)+ipid)*fNAsymCuts + iasym;
}// multi cut analysis
fhEvents = (TH3D *) outputList->At(index++);
- fhCentrality = (TH1D *) outputList->At(index++);
+ if(GetNCentrBin()>1)fhCentrality = (TH1D *) outputList->At(index++);
+ if(GetNCentrBin()>1)fhCentralityNoPair = (TH1D *) outputList->At(index++);
+ if(GetNRPBin()>1) fhEventPlaneAngle = (TH1D *) outputList->At(index++);
+ if(GetNRPBin()>1 && GetNCentrBin()>1)fhEventPlaneResolution = (TH2D *) outputList->At(index++);
fhRealOpeningAngle = (TH2D*) outputList->At(index++);
fhRealCosOpeningAngle = (TH2D*) outputList->At(index++);
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff