ClassImp(AliAnaPi0)
-//________________________________________________________________________________________________________________________________________________
+//______________________________________________________
AliAnaPi0::AliAnaPi0() : AliAnaCaloTrackCorrBaseClass(),
-fDoOwnMix(kFALSE), fEventsList(0x0),
+fEventsList(0x0),
fCalorimeter(""), fNModules(12),
fUseAngleCut(kFALSE), fUseAngleEDepCut(kFALSE), fAngleCut(0), fAngleMaxCut(7.),
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), fFillBadDistHisto(kFALSE), fFillSSCombinations(kFALSE),
-fFillAngleHisto(kFALSE), fFillAsymmetryHisto(kFALSE),
+fMakeInvPtPlots(kFALSE), fSameSM(kFALSE),
+fFillSMCombinations(kFALSE), fCheckConversion(kFALSE),
+fFillBadDistHisto(kFALSE), fFillSSCombinations(kFALSE),
+fFillAngleHisto(kFALSE), fFillAsymmetryHisto(kFALSE), fFillOriginHisto(0), fFillArmenterosThetaStar(0),
//Histograms
fhAverTotECluster(0), fhAverTotECell(0), fhAverTotECellvsCluster(0),
fhEDensityCluster(0), fhEDensityCell(0), fhEDensityCellvsCluster(0),
fhRePtNCellAsymCuts(0x0), fhMiPtNCellAsymCuts(0x0), fhRePtNCellAsymCutsSM(),
fhRePIDBits(0x0), fhRePtMult(0x0), fhReSS(),
fhRePtAsym(0x0), fhRePtAsymPi0(0x0), fhRePtAsymEta(0x0),
-fhEvents(0x0), fhCentrality(0x0), fhCentralityNoPair(0x0),
-fhEventPlaneAngle(0x0), fhEventPlaneResolution(0x0),
+fhEventBin(0), fhEventMixBin(0),
+fhCentrality(0x0), fhCentralityNoPair(0x0),
+fhEventPlaneResolution(0x0),
fhRealOpeningAngle(0x0), fhRealCosOpeningAngle(0x0), fhMixedOpeningAngle(0x0), fhMixedCosOpeningAngle(0x0),
// MC histograms
-fhPrimPi0Pt(0x0), fhPrimPi0AccPt(0x0), fhPrimPi0Y(0x0), fhPrimPi0AccY(0x0),
-fhPrimPi0Phi(0x0), fhPrimPi0AccPhi(0x0), fhPrimPi0OpeningAngle(0x0), fhPrimPi0CosOpeningAngle(0x0),
-fhPrimEtaPt(0x0), fhPrimEtaAccPt(0x0), fhPrimEtaY(0x0), fhPrimEtaAccY(0x0),
-fhPrimEtaPhi(0x0), fhPrimEtaAccPhi(0x0), fhPrimPi0PtOrigin(0x0), fhPrimEtaPtOrigin(0x0),
+fhPrimPi0E(0x0), fhPrimPi0Pt(0x0),
+fhPrimPi0AccE(0x0), fhPrimPi0AccPt(0x0),
+fhPrimPi0Y(0x0), fhPrimPi0AccY(0x0),
+fhPrimPi0Phi(0x0), fhPrimPi0AccPhi(0x0),
+fhPrimPi0OpeningAngle(0x0), fhPrimPi0OpeningAngleAsym(0x0),fhPrimPi0CosOpeningAngle(0x0),
+fhPrimPi0PtCentrality(0), fhPrimPi0PtEventPlane(0),
+fhPrimPi0AccPtCentrality(0), fhPrimPi0AccPtEventPlane(0),
+fhPrimEtaE(0x0), fhPrimEtaPt(0x0),
+fhPrimEtaAccE(0x0), fhPrimEtaAccPt(0x0),
+fhPrimEtaY(0x0), fhPrimEtaAccY(0x0),
+fhPrimEtaPhi(0x0), fhPrimEtaAccPhi(0x0),
+fhPrimEtaOpeningAngle(0x0), fhPrimEtaOpeningAngleAsym(0x0),fhPrimEtaCosOpeningAngle(0x0),
+fhPrimEtaPtCentrality(0), fhPrimEtaPtEventPlane(0),
+fhPrimEtaAccPtCentrality(0), fhPrimEtaAccPtEventPlane(0),
+fhPrimPi0PtOrigin(0x0), fhPrimEtaPtOrigin(0x0),
fhMCOrgMass(), fhMCOrgAsym(), fhMCOrgDeltaEta(), fhMCOrgDeltaPhi(),
fhMCPi0MassPtRec(), fhMCPi0MassPtTrue(), fhMCPi0PtTruePtRec(),
fhMCEtaMassPtRec(), fhMCEtaMassPtTrue(), fhMCEtaPtTruePtRec(),
fhMCPi0PtOrigin(0x0), fhMCEtaPtOrigin(0x0),
-fhReMCFromConversion(0), fhReMCFromNotConversion(0), fhReMCFromMixConversion(0)
+fhReMCFromConversion(0), fhReMCFromNotConversion(0), fhReMCFromMixConversion(0),
+fhCosThStarPrimPi0(0), fhCosThStarPrimEta(0)//,
{
-//Default Ctor
- InitParameters();
+ //Default Ctor
+ InitParameters();
+
+ for(Int_t i = 0; i < 4; i++)
+ {
+ fhArmPrimEta[i] = 0;
+ fhArmPrimPi0[i] = 0;
+ }
}
-//________________________________________________________________________________________________________________________________________________
-AliAnaPi0::~AliAnaPi0() {
+//_____________________
+AliAnaPi0::~AliAnaPi0()
+{
// Remove event containers
- if(fDoOwnMix && fEventsList){
- 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() ;
- delete fEventsList[ic*GetNZvertBin()*GetNRPBin()+iz*GetNRPBin()+irp] ;
+ if(DoOwnMix() && fEventsList)
+ {
+ for(Int_t ic=0; ic<GetNCentrBin(); ic++)
+ {
+ for(Int_t iz=0; iz<GetNZvertBin(); iz++)
+ {
+ for(Int_t irp=0; irp<GetNRPBin(); irp++)
+ {
+ Int_t bin = GetEventMixBin(ic,iz,irp);
+ fEventsList[bin]->Delete() ;
+ delete fEventsList[bin] ;
}
}
}
}
-//________________________________________________________________________________________________________________________________________________
+//______________________________
void AliAnaPi0::InitParameters()
{
//Init parameters when first called the analysis
}
-//________________________________________________________________________________________________________________________________________________
+//_______________________________________
TObjString * AliAnaPi0::GetAnalysisCuts()
{
//Save parameters used for analysis
parList+=onePar ;
snprintf(onePar,buffersize,"Depth of event buffer: %d \n",GetNMaxEvMix()) ;
parList+=onePar ;
- 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 ;
snprintf(onePar,buffersize," Asymmetry cuts: n = %d, asymmetry < ",fNAsymCuts) ;
return new TObjString(parList) ;
}
-//________________________________________________________________________________________________________________________________________________
+//_________________________________________
TList * AliAnaPi0::GetCreateOutputObjects()
{
// Create histograms to be saved in output file and
//create event containers
fEventsList = new TList*[GetNCentrBin()*GetNZvertBin()*GetNRPBin()] ;
-
- 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() ;
- fEventsList[ic*GetNZvertBin()*GetNRPBin()+iz*GetNRPBin()+irp]->SetOwner(kFALSE);
+
+ for(Int_t ic=0; ic<GetNCentrBin(); ic++)
+ {
+ for(Int_t iz=0; iz<GetNZvertBin(); iz++)
+ {
+ for(Int_t irp=0; irp<GetNRPBin(); irp++)
+ {
+ Int_t bin = GetEventMixBin(ic,iz,irp);
+ fEventsList[bin] = new TList() ;
+ fEventsList[bin]->SetOwner(kFALSE);
}
}
}
fhReMod = new TH2F*[fNModules] ;
fhMiMod = new TH2F*[fNModules] ;
- if(fCalorimeter == "PHOS"){
+ if(fCalorimeter == "PHOS")
+ {
fhReDiffPHOSMod = new TH2F*[fNModules] ;
fhMiDiffPHOSMod = new TH2F*[fNModules] ;
}
- else{
+ else
+ {
fhReSameSectorEMCALMod = new TH2F*[fNModules/2] ;
fhReSameSideEMCALMod = new TH2F*[fNModules-2] ;
fhMiSameSectorEMCALMod = new TH2F*[fNModules/2] ;
fhRe1 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
fhMi1 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
- if(fFillBadDistHisto){
+ if(fFillBadDistHisto)
+ {
fhRe2 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
fhRe3 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
fhMi2 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
fhMi3 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
}
- if(fMakeInvPtPlots) {
+ if(fMakeInvPtPlots)
+ {
fhReInvPt1 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
fhMiInvPt1 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
if(fFillBadDistHisto){
Int_t ntrmbins = GetHistogramRanges()->GetHistoTrackMultiplicityBins();
Int_t ntrmmax = GetHistogramRanges()->GetHistoTrackMultiplicityMax();
Int_t ntrmmin = GetHistogramRanges()->GetHistoTrackMultiplicityMin();
-
- 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) ;
-
- }//counting and average histograms
-
- if(fCheckConversion){
+ if(fCheckConversion)
+ {
fhReConv = new TH2F("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})");
fhReConv2->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
outputContainer->Add(fhReConv2) ;
- if(fDoOwnMix){
+ if(DoOwnMix())
+ {
fhMiConv = new TH2F("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})");
}
}
- for(Int_t ic=0; ic<GetNCentrBin(); ic++){
- for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
- for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
+ 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;
//printf("cen %d, pid %d, asy %d, Index %d\n",ic,ipid,iasym,index);
//Distance to bad module 1
//printf("name: %s\n ",fhRe1[index]->GetName());
outputContainer->Add(fhRe1[index]) ;
- if(fFillBadDistHisto){
+ if(fFillBadDistHisto)
+ {
//Distance to bad module 2
snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
}
//Inverse pT
- if(fMakeInvPtPlots){
+ if(fMakeInvPtPlots)
+ {
//Distance to bad module 1
snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
outputContainer->Add(fhReInvPt3[index]) ;
}
}
- if(fDoOwnMix){
+
+ if(DoOwnMix())
+ {
//Distance to bad module 1
snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
fhMi3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
outputContainer->Add(fhMi3[index]) ;
}
+
//Inverse pT
- if(fMakeInvPtPlots){
+ if(fMakeInvPtPlots)
+ {
//Distance to bad module 1
snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
}
}
- if(fFillAsymmetryHisto){
+ if(fFillAsymmetryHisto)
+ {
fhRePtAsym = new TH2F("hRePtAsym","Asymmetry vs pt, for pairs",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
fhRePtAsym->SetXTitle("p_{T} (GeV/c)");
fhRePtAsym->SetYTitle("Asymmetry");
outputContainer->Add(fhRePtAsymEta);
}
- if(fMultiCutAna){
-
+ if(fMultiCutAna)
+ {
fhRePIDBits = new TH2F*[fNPIDBits];
for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
snprintf(key, buffersize,"hRe_pidbit%d",ipid) ;
fhRePtNCellAsymCuts = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
fhMiPtNCellAsymCuts = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
- if(fFillSMCombinations){
+ if(fFillSMCombinations)
for(Int_t iSM = 0; iSM < fNModules; iSM++) fhRePtNCellAsymCutsSM[iSM] = new TH2F*[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++){
+ for(Int_t ipt=0; ipt<fNPtCuts; ipt++)
+ {
+ for(Int_t icell=0; icell<fNCellNCuts; icell++)
+ {
+ for(Int_t iasym=0; iasym<fNAsymCuts; iasym++)
+ {
snprintf(key, buffersize,"hRe_pt%d_cell%d_asym%d",ipt,icell,iasym) ;
snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f ",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
fhMiPtNCellAsymCuts[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
outputContainer->Add(fhMiPtNCellAsymCuts[index]) ;
- if(fFillSMCombinations){
- for(Int_t iSM = 0; iSM < fNModules; iSM++){
+ if(fFillSMCombinations)
+ {
+ for(Int_t iSM = 0; iSM < fNModules; iSM++)
+ {
snprintf(key, buffersize,"hRe_pt%d_cell%d_asym%d_SM%d",ipt,icell,iasym,iSM) ;
snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f, SM %d ",
fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym],iSM) ;
}
}
- if(ntrmbins!=0){
+ if(ntrmbins!=0)
+ {
fhRePtMult = new TH3F*[fNAsymCuts] ;
- for(Int_t iasym = 0; iasym<fNAsymCuts; iasym++){
+ for(Int_t iasym = 0; iasym<fNAsymCuts; iasym++)
+ {
fhRePtMult[iasym] = new TH3F(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)");
outputContainer->Add(fhReSS[2]) ;
}
- fhEvents=new TH3F("hEvents","Number of events",GetNCentrBin(),0.,1.*GetNCentrBin(),
- GetNZvertBin(),0.,1.*GetNZvertBin(),GetNRPBin(),0.,1.*GetNRPBin()) ;
+ if(DoOwnMix())
+ {
+ fhEventBin=new TH1I("hEventBin","Number of real pairs per bin(cen,vz,rp)",
+ GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1,0,
+ GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1) ;
+ fhEventBin->SetXTitle("bin");
+ outputContainer->Add(fhEventBin) ;
+
+
+ fhEventMixBin=new TH1I("hEventMixBin","Number of mixed pairs per bin(cen,vz,rp)",
+ GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1,0,
+ GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1) ;
+ fhEventMixBin->SetXTitle("bin");
+ outputContainer->Add(fhEventMixBin) ;
+ }
- fhEvents->SetXTitle("Centrality bin");
- fhEvents->SetYTitle("Z vertex bin bin");
- fhEvents->SetZTitle("RP bin");
- outputContainer->Add(fhEvents) ;
-
- if(GetNCentrBin()>1){
+ if(GetNCentrBin()>1)
+ {
fhCentrality=new TH1F("hCentralityBin","Number of events in centrality bin",GetNCentrBin(),0.,1.*GetNCentrBin()) ;
fhCentrality->SetXTitle("Centrality bin");
outputContainer->Add(fhCentrality) ;
outputContainer->Add(fhCentralityNoPair) ;
}
- if(GetNRPBin() > 1 ){
-
- fhEventPlaneAngle=new TH1F("hEventPlaneAngleBin","Number of events in centrality bin",100,0.,TMath::TwoPi()) ;
- fhEventPlaneAngle->SetXTitle("EP angle (rad)");
- outputContainer->Add(fhEventPlaneAngle) ;
-
- if(GetNCentrBin()>1){
- fhEventPlaneResolution=new TH2F("hEventPlaneResolution","Event plane resolution",GetNCentrBin(),0,GetNCentrBin(),100,0.,TMath::TwoPi()) ;
- fhEventPlaneResolution->SetYTitle("Resolution");
- fhEventPlaneResolution->SetXTitle("Centrality Bin");
- outputContainer->Add(fhEventPlaneResolution) ;
- }
+ if(GetNRPBin() > 1 && GetNCentrBin()>1 )
+ {
+ fhEventPlaneResolution=new TH2F("hEventPlaneResolution","Event plane resolution",GetNCentrBin(),0,GetNCentrBin(),100,0.,TMath::TwoPi()) ;
+ fhEventPlaneResolution->SetYTitle("Resolution");
+ fhEventPlaneResolution->SetXTitle("Centrality Bin");
+ outputContainer->Add(fhEventPlaneResolution) ;
}
- if(fFillAngleHisto){
+ if(fFillAngleHisto)
+ {
fhRealOpeningAngle = new TH2F
("hRealOpeningAngle","Angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,300,0,TMath::Pi());
fhRealOpeningAngle->SetYTitle("#theta(rad)");
fhRealCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
outputContainer->Add(fhRealCosOpeningAngle) ;
- if(fDoOwnMix){
-
+ if(DoOwnMix())
+ {
fhMixedOpeningAngle = new TH2F
("hMixedOpeningAngle","Angle between all #gamma pair vs E_{#pi^{0}}, Mixed pairs",nptbins,ptmin,ptmax,300,0,TMath::Pi());
fhMixedOpeningAngle->SetYTitle("#theta(rad)");
}
//Histograms filled only if MC data is requested
- if(IsDataMC()){
-
+ if(IsDataMC())
+ {
fhReMCFromConversion = new TH2F("hReMCFromConversion","Invariant mass of 2 clusters originated in conversions",
nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
fhReMCFromConversion->SetXTitle("p_{T} (GeV/c)");
outputContainer->Add(fhReMCFromMixConversion) ;
//Pi0
+
+ fhPrimPi0E = new TH1F("hPrimPi0E","Primary pi0 E, Y<1",nptbins,ptmin,ptmax) ;
+ fhPrimPi0AccE = new TH1F("hPrimPi0AccE","Primary pi0 E with both photons in acceptance",nptbins,ptmin,ptmax) ;
+ fhPrimPi0E ->SetXTitle("E (GeV)");
+ fhPrimPi0AccE->SetXTitle("E (GeV)");
+ outputContainer->Add(fhPrimPi0E) ;
+ outputContainer->Add(fhPrimPi0AccE) ;
+
fhPrimPi0Pt = new TH1F("hPrimPi0Pt","Primary pi0 pt, Y<1",nptbins,ptmin,ptmax) ;
fhPrimPi0AccPt = new TH1F("hPrimPi0AccPt","Primary pi0 pt with both photons in acceptance",nptbins,ptmin,ptmax) ;
fhPrimPi0Pt ->SetXTitle("p_{T} (GeV/c)");
fhPrimPi0AccPhi->SetXTitle("p_{T} (GeV/c)");
outputContainer->Add(fhPrimPi0AccPhi) ;
+ fhPrimPi0PtCentrality = new TH2F("hPrimPi0PtCentrality","Primary pi0 pt vs reco centrality, Y<1",nptbins,ptmin,ptmax, 100, 0, 100) ;
+ fhPrimPi0AccPtCentrality = new TH2F("hPrimPi0AccPtCentrality","Primary pi0 with both photons in acceptance pt vs reco centrality",nptbins,ptmin,ptmax, 100, 0, 100) ;
+ fhPrimPi0PtCentrality ->SetXTitle("p_{T} (GeV/c)");
+ fhPrimPi0AccPtCentrality->SetXTitle("p_{T} (GeV/c)");
+ fhPrimPi0PtCentrality ->SetYTitle("Centrality");
+ fhPrimPi0AccPtCentrality->SetYTitle("Centrality");
+ outputContainer->Add(fhPrimPi0PtCentrality) ;
+ outputContainer->Add(fhPrimPi0AccPtCentrality) ;
+
+ fhPrimPi0PtEventPlane = new TH2F("hPrimPi0PtEventPlane","Primary pi0 pt vs reco event plane angle, Y<1",nptbins,ptmin,ptmax, 100, 0, TMath::Pi()) ;
+ fhPrimPi0AccPtEventPlane = new TH2F("hPrimPi0AccPtEventPlane","Primary pi0 with both photons in acceptance pt vs reco event plane angle",nptbins,ptmin,ptmax, 100, 0, TMath::Pi()) ;
+ fhPrimPi0PtEventPlane ->SetXTitle("p_{T} (GeV/c)");
+ fhPrimPi0AccPtEventPlane->SetXTitle("p_{T} (GeV/c)");
+ fhPrimPi0PtEventPlane ->SetYTitle("Event Plane Angle (rad)");
+ fhPrimPi0AccPtEventPlane->SetYTitle("Event Plane Angle (rad)");
+ outputContainer->Add(fhPrimPi0PtEventPlane) ;
+ outputContainer->Add(fhPrimPi0AccPtEventPlane) ;
+
//Eta
+
+ fhPrimEtaE = new TH1F("hPrimEtaE","Primary eta E",nptbins,ptmin,ptmax) ;
+ fhPrimEtaAccE = new TH1F("hPrimEtaAccE","Primary eta E with both photons in acceptance",nptbins,ptmin,ptmax) ;
+ fhPrimEtaE ->SetXTitle("E (GeV)");
+ fhPrimEtaAccE->SetXTitle("E (GeV)");
+ outputContainer->Add(fhPrimEtaE) ;
+ outputContainer->Add(fhPrimEtaAccE) ;
+
fhPrimEtaPt = new TH1F("hPrimEtaPt","Primary eta pt",nptbins,ptmin,ptmax) ;
fhPrimEtaAccPt = new TH1F("hPrimEtaAccPt","Primary eta pt with both photons in acceptance",nptbins,ptmin,ptmax) ;
fhPrimEtaPt ->SetXTitle("p_{T} (GeV/c)");
fhPrimEtaAccPhi->SetXTitle("p_{T} (GeV/c)");
outputContainer->Add(fhPrimEtaAccPhi) ;
-
- //Prim origin
- //Pi0
- fhPrimPi0PtOrigin = new TH2F("hPrimPi0PtOrigin","Primary pi0 pt vs origin",nptbins,ptmin,ptmax,11,0,11) ;
- fhPrimPi0PtOrigin->SetXTitle("p_{T} (GeV/c)");
- fhPrimPi0PtOrigin->SetYTitle("Origin");
- fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
- fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
- fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances ");
- fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
- fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(5 ,"#rho");
- fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(6 ,"#omega");
- fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(7 ,"K");
- fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(8 ,"Other");
- fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(9 ,"#eta");
- fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(10 ,"#eta prime");
- outputContainer->Add(fhPrimPi0PtOrigin) ;
-
- fhMCPi0PtOrigin = new TH2F("hMCPi0PtOrigin","Reconstructed pair from generated pi0 pt vs origin",nptbins,ptmin,ptmax,11,0,11) ;
- fhMCPi0PtOrigin->SetXTitle("p_{T} (GeV/c)");
- fhMCPi0PtOrigin->SetYTitle("Origin");
- fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
- fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
- fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
- fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
- fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(5 ,"#rho");
- fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(6 ,"#omega");
- fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(7 ,"K");
- fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(8 ,"Other");
- fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(9 ,"#eta");
- fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(10 ,"#eta prime");
- outputContainer->Add(fhMCPi0PtOrigin) ;
-
- //Eta
- fhPrimEtaPtOrigin = new TH2F("hPrimEtaPtOrigin","Primary pi0 pt vs origin",nptbins,ptmin,ptmax,7,0,7) ;
- fhPrimEtaPtOrigin->SetXTitle("p_{T} (GeV/c)");
- fhPrimEtaPtOrigin->SetYTitle("Origin");
- fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
- fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
- fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
- fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
- fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(5 ,"Other");
- fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(6 ,"#eta prime ");
-
- outputContainer->Add(fhPrimEtaPtOrigin) ;
-
- fhMCEtaPtOrigin = new TH2F("hMCEtaPtOrigin","Reconstructed pair from generated pi0 pt vs origin",nptbins,ptmin,ptmax,7,0,7) ;
- fhMCEtaPtOrigin->SetXTitle("p_{T} (GeV/c)");
- fhMCEtaPtOrigin->SetYTitle("Origin");
- fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
- fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
- fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
- fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
- fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(5 ,"Other");
- fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(6 ,"#eta prime");
-
- outputContainer->Add(fhMCEtaPtOrigin) ;
-
- if(fFillAngleHisto){
+ fhPrimEtaPtCentrality = new TH2F("hPrimEtaPtCentrality","Primary eta pt vs reco centrality, Y<1",nptbins,ptmin,ptmax, 100, 0, 100) ;
+ fhPrimEtaAccPtCentrality = new TH2F("hPrimEtaAccPtCentrality","Primary eta with both photons in acceptance pt vs reco centrality",nptbins,ptmin,ptmax, 100, 0, 100) ;
+ fhPrimEtaPtCentrality ->SetXTitle("p_{T} (GeV/c)");
+ fhPrimEtaAccPtCentrality->SetXTitle("p_{T} (GeV/c)");
+ fhPrimEtaPtCentrality ->SetYTitle("Centrality");
+ fhPrimEtaAccPtCentrality->SetYTitle("Centrality");
+ outputContainer->Add(fhPrimEtaPtCentrality) ;
+ outputContainer->Add(fhPrimEtaAccPtCentrality) ;
+
+ fhPrimEtaPtEventPlane = new TH2F("hPrimEtaPtEventPlane","Primary eta pt vs reco event plane angle, Y<1",nptbins,ptmin,ptmax, 100, 0, TMath::Pi()) ;
+ fhPrimEtaAccPtEventPlane = new TH2F("hPrimEtaAccPtEventPlane","Primary eta with both photons in acceptance pt vs reco event plane angle",nptbins,ptmin,ptmax, 100, 0, TMath::Pi()) ;
+ fhPrimEtaPtEventPlane ->SetXTitle("p_{T} (GeV/c)");
+ fhPrimEtaAccPtEventPlane->SetXTitle("p_{T} (GeV/c)");
+ fhPrimEtaPtEventPlane ->SetYTitle("Event Plane Angle (rad)");
+ fhPrimEtaAccPtEventPlane->SetYTitle("Event Plane Angle (rad)");
+ outputContainer->Add(fhPrimEtaPtEventPlane) ;
+ outputContainer->Add(fhPrimEtaAccPtEventPlane) ;
+
+ if(fFillAngleHisto)
+ {
fhPrimPi0OpeningAngle = new TH2F
("hPrimPi0OpeningAngle","Angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,0,0.5);
fhPrimPi0OpeningAngle->SetYTitle("#theta(rad)");
fhPrimPi0OpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
outputContainer->Add(fhPrimPi0OpeningAngle) ;
+ fhPrimPi0OpeningAngleAsym = new TH2F
+ ("hPrimPi0OpeningAngleAsym","Angle between all primary #gamma pair vs Asymmetry, p_{T}>5 GeV/c",100,0,1,100,0,0.5);
+ fhPrimPi0OpeningAngleAsym->SetXTitle("|A|=| (E_{1}-E_{2}) / (E_{1}+E_{2}) |");
+ fhPrimPi0OpeningAngleAsym->SetYTitle("#theta(rad)");
+ outputContainer->Add(fhPrimPi0OpeningAngleAsym) ;
+
fhPrimPi0CosOpeningAngle = new TH2F
("hPrimPi0CosOpeningAngle","Cosinus of angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,-1,1);
fhPrimPi0CosOpeningAngle->SetYTitle("cos (#theta) ");
fhPrimPi0CosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
outputContainer->Add(fhPrimPi0CosOpeningAngle) ;
- }
-
- for(Int_t i = 0; i<13; i++){
- fhMCOrgMass[i] = new TH2F(Form("hMCOrgMass_%d",i),Form("mass vs pt, origin %d",i),nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
- fhMCOrgMass[i]->SetXTitle("p_{T} (GeV/c)");
- fhMCOrgMass[i]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
- outputContainer->Add(fhMCOrgMass[i]) ;
- fhMCOrgAsym[i]= new TH2F(Form("hMCOrgAsym_%d",i),Form("asymmetry vs pt, origin %d",i),nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
- fhMCOrgAsym[i]->SetXTitle("p_{T} (GeV/c)");
- fhMCOrgAsym[i]->SetYTitle("A");
- outputContainer->Add(fhMCOrgAsym[i]) ;
+ fhPrimEtaOpeningAngle = new TH2F
+ ("hPrimEtaOpeningAngle","Angle between all primary #gamma pair vs E_{#eta}",nptbins,ptmin,ptmax,100,0,0.5);
+ fhPrimEtaOpeningAngle->SetYTitle("#theta(rad)");
+ fhPrimEtaOpeningAngle->SetXTitle("E_{#eta} (GeV)");
+ outputContainer->Add(fhPrimEtaOpeningAngle) ;
- fhMCOrgDeltaEta[i] = new TH2F(Form("hMCOrgDeltaEta_%d",i),Form("#Delta #eta of pair vs pt, origin %d",i),nptbins,ptmin,ptmax,netabins,-1.4,1.4) ;
- fhMCOrgDeltaEta[i]->SetXTitle("p_{T} (GeV/c)");
- fhMCOrgDeltaEta[i]->SetYTitle("#Delta #eta");
- outputContainer->Add(fhMCOrgDeltaEta[i]) ;
+ fhPrimEtaOpeningAngleAsym = new TH2F
+ ("hPrimEtaOpeningAngleAsym","Angle between all primary #gamma pair vs Asymmetry, p_{T}>5 GeV/c",100,0,1,100,0,0.5);
+ fhPrimEtaOpeningAngleAsym->SetXTitle("|A|=| (E_{1}-E_{2}) / (E_{1}+E_{2}) |");
+ fhPrimEtaOpeningAngleAsym->SetYTitle("#theta(rad)");
+ outputContainer->Add(fhPrimEtaOpeningAngleAsym) ;
+
- fhMCOrgDeltaPhi[i]= new TH2F(Form("hMCOrgDeltaPhi_%d",i),Form("#Delta #phi of pair vs p_{T}, origin %d",i),nptbins,ptmin,ptmax,nphibins,-0.7,0.7) ;
- fhMCOrgDeltaPhi[i]->SetXTitle("p_{T} (GeV/c)");
- fhMCOrgDeltaPhi[i]->SetYTitle("#Delta #phi (rad)");
- outputContainer->Add(fhMCOrgDeltaPhi[i]) ;
+ fhPrimEtaCosOpeningAngle = new TH2F
+ ("hPrimEtaCosOpeningAngle","Cosinus of angle between all primary #gamma pair vs E_{#eta}",nptbins,ptmin,ptmax,100,-1,1);
+ fhPrimEtaCosOpeningAngle->SetYTitle("cos (#theta) ");
+ fhPrimEtaCosOpeningAngle->SetXTitle("E_{ #eta} (GeV)");
+ outputContainer->Add(fhPrimEtaCosOpeningAngle) ;
+
}
- if(fMultiCutAnaSim){
- fhMCPi0MassPtTrue = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
- fhMCPi0MassPtRec = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
- fhMCPi0PtTruePtRec = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
- fhMCEtaMassPtRec = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
- fhMCEtaMassPtTrue = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
- fhMCEtaPtTruePtRec = new TH2F*[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++){
- Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
-
- fhMCPi0MassPtRec[index] = new TH2F(Form("hMCPi0MassPtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
- Form("Reconstructed Mass vs reconstructed p_T of true #pi^{0} cluster pairs for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
- nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
- fhMCPi0MassPtRec[index]->SetXTitle("p_{T, reconstructed} (GeV/c)");
- fhMCPi0MassPtRec[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
- outputContainer->Add(fhMCPi0MassPtRec[index]) ;
-
- fhMCPi0MassPtTrue[index] = new TH2F(Form("hMCPi0MassPtTrue_pt%d_cell%d_asym%d",ipt,icell,iasym),
- Form("Reconstructed Mass vs generated p_T of true #pi^{0} cluster pairs for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
- nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
- fhMCPi0MassPtTrue[index]->SetXTitle("p_{T, generated} (GeV/c)");
- fhMCPi0MassPtTrue[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
- outputContainer->Add(fhMCPi0MassPtTrue[index]) ;
-
- fhMCPi0PtTruePtRec[index] = new TH2F(Form("hMCPi0PtTruePtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
- Form("Generated vs reconstructed p_T of true #pi^{0} cluster pairs, 0.01 < rec. mass < 0.17 MeV/c^{2} for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
- nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
- fhMCPi0PtTruePtRec[index]->SetXTitle("p_{T, generated} (GeV/c)");
- fhMCPi0PtTruePtRec[index]->SetYTitle("p_{T, reconstructed} (GeV/c)");
- outputContainer->Add(fhMCPi0PtTruePtRec[index]) ;
-
- fhMCEtaMassPtRec[index] = new TH2F(Form("hMCEtaMassPtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
- Form("Reconstructed Mass vs reconstructed p_T of true #eta cluster pairs for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
- nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
- fhMCEtaMassPtRec[index]->SetXTitle("p_{T, generated} (GeV/c)");
- fhMCEtaMassPtRec[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
- outputContainer->Add(fhMCEtaMassPtRec[index]) ;
-
- fhMCEtaMassPtTrue[index] = new TH2F(Form("hMCEtaMassPtTrue_pt%d_cell%d_asym%d",ipt,icell,iasym),
- Form("Reconstructed Mass vs generated p_T of true #eta cluster pairs for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
- nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
- fhMCEtaMassPtTrue[index]->SetXTitle("p_{T, generated} (GeV/c)");
- fhMCEtaMassPtTrue[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
- outputContainer->Add(fhMCEtaMassPtTrue[index]) ;
-
- fhMCEtaPtTruePtRec[index] = new TH2F(Form("hMCEtaPtTruePtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
- Form("Generated vs reconstructed p_T of true #eta cluster pairs, 0.01 < rec. mass < 0.17 MeV/c^{2} for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
- nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
- fhMCEtaPtTruePtRec[index]->SetXTitle("p_{T, generated} (GeV/c)");
- fhMCEtaPtTruePtRec[index]->SetYTitle("p_{T, reconstructed} (GeV/c)");
- outputContainer->Add(fhMCEtaPtTruePtRec[index]) ;
- }
- }
- }
- }//multi cut ana
- else {
- fhMCPi0MassPtTrue = new TH2F*[1];
- fhMCPi0PtTruePtRec = new TH2F*[1];
- fhMCEtaMassPtTrue = new TH2F*[1];
- fhMCEtaPtTruePtRec = new TH2F*[1];
+ if(fFillOriginHisto)
+ {
+ //Prim origin
+ //Pi0
+ fhPrimPi0PtOrigin = new TH2F("hPrimPi0PtOrigin","Primary pi0 pt vs origin",nptbins,ptmin,ptmax,11,0,11) ;
+ fhPrimPi0PtOrigin->SetXTitle("p_{T} (GeV/c)");
+ fhPrimPi0PtOrigin->SetYTitle("Origin");
+ fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
+ fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
+ fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances ");
+ fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
+ fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(5 ,"#rho");
+ fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(6 ,"#omega");
+ fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(7 ,"K");
+ fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(8 ,"Other");
+ fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(9 ,"#eta");
+ fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(10 ,"#eta prime");
+ outputContainer->Add(fhPrimPi0PtOrigin) ;
+
+ fhMCPi0PtOrigin = new TH2F("hMCPi0PtOrigin","Reconstructed pair from generated pi0 pt vs origin",nptbins,ptmin,ptmax,11,0,11) ;
+ fhMCPi0PtOrigin->SetXTitle("p_{T} (GeV/c)");
+ fhMCPi0PtOrigin->SetYTitle("Origin");
+ fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
+ fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
+ fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
+ fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
+ fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(5 ,"#rho");
+ fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(6 ,"#omega");
+ fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(7 ,"K");
+ fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(8 ,"Other");
+ fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(9 ,"#eta");
+ fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(10 ,"#eta prime");
+ outputContainer->Add(fhMCPi0PtOrigin) ;
+
+ //Eta
+ fhPrimEtaPtOrigin = new TH2F("hPrimEtaPtOrigin","Primary pi0 pt vs origin",nptbins,ptmin,ptmax,7,0,7) ;
+ fhPrimEtaPtOrigin->SetXTitle("p_{T} (GeV/c)");
+ fhPrimEtaPtOrigin->SetYTitle("Origin");
+ fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
+ fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
+ fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
+ fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
+ fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(5 ,"Other");
+ fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(6 ,"#eta prime ");
+
+ outputContainer->Add(fhPrimEtaPtOrigin) ;
- fhMCPi0MassPtTrue[0] = new TH2F("hMCPi0MassPtTrue","Reconstructed Mass vs generated p_T of true #pi^{0} cluster pairs",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
- fhMCPi0MassPtTrue[0]->SetXTitle("p_{T, generated} (GeV/c)");
- fhMCPi0MassPtTrue[0]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
- outputContainer->Add(fhMCPi0MassPtTrue[0]) ;
+ fhMCEtaPtOrigin = new TH2F("hMCEtaPtOrigin","Reconstructed pair from generated pi0 pt vs origin",nptbins,ptmin,ptmax,7,0,7) ;
+ fhMCEtaPtOrigin->SetXTitle("p_{T} (GeV/c)");
+ fhMCEtaPtOrigin->SetYTitle("Origin");
+ fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
+ fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
+ fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
+ fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
+ fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(5 ,"Other");
+ fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(6 ,"#eta prime");
- fhMCPi0PtTruePtRec[0]= new TH2F("hMCPi0PtTruePtRec","Generated vs reconstructed p_T of true #pi^{0} cluster pairs, 0.01 < rec. mass < 0.17 MeV/c^{2}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
- fhMCPi0PtTruePtRec[0]->SetXTitle("p_{T, generated} (GeV/c)");
- fhMCPi0PtTruePtRec[0]->SetYTitle("p_{T, reconstructed} (GeV/c)");
- outputContainer->Add(fhMCPi0PtTruePtRec[0]) ;
+ outputContainer->Add(fhMCEtaPtOrigin) ;
- fhMCEtaMassPtTrue[0] = new TH2F("hMCEtaMassPtTrue","Reconstructed Mass vs generated p_T of true #eta cluster pairs",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
- fhMCEtaMassPtTrue[0]->SetXTitle("p_{T, generated} (GeV/c)");
- fhMCEtaMassPtTrue[0]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
- outputContainer->Add(fhMCEtaMassPtTrue[0]) ;
+ for(Int_t i = 0; i<13; i++)
+ {
+ fhMCOrgMass[i] = new TH2F(Form("hMCOrgMass_%d",i),Form("mass vs pt, origin %d",i),nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+ fhMCOrgMass[i]->SetXTitle("p_{T} (GeV/c)");
+ fhMCOrgMass[i]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+ outputContainer->Add(fhMCOrgMass[i]) ;
+
+ fhMCOrgAsym[i]= new TH2F(Form("hMCOrgAsym_%d",i),Form("asymmetry vs pt, origin %d",i),nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
+ fhMCOrgAsym[i]->SetXTitle("p_{T} (GeV/c)");
+ fhMCOrgAsym[i]->SetYTitle("A");
+ outputContainer->Add(fhMCOrgAsym[i]) ;
+
+ fhMCOrgDeltaEta[i] = new TH2F(Form("hMCOrgDeltaEta_%d",i),Form("#Delta #eta of pair vs pt, origin %d",i),nptbins,ptmin,ptmax,netabins,-1.4,1.4) ;
+ fhMCOrgDeltaEta[i]->SetXTitle("p_{T} (GeV/c)");
+ fhMCOrgDeltaEta[i]->SetYTitle("#Delta #eta");
+ outputContainer->Add(fhMCOrgDeltaEta[i]) ;
+
+ fhMCOrgDeltaPhi[i]= new TH2F(Form("hMCOrgDeltaPhi_%d",i),Form("#Delta #phi of pair vs p_{T}, origin %d",i),nptbins,ptmin,ptmax,nphibins,-0.7,0.7) ;
+ fhMCOrgDeltaPhi[i]->SetXTitle("p_{T} (GeV/c)");
+ fhMCOrgDeltaPhi[i]->SetYTitle("#Delta #phi (rad)");
+ outputContainer->Add(fhMCOrgDeltaPhi[i]) ;
+
+ }
- fhMCEtaPtTruePtRec[0]= new TH2F("hMCEtaPtTruePtRec","Generated vs reconstructed p_T of true #eta cluster pairs, 0.01 < rec. mass < 0.17 MeV/c^{2}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
- fhMCEtaPtTruePtRec[0]->SetXTitle("p_{T, generated} (GeV/c)");
- fhMCEtaPtTruePtRec[0]->SetYTitle("p_{T, reconstructed} (GeV/c)");
- outputContainer->Add(fhMCEtaPtTruePtRec[0]) ;
+ if(fMultiCutAnaSim)
+ {
+ fhMCPi0MassPtTrue = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
+ fhMCPi0MassPtRec = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
+ fhMCPi0PtTruePtRec = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
+ fhMCEtaMassPtRec = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
+ fhMCEtaMassPtTrue = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
+ fhMCEtaPtTruePtRec = new TH2F*[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++){
+ Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
+
+ fhMCPi0MassPtRec[index] = new TH2F(Form("hMCPi0MassPtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
+ Form("Reconstructed Mass vs reconstructed p_T of true #pi^{0} cluster pairs for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
+ nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+ fhMCPi0MassPtRec[index]->SetXTitle("p_{T, reconstructed} (GeV/c)");
+ fhMCPi0MassPtRec[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+ outputContainer->Add(fhMCPi0MassPtRec[index]) ;
+
+ fhMCPi0MassPtTrue[index] = new TH2F(Form("hMCPi0MassPtTrue_pt%d_cell%d_asym%d",ipt,icell,iasym),
+ Form("Reconstructed Mass vs generated p_T of true #pi^{0} cluster pairs for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
+ nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+ fhMCPi0MassPtTrue[index]->SetXTitle("p_{T, generated} (GeV/c)");
+ fhMCPi0MassPtTrue[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+ outputContainer->Add(fhMCPi0MassPtTrue[index]) ;
+
+ fhMCPi0PtTruePtRec[index] = new TH2F(Form("hMCPi0PtTruePtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
+ Form("Generated vs reconstructed p_T of true #pi^{0} cluster pairs, 0.01 < rec. mass < 0.17 MeV/c^{2} for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
+ nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
+ fhMCPi0PtTruePtRec[index]->SetXTitle("p_{T, generated} (GeV/c)");
+ fhMCPi0PtTruePtRec[index]->SetYTitle("p_{T, reconstructed} (GeV/c)");
+ outputContainer->Add(fhMCPi0PtTruePtRec[index]) ;
+
+ fhMCEtaMassPtRec[index] = new TH2F(Form("hMCEtaMassPtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
+ Form("Reconstructed Mass vs reconstructed p_T of true #eta cluster pairs for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
+ nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+ fhMCEtaMassPtRec[index]->SetXTitle("p_{T, generated} (GeV/c)");
+ fhMCEtaMassPtRec[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+ outputContainer->Add(fhMCEtaMassPtRec[index]) ;
+
+ fhMCEtaMassPtTrue[index] = new TH2F(Form("hMCEtaMassPtTrue_pt%d_cell%d_asym%d",ipt,icell,iasym),
+ Form("Reconstructed Mass vs generated p_T of true #eta cluster pairs for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
+ nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+ fhMCEtaMassPtTrue[index]->SetXTitle("p_{T, generated} (GeV/c)");
+ fhMCEtaMassPtTrue[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+ outputContainer->Add(fhMCEtaMassPtTrue[index]) ;
+
+ fhMCEtaPtTruePtRec[index] = new TH2F(Form("hMCEtaPtTruePtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
+ Form("Generated vs reconstructed p_T of true #eta cluster pairs, 0.01 < rec. mass < 0.17 MeV/c^{2} for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
+ nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
+ fhMCEtaPtTruePtRec[index]->SetXTitle("p_{T, generated} (GeV/c)");
+ fhMCEtaPtTruePtRec[index]->SetYTitle("p_{T, reconstructed} (GeV/c)");
+ outputContainer->Add(fhMCEtaPtTruePtRec[index]) ;
+ }
+ }
+ }
+ }//multi cut ana
+ else
+ {
+ fhMCPi0MassPtTrue = new TH2F*[1];
+ fhMCPi0PtTruePtRec = new TH2F*[1];
+ fhMCEtaMassPtTrue = new TH2F*[1];
+ fhMCEtaPtTruePtRec = new TH2F*[1];
+
+ fhMCPi0MassPtTrue[0] = new TH2F("hMCPi0MassPtTrue","Reconstructed Mass vs generated p_T of true #pi^{0} cluster pairs",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+ fhMCPi0MassPtTrue[0]->SetXTitle("p_{T, generated} (GeV/c)");
+ fhMCPi0MassPtTrue[0]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+ outputContainer->Add(fhMCPi0MassPtTrue[0]) ;
+
+ fhMCPi0PtTruePtRec[0]= new TH2F("hMCPi0PtTruePtRec","Generated vs reconstructed p_T of true #pi^{0} cluster pairs, 0.01 < rec. mass < 0.17 MeV/c^{2}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
+ fhMCPi0PtTruePtRec[0]->SetXTitle("p_{T, generated} (GeV/c)");
+ fhMCPi0PtTruePtRec[0]->SetYTitle("p_{T, reconstructed} (GeV/c)");
+ outputContainer->Add(fhMCPi0PtTruePtRec[0]) ;
+
+ fhMCEtaMassPtTrue[0] = new TH2F("hMCEtaMassPtTrue","Reconstructed Mass vs generated p_T of true #eta cluster pairs",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+ fhMCEtaMassPtTrue[0]->SetXTitle("p_{T, generated} (GeV/c)");
+ fhMCEtaMassPtTrue[0]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+ outputContainer->Add(fhMCEtaMassPtTrue[0]) ;
+
+ fhMCEtaPtTruePtRec[0]= new TH2F("hMCEtaPtTruePtRec","Generated vs reconstructed p_T of true #eta cluster pairs, 0.01 < rec. mass < 0.17 MeV/c^{2}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
+ fhMCEtaPtTruePtRec[0]->SetXTitle("p_{T, generated} (GeV/c)");
+ fhMCEtaPtTruePtRec[0]->SetYTitle("p_{T, reconstructed} (GeV/c)");
+ outputContainer->Add(fhMCEtaPtTruePtRec[0]) ;
+ }
}
}
- if(fFillSMCombinations){
+ 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++){
+ 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]->SetXTitle("p_{T} (GeV/c)");
fhReMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
outputContainer->Add(fhReMod[imod]) ;
- if(fCalorimeter=="PHOS"){
+ 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 TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
outputContainer->Add(fhReDiffPHOSMod[imod]) ;
}
else{//EMCAL
- if(imod<fNModules/2){
+ 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 TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
fhReSameSectorEMCALMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
outputContainer->Add(fhReSameSectorEMCALMod[imod]) ;
}
- if(imod<fNModules-2){
+ 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 TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
}
}//EMCAL
- if(fDoOwnMix){
+ if(DoOwnMix())
+ {
snprintf(key, buffersize,"hMiMod_%d",imod) ;
snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for Module %d",imod) ;
fhMiMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
outputContainer->Add(fhMiDiffPHOSMod[imod]) ;
}//PHOS
else{//EMCAL
- if(imod<fNModules/2){
+ 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 TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
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 TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
}//loop combinations
} // SM combinations
+ if(fFillArmenterosThetaStar && IsDataMC())
+ {
+ TString ebin[] = {"8 < E < 12 GeV","12 < E < 16 GeV", "16 < E < 20 GeV", "E > 20 GeV" };
+ Int_t narmbins = 400;
+ Float_t armmin = 0;
+ Float_t armmax = 0.4;
+
+ for(Int_t i = 0; i < 4; i++)
+ {
+ fhArmPrimPi0[i] = new TH2F(Form("hArmenterosPrimPi0EBin%d",i),
+ Form("Armenteros of primary #pi^{0}, %s",ebin[i].Data()),
+ 200, -1, 1, narmbins,armmin,armmax);
+ fhArmPrimPi0[i]->SetYTitle("p_{T}^{Arm}");
+ fhArmPrimPi0[i]->SetXTitle("#alpha^{Arm}");
+ outputContainer->Add(fhArmPrimPi0[i]) ;
+
+ fhArmPrimEta[i] = new TH2F(Form("hArmenterosPrimEtaEBin%d",i),
+ Form("Armenteros of primary #eta, %s",ebin[i].Data()),
+ 200, -1, 1, narmbins,armmin,armmax);
+ fhArmPrimEta[i]->SetYTitle("p_{T}^{Arm}");
+ fhArmPrimEta[i]->SetXTitle("#alpha^{Arm}");
+ outputContainer->Add(fhArmPrimEta[i]) ;
+
+ }
+
+ // Same as asymmetry ...
+ fhCosThStarPrimPi0 = new TH2F
+ ("hCosThStarPrimPi0","cos(#theta *) for primary #pi^{0}",nptbins,ptmin,ptmax,200,-1,1);
+ fhCosThStarPrimPi0->SetYTitle("cos(#theta *)");
+ fhCosThStarPrimPi0->SetXTitle("E_{ #pi^{0}} (GeV)");
+ outputContainer->Add(fhCosThStarPrimPi0) ;
+
+ fhCosThStarPrimEta = new TH2F
+ ("hCosThStarPrimEta","cos(#theta *) for primary #eta",nptbins,ptmin,ptmax,200,-1,1);
+ fhCosThStarPrimEta->SetYTitle("cos(#theta *)");
+ fhCosThStarPrimEta->SetXTitle("E_{ #eta} (GeV)");
+ outputContainer->Add(fhCosThStarPrimEta) ;
+
+ }
+
// for(Int_t i = 0; i < outputContainer->GetEntries() ; i++){
//
// printf("Histogram %d, name: %s\n ",i, outputContainer->At(i)->GetName());
return outputContainer;
}
-//_________________________________________________________________________________________________________________________________________________
+//___________________________________________________
void AliAnaPi0::Print(const Option_t * /*opt*/) const
{
//Print some relevant parameters set for the analysis
printf("------------------------------------------------------\n") ;
}
-//_____________________________________________________________
-void AliAnaPi0::FillAcceptanceHistograms(){
+//________________________________________
+void AliAnaPi0::FillAcceptanceHistograms()
+{
//Fill acceptance histograms if MC data is available
- if(GetReader()->ReadStack()){
+ Float_t cen = GetEventCentrality();
+ Float_t ep = GetEventPlaneAngle();
+
+ if(GetReader()->ReadStack())
+ {
AliStack * stack = GetMCStack();
- if(stack){
- for(Int_t i=0 ; i<stack->GetNtrack(); i++){
+ if(stack)
+ {
+ for(Int_t i=0 ; i<stack->GetNtrack(); i++)
+ {
+ if(GetReader()->AcceptOnlyHIJINGLabels() && !GetReader()->IsHIJINGLabel(i)) continue ;
+
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){
+ if( pdg == 111 || pdg == 221)
+ {
Double_t pi0Pt = prim->Pt() ;
- 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 pi0E = prim->Energy() ;
+ if(pi0E == TMath::Abs(prim->Pz())) continue ; //Protection against floating point exception
+ Double_t pi0Y = 0.5*TMath::Log((pi0E-prim->Pz())/(pi0E+prim->Pz())) ;
Double_t phi = TMath::RadToDeg()*prim->Phi() ;
- if(pdg == 111){
- if(TMath::Abs(pi0Y) < 1.0){
+ if(pdg == 111)
+ {
+ if(TMath::Abs(pi0Y) < 1.0)
+ {
+ fhPrimPi0E ->Fill(pi0E ) ;
fhPrimPi0Pt ->Fill(pi0Pt) ;
fhPrimPi0Phi->Fill(pi0Pt, phi) ;
+ fhPrimPi0PtCentrality->Fill(pi0Pt,cen) ;
+ fhPrimPi0PtEventPlane->Fill(pi0Pt,ep ) ;
}
fhPrimPi0Y ->Fill(pi0Pt, pi0Y) ;
}
- else if(pdg == 221){
- if(TMath::Abs(pi0Y) < 1.0){
+ else if(pdg == 221)
+ {
+ if(TMath::Abs(pi0Y) < 1.0)
+ {
+ fhPrimEtaE ->Fill(pi0E ) ;
fhPrimEtaPt ->Fill(pi0Pt) ;
fhPrimEtaPhi->Fill(pi0Pt, phi) ;
+ fhPrimEtaPtCentrality->Fill(pi0Pt,cen) ;
+ fhPrimEtaPtEventPlane->Fill(pi0Pt,ep ) ;
}
fhPrimEtaY ->Fill(pi0Pt, pi0Y) ;
}
//Origin of meson
- Int_t momindex = prim->GetFirstMother();
- 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
+ if(fFillOriginHisto)
+ {
+ Int_t momindex = prim->GetFirstMother();
+ 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) 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()){
+ if(iphot1>-1 && iphot1<stack->GetNtrack() && iphot2>-1 && iphot2<stack->GetNtrack())
+ {
TParticle * phot1 = stack->Particle(iphot1) ;
TParticle * phot2 = stack->Particle(iphot2) ;
- if(phot1 && phot2 && phot1->GetPdgCode()==22 && phot2->GetPdgCode()==22){
+ if(phot1 && phot2 && phot1->GetPdgCode()==22 && phot2->GetPdgCode()==22)
+ {
//printf("2 photons: photon 1: pt %2.2f, phi %3.2f, eta %1.2f; photon 2: pt %2.2f, phi %3.2f, eta %1.2f\n",
// phot1->Pt(), phot1->Phi()*180./3.1415, phot1->Eta(), phot2->Pt(), phot2->Phi()*180./3.1415, phot2->Eta());
- TLorentzVector lv1, lv2;
+ TLorentzVector lv1, lv2,lvmeson;
phot1->Momentum(lv1);
phot2->Momentum(lv2);
+ prim ->Momentum(lvmeson);
+
+ if(fFillArmenterosThetaStar) FillArmenterosThetaStar(pdg,lvmeson,lv1,lv2);
+
Bool_t inacceptance = kFALSE;
- if(fCalorimeter == "PHOS"){
- if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet()){
+ if(fCalorimeter == "PHOS")
+ {
+ if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet())
+ {
Int_t mod ;
Double_t x,z ;
if(GetPHOSGeometry()->ImpactOnEmc(phot1,mod,z,x) && GetPHOSGeometry()->ImpactOnEmc(phot2,mod,z,x))
inacceptance = kTRUE;
if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
}
- else{
-
+ else
+ {
if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
inacceptance = kTRUE ;
if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
}
-
}
- else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet()){
- if(GetEMCALGeometry()){
-
+ else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet())
+ {
+ if(GetEMCALGeometry())
+ {
Int_t absID1=0;
Int_t absID2=0;
// inacceptance = kTRUE;
if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
}
- else{
+ else
+ {
if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
inacceptance = kTRUE ;
if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
}
}
- if(inacceptance){
- if(pdg==111){
+ if(inacceptance)
+ {
+ Float_t asym = TMath::Abs((lv1.E()-lv2.E()) / (lv1.E()+lv2.E()));
+ Double_t angle = lv1.Angle(lv2.Vect());
+
+ if(pdg==111)
+ {
+ fhPrimPi0AccE ->Fill(pi0E) ;
fhPrimPi0AccPt ->Fill(pi0Pt) ;
fhPrimPi0AccPhi->Fill(pi0Pt, phi) ;
fhPrimPi0AccY ->Fill(pi0Pt, pi0Y) ;
- if(fFillAngleHisto){
- Double_t angle = lv1.Angle(lv2.Vect());
- fhPrimPi0OpeningAngle ->Fill(pi0Pt,angle);
- fhPrimPi0CosOpeningAngle->Fill(pi0Pt,TMath::Cos(angle));
+ fhPrimPi0AccPtCentrality->Fill(pi0Pt,cen) ;
+ fhPrimPi0AccPtEventPlane->Fill(pi0Pt,ep ) ;
+
+ if(fFillAngleHisto)
+ {
+ fhPrimPi0OpeningAngle ->Fill(pi0Pt,angle);
+ if(pi0Pt > 5)fhPrimPi0OpeningAngleAsym->Fill(asym,angle);
+ fhPrimPi0CosOpeningAngle ->Fill(pi0Pt,TMath::Cos(angle));
}
}
- else if(pdg==221){
+ else if(pdg==221)
+ {
+ fhPrimEtaAccE ->Fill(pi0E ) ;
fhPrimEtaAccPt ->Fill(pi0Pt) ;
fhPrimEtaAccPhi->Fill(pi0Pt, phi) ;
fhPrimEtaAccY ->Fill(pi0Pt, pi0Y) ;
+ fhPrimEtaAccPtCentrality->Fill(pi0Pt,cen) ;
+ fhPrimEtaAccPtEventPlane->Fill(pi0Pt,ep ) ;
+
+ if(fFillAngleHisto)
+ {
+ fhPrimEtaOpeningAngle ->Fill(pi0Pt,angle);
+ if(pi0Pt > 5)fhPrimEtaOpeningAngleAsym->Fill(asym,angle);
+ fhPrimEtaCosOpeningAngle ->Fill(pi0Pt,TMath::Cos(angle));
+ }
}
}//Accepted
}// 2 photons
}//loop on primaries
}//stack exists and data is MC
}//read stack
- else if(GetReader()->ReadAODMCParticles()){
- TClonesArray * mcparticles = GetReader()->GetAODMCParticles(0);
- if(mcparticles){
+ else if(GetReader()->ReadAODMCParticles())
+ {
+ TClonesArray * mcparticles = GetReader()->GetAODMCParticles();
+ if(mcparticles)
+ {
Int_t nprim = mcparticles->GetEntriesFast();
for(Int_t i=0; i < nprim; i++)
{
+ if(GetReader()->AcceptOnlyHIJINGLabels() && !GetReader()->IsHIJINGLabel(i)) continue ;
+
AliAODMCParticle * prim = (AliAODMCParticle *) mcparticles->At(i);
// 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){
+ if( pdg == 111 || pdg == 221)
+ {
Double_t pi0Pt = prim->Pt() ;
+ Double_t pi0E = prim->E() ;
//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
+ if(pi0E == 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) < 1){
- fhPrimPi0Pt->Fill(pi0Pt) ;
+ if(pdg == 111)
+ {
+ if(TMath::Abs(pi0Y) < 1)
+ {
+ fhPrimPi0E ->Fill(pi0E ) ;
+ fhPrimPi0Pt ->Fill(pi0Pt) ;
fhPrimPi0Phi->Fill(pi0Pt, phi) ;
+ fhPrimPi0PtCentrality->Fill(pi0Pt,cen) ;
+ fhPrimPi0PtEventPlane->Fill(pi0Pt,ep ) ;
}
fhPrimPi0Y ->Fill(pi0Pt, pi0Y) ;
}
- else if(pdg == 221){
- if(TMath::Abs(pi0Y) < 1){
- fhPrimEtaPt->Fill(pi0Pt) ;
+ else if(pdg == 221)
+ {
+ if(TMath::Abs(pi0Y) < 1)
+ {
+ fhPrimEtaE ->Fill(pi0E ) ;
+ fhPrimEtaPt ->Fill(pi0Pt) ;
fhPrimEtaPhi->Fill(pi0Pt, phi) ;
+ fhPrimEtaPtCentrality->Fill(pi0Pt,cen) ;
+ fhPrimEtaPtEventPlane->Fill(pi0Pt,ep ) ;
}
fhPrimEtaY ->Fill(pi0Pt, pi0Y) ;
}
//Origin of meson
Int_t momindex = prim->GetMother();
- if(momindex >= 0) {
+ 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 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(fFillOriginHisto)
+ {
+ 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
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){
+ if(iphot1>-1 && iphot1<nprim && iphot2>-1 && iphot2<nprim)
+ {
AliAODMCParticle * phot1 = (AliAODMCParticle *) mcparticles->At(iphot1);
AliAODMCParticle * phot2 = (AliAODMCParticle *) mcparticles->At(iphot2);
- if(phot1 && phot2 && phot1->GetPdgCode()==22 && phot2->GetPdgCode()==22){
- TLorentzVector lv1, lv2;
+ if(phot1 && phot2 && phot1->GetPdgCode()==22 && phot2->GetPdgCode()==22)
+ {
+ TLorentzVector lv1, lv2,lvmeson;
lv1.SetPxPyPzE(phot1->Px(),phot1->Py(),phot1->Pz(),phot1->E());
lv2.SetPxPyPzE(phot2->Px(),phot2->Py(),phot2->Pz(),phot2->E());
+ lvmeson.SetPxPyPzE(prim->Px(),prim->Py(),prim->Pz(),prim->E());
+ if(fFillArmenterosThetaStar) FillArmenterosThetaStar(pdg,lvmeson,lv1,lv2);
+
Bool_t inacceptance = kFALSE;
- if(fCalorimeter == "PHOS"){
- if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet()){
+ if(fCalorimeter == "PHOS")
+ {
+ if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet())
+ {
Int_t mod ;
Double_t x,z ;
Double_t vtx []={phot1->Xv(),phot1->Yv(),phot1->Zv()};
inacceptance = kTRUE;
if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
}
- else{
-
+ else
+ {
if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
inacceptance = kTRUE ;
if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
}
-
}
- else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet()){
- if(GetEMCALGeometry()){
-
+ else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet())
+ {
+ if(GetEMCALGeometry())
+ {
Int_t absID1=0;
Int_t absID2=0;
if( absID1 >= 0 && absID2 >= 0)
inacceptance = kTRUE;
-
if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
}
- else{
+ else
+ {
if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
inacceptance = kTRUE ;
if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
}
}
- if(inacceptance){
- if(pdg==111){
+ if(inacceptance)
+ {
+ Float_t asym = TMath::Abs((lv1.E()-lv2.E()) / (lv1.E()+lv2.E()));
+ Double_t angle = lv1.Angle(lv2.Vect());
+
+ if(pdg==111)
+ {
// printf("ACCEPTED pi0: pt %2.2f, phi %3.2f, eta %1.2f\n",pi0Pt,phi,pi0Y);
+ fhPrimPi0AccE ->Fill(pi0E ) ;
fhPrimPi0AccPt ->Fill(pi0Pt) ;
fhPrimPi0AccPhi->Fill(pi0Pt, phi) ;
fhPrimPi0AccY ->Fill(pi0Pt, pi0Y) ;
- if(fFillAngleHisto){
- Double_t angle = lv1.Angle(lv2.Vect());
- fhPrimPi0OpeningAngle ->Fill(pi0Pt,angle);
- fhPrimPi0CosOpeningAngle->Fill(pi0Pt,TMath::Cos(angle));
+ fhPrimPi0AccPtCentrality->Fill(pi0Pt,cen) ;
+ fhPrimPi0AccPtEventPlane->Fill(pi0Pt,ep ) ;
+
+ if(fFillAngleHisto)
+ {
+ fhPrimPi0OpeningAngle ->Fill(pi0Pt,angle);
+ if(pi0Pt > 5)fhPrimPi0OpeningAngleAsym->Fill(asym,angle);
+ fhPrimPi0CosOpeningAngle ->Fill(pi0Pt,TMath::Cos(angle));
}
}
- else if(pdg==221){
+ else if(pdg==221)
+ {
+ fhPrimEtaAccE ->Fill(pi0E ) ;
fhPrimEtaAccPt ->Fill(pi0Pt) ;
fhPrimEtaAccPhi->Fill(pi0Pt, phi) ;
fhPrimEtaAccY ->Fill(pi0Pt, pi0Y) ;
+ fhPrimEtaAccPtCentrality->Fill(pi0Pt,cen) ;
+ fhPrimEtaAccPtEventPlane->Fill(pi0Pt,ep ) ;
+
+ if(fFillAngleHisto)
+ {
+ fhPrimEtaOpeningAngle ->Fill(pi0Pt,angle);
+ if(pi0Pt > 5)fhPrimEtaOpeningAngleAsym->Fill(asym,angle);
+ fhPrimEtaCosOpeningAngle ->Fill(pi0Pt,TMath::Cos(angle));
+ }
}
}//Accepted
}// 2 photons
} // read AOD MC
}
-//_____________________________________________________________
-void AliAnaPi0::FillMCVersusRecDataHistograms(const Int_t index1, const Int_t index2,
- const Float_t pt1, const Float_t pt2,
- const Int_t ncell1, const Int_t ncell2,
- const Double_t mass, const Double_t pt, const Double_t asym,
- const Double_t deta, const Double_t dphi){
+//__________________________________________________________________________________
+void AliAnaPi0::FillArmenterosThetaStar(Int_t pdg, TLorentzVector meson,
+ TLorentzVector daugh1, TLorentzVector daugh2)
+{
+ // Fill armenteros plots
+
+ // Get pTArm and AlphaArm
+ Float_t momentumSquaredMother = meson.P()*meson.P();
+ Float_t momentumDaughter1AlongMother = 0.;
+ Float_t momentumDaughter2AlongMother = 0.;
+
+ if (momentumSquaredMother > 0.)
+ {
+ momentumDaughter1AlongMother = (daugh1.Px()*meson.Px() + daugh1.Py()*meson.Py()+ daugh1.Pz()*meson.Pz()) / sqrt(momentumSquaredMother);
+ momentumDaughter2AlongMother = (daugh2.Px()*meson.Px() + daugh2.Py()*meson.Py()+ daugh2.Pz()*meson.Pz()) / sqrt(momentumSquaredMother);
+ }
+
+ Float_t momentumSquaredDaughter1 = daugh1.P()*daugh1.P();
+ Float_t ptArmSquared = momentumSquaredDaughter1 - momentumDaughter1AlongMother*momentumDaughter1AlongMother;
+
+ Float_t pTArm = 0.;
+ if (ptArmSquared > 0.)
+ pTArm = sqrt(ptArmSquared);
+
+ Float_t alphaArm = 0.;
+ if(momentumDaughter1AlongMother +momentumDaughter2AlongMother > 0)
+ alphaArm = (momentumDaughter1AlongMother -momentumDaughter2AlongMother) / (momentumDaughter1AlongMother + momentumDaughter2AlongMother);
+
+ TLorentzVector daugh1Boost = daugh1;
+ daugh1Boost.Boost(-meson.BoostVector());
+ Float_t cosThStar=TMath::Cos(daugh1Boost.Vect().Angle(meson.Vect()));
+
+ Float_t en = meson.Energy();
+ Int_t ebin = -1;
+ if(en > 8 && en <= 12) ebin = 0;
+ if(en > 12 && en <= 16) ebin = 1;
+ if(en > 16 && en <= 20) ebin = 2;
+ if(en > 20) ebin = 3;
+ if(ebin < 0 || ebin > 3) return ;
+
+
+ if(pdg==111)
+ {
+ fhCosThStarPrimPi0->Fill(en,cosThStar);
+ fhArmPrimPi0[ebin]->Fill(alphaArm,pTArm);
+ }
+ else
+ {
+ fhCosThStarPrimEta->Fill(en,cosThStar);
+ fhArmPrimEta[ebin]->Fill(alphaArm,pTArm);
+ }
+
+ if(GetDebug() > 2 )
+ {
+ Float_t asym = 0;
+ if(daugh1.E()+daugh2.E() > 0) asym = TMath::Abs(daugh1.E()-daugh2.E())/(daugh1.E()+daugh2.E());
+
+ printf("AliAnaPi0::FillArmenterosThetaStar() - E %f, alphaArm %f, pTArm %f, cos(theta*) %f, asymmetry %1.3f\n",
+ en,alphaArm,pTArm,cosThStar,asym);
+ }
+}
+
+//_______________________________________________________________________________________
+void AliAnaPi0::FillMCVersusRecDataHistograms(Int_t index1, Int_t index2,
+ Float_t pt1, Float_t pt2,
+ Int_t ncell1, Int_t ncell2,
+ Double_t mass, Double_t pt, Double_t asym,
+ Double_t deta, Double_t dphi)
+{
//Do some MC checks on the origin of the pair, is there any common ancestor and if there is one, who?
//Adjusted for Pythia, need to see what to do for other generators.
//Array of histograms ordered as follows: 0-Photon, 1-electron, 2-pi0, 3-eta, 4-a-proton, 5-a-neutron, 6-stable particles,
// 7-other decays, 8-string, 9-final parton, 10-initial parton, intermediate, 11-colliding proton, 12-unrelated
+ if(!fFillOriginHisto) return;
+
Int_t ancPDG = 0;
Int_t ancStatus = 0;
TLorentzVector ancMomentum;
fhMCPi0MassPtTrue[0]->Fill(ancMomentum.Pt(),mass);
if(mass < 0.17 && mass > 0.1) {
fhMCPi0PtTruePtRec[0]->Fill(ancMomentum.Pt(),pt);
-
- if(GetReader()->ReadStack()){
- TParticle* ancestor = GetMCStack()->Particle(ancLabel);
- momindex = ancestor->GetFirstMother();
- if(momindex < 0) return;
- TParticle* mother = GetMCStack()->Particle(momindex);
- mompdg = TMath::Abs(mother->GetPdgCode());
- momstatus = mother->GetStatusCode();
+ if(fFillOriginHisto)
+ {
+ if(GetReader()->ReadStack())
+ {
+ TParticle* ancestor = GetMCStack()->Particle(ancLabel);
+ momindex = ancestor->GetFirstMother();
+ if(momindex < 0) return;
+ TParticle* mother = GetMCStack()->Particle(momindex);
+ mompdg = TMath::Abs(mother->GetPdgCode());
+ momstatus = mother->GetStatusCode();
+ }
+ else
+ {
+ TClonesArray * mcparticles = GetReader()->GetAODMCParticles();
+ AliAODMCParticle* ancestor = (AliAODMCParticle *) mcparticles->At(ancLabel);
+ momindex = ancestor->GetMother();
+ if(momindex < 0) return;
+ AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
+ mompdg = TMath::Abs(mother->GetPdgCode());
+ momstatus = mother->GetStatus();
+ }
+
+ if (momstatus == 21) fhMCPi0PtOrigin->Fill(pt,0.5);//parton
+ else if(mompdg < 22 ) fhMCPi0PtOrigin->Fill(pt,1.5);//quark
+ else if(mompdg > 2100 && mompdg < 2210) fhMCPi0PtOrigin->Fill(pt,2.5);// resonances
+ else if(mompdg == 221) fhMCPi0PtOrigin->Fill(pt,8.5);//eta
+ else if(mompdg == 331) fhMCPi0PtOrigin->Fill(pt,9.5);//eta prime
+ else if(mompdg == 213) fhMCPi0PtOrigin->Fill(pt,4.5);//rho
+ else if(mompdg == 223) fhMCPi0PtOrigin->Fill(pt,5.5);//omega
+ else if(mompdg >= 310 && mompdg <= 323) fhMCPi0PtOrigin->Fill(pt,6.5);//k0S, k+-,k*
+ else if(mompdg == 130) fhMCPi0PtOrigin->Fill(pt,6.5);//k0L
+ else if(momstatus == 11 || momstatus == 12 ) fhMCPi0PtOrigin->Fill(pt,3.5);//resonances
+ else fhMCPi0PtOrigin->Fill(pt,7.5);//other?
}
- else {
- TClonesArray * mcparticles = GetReader()->GetAODMCParticles(0);
- AliAODMCParticle* ancestor = (AliAODMCParticle *) mcparticles->At(ancLabel);
- momindex = ancestor->GetMother();
- if(momindex < 0) return;
- AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
- mompdg = TMath::Abs(mother->GetPdgCode());
- momstatus = mother->GetStatus();
- }
-
- if (momstatus == 21) fhMCPi0PtOrigin->Fill(pt,0.5);//parton
- else if(mompdg < 22 ) fhMCPi0PtOrigin->Fill(pt,1.5);//quark
- else if(mompdg > 2100 && mompdg < 2210) fhMCPi0PtOrigin->Fill(pt,2.5);// resonances
- else if(mompdg == 221) fhMCPi0PtOrigin->Fill(pt,8.5);//eta
- else if(mompdg == 331) fhMCPi0PtOrigin->Fill(pt,9.5);//eta prime
- else if(mompdg == 213) fhMCPi0PtOrigin->Fill(pt,4.5);//rho
- else if(mompdg == 223) fhMCPi0PtOrigin->Fill(pt,5.5);//omega
- else if(mompdg >= 310 && mompdg <= 323) fhMCPi0PtOrigin->Fill(pt,6.5);//k0S, k+-,k*
- else if(mompdg == 130) fhMCPi0PtOrigin->Fill(pt,6.5);//k0L
- else if(momstatus == 11 || momstatus == 12 ) fhMCPi0PtOrigin->Fill(pt,3.5);//resonances
- else fhMCPi0PtOrigin->Fill(pt,7.5);//other?
-
}//pi0 mass region
}
fhMCEtaMassPtTrue[0]->Fill(ancMomentum.Pt(),mass);
if(mass < 0.65 && mass > 0.45) fhMCEtaPtTruePtRec[0]->Fill(ancMomentum.Pt(),pt);
- if(GetReader()->ReadStack()){
- TParticle* ancestor = GetMCStack()->Particle(ancLabel);
- momindex = ancestor->GetFirstMother();
- if(momindex < 0) return;
- TParticle* mother = GetMCStack()->Particle(momindex);
- mompdg = TMath::Abs(mother->GetPdgCode());
- momstatus = mother->GetStatusCode();
+ if(fFillOriginHisto)
+ {
+ if(GetReader()->ReadStack())
+ {
+ TParticle* ancestor = GetMCStack()->Particle(ancLabel);
+ momindex = ancestor->GetFirstMother();
+ if(momindex < 0) return;
+ TParticle* mother = GetMCStack()->Particle(momindex);
+ mompdg = TMath::Abs(mother->GetPdgCode());
+ momstatus = mother->GetStatusCode();
+ }
+ else
+ {
+ TClonesArray * mcparticles = GetReader()->GetAODMCParticles();
+ AliAODMCParticle* ancestor = (AliAODMCParticle *) mcparticles->At(ancLabel);
+ momindex = ancestor->GetMother();
+ if(momindex < 0) return;
+ AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
+ mompdg = TMath::Abs(mother->GetPdgCode());
+ momstatus = mother->GetStatus();
+ }
+
+ if (momstatus == 21 ) fhMCEtaPtOrigin->Fill(pt,0.5);//parton
+ else if(mompdg < 22 ) fhMCEtaPtOrigin->Fill(pt,1.5);//quark
+ else if(mompdg > 2100 && mompdg < 2210) fhMCEtaPtOrigin->Fill(pt,2.5);//qq resonances
+ else if(mompdg == 331) fhMCEtaPtOrigin->Fill(pt,5.5);//eta prime
+ else if(momstatus == 11 || momstatus == 12 ) fhMCEtaPtOrigin->Fill(pt,3.5);//resonances
+ else fhMCEtaPtOrigin->Fill(pt,4.5);//stable, conversions?
+ //printf("Other Meson pdg %d, Mother %s, pdg %d, status %d\n",pdg, TDatabasePDG::Instance()->GetParticle(mompdg)->GetName(),mompdg, momstatus );
}
- else {
- TClonesArray * mcparticles = GetReader()->GetAODMCParticles(0);
- AliAODMCParticle* ancestor = (AliAODMCParticle *) mcparticles->At(ancLabel);
- momindex = ancestor->GetMother();
- if(momindex < 0) return;
- AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
- mompdg = TMath::Abs(mother->GetPdgCode());
- momstatus = mother->GetStatus();
- }
-
- if (momstatus == 21 ) fhMCEtaPtOrigin->Fill(pt,0.5);//parton
- else if(mompdg < 22 ) fhMCEtaPtOrigin->Fill(pt,1.5);//quark
- else if(mompdg > 2100 && mompdg < 2210) fhMCEtaPtOrigin->Fill(pt,2.5);//qq resonances
- else if(mompdg == 331) fhMCEtaPtOrigin->Fill(pt,5.5);//eta prime
- else if(momstatus == 11 || momstatus == 12 ) fhMCEtaPtOrigin->Fill(pt,3.5);//resonances
- else fhMCEtaPtOrigin->Fill(pt,4.5);//stable, conversions?
- //printf("Other Meson pdg %d, Mother %s, pdg %d, status %d\n",pdg, TDatabasePDG::Instance()->GetParticle(mompdg)->GetName(),mompdg, momstatus );
}// eta mass region
}
else if(ancPDG==-2212){//AProton
}
}
-//____________________________________________________________________________________________________________________________________________________
-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;
- }// 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;
- }// 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);
-}
-
-//____________________________________________________________________________________________________________________________________________________
+//__________________________________________
void AliAnaPi0::MakeAnalysisFillHistograms()
{
//Process one event and extract photons from AOD branch
//if (GetReader()->GetEventNumber()%10000 == 0)
// printf("--- Event %d ---\n",GetReader()->GetEventNumber());
+ if(!GetInputAODBranch())
+ {
+ printf("AliAnaPi0::MakeAnalysisFillHistograms() - No input aod photons in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
+ abort();
+ }
+
//Init some variables
Int_t nPhot = GetInputAODBranch()->GetEntriesFast() ;
- Int_t nClus = 0;
- Int_t nCell = 0;
- Float_t eClusTot = 0;
- Float_t eCellTot = 0;
- Float_t eDenClus = 0;
- Float_t eDenCell = 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);
return ;
}
+ Int_t ncentr = GetNCentrBin();
+ if(GetNCentrBin()==0) ncentr = 1;
+
//Init variables
Int_t module1 = -1;
Int_t module2 = -1;
Double_t vert[] = {0.0, 0.0, 0.0} ; //vertex
Int_t evtIndex1 = 0 ;
Int_t currentEvtIndex = -1;
- Int_t curCentrBin = 0 ;
- Int_t curRPBin = 0 ;
- Int_t curZvertBin = 0 ;
+ Int_t curCentrBin = GetEventCentralityBin();
+ //Int_t curVzBin = GetEventVzBin();
+ //Int_t curRPBin = GetEventRPBin();
+ Int_t eventbin = GetEventMixBin();
+ if(eventbin > GetNCentrBin()*GetNZvertBin()*GetNRPBin())
+ {
+ printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mix Bin not exepcted: cen bin %d, z bin %d, rp bin %d, total bin %d, Event Centrality %d, z vertex %2.3f, Reaction Plane %2.3f\n",GetEventCentralityBin(),GetEventVzBin(), GetEventRPBin(),eventbin,GetEventCentrality(),vert[2],GetEventPlaneAngle());
+ return;
+ }
+
//Get shower shape information of clusters
TObjArray *clusters = 0;
if (fCalorimeter=="EMCAL") clusters = GetEMCALClusters();
//---------------------------------
//First loop on photons/clusters
//---------------------------------
- for(Int_t i1=0; i1<nPhot-1; i1++){
+ for(Int_t i1=0; i1<nPhot-1; i1++)
+ {
AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
//printf("AliAnaPi0::MakeAnalysisFillHistograms() : cluster1 id %d\n",p1->GetCaloLabel(0));
if(TMath::Abs(vert[2]) > GetZvertexCut()) continue ; //vertex cut
- //----------------------------------------------------------------------------
- // Get the multiplicity bin. Different cases: centrality (PbPb),
- // average cluster multiplicity, average cell multiplicity, track multiplicity
- // default is centrality bins
- //----------------------------------------------------------------------------
- if (evtIndex1 != currentEvtIndex) {
- if(fUseTrackMultBins){ // Track multiplicity bins
- //printf("track mult %d\n",GetTrackMultiplicity());
- curCentrBin = (GetTrackMultiplicity()-1)/5;
- 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);
- 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 * 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*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*GetNCentrBin();
- if(curCentrBin > GetNCentrBin()-1) curCentrBin=GetNCentrBin()-1;
- //printf("cluster Eden average %f, bin %d \n",eDenClus,curCentrBin);
- }
- else { //Event centrality
- // 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 >= GetNCentrBin()){
- if(GetDebug() > 0)
- 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()) ;
-
+ if (evtIndex1 != currentEvtIndex)
+ {
//Fill event bin info
- fhEvents ->Fill(curCentrBin+0.5,curZvertBin+0.5,curRPBin+0.5) ;
- if(GetNCentrBin() > 1) {
+ if(DoOwnMix()) fhEventBin->Fill(eventbin) ;
+ 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);
}
//printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 1 Evt %d Vertex : %f,%f,%f\n",evtIndex1, GetVertex(evtIndex1)[0] ,GetVertex(evtIndex1)[1],GetVertex(evtIndex1)[2]);
Bool_t bFound1 = kFALSE;
Int_t caloLabel1 = p1->GetCaloLabel(0);
Bool_t iclus1 =-1;
- if(clusters){
+ 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) {
+ if(cluster)
+ {
+ if (cluster->GetID()==caloLabel1)
+ {
bFound1 = kTRUE ;
cluster1 = cluster;
iclus1 = iclus;
//---------------------------------
//Second loop on photons/clusters
//---------------------------------
- for(Int_t i2=i1+1; i2<nPhot; i2++){
+ 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
Float_t tof2 = -1;
Float_t l02 = -1;
- if(cluster2 && bFound2){
+ if(cluster2 && bFound2)
+ {
tof2 = cluster2->GetTOF()*1e9;
l02 = cluster2->GetM02();
// else printf("cluster2 not available: calo label %d / %d, cluster ID %d\n",
// p2->GetCaloLabel(0),(GetReader()->GetInputEvent())->GetNumberOfCaloClusters()-1,cluster2->GetID());
- if(clusters){
+ if(clusters)
+ {
Double_t t12diff = tof1-tof2;
if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
}
//-------------------------------------------------------------------------------------------------
//Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
//-------------------------------------------------------------------------------------------------
- if(a < fAsymCuts[0] && fFillSMCombinations){
+ if(a < fAsymCuts[0] && fFillSMCombinations)
+ {
if(module1==module2 && module1 >=0 && module1<fNModules)
fhReMod[module1]->Fill(pt,m) ;
- if(fCalorimeter=="EMCAL"){
-
+ if(fCalorimeter=="EMCAL")
+ {
// Same sector
Int_t j=0;
- for(Int_t i = 0; i < fNModules/2; i++){
+ for(Int_t i = 0; i < fNModules/2; i++)
+ {
j=2*i;
if((module1==j && module2==j+1) || (module1==j+1 && module2==j)) fhReSameSectorEMCALMod[i]->Fill(pt,m) ;
}
//In case we want only pairs in same (super) module, check their origin.
Bool_t ok = kTRUE;
if(fSameSM && module1!=module2) ok=kFALSE;
- if(ok){
-
+ if(ok)
+ {
//Check if one of the clusters comes from a conversion
- if(fCheckConversion){
+ if(fCheckConversion)
+ {
if (p1->IsTagged() && p2->IsTagged()) fhReConv2->Fill(pt,m);
else if(p1->IsTagged() || p2->IsTagged()) fhReConv ->Fill(pt,m);
}
for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
if((p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)) && (p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton))){
for(Int_t iasym=0; iasym < fNAsymCuts; iasym++){
- if(a < fAsymCuts[iasym]){
+ if(a < fAsymCuts[iasym])
+ {
Int_t index = ((curCentrBin*fNPIDBits)+ipid)*fNAsymCuts + iasym;
- //printf("index %d :(cen %d * nPID %d + ipid %d)*nasym %d + iasym %d\n",index,curCentrBin,fNPIDBits,ipid,fNAsymCuts,iasym);
+ //printf("index %d :(cen %d * nPID %d + ipid %d)*nasym %d + iasym %d - max index %d\n",index,curCentrBin,fNPIDBits,ipid,fNAsymCuts,iasym, curCentrBin*fNPIDBits*fNAsymCuts);
+
+ if(index < 0 || index >= ncentr*fNPIDBits*fNAsymCuts) continue ;
+
fhRe1 [index]->Fill(pt,m);
if(fMakeInvPtPlots)fhReInvPt1[index]->Fill(pt,m,1./pt) ;
if(fFillBadDistHisto){
}// pid bit loop
//Fill histograms with opening angle
- if(fFillAngleHisto){
+ if(fFillAngleHisto)
+ {
fhRealOpeningAngle ->Fill(pt,angle);
fhRealCosOpeningAngle->Fill(pt,TMath::Cos(angle));
}
//Fill histograms with pair assymmetry
- if(fFillAsymmetryHisto){
+ if(fFillAsymmetryHisto)
+ {
fhRePtAsym->Fill(pt,a);
if(m > 0.10 && m < 0.17) fhRePtAsymPi0->Fill(pt,a);
if(m > 0.45 && m < 0.65) fhRePtAsymEta->Fill(pt,a);
//-------------------------------------------------------
Int_t ncell1 = 0;
Int_t ncell2 = 0;
- if(fMultiCutAna || (IsDataMC() && fMultiCutAnaSim)){
-
+ if(fMultiCutAna || (IsDataMC() && fMultiCutAnaSim))
+ {
AliVEvent * event = GetReader()->GetInputEvent();
if(event){
- for(Int_t iclus = 0; iclus < event->GetNumberOfCaloClusters(); iclus++){
+ for(Int_t iclus = 0; iclus < event->GetNumberOfCaloClusters(); iclus++)
+ {
AliVCluster *cluster = event->GetCaloCluster(iclus);
Bool_t is = kFALSE;
}
//-----------------------
- //Multi cuts analysis
+ //Multi cuts analysis
//-----------------------
- if(fMultiCutAna){
+ if(fMultiCutAna)
+ {
//Histograms for different PID bits selection
for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
//-------------------------------------------------------------
// Mixing
//-------------------------------------------------------------
- if(fDoOwnMix){
- //printf("Cen bin %d, RP bin %d, e aver %f, mult %d\n",curCentrBin,curRPBin, eClusTot, nPhot);
+ if(DoOwnMix())
+ {
//Recover events in with same characteristics as the current event
- TList * evMixList=fEventsList[curCentrBin*GetNZvertBin()*GetNRPBin()+curZvertBin*GetNRPBin()+curRPBin] ;
+
+ //Check that the bin exists, if not (bad determination of RP, centrality or vz bin) do nothing
+ if(eventbin < 0) return ;
+
+ TList * evMixList=fEventsList[eventbin] ;
+
+ if(!evMixList)
+ {
+ printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mix event list not available, bin %d \n",eventbin);
+ return;
+ }
+
Int_t nMixed = evMixList->GetSize() ;
- for(Int_t ii=0; ii<nMixed; ii++){
+ for(Int_t ii=0; ii<nMixed; ii++)
+ {
TClonesArray* ev2= (TClonesArray*) (evMixList->At(ii));
Int_t nPhot2=ev2->GetEntriesFast() ;
Double_t m = -999;
if(GetDebug() > 1)
- printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed event %d photon entries %d, centrality bin %d\n", ii, nPhot2, curCentrBin);
-
+ printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed event %d photon entries %d, centrality bin %d\n", ii, nPhot2, GetEventCentralityBin());
+
+ fhEventMixBin->Fill(eventbin) ;
+
//---------------------------------
//First loop on photons/clusters
//---------------------------------
}
//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))){
- for(Int_t iasym=0; iasym < fNAsymCuts; iasym++){
- if(a < fAsymCuts[iasym]){
+ if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton)))
+ {
+ for(Int_t iasym=0; iasym < fNAsymCuts; iasym++)
+ {
+ if(a < fAsymCuts[iasym])
+ {
Int_t index = ((curCentrBin*fNPIDBits)+ipid)*fNAsymCuts + iasym;
+
+ if(index < 0 || index >= ncentr*fNPIDBits*fNAsymCuts) continue ;
+
fhMi1 [index]->Fill(pt,m) ;
if(fMakeInvPtPlots)fhMiInvPt1[index]->Fill(pt,m,1./pt) ;
- if(fFillBadDistHisto){
- if(p1->DistToBad()>0 && p2->DistToBad()>0){
+ if(fFillBadDistHisto)
+ {
+ if(p1->DistToBad()>0 && p2->DistToBad()>0)
+ {
fhMi2 [index]->Fill(pt,m) ;
if(fMakeInvPtPlots)fhMiInvPt2[index]->Fill(pt,m,1./pt) ;
- if(p1->DistToBad()>1 && p2->DistToBad()>1){
+ if(p1->DistToBad()>1 && p2->DistToBad()>1)
+ {
fhMi3 [index]->Fill(pt,m) ;
if(fMakeInvPtPlots)fhMiInvPt3[index]->Fill(pt,m,1./pt) ;
}
}
-//____________________________________________________________________________________________________________________________________________________
+//________________________________________________________________________
Int_t AliAnaPi0::GetEventIndex(AliAODPWG4Particle * part, Double_t * vert)
{
// retieves the event index and checks the vertex
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff