/************************************************************************** * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * * * * Author: The ALICE Off-line Project. * * Contributors are mentioned in the code where appropriate. * * * * Permission to use, copy, modify and distribute this software and its * * documentation strictly for non-commercial purposes is hereby granted * * without fee, provided that the above copyright notice appears in all * * copies and that both the copyright notice and this permission notice * * appear in the supporting documentation. The authors make no claims * * about the suitability of this software for any purpose. It is * * provided "as is" without express or implied warranty. * **************************************************************************/ //_________________________________________________________________________ // Class to collect two-photon invariant mass distributions for // extracting raw pi0 yield. // Input is produced by AliAnaPhoton (or any other analysis producing output AliAODPWG4Particles), // it will do nothing if executed alone // //-- Author: Dmitri Peressounko (RRC "KI") //-- Adapted to CaloTrackCorr frame by Lamia Benhabib (SUBATECH) //-- and Gustavo Conesa (INFN-Frascati) //_________________________________________________________________________ // --- ROOT system --- #include "TH3.h" #include "TH2F.h" //#include "Riostream.h" #include "TCanvas.h" #include "TPad.h" #include "TROOT.h" #include "TClonesArray.h" #include "TObjString.h" #include "TDatabasePDG.h" //---- AliRoot system ---- #include "AliAnaPi0.h" #include "AliCaloTrackReader.h" #include "AliCaloPID.h" #include "AliStack.h" #include "AliFiducialCut.h" #include "TParticle.h" #include "AliVEvent.h" #include "AliESDCaloCluster.h" #include "AliESDEvent.h" #include "AliAODEvent.h" #include "AliNeutralMesonSelection.h" #include "AliMixedEvent.h" #include "AliAODMCParticle.h" // --- Detectors --- #include "AliPHOSGeoUtils.h" #include "AliEMCALGeometry.h" ClassImp(AliAnaPi0) //______________________________________________________ AliAnaPi0::AliAnaPi0() : AliAnaCaloTrackCorrBaseClass(), fEventsList(0x0), fNModules(22), 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), fFillBadDistHisto(kFALSE), fFillSSCombinations(kFALSE), fFillAngleHisto(kFALSE), fFillAsymmetryHisto(kFALSE), fFillOriginHisto(0), fFillArmenterosThetaStar(0), fCheckAccInSector(kFALSE), fPhotonMom1(), fPhotonMom1Boost(), fPhotonMom2(), fPi0Mom(), fProdVertex(), //Histograms fhAverTotECluster(0), fhAverTotECell(0), fhAverTotECellvsCluster(0), fhEDensityCluster(0), fhEDensityCell(0), fhEDensityCellvsCluster(0), fhReMod(0x0), fhReSameSideEMCALMod(0x0), fhReSameSectorEMCALMod(0x0), fhReDiffPHOSMod(0x0), fhMiMod(0x0), fhMiSameSideEMCALMod(0x0), fhMiSameSectorEMCALMod(0x0), fhMiDiffPHOSMod(0x0), fhReConv(0x0), fhMiConv(0x0), fhReConv2(0x0), fhMiConv2(0x0), fhRe1(0x0), fhMi1(0x0), fhRe2(0x0), fhMi2(0x0), fhRe3(0x0), fhMi3(0x0), fhReInvPt1(0x0), fhMiInvPt1(0x0), fhReInvPt2(0x0), fhMiInvPt2(0x0), fhReInvPt3(0x0), fhMiInvPt3(0x0), fhRePtNCellAsymCuts(0x0), fhMiPtNCellAsymCuts(0x0), fhRePtNCellAsymCutsSM(), fhRePIDBits(0x0), fhRePtMult(0x0), fhReSS(), fhRePtAsym(0x0), fhRePtAsymPi0(0x0), fhRePtAsymEta(0x0), fhEventBin(0), fhEventMixBin(0), fhCentrality(0x0), fhCentralityNoPair(0x0), fhEventPlaneResolution(0x0), fhRealOpeningAngle(0x0), fhRealCosOpeningAngle(0x0), fhMixedOpeningAngle(0x0), fhMixedCosOpeningAngle(0x0), // MC histograms fhPrimPi0E(0x0), fhPrimPi0Pt(0x0), fhPrimPi0AccE(0x0), fhPrimPi0AccPt(0x0), fhPrimPi0Y(0x0), fhPrimPi0AccY(0x0), fhPrimPi0Yeta(0x0), fhPrimPi0YetaYcut(0x0), fhPrimPi0AccYeta(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), fhPrimEtaYeta(0x0), fhPrimEtaYetaYcut(0x0), fhPrimEtaAccYeta(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), fhMCPi0ProdVertex(0), fhMCEtaProdVertex(0), fhPrimPi0ProdVertex(0), fhPrimEtaProdVertex(0), fhReMCFromConversion(0), fhReMCFromNotConversion(0), fhReMCFromMixConversion(0), fhCosThStarPrimPi0(0), fhCosThStarPrimEta(0)//, { //Default Ctor InitParameters(); for(Int_t i = 0; i < 4; i++) { fhArmPrimEta[i] = 0; fhArmPrimPi0[i] = 0; } } //_____________________ AliAnaPi0::~AliAnaPi0() { // Remove event containers if(DoOwnMix() && fEventsList) { for(Int_t ic=0; icDelete() ; delete fEventsList[bin] ; } } } delete[] fEventsList; } } //______________________________ void AliAnaPi0::InitParameters() { //Init parameters when first called the analysis //Set default parameters SetInputAODName("PWG4Particle"); AddToHistogramsName("AnaPi0_"); fUseAngleCut = kFALSE; fUseAngleEDepCut = kFALSE; fAngleCut = 0.; fAngleMaxCut = TMath::Pi(); fMultiCutAna = kFALSE; fNPtCuts = 3; fPtCuts[0] = 0.; fPtCuts[1] = 0.3; fPtCuts[2] = 0.5; for(Int_t i = fNPtCuts; i < 10; i++)fPtCuts[i] = 0.; 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; fCellNCuts[0] = 0; fCellNCuts[1] = 1; fCellNCuts[2] = 2; for(Int_t i = fNCellNCuts; i < 10; i++)fCellNCuts[i] = 0; fNPIDBits = 2; 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; } //_______________________________________ TObjString * AliAnaPi0::GetAnalysisCuts() { //Save parameters used for analysis TString parList ; //this will be list of parameters used for this analysis. const Int_t buffersize = 255; char onePar[buffersize] ; snprintf(onePar,buffersize,"--- AliAnaPi0 ---\n") ; parList+=onePar ; 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",GetNMaxEvMix()) ; 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) ; for(Int_t i = 0; i < fNAsymCuts; i++) snprintf(onePar,buffersize,"%s %2.2f;",onePar,fAsymCuts[i]); parList+=onePar ; snprintf(onePar,buffersize," PID selection bits: n = %d, PID bit =\n",fNPIDBits) ; for(Int_t i = 0; i < fNPIDBits; i++) snprintf(onePar,buffersize,"%s %d;",onePar,fPIDBits[i]); parList+=onePar ; snprintf(onePar,buffersize,"Cuts: \n") ; parList+=onePar ; snprintf(onePar,buffersize,"Z vertex position: -%f < z < %f \n",GetZvertexCut(),GetZvertexCut()) ; parList+=onePar ; snprintf(onePar,buffersize,"Calorimeter: %s \n",GetCalorimeterString().Data()) ; parList+=onePar ; snprintf(onePar,buffersize,"Number of modules: %d \n",fNModules) ; parList+=onePar ; if(fMultiCutAna){ snprintf(onePar, buffersize," pT cuts: n = %d, pt > ",fNPtCuts) ; for(Int_t i = 0; i < fNPtCuts; i++) snprintf(onePar,buffersize,"%s %2.2f;",onePar,fPtCuts[i]); parList+=onePar ; snprintf(onePar,buffersize, " N cell in cluster cuts: n = %d, nCell > ",fNCellNCuts) ; for(Int_t i = 0; i < fNCellNCuts; i++) snprintf(onePar,buffersize,"%s %d;",onePar,fCellNCuts[i]); parList+=onePar ; } return new TObjString(parList) ; } //_________________________________________ TList * AliAnaPi0::GetCreateOutputObjects() { // Create histograms to be saved in output file and // store them in fOutputContainer // Init the number of modules, set in the class AliCalorimeterUtils fNModules = GetCaloUtils()->GetNumberOfSuperModulesUsed(); if(GetCalorimeter()==kPHOS && fNModules > 4) fNModules = 4; //create event containers fEventsList = new TList*[GetNCentrBin()*GetNZvertBin()*GetNRPBin()] ; for(Int_t ic=0; icSetOwner(kFALSE); } } } TList * outputContainer = new TList() ; outputContainer->SetName(GetName()); fhReMod = new TH2F*[fNModules] ; fhMiMod = new TH2F*[fNModules] ; if(GetCalorimeter() == kPHOS) { fhReDiffPHOSMod = new TH2F*[fNModules] ; fhMiDiffPHOSMod = new TH2F*[fNModules] ; } else { fhReSameSectorEMCALMod = new TH2F*[fNModules/2] ; fhReSameSideEMCALMod = new TH2F*[fNModules-2] ; fhMiSameSectorEMCALMod = new TH2F*[fNModules/2] ; fhMiSameSideEMCALMod = new TH2F*[fNModules-2] ; } fhRe1 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ; fhMi1 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ; 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) { fhReInvPt1 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ; fhMiInvPt1 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ; if(fFillBadDistHisto){ fhReInvPt2 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ; fhReInvPt3 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ; fhMiInvPt2 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ; fhMiInvPt3 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ; } } const Int_t buffersize = 255; char key[buffersize] ; char title[buffersize] ; Int_t nptbins = GetHistogramRanges()->GetHistoPtBins(); Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins(); Int_t netabins = GetHistogramRanges()->GetHistoEtaBins(); Float_t ptmax = GetHistogramRanges()->GetHistoPtMax(); Float_t phimax = GetHistogramRanges()->GetHistoPhiMax(); Float_t etamax = GetHistogramRanges()->GetHistoEtaMax(); Float_t ptmin = GetHistogramRanges()->GetHistoPtMin(); Float_t phimin = GetHistogramRanges()->GetHistoPhiMin(); Float_t etamin = GetHistogramRanges()->GetHistoEtaMin(); Int_t nmassbins = GetHistogramRanges()->GetHistoMassBins(); Int_t nasymbins = GetHistogramRanges()->GetHistoAsymmetryBins(); Float_t massmax = GetHistogramRanges()->GetHistoMassMax(); Float_t asymmax = GetHistogramRanges()->GetHistoAsymmetryMax(); Float_t massmin = GetHistogramRanges()->GetHistoMassMin(); Float_t asymmin = GetHistogramRanges()->GetHistoAsymmetryMin(); Int_t ntrmbins = GetHistogramRanges()->GetHistoTrackMultiplicityBins(); Int_t ntrmmax = GetHistogramRanges()->GetHistoTrackMultiplicityMax(); Int_t ntrmmin = GetHistogramRanges()->GetHistoTrackMultiplicityMin(); if(fCheckConversion) { fhReConv = new TH2F("hReConv","Real Pair with one recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ; fhReConv->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhReConv->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})"); outputContainer->Add(fhReConv) ; fhReConv2 = new TH2F("hReConv2","Real Pair with 2 recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ; fhReConv2->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhReConv2->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})"); outputContainer->Add(fhReConv2) ; if(DoOwnMix()) { fhMiConv = new TH2F("hMiConv","Mixed Pair with one recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ; fhMiConv->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhMiConv->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})"); outputContainer->Add(fhMiConv) ; fhMiConv2 = new TH2F("hMiConv2","Mixed Pair with 2 recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ; fhMiConv2->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhMiConv2->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})"); outputContainer->Add(fhMiConv2) ; } } for(Int_t ic=0; icSetXTitle("#it{p}_{T} (GeV/#it{c})"); fhRe1[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})"); //printf("name: %s\n ",fhRe1[index]->GetName()); outputContainer->Add(fhRe1[index]) ; 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 #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2", ic,fPIDBits[ipid], fAsymCuts[iasym]) ; fhRe2[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ; fhRe2[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhRe2[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})"); outputContainer->Add(fhRe2[index]) ; //Distance to bad module 3 snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ; snprintf(title, buffersize,"Real #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3", ic,fPIDBits[ipid], fAsymCuts[iasym]) ; fhRe3[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ; fhRe3[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhRe3[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})"); outputContainer->Add(fhRe3[index]) ; } //Inverse pT 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 #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1", ic,fPIDBits[ipid], fAsymCuts[iasym]) ; fhReInvPt1[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ; fhReInvPt1[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhReInvPt1[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})"); outputContainer->Add(fhReInvPt1[index]) ; if(fFillBadDistHisto){ //Distance to bad module 2 snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ; snprintf(title, buffersize,"Real #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2", ic,fPIDBits[ipid], fAsymCuts[iasym]) ; fhReInvPt2[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ; fhReInvPt2[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhReInvPt2[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})"); outputContainer->Add(fhReInvPt2[index]) ; //Distance to bad module 3 snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ; snprintf(title, buffersize,"Real #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3", ic,fPIDBits[ipid], fAsymCuts[iasym]) ; fhReInvPt3[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ; fhReInvPt3[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhReInvPt3[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})"); outputContainer->Add(fhReInvPt3[index]) ; } } 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 #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1", ic,fPIDBits[ipid], fAsymCuts[iasym]) ; fhMi1[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ; fhMi1[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhMi1[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})"); outputContainer->Add(fhMi1[index]) ; if(fFillBadDistHisto){ //Distance to bad module 2 snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ; snprintf(title, buffersize,"Mixed #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2", ic,fPIDBits[ipid], fAsymCuts[iasym]) ; fhMi2[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ; fhMi2[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhMi2[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})"); outputContainer->Add(fhMi2[index]) ; //Distance to bad module 3 snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ; snprintf(title, buffersize,"Mixed #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3", ic,fPIDBits[ipid], fAsymCuts[iasym]) ; fhMi3[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ; fhMi3[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhMi3[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})"); outputContainer->Add(fhMi3[index]) ; } //Inverse pT 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 #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1", ic,fPIDBits[ipid], fAsymCuts[iasym]) ; fhMiInvPt1[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ; fhMiInvPt1[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhMiInvPt1[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})"); outputContainer->Add(fhMiInvPt1[index]) ; if(fFillBadDistHisto){ //Distance to bad module 2 snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ; snprintf(title, buffersize,"Mixed #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2", ic,fPIDBits[ipid], fAsymCuts[iasym]) ; fhMiInvPt2[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ; fhMiInvPt2[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhMiInvPt2[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})"); outputContainer->Add(fhMiInvPt2[index]) ; //Distance to bad module 3 snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ; snprintf(title, buffersize,"Mixed #it{M}_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f,dist bad 3", ic,fPIDBits[ipid], fAsymCuts[iasym]) ; fhMiInvPt3[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ; fhMiInvPt3[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhMiInvPt3[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})"); outputContainer->Add(fhMiInvPt3[index]) ; } } } } } } if(fFillAsymmetryHisto) { fhRePtAsym = new TH2F("hRePtAsym","#it{Asymmetry} vs #it{p}_{T} , for pairs",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ; fhRePtAsym->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhRePtAsym->SetYTitle("#it{Asymmetry}"); outputContainer->Add(fhRePtAsym); fhRePtAsymPi0 = new TH2F("hRePtAsymPi0","#it{Asymmetry} vs #it{p}_{T} , for pairs close to #pi^{0} mass",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ; fhRePtAsymPi0->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhRePtAsymPi0->SetYTitle("Asymmetry"); outputContainer->Add(fhRePtAsymPi0); fhRePtAsymEta = new TH2F("hRePtAsymEta","#it{Asymmetry} vs #it{p}_{T} , for pairs close to #eta mass",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ; fhRePtAsymEta->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhRePtAsymEta->SetYTitle("#it{Asymmetry}"); outputContainer->Add(fhRePtAsymEta); } if(fMultiCutAna) { fhRePIDBits = new TH2F*[fNPIDBits]; for(Int_t ipid=0; ipidSetXTitle("#it{p}_{T} (GeV/#it{c})"); fhRePIDBits[ipid]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})"); outputContainer->Add(fhRePIDBits[ipid]) ; }// pid bit loop fhRePtNCellAsymCuts = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts]; fhMiPtNCellAsymCuts = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts]; if(fFillSMCombinations) for(Int_t iSM = 0; iSM < fNModules; iSM++) fhRePtNCellAsymCutsSM[iSM] = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts]; for(Int_t ipt=0; ipt%2.2f, ncell>%d and asym >%1.2f ",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ; Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym; //printf("ipt %d, icell %d, iassym %d, index %d\n",ipt, icell, iasym, index); fhRePtNCellAsymCuts[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ; fhRePtNCellAsymCuts[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhRePtNCellAsymCuts[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})"); outputContainer->Add(fhRePtNCellAsymCuts[index]) ; snprintf(key, buffersize,"hMi_pt%d_cell%d_asym%d",ipt,icell,iasym) ; snprintf(title, buffersize,"Mixed #it{M}_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ; fhMiPtNCellAsymCuts[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ; fhMiPtNCellAsymCuts[index]->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhMiPtNCellAsymCuts[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})"); outputContainer->Add(fhMiPtNCellAsymCuts[index]) ; 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 #it{M}_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f, SM %d ", fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym],iSM) ; fhRePtNCellAsymCutsSM[iSM][index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ; fhRePtNCellAsymCutsSM[iSM][index]->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhRePtNCellAsymCutsSM[iSM][index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})"); outputContainer->Add(fhRePtNCellAsymCutsSM[iSM][index]) ; } } } } } if(ntrmbins!=0) { fhRePtMult = new TH3F*[fNAsymCuts] ; for(Int_t iasym = 0; iasymSetXTitle("#it{p}_{T} (GeV/#it{c})"); fhRePtMult[iasym]->SetYTitle("Track multiplicity"); fhRePtMult[iasym]->SetZTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})"); outputContainer->Add(fhRePtMult[iasym]) ; } } }// multi cuts analysis if(fFillSSCombinations) { fhReSS[0] = new TH2F("hRe_SS_Tight"," 0.01 < #lambda_{0}^{2} < 0.4", nptbins,ptmin,ptmax,nmassbins,massmin,massmax); fhReSS[0]->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhReSS[0]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})"); outputContainer->Add(fhReSS[0]) ; fhReSS[1] = new TH2F("hRe_SS_Loose"," #lambda_{0}^{2} > 0.4", nptbins,ptmin,ptmax,nmassbins,massmin,massmax); fhReSS[1]->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhReSS[1]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})"); outputContainer->Add(fhReSS[1]) ; fhReSS[2] = new TH2F("hRe_SS_Both"," cluster_{1} #lambda_{0}^{2} > 0.4; cluster_{2} 0.01 < #lambda_{0}^{2} < 0.4", nptbins,ptmin,ptmax,nmassbins,massmin,massmax); fhReSS[2]->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhReSS[2]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})"); outputContainer->Add(fhReSS[2]) ; } 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) ; } if(GetNCentrBin()>1) { fhCentrality=new TH1F("hCentralityBin","Number of events in centrality bin",GetNCentrBin(),0.,1.*GetNCentrBin()) ; fhCentrality->SetXTitle("Centrality bin"); outputContainer->Add(fhCentrality) ; fhCentralityNoPair=new TH1F("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 && 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) { fhRealOpeningAngle = new TH2F ("hRealOpeningAngle","Angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,300,0,TMath::Pi()); fhRealOpeningAngle->SetYTitle("#theta(rad)"); fhRealOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)"); outputContainer->Add(fhRealOpeningAngle) ; fhRealCosOpeningAngle = new TH2F ("hRealCosOpeningAngle","Cosinus of angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,0,1); fhRealCosOpeningAngle->SetYTitle("cos (#theta) "); fhRealCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)"); outputContainer->Add(fhRealCosOpeningAngle) ; 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)"); fhMixedOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)"); outputContainer->Add(fhMixedOpeningAngle) ; fhMixedCosOpeningAngle = new TH2F ("hMixedCosOpeningAngle","Cosinus of angle between all #gamma pair vs E_{#pi^{0}}, Mixed pairs",nptbins,ptmin,ptmax,100,0,1); fhMixedCosOpeningAngle->SetYTitle("cos (#theta) "); fhMixedCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)"); outputContainer->Add(fhMixedCosOpeningAngle) ; } } //Histograms filled only if MC data is requested if(IsDataMC()) { fhReMCFromConversion = new TH2F("hReMCFromConversion","Invariant mass of 2 clusters originated in conversions", nptbins,ptmin,ptmax,nmassbins,massmin,massmax); fhReMCFromConversion->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhReMCFromConversion->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})"); outputContainer->Add(fhReMCFromConversion) ; fhReMCFromNotConversion = new TH2F("hReMCNotFromConversion","Invariant mass of 2 clusters not originated in conversions", nptbins,ptmin,ptmax,nmassbins,massmin,massmax); fhReMCFromNotConversion->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhReMCFromNotConversion->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})"); outputContainer->Add(fhReMCFromNotConversion) ; fhReMCFromMixConversion = new TH2F("hReMCFromMixConversion","Invariant mass of 2 clusters one from conversion and the other not", nptbins,ptmin,ptmax,nmassbins,massmin,massmax); fhReMCFromMixConversion->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhReMCFromMixConversion->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})"); outputContainer->Add(fhReMCFromMixConversion) ; //Pi0 fhPrimPi0E = new TH1F("hPrimPi0E","Primary #pi^{0} E, |#it{Y}|<1",nptbins,ptmin,ptmax) ; fhPrimPi0AccE = new TH1F("hPrimPi0AccE","Primary #pi^{0} #it{E} with both photons in acceptance",nptbins,ptmin,ptmax) ; fhPrimPi0E ->SetXTitle("#it{E} (GeV)"); fhPrimPi0AccE->SetXTitle("#it{E} (GeV)"); outputContainer->Add(fhPrimPi0E) ; outputContainer->Add(fhPrimPi0AccE) ; fhPrimPi0Pt = new TH1F("hPrimPi0Pt","Primary #pi^{0} #it{p}_{T} , |#it{Y}|<1",nptbins,ptmin,ptmax) ; fhPrimPi0AccPt = new TH1F("hPrimPi0AccPt","Primary #pi^{0} #it{p}_{T} with both photons in acceptance",nptbins,ptmin,ptmax) ; fhPrimPi0Pt ->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhPrimPi0AccPt->SetXTitle("#it{p}_{T} (GeV/#it{c})"); outputContainer->Add(fhPrimPi0Pt) ; outputContainer->Add(fhPrimPi0AccPt) ; Int_t netabinsopen = TMath::Nint(netabins*4/(etamax-etamin)); fhPrimPi0Y = new TH2F("hPrimPi0Rapidity","Rapidity of primary #pi^{0}",nptbins,ptmin,ptmax,netabinsopen,-2, 2) ; fhPrimPi0Y ->SetYTitle("#it{Rapidity}"); fhPrimPi0Y ->SetXTitle("#it{p}_{T} (GeV/#it{c})"); outputContainer->Add(fhPrimPi0Y) ; fhPrimPi0Yeta = new TH2F("hPrimPi0PseudoRapidity","PseudoRapidity of primary #pi^{0}",nptbins,ptmin,ptmax,netabinsopen,-2, 2) ; fhPrimPi0Yeta ->SetYTitle("#eta"); fhPrimPi0Yeta ->SetXTitle("#it{p}_{T} (GeV/#it{c})"); outputContainer->Add(fhPrimPi0Yeta) ; fhPrimPi0YetaYcut = new TH2F("hPrimPi0PseudoRapidityYcut","PseudoRapidity of primary #pi^{0}, |#it{Y}|<1",nptbins,ptmin,ptmax,netabinsopen,-2, 2) ; fhPrimPi0YetaYcut ->SetYTitle("#eta"); fhPrimPi0YetaYcut ->SetXTitle("#it{p}_{T} (GeV/#it{c})"); outputContainer->Add(fhPrimPi0YetaYcut) ; fhPrimPi0AccY = new TH2F("hPrimPi0AccRapidity","Rapidity of primary #pi^{0} with accepted daughters",nptbins,ptmin,ptmax,netabins,etamin,etamax) ; fhPrimPi0AccY->SetYTitle("Rapidity"); fhPrimPi0AccY->SetXTitle("#it{p}_{T} (GeV/#it{c})"); outputContainer->Add(fhPrimPi0AccY) ; fhPrimPi0AccYeta = new TH2F("hPrimPi0AccPseudoRapidity","PseudoRapidity of primary #pi^{0} with accepted daughters",nptbins,ptmin,ptmax,netabins,etamin,etamax) ; fhPrimPi0AccYeta ->SetYTitle("#eta"); fhPrimPi0AccYeta ->SetXTitle("#it{p}_{T} (GeV/#it{c})"); outputContainer->Add(fhPrimPi0AccYeta) ; Int_t nphibinsopen = TMath::Nint(nphibins*TMath::TwoPi()/(phimax-phimin)); fhPrimPi0Phi = new TH2F("hPrimPi0Phi","#phi of primary #pi^{0}, |#it{Y}|<1",nptbins,ptmin,ptmax,nphibinsopen,0,360) ; fhPrimPi0Phi->SetYTitle("#phi (deg)"); fhPrimPi0Phi->SetXTitle("#it{p}_{T} (GeV/#it{c})"); outputContainer->Add(fhPrimPi0Phi) ; fhPrimPi0AccPhi = new TH2F("hPrimPi0AccPhi","#phi of primary #pi^{0} with accepted daughters",nptbins,ptmin,ptmax, nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ; fhPrimPi0AccPhi->SetYTitle("#phi (deg)"); fhPrimPi0AccPhi->SetXTitle("#it{p}_{T} (GeV/#it{c})"); outputContainer->Add(fhPrimPi0AccPhi) ; fhPrimPi0PtCentrality = new TH2F("hPrimPi0PtCentrality","Primary #pi^{0} #it{p}_{T} vs reco centrality, |#it{Y}|<1", nptbins,ptmin,ptmax, 100, 0, 100) ; fhPrimPi0AccPtCentrality = new TH2F("hPrimPi0AccPtCentrality","Primary #pi^{0} with both photons in acceptance #it{p}_{T} vs reco centrality", nptbins,ptmin,ptmax, 100, 0, 100) ; fhPrimPi0PtCentrality ->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhPrimPi0AccPtCentrality->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhPrimPi0PtCentrality ->SetYTitle("Centrality"); fhPrimPi0AccPtCentrality->SetYTitle("Centrality"); outputContainer->Add(fhPrimPi0PtCentrality) ; outputContainer->Add(fhPrimPi0AccPtCentrality) ; fhPrimPi0PtEventPlane = new TH2F("hPrimPi0PtEventPlane","Primary #pi^{0} #it{p}_{T} vs reco event plane angle, |#it{Y}|<1", nptbins,ptmin,ptmax, 100, 0, TMath::Pi()) ; fhPrimPi0AccPtEventPlane = new TH2F("hPrimPi0AccPtEventPlane","Primary #pi^{0} with both photons in acceptance #it{p}_{T} vs reco event plane angle", nptbins,ptmin,ptmax, 100, 0, TMath::Pi()) ; fhPrimPi0PtEventPlane ->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhPrimPi0AccPtEventPlane->SetXTitle("#it{p}_{T} (GeV/#it{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 #it{E} with both photons in acceptance",nptbins,ptmin,ptmax) ; fhPrimEtaE ->SetXTitle("#it{E} (GeV)"); fhPrimEtaAccE->SetXTitle("#it{E} (GeV)"); outputContainer->Add(fhPrimEtaE) ; outputContainer->Add(fhPrimEtaAccE) ; fhPrimEtaPt = new TH1F("hPrimEtaPt","Primary #eta #it{p}_{T}",nptbins,ptmin,ptmax) ; fhPrimEtaAccPt = new TH1F("hPrimEtaAccPt","Primary eta #it{p}_{T} with both photons in acceptance",nptbins,ptmin,ptmax) ; fhPrimEtaPt ->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhPrimEtaAccPt->SetXTitle("#it{p}_{T} (GeV/#it{c})"); outputContainer->Add(fhPrimEtaPt) ; outputContainer->Add(fhPrimEtaAccPt) ; fhPrimEtaY = new TH2F("hPrimEtaRapidity","Rapidity of primary #eta",nptbins,ptmin,ptmax,netabinsopen,-2, 2) ; fhPrimEtaY->SetYTitle("#it{Rapidity}"); fhPrimEtaY->SetXTitle("#it{p}_{T} (GeV/#it{c})"); outputContainer->Add(fhPrimEtaY) ; fhPrimEtaYeta = new TH2F("hPrimEtaPseudoRapidityEta","PseudoRapidity of primary #eta",nptbins,ptmin,ptmax,netabinsopen,-2, 2) ; fhPrimEtaYeta->SetYTitle("#it{Rapidity}"); fhPrimEtaYeta->SetXTitle("#it{p}_{T} (GeV/#it{c})"); outputContainer->Add(fhPrimEtaYeta) ; fhPrimEtaYetaYcut = new TH2F("hPrimEtaPseudoRapidityEtaYcut","PseudoRapidity of primary #eta, |#it{Y}|<1",nptbins,ptmin,ptmax,netabinsopen,-2, 2) ; fhPrimEtaYetaYcut->SetYTitle("#it{Pseudorapidity}"); fhPrimEtaYetaYcut->SetXTitle("#it{p}_{T} (GeV/#it{c})"); outputContainer->Add(fhPrimEtaYetaYcut) ; fhPrimEtaAccY = new TH2F("hPrimEtaAccRapidity","Rapidity of primary #eta",nptbins,ptmin,ptmax, netabins,etamin,etamax) ; fhPrimEtaAccY->SetYTitle("#it{Rapidity}"); fhPrimEtaAccY->SetXTitle("#it{p}_{T} (GeV/#it{c})"); outputContainer->Add(fhPrimEtaAccY) ; fhPrimEtaAccYeta = new TH2F("hPrimEtaAccPseudoRapidity","PseudoRapidity of primary #eta",nptbins,ptmin,ptmax, netabins,etamin,etamax) ; fhPrimEtaAccYeta->SetYTitle("#it{Pseudorapidity}"); fhPrimEtaAccYeta->SetXTitle("#it{p}_{T} (GeV/#it{c})"); outputContainer->Add(fhPrimEtaAccYeta) ; fhPrimEtaPhi = new TH2F("hPrimEtaPhi","Azimuthal of primary #eta",nptbins,ptmin,ptmax, nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ; fhPrimEtaPhi->SetYTitle("#phi (deg)"); fhPrimEtaPhi->SetXTitle("#it{p}_{T} (GeV/#it{c})"); outputContainer->Add(fhPrimEtaPhi) ; fhPrimEtaAccPhi = new TH2F("hPrimEtaAccPhi","Azimuthal of primary #eta with accepted daughters",nptbins,ptmin,ptmax, nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ; fhPrimEtaAccPhi->SetYTitle("#phi (deg)"); fhPrimEtaAccPhi->SetXTitle("#it{p}_{T} (GeV/#it{c})"); outputContainer->Add(fhPrimEtaAccPhi) ; fhPrimEtaPtCentrality = new TH2F("hPrimEtaPtCentrality","Primary #eta #it{p}_{T} vs reco centrality, |#it{Y}|<1",nptbins,ptmin,ptmax, 100, 0, 100) ; fhPrimEtaAccPtCentrality = new TH2F("hPrimEtaAccPtCentrality","Primary #eta with both photons in acceptance #it{p}_{T} vs reco centrality",nptbins,ptmin,ptmax, 100, 0, 100) ; fhPrimEtaPtCentrality ->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhPrimEtaAccPtCentrality->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhPrimEtaPtCentrality ->SetYTitle("Centrality"); fhPrimEtaAccPtCentrality->SetYTitle("Centrality"); outputContainer->Add(fhPrimEtaPtCentrality) ; outputContainer->Add(fhPrimEtaAccPtCentrality) ; fhPrimEtaPtEventPlane = new TH2F("hPrimEtaPtEventPlane","Primary #eta #it{p}_{T} vs reco event plane angle, |#it{Y}|<1",nptbins,ptmin,ptmax, 100, 0, TMath::Pi()) ; fhPrimEtaAccPtEventPlane = new TH2F("hPrimEtaAccPtEventPlane","Primary #eta with both #gamma_{decay} in acceptance #it{p}_{T} vs reco event plane angle",nptbins,ptmin,ptmax, 100, 0, TMath::Pi()) ; fhPrimEtaPtEventPlane ->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhPrimEtaAccPtEventPlane->SetXTitle("#it{p}_{T} (GeV/#it{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 #it{Asymmetry}, #it{p}_{T}>5 GeV/#it{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) ; 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) ; fhPrimEtaOpeningAngleAsym = new TH2F ("hPrimEtaOpeningAngleAsym","Angle between all primary #gamma pair vs #it{Asymmetry}, #it{p}_{T}>5 GeV/#it{c}",100,0,1,100,0,0.5); fhPrimEtaOpeningAngleAsym->SetXTitle("|#it{A}|=| (#it{E}_{1}-#it{E}_{2}) / (#it{E}_{1}+#it{E}_{2}) |"); fhPrimEtaOpeningAngleAsym->SetYTitle("#theta(rad)"); outputContainer->Add(fhPrimEtaOpeningAngleAsym) ; 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("#it{E}_{ #eta} (GeV)"); outputContainer->Add(fhPrimEtaCosOpeningAngle) ; } if(fFillOriginHisto) { //Prim origin //Pi0 fhPrimPi0PtOrigin = new TH2F("hPrimPi0PtOrigin","Primary #pi^{0} #it{p}_{T} vs origin",nptbins,ptmin,ptmax,11,0,11) ; fhPrimPi0PtOrigin->SetXTitle("#it{p}_{T} (GeV/#it{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 #pi^{0} #it{p}_{T} vs origin",nptbins,ptmin,ptmax,11,0,11) ; fhMCPi0PtOrigin->SetXTitle("#it{p}_{T} (GeV/#it{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 #pi^{0} #it{p}_{T} vs origin",nptbins,ptmin,ptmax,7,0,7) ; fhPrimEtaPtOrigin->SetXTitle("#it{p}_{T} (GeV/#it{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 #pi^{0} #it{p}_{T} vs origin",nptbins,ptmin,ptmax,7,0,7) ; fhMCEtaPtOrigin->SetXTitle("#it{p}_{T} (GeV/#it{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) ; fhMCPi0ProdVertex = new TH2F("hMCPi0ProdVertex","Selected reco pair from generated #pi^{0} #it{p}_{T} vs production vertex", 200,0.,20.,5000,0,500) ; fhMCPi0ProdVertex->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhMCPi0ProdVertex->SetYTitle("#it{R} (cm)"); outputContainer->Add(fhMCPi0ProdVertex) ; fhMCEtaProdVertex = new TH2F("hMCEtaProdVertex","Selected reco pair from generated #eta #it{p}_{T} vs production vertex", 200,0.,20.,5000,0,500) ; fhMCEtaProdVertex->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhMCEtaProdVertex->SetYTitle("#it{R} (cm)"); outputContainer->Add(fhMCEtaProdVertex) ; fhPrimPi0ProdVertex = new TH2F("hPrimPi0ProdVertex","generated #pi^{0} #it{p}_{T} vs production vertex", 200,0.,20.,5000,0,500) ; fhPrimPi0ProdVertex->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhPrimPi0ProdVertex->SetYTitle("#it{R} (cm)"); outputContainer->Add(fhPrimPi0ProdVertex) ; fhPrimEtaProdVertex = new TH2F("hPrimEtaProdVertex","generated #eta #it{p}_{T} vs production vertex", 200,0.,20.,5000,0,500) ; fhPrimEtaProdVertex->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhPrimEtaProdVertex->SetYTitle("#it{R} (cm)"); outputContainer->Add(fhPrimEtaProdVertex) ; for(Int_t i = 0; i<13; i++) { fhMCOrgMass[i] = new TH2F(Form("hMCOrgMass_%d",i),Form("#it{M} vs #it{p}_{T}, origin %d",i),nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ; fhMCOrgMass[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhMCOrgMass[i]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})"); outputContainer->Add(fhMCOrgMass[i]) ; fhMCOrgAsym[i]= new TH2F(Form("hMCOrgAsym_%d",i),Form("#it{Asymmetry} vs #it{p}_{T}, origin %d",i),nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ; fhMCOrgAsym[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhMCOrgAsym[i]->SetYTitle("A"); outputContainer->Add(fhMCOrgAsym[i]) ; fhMCOrgDeltaEta[i] = new TH2F(Form("hMCOrgDeltaEta_%d",i),Form("#Delta #eta of pair vs #it{p}_{T}, origin %d",i),nptbins,ptmin,ptmax,netabins,-1.4,1.4) ; fhMCOrgDeltaEta[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhMCOrgDeltaEta[i]->SetYTitle("#Delta #eta"); outputContainer->Add(fhMCOrgDeltaEta[i]) ; fhMCOrgDeltaPhi[i]= new TH2F(Form("hMCOrgDeltaPhi_%d",i),Form("#Delta #phi of pair vs #it{p}_{T}, origin %d",i),nptbins,ptmin,ptmax,nphibins,-0.7,0.7) ; fhMCOrgDeltaPhi[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhMCOrgDeltaPhi[i]->SetYTitle("#Delta #phi (rad)"); outputContainer->Add(fhMCOrgDeltaPhi[i]) ; } 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%2.2f, #it{N}^{cluster}_{cell}>%d and |#it{A}|>%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]), nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ; fhMCPi0MassPtRec[index]->SetXTitle("#it{p}_{T, reconstructed} (GeV/#it{c})"); fhMCPi0MassPtRec[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})"); outputContainer->Add(fhMCPi0MassPtRec[index]) ; fhMCPi0MassPtTrue[index] = new TH2F(Form("hMCPi0MassPtTrue_pt%d_cell%d_asym%d",ipt,icell,iasym), Form("Reconstructed #it{M} vs generated #it{p}_{T} of true #pi^{0} cluster pairs for #it{p}_{T} >%2.2f, #it{N}^{cluster}_{cell}>%d and |#it{A}|>%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]), nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ; fhMCPi0MassPtTrue[index]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})"); fhMCPi0MassPtTrue[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{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 #it{p}_{T} of true #pi^{0} cluster pairs, 0.01 < rec. mass < 0.17 MeV/#it{c}^{2} for #it{p}_{T} >%2.2f, #it{N}^{cluster}_{cell}>%d and |#it{A}|>%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]), nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ; fhMCPi0PtTruePtRec[index]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})"); fhMCPi0PtTruePtRec[index]->SetYTitle("#it{p}_{T, reconstructed} (GeV/#it{c})"); outputContainer->Add(fhMCPi0PtTruePtRec[index]) ; fhMCEtaMassPtRec[index] = new TH2F(Form("hMCEtaMassPtRec_pt%d_cell%d_asym%d",ipt,icell,iasym), Form("Reconstructed #it{M} vs reconstructed #it{p}_{T} of true #eta cluster pairs for #it{p}_{T} >%2.2f, #it{N}^{cluster}_{cell}>%d and |#it{A}|>%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]), nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ; fhMCEtaMassPtRec[index]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})"); fhMCEtaMassPtRec[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})"); outputContainer->Add(fhMCEtaMassPtRec[index]) ; fhMCEtaMassPtTrue[index] = new TH2F(Form("hMCEtaMassPtTrue_pt%d_cell%d_asym%d",ipt,icell,iasym), Form("Reconstructed #it{M} vs generated #it{p}_{T} of true #eta cluster pairs for #it{p}_{T} >%2.2f, #it{N}^{cluster}_{cell}>%d and |#it{A}|>%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]), nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ; fhMCEtaMassPtTrue[index]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})"); fhMCEtaMassPtTrue[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{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 #it{p}_{T} of true #eta cluster pairs, 0.01 < rec. mass < 0.17 MeV/#it{c}^{2} for #it{p}_{T} >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]), nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ; fhMCEtaPtTruePtRec[index]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})"); fhMCEtaPtTruePtRec[index]->SetYTitle("#it{p}_{T, reconstructed} (GeV/#it{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("#it{p}_{T, generated} (GeV/#it{c})"); fhMCPi0MassPtTrue[0]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{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/#it{c}^{2}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ; fhMCPi0PtTruePtRec[0]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})"); fhMCPi0PtTruePtRec[0]->SetYTitle("#it{p}_{T, reconstructed} (GeV/#it{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("#it{p}_{T, generated} (GeV/#it{c})"); fhMCEtaMassPtTrue[0]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{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/#it{c}^{2}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ; fhMCEtaPtTruePtRec[0]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})"); fhMCEtaPtTruePtRec[0]->SetYTitle("#it{p}_{T, reconstructed} (GeV/#it{c})"); outputContainer->Add(fhMCEtaPtTruePtRec[0]) ; } } } 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; imodSetXTitle("#it{p}_{T} (GeV/#it{c})"); fhReMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})"); outputContainer->Add(fhReMod[imod]) ; if(GetCalorimeter()==kPHOS) { 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) ; fhReDiffPHOSMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhReDiffPHOSMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})"); outputContainer->Add(fhReDiffPHOSMod[imod]) ; } else{//EMCAL if(imodSetXTitle("#it{p}_{T} (GeV/#it{c})"); fhReSameSectorEMCALMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})"); outputContainer->Add(fhReSameSectorEMCALMod[imod]) ; } if(imodSetXTitle("#it{p}_{T} (GeV/#it{c})"); fhReSameSideEMCALMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})"); outputContainer->Add(fhReSameSideEMCALMod[imod]) ; } }//EMCAL if(DoOwnMix()) { snprintf(key, buffersize,"hMiMod_%d",imod) ; snprintf(title, buffersize,"Mixed #it{M}_{#gamma#gamma} distr. for Module %d",imod) ; fhMiMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ; fhMiMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhMiMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})"); outputContainer->Add(fhMiMod[imod]) ; if(GetCalorimeter()==kPHOS){ snprintf(key, buffersize,"hMiDiffPHOSMod_%d",imod) ; snprintf(title, buffersize,"Mixed pairs PHOS, clusters in different Modules: %s",(pairnamePHOS[imod]).Data()) ; fhMiDiffPHOSMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ; fhMiDiffPHOSMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhMiDiffPHOSMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})"); outputContainer->Add(fhMiDiffPHOSMod[imod]) ; }//PHOS else{//EMCAL if(imodSetXTitle("#it{p}_{T} (GeV/#it{c})"); fhMiSameSectorEMCALMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})"); outputContainer->Add(fhMiSameSectorEMCALMod[imod]) ; } if(imodSetXTitle("#it{p}_{T} (GeV/#it{c})"); fhMiSameSideEMCALMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})"); outputContainer->Add(fhMiSameSideEMCALMod[imod]) ; } }//EMCAL }// own mix }//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("#it{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("#it{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("**** Print %s %s ****\n", GetName(), GetTitle() ) ; AliAnaCaloTrackCorrBaseClass::Print(" "); 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",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? printf("Number of modules: %d \n",fNModules) ; printf("Select pairs with their angle: %d, edep %d, min angle %2.3f, max angle %2.3f \n",fUseAngleCut, fUseAngleEDepCut, fAngleCut, fAngleMaxCut) ; printf("Asymmetry cuts: n = %d, \n",fNAsymCuts) ; printf("\tasymmetry < "); for(Int_t i = 0; i < fNAsymCuts; i++) printf("%2.2f ",fAsymCuts[i]); printf("\n"); printf("PID selection bits: n = %d, \n",fNPIDBits) ; printf("\tPID bit = "); for(Int_t i = 0; i < fNPIDBits; i++) printf("%d ",fPIDBits[i]); printf("\n"); if(fMultiCutAna){ printf("pT cuts: n = %d, \n",fNPtCuts) ; printf("\tpT > "); for(Int_t i = 0; i < fNPtCuts; i++) printf("%2.2f ",fPtCuts[i]); printf("GeV/c\n"); printf("N cell in cluster cuts: n = %d, \n",fNCellNCuts) ; printf("\tnCell > "); for(Int_t i = 0; i < fNCellNCuts; i++) printf("%d ",fCellNCuts[i]); printf("\n"); } printf("------------------------------------------------------\n") ; } //________________________________________ void AliAnaPi0::FillAcceptanceHistograms() { //Fill acceptance histograms if MC data is available Double_t mesonY = -100 ; Double_t mesonE = -1 ; Double_t mesonPt = -1 ; Double_t mesonPhi = 100 ; Double_t mesonYeta= -1 ; Int_t pdg = 0 ; Int_t nprim = 0 ; Int_t nDaught = 0 ; Int_t iphot1 = 0 ; Int_t iphot2 = 0 ; Float_t cen = GetEventCentrality(); Float_t ep = GetEventPlaneAngle(); TParticle * primStack = 0; AliAODMCParticle * primAOD = 0; // Get the ESD MC particles container AliStack * stack = 0; if( GetReader()->ReadStack() ) { stack = GetMCStack(); if(!stack ) return; nprim = stack->GetNtrack(); } // Get the AOD MC particles container TClonesArray * mcparticles = 0; if( GetReader()->ReadAODMCParticles() ) { mcparticles = GetReader()->GetAODMCParticles(); if( !mcparticles ) return; nprim = mcparticles->GetEntriesFast(); } for(Int_t i=0 ; i < nprim; i++) { if(GetReader()->AcceptOnlyHIJINGLabels() && !GetReader()->IsHIJINGLabel(i)) continue ; if(GetReader()->ReadStack()) { primStack = stack->Particle(i) ; if(!primStack) { printf("AliAnaPi0::FillAcceptanceHistograms() - ESD primaries pointer not available!!\n"); continue; } // If too small skip if( primStack->Energy() < 0.4 ) continue; pdg = primStack->GetPdgCode(); nDaught = primStack->GetNDaughters(); iphot1 = primStack->GetDaughter(0) ; iphot2 = primStack->GetDaughter(1) ; if(primStack->Energy() == TMath::Abs(primStack->Pz())) continue ; //Protection against floating point exception //printf("i %d, %s %d %s %d \n",i, stack->Particle(i)->GetName(), stack->Particle(i)->GetPdgCode(), // prim->GetName(), prim->GetPdgCode()); //Photon kinematics primStack->Momentum(fPi0Mom); mesonY = 0.5*TMath::Log((primStack->Energy()+primStack->Pz())/(primStack->Energy()-primStack->Pz())) ; } else { primAOD = (AliAODMCParticle *) mcparticles->At(i); if(!primAOD) { printf("AliAnaPi0::FillAcceptanceHistograms() - AOD primaries pointer not available!!\n"); continue; } // If too small skip if( primAOD->E() < 0.4 ) continue; pdg = primAOD->GetPdgCode(); nDaught = primAOD->GetNDaughters(); iphot1 = primAOD->GetFirstDaughter() ; iphot2 = primAOD->GetLastDaughter() ; if(primAOD->E() == TMath::Abs(primAOD->Pz())) continue ; //Protection against floating point exception //Photon kinematics fPi0Mom.SetPxPyPzE(primAOD->Px(),primAOD->Py(),primAOD->Pz(),primAOD->E()); mesonY = 0.5*TMath::Log((primAOD->E()+primAOD->Pz())/(primAOD->E()-primAOD->Pz())) ; } // Select only pi0 or eta if( pdg != 111 && pdg != 221) continue ; mesonPt = fPi0Mom.Pt () ; mesonE = fPi0Mom.E () ; mesonYeta= fPi0Mom.Eta() ; mesonPhi = fPi0Mom.Phi() ; if( mesonPhi < 0 ) mesonPhi+=TMath::TwoPi(); mesonPhi *= TMath::RadToDeg(); if(pdg == 111) { if(TMath::Abs(mesonY) < 1.0) { fhPrimPi0E ->Fill(mesonE ) ; fhPrimPi0Pt ->Fill(mesonPt) ; fhPrimPi0Phi->Fill(mesonPt, mesonPhi) ; fhPrimPi0YetaYcut ->Fill(mesonPt,mesonYeta) ; fhPrimPi0PtCentrality->Fill(mesonPt,cen) ; fhPrimPi0PtEventPlane->Fill(mesonPt,ep ) ; } fhPrimPi0Y ->Fill(mesonPt, mesonY) ; fhPrimPi0Yeta->Fill(mesonPt, mesonYeta) ; } else if(pdg == 221) { if(TMath::Abs(mesonY) < 1.0) { fhPrimEtaE ->Fill(mesonE ) ; fhPrimEtaPt ->Fill(mesonPt) ; fhPrimEtaPhi->Fill(mesonPt, mesonPhi) ; fhPrimEtaYetaYcut ->Fill(mesonPt,mesonYeta) ; fhPrimEtaPtCentrality->Fill(mesonPt,cen) ; fhPrimEtaPtEventPlane->Fill(mesonPt,ep ) ; } fhPrimEtaY ->Fill(mesonPt, mesonY) ; fhPrimEtaYeta->Fill(mesonPt, mesonYeta) ; } //Origin of meson if(fFillOriginHisto && TMath::Abs(mesonY) < 0.7) { Int_t momindex = -1; Int_t mompdg = -1; Int_t momstatus = -1; Float_t momR = 0; if(GetReader()->ReadStack()) momindex = primStack->GetFirstMother(); if(GetReader()->ReadAODMCParticles()) momindex = primAOD ->GetMother(); if(momindex >= 0 && momindex < nprim) { if(GetReader()->ReadStack()) { TParticle* mother = stack->Particle(momindex); mompdg = TMath::Abs(mother->GetPdgCode()); momstatus = mother->GetStatusCode(); momR = mother->R(); } if(GetReader()->ReadAODMCParticles()) { AliAODMCParticle* mother = (AliAODMCParticle*) mcparticles->At(momindex); mompdg = TMath::Abs(mother->GetPdgCode()); momstatus = mother->GetStatus(); momR = TMath::Sqrt(mother->Xv()*mother->Xv()+mother->Yv()*mother->Yv()); } if(pdg == 111) { if (momstatus == 21)fhPrimPi0PtOrigin->Fill(mesonPt,0.5);//parton else if(mompdg < 22 ) fhPrimPi0PtOrigin->Fill(mesonPt,1.5);//quark else if(mompdg > 2100 && mompdg < 2210) fhPrimPi0PtOrigin->Fill(mesonPt,2.5);// resonances else if(mompdg == 221) fhPrimPi0PtOrigin->Fill(mesonPt,8.5);//eta else if(mompdg == 331) fhPrimPi0PtOrigin->Fill(mesonPt,9.5);//eta prime else if(mompdg == 213) fhPrimPi0PtOrigin->Fill(mesonPt,4.5);//rho else if(mompdg == 223) fhPrimPi0PtOrigin->Fill(mesonPt,5.5);//omega else if(mompdg >= 310 && mompdg <= 323) fhPrimPi0PtOrigin->Fill(mesonPt,6.5);//k0S, k+-,k* else if(mompdg == 130) fhPrimPi0PtOrigin->Fill(mesonPt,6.5);//k0L else if(momstatus == 11 || momstatus == 12 ) fhPrimPi0PtOrigin->Fill(mesonPt,3.5);//resonances else fhPrimPi0PtOrigin->Fill(mesonPt,7.5);//other? //printf("Prim Pi0: index %d, pt %2.2f Prod vertex %3.3f, origin pdg %d, origin status %d, origin UI %d\n", // momindex, mesonPt,mother->R(),mompdg,momstatus,mother->GetUniqueID()); fhPrimPi0ProdVertex->Fill(mesonPt,momR); }//pi0 else { if (momstatus == 21 ) fhPrimEtaPtOrigin->Fill(mesonPt,0.5);//parton else if(mompdg < 22 ) fhPrimEtaPtOrigin->Fill(mesonPt,1.5);//quark else if(mompdg > 2100 && mompdg < 2210) fhPrimEtaPtOrigin->Fill(mesonPt,2.5);//qq resonances else if(mompdg == 331) fhPrimEtaPtOrigin->Fill(mesonPt,5.5);//eta prime else if(momstatus == 11 || momstatus == 12 ) fhPrimEtaPtOrigin->Fill(mesonPt,3.5);//resonances else fhPrimEtaPtOrigin->Fill(mesonPt,4.5);//stable, conversions? //printf("Other Meson pdg %d, Mother %s, pdg %d, status %d\n",pdg, TDatabasePDG::Instance()->GetParticle(mompdg)->GetName(),mompdg, momstatus ); fhPrimEtaProdVertex->Fill(mesonPt,momR); } } // pi0 has mother } //Check if both photons hit Calorimeter if(nDaught != 2 ) continue; //Only interested in 2 gamma decay if(iphot1 < 0 || iphot1 >= nprim || iphot2 < 0 || iphot2 >= nprim) continue ; Int_t pdg1 = 0; Int_t pdg2 = 0; Bool_t inacceptance1 = kTRUE; Bool_t inacceptance2 = kTRUE; if(GetReader()->ReadStack()) { TParticle * phot1 = stack->Particle(iphot1) ; TParticle * phot2 = stack->Particle(iphot2) ; if(!phot1 || !phot2) continue ; pdg1 = phot1->GetPdgCode(); pdg2 = phot2->GetPdgCode(); phot1->Momentum(fPhotonMom1); phot2->Momentum(fPhotonMom2); // Check if photons hit the Calorimeter acceptance if(IsRealCaloAcceptanceOn()) { if( !GetCaloUtils()->IsMCParticleInCalorimeterAcceptance( GetCalorimeter(), phot1 )) inacceptance1 = kFALSE ; if( !GetCaloUtils()->IsMCParticleInCalorimeterAcceptance( GetCalorimeter(), phot2 )) inacceptance2 = kFALSE ; } } if(GetReader()->ReadAODMCParticles()) { AliAODMCParticle * phot1 = (AliAODMCParticle *) mcparticles->At(iphot1) ; AliAODMCParticle * phot2 = (AliAODMCParticle *) mcparticles->At(iphot2) ; if(!phot1 || !phot2) continue ; pdg1 = phot1->GetPdgCode(); pdg2 = phot2->GetPdgCode(); fPhotonMom1.SetPxPyPzE(phot1->Px(),phot1->Py(),phot1->Pz(),phot1->E()); fPhotonMom2.SetPxPyPzE(phot2->Px(),phot2->Py(),phot2->Pz(),phot2->E()); // Check if photons hit the Calorimeter acceptance if(IsRealCaloAcceptanceOn()) { if( !GetCaloUtils()->IsMCParticleInCalorimeterAcceptance( GetCalorimeter(), phot1 )) inacceptance1 = kFALSE ; if( !GetCaloUtils()->IsMCParticleInCalorimeterAcceptance( GetCalorimeter(), phot2 )) inacceptance2 = kFALSE ; } } if( pdg1 != 22 || pdg2 !=22) continue ; // Check if photons hit desired acceptance in the fidutial borders fixed in the analysis if(IsFiducialCutOn()) { if( inacceptance1 && !GetFiducialCut()->IsInFiducialCut(fPhotonMom1.Eta(), fPhotonMom1.Phi(), GetCalorimeter()) ) inacceptance1 = kFALSE ; if( inacceptance2 && !GetFiducialCut()->IsInFiducialCut(fPhotonMom2.Eta(), fPhotonMom2.Phi(), GetCalorimeter()) ) inacceptance2 = kFALSE ; } if(fFillArmenterosThetaStar) FillArmenterosThetaStar(pdg); if(GetCalorimeter()==kEMCAL && inacceptance1 && inacceptance2 && fCheckAccInSector) { Int_t absID1=0; Int_t absID2=0; Float_t photonPhi1 = fPhotonMom1.Phi(); Float_t photonPhi2 = fPhotonMom2.Phi(); if(photonPhi1 < 0) photonPhi1+=TMath::TwoPi(); if(photonPhi2 < 0) photonPhi2+=TMath::TwoPi(); GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(fPhotonMom1.Eta(),photonPhi1,absID1); GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(fPhotonMom2.Eta(),photonPhi2,absID2); Int_t sm1 = GetEMCALGeometry()->GetSuperModuleNumber(absID1); Int_t sm2 = GetEMCALGeometry()->GetSuperModuleNumber(absID2); Int_t j=0; Bool_t sameSector = kFALSE; for(Int_t isector = 0; isector < fNModules/2; isector++) { j=2*isector; if((sm1==j && sm2==j+1) || (sm1==j+1 && sm2==j)) sameSector = kTRUE; } if(sm1!=sm2 && !sameSector) { inacceptance1 = kFALSE; inacceptance2 = kFALSE; } //if(sm1!=sm2)printf("sm1 %d, sm2 %d, same sector %d, in acceptance %d\n",sm1,sm2,sameSector,inacceptance); // if(GetEMCALGeometry()->Impact(phot1) && GetEMCALGeometry()->Impact(phot2)) // inacceptance = kTRUE; } if(GetDebug() > 2) printf("Accepted in %s?: m (%2.2f,%2.2f,%2.2f), p1 (%2.2f,%2.2f,%2.2f), p2 (%2.2f,%2.2f,%2.2f) : in1 %d, in2 %d\n", GetCalorimeterString().Data(), mesonPt,mesonYeta,mesonPhi, fPhotonMom1.Pt(),fPhotonMom1.Eta(),fPhotonMom1.Phi()*TMath::RadToDeg(), fPhotonMom2.Pt(),fPhotonMom2.Eta(),fPhotonMom2.Phi()*TMath::RadToDeg(), inacceptance1, inacceptance2); if(inacceptance1 && inacceptance2) { Float_t asym = TMath::Abs((fPhotonMom1.E()-fPhotonMom2.E()) / (fPhotonMom1.E()+fPhotonMom2.E())); Double_t angle = fPhotonMom1.Angle(fPhotonMom2.Vect()); if(GetDebug() > 2) printf("\t ACCEPTED pdg %d: pt %2.2f, phi %2.2f, eta %2.2f\n",pdg,mesonPt,mesonPhi,mesonYeta); if(pdg==111) { fhPrimPi0AccE ->Fill(mesonE) ; fhPrimPi0AccPt ->Fill(mesonPt) ; fhPrimPi0AccPhi ->Fill(mesonPt, mesonPhi) ; fhPrimPi0AccY ->Fill(mesonPt, mesonY) ; fhPrimPi0AccYeta->Fill(mesonPt, mesonYeta) ; fhPrimPi0AccPtCentrality->Fill(mesonPt,cen) ; fhPrimPi0AccPtEventPlane->Fill(mesonPt,ep ) ; if(fFillAngleHisto) { fhPrimPi0OpeningAngle ->Fill(mesonPt,angle); if(mesonPt > 5)fhPrimPi0OpeningAngleAsym->Fill(asym,angle); fhPrimPi0CosOpeningAngle ->Fill(mesonPt,TMath::Cos(angle)); } } else if(pdg==221) { fhPrimEtaAccE ->Fill(mesonE ) ; fhPrimEtaAccPt ->Fill(mesonPt) ; fhPrimEtaAccPhi ->Fill(mesonPt, mesonPhi) ; fhPrimEtaAccY ->Fill(mesonPt, mesonY) ; fhPrimEtaAccYeta->Fill(mesonPt, mesonYeta) ; fhPrimEtaAccPtCentrality->Fill(mesonPt,cen) ; fhPrimEtaAccPtEventPlane->Fill(mesonPt,ep ) ; if(fFillAngleHisto) { fhPrimEtaOpeningAngle ->Fill(mesonPt,angle); if(mesonPt > 5)fhPrimEtaOpeningAngleAsym->Fill(asym,angle); fhPrimEtaCosOpeningAngle ->Fill(mesonPt,TMath::Cos(angle)); } } }//Accepted }//loop on primaries } //________________________________________________ void AliAnaPi0::FillArmenterosThetaStar(Int_t pdg) { // Fill armenteros plots // Get pTArm and AlphaArm Float_t momentumSquaredMother = fPi0Mom.P()*fPi0Mom.P(); Float_t momentumDaughter1AlongMother = 0.; Float_t momentumDaughter2AlongMother = 0.; if (momentumSquaredMother > 0.) { momentumDaughter1AlongMother = (fPhotonMom1.Px()*fPi0Mom.Px() + fPhotonMom1.Py()*fPi0Mom.Py()+ fPhotonMom1.Pz()*fPi0Mom.Pz()) / sqrt(momentumSquaredMother); momentumDaughter2AlongMother = (fPhotonMom2.Px()*fPi0Mom.Px() + fPhotonMom2.Py()*fPi0Mom.Py()+ fPhotonMom2.Pz()*fPi0Mom.Pz()) / sqrt(momentumSquaredMother); } Float_t momentumSquaredDaughter1 = fPhotonMom1.P()*fPhotonMom1.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); fPhotonMom1Boost = fPhotonMom1; fPhotonMom1Boost.Boost(-fPi0Mom.BoostVector()); Float_t cosThStar=TMath::Cos(fPhotonMom1Boost.Vect().Angle(fPi0Mom.Vect())); Float_t en = fPi0Mom.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(fPhotonMom1.E()+fPhotonMom2.E() > 0) asym = TMath::Abs(fPhotonMom1.E()-fPhotonMom2.E())/(fPhotonMom1.E()+fPhotonMom2.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 Int_t ancPDG = 0; Int_t ancStatus = 0; Int_t ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(index1, index2, GetReader(), ancPDG, ancStatus,fPi0Mom, fProdVertex); Int_t momindex = -1; Int_t mompdg = -1; Int_t momstatus = -1; if(GetDebug() > 1 ) { if(ancLabel >= 0) printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Common ancestor label %d, pdg %d, name %s, status %d; \n", ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus); else printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Common ancestor not found \n"); } Float_t prodR = -1; Int_t mcIndex = -1; if(ancLabel > -1) { if(ancPDG==22) {//gamma mcIndex = 0; } else if(TMath::Abs(ancPDG)==11) {//e mcIndex = 1; } else if(ancPDG==111) {//Pi0 mcIndex = 2; if(fMultiCutAnaSim) { for(Int_t ipt=0; ipt fPtCuts[ipt] && pt2 > fPtCuts[ipt] && asym < fAsymCuts[iasym] && ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell]) { fhMCPi0MassPtRec [index]->Fill(pt,mass); fhMCPi0MassPtTrue[index]->Fill(fPi0Mom.Pt(),mass); if(mass < 0.17 && mass > 0.1) fhMCPi0PtTruePtRec[index]->Fill(fPi0Mom.Pt(),pt); }//pass the different cuts }// pid bit cut loop }// icell loop }// pt cut loop }//Multi cut ana sim else { fhMCPi0MassPtTrue[0]->Fill(fPi0Mom.Pt(),mass); if(mass < 0.17 && mass > 0.1) { fhMCPi0PtTruePtRec[0]->Fill(fPi0Mom.Pt(),pt); //Int_t uniqueId = -1; 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(); prodR = mother->R(); //uniqueId = mother->GetUniqueID(); } 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(); prodR = TMath::Sqrt(mother->Xv()*mother->Xv()+mother->Yv()*mother->Yv()); //uniqueId = mother->GetUniqueID(); } // printf("Reco Pi0: pt %2.2f Prod vertex %3.3f, origin pdg %d, origin status %d, origin UI %d\n", // pt,prodR,mompdg,momstatus,uniqueId); fhMCPi0ProdVertex->Fill(pt,prodR); 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 } } else if(ancPDG==221) {//Eta mcIndex = 3; if(fMultiCutAnaSim) { for(Int_t ipt=0; ipt fPtCuts[ipt] && pt2 > fPtCuts[ipt] && asym < fAsymCuts[iasym] && ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell]) { fhMCEtaMassPtRec [index]->Fill(pt,mass); fhMCEtaMassPtTrue[index]->Fill(fPi0Mom.Pt(),mass); if(mass < 0.65 && mass > 0.45) fhMCEtaPtTruePtRec[index]->Fill(fPi0Mom.Pt(),pt); }//pass the different cuts }// pid bit cut loop }// icell loop }// pt cut loop } //Multi cut ana sim else { fhMCEtaMassPtTrue[0]->Fill(fPi0Mom.Pt(),mass); if(mass < 0.65 && mass > 0.45) fhMCEtaPtTruePtRec[0]->Fill(fPi0Mom.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(); prodR = mother->R(); } 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(); prodR = TMath::Sqrt(mother->Xv()*mother->Xv()+mother->Yv()*mother->Yv()); } fhMCEtaProdVertex->Fill(pt,prodR); 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 mcIndex = 4; } else if(ancPDG==-2112){//ANeutron mcIndex = 5; } else if(TMath::Abs(ancPDG)==13){//muons mcIndex = 6; } else if (TMath::Abs(ancPDG) > 100 && ancLabel > 7) { if(ancStatus==1) {//Stable particles, converted? not decayed resonances mcIndex = 6; } else {//resonances and other decays, more hadron conversions? mcIndex = 7; } } else {//Partons, colliding protons, strings, intermediate corrections if(ancStatus==11 || ancStatus==12) {//String fragmentation mcIndex = 8; } else if (ancStatus==21){ if(ancLabel < 2) {//Colliding protons mcIndex = 11; }//colliding protons else if(ancLabel < 6) {//partonic initial states interactions mcIndex = 9; } else if(ancLabel < 8) {//Final state partons radiations? mcIndex = 10; } // 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); // } }////Partons, colliding protons, strings, intermediate corrections }//ancLabel > -1 else { //ancLabel <= -1 //printf("Not related at all label = %d\n",ancLabel); mcIndex = 12; } if(mcIndex >=0 && mcIndex < 13) { fhMCOrgMass[mcIndex]->Fill(pt,mass); fhMCOrgAsym[mcIndex]->Fill(pt,asym); fhMCOrgDeltaEta[mcIndex]->Fill(pt,deta); fhMCOrgDeltaPhi[mcIndex]->Fill(pt,dphi); } } //__________________________________________ void AliAnaPi0::MakeAnalysisFillHistograms() { //Process one event and extract photons from AOD branch // filled with AliAnaPhoton and fill histos with invariant mass //In case of simulated data, fill acceptance histograms if(IsDataMC())FillAcceptanceHistograms(); //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() ; 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 ) { 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 ; } 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 = 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 (GetCalorimeter()==kEMCAL) clusters = GetEMCALClusters(); else if(GetCalorimeter()==kPHOS ) clusters = GetPHOSClusters() ; //--------------------------------- //First loop on photons/clusters //--------------------------------- for(Int_t i1=0; i1At(i1)) ; //printf("AliAnaPi0::MakeAnalysisFillHistograms() : cluster1 id %d/%d\n",i1,nPhot-1); // get the event index in the mixed buffer where the photon comes from // in case of mixing with analysis frame, not own mixing evtIndex1 = GetEventIndex(p1, vert) ; if ( evtIndex1 == -1 ) return ; if ( evtIndex1 == -2 ) continue ; // Only effective in case of mixed event frame if(TMath::Abs(vert[2]) > GetZvertexCut()) continue ; //vertex cut if (evtIndex1 != currentEvtIndex) { //Fill event bin info 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 ; } //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 1 Evt %d Vertex : %f,%f,%f\n",evtIndex1, GetVertex(evtIndex1)[0] ,GetVertex(evtIndex1)[1],GetVertex(evtIndex1)[2]); //Get the momentum of this cluster fPhotonMom1.SetPxPyPzE(p1->Px(),p1->Py(),p1->Pz(),p1->E()); //Get (Super)Module number of this cluster module1 = GetModuleNumber(p1); //------------------------------------------ // Recover original cluster Int_t iclus1 = -1 ; AliVCluster * cluster1 = FindCluster(clusters,p1->GetCaloLabel(0),iclus1); if(!cluster1) printf("AliAnaPi0 - Cluster1 not found!\n"); //--------------------------------- //Second loop on photons/clusters //--------------------------------- for(Int_t i2=i1+1; i2At(i2)) ; //printf("AliAnaPi0::MakeAnalysisFillHistograms() : cluster2 i %d/%d\n",i2,nPhot); //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 ) return ; if ( evtIndex2 == -2 ) continue ; if (GetMixedEvent() && (evtIndex1 == evtIndex2)) continue ; //------------------------------------------ // Recover original cluster Int_t iclus2 = -1; AliVCluster * cluster2 = FindCluster(clusters,p2->GetCaloLabel(0),iclus2,iclus1+1); // start new loop from iclus1+1 to gain some time if(!cluster2) printf("AliAnaPi0 - Cluster2 not found!\n"); // Get the TOF,l0 and ncells from the clusters Float_t tof1 = -1; Float_t l01 = -1; Int_t ncell1 = 0; if(cluster1) { tof1 = cluster1->GetTOF()*1e9; l01 = cluster1->GetM02(); ncell1 = cluster1->GetNCells(); //printf("cluster1: E %2.2f (%2.2f), l0 %2.2f, tof %2.2f\n",cluster1->E(),p1->E(),l01,tof1); } //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; Float_t l02 = -1; Int_t ncell2 = 0; if(cluster2) { tof2 = cluster2->GetTOF()*1e9; l02 = cluster2->GetM02(); ncell2 = cluster2->GetNCells(); //printf("cluster2: E %2.2f (%2.2f), l0 %2.2f, tof %2.2f\n",cluster2->E(),p2->E(),l02,tof2); } //else printf("cluster2 not available: calo label %d / %d, cluster ID %d\n", // p2->GetCaloLabel(0),(GetReader()->GetInputEvent())->GetNumberOfCaloClusters()-1,cluster2->GetID()); if(cluster1 && cluster2) { 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 fPhotonMom2.SetPxPyPzE(p2->Px(),p2->Py(),p2->Pz(),p2->E()); //Get module number module2 = GetModuleNumber(p2); //--------------------------------- // Get pair kinematics //--------------------------------- Double_t m = (fPhotonMom1 + fPhotonMom2).M() ; Double_t pt = (fPhotonMom1 + fPhotonMom2).Pt(); Double_t deta = fPhotonMom1.Eta() - fPhotonMom2.Eta(); Double_t dphi = fPhotonMom1.Phi() - fPhotonMom2.Phi(); Double_t a = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ; if(GetDebug() > 2) printf(" E: fPhotonMom1 %f, fPhotonMom2 %f; Pair: pT %f, mass %f, a %f\n", p1->E(), p2->E(), (fPhotonMom1 + fPhotonMom2).E(),m,a); //-------------------------------- // Opening angle selection //-------------------------------- //Check if opening angle is too large or too small compared to what is expected Double_t angle = fPhotonMom1.Angle(fPhotonMom2.Vect()); if(fUseAngleEDepCut && !GetNeutralMesonSelection()->IsAngleInWindow((fPhotonMom1+fPhotonMom2).E(),angle+0.05)) { if(GetDebug() > 2) printf("AliAnaPi0::MakeAnalysisFillHistograms() -Real pair angle %f not in E %f window\n",angle, (fPhotonMom1+fPhotonMom2).E()); continue; } if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut)) { if(GetDebug() > 2) printf("AliAnaPi0::MakeAnalysisFillHistograms() - Real pair cut %f < angle %f < cut %f\n",fAngleCut, angle, fAngleMaxCut); continue; } //------------------------------------------------------------------------------------------------- //Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array //------------------------------------------------------------------------------------------------- if(a < fAsymCuts[0] && fFillSMCombinations) { if(module1==module2 && module1 >=0 && module1Fill(pt,m) ; if(GetCalorimeter()==kEMCAL) { // Same sector Int_t j=0; 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) ; } // Same side for(Int_t i = 0; i < fNModules-2; i++){ if((module1==i && module2==i+2) || (module1==i+2 && module2==i)) fhReSameSideEMCALMod[i]->Fill(pt,m); } }//EMCAL else {//PHOS if((module1==0 && module2==1) || (module1==1 && module2==0)) fhReDiffPHOSMod[0]->Fill(pt,m) ; if((module1==0 && module2==2) || (module1==2 && module2==0)) fhReDiffPHOSMod[1]->Fill(pt,m) ; if((module1==1 && module2==2) || (module1==2 && module2==1)) fhReDiffPHOSMod[2]->Fill(pt,m) ; }//PHOS } //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) { //Check if one of the clusters comes from a conversion if(fCheckConversion) { if (p1->IsTagged() && p2->IsTagged()) fhReConv2->Fill(pt,m); else if(p1->IsTagged() || p2->IsTagged()) fhReConv ->Fill(pt,m); } // Fill shower shape cut histograms if(fFillSSCombinations) { if ( l01 > 0.01 && l01 < 0.4 && l02 > 0.01 && l02 < 0.4 ) fhReSS[0]->Fill(pt,m); // Tight else if( l01 > 0.4 && l02 > 0.4 ) fhReSS[1]->Fill(pt,m); // Loose else if( l01 > 0.01 && l01 < 0.4 && l02 > 0.4 ) fhReSS[2]->Fill(pt,m); // Both else if( l02 > 0.01 && l02 < 0.4 && l01 > 0.4 ) fhReSS[2]->Fill(pt,m); // Both } //Fill histograms for different bad channel distance, centrality, assymmetry cut and pid bit for(Int_t ipid=0; ipidIsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)) && (p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton))) { for(Int_t iasym=0; iasym < fNAsymCuts; 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 - 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) { if(p1->DistToBad()>0 && p2->DistToBad()>0) { fhRe2 [index]->Fill(pt,m) ; if(fMakeInvPtPlots)fhReInvPt2[index]->Fill(pt,m,1./pt) ; if(p1->DistToBad()>1 && p2->DistToBad()>1) { fhRe3 [index]->Fill(pt,m) ; if(fMakeInvPtPlots)fhReInvPt3[index]->Fill(pt,m,1./pt) ; }// bad 3 }// bad2 }// Fill bad dist histos }//assymetry cut }// asymmetry cut loop }// bad 1 }// pid bit loop //Fill histograms with opening angle if(fFillAngleHisto) { fhRealOpeningAngle ->Fill(pt,angle); fhRealCosOpeningAngle->Fill(pt,TMath::Cos(angle)); } //Fill histograms with pair assymmetry 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); } //--------- // MC data //--------- //Do some MC checks on the origin of the pair, is there any common ancestor and if there is one, who? if(IsDataMC()) { if(GetMCAnalysisUtils()->CheckTagBit(p1->GetTag(),AliMCAnalysisUtils::kMCConversion) && GetMCAnalysisUtils()->CheckTagBit(p2->GetTag(),AliMCAnalysisUtils::kMCConversion)) { fhReMCFromConversion->Fill(pt,m); } else if(!GetMCAnalysisUtils()->CheckTagBit(p1->GetTag(),AliMCAnalysisUtils::kMCConversion) && !GetMCAnalysisUtils()->CheckTagBit(p2->GetTag(),AliMCAnalysisUtils::kMCConversion)) { fhReMCFromNotConversion->Fill(pt,m); } else { fhReMCFromMixConversion->Fill(pt,m); } if(fFillOriginHisto) FillMCVersusRecDataHistograms(p1->GetLabel(), p2->GetLabel(),p1->Pt(), p2->Pt(),ncell1, ncell2, m, pt, a,deta, dphi); } //----------------------- //Multi cuts analysis //----------------------- if(fMultiCutAna) { //Histograms for different PID bits selection for(Int_t ipid=0; ipidIsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton) && p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)) fhRePIDBits[ipid]->Fill(pt,m) ; //printf("ipt %d, ipid%d, name %s\n",ipt, ipid, fhRePtPIDCuts[ipt*fNPIDBitsBits+ipid]->GetName()); } // pid bit cut loop //Several pt,ncell and asymmetry cuts for(Int_t ipt=0; iptE() > fPtCuts[ipt] && p2->E() > fPtCuts[ipt] && a < fAsymCuts[iasym] && 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(fFillSMCombinations && module1==module2) { fhRePtNCellAsymCutsSM[module1][index]->Fill(pt,m) ; } } }// pid bit cut loop }// icell loop }// pt cut loop if(GetHistogramRanges()->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 }// second same event particle }// first cluster //------------------------------------------------------------- // Mixing //------------------------------------------------------------- if(DoOwnMix()) { //Recover events in with same characteristics as the current event //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; iiAt(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, GetEventCentralityBin()); fhEventMixBin->Fill(eventbin) ; //--------------------------------- //First loop on photons/clusters //--------------------------------- for(Int_t i1=0; i1At(i1)) ; if(fSameSM && GetModuleNumber(p1)!=module1) continue; //Get kinematics of cluster and (super) module of this cluster fPhotonMom1.SetPxPyPzE(p1->Px(),p1->Py(),p1->Pz(),p1->E()); module1 = GetModuleNumber(p1); //--------------------------------- //First loop on photons/clusters //--------------------------------- for(Int_t i2=0; i2At(i2)) ; //Get kinematics of second cluster and calculate those of the pair fPhotonMom2.SetPxPyPzE(p2->Px(),p2->Py(),p2->Pz(),p2->E()); m = (fPhotonMom1+fPhotonMom2).M() ; Double_t pt = (fPhotonMom1 + fPhotonMom2).Pt(); Double_t a = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ; //Check if opening angle is too large or too small compared to what is expected Double_t angle = fPhotonMom1.Angle(fPhotonMom2.Vect()); if(fUseAngleEDepCut && !GetNeutralMesonSelection()->IsAngleInWindow((fPhotonMom1+fPhotonMom2).E(),angle+0.05)) { if(GetDebug() > 2) printf("AliAnaPi0::MakeAnalysisFillHistograms() -Mix pair angle %f not in E %f window\n",angle, (fPhotonMom1+fPhotonMom2).E()); continue; } if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut)) { if(GetDebug() > 2) printf("AliAnaPi0::MakeAnalysisFillHistograms() -Mix pair angle %f < cut %f\n",angle,fAngleCut); continue; } if(GetDebug() > 2) printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed Event: pT: fPhotonMom1 %2.2f, fPhotonMom2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n", p1->Pt(), p2->Pt(), pt,m,a); //In case we want only pairs in same (super) module, check their origin. module2 = GetModuleNumber(p2); //------------------------------------------------------------------------------------------------- //Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array //------------------------------------------------------------------------------------------------- if(a < fAsymCuts[0] && fFillSMCombinations){ if(module1==module2 && module1 >=0 && module1Fill(pt,m) ; if(GetCalorimeter()==kEMCAL) { // Same sector Int_t j=0; for(Int_t i = 0; i < fNModules/2; i++) { j=2*i; if((module1==j && module2==j+1) || (module1==j+1 && module2==j)) fhMiSameSectorEMCALMod[i]->Fill(pt,m) ; } // Same side for(Int_t i = 0; i < fNModules-2; i++) { if((module1==i && module2==i+2) || (module1==i+2 && module2==i)) fhMiSameSideEMCALMod[i]->Fill(pt,m); } }//EMCAL else {//PHOS if((module1==0 && module2==1) || (module1==1 && module2==0)) fhMiDiffPHOSMod[0]->Fill(pt,m) ; if((module1==0 && module2==2) || (module1==2 && module2==0)) fhMiDiffPHOSMod[1]->Fill(pt,m) ; if((module1==1 && module2==2) || (module1==2 && module2==1)) fhMiDiffPHOSMod[2]->Fill(pt,m) ; }//PHOS } Bool_t ok = kTRUE; if(fSameSM && module1!=module2) ok=kFALSE; if(ok){ //Check if one of the clusters comes from a conversion 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; ipidIsPIDOK(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) { fhMi2 [index]->Fill(pt,m) ; if(fMakeInvPtPlots)fhMiInvPt2[index]->Fill(pt,m,1./pt) ; if(p1->DistToBad()>1 && p2->DistToBad()>1) { fhMi3 [index]->Fill(pt,m) ; if(fMakeInvPtPlots)fhMiInvPt3[index]->Fill(pt,m,1./pt) ; } } }// Fill bad dist histo }//Asymmetry cut }// Asymmetry loop }//PID cut }//loop for histograms //----------------------- //Multi cuts analysis //----------------------- if(fMultiCutAna){ //Several pt,ncell and asymmetry cuts for(Int_t ipt=0; iptPt() > fPtCuts[ipt] && p2->Pt() > fPtCuts[ipt] && a < fAsymCuts[iasym] //&& //p1->GetBtag() >= fCellNCuts[icell] && p2->GetBtag() >= fCellNCuts[icell] // trick, correct it. ) { fhMiPtNCellAsymCuts[index]->Fill(pt,m) ; //printf("ipt %d, icell%d, iasym %d, name %s\n",ipt, icell, iasym, fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->GetName()); } }// pid bit cut loop }// icell loop }// pt cut loop } // Multi cut ana //Fill histograms with opening angle if(fFillAngleHisto) { fhMixedOpeningAngle ->Fill(pt,angle); fhMixedCosOpeningAngle->Fill(pt,TMath::Cos(angle)); } }//ok }// second cluster loop }//first cluster loop }//loop on mixed events //-------------------------------------------------------- // Add the current event to the list of events for mixing //-------------------------------------------------------- TClonesArray *currentEvent = new TClonesArray(*GetInputAODBranch()); //Add current event to buffer and Remove redundant events 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() >= GetNMaxEvMix() ) { TClonesArray * tmp = (TClonesArray*) (evMixList->Last()) ; evMixList->RemoveLast() ; delete tmp ; } } else { //empty event delete currentEvent ; currentEvent=0 ; } }// DoOwnMix if(GetDebug() > 0) printf("AliAnaPi0::MakeAnalysisFillHistograms() - End fill histograms\n"); } //________________________________________________________________________ Int_t AliAnaPi0::GetEventIndex(AliAODPWG4Particle * part, Double_t * vert) { // retieves the event index and checks the vertex // in the mixed buffer returns -2 if vertex NOK // for normal events returns 0 if vertex OK and -1 if vertex NOK Int_t evtIndex = -1 ; if(GetReader()->GetDataType()!=AliCaloTrackReader::kMC) { if (GetMixedEvent()) { evtIndex = GetMixedEvent()->EventIndexForCaloCluster(part->GetCaloLabel(0)) ; GetVertex(vert,evtIndex); if(TMath::Abs(vert[2])> GetZvertexCut()) evtIndex = -2 ; //Event can not be used (vertex, centrality,... cuts not fulfilled) } else { // Single event GetVertex(vert); if(TMath::Abs(vert[2])> GetZvertexCut()) evtIndex = -1 ; //Event can not be used (vertex, centrality,... cuts not fulfilled) else evtIndex = 0 ; } }//No MC reader else { evtIndex = 0; vert[0] = 0. ; vert[1] = 0. ; vert[2] = 0. ; } return evtIndex ; }