1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
16 //_________________________________________________________________________
17 // Class for analysis of particle isolation
18 // Input is selected particles put in AOD branch (AliAODPWG4ParticleCorrelation)
20 // Class created from old AliPHOSGammaJet
21 // (see AliRoot versions previous Release 4-09)
23 // -- Author: Gustavo Conesa (LNF-INFN)
25 //-Yaxian Mao (add the possibility for different IC method with different pt range, 01/10/2010)
26 //////////////////////////////////////////////////////////////////////////////
28 // --- ROOT system ---
29 #include <TClonesArray.h>
31 #include <TObjString.h>
35 #include "TParticle.h"
36 #include "TDatabasePDG.h"
38 // --- Analysis system ---
39 #include "AliAnaParticleIsolation.h"
40 #include "AliCaloTrackReader.h"
42 #include "AliIsolationCut.h"
43 #include "AliFiducialCut.h"
44 #include "AliMCAnalysisUtils.h"
45 #include "AliNeutralMesonSelection.h"
46 #include "AliAODMCParticle.h"
47 #include "AliAODPWG4ParticleCorrelation.h"
48 #include "AliMCAnalysisUtils.h"
49 #include "AliVTrack.h"
50 #include "AliVCluster.h"
51 #include "AliESDEvent.h"
52 #include "AliAODEvent.h"
54 #include "AliEMCALGeometry.h"
55 #include "AliPHOSGeoUtils.h"
57 ClassImp(AliAnaParticleIsolation)
59 //______________________________________________________________________________
60 AliAnaParticleIsolation::AliAnaParticleIsolation() :
61 AliAnaCaloTrackCorrBaseClass(),
62 fCalorimeter(""), fIsoDetector(""),
63 fReMakeIC(0), fMakeSeveralIC(0),
64 fFillPileUpHistograms(0),
65 fFillTMHisto(0), fFillSSHisto(1),
66 fFillUEBandSubtractHistograms(1), fFillCellHistograms(0),
67 fFillHighMultHistograms(0), fFillTaggedDecayHistograms(0),
68 fNDecayBits(0), fDecayBits(),
69 fFillNLMHistograms(0),
70 fLeadingOnly(0), fCheckLeadingWithNeutralClusters(0),
71 fSelectPrimariesInCone(0), fMakePrimaryPi0DecayStudy(0),
72 fFillBackgroundBinHistograms(0), fNBkgBin(0),
74 fNCones(0), fNPtThresFrac(0),
75 fConeSizes(), fPtThresholds(),
76 fPtFractions(), fSumPtThresholds(),
78 fhEIso(0), fhPtIso(0),
79 fhPtCentralityIso(0), fhPtEventPlaneIso(0),
81 fhPhiIso(0), fhEtaIso(0), fhEtaPhiIso(0),
83 fhENoIso(0), fhPtNoIso(0), fhPtNLocMaxNoIso(0),
85 fhPtClusterInCone(0), fhPtCellInCone(0), fhPtTrackInCone(0),
86 fhPtTrackInConeOtherBC(0), fhPtTrackInConeOtherBCPileUpSPD(0),
87 fhPtTrackInConeBC0(0), fhPtTrackInConeVtxBC0(0),
88 fhPtTrackInConeBC0PileUpSPD(0),
89 fhPtInConePileUp(), fhPtInConeCent(0),
90 fhPerpConeSumPt(0), fhPtInPerpCone(0),
91 fhEtaPhiInConeCluster(0), fhEtaPhiCluster(0),
92 fhEtaPhiInConeTrack(0), fhEtaPhiTrack(0),
93 fhEtaBandCluster(0), fhPhiBandCluster(0),
94 fhEtaBandTrack(0), fhPhiBandTrack(0),
95 fhEtaBandCell(0), fhPhiBandCell(0),
96 fhConePtLead(0), fhConePtLeadCluster(0), fhConePtLeadTrack(0),
97 fhConePtLeadClustervsTrack(0), fhConePtLeadClusterTrackFrac(0),
98 fhConeSumPt(0), fhConeSumPtCellTrack(0),
99 fhConeSumPtCell(0), fhConeSumPtCluster(0), fhConeSumPtTrack(0),
100 fhConeSumPtEtaBandUECluster(0), fhConeSumPtPhiBandUECluster(0),
101 fhConeSumPtEtaBandUETrack(0), fhConeSumPtPhiBandUETrack(0),
102 fhConeSumPtEtaBandUECell(0), fhConeSumPtPhiBandUECell(0),
103 fhConeSumPtTrigEtaPhi(0),
104 fhConeSumPtCellTrackTrigEtaPhi(0),
105 fhConeSumPtEtaBandUEClusterTrigEtaPhi(0), fhConeSumPtPhiBandUEClusterTrigEtaPhi(0),
106 fhConeSumPtEtaBandUETrackTrigEtaPhi(0), fhConeSumPtPhiBandUETrackTrigEtaPhi(0),
107 fhConeSumPtEtaBandUECellTrigEtaPhi(0), fhConeSumPtPhiBandUECellTrigEtaPhi(0),
108 fhConeSumPtEtaUESub(0), fhConeSumPtPhiUESub(0),
109 fhConeSumPtEtaUESubTrigEtaPhi(0), fhConeSumPtPhiUESubTrigEtaPhi(0),
110 fhConeSumPtEtaUESubTrackCell(0), fhConeSumPtPhiUESubTrackCell(0),
111 fhConeSumPtEtaUESubTrackCellTrigEtaPhi(0), fhConeSumPtPhiUESubTrackCellTrigEtaPhi(0),
112 fhConeSumPtEtaUESubCluster(0), fhConeSumPtPhiUESubCluster(0),
113 fhConeSumPtEtaUESubClusterTrigEtaPhi(0), fhConeSumPtPhiUESubClusterTrigEtaPhi(0),
114 fhConeSumPtEtaUESubCell(0), fhConeSumPtPhiUESubCell(0),
115 fhConeSumPtEtaUESubCellTrigEtaPhi(0), fhConeSumPtPhiUESubCellTrigEtaPhi(0),
116 fhConeSumPtEtaUESubTrack(0), fhConeSumPtPhiUESubTrack(0),
117 fhConeSumPtEtaUESubTrackTrigEtaPhi(0), fhConeSumPtPhiUESubTrackTrigEtaPhi(0),
118 fhFractionTrackOutConeEta(0), fhFractionTrackOutConeEtaTrigEtaPhi(0),
119 fhFractionClusterOutConeEta(0), fhFractionClusterOutConeEtaTrigEtaPhi(0),
120 fhFractionClusterOutConePhi(0), fhFractionClusterOutConePhiTrigEtaPhi(0),
121 fhFractionCellOutConeEta(0), fhFractionCellOutConeEtaTrigEtaPhi(0),
122 fhFractionCellOutConePhi(0), fhFractionCellOutConePhiTrigEtaPhi(0),
123 fhConeSumPtClustervsTrack(0), fhConeSumPtClusterTrackFrac(0),
124 fhConeSumPtEtaUESubClustervsTrack(0), fhConeSumPtPhiUESubClustervsTrack(0),
125 fhConeSumPtCellvsTrack(0),
126 fhConeSumPtEtaUESubCellvsTrack(0), fhConeSumPtPhiUESubCellvsTrack(0),
127 fhEtaBandClustervsTrack(0), fhPhiBandClustervsTrack(0),
128 fhEtaBandNormClustervsTrack(0), fhPhiBandNormClustervsTrack(0),
129 fhEtaBandCellvsTrack(0), fhPhiBandCellvsTrack(0),
130 fhEtaBandNormCellvsTrack(0), fhPhiBandNormCellvsTrack(0),
131 fhConeSumPtSubvsConeSumPtTotPhiTrack(0), fhConeSumPtSubNormvsConeSumPtTotPhiTrack(0),
132 fhConeSumPtSubvsConeSumPtTotEtaTrack(0), fhConeSumPtSubNormvsConeSumPtTotEtaTrack(0),
133 fhConeSumPtSubvsConeSumPtTotPhiCluster(0), fhConeSumPtSubNormvsConeSumPtTotPhiCluster(0),
134 fhConeSumPtSubvsConeSumPtTotEtaCluster(0), fhConeSumPtSubNormvsConeSumPtTotEtaCluster(0),
135 fhConeSumPtSubvsConeSumPtTotPhiCell(0), fhConeSumPtSubNormvsConeSumPtTotPhiCell(0),
136 fhConeSumPtSubvsConeSumPtTotEtaCell(0), fhConeSumPtSubNormvsConeSumPtTotEtaCell(0),
137 fhConeSumPtVSUETracksEtaBand(0), fhConeSumPtVSUETracksPhiBand(0),
138 fhConeSumPtVSUEClusterEtaBand(0), fhConeSumPtVSUEClusterPhiBand(0),
139 fhPtPrimMCPi0DecayPairOutOfCone(0), fhPtPrimMCPi0DecayPairOutOfAcceptance(0),
140 fhPtPrimMCPi0DecayPairAcceptInConeLowPt(0), fhPtPrimMCPi0DecayPairAcceptInConeLowPtNoOverlap(0), fhPtPrimMCPi0DecayPairAcceptInConeLowPtNoOverlapCaloE(0),
141 fhPtPrimMCPi0DecayPairNoOverlap(0),
142 fhPtPrimMCPi0DecayIsoPairOutOfCone(0), fhPtPrimMCPi0DecayIsoPairOutOfAcceptance(0),
143 fhPtPrimMCPi0DecayIsoPairAcceptInConeLowPt(0),fhPtPrimMCPi0DecayIsoPairAcceptInConeLowPtNoOverlap(0), fhPtPrimMCPi0DecayIsoPairAcceptInConeLowPtNoOverlapCaloE(0),
144 fhPtPrimMCPi0DecayIsoPairNoOverlap(0),
145 fhPtPrimMCPi0Overlap(0), fhPtPrimMCPi0IsoOverlap(0),
146 fhPtLeadConeBinLambda0(0), fhSumPtConeBinLambda0(0),
147 fhPtLeadConeBinLambda0MC(0), fhSumPtConeBinLambda0MC(0),
148 // Number of local maxima in cluster
150 fhELambda0LocMax1(), fhELambda1LocMax1(),
151 fhELambda0LocMax2(), fhELambda1LocMax2(),
152 fhELambda0LocMaxN(), fhELambda1LocMaxN(),
154 fhEIsoPileUp(), fhPtIsoPileUp(),
155 fhENoIsoPileUp(), fhPtNoIsoPileUp(),
156 fhTimeENoCut(0), fhTimeESPD(0), fhTimeESPDMulti(0),
157 fhTimeNPileUpVertSPD(0), fhTimeNPileUpVertTrack(0),
158 fhTimeNPileUpVertContributors(0),
159 fhTimePileUpMainVertexZDistance(0), fhTimePileUpMainVertexZDiamond(0)
163 //Initialize parameters
166 for(Int_t i = 0; i < 5 ; i++)
170 for(Int_t imc = 0; imc < fgkNmcTypes; imc++)
171 fhSumPtLeadingPtMC[imc][i] = 0 ;
173 for(Int_t j = 0; j < 5 ; j++)
175 fhPtThresIsolated [i][j] = 0 ;
176 fhPtFracIsolated [i][j] = 0 ;
177 fhSumPtIsolated [i][j] = 0 ;
179 fhEtaPhiPtThresIso [i][j] = 0 ;
180 fhEtaPhiPtThresDecayIso [i][j] = 0 ;
181 fhPtPtThresDecayIso [i][j] = 0 ;
183 fhEtaPhiPtFracIso [i][j] = 0 ;
184 fhEtaPhiPtFracDecayIso [i][j] = 0 ;
185 fhPtPtFracDecayIso [i][j] = 0 ;
186 fhPtPtSumDecayIso [i][j] = 0 ;
187 fhPtSumDensityIso [i][j] = 0 ;
188 fhPtSumDensityDecayIso [i][j] = 0 ;
189 fhEtaPhiSumDensityIso [i][j] = 0 ;
190 fhEtaPhiSumDensityDecayIso [i][j] = 0 ;
191 fhPtFracPtSumIso [i][j] = 0 ;
192 fhPtFracPtSumDecayIso [i][j] = 0 ;
193 fhEtaPhiFracPtSumIso [i][j] = 0 ;
194 fhEtaPhiFracPtSumDecayIso [i][j] = 0 ;
196 for(Int_t imc = 0; imc < fgkNmcTypes; imc++)
198 fhPtThresIsolatedMC[imc][i][j] = 0 ;
199 fhPtFracIsolatedMC [imc][i][j] = 0 ;
200 fhSumPtIsolatedMC [imc][i][j] = 0 ;
206 for(Int_t ibit =0; ibit< 4; ibit++)
208 fhPtDecayIso [ibit] = 0;
209 fhPtLambda0Decay[0][ibit] = 0;
210 fhPtLambda0Decay[1][ibit] = 0;
211 fhPtDecayNoIso [ibit] = 0;
212 fhEtaPhiDecayIso [ibit] = 0;
213 fhEtaPhiDecayNoIso [ibit] = 0;
214 for(Int_t imc = 0; imc < fgkNmcTypes; imc++)
216 fhPtDecayIsoMC [ibit][imc] = 0;
217 fhPtDecayNoIsoMC[ibit][imc] = 0;
221 for(Int_t i = 0; i < 5 ; i++)
223 fPtFractions [i] = 0 ;
224 fPtThresholds [i] = 0 ;
225 fSumPtThresholds[i] = 0 ;
227 fhSumPtLeadingPt [i] = 0 ;
228 fhPtLeadingPt [i] = 0 ;
229 fhPerpSumPtLeadingPt[i] = 0 ;
230 fhPerpPtLeadingPt [i] = 0 ;
233 for(Int_t imc = 0; imc < fgkNmcTypes; imc++)
235 fhPtNoIsoMC [imc] = 0;
237 fhPhiIsoMC [imc] = 0;
238 fhEtaIsoMC [imc] = 0;
239 fhPtLambda0MC[imc][0] = 0;
240 fhPtLambda0MC[imc][1] = 0;
243 for(Int_t i = 0; i < 2 ; i++)
245 fhTrackMatchedDEta[i] = 0 ; fhTrackMatchedDPhi[i] = 0 ; fhTrackMatchedDEtaDPhi [i] = 0 ;
246 fhdEdx [i] = 0 ; fhEOverP [i] = 0 ; fhTrackMatchedMCParticle[i] = 0 ;
247 fhELambda0 [i] = 0 ; fhELambda1 [i] = 0 ; fhPtLambda0 [i] = 0 ;
248 fhELambda0TRD [i] = 0 ; fhELambda1TRD [i] = 0 ; fhPtLambda0TRD [i] = 0 ;
250 // Number of local maxima in cluster
252 fhELambda0LocMax1[i] = 0 ; fhELambda1LocMax1[i] = 0 ;
253 fhELambda0LocMax2[i] = 0 ; fhELambda1LocMax2[i] = 0 ;
254 fhELambda0LocMaxN[i] = 0 ; fhELambda1LocMaxN[i] = 0 ;
258 for(Int_t i = 0; i < fgkNmcPrimTypes; i++)
260 fhPtPrimMCiso[i] = 0;
269 for(Int_t i = 0 ; i < 7 ; i++)
271 fhPtInConePileUp[i] = 0 ;
272 fhEIsoPileUp [i] = 0 ;
273 fhPtIsoPileUp [i] = 0 ;
274 fhENoIsoPileUp [i] = 0 ;
275 fhPtNoIsoPileUp [i] = 0 ;
279 //_______________________________________________________________________________________________
280 void AliAnaParticleIsolation::CalculateCaloUEBand(AliAODPWG4ParticleCorrelation * pCandidate,
281 Float_t & etaBandPtSum, Float_t & phiBandPtSum)
283 // Get the clusters pT or sum of pT in phi/eta bands or at 45 degrees from trigger
285 if( GetIsolationCut()->GetParticleTypeInCone()==AliIsolationCut::kOnlyCharged ) return ;
287 Float_t conesize = GetIsolationCut()->GetConeSize();
290 //Select the Calorimeter
291 TObjArray * pl = 0x0;
292 if (fCalorimeter == "PHOS" )
293 pl = GetPHOSClusters();
294 else if (fCalorimeter == "EMCAL")
295 pl = GetEMCALClusters();
299 //Get vertex for cluster momentum calculation
300 Double_t vertex[] = {0,0,0} ; //vertex ;
301 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
302 GetReader()->GetVertex(vertex);
304 Float_t ptTrig = pCandidate->Pt() ;
305 Float_t phiTrig = pCandidate->Phi();
306 Float_t etaTrig = pCandidate->Eta();
308 for(Int_t icluster=0; icluster < pl->GetEntriesFast(); icluster++)
310 AliVCluster* cluster = (AliVCluster *) pl->At(icluster);
314 printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Cluster not available?");
318 //Do not count the candidate (photon or pi0) or the daughters of the candidate
319 if(cluster->GetID() == pCandidate->GetCaloLabel(0) ||
320 cluster->GetID() == pCandidate->GetCaloLabel(1) ) continue ;
322 //Remove matched clusters to tracks if Neutral and Track info is used
323 if( GetIsolationCut()->GetParticleTypeInCone()==AliIsolationCut::kNeutralAndCharged &&
324 IsTrackMatched(cluster,GetReader()->GetInputEvent())) continue ;
326 cluster->GetMomentum(mom,vertex) ;//Assume that come from vertex in straight line
328 //exclude particles in cone
329 Float_t rad = GetIsolationCut()->Radius(etaTrig, phiTrig, mom.Eta(), mom.Phi());
331 // histo of eta and phi for all clusters
332 fhEtaPhiCluster->Fill(mom.Eta(), mom.Phi());
334 // histos for all clusters in cone
335 fhEtaPhiInConeCluster->Fill(mom.Eta(), mom.Phi());
338 //fill histogram for UE in phi band in EMCal acceptance
339 if(mom.Eta() > (etaTrig-conesize) && mom.Eta() < (etaTrig+conesize))
341 phiBandPtSum+=mom.Pt();
342 fhPhiBandCluster->Fill(mom.Eta(),mom.Phi());
346 //fill histogram for UE in eta band in EMCal acceptance
347 if(mom.Phi() > (phiTrig-conesize) && mom.Phi() < (phiTrig+conesize))
349 etaBandPtSum+=mom.Pt();
350 fhEtaBandCluster->Fill(mom.Eta(),mom.Phi());
354 fhConeSumPtEtaBandUECluster ->Fill(ptTrig , etaBandPtSum);
355 fhConeSumPtPhiBandUECluster ->Fill(ptTrig , phiBandPtSum);
356 fhConeSumPtEtaBandUEClusterTrigEtaPhi->Fill(etaTrig,phiTrig,etaBandPtSum);
357 fhConeSumPtPhiBandUEClusterTrigEtaPhi->Fill(etaTrig,phiTrig,phiBandPtSum);
361 //________________________________________________________________________________________________
362 void AliAnaParticleIsolation::CalculateCaloCellUEBand(AliAODPWG4ParticleCorrelation * pCandidate,
363 Float_t & etaBandPtSumCells, Float_t & phiBandPtSumCells)
365 // Get the cells amplitude or sum of amplitude in phi/eta bands or at 45 degrees from trigger
367 if( GetIsolationCut()->GetParticleTypeInCone()==AliIsolationCut::kOnlyCharged ) return ;
369 Float_t conesize = GetIsolationCut()->GetConeSize();
371 Float_t phiTrig = pCandidate->Phi();
372 if(phiTrig<0) phiTrig += TMath::TwoPi();
373 Float_t etaTrig = pCandidate->Eta();
375 if(pCandidate->GetDetector()=="EMCAL")
377 AliEMCALGeometry* eGeom = AliEMCALGeometry::GetInstance();
380 if (eGeom->GetAbsCellIdFromEtaPhi(etaTrig,phiTrig,absId))
382 if(!eGeom->CheckAbsCellId(absId)) return ;
384 // Get absolute (col,row) of trigger particle
385 Int_t nSupMod = eGeom->GetSuperModuleNumber(absId);
387 Int_t imEta=-1, imPhi=-1;
388 Int_t ieta =-1, iphi =-1;
390 if (eGeom->GetCellIndex(absId,nSupMod,nModule,imPhi,imEta))
392 eGeom->GetCellPhiEtaIndexInSModule(nSupMod,nModule,imPhi,imEta,iphi,ieta);
394 Int_t colTrig = ieta;
395 if (nSupMod % 2) colTrig = AliEMCALGeoParams::fgkEMCALCols + ieta ;
396 Int_t rowTrig = iphi + AliEMCALGeoParams::fgkEMCALRows*int(nSupMod/2);
398 Int_t sqrSize = int(conesize/0.0143);
400 AliVCaloCells * cells = GetEMCALCells();
402 Int_t nTotalRows = AliEMCALGeoParams::fgkEMCALRows*16/3 ; // 24*(16/3) 5 full-size Sectors (2 SM) + 1 one-third Sector (2 SM)
403 Int_t nTotalCols = 2*AliEMCALGeoParams::fgkEMCALCols;
404 // printf("nTotalRows %i, nTotalCols %i\n",nTotalRows,nTotalCols);
405 // Loop on cells in eta band
407 Int_t irowmin = rowTrig-sqrSize;
408 if(irowmin<0) irowmin=0;
409 Int_t irowmax = rowTrig+sqrSize;
410 if(irowmax>AliEMCALGeoParams::fgkEMCALRows) irowmax=AliEMCALGeoParams::fgkEMCALRows;
413 for(Int_t irow = irowmin; irow <irowmax; irow++)
415 for(Int_t icol = 0; icol < nTotalCols; icol++)
417 Int_t inSector = int(irow/AliEMCALGeoParams::fgkEMCALRows);
418 if(inSector==5) continue;
421 if(icol < AliEMCALGeoParams::fgkEMCALCols)
423 inSupMod = 2*inSector + 1;
426 else if(icol > AliEMCALGeoParams::fgkEMCALCols - 1)
428 inSupMod = 2*inSector;
429 icolLoc = icol-AliEMCALGeoParams::fgkEMCALCols;
432 Int_t irowLoc = irow - AliEMCALGeoParams::fgkEMCALRows*inSector ;
434 // Exclude cells in cone
435 if(TMath::Abs(icol-colTrig) < sqrSize || TMath::Abs(irow-rowTrig) < sqrSize){
438 Int_t iabsId = eGeom->GetAbsCellIdFromCellIndexes(inSupMod,irowLoc,icolLoc);
439 if(!eGeom->CheckAbsCellId(iabsId)) continue;
440 etaBandPtSumCells += cells->GetCellAmplitude(iabsId);
441 fhEtaBandCell->Fill(colTrig,rowTrig);
443 // printf("ETA inSupMod %i,irowLoc %i,icolLoc %i, iabsId %i, etaBandPtSumCells %f\n",nSupMod,irowLoc,icolLoc,iabsId,etaBandPtSumCells);
446 Int_t icolmin = colTrig-sqrSize;
447 if(icolmin<0) icolmin=0;
448 Int_t icolmax = colTrig+sqrSize;
449 if(icolmax>AliEMCALGeoParams::fgkEMCALCols) icolmax=AliEMCALGeoParams::fgkEMCALCols;
451 // Loop on cells in phi band
452 for(Int_t icol = icolmin; icol < icolmax; icol++)
454 for(Int_t irow = 0; irow < nTotalRows; irow++)
456 Int_t inSector = int(irow/AliEMCALGeoParams::fgkEMCALRows);
457 if(inSector==5) continue;
460 // printf("icol %i, irow %i, inSector %i\n",icol,irow ,inSector);
461 if(icol < AliEMCALGeoParams::fgkEMCALCols)
463 // printf("icol < AliEMCALGeoParams::fgkEMCALCols %i\n",AliEMCALGeoParams::fgkEMCALCols );
464 inSupMod = 2*inSector + 1;
467 else if(icol > AliEMCALGeoParams::fgkEMCALCols - 1)
469 // printf("icol > AliEMCALGeoParams::fgkEMCALCols -1 %i\n",AliEMCALGeoParams::fgkEMCALCols -1 );
470 inSupMod = 2*inSector;
471 icolLoc = icol-AliEMCALGeoParams::fgkEMCALCols;
474 Int_t irowLoc = irow - AliEMCALGeoParams::fgkEMCALRows*inSector ; // Stesso problema di sopra //
476 // Exclude cells in cone
477 if(TMath::Abs(icol-colTrig) < sqrSize) {
478 //printf("TMath::Abs(icol-colTrig) %i < sqrSize %i\n",TMath::Abs(icol-colTrig) ,sqrSize);continue ;
480 if(TMath::Abs(irow-rowTrig) < sqrSize) {
481 //printf("TMath::Abs(irow-rowTrig) %i < sqrSize %i\n",TMath::Abs(irow-rowTrig) ,sqrSize);continue ;
484 Int_t iabsId = eGeom->GetAbsCellIdFromCellIndexes(inSupMod,irowLoc,icolLoc);
485 if(!eGeom->CheckAbsCellId(iabsId)) {printf("!eGeom->CheckAbsCellId(iabsId=%i) inSupMod %i irowLoc %i icolLoc %i \n",iabsId,inSupMod, irowLoc, icolLoc);continue;}
486 phiBandPtSumCells += cells->GetCellAmplitude(iabsId);
487 fhPhiBandCell->Fill(colTrig,rowTrig);
488 //printf("inSupMod %i,irowLoc %i,icolLoc %i, iabsId %i, phiBandPtSumCells %f\n",nSupMod,irowLoc,icolLoc,iabsId,phiBandPtSumCells);
495 Float_t ptTrig = pCandidate->Pt();
497 fhConeSumPtEtaBandUECell ->Fill(ptTrig , etaBandPtSumCells);
498 fhConeSumPtPhiBandUECell ->Fill(ptTrig , phiBandPtSumCells);
499 fhConeSumPtEtaBandUECellTrigEtaPhi->Fill(etaTrig,phiTrig,etaBandPtSumCells);
500 fhConeSumPtPhiBandUECellTrigEtaPhi->Fill(etaTrig,phiTrig,phiBandPtSumCells);
504 //________________________________________________________________________________________________
505 void AliAnaParticleIsolation::CalculateTrackUEBand(AliAODPWG4ParticleCorrelation * pCandidate,
506 Float_t & etaBandPtSum, Float_t & phiBandPtSum)
508 // Get the track pT or sum of pT in phi/eta bands or at 45 degrees from trigger
510 if( GetIsolationCut()->GetParticleTypeInCone()==AliIsolationCut::kOnlyNeutral ) return ;
512 Float_t conesize = GetIsolationCut()->GetConeSize();
514 Double_t sumptPerp= 0. ;
515 Float_t ptTrig = pCandidate->Pt() ;
516 Float_t phiTrig = pCandidate->Phi();
517 Float_t etaTrig = pCandidate->Eta();
519 TObjArray * trackList = GetCTSTracks() ;
520 for(Int_t itrack=0; itrack < trackList->GetEntriesFast(); itrack++)
522 AliVTrack* track = (AliVTrack *) trackList->At(itrack);
526 printf("AliAnaParticleIsolation::CalculateTrackUEBand() - Track not available?");
530 //Do not count the candidate (pion, conversion photon) or the daughters of the candidate
531 if(track->GetID() == pCandidate->GetTrackLabel(0) || track->GetID() == pCandidate->GetTrackLabel(1) ||
532 track->GetID() == pCandidate->GetTrackLabel(2) || track->GetID() == pCandidate->GetTrackLabel(3) ) continue ;
534 // histo of eta:phi for all tracks
535 fhEtaPhiTrack->Fill(track->Eta(),track->Phi());
537 //exclude particles in cone
538 Float_t rad = GetIsolationCut()->Radius(etaTrig, phiTrig, track->Eta(), track->Phi());
540 // histo of eta:phi for all tracks in cone
541 fhEtaPhiInConeTrack->Fill(track->Eta(),track->Phi());
545 //fill histogram for UE in phi band
546 if(track->Eta() > (etaTrig-conesize) && track->Eta() < (etaTrig+conesize))
548 phiBandPtSum+=track->Pt();
549 fhPhiBandTrack->Fill(track->Eta(),track->Phi());
552 //fill histogram for UE in eta band in EMCal acceptance
553 if(track->Phi() > (phiTrig-conesize) && track->Phi() < (phiTrig+conesize))
555 etaBandPtSum+=track->Pt();
556 fhEtaBandTrack->Fill(track->Eta(),track->Phi());
559 //fill the histograms at +-45 degrees in phi from trigger particle, perpedicular to trigger axis in phi
560 Double_t dPhi = phiTrig - track->Phi() + TMath::PiOver2();
561 Double_t dEta = etaTrig - track->Eta();
562 Double_t arg = dPhi*dPhi + dEta*dEta;
563 if(TMath::Sqrt(arg) < conesize)
565 fhPtInPerpCone->Fill(ptTrig,TMath::Sqrt(track->Px()*track->Px()+track->Py()*track->Py()));
566 sumptPerp+=track->Pt();
569 dPhi = phiTrig - track->Phi() - TMath::PiOver2();
570 arg = dPhi*dPhi + dEta*dEta;
571 if(TMath::Sqrt(arg) < conesize)
573 fhPtInPerpCone->Fill(ptTrig,TMath::Sqrt(track->Px()*track->Px()+track->Py()*track->Py()));
574 sumptPerp+=track->Pt();
578 fhPerpConeSumPt ->Fill(ptTrig , sumptPerp );
579 fhConeSumPtEtaBandUETrack ->Fill(ptTrig , etaBandPtSum);
580 fhConeSumPtPhiBandUETrack ->Fill(ptTrig , phiBandPtSum);
581 fhConeSumPtEtaBandUETrackTrigEtaPhi->Fill(etaTrig,phiTrig,etaBandPtSum);
582 fhConeSumPtPhiBandUETrackTrigEtaPhi->Fill(etaTrig,phiTrig,phiBandPtSum);
588 //_____________________________________________________________________________________________________________________________________
589 void AliAnaParticleIsolation::CalculateNormalizeUEBandPerUnitArea(AliAODPWG4ParticleCorrelation * pCandidate, Float_t coneptsumCluster,
590 Float_t coneptsumCell, Float_t coneptsumTrack,
591 Float_t &etaBandptsumTrackNorm, Float_t &etaBandptsumClusterNorm)
593 //normalize phi/eta band per area unit
595 Float_t etaUEptsumTrack = 0 ;
596 Float_t phiUEptsumTrack = 0 ;
597 Float_t etaUEptsumCluster = 0 ;
598 Float_t phiUEptsumCluster = 0 ;
599 Float_t etaUEptsumCell = 0 ;
600 Float_t phiUEptsumCell = 0 ;
602 Int_t partTypeInCone = GetIsolationCut()->GetParticleTypeInCone();
604 // Do the normalization
606 Float_t conesize = GetIsolationCut()->GetConeSize();
607 Float_t coneA = conesize*conesize*TMath::Pi(); // A = pi R^2, isolation cone area
608 Float_t ptTrig = pCandidate->Pt() ;
609 Float_t phiTrig = pCandidate->Phi();
610 Float_t etaTrig = pCandidate->Eta();
616 Float_t phiUEptsumTrackNorm = 0 ;
617 Float_t etaUEptsumTrackNorm = 0 ;
618 Float_t coneptsumTrackSubPhi = 0 ;
619 Float_t coneptsumTrackSubEta = 0 ;
620 Float_t coneptsumTrackSubPhiNorm = 0 ;
621 Float_t coneptsumTrackSubEtaNorm = 0 ;
622 etaBandptsumTrackNorm = 0 ;
624 if( partTypeInCone!=AliIsolationCut::kOnlyNeutral )
626 // Sum the pT in the phi or eta band for clusters or tracks
627 CalculateTrackUEBand (pCandidate,etaUEptsumTrack ,phiUEptsumTrack );// rajouter ici l'histo eta phi
630 fhConeSumPtVSUETracksEtaBand->Fill(coneptsumTrack,etaUEptsumTrack);
631 fhConeSumPtVSUETracksPhiBand->Fill(coneptsumTrack,phiUEptsumTrack);
634 Float_t correctConeSumTrack = 1;
635 Float_t correctConeSumTrackPhi = 1;
637 GetIsolationCut()->CalculateUEBandTrackNormalization(GetReader(),etaTrig, phiTrig,
638 phiUEptsumTrack,etaUEptsumTrack,
639 phiUEptsumTrackNorm,etaUEptsumTrackNorm,
640 correctConeSumTrack,correctConeSumTrackPhi);
642 coneptsumTrackSubPhi = coneptsumTrack - phiUEptsumTrackNorm;
643 coneptsumTrackSubEta = coneptsumTrack - etaUEptsumTrackNorm;
645 etaBandptsumTrackNorm = etaUEptsumTrackNorm;
647 fhConeSumPtPhiUESubTrack ->Fill(ptTrig , coneptsumTrackSubPhi);
648 fhConeSumPtPhiUESubTrackTrigEtaPhi ->Fill(etaTrig, phiTrig, coneptsumTrackSubPhi);
649 fhConeSumPtEtaUESubTrack ->Fill(ptTrig , coneptsumTrackSubEta);
650 fhConeSumPtEtaUESubTrackTrigEtaPhi ->Fill(etaTrig, phiTrig, coneptsumTrackSubEta);
652 fhFractionTrackOutConeEta ->Fill(ptTrig , correctConeSumTrack-1);
653 fhFractionTrackOutConeEtaTrigEtaPhi->Fill(etaTrig, phiTrig,correctConeSumTrack-1);
655 if(coneptsumTrack > 0)
657 coneptsumTrackSubPhiNorm = coneptsumTrackSubPhi/coneptsumTrack;
658 coneptsumTrackSubEtaNorm = coneptsumTrackSubEta/coneptsumTrack;
661 fhConeSumPtSubvsConeSumPtTotPhiTrack ->Fill(coneptsumTrack,coneptsumTrackSubPhi);
662 fhConeSumPtSubNormvsConeSumPtTotPhiTrack->Fill(coneptsumTrack,coneptsumTrackSubPhiNorm);
663 fhConeSumPtSubvsConeSumPtTotEtaTrack ->Fill(coneptsumTrack,coneptsumTrackSubEta);
664 fhConeSumPtSubNormvsConeSumPtTotEtaTrack->Fill(coneptsumTrack,coneptsumTrackSubEtaNorm);
668 // ------------------------ //
669 // EMCal Clusters and cells //
670 // ------------------------ //
671 Float_t phiUEptsumClusterNorm = 0 ;
672 Float_t etaUEptsumClusterNorm = 0 ;
673 Float_t coneptsumClusterSubPhi = 0 ;
674 Float_t coneptsumClusterSubEta = 0 ;
675 Float_t coneptsumClusterSubPhiNorm = 0 ;
676 Float_t coneptsumClusterSubEtaNorm = 0 ;
677 Float_t phiUEptsumCellNorm = 0 ;
678 Float_t etaUEptsumCellNorm = 0 ;
679 Float_t coneptsumCellSubPhi = 0 ;
680 Float_t coneptsumCellSubEta = 0 ;
681 Float_t coneptsumCellSubPhiNorm = 0 ;
682 Float_t coneptsumCellSubEtaNorm = 0 ;
683 etaBandptsumClusterNorm = 0;
685 if( partTypeInCone!=AliIsolationCut::kOnlyCharged )
692 // Sum the pT in the phi or eta band for clusters or tracks
693 CalculateCaloUEBand (pCandidate,etaUEptsumCluster,phiUEptsumCluster);// rajouter ici l'histo eta phi
696 fhConeSumPtVSUEClusterEtaBand->Fill(coneptsumCluster,etaUEptsumCluster);
697 fhConeSumPtVSUEClusterPhiBand->Fill(coneptsumCluster,phiUEptsumCluster);
700 Float_t correctConeSumClusterEta = 1;
701 Float_t correctConeSumClusterPhi = 1;
703 GetIsolationCut()->CalculateUEBandClusterNormalization(GetReader(),etaTrig, phiTrig,
704 phiUEptsumCluster,etaUEptsumCluster,
705 phiUEptsumClusterNorm,etaUEptsumClusterNorm,
706 correctConeSumClusterEta,correctConeSumClusterPhi);
708 // In case that cone is out of eta and phi side, we are over correcting, not too often with the current cuts ...
709 // Comment if not used
710 // Float_t coneBadCellsCoeff =1;
711 // Float_t etaBandBadCellsCoeff=1;
712 // Float_t phiBandBadCellsCoeff=1;
713 // GetIsolationCut()->GetCoeffNormBadCell(pCandidate, GetReader(),coneBadCellsCoeff,etaBandBadCellsCoeff,phiBandBadCellsCoeff) ;
715 //coneptsumCluster=coneptsumCluster*coneBadCellsCoeff*correctConeSumClusterEta*correctConeSumClusterPhi;
717 coneptsumClusterSubPhi = coneptsumCluster - phiUEptsumClusterNorm;
718 coneptsumClusterSubEta = coneptsumCluster - etaUEptsumClusterNorm;
720 etaBandptsumClusterNorm = etaUEptsumClusterNorm;
722 fhConeSumPtPhiUESubCluster ->Fill(ptTrig , coneptsumClusterSubPhi);
723 fhConeSumPtPhiUESubClusterTrigEtaPhi ->Fill(etaTrig, phiTrig, coneptsumClusterSubPhi);
724 fhConeSumPtEtaUESubCluster ->Fill(ptTrig , coneptsumClusterSubEta);
725 fhConeSumPtEtaUESubClusterTrigEtaPhi ->Fill(etaTrig, phiTrig, coneptsumClusterSubEta);
727 fhFractionClusterOutConeEta ->Fill(ptTrig , correctConeSumClusterEta-1);
728 fhFractionClusterOutConeEtaTrigEtaPhi->Fill(etaTrig, phiTrig, correctConeSumClusterEta-1);
729 fhFractionClusterOutConePhi ->Fill(ptTrig , correctConeSumClusterPhi-1);
730 fhFractionClusterOutConePhiTrigEtaPhi->Fill(etaTrig, phiTrig, correctConeSumClusterPhi-1);
732 if(coneptsumCluster!=0)
734 coneptsumClusterSubPhiNorm = coneptsumClusterSubPhi/coneptsumCluster;
735 coneptsumClusterSubEtaNorm = coneptsumClusterSubEta/coneptsumCluster;
738 fhConeSumPtSubvsConeSumPtTotPhiCluster ->Fill(coneptsumCluster,coneptsumClusterSubPhi);
739 fhConeSumPtSubNormvsConeSumPtTotPhiCluster->Fill(coneptsumCluster,coneptsumClusterSubPhiNorm);
740 fhConeSumPtSubvsConeSumPtTotEtaCluster ->Fill(coneptsumCluster,coneptsumClusterSubEta);
741 fhConeSumPtSubNormvsConeSumPtTotEtaCluster->Fill(coneptsumCluster,coneptsumClusterSubEtaNorm);
747 if(fFillCellHistograms)
749 // Sum the pT in the phi or eta band for clusters or tracks
750 CalculateCaloCellUEBand(pCandidate,etaUEptsumCell ,phiUEptsumCell );
752 // Move to AliIsolationCut the calculation not the histograms??
754 //Careful here if EMCal limits changed .. 2010 (4 SM) to 2011-12 (10 SM), for the moment consider 100 deg in phi
755 Float_t emcEtaSize = 0.7*2; // TO FIX
756 Float_t emcPhiSize = TMath::DegToRad()*100.; // TO FIX
758 if(((2*conesize*emcPhiSize)-coneA)!=0)phiUEptsumCellNorm = phiUEptsumCell*(coneA / ((2*conesize*emcPhiSize)-coneA));
759 if(((2*conesize*emcEtaSize)-coneA)!=0)etaUEptsumCellNorm = etaUEptsumCell*(coneA / ((2*conesize*emcEtaSize)-coneA));
761 // Need to correct coneptsumCluster by the fraction of the cone out of the calorimeter cut acceptance!
763 Float_t correctConeSumCellEta = 1;
764 if(TMath::Abs(etaTrig)+conesize > emcEtaSize/2.)
766 Float_t excess = TMath::Abs(etaTrig) + conesize - emcEtaSize/2.;
767 correctConeSumCellEta = GetIsolationCut()->CalculateExcessAreaFraction(excess);
768 //printf("Excess EMC-Eta %2.3f, coneA %2.2f, excessA %2.2f, angle %2.2f,factor %2.2f\n",excess,coneA, excessA, angle*TMath::RadToDeg(), correctConeSumClusterEta);
769 // Need to correct phi band surface if part of the cone falls out of track cut acceptance!
770 if(((2*(conesize-excess)*emcPhiSize)-(coneA-correctConeSumCellEta))!=0)phiUEptsumCellNorm = phiUEptsumCell*(coneA / ((2*(conesize-excess)*emcPhiSize)-(coneA-correctConeSumCellEta)));
773 Float_t correctConeSumCellPhi = 1;
774 //printf("EMCPhiTrig %2.2f, conesize %2.2f, sum %2.2f, rest %2.2f \n",phiTrig*TMath::RadToDeg(),conesize*TMath::RadToDeg(),(phiTrig+conesize)*TMath::RadToDeg(),(phiTrig-conesize)*TMath::RadToDeg() );
775 if((phiTrig+conesize > 180*TMath::DegToRad()) ||
776 (phiTrig-conesize < 80*TMath::DegToRad()))
779 if( phiTrig+conesize > 180*TMath::DegToRad() ) excess = conesize + phiTrig - 180*TMath::DegToRad() ;
780 else excess = conesize - phiTrig + 80*TMath::DegToRad() ;
782 correctConeSumCellPhi = GetIsolationCut()->CalculateExcessAreaFraction(excess);
783 //printf("Excess EMC-Phi %2.3f, coneA %2.2f, excessA %2.2f, angle %2.2f,factor %2.2f\n",excess,coneA, excessA, angle*TMath::RadToDeg(), correctConeSumClusterPhi);
785 // Need to correct eta band surface if part of the cone falls out of track cut acceptance!
786 if(((2*(conesize-excess)*emcEtaSize)-(coneA-correctConeSumCellPhi))!=0)etaUEptsumCellNorm = etaUEptsumCell*(coneA / ((2*(conesize-excess)*emcEtaSize)-(coneA-correctConeSumCellPhi)));
790 // In case that cone is out of eta and phi side, we are over correcting, not too often with the current cuts ...
791 coneptsumCellSubPhi = coneptsumCell*correctConeSumCellEta*correctConeSumCellPhi - phiUEptsumCellNorm;
792 coneptsumCellSubEta = coneptsumCell*correctConeSumCellEta*correctConeSumCellPhi - etaUEptsumCellNorm;
794 fhConeSumPtPhiUESubCell ->Fill(ptTrig , coneptsumCellSubPhi);
795 fhConeSumPtPhiUESubCellTrigEtaPhi ->Fill(etaTrig, phiTrig, coneptsumCellSubPhi);
796 fhConeSumPtEtaUESubCell ->Fill(ptTrig , coneptsumCellSubEta);
797 fhConeSumPtEtaUESubCellTrigEtaPhi ->Fill(etaTrig, phiTrig, coneptsumCellSubEta);
799 fhFractionCellOutConeEta ->Fill(ptTrig , correctConeSumCellEta-1);
800 fhFractionCellOutConeEtaTrigEtaPhi->Fill(etaTrig, phiTrig, correctConeSumCellEta-1);
801 fhFractionCellOutConePhi ->Fill(ptTrig , correctConeSumCellPhi-1);
802 fhFractionCellOutConePhiTrigEtaPhi->Fill(etaTrig, phiTrig, correctConeSumCellPhi-1);
805 coneptsumCellSubPhiNorm = coneptsumCellSubPhi/coneptsumCell;
806 coneptsumCellSubEtaNorm = coneptsumCellSubEta/coneptsumCell;
809 fhConeSumPtSubvsConeSumPtTotPhiCell ->Fill(coneptsumCell,coneptsumCellSubPhi);
810 fhConeSumPtSubNormvsConeSumPtTotPhiCell->Fill(coneptsumCell,coneptsumCellSubPhiNorm);
811 fhConeSumPtSubvsConeSumPtTotEtaCell ->Fill(coneptsumCell,coneptsumCellSubEta);
812 fhConeSumPtSubNormvsConeSumPtTotEtaCell->Fill(coneptsumCell,coneptsumCellSubEtaNorm);
816 if( partTypeInCone==AliIsolationCut::kNeutralAndCharged )
818 // --------------------------- //
819 // Tracks and clusters in cone //
820 // --------------------------- //
822 Double_t sumPhiUESub = coneptsumClusterSubPhi + coneptsumTrackSubPhi;
823 Double_t sumEtaUESub = coneptsumClusterSubEta + coneptsumTrackSubEta;
825 fhConeSumPtPhiUESub ->Fill(ptTrig , sumPhiUESub);
826 fhConeSumPtPhiUESubTrigEtaPhi->Fill(etaTrig, phiTrig, sumPhiUESub);
827 fhConeSumPtEtaUESub ->Fill(ptTrig , sumEtaUESub);
828 fhConeSumPtEtaUESubTrigEtaPhi->Fill(etaTrig, phiTrig, sumEtaUESub);
830 fhEtaBandClustervsTrack ->Fill(etaUEptsumCluster ,etaUEptsumTrack );
831 fhPhiBandClustervsTrack ->Fill(phiUEptsumCluster ,phiUEptsumTrack );
832 fhEtaBandNormClustervsTrack->Fill(etaUEptsumClusterNorm,etaUEptsumTrackNorm);
833 fhPhiBandNormClustervsTrack->Fill(phiUEptsumClusterNorm,phiUEptsumTrackNorm);
835 fhConeSumPtEtaUESubClustervsTrack->Fill(coneptsumClusterSubEta,coneptsumTrackSubEta);
836 fhConeSumPtPhiUESubClustervsTrack->Fill(coneptsumClusterSubPhi,coneptsumTrackSubPhi);
838 // ------------------------ //
839 // Tracks and cells in cone //
840 // ------------------------ //
842 if(fFillCellHistograms)
844 Double_t sumPhiUESubTrackCell = coneptsumCellSubPhi + coneptsumTrackSubPhi;
845 Double_t sumEtaUESubTrackCell = coneptsumCellSubEta + coneptsumTrackSubEta;
847 fhConeSumPtPhiUESubTrackCell ->Fill(ptTrig , sumPhiUESubTrackCell);
848 fhConeSumPtPhiUESubTrackCellTrigEtaPhi->Fill(etaTrig, phiTrig, sumPhiUESubTrackCell);
849 fhConeSumPtEtaUESubTrackCell ->Fill(ptTrig , sumEtaUESubTrackCell);
850 fhConeSumPtEtaUESubTrackCellTrigEtaPhi->Fill(etaTrig, phiTrig, sumEtaUESubTrackCell);
852 fhEtaBandCellvsTrack ->Fill(etaUEptsumCell ,etaUEptsumTrack );
853 fhPhiBandCellvsTrack ->Fill(phiUEptsumCell ,phiUEptsumTrack );
854 fhEtaBandNormCellvsTrack->Fill(etaUEptsumCellNorm,etaUEptsumTrackNorm);
855 fhPhiBandNormCellvsTrack->Fill(phiUEptsumCellNorm,phiUEptsumTrackNorm);
857 fhConeSumPtEtaUESubCellvsTrack->Fill(coneptsumCellSubEta,coneptsumTrackSubEta);
858 fhConeSumPtPhiUESubCellvsTrack->Fill(coneptsumCellSubPhi,coneptsumTrackSubPhi);
865 //______________________________________________________________________________________________________________
866 void AliAnaParticleIsolation::CalculateCaloSignalInCone(AliAODPWG4ParticleCorrelation * aodParticle,
867 Float_t & coneptsumCluster, Float_t & coneptLeadCluster)
869 // Get the cluster pT or sum of pT in isolation cone
870 coneptLeadCluster = 0;
871 coneptsumCluster = 0;
873 if( GetIsolationCut()->GetParticleTypeInCone()==AliIsolationCut::kOnlyCharged ) return ;
875 //Recover reference arrays with clusters and tracks
876 TObjArray * refclusters = aodParticle->GetObjArray(GetAODObjArrayName()+"Clusters");
877 if(!refclusters) return ;
879 Float_t ptTrig = aodParticle->Pt();
881 //Get vertex for cluster momentum calculation
882 Double_t vertex[] = {0,0,0} ; //vertex ;
883 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
884 GetReader()->GetVertex(vertex);
887 for(Int_t icalo=0; icalo < refclusters->GetEntriesFast(); icalo++)
889 AliVCluster* calo = (AliVCluster *) refclusters->At(icalo);
890 calo->GetMomentum(mom,vertex) ;//Assume that come from vertex in straight line
892 fhPtInCone ->Fill(ptTrig, mom.Pt());
893 fhPtClusterInCone->Fill(ptTrig, mom.Pt());
895 if(fFillPileUpHistograms)
897 if(GetReader()->IsPileUpFromSPD()) fhPtInConePileUp[0]->Fill(ptTrig,mom.Pt());
898 if(GetReader()->IsPileUpFromEMCal()) fhPtInConePileUp[1]->Fill(ptTrig,mom.Pt());
899 if(GetReader()->IsPileUpFromSPDOrEMCal()) fhPtInConePileUp[2]->Fill(ptTrig,mom.Pt());
900 if(GetReader()->IsPileUpFromSPDAndEMCal()) fhPtInConePileUp[3]->Fill(ptTrig,mom.Pt());
901 if(GetReader()->IsPileUpFromSPDAndNotEMCal()) fhPtInConePileUp[4]->Fill(ptTrig,mom.Pt());
902 if(GetReader()->IsPileUpFromEMCalAndNotSPD()) fhPtInConePileUp[5]->Fill(ptTrig,mom.Pt());
903 if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) fhPtInConePileUp[6]->Fill(ptTrig,mom.Pt());
906 if(fFillHighMultHistograms) fhPtInConeCent->Fill(GetEventCentrality(),mom.Pt());
908 coneptsumCluster+=mom.Pt();
909 if(mom.Pt() > coneptLeadCluster) coneptLeadCluster = mom.Pt();
912 fhConeSumPtCluster ->Fill(ptTrig, coneptsumCluster );
913 fhConePtLeadCluster->Fill(ptTrig, coneptLeadCluster);
916 //______________________________________________________________________________________________________
917 void AliAnaParticleIsolation::CalculateCaloCellSignalInCone(AliAODPWG4ParticleCorrelation * aodParticle,
918 Float_t & coneptsumCell)
920 // Get the cell amplityde or sum of amplitudes in isolation cone
921 // Mising: Remove signal cells in cone in case the trigger is a cluster!
923 if( GetIsolationCut()->GetParticleTypeInCone()==AliIsolationCut::kOnlyCharged ) return ;
925 Float_t conesize = GetIsolationCut()->GetConeSize();
927 Float_t ptTrig = aodParticle->Pt();
928 Float_t phiTrig = aodParticle->Phi();
929 if(phiTrig<0) phiTrig += TMath::TwoPi();
930 Float_t etaTrig = aodParticle->Eta();
932 if(aodParticle->GetDetector()=="EMCAL")
934 AliEMCALGeometry* eGeom = AliEMCALGeometry::GetInstance();
937 if (eGeom->GetAbsCellIdFromEtaPhi(etaTrig,phiTrig,absId))
939 if(!eGeom->CheckAbsCellId(absId)) return ;
941 // Get absolute (col,row) of trigger particle
942 Int_t nSupMod = eGeom->GetSuperModuleNumber(absId);
944 Int_t imEta=-1, imPhi=-1;
945 Int_t ieta =-1, iphi =-1;
947 if (eGeom->GetCellIndex(absId,nSupMod,nModule,imPhi,imEta))
949 Int_t iEta=-1, iPhi=-1;
950 eGeom->GetCellPhiEtaIndexInSModule(nSupMod,nModule,imPhi,imEta,iphi,ieta);
952 Int_t colTrig = iEta;
953 if (nSupMod % 2) colTrig = AliEMCALGeoParams::fgkEMCALCols + iEta ;
954 Int_t rowTrig = iPhi + AliEMCALGeoParams::fgkEMCALRows*int(nSupMod/2);
956 Int_t sqrSize = int(conesize/0.0143);
958 AliVCaloCells * cells = GetEMCALCells();
960 // Loop on cells in cone
961 for(Int_t irow = rowTrig-sqrSize; irow < rowTrig+sqrSize; irow++)
963 for(Int_t icol = colTrig-sqrSize; icol < colTrig+sqrSize; icol++)
965 Int_t inSector = int(irow/AliEMCALGeoParams::fgkEMCALRows);
966 if(inSector==5) continue;
970 if(icol < AliEMCALGeoParams::fgkEMCALCols)
972 inSupMod = 2*inSector + 1;
975 else if(icol > AliEMCALGeoParams::fgkEMCALCols - 1)
977 inSupMod = 2*inSector;
978 icolLoc = icol-AliEMCALGeoParams::fgkEMCALCols;
981 Int_t irowLoc = irow - AliEMCALGeoParams::fgkEMCALRows*inSector ;
983 Int_t iabsId = eGeom->GetAbsCellIdFromCellIndexes(inSupMod,irowLoc,icolLoc);
984 if(!eGeom->CheckAbsCellId(iabsId)) continue;
986 fhPtCellInCone->Fill(ptTrig, cells->GetCellAmplitude(iabsId));
987 coneptsumCell += cells->GetCellAmplitude(iabsId);
994 fhConeSumPtCell->Fill(ptTrig,coneptsumCell);
998 //___________________________________________________________________________________________________________
999 void AliAnaParticleIsolation::CalculateTrackSignalInCone(AliAODPWG4ParticleCorrelation * aodParticle,
1000 Float_t & coneptsumTrack, Float_t & coneptLeadTrack)
1002 // Get the track pT or sum of pT in isolation cone
1004 if( GetIsolationCut()->GetParticleTypeInCone()==AliIsolationCut::kOnlyNeutral ) return ;
1006 //Recover reference arrays with clusters and tracks
1007 TObjArray * reftracks = aodParticle->GetObjArray(GetAODObjArrayName()+"Tracks");
1008 if(!reftracks) return ;
1010 Float_t ptTrig = aodParticle->Pt();
1011 Double_t bz = GetReader()->GetInputEvent()->GetMagneticField();
1013 for(Int_t itrack=0; itrack < reftracks->GetEntriesFast(); itrack++)
1015 AliVTrack* track = (AliVTrack *) reftracks->At(itrack);
1016 Float_t pTtrack = track->Pt();
1018 fhPtInCone ->Fill(ptTrig,pTtrack);
1019 fhPtTrackInCone->Fill(ptTrig,pTtrack);
1021 if(fFillPileUpHistograms)
1023 ULong_t status = track->GetStatus();
1024 Bool_t okTOF = ( (status & AliVTrack::kTOFout) == AliVTrack::kTOFout ) ;
1025 //Double32_t tof = track->GetTOFsignal()*1e-3;
1026 Int_t trackBC = track->GetTOFBunchCrossing(bz);
1028 if ( okTOF && trackBC!=0 ) fhPtTrackInConeOtherBC->Fill(ptTrig,pTtrack);
1029 else if( okTOF && trackBC==0 ) fhPtTrackInConeBC0 ->Fill(ptTrig,pTtrack);
1031 Int_t vtxBC = GetReader()->GetVertexBC();
1032 if(vtxBC == 0 || vtxBC==AliVTrack::kTOFBCNA) fhPtTrackInConeVtxBC0->Fill(ptTrig,pTtrack);
1034 if(GetReader()->IsPileUpFromSPD()) { fhPtInConePileUp[0]->Fill(ptTrig,pTtrack);
1035 if(okTOF && trackBC!=0 ) fhPtTrackInConeOtherBCPileUpSPD->Fill(ptTrig,pTtrack);
1036 if(okTOF && trackBC==0 ) fhPtTrackInConeBC0PileUpSPD ->Fill(ptTrig,pTtrack); }
1037 if(GetReader()->IsPileUpFromEMCal()) fhPtInConePileUp[1]->Fill(ptTrig,pTtrack);
1038 if(GetReader()->IsPileUpFromSPDOrEMCal()) fhPtInConePileUp[2]->Fill(ptTrig,pTtrack);
1039 if(GetReader()->IsPileUpFromSPDAndEMCal()) fhPtInConePileUp[3]->Fill(ptTrig,pTtrack);
1040 if(GetReader()->IsPileUpFromSPDAndNotEMCal()) fhPtInConePileUp[4]->Fill(ptTrig,pTtrack);
1041 if(GetReader()->IsPileUpFromEMCalAndNotSPD()) fhPtInConePileUp[5]->Fill(ptTrig,pTtrack);
1042 if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) fhPtInConePileUp[6]->Fill(ptTrig,pTtrack);
1045 if(fFillHighMultHistograms) fhPtInConeCent->Fill(GetEventCentrality(),pTtrack);
1047 coneptsumTrack+=pTtrack;
1048 if(pTtrack > coneptLeadTrack) coneptLeadTrack = pTtrack;
1051 fhConeSumPtTrack ->Fill(ptTrig, coneptsumTrack );
1052 fhConePtLeadTrack->Fill(ptTrig, coneptLeadTrack);
1056 //_________________________________________________________________
1057 void AliAnaParticleIsolation::FillPileUpHistograms(Int_t clusterID)
1059 // Fill some histograms to understand pile-up
1063 printf("AliAnaParticleIsolation::FillPileUpHistograms(), ID of cluster = %d, not possible! ", clusterID);
1068 TObjArray* clusters = 0x0;
1069 if (fCalorimeter == "EMCAL") clusters = GetEMCALClusters();
1070 else if(fCalorimeter == "PHOS" ) clusters = GetPHOSClusters();
1073 Float_t time = -1000;
1077 AliVCluster *cluster = FindCluster(clusters,clusterID,iclus);
1078 energy = cluster->E();
1079 time = cluster->GetTOF()*1e9;
1082 //printf("E %f, time %f\n",energy,time);
1083 AliVEvent * event = GetReader()->GetInputEvent();
1085 fhTimeENoCut->Fill(energy,time);
1086 if(GetReader()->IsPileUpFromSPD()) fhTimeESPD ->Fill(energy,time);
1087 if(event->IsPileupFromSPDInMultBins()) fhTimeESPDMulti->Fill(energy,time);
1089 if(energy < 8) return; // Fill time figures for high energy clusters not too close to trigger threshold
1091 AliESDEvent* esdEv = dynamic_cast<AliESDEvent*> (event);
1092 AliAODEvent* aodEv = dynamic_cast<AliAODEvent*> (event);
1094 // N pile up vertices
1095 Int_t nVerticesSPD = -1;
1096 Int_t nVerticesTracks = -1;
1100 nVerticesSPD = esdEv->GetNumberOfPileupVerticesSPD();
1101 nVerticesTracks = esdEv->GetNumberOfPileupVerticesTracks();
1106 nVerticesSPD = aodEv->GetNumberOfPileupVerticesSPD();
1107 nVerticesTracks = aodEv->GetNumberOfPileupVerticesTracks();
1110 fhTimeNPileUpVertSPD ->Fill(time,nVerticesSPD);
1111 fhTimeNPileUpVertTrack->Fill(time,nVerticesTracks);
1113 //printf("Is SPD %d, Is SPD Multi %d, n spd %d, n track %d\n",
1114 // GetReader()->IsPileUpFromSPD(),event->IsPileupFromSPDInMultBins(),nVerticesSPD,nVerticesTracks);
1117 Float_t z1 = -1, z2 = -1;
1119 for(Int_t iVert=0; iVert<nVerticesSPD;iVert++)
1123 const AliESDVertex* pv=esdEv->GetPileupVertexSPD(iVert);
1124 ncont=pv->GetNContributors();
1125 z1 = esdEv->GetPrimaryVertexSPD()->GetZ();
1127 diamZ = esdEv->GetDiamondZ();
1131 AliAODVertex *pv=aodEv->GetVertex(iVert);
1132 if(pv->GetType()!=AliAODVertex::kPileupSPD) continue;
1133 ncont=pv->GetNContributors();
1134 z1=aodEv->GetPrimaryVertexSPD()->GetZ();
1136 diamZ = aodEv->GetDiamondZ();
1139 Double_t distZ = TMath::Abs(z2-z1);
1140 diamZ = TMath::Abs(z2-diamZ);
1142 fhTimeNPileUpVertContributors ->Fill(time,ncont);
1143 fhTimePileUpMainVertexZDistance->Fill(time,distZ);
1144 fhTimePileUpMainVertexZDiamond ->Fill(time,diamZ);
1149 //_____________________________________________________________________________________________________________________
1150 void AliAnaParticleIsolation::FillTrackMatchingShowerShapeControlHistograms(AliAODPWG4ParticleCorrelation *pCandidate,
1151 Float_t coneptsum, Float_t coneleadpt,
1154 // Fill Track matching and Shower Shape control histograms
1155 if(!fFillTMHisto && !fFillSSHisto && !fFillBackgroundBinHistograms) return;
1157 Int_t clusterID = pCandidate->GetCaloLabel(0) ;
1158 Int_t nMaxima = pCandidate->GetFiducialArea(); // bad name, just place holder for the moment
1159 Int_t mcTag = pCandidate->GetTag() ;
1160 Bool_t isolated = pCandidate->IsIsolated();
1164 printf("AliAnaParticleIsolation::FillTrackMatchingShowerShapeControlHistograms(), ID of cluster = %d, not possible! \n", clusterID);
1169 TObjArray* clusters = 0x0;
1170 if (fCalorimeter == "EMCAL") clusters = GetEMCALClusters();
1171 else if(fCalorimeter == "PHOS" ) clusters = GetPHOSClusters();
1173 if(!clusters) return;
1175 AliVCluster *cluster = FindCluster(clusters,clusterID,iclus);
1177 Float_t m02 = cluster->GetM02() ;
1178 Float_t energy = pCandidate->E();
1179 Float_t pt = pCandidate->Pt();
1181 // Get the max pt leading in cone or the sum of pt in cone
1182 // assign a bin to the candidate, depending on both quantities
1183 // see the shower shape in those bins.
1184 if(fFillBackgroundBinHistograms)
1186 // Get the background bin for this cone and trigger
1187 Int_t ptsumBin = -1;
1188 Int_t leadptBin = -1;
1190 if( GetDebug() > 1 )
1191 printf("AliAnaParticleIsolation::FillTrackMatchingShowerShapeControlHistograms() - pT cand: %2.2f, In cone pT: Sum %2.2f, Lead %2.2f, n bins %d\n",
1192 pt,coneptsum,coneleadpt,fNBkgBin);
1194 for(Int_t ibin = 0; ibin < fNBkgBin; ibin++)
1196 if( coneptsum >= fBkgBinLimit[ibin] && coneptsum < fBkgBinLimit[ibin+1]) ptsumBin = ibin;
1197 if( coneleadpt >= fBkgBinLimit[ibin] && coneleadpt < fBkgBinLimit[ibin+1]) leadptBin = ibin;
1200 // Fill the histograms per bin of pt lead or pt sum
1201 if( GetDebug() > 1 && ptsumBin >=0 ) printf("\t Sum bin %d [%2.2f,%2.2f]\n" , ptsumBin ,fBkgBinLimit[ptsumBin] ,fBkgBinLimit[ptsumBin +1]);
1202 if( GetDebug() > 1 && leadptBin >=0 ) printf("\t Lead bin %d [%2.2f,%2.2f]\n", leadptBin,fBkgBinLimit[leadptBin],fBkgBinLimit[leadptBin+1]);
1204 if( leadptBin >=0 ) fhPtLeadConeBinLambda0[leadptBin]->Fill(pt,m02);
1205 if( ptsumBin >=0 ) fhSumPtConeBinLambda0 [ ptsumBin]->Fill(pt,m02);
1207 if( GetDebug() > 1 && leadptBin == 0 )
1208 printf("No track/clusters in isolation cone: cand pt %2.2f GeV/c, track multiplicity %d, N clusters %d\n",
1209 pt, GetTrackMultiplicity(),GetEMCALClusters()->GetEntriesFast());
1213 Int_t leadptBinMC = leadptBin+mcIndex*fNBkgBin;
1214 Int_t ptsumBinMC = ptsumBin+mcIndex*fNBkgBin;
1215 if( leadptBin >=0 ) fhPtLeadConeBinLambda0MC[leadptBinMC]->Fill(pt,m02);
1216 if( ptsumBin >=0 ) fhSumPtConeBinLambda0MC [ ptsumBinMC]->Fill(pt,m02);
1217 if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton))
1219 leadptBinMC = leadptBin+kmcPhoton*fNBkgBin;
1220 ptsumBinMC = ptsumBin+kmcPhoton*fNBkgBin;
1221 if( leadptBin >=0 ) fhPtLeadConeBinLambda0MC[leadptBinMC]->Fill(pt,m02);
1222 if( ptsumBin >=0 ) fhSumPtConeBinLambda0MC [ ptsumBinMC]->Fill(pt,m02);
1229 fhELambda0 [isolated]->Fill(energy, m02);
1230 fhPtLambda0[isolated]->Fill(pt, m02);
1231 fhELambda1 [isolated]->Fill(energy, m02);
1232 if(fFillTaggedDecayHistograms)
1234 Int_t decayTag = pCandidate->GetBtag(); // temporary
1235 if(decayTag < 0) decayTag = 0; // temporary
1236 for(Int_t ibit = 0; ibit < fNDecayBits; ibit++)
1238 if(GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[ibit]))
1239 fhPtLambda0Decay[isolated][ibit]->Fill(pt,m02);
1245 if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton))
1246 fhPtLambda0MC[kmcPhoton][isolated]->Fill(pt,m02);
1248 if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCDecayPairLost) && mcIndex==kmcPi0Decay )
1249 fhPtLambda0MC[kmcPi0DecayLostPair][isolated]->Fill(pt,m02);
1251 fhPtLambda0MC[mcIndex][isolated]->Fill(pt,m02);
1254 if(fCalorimeter == "EMCAL" && GetFirstSMCoveredByTRD() >= 0 &&
1255 GetModuleNumber(cluster) >= GetFirstSMCoveredByTRD() )
1257 fhELambda0TRD [isolated]->Fill(energy, m02 );
1258 fhPtLambda0TRD[isolated]->Fill(pt , m02 );
1259 fhELambda1TRD [isolated]->Fill(energy, m02 );
1262 if(fFillNLMHistograms)
1264 fhNLocMax[isolated]->Fill(energy,nMaxima);
1265 if (nMaxima==1) { fhELambda0LocMax1[isolated]->Fill(energy,m02); fhELambda1LocMax1[isolated]->Fill(energy,m02); }
1266 else if(nMaxima==2) { fhELambda0LocMax2[isolated]->Fill(energy,m02); fhELambda1LocMax2[isolated]->Fill(energy,m02); }
1267 else { fhELambda0LocMaxN[isolated]->Fill(energy,m02); fhELambda1LocMaxN[isolated]->Fill(energy,m02); }
1273 Float_t dZ = cluster->GetTrackDz();
1274 Float_t dR = cluster->GetTrackDx();
1276 if(cluster->IsEMCAL() && GetCaloUtils()->IsRecalculationOfClusterTrackMatchingOn())
1278 dR = 2000., dZ = 2000.;
1279 GetCaloUtils()->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
1282 //printf("ParticleIsolation: dPhi %f, dEta %f\n",dR,dZ);
1283 if(fhTrackMatchedDEta[isolated] && TMath::Abs(dR) < 999)
1285 fhTrackMatchedDEta[isolated]->Fill(energy,dZ);
1286 fhTrackMatchedDPhi[isolated]->Fill(energy,dR);
1287 if(energy > 0.5) fhTrackMatchedDEtaDPhi[isolated]->Fill(dZ,dR);
1290 // Check dEdx and E/p of matched clusters
1292 if(TMath::Abs(dZ) < 0.05 && TMath::Abs(dR) < 0.05)
1295 AliVTrack *track = GetCaloUtils()->GetMatchedTrack(cluster, GetReader()->GetInputEvent());
1299 Float_t dEdx = track->GetTPCsignal();
1300 fhdEdx[isolated]->Fill(cluster->E(), dEdx);
1302 Float_t eOverp = cluster->E()/track->P();
1303 fhEOverP[isolated]->Fill(cluster->E(), eOverp);
1306 // printf("AliAnaParticleIsolation::FillTrackMatchingShowerShapeHistograms() - Residual OK but (dR, dZ)= (%2.4f,%2.4f) no track associated WHAT? \n", dR,dZ);
1311 if ( !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion) )
1313 if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0) ||
1314 GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta) ) fhTrackMatchedMCParticle[isolated]->Fill(energy, 2.5 );
1315 else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton) ) fhTrackMatchedMCParticle[isolated]->Fill(energy, 0.5 );
1316 else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron) ) fhTrackMatchedMCParticle[isolated]->Fill(energy, 1.5 );
1317 else fhTrackMatchedMCParticle[isolated]->Fill(energy, 3.5 );
1322 if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0) ||
1323 GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta) ) fhTrackMatchedMCParticle[isolated]->Fill(energy, 6.5 );
1324 else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton) ) fhTrackMatchedMCParticle[isolated]->Fill(energy, 4.5 );
1325 else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron) ) fhTrackMatchedMCParticle[isolated]->Fill(energy, 5.5 );
1326 else fhTrackMatchedMCParticle[isolated]->Fill(energy, 7.5 );
1336 //______________________________________________________
1337 TObjString * AliAnaParticleIsolation::GetAnalysisCuts()
1339 //Save parameters used for analysis
1340 TString parList ; //this will be list of parameters used for this analysis.
1341 const Int_t buffersize = 255;
1342 char onePar[buffersize] ;
1344 snprintf(onePar, buffersize,"--- AliAnaParticleIsolation ---\n") ;
1346 snprintf(onePar, buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
1348 snprintf(onePar, buffersize,"Isolation Cand Detector: %s\n",fIsoDetector.Data()) ;
1350 snprintf(onePar, buffersize,"fReMakeIC =%d (Flag for reisolation during histogram filling) \n",fReMakeIC) ;
1352 snprintf(onePar, buffersize,"fMakeSeveralIC=%d (Flag for isolation with several cuts at the same time ) \n",fMakeSeveralIC) ;
1354 snprintf(onePar, buffersize,"fFillTMHisto=%d (Flag for track matching histograms) \n",fFillTMHisto) ;
1356 snprintf(onePar, buffersize,"fFillSSHisto=%d (Flag for shower shape histograms) \n",fFillSSHisto) ;
1361 snprintf(onePar, buffersize,"fNCones =%d (Number of cone sizes) \n",fNCones) ;
1363 snprintf(onePar, buffersize,"fNPtThresFrac=%d (Flag for isolation with several cuts at the same time ) \n",fNPtThresFrac) ;
1366 for(Int_t icone = 0; icone < fNCones ; icone++)
1368 snprintf(onePar, buffersize,"fConeSizes[%d]=%1.2f (isolation cone size) \n",icone, fConeSizes[icone]) ;
1371 for(Int_t ipt = 0; ipt < fNPtThresFrac ; ipt++)
1373 snprintf(onePar, buffersize,"fPtThresholds[%d]=%1.2f (isolation pt threshold) \n",ipt, fPtThresholds[ipt]) ;
1376 for(Int_t ipt = 0; ipt < fNPtThresFrac ; ipt++)
1378 snprintf(onePar, buffersize,"fPtFractions[%d]=%1.2f (isolation pt fraction threshold) \n",ipt, fPtFractions[ipt]) ;
1381 for(Int_t ipt = 0; ipt < fNPtThresFrac ; ipt++)
1383 snprintf(onePar, buffersize,"fSumPtThresholds[%d]=%1.2f (isolation sum pt threshold) \n",ipt, fSumPtThresholds[ipt]) ;
1388 //Get parameters set in base class.
1389 parList += GetBaseParametersList() ;
1391 //Get parameters set in IC class.
1392 if(!fMakeSeveralIC)parList += GetIsolationCut()->GetICParametersList() ;
1394 return new TObjString(parList) ;
1398 //________________________________________________________
1399 TList * AliAnaParticleIsolation::GetCreateOutputObjects()
1401 // Create histograms to be saved in output file and
1402 // store them in outputContainer
1403 TList * outputContainer = new TList() ;
1404 outputContainer->SetName("IsolatedParticleHistos") ;
1406 Int_t nptbins = GetHistogramRanges()->GetHistoPtBins();
1407 Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins();
1408 Int_t netabins = GetHistogramRanges()->GetHistoEtaBins();
1409 Float_t ptmax = GetHistogramRanges()->GetHistoPtMax();
1410 Float_t phimax = GetHistogramRanges()->GetHistoPhiMax();
1411 Float_t etamax = GetHistogramRanges()->GetHistoEtaMax();
1412 Float_t ptmin = GetHistogramRanges()->GetHistoPtMin();
1413 Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();
1414 Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
1415 Int_t ssbins = GetHistogramRanges()->GetHistoShowerShapeBins();
1416 Float_t ssmax = GetHistogramRanges()->GetHistoShowerShapeMax();
1417 Float_t ssmin = GetHistogramRanges()->GetHistoShowerShapeMin();
1418 Int_t ntimebins= GetHistogramRanges()->GetHistoTimeBins();
1419 Float_t timemax = GetHistogramRanges()->GetHistoTimeMax();
1420 Float_t timemin = GetHistogramRanges()->GetHistoTimeMin();
1422 Int_t nresetabins = GetHistogramRanges()->GetHistoTrackResidualEtaBins();
1423 Float_t resetamax = GetHistogramRanges()->GetHistoTrackResidualEtaMax();
1424 Float_t resetamin = GetHistogramRanges()->GetHistoTrackResidualEtaMin();
1425 Int_t nresphibins = GetHistogramRanges()->GetHistoTrackResidualPhiBins();
1426 Float_t resphimax = GetHistogramRanges()->GetHistoTrackResidualPhiMax();
1427 Float_t resphimin = GetHistogramRanges()->GetHistoTrackResidualPhiMin();
1429 Int_t ndedxbins = GetHistogramRanges()->GetHistodEdxBins();
1430 Float_t dedxmax = GetHistogramRanges()->GetHistodEdxMax();
1431 Float_t dedxmin = GetHistogramRanges()->GetHistodEdxMin();
1432 Int_t nPoverEbins = GetHistogramRanges()->GetHistoPOverEBins();
1433 Float_t pOverEmax = GetHistogramRanges()->GetHistoPOverEMax();
1434 Float_t pOverEmin = GetHistogramRanges()->GetHistoPOverEMin();
1436 Int_t nptsumbins = GetHistogramRanges()->GetHistoNPtSumBins();
1437 Float_t ptsummax = GetHistogramRanges()->GetHistoPtSumMax();
1438 Float_t ptsummin = GetHistogramRanges()->GetHistoPtSumMin();
1439 Int_t nptinconebins = GetHistogramRanges()->GetHistoNPtInConeBins();
1440 Float_t ptinconemax = GetHistogramRanges()->GetHistoPtInConeMax();
1441 Float_t ptinconemin = GetHistogramRanges()->GetHistoPtInConeMin();
1443 //Float_t ptthre = GetIsolationCut()->GetPtThreshold();
1444 //Float_t ptsumthre = GetIsolationCut()->GetSumPtThreshold();
1445 //Float_t ptfrac = GetIsolationCut()->GetPtFraction();
1446 Float_t r = GetIsolationCut()->GetConeSize();
1447 Int_t method = GetIsolationCut()->GetICMethod() ;
1448 Int_t particle = GetIsolationCut()->GetParticleTypeInCone() ;
1450 TString sThreshold = "";
1451 if ( method == AliIsolationCut::kSumPtIC )
1453 sThreshold = Form(", %2.2f < #Sigma #it{p}_{T}^{in cone} < %2.2f GeV/#it{c}",
1454 GetIsolationCut()->GetSumPtThreshold(), GetIsolationCut()->GetSumPtThresholdMax());
1455 if(GetIsolationCut()->GetSumPtThresholdMax() > 200)
1456 sThreshold = Form(", #Sigma #it{p}_{T}^{in cone} = %2.2f GeV/#it{c}",
1457 GetIsolationCut()->GetSumPtThreshold());
1459 else if ( method == AliIsolationCut::kPtThresIC)
1461 sThreshold = Form(", %2.2f < #it{p}_{T}^{th} < %2.2f GeV/#it{c}",
1462 GetIsolationCut()->GetPtThreshold(),GetIsolationCut()->GetPtThresholdMax());
1463 if(GetIsolationCut()->GetSumPtThreshold() > 200)
1464 sThreshold = Form(", #it{p}_{T}^{th} = %2.2f GeV/#it{c}",
1465 GetIsolationCut()->GetPtThreshold());
1467 else if ( method == AliIsolationCut::kPtFracIC)
1468 sThreshold = Form(", #Sigma #it{p}_{T}^{in cone}/#it{p}_{T}^{trig} = %2.2f" ,
1469 GetIsolationCut()->GetPtFraction());
1471 TString sParticle = ", x^{0,#pm}";
1472 if ( particle == AliIsolationCut::kOnlyNeutral ) sParticle = ", x^{0}";
1473 else if ( particle == AliIsolationCut::kOnlyCharged ) sParticle = ", x^{#pm}";
1475 TString parTitle = Form("#it{R} = %2.2f%s%s",GetIsolationCut()->GetConeSize(), sThreshold.Data(),sParticle.Data());
1477 TString pileUpName[] = {"SPD","EMCAL","SPDOrEMCAL","SPDAndEMCAL","SPDAndNotEMCAL","EMCALAndNotSPD","NotSPDAndNotEMCAL"} ;
1479 // MC histograms title and name
1480 TString mcPartType[] = { "#gamma", "#gamma_{prompt}", "#gamma_{fragmentation}",
1481 "#pi^{0} (merged #gamma)","#gamma_{#pi decay}","#gamma_{#pi decay} lost companion",
1482 "#gamma_{#eta decay}","#gamma_{other decay}",
1483 "e^{#pm}","hadrons?"} ;
1485 TString mcPartName[] = { "Photon","PhotonPrompt","PhotonFrag",
1486 "Pi0","Pi0Decay","Pi0DecayLostPair","EtaDecay","OtherDecay",
1487 "Electron","Hadron"} ;
1489 // Primary MC histograms title and name
1490 TString pptype[] = { "#gamma", "#gamma_{#pi decay}","#gamma_{other decay}",
1491 "#gamma_{prompt}","#gamma_{fragmentation}","#gamma_{ISR}","#pi^{0}"} ;
1493 TString ppname[] = { "Photon","PhotonPi0Decay","PhotonOtherDecay",
1494 "PhotonPrompt","PhotonFrag","PhotonISR","Pi0"} ;
1496 // Not Isolated histograms, reference histograms
1498 fhENoIso = new TH1F("hENoIso",
1499 Form("Number of not isolated leading particles vs #it{p}_{T}, %s",parTitle.Data()),
1500 nptbins,ptmin,ptmax);
1501 fhENoIso->SetYTitle("#it{counts}");
1502 fhENoIso->SetXTitle("E (GeV/#it{c})");
1503 outputContainer->Add(fhENoIso) ;
1505 fhPtNoIso = new TH1F("hPtNoIso",
1506 Form("Number of not isolated leading particles vs #it{p}_{T}, %s",parTitle.Data()),
1507 nptbins,ptmin,ptmax);
1508 fhPtNoIso->SetYTitle("#it{counts}");
1509 fhPtNoIso->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1510 outputContainer->Add(fhPtNoIso) ;
1512 fhEtaPhiNoIso = new TH2F("hEtaPhiNoIso",
1513 Form("Number of not isolated leading particles #eta vs #phi, %s",parTitle.Data()),
1514 netabins,etamin,etamax,nphibins,phimin,phimax);
1515 fhEtaPhiNoIso->SetXTitle("#eta");
1516 fhEtaPhiNoIso->SetYTitle("#phi");
1517 outputContainer->Add(fhEtaPhiNoIso) ;
1521 // For histograms in arrays, index in the array, corresponding to any particle origin
1523 for(Int_t imc = 0; imc < fgkNmcTypes; imc++)
1525 fhPtNoIsoMC[imc] = new TH1F(Form("hPtNoIsoMC%s",mcPartName[imc].Data()),
1526 Form("#it{p}_{T} of NOT isolated %s, %s",mcPartType[imc].Data(),parTitle.Data()),
1527 nptbins,ptmin,ptmax);
1528 fhPtNoIsoMC[imc]->SetYTitle("#it{counts}");
1529 fhPtNoIsoMC[imc]->SetXTitle("#it{p}_{T}(GeV/#it{c})");
1530 outputContainer->Add(fhPtNoIsoMC[imc]) ;
1532 fhPtIsoMC[imc] = new TH1F(Form("hPtMC%s",mcPartName[imc].Data()),
1533 Form("#it{p}_{T} of isolated %s, %s",mcPartType[imc].Data(),parTitle.Data()),
1534 nptbins,ptmin,ptmax);
1535 fhPtIsoMC[imc]->SetYTitle("#it{counts}");
1536 fhPtIsoMC[imc]->SetXTitle("#it{p}_{T}(GeV/#it{c})");
1537 outputContainer->Add(fhPtIsoMC[imc]) ;
1539 fhPhiIsoMC[imc] = new TH2F(Form("hPhiMC%s",mcPartName[imc].Data()),
1540 Form("#phi vs #it{p}_{T} of isolated %s, %s",mcPartType[imc].Data(),parTitle.Data()),
1541 nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1542 fhPhiIsoMC[imc]->SetYTitle("#phi");
1543 fhPhiIsoMC[imc]->SetXTitle("#it{p}_{T}(GeV/#it{c})");
1544 outputContainer->Add(fhPhiIsoMC[imc]) ;
1546 fhEtaIsoMC[imc] = new TH2F(Form("hEtaMC%s",mcPartName[imc].Data()),
1547 Form("#phi vs #it{p}_{T} of isolated %s, %s",mcPartType[imc].Data(),parTitle.Data()),
1548 nptbins,ptmin,ptmax,netabins,etamin,etamax);
1549 fhEtaIsoMC[imc]->SetYTitle("#eta");
1550 fhEtaIsoMC[imc]->SetXTitle("#it{p}_{T}(GeV/#it{c})");
1551 outputContainer->Add(fhEtaIsoMC[imc]) ;
1555 // Histograms for tagged candidates as decay
1556 if(fFillTaggedDecayHistograms)
1558 for(Int_t ibit = 0; ibit < fNDecayBits; ibit++)
1560 fhPtDecayNoIso[ibit] =
1561 new TH1F(Form("hPtDecayNoIso_bit%d",fDecayBits[ibit]),
1562 Form("Number of not isolated leading pi0 decay particles vs #it{p}_{T}, bit %d, %s",fDecayBits[ibit],parTitle.Data()),
1563 nptbins,ptmin,ptmax);
1564 fhPtDecayNoIso[ibit]->SetYTitle("#it{counts}");
1565 fhPtDecayNoIso[ibit]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1566 outputContainer->Add(fhPtDecayNoIso[ibit]) ;
1568 fhEtaPhiDecayNoIso[ibit] =
1569 new TH2F(Form("hEtaPhiDecayNoIso_bit%d",fDecayBits[ibit]),
1570 Form("Number of not isolated leading Pi0 decay particles #eta vs #phi, bit %d, %s",fDecayBits[ibit],parTitle.Data()),
1571 netabins,etamin,etamax,nphibins,phimin,phimax);
1572 fhEtaPhiDecayNoIso[ibit]->SetXTitle("#eta");
1573 fhEtaPhiDecayNoIso[ibit]->SetYTitle("#phi");
1574 outputContainer->Add(fhEtaPhiDecayNoIso[ibit]) ;
1578 fhPtDecayIso[ibit] =
1579 new TH1F(Form("hPtDecayIso_bit%d",fDecayBits[ibit]),
1580 Form("Number of isolated #pi^{0} decay particles vs #it{p}_{T}, bit %d, %s",fDecayBits[ibit],parTitle.Data()),
1581 nptbins,ptmin,ptmax);
1582 fhPtDecayIso[ibit]->SetYTitle("#it{counts}");
1583 fhPtDecayIso[ibit]->SetXTitle("#it{p}_{T}(GeV/#it{c})");
1584 outputContainer->Add(fhPtDecayIso[ibit]) ;
1586 fhEtaPhiDecayIso[ibit] =
1587 new TH2F(Form("hEtaPhiDecayIso_bit%d",fDecayBits[ibit]),
1588 Form("Number of isolated Pi0 decay particles #eta vs #phi, bit %d, %s",fDecayBits[ibit],parTitle.Data()),
1589 netabins,etamin,etamax,nphibins,phimin,phimax);
1590 fhEtaPhiDecayIso[ibit]->SetXTitle("#eta");
1591 fhEtaPhiDecayIso[ibit]->SetYTitle("#phi");
1592 outputContainer->Add(fhEtaPhiDecayIso[ibit]) ;
1597 for(Int_t imc = 0; imc < fgkNmcTypes; imc++)
1600 fhPtDecayNoIsoMC[ibit][imc] =
1601 new TH1F(Form("hPtDecayNoIso_bit%d_MC%s",fDecayBits[ibit],mcPartName[imc].Data()),
1602 Form("#it{p}_{T} of NOT isolated, decay bit %d, %s, %s",fDecayBits[ibit],mcPartType[imc].Data(),parTitle.Data()),
1603 nptbins,ptmin,ptmax);
1604 fhPtDecayNoIsoMC[ibit][imc]->SetYTitle("#it{counts}");
1605 fhPtDecayNoIsoMC[ibit][imc]->SetXTitle("#it{p}_{T}(GeV/#it{c})");
1606 outputContainer->Add(fhPtDecayNoIsoMC[ibit][imc]) ;
1610 fhPtDecayIsoMC[ibit][imc] =
1611 new TH1F(Form("hPtDecay_bit%d_MC%s",fDecayBits[ibit],mcPartName[imc].Data()),
1612 Form("#it{p}_{T} of isolated %s, decay bit %d, %s",mcPartType[imc].Data(),fDecayBits[ibit],parTitle.Data()),
1613 nptbins,ptmin,ptmax);
1614 fhPtDecayIsoMC[ibit][imc]->SetYTitle("#it{counts}");
1615 fhPtDecayIsoMC[ibit][imc]->SetXTitle("#it{p}_{T}(GeV/#it{c})");
1616 outputContainer->Add(fhPtDecayIsoMC[ibit][imc]) ;
1618 }// MC particle loop
1625 TString isoName [] = {"NoIso",""};
1626 TString isoTitle[] = {"Not isolated" ,"isolated"};
1628 fhEIso = new TH1F("hE",
1629 Form("Number of isolated particles vs E, %s",parTitle.Data()),
1630 nptbins,ptmin,ptmax);
1631 fhEIso->SetYTitle("d#it{N} / d#it{E}");
1632 fhEIso->SetXTitle("#it{E} (GeV/#it{c})");
1633 outputContainer->Add(fhEIso) ;
1635 fhPtIso = new TH1F("hPt",
1636 Form("Number of isolated particles vs #it{p}_{T}, %s",parTitle.Data()),
1637 nptbins,ptmin,ptmax);
1638 fhPtIso->SetYTitle("d#it{N} / #it{p}_{T}");
1639 fhPtIso->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1640 outputContainer->Add(fhPtIso) ;
1642 fhPhiIso = new TH2F("hPhi",
1643 Form("Number of isolated particles vs #phi, %s",parTitle.Data()),
1644 nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1645 fhPhiIso->SetYTitle("#phi");
1646 fhPhiIso->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1647 outputContainer->Add(fhPhiIso) ;
1649 fhEtaIso = new TH2F("hEta",
1650 Form("Number of isolated particles vs #eta, %s",parTitle.Data()),
1651 nptbins,ptmin,ptmax,netabins,etamin,etamax);
1652 fhEtaIso->SetYTitle("#eta");
1653 fhEtaIso->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1654 outputContainer->Add(fhEtaIso) ;
1656 fhEtaPhiIso = new TH2F("hEtaPhiIso",
1657 Form("Number of isolated particles #eta vs #phi, %s",parTitle.Data()),
1658 netabins,etamin,etamax,nphibins,phimin,phimax);
1659 fhEtaPhiIso->SetXTitle("#eta");
1660 fhEtaPhiIso->SetYTitle("#phi");
1661 outputContainer->Add(fhEtaPhiIso) ;
1663 if(fFillHighMultHistograms)
1665 fhPtCentralityIso = new TH2F("hPtCentrality",
1666 Form("centrality vs #it{p}_{T} for isolated particles, %s",parTitle.Data()),
1667 nptbins,ptmin,ptmax, 100,0,100);
1668 fhPtCentralityIso->SetYTitle("centrality");
1669 fhPtCentralityIso->SetXTitle("#it{p}_{T}(GeV/#it{c})");
1670 outputContainer->Add(fhPtCentralityIso) ;
1672 fhPtEventPlaneIso = new TH2F("hPtEventPlane",
1673 Form("event plane angle vs #it{p}_{T} for isolated particles, %s",parTitle.Data()),
1674 nptbins,ptmin,ptmax, 100,0,TMath::Pi());
1675 fhPtEventPlaneIso->SetYTitle("Event plane angle (rad)");
1676 fhPtEventPlaneIso->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1677 outputContainer->Add(fhPtEventPlaneIso) ;
1680 if(fFillNLMHistograms)
1682 fhPtNLocMaxIso = new TH2F("hPtNLocMax",
1683 Form("Number of isolated particles vs #it{p}_{T}, %s",parTitle.Data()),
1684 nptbins,ptmin,ptmax,10,0,10);
1685 fhPtNLocMaxIso->SetYTitle("#it{NLM}");
1686 fhPtNLocMaxIso->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1688 fhPtNLocMaxNoIso = new TH2F("hPtNLocMaxNoIso",
1689 Form("Number of not isolated particles vs #it{p}_{T}, %s",parTitle.Data()),
1690 nptbins,ptmin,ptmax,10,0,10);
1691 fhPtNLocMaxNoIso->SetYTitle("#it{NLM}");
1692 fhPtNLocMaxNoIso->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1693 outputContainer->Add(fhPtNLocMaxNoIso) ;
1696 fhConePtLead = new TH2F("hConePtLead",
1697 Form("Track or Cluster leading #it{p}_{T} in isolation cone for #it{R} = %2.2f",r),
1698 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1699 fhConePtLead->SetYTitle("#it{p}_{T, leading} (GeV/#it{c})");
1700 fhConePtLead->SetXTitle("#it{p}_{T, trigger} (GeV/#it{c})");
1701 outputContainer->Add(fhConePtLead) ;
1703 fhConeSumPt = new TH2F("hConePtSum",
1704 Form("Track and Cluster #Sigma #it{p}_{T} in isolation cone for #it{R} = %2.2f",r),
1705 nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
1706 fhConeSumPt->SetYTitle("#Sigma #it{p}_{T}");
1707 fhConeSumPt->SetXTitle("#it{p}_{T, trigger} (GeV/#it{c})");
1708 outputContainer->Add(fhConeSumPt) ;
1710 fhConeSumPtTrigEtaPhi = new TH2F("hConePtSumTrigEtaPhi",
1711 Form("Trigger #eta vs #phi, #Sigma #it{p}_{T} in isolation cone for #it{R} = %2.2f",r),
1712 netabins,etamin,etamax,nphibins,phimin,phimax);
1713 fhConeSumPtTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
1714 fhConeSumPtTrigEtaPhi->SetXTitle("#eta_{trigger}");
1715 fhConeSumPtTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
1716 outputContainer->Add(fhConeSumPtTrigEtaPhi) ;
1718 fhPtInCone = new TH2F("hPtInCone",
1719 Form("#it{p}_{T} of clusters and tracks in isolation cone for #it{R} = %2.2f",r),
1720 nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
1721 fhPtInCone->SetYTitle("#it{p}_{T in cone} (GeV/#it{c})");
1722 fhPtInCone->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1723 outputContainer->Add(fhPtInCone) ;
1725 if(fFillBackgroundBinHistograms)
1727 fhPtLeadConeBinLambda0 = new TH2F*[fNBkgBin];
1728 fhSumPtConeBinLambda0 = new TH2F*[fNBkgBin];
1732 fhPtLeadConeBinLambda0MC = new TH2F*[fNBkgBin*fgkNmcTypes];
1733 fhSumPtConeBinLambda0MC = new TH2F*[fNBkgBin*fgkNmcTypes];
1736 for(Int_t ibin = 0; ibin < fNBkgBin; ibin++)
1738 fhPtLeadConeBinLambda0[ibin] = new TH2F
1739 (Form("hPtLeadConeLambda0_Bin%d",ibin),
1740 Form("#lambda_{0}, in cone %2.2f<#it{p}_{T}^{leading}<%2.2f (GeV/#it{c}), %s",
1741 fBkgBinLimit[ibin],fBkgBinLimit[ibin+1], parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1742 fhPtLeadConeBinLambda0[ibin]->SetYTitle("#lambda_{0}^{2}");
1743 fhPtLeadConeBinLambda0[ibin]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1744 outputContainer->Add(fhPtLeadConeBinLambda0[ibin]) ;
1746 fhSumPtConeBinLambda0[ibin] = new TH2F
1747 (Form("hSumPtConeLambda0_Bin%d",ibin),
1748 Form("#lambda_{0}, in cone %2.2f <#Sigma #it{p}_{T}< %2.2f (GeV/#it{c}), %s",
1749 fBkgBinLimit[ibin],fBkgBinLimit[ibin+1], parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1750 fhSumPtConeBinLambda0[ibin]->SetYTitle("#lambda_{0}^{2}");
1751 fhSumPtConeBinLambda0[ibin]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1752 outputContainer->Add(fhSumPtConeBinLambda0[ibin]) ;
1756 for(Int_t imc = 0; imc < fgkNmcTypes; imc++)
1758 Int_t binmc = ibin+imc*fNBkgBin;
1759 fhPtLeadConeBinLambda0MC[binmc] = new TH2F
1760 (Form("hPtLeadConeLambda0_Bin%d_MC%s",ibin, mcPartName[imc].Data()),
1761 Form("#lambda_{0}, in cone %2.2f<#it{p}_{T}^{leading}<%2.2f (GeV/#it{c}), MC %s, %s",
1762 fBkgBinLimit[ibin],fBkgBinLimit[ibin+1], mcPartType[imc].Data(), parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1763 fhPtLeadConeBinLambda0MC[binmc]->SetYTitle("#lambda_{0}^{2}");
1764 fhPtLeadConeBinLambda0MC[binmc]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1765 outputContainer->Add(fhPtLeadConeBinLambda0MC[binmc]) ;
1767 fhSumPtConeBinLambda0MC[binmc] = new TH2F
1768 (Form("hSumPtConeLambda0_Bin%d_MC%s",ibin,mcPartName[imc].Data()),
1769 Form("#lambda_{0}, in cone %2.2f <#Sigma #it{p}_{T}< %2.2f (GeV/#it{c}), MC %s, %s",
1770 fBkgBinLimit[ibin],fBkgBinLimit[ibin+1], mcPartType[imc].Data(), parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1771 fhSumPtConeBinLambda0MC[binmc]->SetYTitle("#lambda_{0}^{2}");
1772 fhSumPtConeBinLambda0MC[binmc]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1773 outputContainer->Add(fhSumPtConeBinLambda0MC[binmc]) ;
1774 } // MC particle loop
1778 } // bkg cone pt bin histograms
1780 if(fFillHighMultHistograms)
1782 fhPtInConeCent = new TH2F("hPtInConeCent",
1783 Form("#it{p}_{T} in isolation cone for #it{R} = %2.2f",r),
1784 100,0,100,nptinconebins,ptinconemin,ptinconemax);
1785 fhPtInConeCent->SetYTitle("#it{p}_{T in cone} (GeV/#it{c})");
1786 fhPtInConeCent->SetXTitle("centrality");
1787 outputContainer->Add(fhPtInConeCent) ;
1790 // Cluster only histograms
1791 if(GetIsolationCut()->GetParticleTypeInCone()!=AliIsolationCut::kOnlyCharged)
1793 fhConeSumPtCluster = new TH2F("hConePtSumCluster",
1794 Form("Cluster #Sigma #it{p}_{T} in isolation cone for #it{R} = %2.2f",r),
1795 nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
1796 fhConeSumPtCluster->SetYTitle("#Sigma #it{p}_{T}");
1797 fhConeSumPtCluster->SetXTitle("#it{p}_{T, trigger} (GeV/#it{c})");
1798 outputContainer->Add(fhConeSumPtCluster) ;
1800 fhConePtLeadCluster = new TH2F("hConeLeadPtCluster",
1801 Form("Cluster leading in isolation cone for #it{R} = %2.2f",r),
1802 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1803 fhConePtLeadCluster->SetYTitle("#it{p}_{T, leading} (GeV/#it{c})");
1804 fhConePtLeadCluster->SetXTitle("#it{p}_{T, trigger} (GeV/#it{c})");
1805 outputContainer->Add(fhConePtLeadCluster) ;
1808 if(fFillCellHistograms)
1810 fhConeSumPtCell = new TH2F("hConePtSumCell",
1811 Form("Cell #Sigma #it{p}_{T} in isolation cone for #it{R} = %2.2f",r),
1812 nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
1813 fhConeSumPtCell->SetYTitle("#Sigma #it{p}_{T}");
1814 fhConeSumPtCell->SetXTitle("#it{p}_{T, trigger} (GeV/#it{c})");
1815 outputContainer->Add(fhConeSumPtCell) ;
1818 if(fFillUEBandSubtractHistograms)
1820 fhConeSumPtEtaBandUECluster = new TH2F("hConePtSumEtaBandUECluster",
1821 "#Sigma cluster #it{p}_{T} in UE Eta Band",
1822 nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
1823 fhConeSumPtEtaBandUECluster->SetYTitle("#Sigma #it{p}_{T}");
1824 fhConeSumPtEtaBandUECluster->SetXTitle("#it{p}_{T, trigger} (GeV/#it{c})");
1825 outputContainer->Add(fhConeSumPtEtaBandUECluster) ;
1827 fhConeSumPtPhiBandUECluster = new TH2F("hConePtSumPhiBandUECluster",
1828 "#Sigma cluster #it{p}_{T} UE Phi Band",
1829 nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
1830 fhConeSumPtPhiBandUECluster->SetYTitle("#Sigma #it{p}_{T}");
1831 fhConeSumPtPhiBandUECluster->SetXTitle("#it{p}_{T, trigger} (GeV/#it{c})");
1832 outputContainer->Add(fhConeSumPtPhiBandUECluster) ;
1834 fhConeSumPtEtaBandUEClusterTrigEtaPhi = new TH2F("hConePtSumEtaBandUEClusterTrigEtaPhi",
1835 "Trigger #eta vs #phi, #Sigma cluster #it{p}_{T} in UE Eta Band",
1836 netabins,etamin,etamax,nphibins,phimin,phimax);
1837 fhConeSumPtEtaBandUEClusterTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
1838 fhConeSumPtEtaBandUEClusterTrigEtaPhi->SetXTitle("#eta_{trigger}");
1839 fhConeSumPtEtaBandUEClusterTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
1840 outputContainer->Add(fhConeSumPtEtaBandUEClusterTrigEtaPhi) ;
1842 fhConeSumPtPhiBandUEClusterTrigEtaPhi = new TH2F("hConePtSumPhiBandUEClusterTrigEtaPhi",
1843 "Trigger #eta vs #phi, #Sigma cluster #it{p}_{T} UE Phi Band",
1844 netabins,etamin,etamax,nphibins,phimin,phimax);
1845 fhConeSumPtPhiBandUEClusterTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
1846 fhConeSumPtPhiBandUEClusterTrigEtaPhi->SetXTitle("#eta_{trigger}");
1847 fhConeSumPtPhiBandUEClusterTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
1848 outputContainer->Add(fhConeSumPtPhiBandUEClusterTrigEtaPhi) ;
1849 if(fFillCellHistograms)
1852 fhConeSumPtEtaBandUECell = new TH2F("hConePtSumEtaBandUECell",
1853 "#Sigma cell #it{p}_{T} in UE Eta Band",
1854 nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
1855 fhConeSumPtEtaBandUECell->SetYTitle("#Sigma #it{p}_{T}");
1856 fhConeSumPtEtaBandUECell->SetXTitle("#it{p}_{T, trigger} (GeV/#it{c})");
1857 outputContainer->Add(fhConeSumPtEtaBandUECell) ;
1859 fhConeSumPtPhiBandUECell = new TH2F("hConePtSumPhiBandUECell",
1860 "#Sigma cell #it{p}_{T} UE Phi Band",
1861 nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
1862 fhConeSumPtPhiBandUECell->SetYTitle("#Sigma #it{p}_{T}");
1863 fhConeSumPtPhiBandUECell->SetXTitle("#it{p}_{T, trigger} (GeV/#it{c})");
1864 outputContainer->Add(fhConeSumPtPhiBandUECell) ;
1866 fhConeSumPtEtaBandUECellTrigEtaPhi = new TH2F("hConePtSumEtaBandUECellTrigEtaPhi",
1867 "Trigger #eta vs #phi, #Sigma cell #it{p}_{T} in UE Eta Band",
1868 netabins,etamin,etamax,nphibins,phimin,phimax);
1869 fhConeSumPtEtaBandUECellTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
1870 fhConeSumPtEtaBandUECellTrigEtaPhi->SetXTitle("#eta_{trigger}");
1871 fhConeSumPtEtaBandUECellTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
1872 outputContainer->Add(fhConeSumPtEtaBandUECellTrigEtaPhi) ;
1874 fhConeSumPtPhiBandUECellTrigEtaPhi = new TH2F("hConePtSumPhiBandUECellTrigEtaPhi",
1875 "Trigger #eta vs #phi, #Sigma cell #it{p}_{T} UE Phi Band",
1876 netabins,etamin,etamax,nphibins,phimin,phimax);
1877 fhConeSumPtPhiBandUECellTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
1878 fhConeSumPtPhiBandUECellTrigEtaPhi->SetXTitle("#eta_{trigger}");
1879 fhConeSumPtPhiBandUECellTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
1880 outputContainer->Add(fhConeSumPtPhiBandUECellTrigEtaPhi) ;
1883 fhEtaBandCluster = new TH2F("hEtaBandCluster",
1884 Form("#eta vs #phi of clusters in #eta band isolation cone for #it{R} = %2.2f",r),
1885 netabins,-1,1,nphibins,0,TMath::TwoPi());
1886 fhEtaBandCluster->SetXTitle("#eta");
1887 fhEtaBandCluster->SetYTitle("#phi");
1888 outputContainer->Add(fhEtaBandCluster) ;
1890 fhPhiBandCluster = new TH2F("hPhiBandCluster",
1891 Form("#eta vs #phi of clusters in #phi band isolation cone for #it{R} = %2.2f",r),
1892 netabins,-1,1,nphibins,0,TMath::TwoPi());
1893 fhPhiBandCluster->SetXTitle("#eta");
1894 fhPhiBandCluster->SetYTitle("#phi");
1895 outputContainer->Add(fhPhiBandCluster) ;
1897 fhEtaPhiInConeCluster= new TH2F("hEtaPhiInConeCluster",
1898 Form("#eta vs #phi of clusters in cone for #it{R} = %2.2f",r),
1899 netabins,-1,1,nphibins,0,TMath::TwoPi());
1900 fhEtaPhiInConeCluster->SetXTitle("#eta");
1901 fhEtaPhiInConeCluster->SetYTitle("#phi");
1902 outputContainer->Add(fhEtaPhiInConeCluster) ;
1904 fhEtaPhiCluster= new TH2F("hEtaPhiCluster",
1905 Form("#eta vs #phi of all clusters"),
1906 netabins,-1,1,nphibins,0,TMath::TwoPi());
1907 fhEtaPhiCluster->SetXTitle("#eta");
1908 fhEtaPhiCluster->SetYTitle("#phi");
1909 outputContainer->Add(fhEtaPhiCluster) ;
1913 fhPtClusterInCone = new TH2F("hPtClusterInCone",
1914 Form("#it{p}_{T} of clusters in isolation cone for #it{R} = %2.2f",r),
1915 nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
1916 fhPtClusterInCone->SetYTitle("#it{p}_{T in cone} (GeV/#it{c})");
1917 fhPtClusterInCone->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1918 outputContainer->Add(fhPtClusterInCone) ;
1920 if(fFillCellHistograms)
1922 fhPtCellInCone = new TH2F("hPtCellInCone",
1923 Form("#it{p}_{T} of cells in isolation cone for #it{R} = %2.2f",r),
1924 nptbins,ptmin,ptmax,1000,0,50);
1925 fhPtCellInCone->SetYTitle("#it{p}_{T in cone} (GeV/#it{c})");
1926 fhPtCellInCone->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1927 outputContainer->Add(fhPtCellInCone) ;
1929 fhEtaBandCell = new TH2F("hEtaBandCell",
1930 Form("#col vs #row of cells in #eta band isolation cone for #it{R} = %2.2f",r),
1932 fhEtaBandCell->SetXTitle("#col");
1933 fhEtaBandCell->SetYTitle("#row");
1934 outputContainer->Add(fhEtaBandCell) ;
1936 fhPhiBandCell = new TH2F("hPhiBandCell",
1937 Form("#col vs #row of cells in #phi band isolation cone for #it{R} = %2.2f",r),
1939 fhPhiBandCell->SetXTitle("#col");
1940 fhPhiBandCell->SetYTitle("#row");
1941 outputContainer->Add(fhPhiBandCell) ;
1944 if(fFillUEBandSubtractHistograms)
1946 fhConeSumPtEtaUESubCluster = new TH2F("hConeSumPtEtaUESubCluster",
1947 Form("Clusters #Sigma #it{p}_{T} after bkg subtraction from eta band in the isolation cone for #it{R} = %2.2f",r),
1948 nptbins,ptmin,ptmax,2*nptsumbins,-ptsummax,ptsummax);
1949 fhConeSumPtEtaUESubCluster->SetYTitle("#Sigma #it{p}_{T}");
1950 fhConeSumPtEtaUESubCluster->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1951 outputContainer->Add(fhConeSumPtEtaUESubCluster) ;
1953 fhConeSumPtPhiUESubCluster = new TH2F("hConeSumPtPhiUESubCluster",
1954 Form("Clusters #Sigma #it{p}_{T} after bkg subtraction from phi band in the isolation cone for #it{R} = %2.2f",r),
1955 nptbins,ptmin,ptmax,2*nptsumbins,-ptsummax,ptsummax);
1956 fhConeSumPtPhiUESubCluster->SetYTitle("#Sigma #it{p}_{T}");
1957 fhConeSumPtPhiUESubCluster->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1958 outputContainer->Add(fhConeSumPtPhiUESubCluster) ;
1960 fhConeSumPtEtaUESubClusterTrigEtaPhi = new TH2F("hConeSumPtEtaUESubClusterTrigEtaPhi",
1961 Form("Trigger #eta vs #phi, Clusters #Sigma #it{p}_{T} after bkg subtraction from eta band in the isolation cone for #it{R} = %2.2f",r),
1962 netabins,etamin,etamax,nphibins,phimin,phimax);
1963 fhConeSumPtEtaUESubClusterTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
1964 fhConeSumPtEtaUESubClusterTrigEtaPhi->SetXTitle("#eta_{trigger}");
1965 fhConeSumPtEtaUESubClusterTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
1966 outputContainer->Add(fhConeSumPtEtaUESubClusterTrigEtaPhi) ;
1968 fhConeSumPtPhiUESubClusterTrigEtaPhi = new TH2F("hConeSumPtPhiUESubClusterTrigEtaPhi",
1969 Form("Trigger #eta vs #phi, Clusters #Sigma #it{p}_{T} after bkg subtraction from phi band in the isolation cone for #it{R} = %2.2f",r),
1970 netabins,etamin,etamax,nphibins,phimin,phimax);
1971 fhConeSumPtPhiUESubClusterTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
1972 fhConeSumPtPhiUESubClusterTrigEtaPhi->SetXTitle("#eta_{trigger}");
1973 fhConeSumPtPhiUESubClusterTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
1974 outputContainer->Add(fhConeSumPtPhiUESubClusterTrigEtaPhi) ;
1976 if(fFillCellHistograms)
1978 fhConeSumPtEtaUESubCell = new TH2F("hConeSumPtEtaUESubCell",
1979 Form("Cells #Sigma #it{p}_{T} after bkg subtraction from eta band in the isolation cone for #it{R} = %2.2f",r),
1980 nptbins,ptmin,ptmax,2*nptsumbins,-ptsummax,ptsummax);
1981 fhConeSumPtEtaUESubCell->SetYTitle("#Sigma #it{p}_{T}");
1982 fhConeSumPtEtaUESubCell->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1983 outputContainer->Add(fhConeSumPtEtaUESubCell) ;
1985 fhConeSumPtPhiUESubCell = new TH2F("hConeSumPtPhiUESubCell",
1986 Form("Cells #Sigma #it{p}_{T} after bkg subtraction from phi band in the isolation cone for #it{R} = %2.2f",r),
1987 nptbins,ptmin,ptmax,2*nptsumbins,-ptsummax,ptsummax);
1988 fhConeSumPtPhiUESubCell->SetYTitle("#Sigma #it{p}_{T}");
1989 fhConeSumPtPhiUESubCell->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1990 outputContainer->Add(fhConeSumPtPhiUESubCell) ;
1992 fhConeSumPtEtaUESubCellTrigEtaPhi = new TH2F("hConeSumPtEtaUESubCellTrigEtaPhi",
1993 Form("Trigger #eta vs #phi, Cells #Sigma #it{p}_{T} after bkg subtraction from eta band in the isolation cone for #it{R} = %2.2f",r),
1994 netabins,etamin,etamax,nphibins,phimin,phimax);
1995 fhConeSumPtEtaUESubCellTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
1996 fhConeSumPtEtaUESubCellTrigEtaPhi->SetXTitle("#eta_{trigger}");
1997 fhConeSumPtEtaUESubCellTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
1998 outputContainer->Add(fhConeSumPtEtaUESubCellTrigEtaPhi) ;
2000 fhConeSumPtPhiUESubCellTrigEtaPhi = new TH2F("hConeSumPtPhiUESubCellTrigEtaPhi",
2001 Form("Trigger #eta vs #phi, Cells #Sigma #it{p}_{T} after bkg subtraction from phi band in the isolation cone for #it{R} = %2.2f",r),
2002 netabins,etamin,etamax,nphibins,phimin,phimax);
2003 fhConeSumPtPhiUESubCellTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
2004 fhConeSumPtPhiUESubCellTrigEtaPhi->SetXTitle("#eta_{trigger}");
2005 fhConeSumPtPhiUESubCellTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2006 outputContainer->Add(fhConeSumPtPhiUESubCellTrigEtaPhi) ;
2009 fhFractionClusterOutConeEta = new TH2F("hFractionClusterOutConeEta",
2010 Form("Fraction of the isolation cone #it{R} = %2.2f, out of clusters #eta acceptance",r),
2011 nptbins,ptmin,ptmax,100,0,1);
2012 fhFractionClusterOutConeEta->SetYTitle("#it{fraction}");
2013 fhFractionClusterOutConeEta->SetXTitle("#it{p}_{T,trigger} (GeV/#it{c})");
2014 outputContainer->Add(fhFractionClusterOutConeEta) ;
2016 fhFractionClusterOutConeEtaTrigEtaPhi = new TH2F("hFractionClusterOutConeEtaTrigEtaPhi",
2017 Form("Fraction of the isolation cone #it{R} = %2.2f, out of clusters #eta acceptance, in trigger #eta-#phi ",r),
2018 netabins,etamin,etamax,nphibins,phimin,phimax);
2019 fhFractionClusterOutConeEtaTrigEtaPhi->SetZTitle("#it{fraction}");
2020 fhFractionClusterOutConeEtaTrigEtaPhi->SetXTitle("#eta_{trigger}");
2021 fhFractionClusterOutConeEtaTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2022 outputContainer->Add(fhFractionClusterOutConeEtaTrigEtaPhi) ;
2024 fhFractionClusterOutConePhi = new TH2F("hFractionClusterOutConePhi",
2025 Form("Fraction of the isolation cone #it{R} = %2.2f, out of clusters #phi acceptance",r),
2026 nptbins,ptmin,ptmax,100,0,1);
2027 fhFractionClusterOutConePhi->SetYTitle("#it{fraction}");
2028 fhFractionClusterOutConePhi->SetXTitle("#it{p}_{T,trigger} (GeV/#it{c})");
2029 outputContainer->Add(fhFractionClusterOutConePhi) ;
2031 fhFractionClusterOutConePhiTrigEtaPhi = new TH2F("hFractionClusterOutConePhiTrigEtaPhi",
2032 Form("Fraction of the isolation cone #it{R} = %2.2f, out of clusters #phi acceptance, in trigger #eta-#phi ",r),
2033 netabins,etamin,etamax,nphibins,phimin,phimax);
2034 fhFractionClusterOutConePhiTrigEtaPhi->SetZTitle("#it{fraction}");
2035 fhFractionClusterOutConePhiTrigEtaPhi->SetXTitle("#eta_{trigger}");
2036 fhFractionClusterOutConePhiTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2037 outputContainer->Add(fhFractionClusterOutConePhiTrigEtaPhi) ;
2039 fhConeSumPtSubvsConeSumPtTotPhiCluster = new TH2F("hConeSumPtSubvsConeSumPtTotPhiCluster",
2040 Form("#Sigma #it{p}_{T} in cone after bkg sub from #phi band vs #Sigma #it{p}_{T} in cone before bkg sub, R=%2.2f",r),
2041 nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2042 fhConeSumPtSubvsConeSumPtTotPhiCluster->SetXTitle("#Sigma #it{p}_{T, tot} (GeV/#it{c})");
2043 fhConeSumPtSubvsConeSumPtTotPhiCluster->SetYTitle("#Sigma #it{p}_{T, sub} (GeV/#it{c})");
2044 outputContainer->Add(fhConeSumPtSubvsConeSumPtTotPhiCluster);
2046 fhConeSumPtSubNormvsConeSumPtTotPhiCluster = new TH2F("hConeSumPtSubNormvsConeSumPtTotPhiCluster",
2047 Form("#Sigma #it{p}_{T, norm} in cone after bkg sub from #phi band vs #Sigma #it{p}_{T} in cone before bkg sub, R=%2.2f",r),
2048 nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2049 fhConeSumPtSubNormvsConeSumPtTotPhiCluster->SetXTitle("#Sigma #it{p}_{T, tot} (GeV/#it{c})");
2050 fhConeSumPtSubNormvsConeSumPtTotPhiCluster->SetYTitle("#Sigma #it{p}_{T, sub norm} (GeV/#it{c})");
2051 outputContainer->Add(fhConeSumPtSubNormvsConeSumPtTotPhiCluster);
2053 fhConeSumPtSubvsConeSumPtTotEtaCluster = new TH2F("hConeSumPtSubvsConeSumPtTotEtaCluster",
2054 Form("#Sigma #it{p}_{T} in cone after bkg sub from #eta band vs #Sigma #it{p}_{T} in cone before bkg sub, R=%2.2f",r),
2055 nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2056 fhConeSumPtSubvsConeSumPtTotEtaCluster->SetXTitle("#Sigma #it{p}_{T, tot} (GeV/#it{c})");
2057 fhConeSumPtSubvsConeSumPtTotEtaCluster->SetYTitle("#Sigma #it{p}_{T, sub} (GeV/#it{c})");
2058 outputContainer->Add(fhConeSumPtSubvsConeSumPtTotEtaCluster);
2060 fhConeSumPtSubNormvsConeSumPtTotEtaCluster = new TH2F("hConeSumPtSubNormvsConeSumPtTotEtaCluster",
2061 Form("#Sigma #it{p}_{T, norm} in cone after bkg sub from #eta band vs #Sigma #it{p}_{T} in cone before bkg sub, R=%2.2f",r),
2062 nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2063 fhConeSumPtSubNormvsConeSumPtTotEtaCluster->SetXTitle("#Sigma #it{p}_{T, tot} (GeV/#it{c})");
2064 fhConeSumPtSubNormvsConeSumPtTotEtaCluster->SetYTitle("#Sigma #it{p}_{T, sub norm} (GeV/#it{c})");
2065 outputContainer->Add(fhConeSumPtSubNormvsConeSumPtTotEtaCluster);
2067 fhConeSumPtVSUEClusterEtaBand = new TH2F("hConeSumPtVSUEClusterEtaBand",
2068 Form("#Sigma #it{p}_{T} in cone versus #Sigma #it{p}_{T} in eta band for cluster (before normalization), R=%2.2f",r),
2069 nptsumbins,ptsummin,ptsummax,2*nptsumbins,ptsummin,2*ptsummax);
2070 fhConeSumPtVSUEClusterEtaBand->SetXTitle("#Sigma #it{p}_{T} cone (GeV/#it{c})");
2071 fhConeSumPtVSUEClusterEtaBand->SetYTitle("#Sigma #it{p}_{T} UE (GeV/#it{c})");
2072 outputContainer->Add(fhConeSumPtVSUEClusterEtaBand);
2074 fhConeSumPtVSUEClusterPhiBand = new TH2F("hConeSumPtVSUEClusterPhiBand",
2075 Form("#Sigma #it{p}_{T} in cone versus #Sigma #it{p}_{T} in phi band for cluster (before normalization), R=%2.2f",r),
2076 nptsumbins,ptsummin,ptsummax,8*nptsumbins,ptsummin,8*ptsummax);
2077 fhConeSumPtVSUEClusterPhiBand->SetXTitle("#Sigma #it{p}_{T} cone (GeV/#it{c})");
2078 fhConeSumPtVSUEClusterPhiBand->SetYTitle("#Sigma #it{p}_{T} UE (GeV/#it{c})");
2079 outputContainer->Add(fhConeSumPtVSUEClusterPhiBand);
2081 if(fFillCellHistograms)
2083 fhFractionCellOutConeEta = new TH2F("hFractionCellOutConeEta",
2084 Form("Fraction of the isolation cone #it{R} = %2.2f, out of cells #eta acceptance",r),
2085 nptbins,ptmin,ptmax,100,0,1);
2086 fhFractionCellOutConeEta->SetYTitle("#it{fraction}");
2087 fhFractionCellOutConeEta->SetXTitle("#it{p}_{T,trigger} (GeV/#it{c})");
2088 outputContainer->Add(fhFractionCellOutConeEta) ;
2090 fhFractionCellOutConeEtaTrigEtaPhi = new TH2F("hFractionCellOutConeEtaTrigEtaPhi",
2091 Form("Fraction of the isolation cone #it{R} = %2.2f, out of cells #eta acceptance, in trigger #eta-#phi ",r),
2092 netabins,etamin,etamax,nphibins,phimin,phimax);
2093 fhFractionCellOutConeEtaTrigEtaPhi->SetZTitle("#it{fraction}");
2094 fhFractionCellOutConeEtaTrigEtaPhi->SetXTitle("#eta_{trigger}");
2095 fhFractionCellOutConeEtaTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2096 outputContainer->Add(fhFractionCellOutConeEtaTrigEtaPhi) ;
2098 fhFractionCellOutConePhi = new TH2F("hFractionCellOutConePhi",
2099 Form("Fraction of the isolation cone #it{R} = %2.2f, out of cells #phi acceptance",r),
2100 nptbins,ptmin,ptmax,100,0,1);
2101 fhFractionCellOutConePhi->SetYTitle("#it{fraction}");
2102 fhFractionCellOutConePhi->SetXTitle("#it{p}_{T,trigger} (GeV/#it{c})");
2103 outputContainer->Add(fhFractionCellOutConePhi) ;
2105 fhFractionCellOutConePhiTrigEtaPhi = new TH2F("hFractionCellOutConePhiTrigEtaPhi",
2106 Form("Fraction of the isolation cone #it{R} = %2.2f, out of cells #phi acceptance, in trigger #eta-#phi ",r),
2107 netabins,etamin,etamax,nphibins,phimin,phimax);
2108 fhFractionCellOutConePhiTrigEtaPhi->SetZTitle("#it{fraction}");
2109 fhFractionCellOutConePhiTrigEtaPhi->SetXTitle("#eta_{trigger}");
2110 fhFractionCellOutConePhiTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2111 outputContainer->Add(fhFractionCellOutConePhiTrigEtaPhi) ;
2114 fhConeSumPtSubvsConeSumPtTotPhiCell = new TH2F("hConeSumPtSubvsConeSumPtTotPhiCell",
2115 Form("#Sigma #it{p}_{T} in cone after bkg sub from #phi band vs #Sigma #it{p}_{T} in cone before bkg sub, R=%2.2f",r),
2116 nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2117 fhConeSumPtSubvsConeSumPtTotPhiCell->SetXTitle("#Sigma #it{p}_{T, tot} (GeV/#it{c})");
2118 fhConeSumPtSubvsConeSumPtTotPhiCell->SetYTitle("#Sigma #it{p}_{T, sub} (GeV/#it{c})");
2119 outputContainer->Add(fhConeSumPtSubvsConeSumPtTotPhiCell);
2121 fhConeSumPtSubNormvsConeSumPtTotPhiCell = new TH2F("hConeSumPtSubNormvsConeSumPtTotPhiCell",
2122 Form("#Sigma #it{p}_{T, norm} in cone after bkg sub from #phi band vs #Sigma #it{p}_{T} in cone before bkg sub, R=%2.2f",r),
2123 nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2124 fhConeSumPtSubNormvsConeSumPtTotPhiCell->SetXTitle("#Sigma #it{p}_{T, tot} (GeV/#it{c})");
2125 fhConeSumPtSubNormvsConeSumPtTotPhiCell->SetYTitle("#Sigma #it{p}_{T, sub norm} (GeV/#it{c})");
2126 outputContainer->Add(fhConeSumPtSubNormvsConeSumPtTotPhiCell);
2128 fhConeSumPtSubvsConeSumPtTotEtaCell = new TH2F("hConeSumPtSubvsConeSumPtTotEtaCell",
2129 Form("#Sigma #it{p}_{T} in cone after bkg sub from #eta band vs #Sigma #it{p}_{T} in cone before bkg sub, R=%2.2f",r),
2130 nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2131 fhConeSumPtSubvsConeSumPtTotEtaCell->SetXTitle("#Sigma #it{p}_{T, tot} (GeV/#it{c})");
2132 fhConeSumPtSubvsConeSumPtTotEtaCell->SetYTitle("#Sigma #it{p}_{T, sub} (GeV/#it{c})");
2133 outputContainer->Add(fhConeSumPtSubvsConeSumPtTotEtaCell);
2135 fhConeSumPtSubNormvsConeSumPtTotEtaCell = new TH2F("hConeSumPtSubNormvsConeSumPtTotEtaCell",
2136 Form("#Sigma #it{p}_{T, norm} in cone after bkg sub from #eta band vs #Sigma #it{p}_{T} in cone before bkg sub, R=%2.2f",r),
2137 nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2138 fhConeSumPtSubNormvsConeSumPtTotEtaCell->SetXTitle("#Sigma #it{p}_{T, tot} (GeV/#it{c})");
2139 fhConeSumPtSubNormvsConeSumPtTotEtaCell->SetYTitle("#Sigma #it{p}_{T, sub norm} (GeV/#it{c})");
2140 outputContainer->Add(fhConeSumPtSubNormvsConeSumPtTotEtaCell);
2145 // Track only histograms
2146 if(GetIsolationCut()->GetParticleTypeInCone()!=AliIsolationCut::kOnlyNeutral)
2148 fhConeSumPtTrack = new TH2F("hConePtSumTrack",
2149 Form("Track #Sigma #it{p}_{T} in isolation cone for #it{R} = %2.2f",r),
2150 nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
2151 fhConeSumPtTrack->SetYTitle("#Sigma #it{p}_{T}");
2152 fhConeSumPtTrack->SetXTitle("#it{p}_{T, trigger} (GeV/#it{c})");
2153 outputContainer->Add(fhConeSumPtTrack) ;
2155 fhConePtLeadTrack = new TH2F("hConeLeadPtTrack",
2156 Form("Track leading in isolation cone for #it{R} = %2.2f",r),
2157 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2158 fhConePtLeadTrack->SetYTitle("#it{p}_{T, leading} (GeV/#it{c})");
2159 fhConePtLeadTrack->SetXTitle("#it{p}_{T, trigger} (GeV/#it{c})");
2160 outputContainer->Add(fhConePtLeadTrack) ;
2162 fhPtTrackInCone = new TH2F("hPtTrackInCone",
2163 Form("#it{p}_{T} of tracks in isolation cone for #it{R} = %2.2f",r),
2164 nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
2165 fhPtTrackInCone->SetYTitle("#it{p}_{T in cone} (GeV/#it{c})");
2166 fhPtTrackInCone->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2167 outputContainer->Add(fhPtTrackInCone) ;
2170 if(fFillUEBandSubtractHistograms)
2172 fhConeSumPtEtaBandUETrack = new TH2F("hConePtSumEtaBandUETrack",
2173 "#Sigma track #it{p}_{T} in UE Eta Band",
2174 nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
2175 fhConeSumPtEtaBandUETrack->SetYTitle("#Sigma #it{p}_{T}");
2176 fhConeSumPtEtaBandUETrack->SetXTitle("#it{p}_{T, trigger} (GeV/#it{c})");
2177 outputContainer->Add(fhConeSumPtEtaBandUETrack) ;
2179 fhConeSumPtPhiBandUETrack = new TH2F("hConePtSumPhiBandUETrack",
2180 "#Sigma track #it{p}_{T} in UE Phi Band",
2181 nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax*8);
2182 fhConeSumPtPhiBandUETrack->SetYTitle("#Sigma #it{p}_{T}");
2183 fhConeSumPtPhiBandUETrack->SetXTitle("#it{p}_{T, trigger} (GeV/#it{c})");
2184 outputContainer->Add(fhConeSumPtPhiBandUETrack) ;
2187 fhConeSumPtEtaBandUETrackTrigEtaPhi = new TH2F("hConePtSumEtaBandUETrackTrigEtaPhi",
2188 "Trigger #eta vs #phi, #Sigma track #it{p}_{T} in UE Eta Band",
2189 netabins,etamin,etamax,nphibins,phimin,phimax);
2190 fhConeSumPtEtaBandUETrackTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
2191 fhConeSumPtEtaBandUETrackTrigEtaPhi->SetXTitle("#eta_{trigger}");
2192 fhConeSumPtEtaBandUETrackTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2193 outputContainer->Add(fhConeSumPtEtaBandUETrackTrigEtaPhi) ;
2195 fhConeSumPtPhiBandUETrackTrigEtaPhi = new TH2F("hConePtSumPhiBandUETrackTrigEtaPhi",
2196 "Trigger #eta vs #phi, #Sigma track #it{p}_{T} in UE Phi Band",
2197 netabins,etamin,etamax,nphibins,phimin,phimax);
2198 fhConeSumPtPhiBandUETrackTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
2199 fhConeSumPtPhiBandUETrackTrigEtaPhi->SetXTitle("#eta_{trigger}");
2200 fhConeSumPtPhiBandUETrackTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2201 outputContainer->Add(fhConeSumPtPhiBandUETrackTrigEtaPhi) ;
2203 fhEtaBandTrack = new TH2F("hEtaBandTrack",
2204 Form("#eta vs #phi of tracks in #eta band isolation cone for #it{R} = %2.2f",r),
2205 netabins,-1,1,nphibins,0,TMath::TwoPi());
2206 fhEtaBandTrack->SetXTitle("#eta");
2207 fhEtaBandTrack->SetYTitle("#phi");
2208 outputContainer->Add(fhEtaBandTrack) ;
2210 fhPhiBandTrack = new TH2F("hPhiBandTrack",
2211 Form("#eta vs #phi of tracks in #phi band isolation cone for #it{R} = %2.2f",r),
2212 netabins,-1,1,nphibins,0,TMath::TwoPi());
2213 fhPhiBandTrack->SetXTitle("#eta");
2214 fhPhiBandTrack->SetYTitle("#phi");
2215 outputContainer->Add(fhPhiBandTrack) ;
2217 fhConeSumPtEtaUESubTrack = new TH2F("hConeSumPtEtaUESubTrack",
2218 Form("Tracks #Sigma #it{p}_{T} after bkg subtraction from eta band in the isolation cone for #it{R} = %2.2f",r),
2219 nptbins,ptmin,ptmax,2*nptsumbins,-ptsummax,ptsummax);
2220 fhConeSumPtEtaUESubTrack->SetYTitle("#Sigma #it{p}_{T}");
2221 fhConeSumPtEtaUESubTrack->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2222 outputContainer->Add(fhConeSumPtEtaUESubTrack) ;
2224 fhConeSumPtPhiUESubTrack = new TH2F("hConeSumPtPhiUESubTrack",
2225 Form("Tracks #Sigma #it{p}_{T} after bkg subtraction from phi band in the isolation cone for #it{R} = %2.2f",r),
2226 nptbins,ptmin,ptmax,2*nptsumbins,-ptsummax,ptsummax);
2227 fhConeSumPtPhiUESubTrack->SetYTitle("#Sigma #it{p}_{T}");
2228 fhConeSumPtPhiUESubTrack->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2229 outputContainer->Add(fhConeSumPtPhiUESubTrack) ;
2231 fhConeSumPtEtaUESubTrackTrigEtaPhi = new TH2F("hConeSumPtEtaUESubTrackTrigEtaPhi",
2232 Form("Trigger #eta vs #phi, Tracks #Sigma #it{p}_{T} after bkg subtraction from eta band in the isolation cone for #it{R} = %2.2f",r),
2233 netabins,etamin,etamax,nphibins,phimin,phimax);
2234 fhConeSumPtEtaUESubTrackTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
2235 fhConeSumPtEtaUESubTrackTrigEtaPhi->SetXTitle("#eta_{trigger}");
2236 fhConeSumPtEtaUESubTrackTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2237 outputContainer->Add(fhConeSumPtEtaUESubTrackTrigEtaPhi) ;
2239 fhConeSumPtPhiUESubTrackTrigEtaPhi = new TH2F("hConeSumPtPhiUESubTrackTrigEtaPhi",
2240 Form("Trigger #eta vs #phi, Tracks #Sigma #it{p}_{T} after bkg subtraction from phi band in the isolation cone for #it{R} = %2.2f",r),
2241 netabins,etamin,etamax,nphibins,phimin,phimax);
2242 fhConeSumPtPhiUESubTrackTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
2243 fhConeSumPtPhiUESubTrackTrigEtaPhi->SetXTitle("#eta_{trigger}");
2244 fhConeSumPtPhiUESubTrackTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2245 outputContainer->Add(fhConeSumPtPhiUESubTrackTrigEtaPhi) ;
2247 fhFractionTrackOutConeEta = new TH2F("hFractionTrackOutConeEta",
2248 Form("Fraction of the isolation cone #it{R} = %2.2f, out of tracks #eta acceptance",r),
2249 nptbins,ptmin,ptmax,100,0,1);
2250 fhFractionTrackOutConeEta->SetYTitle("#it{fraction}");
2251 fhFractionTrackOutConeEta->SetXTitle("#it{p}_{T,trigger} (GeV/#it{c})");
2252 outputContainer->Add(fhFractionTrackOutConeEta) ;
2254 fhFractionTrackOutConeEtaTrigEtaPhi = new TH2F("hFractionTrackOutConeEtaTrigEtaPhi",
2255 Form("Fraction of the isolation cone #it{R} = %2.2f, out of tracks #eta acceptance, in trigger #eta-#phi ",r),
2256 netabins,etamin,etamax,nphibins,phimin,phimax);
2257 fhFractionTrackOutConeEtaTrigEtaPhi->SetZTitle("#it{fraction}");
2258 fhFractionTrackOutConeEtaTrigEtaPhi->SetXTitle("#eta_{trigger}");
2259 fhFractionTrackOutConeEtaTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2260 outputContainer->Add(fhFractionTrackOutConeEtaTrigEtaPhi) ;
2262 fhConeSumPtSubvsConeSumPtTotPhiTrack = new TH2F("hConeSumPtSubvsConeSumPtTotPhiTrack",
2263 Form("#Sigma #it{p}_{T} in cone after bkg sub from #phi band vs #Sigma #it{p}_{T} in cone before bkg sub, R=%2.2f",r),
2264 nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2265 fhConeSumPtSubvsConeSumPtTotPhiTrack->SetXTitle("#Sigma #it{p}_{T, tot} (GeV/#it{c})");
2266 fhConeSumPtSubvsConeSumPtTotPhiTrack->SetYTitle("#Sigma #it{p}_{T, sub} (GeV/#it{c})");
2267 outputContainer->Add(fhConeSumPtSubvsConeSumPtTotPhiTrack);
2269 fhConeSumPtSubNormvsConeSumPtTotPhiTrack = new TH2F("hConeSumPtSubNormvsConeSumPtTotPhiTrack",
2270 Form("#Sigma #it{p}_{T, norm} in cone after bkg sub from #phi band vs #Sigma #it{p}_{T} in cone before bkg sub, R=%2.2f",r),
2271 nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2272 fhConeSumPtSubNormvsConeSumPtTotPhiTrack->SetXTitle("#Sigma #it{p}_{T, tot} (GeV/#it{c})");
2273 fhConeSumPtSubNormvsConeSumPtTotPhiTrack->SetYTitle("#Sigma #it{p}_{T, sub norm} (GeV/#it{c})");
2274 outputContainer->Add(fhConeSumPtSubNormvsConeSumPtTotPhiTrack);
2276 fhConeSumPtSubvsConeSumPtTotEtaTrack = new TH2F("hConeSumPtSubvsConeSumPtTotEtaTrack",
2277 Form("#Sigma #it{p}_{T} in cone after bkg sub from #eta band vs #Sigma #it{p}_{T} in cone before bkg sub, R=%2.2f",r),
2278 nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2279 fhConeSumPtSubvsConeSumPtTotEtaTrack->SetXTitle("#Sigma #it{p}_{T, tot} (GeV/#it{c})");
2280 fhConeSumPtSubvsConeSumPtTotEtaTrack->SetYTitle("#Sigma #it{p}_{T, sub} (GeV/#it{c})");
2281 outputContainer->Add(fhConeSumPtSubvsConeSumPtTotEtaTrack);
2283 fhConeSumPtSubNormvsConeSumPtTotEtaTrack = new TH2F("hConeSumPtSubNormvsConeSumPtTotEtaTrack",
2284 Form("#Sigma #it{p}_{T, norm} in cone after bkg sub from #eta band vs #Sigma #it{p}_{T} in cone before bkg sub, R=%2.2f",r),
2285 nptsumbins,ptsummin,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2286 fhConeSumPtSubNormvsConeSumPtTotEtaTrack->SetXTitle("#Sigma #it{p}_{T, tot} (GeV/#it{c})");
2287 fhConeSumPtSubNormvsConeSumPtTotEtaTrack->SetYTitle("#Sigma #it{p}_{T, sub norm} (GeV/#it{c})");
2288 outputContainer->Add(fhConeSumPtSubNormvsConeSumPtTotEtaTrack);
2291 // UE in perpendicular cone
2292 fhPerpConeSumPt = new TH2F("hPerpConePtSum",
2293 Form("#Sigma #it{p}_{T} in isolation cone at #pm 45 degree phi from trigger particle, #it{R} = %2.2f",r),
2294 nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
2295 fhPerpConeSumPt->SetYTitle("#Sigma #it{p}_{T}");
2296 fhPerpConeSumPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2297 outputContainer->Add(fhPerpConeSumPt) ;
2299 fhPtInPerpCone = new TH2F("hPtInPerpCone",
2300 Form("#it{p}_{T} in isolation cone at #pm 45 degree phi from trigger particle, #it{R} = %2.2f",r),
2301 nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
2302 fhPtInPerpCone->SetYTitle("#it{p}_{T in cone} (GeV/#it{c})");
2303 fhPtInPerpCone->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2304 outputContainer->Add(fhPtInPerpCone) ;
2306 fhEtaPhiTrack= new TH2F("hEtaPhiTrack",
2307 Form("#eta vs #phi of all Tracks"),
2308 netabins,-1,1,nphibins,0,TMath::TwoPi());
2309 fhEtaPhiTrack->SetXTitle("#eta");
2310 fhEtaPhiTrack->SetYTitle("#phi");
2311 outputContainer->Add(fhEtaPhiTrack) ;
2313 fhEtaPhiInConeTrack= new TH2F("hEtaPhiInConeTrack",
2314 Form("#eta vs #phi of Tracks in cone for #it{R} = %2.2f",r),
2315 netabins,-1,1,nphibins,0,TMath::TwoPi());
2316 fhEtaPhiInConeTrack->SetXTitle("#eta");
2317 fhEtaPhiInConeTrack->SetYTitle("#phi");
2318 outputContainer->Add(fhEtaPhiInConeTrack) ;
2320 fhConeSumPtVSUETracksEtaBand = new TH2F("hConeSumPtVSUETracksEtaBand",
2321 Form("#Sigma #it{p}_{T} in cone versus #Sigma #it{p}_{T} in eta band for tracks (before normalization), R=%2.2f",r),
2322 nptsumbins,ptsummin,ptsummax,2*nptsumbins,ptsummin,2*ptsummax);
2323 fhConeSumPtVSUETracksEtaBand->SetXTitle("#Sigma #it{p}_{T} cone (GeV/#it{c})");
2324 fhConeSumPtVSUETracksEtaBand->SetYTitle("#Sigma #it{p}_{T} UE (GeV/#it{c})");
2325 outputContainer->Add(fhConeSumPtVSUETracksEtaBand);
2327 fhConeSumPtVSUETracksPhiBand = new TH2F("hConeSumPtVSUETracksPhiBand",
2328 Form("#Sigma #it{p}_{T} in cone versus #Sigma #it{p}_{T} in phi band for tracks (before normalization), R=%2.2f",r),
2329 nptsumbins,ptsummin,ptsummax,8*nptsumbins,ptsummin,8*ptsummax);
2330 fhConeSumPtVSUETracksPhiBand->SetXTitle("#Sigma #it{p}_{T} cone (GeV/#it{c})");
2331 fhConeSumPtVSUETracksPhiBand->SetYTitle("#Sigma #it{p}_{T} UE (GeV/#it{c})");
2332 outputContainer->Add(fhConeSumPtVSUETracksPhiBand);
2336 if(GetIsolationCut()->GetParticleTypeInCone()==AliIsolationCut::kNeutralAndCharged )
2338 fhConeSumPtClustervsTrack = new TH2F("hConePtSumClustervsTrack",
2339 Form("Track vs Cluster #Sigma #it{p}_{T} in isolation cone for #it{R} = %2.2f",r),
2340 nptsumbins,ptsummin,ptsummax,nptsumbins,ptsummin,ptsummax);
2341 fhConeSumPtClustervsTrack->SetXTitle("#Sigma #it{p}_{T}^{cluster} (GeV/#it{c})");
2342 fhConeSumPtClustervsTrack->SetYTitle("#Sigma #it{p}_{T}^{track} (GeV/#it{c})");
2343 outputContainer->Add(fhConeSumPtClustervsTrack) ;
2345 fhConeSumPtClusterTrackFrac = new TH2F("hConePtSumClusterTrackFraction",
2346 Form("#Sigma #it{p}_{T}^{cluster}/#Sigma #it{p}_{T}^{track} in isolation cone for #it{R} = %2.2f",r),
2347 nptbins,ptmin,ptmax,200,0,5);
2348 fhConeSumPtClusterTrackFrac->SetYTitle("#Sigma #it{p}^{cluster}_{T} /#Sigma #it{p}_{T}^{track}");
2349 fhConeSumPtClusterTrackFrac->SetXTitle("#it{p}^{trigger}_{T} (GeV/#it{c})");
2350 outputContainer->Add(fhConeSumPtClusterTrackFrac) ;
2353 fhConePtLeadClustervsTrack = new TH2F("hConePtLeadClustervsTrack",
2354 Form("Track vs Cluster lead #it{p}_{T} in isolation cone for #it{R} = %2.2f",r),
2355 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2356 fhConePtLeadClustervsTrack->SetXTitle("#it{p}^{leading cluster}_{T} (GeV/#it{c})");
2357 fhConePtLeadClustervsTrack->SetYTitle("#it{p}^{leading track}_{T} (GeV/#it{c})");
2358 outputContainer->Add(fhConePtLeadClustervsTrack) ;
2360 fhConePtLeadClusterTrackFrac = new TH2F("hConePtLeadClusterTrackFraction",
2361 Form(" #it{p}^{leading cluster}_{T}/#it{p}^{leading track}_{T} in isolation cone for #it{R} = %2.2f",r),
2362 nptbins,ptmin,ptmax,200,0,5);
2363 fhConePtLeadClusterTrackFrac->SetYTitle("#it{p}^{leading cluster}_{T}/ #it{p}^{leading track}_{T}");
2364 fhConePtLeadClusterTrackFrac->SetXTitle("#it{p}^{trigger}_{T} (GeV/#it{c})");
2365 outputContainer->Add(fhConePtLeadClusterTrackFrac) ;
2368 if(fFillCellHistograms)
2370 fhConeSumPtCellvsTrack = new TH2F("hConePtSumCellvsTrack",
2371 Form("Track vs cell #Sigma #it{p}_{T} in isolation cone for #it{R} = %2.2f",r),
2372 nptsumbins,ptsummin,ptsummax,nptsumbins,ptsummin,ptsummax);
2373 fhConeSumPtCellvsTrack->SetXTitle("#Sigma #it{p}_{T} cell");
2374 fhConeSumPtCellvsTrack->SetYTitle("#Sigma #it{p}_{T} track");
2375 outputContainer->Add(fhConeSumPtCellvsTrack) ;
2377 fhConeSumPtCellTrack = new TH2F("hConePtSumCellTrack",
2378 Form("Track and Cell #Sigma #it{p}_{T} in isolation cone for #it{R} = %2.2f",r),
2379 nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
2380 fhConeSumPtCellTrack->SetYTitle("#Sigma #it{p}_{T}");
2381 fhConeSumPtCellTrack->SetXTitle("#it{p}_{T, trigger} (GeV/#it{c})");
2382 outputContainer->Add(fhConeSumPtCellTrack) ;
2384 fhConeSumPtCellTrackTrigEtaPhi = new TH2F("hConePtSumCellTrackTrigEtaPhi",
2385 Form("Trigger #eta vs #phi, #Sigma #it{p}_{T} in isolation cone for #it{R} = %2.2f",r),
2386 netabins,etamin,etamax,nphibins,phimin,phimax);
2387 fhConeSumPtCellTrackTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
2388 fhConeSumPtCellTrackTrigEtaPhi->SetXTitle("#eta_{trigger}");
2389 fhConeSumPtCellTrackTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2390 outputContainer->Add(fhConeSumPtCellTrackTrigEtaPhi) ;
2394 if(fFillUEBandSubtractHistograms)
2396 fhConeSumPtEtaUESub = new TH2F("hConeSumPtEtaUESub",
2397 Form("#Sigma #it{p}_{T} after bkg subtraction from eta band in the isolation cone for #it{R} = %2.2f",r),
2398 nptbins,ptmin,ptmax,2*nptsumbins,-ptsummax,ptsummax);
2399 fhConeSumPtEtaUESub->SetYTitle("#Sigma #it{p}_{T}");
2400 fhConeSumPtEtaUESub->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2401 outputContainer->Add(fhConeSumPtEtaUESub) ;
2403 fhConeSumPtPhiUESub = new TH2F("hConeSumPtPhiUESub",
2404 Form("#Sigma #it{p}_{T} after bkg subtraction from phi band in the isolation cone for #it{R} = %2.2f",r),
2405 nptbins,ptmin,ptmax,2*nptsumbins,-ptsummax,ptsummax);
2406 fhConeSumPtPhiUESub->SetYTitle("#Sigma #it{p}_{T}");
2407 fhConeSumPtPhiUESub->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2408 outputContainer->Add(fhConeSumPtPhiUESub) ;
2410 fhConeSumPtEtaUESubTrigEtaPhi = new TH2F("hConeSumPtEtaUESubTrigEtaPhi",
2411 Form("Trigger #eta vs #phi, #Sigma #it{p}_{T} after bkg subtraction from eta band in the isolation cone for #it{R} = %2.2f",r),
2412 netabins,etamin,etamax,nphibins,phimin,phimax);
2413 fhConeSumPtEtaUESubTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
2414 fhConeSumPtEtaUESubTrigEtaPhi->SetXTitle("#eta_{trigger}");
2415 fhConeSumPtEtaUESubTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2416 outputContainer->Add(fhConeSumPtEtaUESubTrigEtaPhi) ;
2418 fhConeSumPtPhiUESubTrigEtaPhi = new TH2F("hConeSumPtPhiUESubTrigEtaPhi",
2419 Form("Trigger #eta vs #phi, #Sigma #it{p}_{T} after bkg subtraction from phi band in the isolation cone for #it{R} = %2.2f",r),
2420 netabins,etamin,etamax,nphibins,phimin,phimax);
2421 fhConeSumPtPhiUESubTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
2422 fhConeSumPtPhiUESubTrigEtaPhi->SetXTitle("#eta_{trigger}");
2423 fhConeSumPtPhiUESubTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2424 outputContainer->Add(fhConeSumPtPhiUESubTrigEtaPhi) ;
2426 fhConeSumPtEtaUESubClustervsTrack = new TH2F("hConePtSumEtaUESubClustervsTrack",
2427 Form("Track vs Cluster #Sigma #it{p}_{T} UE sub eta band in isolation cone for #it{R} = %2.2f",r),
2428 2*nptsumbins,-ptsummax,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2429 fhConeSumPtEtaUESubClustervsTrack->SetXTitle("#Sigma #it{p}_{T} cluster");
2430 fhConeSumPtEtaUESubClustervsTrack->SetYTitle("#Sigma #it{p}_{T} track");
2431 outputContainer->Add(fhConeSumPtEtaUESubClustervsTrack) ;
2433 fhConeSumPtPhiUESubClustervsTrack = new TH2F("hConePhiUESubPtSumClustervsTrack",
2434 Form("Track vs Cluster #Sigma #it{p}_{T} UE sub phi band in isolation cone for #it{R} = %2.2f",r),
2435 2*nptsumbins,-ptsummax,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2436 fhConeSumPtPhiUESubClustervsTrack->SetXTitle("#Sigma #it{p}_{T} cluster");
2437 fhConeSumPtPhiUESubClustervsTrack->SetYTitle("#Sigma #it{p}_{T} track");
2438 outputContainer->Add(fhConeSumPtPhiUESubClustervsTrack) ;
2440 fhEtaBandClustervsTrack = new TH2F("hEtaBandClustervsTrack",
2441 Form("Track vs Cluster #Sigma #it{p}_{T} in Eta band in isolation cone for #it{R} = %2.2f",r),
2442 nptsumbins,ptsummin,ptsummax,nptsumbins,ptsummin,ptsummax);
2443 fhEtaBandClustervsTrack->SetXTitle("#Sigma #it{p}_{T} cluster");
2444 fhEtaBandClustervsTrack->SetYTitle("#Sigma #it{p}_{T} track");
2445 outputContainer->Add(fhEtaBandClustervsTrack) ;
2447 fhPhiBandClustervsTrack = new TH2F("hPhiBandClustervsTrack",
2448 Form("Track vs Cluster #Sigma #it{p}_{T} in Phi band in isolation cone for #it{R} = %2.2f",r),
2449 nptsumbins,ptsummin,ptsummax*4,nptsumbins,ptsummin,ptsummax*8);
2450 fhPhiBandClustervsTrack->SetXTitle("#Sigma #it{p}_{T} cluster");
2451 fhPhiBandClustervsTrack->SetYTitle("#Sigma #it{p}_{T} track");
2452 outputContainer->Add(fhPhiBandClustervsTrack) ;
2454 fhEtaBandNormClustervsTrack = new TH2F("hEtaBandNormClustervsTrack",
2455 Form("Track vs Cluster #Sigma #it{p}_{T} in Eta band in isolation cone for #it{R} = %2.2f",r),
2456 nptsumbins,ptsummin,ptsummax,nptsumbins,ptsummin,ptsummax);
2457 fhEtaBandNormClustervsTrack->SetXTitle("#Sigma #it{p}_{T} cluster");
2458 fhEtaBandNormClustervsTrack->SetYTitle("#Sigma #it{p}_{T} track");
2459 outputContainer->Add(fhEtaBandNormClustervsTrack) ;
2461 fhPhiBandNormClustervsTrack = new TH2F("hPhiBandNormClustervsTrack",
2462 Form("Track vs Cluster #Sigma #it{p}_{T} in Phi band in isolation cone for #it{R} = %2.2f",r),
2463 nptsumbins,ptsummin,ptsummax,nptsumbins,ptsummin,ptsummax);
2464 fhPhiBandNormClustervsTrack->SetXTitle("#Sigma #it{p}_{T} cluster");
2465 fhPhiBandNormClustervsTrack->SetYTitle("#Sigma #it{p}_{T} track");
2466 outputContainer->Add(fhPhiBandNormClustervsTrack) ;
2468 fhConeSumPtEtaUESubClustervsTrack = new TH2F("hConePtSumEtaUESubClustervsTrack",
2469 Form("Track vs Cluster #Sigma #it{p}_{T} UE sub eta band in isolation cone for #it{R} = %2.2f",r),
2470 2*nptsumbins,-ptsummax,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2471 fhConeSumPtEtaUESubClustervsTrack->SetXTitle("#Sigma #it{p}_{T} cluster");
2472 fhConeSumPtEtaUESubClustervsTrack->SetYTitle("#Sigma #it{p}_{T} track");
2473 outputContainer->Add(fhConeSumPtEtaUESubClustervsTrack) ;
2475 fhConeSumPtPhiUESubClustervsTrack = new TH2F("hConePhiUESubPtSumClustervsTrack",
2476 Form("Track vs Cluster #Sigma #it{p}_{T} UE sub phi band in isolation cone for #it{R} = %2.2f",r),
2477 2*nptsumbins,-ptsummax,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2478 fhConeSumPtPhiUESubClustervsTrack->SetXTitle("#Sigma #it{p}_{T} cluster");
2479 fhConeSumPtPhiUESubClustervsTrack->SetYTitle("#Sigma #it{p}_{T} track");
2480 outputContainer->Add(fhConeSumPtPhiUESubClustervsTrack) ;
2482 if(fFillCellHistograms)
2485 fhConeSumPtEtaUESubCellvsTrack = new TH2F("hConePtSumEtaUESubCellvsTrack",
2486 Form("Track vs Cell #Sigma #it{p}_{T} UE sub eta band in isolation cone for #it{R} = %2.2f",r),
2487 2*nptsumbins,-ptsummax,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2488 fhConeSumPtEtaUESubCellvsTrack->SetXTitle("#Sigma #it{p}_{T} cell");
2489 fhConeSumPtEtaUESubCellvsTrack->SetYTitle("#Sigma #it{p}_{T} track");
2490 outputContainer->Add(fhConeSumPtEtaUESubCellvsTrack) ;
2492 fhConeSumPtPhiUESubCellvsTrack = new TH2F("hConePhiUESubPtSumCellvsTrack",
2493 Form("Track vs Cell #Sigma #it{p}_{T} UE sub phi band in isolation cone for #it{R} = %2.2f",r),
2494 2*nptsumbins,-ptsummax,ptsummax,2*nptsumbins,-ptsummax,ptsummax);
2495 fhConeSumPtPhiUESubCellvsTrack->SetXTitle("#Sigma #it{p}_{T} cell");
2496 fhConeSumPtPhiUESubCellvsTrack->SetYTitle("#Sigma #it{p}_{T} track");
2497 outputContainer->Add(fhConeSumPtPhiUESubCellvsTrack) ;
2499 fhEtaBandCellvsTrack = new TH2F("hEtaBandCellvsTrack",
2500 Form("Track vs Cell #Sigma #it{p}_{T} in Eta band in isolation cone for #it{R} = %2.2f",r),
2501 nptsumbins,ptsummin,ptsummax,nptsumbins,ptsummin,ptsummax);
2502 fhEtaBandCellvsTrack->SetXTitle("#Sigma #it{p}_{T} cell");
2503 fhEtaBandCellvsTrack->SetYTitle("#Sigma #it{p}_{T} track");
2504 outputContainer->Add(fhEtaBandCellvsTrack) ;
2506 fhPhiBandCellvsTrack = new TH2F("hPhiBandCellvsTrack",
2507 Form("Track vs Cell #Sigma #it{p}_{T} in Phi band in isolation cone for #it{R} = %2.2f",r),
2508 nptsumbins,ptsummin,ptsummax*4,nptsumbins,ptsummin,ptsummax*8);
2509 fhPhiBandCellvsTrack->SetXTitle("#Sigma #it{p}_{T} cell");
2510 fhPhiBandCellvsTrack->SetYTitle("#Sigma #it{p}_{T} track");
2511 outputContainer->Add(fhPhiBandCellvsTrack) ;
2513 fhEtaBandNormCellvsTrack = new TH2F("hEtaBandNormCellvsTrack",
2514 Form("Track vs Cell #Sigma #it{p}_{T} in Eta band in isolation cone for #it{R} = %2.2f",r),
2515 nptsumbins,ptsummin,ptsummax,nptsumbins,ptsummin,ptsummax);
2516 fhEtaBandNormCellvsTrack->SetXTitle("#Sigma #it{p}_{T} cell");
2517 fhEtaBandNormCellvsTrack->SetYTitle("#Sigma #it{p}_{T} track");
2518 outputContainer->Add(fhEtaBandNormCellvsTrack) ;
2520 fhPhiBandNormCellvsTrack = new TH2F("hPhiBandNormCellvsTrack",
2521 Form("Track vs Cell #Sigma #it{p}_{T} in Phi band in isolation cone for #it{R} = %2.2f",r),
2522 nptsumbins,ptsummin,ptsummax,nptsumbins,ptsummin,ptsummax);
2523 fhPhiBandNormCellvsTrack->SetXTitle("#Sigma #it{p}_{T} cell");
2524 fhPhiBandNormCellvsTrack->SetYTitle("#Sigma #it{p}_{T} track");
2525 outputContainer->Add(fhPhiBandNormCellvsTrack) ;
2527 fhConeSumPtEtaUESubTrackCell = new TH2F("hConeSumPtEtaUESubTrackCell",
2528 Form("Tracks #Sigma #it{p}_{T} after bkg subtraction from eta band in the isolation cone for #it{R} = %2.2f",r),
2529 nptbins,ptmin,ptmax,2*nptsumbins,-ptsummax,ptsummax);
2530 fhConeSumPtEtaUESubTrackCell->SetYTitle("#Sigma #it{p}_{T}");
2531 fhConeSumPtEtaUESubTrackCell->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2532 outputContainer->Add(fhConeSumPtEtaUESubTrackCell) ;
2534 fhConeSumPtPhiUESubTrackCell = new TH2F("hConeSumPtPhiUESubTrackCell",
2535 Form("Tracks #Sigma #it{p}_{T} after bkg subtraction from phi band in the isolation cone for #it{R} = %2.2f",r),
2536 nptbins,ptmin,ptmax,2*nptsumbins,-ptsummax,ptsummax);
2537 fhConeSumPtPhiUESubTrackCell->SetYTitle("#Sigma #it{p}_{T}");
2538 fhConeSumPtPhiUESubTrackCell->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2539 outputContainer->Add(fhConeSumPtPhiUESubTrackCell) ;
2541 fhConeSumPtEtaUESubTrackCellTrigEtaPhi = new TH2F("hConeSumPtEtaUESubTrackCellTrigEtaPhi",
2542 Form("Trigger #eta vs #phi, Tracks #Sigma #it{p}_{T} after bkg subtraction from eta band in the isolation cone for #it{R} = %2.2f",r),
2543 netabins,etamin,etamax,nphibins,phimin,phimax);
2544 fhConeSumPtEtaUESubTrackCellTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
2545 fhConeSumPtEtaUESubTrackCellTrigEtaPhi->SetXTitle("#eta_{trigger}");
2546 fhConeSumPtEtaUESubTrackCellTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2547 outputContainer->Add(fhConeSumPtEtaUESubTrackCellTrigEtaPhi) ;
2549 fhConeSumPtPhiUESubTrackCellTrigEtaPhi = new TH2F("hConeSumPtPhiUESubTrackCellTrigEtaPhi",
2550 Form("Trigger #eta vs #phi, Tracks #Sigma #it{p}_{T} after bkg subtraction from phi band in the isolation cone for #it{R} = %2.2f",r),
2551 netabins,etamin,etamax,nphibins,phimin,phimax);
2552 fhConeSumPtPhiUESubTrackCellTrigEtaPhi->SetZTitle("#Sigma #it{p}_{T}");
2553 fhConeSumPtPhiUESubTrackCellTrigEtaPhi->SetXTitle("#eta_{trigger}");
2554 fhConeSumPtPhiUESubTrackCellTrigEtaPhi->SetYTitle("#phi_{trigger} (rad)");
2555 outputContainer->Add(fhConeSumPtPhiUESubTrackCellTrigEtaPhi) ;
2560 for(Int_t iso = 0; iso < 2; iso++)
2564 fhTrackMatchedDEta[iso] = new TH2F
2565 (Form("hTrackMatchedDEta%s",isoName[iso].Data()),
2566 Form("%s - d#eta of cluster-track vs cluster energy, %s",isoTitle[iso].Data(),parTitle.Data()),
2567 nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2568 fhTrackMatchedDEta[iso]->SetYTitle("d#eta");
2569 fhTrackMatchedDEta[iso]->SetXTitle("E_{cluster} (GeV)");
2571 fhTrackMatchedDPhi[iso] = new TH2F
2572 (Form("hTrackMatchedDPhi%s",isoName[iso].Data()),
2573 Form("%s - d#phi of cluster-track vs cluster energy, %s",isoTitle[iso].Data(),parTitle.Data()),
2574 nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2575 fhTrackMatchedDPhi[iso]->SetYTitle("d#phi (rad)");
2576 fhTrackMatchedDPhi[iso]->SetXTitle("E_{cluster} (GeV)");
2578 fhTrackMatchedDEtaDPhi[iso] = new TH2F
2579 (Form("hTrackMatchedDEtaDPhi%s",isoName[iso].Data()),
2580 Form("%s - d#eta vs d#phi of cluster-track, %s",isoTitle[iso].Data(),parTitle.Data()),
2581 nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
2582 fhTrackMatchedDEtaDPhi[iso]->SetYTitle("d#phi (rad)");
2583 fhTrackMatchedDEtaDPhi[iso]->SetXTitle("d#eta");
2585 outputContainer->Add(fhTrackMatchedDEta[iso]) ;
2586 outputContainer->Add(fhTrackMatchedDPhi[iso]) ;
2587 outputContainer->Add(fhTrackMatchedDEtaDPhi[iso]) ;
2589 fhdEdx[iso] = new TH2F
2590 (Form("hdEdx%s",isoName[iso].Data()),
2591 Form("%s - Matched track <d#it{E}/d#it{x}> vs cluster #it{E}, %s",isoTitle[iso].Data(),parTitle.Data()),
2592 nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
2593 fhdEdx[iso]->SetXTitle("#it{E} (GeV)");
2594 fhdEdx[iso]->SetYTitle("<d#it{E}/d#it{x}>");
2595 outputContainer->Add(fhdEdx[iso]);
2597 fhEOverP[iso] = new TH2F
2598 (Form("hEOverP%s",isoName[iso].Data()),
2599 Form("%s - Matched track #it{E}/#it{p} vs cluster, %s",isoTitle[iso].Data(),parTitle.Data()),
2600 nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
2601 fhEOverP[iso]->SetXTitle("#it{E} (GeV)");
2602 fhEOverP[iso]->SetYTitle("#it{E}/#it{p}");
2603 outputContainer->Add(fhEOverP[iso]);
2607 fhTrackMatchedMCParticle[iso] = new TH2F
2608 (Form("hTrackMatchedMCParticle%s",isoName[iso].Data()),
2609 Form("%s - Origin of particle vs cluster #it{E}, %s",isoTitle[iso].Data(),parTitle.Data()),
2610 nptbins,ptmin,ptmax,8,0,8);
2611 fhTrackMatchedMCParticle[iso]->SetXTitle("#it{E} (GeV)");
2612 //fhTrackMatchedMCParticle[iso]->SetYTitle("Particle type");
2614 fhTrackMatchedMCParticle[iso]->GetYaxis()->SetBinLabel(1 ,"Photon");
2615 fhTrackMatchedMCParticle[iso]->GetYaxis()->SetBinLabel(2 ,"Electron");
2616 fhTrackMatchedMCParticle[iso]->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
2617 fhTrackMatchedMCParticle[iso]->GetYaxis()->SetBinLabel(4 ,"Rest");
2618 fhTrackMatchedMCParticle[iso]->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
2619 fhTrackMatchedMCParticle[iso]->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
2620 fhTrackMatchedMCParticle[iso]->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
2621 fhTrackMatchedMCParticle[iso]->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
2623 outputContainer->Add(fhTrackMatchedMCParticle[iso]);
2629 fhELambda0[iso] = new TH2F
2630 (Form("hELambda0%s",isoName[iso].Data()),
2631 Form("%s cluster : #it{E} vs #lambda_{0}, %s",isoTitle[iso].Data(),parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2632 fhELambda0[iso]->SetYTitle("#lambda_{0}^{2}");
2633 fhELambda0[iso]->SetXTitle("#it{E} (GeV)");
2634 outputContainer->Add(fhELambda0[iso]) ;
2636 fhELambda1[iso] = new TH2F
2637 (Form("hELambda1%s",isoName[iso].Data()),
2638 Form("%s cluster: #it{E} vs #lambda_{1}, %s",isoTitle[iso].Data(),parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2639 fhELambda1[iso]->SetYTitle("#lambda_{1}^{2}");
2640 fhELambda1[iso]->SetXTitle("#it{E} (GeV)");
2641 outputContainer->Add(fhELambda1[iso]) ;
2643 fhPtLambda0[iso] = new TH2F
2644 (Form("hPtLambda0%s",isoName[iso].Data()),
2645 Form("%s cluster : #it{p}_{T} vs #lambda_{0}, %s",isoTitle[iso].Data(), parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2646 fhPtLambda0[iso]->SetYTitle("#lambda_{0}^{2}");
2647 fhPtLambda0[iso]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2648 outputContainer->Add(fhPtLambda0[iso]) ;
2650 if(fFillTaggedDecayHistograms)
2652 for(Int_t ibit = 0; ibit < fNDecayBits; ibit++)
2654 fhPtLambda0Decay[iso][ibit] = new TH2F
2655 (Form("hPtLambda0Decay%s_bit%d",isoName[iso].Data(),fDecayBits[ibit]),
2656 Form("%s cluster : #it{p}_{T} vs #lambda_{0}, decay bit %d, %s",isoTitle[iso].Data(), fDecayBits[ibit], parTitle.Data()),
2657 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2658 fhPtLambda0Decay[iso][ibit]->SetYTitle("#lambda_{0}^{2}");
2659 fhPtLambda0Decay[iso][ibit]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2660 outputContainer->Add(fhPtLambda0Decay[iso][ibit]) ;
2666 for(Int_t imc = 0; imc < fgkNmcTypes; imc++)
2668 fhPtLambda0MC[imc][iso] = new TH2F(Form("hPtLambda0%s_MC%s",isoName[iso].Data(),mcPartName[imc].Data()),
2669 Form("%s cluster : #it{p}_{T} vs #lambda_{0}: %s %s",isoTitle[iso].Data(),mcPartType[imc].Data(),parTitle.Data()),
2670 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2671 fhPtLambda0MC[imc][iso]->SetYTitle("#lambda_{0}^{2}");
2672 fhPtLambda0MC[imc][iso]->SetXTitle("#it{p}_{T}(GeV/#it{c})");
2673 outputContainer->Add( fhPtLambda0MC[imc][iso]) ;
2677 if(fIsoDetector=="EMCAL" && GetFirstSMCoveredByTRD() >= 0)
2679 fhPtLambda0TRD[iso] = new TH2F
2680 (Form("hPtLambda0TRD%s",isoName[iso].Data()),
2681 Form("%s cluster: #it{p}_{T} vs #lambda_{0}, SM behind TRD, %s",isoTitle[iso].Data(),parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2682 fhPtLambda0TRD[iso]->SetYTitle("#lambda_{0}^{2}");
2683 fhPtLambda0TRD[iso]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2684 outputContainer->Add(fhPtLambda0TRD[iso]) ;
2686 fhELambda0TRD[iso] = new TH2F
2687 (Form("hELambda0TRD%s",isoName[iso].Data()),
2688 Form("%s cluster: #it{E} vs #lambda_{0}, SM behind TRD, %s",isoTitle[iso].Data(),parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2689 fhELambda0TRD[iso]->SetYTitle("#lambda_{0}^{2}");
2690 fhELambda0TRD[iso]->SetXTitle("#it{E} (GeV)");
2691 outputContainer->Add(fhELambda0TRD[iso]) ;
2693 fhELambda1TRD[iso] = new TH2F
2694 (Form("hELambda1TRD%s",isoName[iso].Data()),
2695 Form("%s cluster: #it{E} vs #lambda_{1}, SM behind TRD, %s",isoTitle[iso].Data(),parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2696 fhELambda1TRD[iso]->SetYTitle("#lambda_{1}^{2}");
2697 fhELambda1TRD[iso]->SetXTitle("#it{E} (GeV)");
2698 outputContainer->Add(fhELambda1TRD[iso]) ;
2701 if(fFillNLMHistograms)
2703 fhNLocMax[iso] = new TH2F
2704 (Form("hNLocMax%s",isoName[iso].Data()),
2705 Form("%s - Number of local maxima in cluster, %s",isoTitle[iso].Data(),parTitle.Data()),
2706 nptbins,ptmin,ptmax,10,0,10);
2707 fhNLocMax[iso]->SetYTitle("#it{NLM}");
2708 fhNLocMax[iso]->SetXTitle("#it{E} (GeV)");
2709 outputContainer->Add(fhNLocMax[iso]) ;
2711 fhELambda0LocMax1[iso] = new TH2F
2712 (Form("hELambda0LocMax1%s",isoName[iso].Data()),
2713 Form("%s cluster (#eta) pairs: #it{E} vs #lambda_{0}, #it{NLM}=1, %s",isoTitle[iso].Data(),parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2714 fhELambda0LocMax1[iso]->SetYTitle("#lambda_{0}^{2}");
2715 fhELambda0LocMax1[iso]->SetXTitle("#it{E} (GeV)");
2716 outputContainer->Add(fhELambda0LocMax1[iso]) ;
2718 fhELambda1LocMax1[iso] = new TH2F
2719 (Form("hELambda1LocMax1%s",isoName[iso].Data()),
2720 Form("%s cluster (#eta) pairs: #it{E} vs #lambda_{1}, #it{NLM}=1, %s",isoTitle[iso].Data(),parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2721 fhELambda1LocMax1[iso]->SetYTitle("#lambda_{1}^{2}");
2722 fhELambda1LocMax1[iso]->SetXTitle("#it{E} (GeV)");
2723 outputContainer->Add(fhELambda1LocMax1[iso]) ;
2725 fhELambda0LocMax2[iso] = new TH2F
2726 (Form("hELambda0LocMax2%s",isoName[iso].Data()),
2727 Form("%s cluster (#eta) pairs: #it{E} vs #lambda_{0}, #it{NLM}=2, %s",isoTitle[iso].Data(),parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2728 fhELambda0LocMax2[iso]->SetYTitle("#lambda_{0}^{2}");
2729 fhELambda0LocMax2[iso]->SetXTitle("#it{E} (GeV)");
2730 outputContainer->Add(fhELambda0LocMax2[iso]) ;
2732 fhELambda1LocMax2[iso] = new TH2F
2733 (Form("hELambda1LocMax2%s",isoName[iso].Data()),
2734 Form("%s cluster (#eta) pairs: #it{E} vs #lambda_{1}, #it{NLM}=2, %s",isoTitle[iso].Data(),parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2735 fhELambda1LocMax2[iso]->SetYTitle("#lambda_{1}^{2}");
2736 fhELambda1LocMax2[iso]->SetXTitle("#it{E} (GeV)");
2737 outputContainer->Add(fhELambda1LocMax2[iso]) ;
2739 fhELambda0LocMaxN[iso] = new TH2F
2740 ( Form("hELambda0LocMaxN%s",isoName[iso].Data()),
2741 Form("%s cluster (#eta) pairs: #it{E} vs #lambda_{0}, #it{NLM}>2, %s",isoTitle[iso].Data(),parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2742 fhELambda0LocMaxN[iso]->SetYTitle("#lambda_{0}^{2}");
2743 fhELambda0LocMaxN[iso]->SetXTitle("#it{E} (GeV)");
2744 outputContainer->Add(fhELambda0LocMaxN[iso]) ;
2746 fhELambda1LocMaxN[iso] = new TH2F
2747 (Form("hELambda1LocMaxN%s",isoName[iso].Data()),
2748 Form("%s cluster (#eta) pairs: #it{E} vs #lambda_{1}, #it{NLM}>2, %s",isoTitle[iso].Data(),parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2749 fhELambda1LocMaxN[iso]->SetYTitle("#lambda_{1}^{2}");
2750 fhELambda1LocMaxN[iso]->SetXTitle("#it{E} (GeV)");
2751 outputContainer->Add(fhELambda1LocMaxN[iso]) ;
2754 } // control histograms for isolated and non isolated objects
2757 if(fFillPileUpHistograms)
2759 fhPtTrackInConeOtherBC = new TH2F("hPtTrackInConeOtherBC",
2760 Form("#it{p}_{T} of tracks in isolation cone for #it{R} = %2.2f, TOF from BC!=0",r),
2761 nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
2762 fhPtTrackInConeOtherBC->SetYTitle("#it{p}_{T in cone} (GeV/#it{c})");
2763 fhPtTrackInConeOtherBC->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2764 outputContainer->Add(fhPtTrackInConeOtherBC) ;
2766 fhPtTrackInConeOtherBCPileUpSPD = new TH2F("hPtTrackInConeOtherBCPileUpSPD",
2767 Form("#it{p}_{T} of tracks in isolation cone for #it{R} = %2.2f, TOF from BC!=0, pile-up from SPD",r),
2768 nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
2769 fhPtTrackInConeOtherBCPileUpSPD->SetYTitle("#it{p}_{T in cone} (GeV/#it{c})");
2770 fhPtTrackInConeOtherBCPileUpSPD->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2771 outputContainer->Add(fhPtTrackInConeOtherBCPileUpSPD) ;
2773 fhPtTrackInConeBC0 = new TH2F("hPtTrackInConeBC0",
2774 Form("#it{p}_{T} of tracks in isolation cone for #it{R} = %2.2f, TOF from BC==0",r),
2775 nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
2776 fhPtTrackInConeBC0->SetYTitle("#it{p}_{T in cone} (GeV/#it{c})");
2777 fhPtTrackInConeBC0->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2778 outputContainer->Add(fhPtTrackInConeBC0) ;
2780 fhPtTrackInConeVtxBC0 = new TH2F("hPtTrackInConeVtxBC0",
2781 Form("#it{p}_{T} of tracks in isolation cone for #it{R} = %2.2f, TOF from BC==0",r),
2782 nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
2783 fhPtTrackInConeVtxBC0->SetYTitle("#it{p}_{T in cone} (GeV/#it{c})");
2784 fhPtTrackInConeVtxBC0->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2785 outputContainer->Add(fhPtTrackInConeVtxBC0) ;
2788 fhPtTrackInConeBC0PileUpSPD = new TH2F("hPtTrackInConeBC0PileUpSPD",
2789 Form("#it{p}_{T} of tracks in isolation cone for #it{R} = %2.2f, TOF from BC==0, pile-up from SPD",r),
2790 nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
2791 fhPtTrackInConeBC0PileUpSPD->SetYTitle("#it{p}_{T in cone} (GeV/#it{c})");
2792 fhPtTrackInConeBC0PileUpSPD->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2793 outputContainer->Add(fhPtTrackInConeBC0PileUpSPD) ;
2796 for (Int_t i = 0; i < 7 ; i++)
2798 fhPtInConePileUp[i] = new TH2F(Form("hPtInConePileUp%s",pileUpName[i].Data()),
2799 Form("#it{p}_{T} in isolation cone for #it{R} = %2.2f, from pile-up (%s)",r,pileUpName[i].Data()),
2800 nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
2801 fhPtInConePileUp[i]->SetYTitle("#it{p}_{T in cone} (GeV/#it{c})");
2802 fhPtInConePileUp[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2803 outputContainer->Add(fhPtInConePileUp[i]) ;
2809 // For histograms in arrays, index in the array, corresponding to any particle origin
2811 for(Int_t i = 0; i < fgkNmcPrimTypes; i++)
2813 fhEPrimMC[i] = new TH1F(Form("hEPrim_MC%s",ppname[i].Data()),
2814 Form("primary photon %s : #it{E}, %s",pptype[i].Data(),parTitle.Data()),
2815 nptbins,ptmin,ptmax);
2816 fhEPrimMC[i]->SetXTitle("#it{E} (GeV)");
2817 outputContainer->Add(fhEPrimMC[i]) ;
2819 fhPtPrimMC[i] = new TH1F(Form("hPtPrim_MC%s",ppname[i].Data()),
2820 Form("primary photon %s : #it{p}_{T}, %s",pptype[i].Data(),parTitle.Data()),
2821 nptbins,ptmin,ptmax);
2822 fhPtPrimMC[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2823 outputContainer->Add(fhPtPrimMC[i]) ;
2825 fhPtPrimMCiso[i] = new TH1F(Form("hPtPrim_MCiso%s",ppname[i].Data()),
2826 Form("primary isolated photon %s : #it{p}_{T}, %s",pptype[i].Data(),parTitle.Data()),
2827 nptbins,ptmin,ptmax);
2828 fhPtPrimMCiso[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2829 outputContainer->Add(fhPtPrimMCiso[i]) ;
2831 fhEtaPrimMC[i] = new TH2F(Form("hEtaPrim_MC%s",ppname[i].Data()),
2832 Form("primary photon %s : #eta vs #it{p}_{T}, %s",pptype[i].Data(),parTitle.Data()),
2833 nptbins,ptmin,ptmax,200,-2,2);
2834 fhEtaPrimMC[i]->SetYTitle("#eta");
2835 fhEtaPrimMC[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2836 outputContainer->Add(fhEtaPrimMC[i]) ;
2838 fhPhiPrimMC[i] = new TH2F(Form("hPhiPrim_MC%s",ppname[i].Data()),
2839 Form("primary photon %s : #phi vs #it{p}_{T}, %s",pptype[i].Data(),parTitle.Data()),
2840 nptbins,ptmin,ptmax,200,0.,TMath::TwoPi());
2841 fhPhiPrimMC[i]->SetYTitle("#phi");
2842 fhPhiPrimMC[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2843 outputContainer->Add(fhPhiPrimMC[i]) ;
2846 if(fMakePrimaryPi0DecayStudy)
2848 fhPtPrimMCPi0DecayPairAcceptInConeLowPt = new TH1F("hPtPrim_MCPhotonPi0DecayPairAcceptInConeLowPt",
2849 Form("primary photon %s : #it{p}_{T}, pair in cone, %s",pptype[kmcPrimPi0Decay].Data(),parTitle.Data()),
2850 nptbins,ptmin,ptmax);
2851 fhPtPrimMCPi0DecayPairAcceptInConeLowPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2852 outputContainer->Add(fhPtPrimMCPi0DecayPairAcceptInConeLowPt) ;
2854 fhPtPrimMCPi0DecayIsoPairAcceptInConeLowPt = new TH1F("hPtPrim_MCisoPhotonPi0DecayPairAcceptInConeLowPt",
2855 Form("isolated primary photon %s, pair in cone : #it{p}_{T}, %s",
2856 pptype[kmcPrimPi0Decay].Data(),parTitle.Data()),
2857 nptbins,ptmin,ptmax);
2858 fhPtPrimMCPi0DecayIsoPairAcceptInConeLowPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2859 outputContainer->Add(fhPtPrimMCPi0DecayIsoPairAcceptInConeLowPt) ;
2861 fhPtPrimMCPi0DecayPairAcceptInConeLowPtNoOverlap = new TH1F("hPtPrim_MCPhotonPi0DecayPairAcceptInConeLowPtNoOverlap",
2862 Form("primary photon %s, no overlap, pair in cone : #it{p}_{T}, %s",
2863 pptype[kmcPrimPi0Decay].Data(),parTitle.Data()),
2864 nptbins,ptmin,ptmax);
2865 fhPtPrimMCPi0DecayPairAcceptInConeLowPtNoOverlap->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2866 outputContainer->Add(fhPtPrimMCPi0DecayPairAcceptInConeLowPtNoOverlap) ;
2868 fhPtPrimMCPi0DecayIsoPairAcceptInConeLowPtNoOverlap = new TH1F("hPtPrim_MCisoPhotonPi0DecayPairAcceptInConeLowPtNoOverlap",
2869 Form("isolated primary photon %s, pair in cone,no overlap : #it{p}_{T}, %s",
2870 pptype[kmcPrimPi0Decay].Data(),parTitle.Data()),
2871 nptbins,ptmin,ptmax);
2872 fhPtPrimMCPi0DecayIsoPairAcceptInConeLowPtNoOverlap->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2873 outputContainer->Add(fhPtPrimMCPi0DecayIsoPairAcceptInConeLowPtNoOverlap) ;
2875 fhPtPrimMCPi0DecayPairAcceptInConeLowPtNoOverlapCaloE = new TH1F("hPtPrim_MCPhotonPi0DecayPairAcceptInConeLowPtNoOverlapCaloE",
2876 Form("primary photon %s, no overlap, pair in cone, E > calo min: #it{p}_{T}, %s",
2877 pptype[kmcPrimPi0Decay].Data(),parTitle.Data()),
2878 nptbins,ptmin,ptmax);
2879 fhPtPrimMCPi0DecayPairAcceptInConeLowPtNoOverlapCaloE->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2880 outputContainer->Add(fhPtPrimMCPi0DecayPairAcceptInConeLowPtNoOverlapCaloE) ;
2882 fhPtPrimMCPi0DecayIsoPairAcceptInConeLowPtNoOverlapCaloE = new TH1F("hPtPrim_MCisoPhotonPi0DecayPairAcceptInConeLowPtNoOverlapCaloE",
2883 Form("isolated primary photon %s, pair in cone,no overlap, E > calo min: #it{p}_{T}, %s",
2884 pptype[kmcPrimPi0Decay].Data(),parTitle.Data()),
2885 nptbins,ptmin,ptmax);
2886 fhPtPrimMCPi0DecayIsoPairAcceptInConeLowPtNoOverlapCaloE->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2887 outputContainer->Add(fhPtPrimMCPi0DecayIsoPairAcceptInConeLowPtNoOverlapCaloE) ;
2890 fhPtPrimMCPi0DecayPairNoOverlap = new TH1F("hPtPrim_MCPhotonPi0DecayPairNoOverlap",
2891 Form("primary photon %s, no overlap: #it{p}_{T}, %s",
2892 pptype[kmcPrimPi0Decay].Data(),parTitle.Data()),
2893 nptbins,ptmin,ptmax);
2894 fhPtPrimMCPi0DecayPairNoOverlap->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2895 outputContainer->Add(fhPtPrimMCPi0DecayPairNoOverlap) ;
2897 fhPtPrimMCPi0DecayIsoPairNoOverlap = new TH1F("hPtPrim_MCisoPhotonPi0DecayPairNoOverlap",
2898 Form("isolated primary photon %s, no overlap: #it{p}_{T}, %s",
2899 pptype[kmcPrimPi0Decay].Data(),parTitle.Data()),
2900 nptbins,ptmin,ptmax);
2901 fhPtPrimMCPi0DecayIsoPairNoOverlap->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2902 outputContainer->Add(fhPtPrimMCPi0DecayIsoPairNoOverlap) ;
2904 fhPtPrimMCPi0DecayPairOutOfCone = new TH1F("hPtPrim_MCPhotonPi0DecayPairOutOfCone",
2905 Form("primary photon %s : #it{p}_{T}, pair out of cone, %s",pptype[kmcPrimPi0Decay].Data(),parTitle.Data()),
2906 nptbins,ptmin,ptmax);
2907 fhPtPrimMCPi0DecayPairOutOfCone->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2908 outputContainer->Add(fhPtPrimMCPi0DecayPairOutOfCone) ;
2910 fhPtPrimMCPi0DecayIsoPairOutOfCone = new TH1F("hPtPrim_MCisoPhotonPi0DecayPairOutOfCone",
2911 Form("isolated primary photon %s, pair out of cone : #it{p}_{T}, %s",
2912 pptype[kmcPrimPi0Decay].Data(),parTitle.Data()),
2913 nptbins,ptmin,ptmax);
2914 fhPtPrimMCPi0DecayIsoPairOutOfCone->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2915 outputContainer->Add(fhPtPrimMCPi0DecayIsoPairOutOfCone) ;
2917 fhPtPrimMCPi0DecayPairOutOfAcceptance = new TH1F("hPtPrim_MCPhotonPi0DecayPairOutOfAcceptance",
2918 Form("primary photon %s : #it{p}_{T}, pair out of acceptance, %s",pptype[kmcPrimPi0Decay].Data(),parTitle.Data()),
2919 nptbins,ptmin,ptmax);
2920 fhPtPrimMCPi0DecayPairOutOfAcceptance->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2921 outputContainer->Add(fhPtPrimMCPi0DecayPairOutOfAcceptance) ;
2923 fhPtPrimMCPi0DecayIsoPairOutOfAcceptance = new TH1F("hPtPrim_MCisoPhotonPi0DecayPairOutOfAcceptance",
2924 Form("isolated primary photon %s, pair out of acceptance : #it{p}_{T}, %s",
2925 pptype[kmcPrimPi0Decay].Data(),parTitle.Data()),
2926 nptbins,ptmin,ptmax);
2927 fhPtPrimMCPi0DecayIsoPairOutOfAcceptance->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2928 outputContainer->Add(fhPtPrimMCPi0DecayIsoPairOutOfAcceptance) ;
2930 fhPtPrimMCPi0Overlap = new TH1F("hPtPrim_MCPi0Overlap",
2931 Form("primary %s, overlap: #it{p}_{T}, %s",
2932 pptype[kmcPrimPi0].Data(),parTitle.Data()),
2933 nptbins,ptmin,ptmax);
2934 fhPtPrimMCPi0Overlap->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2935 outputContainer->Add(fhPtPrimMCPi0Overlap) ;
2937 fhPtPrimMCPi0IsoOverlap = new TH1F("hPtPrim_MCisoPi0Overlap",
2938 Form("primary %s, overlap: #it{p}_{T}, %s",
2939 pptype[kmcPrimPi0].Data(),parTitle.Data()),
2940 nptbins,ptmin,ptmax);
2941 fhPtPrimMCPi0IsoOverlap->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2942 outputContainer->Add(fhPtPrimMCPi0IsoOverlap) ;
2952 const Int_t buffersize = 255;
2953 char name[buffersize];
2954 char title[buffersize];
2955 for(Int_t icone = 0; icone<fNCones; icone++)
2957 // sum pt in cone vs. pt leading
2958 snprintf(name, buffersize,"hSumPtLeadingPt_Cone_%d",icone);
2959 snprintf(title, buffersize,"#Sigma #it{p}_{T} in isolation cone for #it{R} = %2.2f",fConeSizes[icone]);
2960 fhSumPtLeadingPt[icone] = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
2961 fhSumPtLeadingPt[icone] ->SetYTitle("#sum_{cone}#it{p}_{T} (GeV/#it{c})");//#Sigma #it{p}_{T}
2962 fhSumPtLeadingPt[icone] ->SetXTitle("#it{p}_{T}^{leading} (GeV/#it{c})");
2963 outputContainer->Add(fhSumPtLeadingPt[icone]) ;
2965 // pt in cone vs. pt leading
2966 snprintf(name, buffersize,"hPtLeadingPt_Cone_%d",icone);
2967 snprintf(title, buffersize,"#it{p}_{T} in isolation cone for #it{R} = %2.2f",fConeSizes[icone]);
2968 fhPtLeadingPt[icone] = new TH2F(name, title, nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
2969 fhPtLeadingPt[icone] ->SetYTitle("#it{p}_{T}^{cone} (GeV/#it{c})");
2970 fhPtLeadingPt[icone] ->SetXTitle("#it{p}_{T}^{leading} (GeV/#it{c})");
2971 outputContainer->Add(fhPtLeadingPt[icone]) ;
2973 // sum pt in cone vs. pt leading in the forward region (for background subtraction studies)
2974 snprintf(name, buffersize,"hPerpSumPtLeadingPt_Cone_%d",icone);
2975 snprintf(title, buffersize,"#Sigma #it{p}_{T} in isolation cone for #it{R} = %2.2f",fConeSizes[icone]);
2976 fhPerpSumPtLeadingPt[icone] = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
2977 fhPerpSumPtLeadingPt[icone] ->SetYTitle("#sum_{cone}#it{p}_{T} (GeV/#it{c})");//#Sigma #it{p}_{T}
2978 fhPerpSumPtLeadingPt[icone] ->SetXTitle("#it{p}_{T}^{leading} (GeV/#it{c})");
2979 outputContainer->Add(fhPerpSumPtLeadingPt[icone]) ;
2981 // pt in cone vs. pt leading in the forward region (for background subtraction studies)
2982 snprintf(name, buffersize,"hPerpPtLeadingPt_Cone_%d",icone);
2983 snprintf(title, buffersize,"#it{p}_{T} in isolation cone for #it{R} = %2.2f",fConeSizes[icone]);
2984 fhPerpPtLeadingPt[icone] = new TH2F(name, title, nptbins,ptmin,ptmax,nptinconebins,ptinconemin,ptinconemax);
2985 fhPerpPtLeadingPt[icone] ->SetYTitle("#it{p}_{T}^{cone} (GeV/#it{c})");
2986 fhPerpPtLeadingPt[icone] ->SetXTitle("#it{p}_{T}^{leading} (GeV/#it{c})");
2987 outputContainer->Add(fhPerpPtLeadingPt[icone]) ;
2991 for(Int_t imc = 0; imc < fgkNmcTypes; imc++)
2993 snprintf(name , buffersize,"hSumPtLeadingPt_MC%s_Cone_%d",mcPartName[imc].Data(),icone);
2994 snprintf(title, buffersize,"Candidate %s #it{p}_{T} vs cone #Sigma #it{p}_{T} for #it{R}=%2.2f",mcPartType[imc].Data(),fConeSizes[icone]);
2995 fhSumPtLeadingPtMC[imc][icone] = new TH2F(name, title,nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
2996 fhSumPtLeadingPtMC[imc][icone]->SetYTitle("#Sigma #it{p}_{T} (GeV/#it{c})");
2997 fhSumPtLeadingPtMC[imc][icone]->SetXTitle("#it{p}_{T}(GeV/#it{c})");
2998 outputContainer->Add(fhSumPtLeadingPtMC[imc][icone]) ;
3002 for(Int_t ipt = 0; ipt<fNPtThresFrac;ipt++)
3004 snprintf(name, buffersize,"hPtThres_Cone_%d_Pt%d",icone,ipt);
3005 snprintf(title, buffersize,"Isolated candidate #it{p}_{T} distribution for #it{R} = %2.2f and #it{p}_{T}^{th} = %2.2f GeV/#it{c}",fConeSizes[icone],fPtThresholds[ipt]);
3006 fhPtThresIsolated[icone][ipt] = new TH1F(name, title,nptbins,ptmin,ptmax);
3007 fhPtThresIsolated[icone][ipt]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3008 outputContainer->Add(fhPtThresIsolated[icone][ipt]) ;
3010 snprintf(name, buffersize,"hPtFrac_Cone_%d_Pt%d",icone,ipt);
3011 snprintf(title, buffersize,"Isolated candidate #it{p}_{T} distribution for #it{R} = %2.2f and #it{p}_{T}^{fr} = %2.2f GeV/#it{c}",fConeSizes[icone],fPtFractions[ipt]);
3012 fhPtFracIsolated[icone][ipt] = new TH1F(name, title,nptbins,ptmin,ptmax);
3013 fhPtFracIsolated[icone][ipt]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3014 outputContainer->Add(fhPtFracIsolated[icone][ipt]) ;
3016 snprintf(name, buffersize,"hSumPt_Cone_%d_Pt%d",icone,ipt);
3017 snprintf(title, buffersize,"Isolated candidate #it{p}_{T} distribution for #it{R} = %2.2f and #it{p}_{T}^{sum} = %2.2f GeV/#it{c}",fConeSizes[icone],fSumPtThresholds[ipt]);
3018 fhSumPtIsolated[icone][ipt] = new TH1F(name, title,nptbins,ptmin,ptmax);
3019 // fhSumPtIsolated[icone][ipt]->SetYTitle("#Sigma #it{p}_{T} (GeV/#it{c})");
3020 fhSumPtIsolated[icone][ipt]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3021 outputContainer->Add(fhSumPtIsolated[icone][ipt]) ;
3023 snprintf(name, buffersize,"hPtSumDensity_Cone_%d_Pt%d",icone,ipt);
3024 snprintf(title, buffersize,"Isolated candidate #it{p}_{T} distribution for density in #it{R} = %2.2f and #it{p}_{T}^{sum} = %2.2f GeV/#it{c}",fConeSizes[icone],fSumPtThresholds[ipt]);
3025 fhPtSumDensityIso[icone][ipt] = new TH1F(name, title,nptbins,ptmin,ptmax);//,nptsumbins,ptsummin,ptsummax);
3026 //fhPtSumIsolated[icone][ipt]->SetYTitle("#Sigma #it{p}_{T} (GeV/#it{c})");
3027 fhPtSumDensityIso[icone][ipt]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3028 outputContainer->Add(fhPtSumDensityIso[icone][ipt]) ;
3030 snprintf(name, buffersize,"hPtFracPtSum_Cone_%d_Pt%d",icone,ipt);
3031 snprintf(title, buffersize,"Isolated candidate #it{p}_{T} distribution for PtFracPtSum in #it{R} = %2.2f and #it{p}_{T}^{fr} = %2.2f GeV/#it{c}",fConeSizes[icone],fPtFractions[ipt]);
3032 fhPtFracPtSumIso[icone][ipt] = new TH1F(name, title,nptbins,ptmin,ptmax);//,nptsumbins,ptsummin,ptsummax);
3033 //fhPtSumIsolated[icone][ipt]->SetYTitle("#Sigma #it{p}_{T} (GeV/#it{c})");
3034 fhPtFracPtSumIso[icone][ipt]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3035 outputContainer->Add(fhPtFracPtSumIso[icone][ipt]) ;
3038 snprintf(name, buffersize,"hEtaPhiPtThres_Cone_%d_Pt%d",icone,ipt);
3039 snprintf(title, buffersize,"Isolated candidate #eta:#phi distribution for #it{R} = %2.2f and #it{p}_{T}^{th} = %2.2f GeV/#it{c}",fConeSizes[icone],fPtThresholds[ipt]);
3040 fhEtaPhiPtThresIso[icone][ipt] = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
3041 fhEtaPhiPtThresIso[icone][ipt]->SetXTitle("#eta");
3042 fhEtaPhiPtThresIso[icone][ipt]->SetYTitle("#phi");
3043 outputContainer->Add(fhEtaPhiPtThresIso[icone][ipt]) ;
3045 snprintf(name, buffersize,"hEtaPhiPtFrac_Cone_%d_Pt%d",icone,ipt);
3046 snprintf(title, buffersize,"Isolated candidate #eta:#phi distribution for #it{R} = %2.2f and #it{p}_{T}^{fr} = %2.2f GeV/#it{c}",fConeSizes[icone],fPtFractions[ipt]);
3047 fhEtaPhiPtFracIso[icone][ipt] = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
3048 fhEtaPhiPtFracIso[icone][ipt]->SetXTitle("#eta");
3049 fhEtaPhiPtFracIso[icone][ipt]->SetYTitle("#phi");
3050 outputContainer->Add(fhEtaPhiPtFracIso[icone][ipt]) ;
3052 snprintf(name, buffersize,"hEtaPhiPtSum_Cone_%d_Pt%d",icone,ipt);
3053 snprintf(title, buffersize,"Isolated candidate #eta:#phi distribution for #it{R} = %2.2f and #it{p}_{T}^{sum} = %2.2f GeV/#it{c}",fConeSizes[icone],fSumPtThresholds[ipt]);
3054 fhEtaPhiPtSumIso[icone][ipt] = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
3055 fhEtaPhiPtSumIso[icone][ipt]->SetXTitle("#eta");
3056 fhEtaPhiPtSumIso[icone][ipt]->SetYTitle("#phi");
3057 outputContainer->Add(fhEtaPhiPtSumIso[icone][ipt]) ;
3059 snprintf(name, buffersize,"hEtaPhiSumDensity_Cone_%d_Pt%d",icone,ipt);
3060 snprintf(title, buffersize,"Isolated candidate #eta:#phi distribution for density #it{R} = %2.2f and #it{p}_{T}^{sum} = %2.2f GeV/#it{c}",fConeSizes[icone],fSumPtThresholds[ipt]);
3061 fhEtaPhiSumDensityIso[icone][ipt] = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
3062 fhEtaPhiSumDensityIso[icone][ipt]->SetXTitle("#eta");
3063 fhEtaPhiSumDensityIso[icone][ipt]->SetYTitle("#phi");
3064 outputContainer->Add(fhEtaPhiSumDensityIso[icone][ipt]) ;
3066 snprintf(name, buffersize,"hEtaPhiFracPtSum_Cone_%d_Pt%d",icone,ipt);
3067 snprintf(title, buffersize,"Isolated candidate #eta:#phi distribution for FracPtSum #it{R} = %2.2f and #it{p}_{T}^{fr} = %2.2f GeV/#it{c}",fConeSizes[icone],fPtFractions[ipt]);
3068 fhEtaPhiFracPtSumIso[icone][ipt] = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
3069 fhEtaPhiFracPtSumIso[icone][ipt]->SetXTitle("#eta");
3070 fhEtaPhiFracPtSumIso[icone][ipt]->SetYTitle("#phi");
3071 outputContainer->Add(fhEtaPhiFracPtSumIso[icone][ipt]) ;
3073 if(fFillTaggedDecayHistograms)
3075 // pt decays isolated
3076 snprintf(name, buffersize,"hPtThres_Decay_Cone_%d_Pt%d",icone,ipt);
3077 snprintf(title, buffersize,"Isolated decay candidate #it{p}_{T} distribution for #it{R} = %2.2f and #it{p}_{T}^{th} = %2.2f GeV/#it{c}",fConeSizes[icone],fPtThresholds[ipt]);
3078 fhPtPtThresDecayIso[icone][ipt] = new TH1F(name, title,nptbins,ptmin,ptmax);
3079 fhPtPtThresDecayIso[icone][ipt]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3080 outputContainer->Add(fhPtPtThresDecayIso[icone][ipt]) ;
3082 snprintf(name, buffersize,"hPtFrac_Decay_Cone_%d_Pt%d",icone,ipt);
3083 snprintf(title, buffersize,"Isolated decay candidate #it{p}_{T} distribution for #it{R} = %2.2f and #it{p}_{T}^{fr} = %2.2f GeV/#it{c}",fConeSizes[icone],fPtFractions[ipt]);
3084 fhPtPtFracDecayIso[icone][ipt] = new TH1F(name, title,nptbins,ptmin,ptmax);
3085 fhPtPtFracDecayIso[icone][ipt]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3086 outputContainer->Add(fhPtPtFracDecayIso[icone][ipt]) ;
3088 snprintf(name, buffersize,"hPtSum_Decay_Cone_%d_Pt%d",icone,ipt);
3089 snprintf(title, buffersize,"Isolated decay candidate #it{p}_{T} distribution for #it{R} = %2.2f and #it{p}_{T}^{sum} = %2.2f GeV/#it{c}",fConeSizes[icone],fSumPtThresholds[ipt]);
3090 fhPtPtSumDecayIso[icone][ipt] = new TH1F(name, title,nptbins,ptmin,ptmax);//,nptsumbins,ptsummin,ptsummax);
3091 // fhPtPtSumDecayIso[icone]->SetYTitle("#Sigma #it{p}_{T} (GeV/#it{c})");
3092 fhPtPtSumDecayIso[icone][ipt]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3093 outputContainer->Add(fhPtPtSumDecayIso[icone][ipt]) ;
3095 snprintf(name, buffersize,"hPtSumDensity_Decay_Cone_%d_Pt%d",icone,ipt);
3096 snprintf(title, buffersize,"Isolated decay candidate #it{p}_{T} distribution for density in #it{R} = %2.2f and #it{p}_{T}^{sum} = %2.2f GeV/#it{c}",fConeSizes[icone],fSumPtThresholds[ipt]);
3097 fhPtSumDensityDecayIso[icone][ipt] = new TH1F(name, title,nptbins,ptmin,ptmax);//,nptsumbins,ptsummin,ptsummax);
3098 // fhPtPtSumDecayIso[icone]->SetYTitle("#Sigma #it{p}_{T} (GeV/#it{c})");
3099 fhPtSumDensityDecayIso[icone][ipt]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3100 outputContainer->Add(fhPtSumDensityDecayIso[icone][ipt]) ;
3102 snprintf(name, buffersize,"hPtFracPtSum_Decay_Cone_%d_Pt%d",icone,ipt);
3103 snprintf(title, buffersize,"Isolated decay candidate #it{p}_{T} distribution for PtFracPtSum in #it{R} = %2.2f and #it{p}_{T}^{fr} = %2.2f GeV/#it{c}",fConeSizes[icone],fPtFractions[ipt]);
3104 fhPtFracPtSumDecayIso[icone][ipt] = new TH1F(name, title,nptbins,ptmin,ptmax);//,nptsumbins,ptsummin,ptsummax);
3105 // fhPtPtSumDecayIso[icone]->SetYTitle("#Sigma #it{p}_{T} (GeV/#it{c})");
3106 fhPtFracPtSumDecayIso[icone][ipt]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3107 outputContainer->Add(fhPtFracPtSumDecayIso[icone][ipt]) ;
3110 snprintf(name, buffersize,"hEtaPhiPtThres_Decay_Cone_%d_Pt%d",icone,ipt);
3111 snprintf(title, buffersize,"Isolated decay candidate #eta:#phi distribution for #it{R} = %2.2f and #it{p}_{T}^{th} = %2.2f GeV/#it{c}",fConeSizes[icone],fPtThresholds[ipt]);
3112 fhEtaPhiPtThresDecayIso[icone][ipt] = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
3113 fhEtaPhiPtThresDecayIso[icone][ipt]->SetXTitle("#eta");
3114 fhEtaPhiPtThresDecayIso[icone][ipt]->SetYTitle("#phi");
3115 outputContainer->Add(fhEtaPhiPtThresDecayIso[icone][ipt]) ;
3117 snprintf(name, buffersize,"hEtaPhiPtFrac_Decay_Cone_%d_Pt%d",icone,ipt);
3118 snprintf(title, buffersize,"Isolated decay candidate #eta:#phi distribution for #it{R} = %2.2f and #it{p}_{T}^{fr} = %2.2f GeV/#it{c}",fConeSizes[icone],fPtFractions[ipt]);
3119 fhEtaPhiPtFracDecayIso[icone][ipt] = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
3120 fhEtaPhiPtFracDecayIso[icone][ipt]->SetXTitle("#eta");
3121 fhEtaPhiPtFracDecayIso[icone][ipt]->SetYTitle("#phi");
3122 outputContainer->Add(fhEtaPhiPtFracDecayIso[icone][ipt]) ;
3125 snprintf(name, buffersize,"hEtaPhiPtSum_Decay_Cone_%d_Pt%d",icone,ipt);
3126 snprintf(title, buffersize,"Isolated decay candidate #eta:#phi distribution for #it{R} = %2.2f and #it{p}_{T}^{sum} = %2.2f GeV/#it{c}",fConeSizes[icone],fSumPtThresholds[ipt]);
3127 fhEtaPhiPtSumDecayIso[icone][ipt] = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
3128 fhEtaPhiPtSumDecayIso[icone][ipt]->SetXTitle("#eta");
3129 fhEtaPhiPtSumDecayIso[icone][ipt]->SetYTitle("#phi");
3130 outputContainer->Add(fhEtaPhiPtSumDecayIso[icone][ipt]) ;
3132 snprintf(name, buffersize,"hEtaPhiSumDensity_Decay_Cone_%d_Pt%d",icone,ipt);
3133 snprintf(title, buffersize,"Isolated decay candidate #eta:#phi distribution for density #it{R} = %2.2f and #it{p}_{T}^{sum} = %2.2f GeV/#it{c}",fConeSizes[icone],fSumPtThresholds[ipt]);
3134 fhEtaPhiSumDensityDecayIso[icone][ipt] = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
3135 fhEtaPhiSumDensityDecayIso[icone][ipt]->SetXTitle("#eta");
3136 fhEtaPhiSumDensityDecayIso[icone][ipt]->SetYTitle("#phi");
3137 outputContainer->Add(fhEtaPhiSumDensityDecayIso[icone][ipt]) ;
3139 snprintf(name, buffersize,"hEtaPhiFracPtSum_Decay_Cone_%d_Pt%d",icone,ipt);
3140 snprintf(title, buffersize,"Isolated decay candidate #eta:#phi distribution for FracPtSum #it{R} = %2.2f and #it{p}_{T}^{fr} = %2.2f GeV/#it{c}",fConeSizes[icone],fPtFractions[ipt]);
3141 fhEtaPhiFracPtSumDecayIso[icone][ipt] = new TH2F(name, title,netabins,etamin,etamax,nphibins,phimin,phimax);
3142 fhEtaPhiFracPtSumDecayIso[icone][ipt]->SetXTitle("#eta");
3143 fhEtaPhiFracPtSumDecayIso[icone][ipt]->SetYTitle("#phi");
3144 outputContainer->Add(fhEtaPhiFracPtSumDecayIso[icone][ipt]) ;
3150 for(Int_t imc = 0; imc < fgkNmcTypes; imc++)
3152 snprintf(name , buffersize,"hPtThreshMC%s_Cone_%d_Pt%d",mcPartName[imc].Data(),icone,ipt);
3153 snprintf(title, buffersize,"Isolated %s #it{p}_{T} for #it{R}=%2.2f and #it{p}_{T}^{th}=%2.2f",
3154 mcPartType[imc].Data(),fConeSizes[icone], fPtThresholds[ipt]);
3155 fhPtThresIsolatedMC[imc][icone][ipt] = new TH1F(name, title,nptbins,ptmin,ptmax);
3156 fhPtThresIsolatedMC[imc][icone][ipt]->SetYTitle("#it{counts}");
3157 fhPtThresIsolatedMC[imc][icone][ipt]->SetXTitle("#it{p}_{T}(GeV/#it{c})");
3158 outputContainer->Add(fhPtThresIsolatedMC[imc][icone][ipt]) ;
3161 snprintf(name , buffersize,"hPtFracMC%s_Cone_%d_Pt%d",mcPartName[imc].Data(),icone,ipt);
3162 snprintf(title, buffersize,"Isolated %s #it{p}_{T} for #it{R}=%2.2f and #Sigma #it{p}_{T}^{in cone}/#it{p}_{T}^{trig}=%2.2f",
3163 mcPartType[imc].Data(),fConeSizes[icone], fPtFractions[ipt]);
3164 fhPtFracIsolatedMC[imc][icone][ipt] = new TH1F(name, title,nptbins,ptmin,ptmax);
3165 fhPtFracIsolatedMC[imc][icone][ipt]->SetYTitle("#it{counts}");
3166 fhPtFracIsolatedMC[imc][icone][ipt]->SetXTitle("#it{p}_{T}(GeV/#it{c})");
3167 outputContainer->Add(fhPtFracIsolatedMC[imc][icone][ipt]) ;
3169 snprintf(name , buffersize,"hSumPtMC%s_Cone_%d_Pt%d",mcPartName[imc].Data(),icone,ipt);
3170 snprintf(title, buffersize,"Isolated %s #it{p}_{T} for #it{R}=%2.2f and #Sigma #it{p}_{T}^{in cone}=%2.2f",
3171 mcPartType[imc].Data(),fConeSizes[icone], fSumPtThresholds[ipt]);
3172 fhSumPtIsolatedMC[imc][icone][ipt] = new TH1F(name, title,nptbins,ptmin,ptmax);
3173 fhSumPtIsolatedMC[imc][icone][ipt]->SetYTitle("#it{counts}");
3174 fhSumPtIsolatedMC[imc][icone][ipt]->SetXTitle("#it{p}_{T}(GeV/#it{c})");
3175 outputContainer->Add(fhSumPtIsolatedMC[imc][icone][ipt]) ;
3182 if(fFillPileUpHistograms)
3184 for (Int_t i = 0; i < 7 ; i++)
3186 fhEIsoPileUp[i] = new TH1F(Form("hEPileUp%s",pileUpName[i].Data()),
3187 Form("Number of isolated particles vs E, %s, pile-up event by %s",parTitle.Data(),pileUpName[i].Data()),
3188 nptbins,ptmin,ptmax);
3189 fhEIsoPileUp[i]->SetYTitle("d#it{N} / d#it{E}");
3190 fhEIsoPileUp[i]->SetXTitle("#it{E} (GeV)");
3191 outputContainer->Add(fhEIsoPileUp[i]) ;
3193 fhPtIsoPileUp[i] = new TH1F(Form("hPtPileUp%s",pileUpName[i].Data()),
3194 Form("Number of isolated particles vs #it{p}_{T}, %s, pile-up event by %s",parTitle.Data(),pileUpName[i].Data()),
3195 nptbins,ptmin,ptmax);
3196 fhPtIsoPileUp[i]->SetYTitle("d#it{N} / #it{p}_{T}");
3197 fhPtIsoPileUp[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3198 outputContainer->Add(fhPtIsoPileUp[i]) ;
3200 fhENoIsoPileUp[i] = new TH1F(Form("hENoIsoPileUp%s",pileUpName[i].Data()),
3201 Form("Number of not isolated particles vs E, %s, pile-up event by %s",parTitle.Data(),pileUpName[i].Data()),
3202 nptbins,ptmin,ptmax);
3203 fhENoIsoPileUp[i]->SetYTitle("d#it{N} / dE");
3204 fhENoIsoPileUp[i]->SetXTitle("#it{E} (GeV)");
3205 outputContainer->Add(fhENoIsoPileUp[i]) ;
3207 fhPtNoIsoPileUp[i] = new TH1F(Form("hPtNoIsoPileUp%s",pileUpName[i].Data()),
3208 Form("Number of not isolated particles vs #it{p}_{T}, %s, pile-up event by %s",parTitle.Data(),pileUpName[i].Data()),
3209 nptbins,ptmin,ptmax);
3210 fhPtNoIsoPileUp[i]->SetYTitle("d#it{N} / #it{p}_{T}");
3211 fhPtNoIsoPileUp[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3212 outputContainer->Add(fhPtNoIsoPileUp[i]) ;
3215 fhTimeENoCut = new TH2F ("hTimeE_NoCut","time of cluster vs E of clusters, no cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
3216 fhTimeENoCut->SetXTitle("#it{E} (GeV)");
3217 fhTimeENoCut->SetYTitle("#it{time} (ns)");
3218 outputContainer->Add(fhTimeENoCut);
3220 fhTimeESPD = new TH2F ("hTimeE_SPD","time of cluster vs E of clusters, SPD cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
3221 fhTimeESPD->SetXTitle("#it{E} (GeV)");
3222 fhTimeESPD->SetYTitle("#it{time} (ns)");
3223 outputContainer->Add(fhTimeESPD);
3225 fhTimeESPDMulti = new TH2F ("hTimeE_SPDMulti","time of cluster vs E of clusters, SPD multi cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
3226 fhTimeESPDMulti->SetXTitle("#it{E} (GeV)");
3227 fhTimeESPDMulti->SetYTitle("#it{time} (ns)");
3228 outputContainer->Add(fhTimeESPDMulti);
3230 fhTimeNPileUpVertSPD = new TH2F ("hTime_NPileUpVertSPD","time of cluster vs N pile-up SPD vertex", ntimebins,timemin,timemax,50,0,50);
3231 fhTimeNPileUpVertSPD->SetYTitle("# vertex ");
3232 fhTimeNPileUpVertSPD->SetXTitle("#it{time} (ns)");
3233 outputContainer->Add(fhTimeNPileUpVertSPD);
3235 fhTimeNPileUpVertTrack = new TH2F ("hTime_NPileUpVertTracks","time of cluster vs N pile-up Tracks vertex", ntimebins,timemin,timemax, 50,0,50 );
3236 fhTimeNPileUpVertTrack->SetYTitle("# vertex ");
3237 fhTimeNPileUpVertTrack->SetXTitle("#it{time} (ns)");
3238 outputContainer->Add(fhTimeNPileUpVertTrack);
3240 fhTimeNPileUpVertContributors = new TH2F ("hTime_NPileUpVertContributors","time of cluster vs N constributors to pile-up SPD vertex", ntimebins,timemin,timemax,50,0,50);
3241 fhTimeNPileUpVertContributors->SetYTitle("# vertex ");
3242 fhTimeNPileUpVertContributors->SetXTitle("#it{time} (ns)");
3243 outputContainer->Add(fhTimeNPileUpVertContributors);
3245 fhTimePileUpMainVertexZDistance = new TH2F ("hTime_PileUpMainVertexZDistance","time of cluster vs distance in Z pile-up SPD vertex - main SPD vertex",ntimebins,timemin,timemax,100,0,50);
3246 fhTimePileUpMainVertexZDistance->SetYTitle("distance #it{z} (cm) ");
3247 fhTimePileUpMainVertexZDistance->SetXTitle("#it{time} (ns)");
3248 outputContainer->Add(fhTimePileUpMainVertexZDistance);
3250 fhTimePileUpMainVertexZDiamond = new TH2F ("hTime_PileUpMainVertexZDiamond","time of cluster vs distance in Z pile-up SPD vertex - z diamond",ntimebins,timemin,timemax,100,0,50);
3251 fhTimePileUpMainVertexZDiamond->SetYTitle("diamond distance #it{z} (cm) ");
3252 fhTimePileUpMainVertexZDiamond->SetXTitle("#it{time} (ns)");
3253 outputContainer->Add(fhTimePileUpMainVertexZDiamond);
3256 return outputContainer ;
3260 //____________________________________________________
3261 Int_t AliAnaParticleIsolation::GetMCIndex(Int_t mcTag)
3263 // Histogram index depending on origin of candidate
3265 if(!IsDataMC()) return -1;
3267 if (GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPrompt))
3271 else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCFragmentation) ||
3272 GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCISR))
3276 else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0))
3280 else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0Decay))
3284 else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEtaDecay))
3288 else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCOtherDecay))
3290 return kmcOtherDecay;
3292 else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron))
3296 else // anything else
3298 // careful can contain also other decays, to be checked.
3303 //__________________________________
3304 void AliAnaParticleIsolation::Init()
3306 // Do some checks and init stuff
3308 // In case of several cone and thresholds analysis, open the cuts for the filling of the
3309 // track and cluster reference arrays in cone when done in the MakeAnalysisFillAOD().
3310 // The different cones, thresholds are tested for this list of tracks, clusters.
3313 printf("AliAnaParticleIsolation::Init() - Open default isolation cuts for multiple Isolation analysis\n");
3314 GetIsolationCut()->SetPtThreshold(100);
3315 GetIsolationCut()->SetPtFraction(100);
3316 GetIsolationCut()->SetConeSize(1);
3319 if(!GetReader()->IsCTSSwitchedOn() && GetIsolationCut()->GetParticleTypeInCone()!=AliIsolationCut::kOnlyNeutral)
3320 AliFatal("STOP!: You want to use CTS tracks in analysis but not read!! \n!!Check the configuration file!!\n");
3324 //____________________________________________
3325 void AliAnaParticleIsolation::InitParameters()
3328 //Initialize the parameters of the analysis.
3329 SetInputAODName("PWG4Particle");
3330 SetAODObjArrayName("IsolationCone");
3331 AddToHistogramsName("AnaIsolation_");
3333 fCalorimeter = "EMCAL" ;
3334 fIsoDetector = "EMCAL" ;
3336 fReMakeIC = kFALSE ;
3337 fMakeSeveralIC = kFALSE ;
3339 fLeadingOnly = kTRUE;
3340 fCheckLeadingWithNeutralClusters = kTRUE;
3343 fDecayBits[0] = AliNeutralMesonSelection::kPi0;
3344 fDecayBits[1] = AliNeutralMesonSelection::kEta;
3345 fDecayBits[2] = AliNeutralMesonSelection::kPi0Side;
3346 fDecayBits[3] = AliNeutralMesonSelection::kEtaSide;
3349 fBkgBinLimit[ 0] = 00.0; fBkgBinLimit[ 1] = 00.2; fBkgBinLimit[ 2] = 00.3; fBkgBinLimit[ 3] = 00.4; fBkgBinLimit[ 4] = 00.5;
3350 fBkgBinLimit[ 5] = 01.0; fBkgBinLimit[ 6] = 01.5; fBkgBinLimit[ 7] = 02.0; fBkgBinLimit[ 8] = 03.0; fBkgBinLimit[ 9] = 05.0;
3351 fBkgBinLimit[10] = 10.0; fBkgBinLimit[11] = 100.;
3352 for(Int_t ibin = 12; ibin < 20; ibin++) fBkgBinLimit[ibin] = 00.0;
3354 //----------- Several IC-----------------
3357 fConeSizes [0] = 0.1; fConeSizes [1] = 0.2; fConeSizes [2] = 0.3; fConeSizes [3] = 0.4; fConeSizes [4] = 0.5;
3358 fPtThresholds [0] = 1.; fPtThresholds [1] = 2.; fPtThresholds [2] = 3.; fPtThresholds [3] = 4.; fPtThresholds [4] = 5.;
3359 fPtFractions [0] = 0.05; fPtFractions [1] = 0.075; fPtFractions [2] = 0.1; fPtFractions [3] = 1.25; fPtFractions [4] = 1.5;
3360 fSumPtThresholds[0] = 1.; fSumPtThresholds[1] = 2.; fSumPtThresholds[2] = 3.; fSumPtThresholds[3] = 4.; fSumPtThresholds[4] = 5.;
3364 //_________________________________________________________________________________________
3365 Bool_t AliAnaParticleIsolation::IsTriggerTheNearSideEventLeadingParticle(Int_t & idLeading)
3367 // Check if the what of the selected isolation candidates is leading particle in the same hemisphere
3368 // comparing with all the candidates, all the tracks or all the clusters.
3370 Double_t ptTrig = GetMinPt();
3371 Double_t phiTrig = 0 ;
3373 AliAODPWG4ParticleCorrelation* pLeading = 0;
3375 // Loop on stored AOD particles, find leading trigger on the selected list, with at least min pT.
3377 for(Int_t iaod = 0; iaod < GetInputAODBranch()->GetEntriesFast() ; iaod++)
3379 AliAODPWG4ParticleCorrelation* particle = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
3380 particle->SetLeadingParticle(kFALSE); // set it later
3382 // Vertex cut in case of mixing
3385 Int_t check = CheckMixedEventVertex(particle->GetCaloLabel(0), particle->GetTrackLabel(0));
3386 if(check == 0) continue;
3387 if(check == -1) return kFALSE; // not sure if it is correct.
3390 //check if it is low pt trigger particle
3391 if((particle->Pt() < GetIsolationCut()->GetPtThreshold() ||
3392 particle->Pt() < GetIsolationCut()->GetSumPtThreshold()) &&
3395 continue ; //trigger should not come from underlying event
3398 // find the leading particles with highest momentum
3399 if (particle->Pt() > ptTrig)
3401 ptTrig = particle->Pt() ;
3402 phiTrig = particle->Phi();
3404 pLeading = particle ;
3406 }// finish search of leading trigger particle on the AOD branch.
3408 if(index < 0) return kFALSE;
3412 //printf("AOD leading pT %2.2f, ID %d\n", pLeading->Pt(),pLeading->GetCaloLabel(0));
3414 if(phiTrig < 0 ) phiTrig += TMath::TwoPi();
3416 // Compare if it is the leading of all tracks
3419 for(Int_t ipr = 0;ipr < GetCTSTracks()->GetEntriesFast() ; ipr ++ )
3421 AliVTrack * track = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
3423 if(track->GetID() == pLeading->GetTrackLabel(0) || track->GetID() == pLeading->GetTrackLabel(1) ||
3424 track->GetID() == pLeading->GetTrackLabel(2) || track->GetID() == pLeading->GetTrackLabel(3) ) continue ;
3426 Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
3427 p3.SetXYZ(mom[0],mom[1],mom[2]);
3428 Float_t pt = p3.Pt();
3429 Float_t phi = p3.Phi() ;
3430 if(phi < 0) phi+=TMath::TwoPi();
3432 //skip this event if near side associated particle pt larger than trigger
3434 Float_t deltaPhi = phiTrig-phi;
3436 // Calculate deltaPhi shift so that for the particles on the opposite side
3437 // it is defined between 90 and 270 degrees
3438 // Shift [-360,-90] to [0, 270]
3439 // and [270,360] to [-90,0]
3440 if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
3441 if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
3443 if(pt > ptTrig && deltaPhi < TMath::PiOver2()) return kFALSE;
3447 // Compare if it is leading of all calorimeter clusters
3449 if(fCheckLeadingWithNeutralClusters)
3451 // Select the calorimeter cluster list
3452 TObjArray * nePl = 0x0;
3453 if (pLeading->GetDetector() == "PHOS" )
3454 nePl = GetPHOSClusters();
3456 nePl = GetEMCALClusters();
3458 if(!nePl) return kTRUE; // Do the selection just with the tracks if no calorimeter is available.
3461 for(Int_t ipr = 0;ipr < nePl->GetEntriesFast() ; ipr ++ )
3463 AliVCluster * cluster = (AliVCluster *) (nePl->At(ipr)) ;
3465 if(cluster->GetID() == pLeading->GetCaloLabel(0) || cluster->GetID() == pLeading->GetCaloLabel(1) ) continue ;
3467 cluster->GetMomentum(lv,GetVertex(0));
3469 Float_t pt = lv.Pt();
3470 Float_t phi = lv.Phi() ;
3471 if(phi < 0) phi+=TMath::TwoPi();
3473 if(IsTrackMatched(cluster,GetReader()->GetInputEvent())) continue ; // avoid charged clusters, already covered by tracks, or cluster merging with track.
3475 // skip this event if near side associated particle pt larger than trigger
3476 // not really needed for calorimeter, unless DCal is included
3478 Float_t deltaPhi = phiTrig-phi;
3479 if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
3480 if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
3482 if(pt > ptTrig && deltaPhi < TMath::PiOver2()) return kFALSE ;
3485 } // check neutral clusters
3488 pLeading->SetLeadingParticle(kTRUE);
3490 if( GetDebug() > 1 )
3491 printf("AliAnaParticleIsolation::IsTriggerTheNearSideEventLeadingParticle() - Particle AOD with index %d is leading with pT %2.2f\n",
3492 idLeading, pLeading->Pt());
3498 //__________________________________________________
3499 void AliAnaParticleIsolation::MakeAnalysisFillAOD()
3501 // Do analysis and fill aods
3502 // Search for the isolated photon in fCalorimeter with GetMinPt() < pt < GetMaxPt()
3503 // and if the particle is leading in the near side (if requested)
3505 if(!GetInputAODBranch())
3506 AliFatal(Form("No input particles in AOD with name branch < %s >, STOP",GetInputAODName().Data()));
3508 if(strcmp(GetInputAODBranch()->GetClass()->GetName(), "AliAODPWG4ParticleCorrelation"))
3509 AliFatal(Form("Wrong type of AOD object, change AOD class name in input AOD: It should be <AliAODPWG4ParticleCorrelation> and not <%s> \n",GetInputAODBranch()->GetClass()->GetName()));
3511 Int_t n = 0, nfrac = 0;
3512 Bool_t isolated = kFALSE ;
3513 Float_t coneptsum = 0, coneptlead = 0;
3514 TObjArray * pl = 0x0; ;
3516 //Select the calorimeter for candidate isolation with neutral particles
3517 if (fCalorimeter == "PHOS" )
3518 pl = GetPHOSClusters();
3519 else if (fCalorimeter == "EMCAL")
3520 pl = GetEMCALClusters();
3522 //Loop on AOD branch, filled previously in AliAnaPhoton, find leading particle to do isolation only with it
3523 TLorentzVector mom ;
3524 Int_t idLeading = -1 ;
3526 Int_t naod = GetInputAODBranch()->GetEntriesFast();
3529 printf("AliAnaParticleIsolation::MakeAnalysisFillAOD() - Input aod branch entries %d\n", naod);
3531 if(IsLeadingOnlyOn())
3533 Bool_t leading = IsTriggerTheNearSideEventLeadingParticle(idLeading);
3536 if(GetDebug() > 1) printf("AliAnaParticleIsolation::MakeAnalysisFillAOD() - Not leading; End fill AODs \n");
3539 iaod0 = idLeading ; // first entry in particle loop
3540 naod = idLeading+1; // last entry in particle loop
3543 // Check isolation of list of candidate particles or leading particle
3545 for(Int_t iaod = iaod0; iaod < naod; iaod++)
3547 AliAODPWG4ParticleCorrelation * aodinput = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
3549 // Check isolation only of clusters in fiducial region
3551 if(IsFiducialCutOn())
3553 Bool_t in = GetFiducialCut()->IsInFiducialCut(*aodinput->Momentum(),aodinput->GetDetector()) ;
3554 if(! in ) continue ;
3557 //If too small or too large pt, skip
3558 Float_t pt = aodinput->Pt();
3559 if(pt < GetMinPt() || pt > GetMaxPt() ) continue ;
3562 //check if it is low pt trigger particle
3563 if( ( pt < GetIsolationCut()->GetPtThreshold() || pt < GetIsolationCut()->GetSumPtThreshold() ) &&
3566 continue ; //trigger should not come from underlying event
3569 //After cuts, study isolation
3570 n=0; nfrac = 0; isolated = kFALSE; coneptsum = 0; coneptlead = 0;
3571 GetIsolationCut()->MakeIsolationCut(GetCTSTracks(),pl,
3572 GetReader(), GetCaloPID(),
3573 kTRUE, aodinput, GetAODObjArrayName(),
3574 n,nfrac,coneptsum,coneptlead,isolated);
3576 if(!fMakeSeveralIC) aodinput->SetIsolated(isolated);
3580 if(isolated) printf("AliAnaParticleIsolation::MakeAnalysisFillAOD() : Particle %d IS ISOLATED \n",iaod);
3581 printf("AliAnaParticleIsolation::MakeAnalysisFillAOD() - End fill AODs \n");
3584 } // particle isolation loop
3588 //_________________________________________________________
3589 void AliAnaParticleIsolation::MakeAnalysisFillHistograms()
3591 // Do analysis and fill histograms
3593 // In case of simulated data, fill acceptance histograms
3594 // Not ready for multiple case analysis.
3595 if(IsDataMC() && !fMakeSeveralIC) FillAcceptanceHistograms();
3597 //Loop on stored AOD
3598 Int_t naod = GetInputAODBranch()->GetEntriesFast();
3600 printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Histo aod branch entries %d\n", naod);
3602 for(Int_t iaod = 0; iaod < naod ; iaod++)
3604 AliAODPWG4ParticleCorrelation* aod = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
3606 if(IsLeadingOnlyOn() && !aod->IsLeadingParticle()) continue; // Try to isolate only leading cluster or track
3608 // Check isolation only of clusters in fiducial region
3609 if(IsFiducialCutOn())
3611 Bool_t in = GetFiducialCut()->IsInFiducialCut(*aod->Momentum(),aod->GetDetector()) ;
3612 if(! in ) continue ;
3615 Float_t pt = aod->Pt();
3617 //If too small or too large pt, skip
3618 if(pt < GetMinPt() || pt > GetMaxPt() ) continue ;
3620 Int_t mcTag = aod->GetTag() ;
3621 Int_t mcIndex = GetMCIndex(mcTag);
3623 // --- In case of redoing isolation from delta AOD ----
3624 // Not standard case, not used since its implementation
3627 //Analysis of multiple IC at same time
3628 MakeSeveralICAnalysis(aod,mcIndex);
3632 // --- In case of redoing isolation multiple cuts ----
3636 //In case a more strict IC is needed in the produced AOD
3637 Bool_t isolated = kFALSE;
3638 Int_t n = 0, nfrac = 0;
3639 Float_t coneptsum = 0, coneptlead = 0;
3641 //Recover reference arrays with clusters and tracks
3642 TObjArray * refclusters = aod->GetObjArray(GetAODObjArrayName()+"Clusters");
3643 TObjArray * reftracks = aod->GetObjArray(GetAODObjArrayName()+"Tracks");
3645 GetIsolationCut()->MakeIsolationCut(reftracks, refclusters,
3646 GetReader(), GetCaloPID(),
3648 n,nfrac,coneptsum,coneptlead,isolated);
3651 Bool_t isolated = aod->IsIsolated();
3652 Float_t energy = aod->E();
3653 Float_t phi = aod->Phi();
3654 Float_t eta = aod->Eta();
3657 if(fFillTaggedDecayHistograms)
3659 decayTag = aod->GetBtag(); // temporary
3660 if(decayTag < 0) decayTag = 0; // temporary
3664 printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - pt %1.1f, eta %1.1f, phi %1.1f, Isolated %d\n",
3665 pt, eta, phi, isolated);
3667 //---------------------------------------------------------------
3668 // Fill pt/sum pT distribution of particles in cone or in UE band
3669 //---------------------------------------------------------------
3671 Float_t coneptLeadCluster= 0;
3672 Float_t coneptLeadTrack = 0;
3673 Float_t coneptsumCluster = 0;
3674 Float_t coneptsumTrack = 0;
3675 Float_t coneptsumCell = 0;
3676 Float_t etaBandptsumClusterNorm = 0;
3677 Float_t etaBandptsumTrackNorm = 0;
3679 CalculateTrackSignalInCone (aod,coneptsumTrack , coneptLeadTrack );
3680 CalculateCaloSignalInCone (aod,coneptsumCluster, coneptLeadCluster);
3681 if(fFillCellHistograms)
3682 CalculateCaloCellSignalInCone(aod,coneptsumCell );
3684 if(GetIsolationCut()->GetParticleTypeInCone()==AliIsolationCut::kNeutralAndCharged)
3686 fhConeSumPtClustervsTrack ->Fill(coneptsumCluster, coneptsumTrack );
3687 fhConePtLeadClustervsTrack->Fill(coneptLeadCluster,coneptLeadTrack);
3689 if(coneptsumTrack > 0) fhConeSumPtClusterTrackFrac ->Fill(pt, coneptsumCluster /coneptsumTrack );
3690 if(coneptLeadTrack > 0) fhConePtLeadClusterTrackFrac->Fill(pt, coneptLeadCluster/coneptLeadTrack);
3692 if(fFillCellHistograms)
3694 fhConeSumPtCellvsTrack ->Fill(coneptsumCell, coneptsumTrack);
3695 fhConeSumPtCellTrack ->Fill(pt, coneptsumTrack+coneptsumCell);
3696 fhConeSumPtCellTrackTrigEtaPhi->Fill(eta,phi,coneptsumTrack+coneptsumCell);
3700 fhConeSumPt ->Fill(pt, coneptsumTrack+coneptsumCluster);
3701 fhConeSumPtTrigEtaPhi ->Fill(eta,phi,coneptsumTrack+coneptsumCluster);
3703 Float_t coneptLead = coneptLeadTrack;
3704 if(coneptLeadCluster > coneptLeadTrack) coneptLead = coneptLeadCluster;
3705 fhConePtLead->Fill(pt, coneptLead );
3708 printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Particle %d Energy Sum in Isolation Cone %2.2f, Leading pT in cone %2.2f\n",
3709 iaod, coneptsumTrack+coneptsumCluster, coneptLead);
3711 //normalize phi/eta band per area unit
3712 if(fFillUEBandSubtractHistograms)
3713 CalculateNormalizeUEBandPerUnitArea(aod, coneptsumCluster, coneptsumCell, coneptsumTrack, etaBandptsumTrackNorm, etaBandptsumClusterNorm) ;
3715 // printf("Histograms analysis : cluster pt = %f, etaBandTrack = %f, etaBandCluster = %f, isolation = %d\n",aod->Pt(),etaBandptsumTrackNorm,etaBandptsumClusterNorm,aod->IsIsolated());
3717 //---------------------------------------------------------------
3718 // Fill Shower shape and track matching histograms
3719 //---------------------------------------------------------------
3721 FillTrackMatchingShowerShapeControlHistograms(aod, coneptsumTrack+coneptsumCluster, coneptLead, mcIndex);
3723 //---------------------------------------------------------------
3724 // Isolated/ Non isolated histograms
3725 //---------------------------------------------------------------
3730 printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Particle %d ISOLATED: fill histograms\n", iaod);
3732 fhEIso ->Fill(energy);
3734 fhPhiIso ->Fill(pt,phi);
3735 fhEtaIso ->Fill(pt,eta);
3736 fhEtaPhiIso ->Fill(eta,phi);
3740 // For histograms in arrays, index in the array, corresponding to any particle origin
3741 if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton))
3743 fhPtIsoMC [kmcPhoton]->Fill(pt);
3744 fhPhiIsoMC[kmcPhoton]->Fill(pt,phi);
3745 fhEtaIsoMC[kmcPhoton]->Fill(pt,eta);
3748 if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCDecayPairLost) && mcIndex==kmcPi0Decay )
3750 fhPtIsoMC [kmcPi0DecayLostPair]->Fill(pt);
3751 fhPhiIsoMC[kmcPi0DecayLostPair]->Fill(pt,phi);
3752 fhEtaIsoMC[kmcPi0DecayLostPair]->Fill(pt,eta);
3755 fhPtIsoMC [mcIndex]->Fill(pt);
3756 fhPhiIsoMC[mcIndex]->Fill(pt,phi);
3757 fhEtaIsoMC[mcIndex]->Fill(pt,eta);
3758 }//Histograms with MC
3760 // Candidates tagged as decay in another analysis (AliAnaPi0EbE)
3761 if(fFillTaggedDecayHistograms && decayTag > 0)
3763 for(Int_t ibit = 0; ibit < fNDecayBits; ibit++)
3765 if(GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[ibit]))
3767 fhPtDecayIso [ibit]->Fill(pt);
3768 fhEtaPhiDecayIso [ibit]->Fill(eta,phi);
3772 if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton))
3773 fhPtDecayIsoMC[ibit][kmcPhoton]->Fill(pt);
3775 if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCDecayPairLost) && mcIndex==kmcPi0Decay )
3776 fhPtDecayIsoMC[ibit][kmcPi0DecayLostPair]->Fill(pt);
3778 fhPtDecayIsoMC[ibit][mcIndex]->Fill(pt);
3782 } // decay histograms
3784 if(fFillNLMHistograms)
3785 fhPtNLocMaxIso ->Fill(pt,aod->GetFiducialArea()) ; // remember to change method name
3787 if(fFillHighMultHistograms)
3789 fhPtCentralityIso ->Fill(pt,GetEventCentrality()) ;
3790 fhPtEventPlaneIso ->Fill(pt,GetEventPlaneAngle() ) ;
3793 if(fFillPileUpHistograms)
3795 if(GetReader()->IsPileUpFromSPD()) { fhEIsoPileUp[0] ->Fill(energy) ; fhPtIsoPileUp[0]->Fill(pt) ; }
3796 if(GetReader()->IsPileUpFromEMCal()) { fhEIsoPileUp[1] ->Fill(energy) ; fhPtIsoPileUp[1]->Fill(pt) ; }
3797 if(GetReader()->IsPileUpFromSPDOrEMCal()) { fhEIsoPileUp[2] ->Fill(energy) ; fhPtIsoPileUp[2]->Fill(pt) ; }
3798 if(GetReader()->IsPileUpFromSPDAndEMCal()) { fhEIsoPileUp[3] ->Fill(energy) ; fhPtIsoPileUp[3]->Fill(pt) ; }
3799 if(GetReader()->IsPileUpFromSPDAndNotEMCal()) { fhEIsoPileUp[4] ->Fill(energy) ; fhPtIsoPileUp[4]->Fill(pt) ; }
3800 if(GetReader()->IsPileUpFromEMCalAndNotSPD()) { fhEIsoPileUp[5] ->Fill(energy) ; fhPtIsoPileUp[5]->Fill(pt) ; }
3801 if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) { fhEIsoPileUp[6] ->Fill(energy) ; fhPtIsoPileUp[6]->Fill(pt) ; }
3803 // Fill histograms to undertand pile-up before other cuts applied
3804 // Remember to relax time cuts in the reader
3805 FillPileUpHistograms(aod->GetCaloLabel(0));
3808 }//Isolated histograms
3809 else // NON isolated
3812 printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Particle %d NOT ISOLATED, fill histograms\n", iaod);
3814 fhENoIso ->Fill(energy);
3815 fhPtNoIso ->Fill(pt);
3816 fhEtaPhiNoIso ->Fill(eta,phi);
3820 if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton) )
3821 fhPtNoIsoMC[kmcPhoton]->Fill(pt);
3823 if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCDecayPairLost) && mcIndex==kmcPi0Decay )
3824 fhPtNoIsoMC[kmcPi0DecayLostPair]->Fill(pt);
3826 fhPtNoIsoMC[mcIndex]->Fill(pt);
3829 // Candidates tagged as decay in another analysis (AliAnaPi0EbE)
3830 if(fFillTaggedDecayHistograms && decayTag > 0)
3832 for(Int_t ibit = 0; ibit < fNDecayBits; ibit++)
3834 if(GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[ibit]))
3836 fhPtDecayNoIso[ibit] ->Fill(pt);
3837 fhEtaPhiDecayNoIso[ibit]->Fill(eta,phi);
3841 if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton))
3842 fhPtDecayNoIsoMC[ibit][kmcPhoton]->Fill(pt);
3844 if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCDecayPairLost) && mcIndex==kmcPi0Decay )
3845 fhPtDecayNoIsoMC[ibit][kmcPi0DecayLostPair]->Fill(pt);
3847 fhPtDecayNoIsoMC[ibit][mcIndex]->Fill(pt);
3851 } // decay histograms
3853 if(fFillNLMHistograms)
3854 fhPtNLocMaxNoIso ->Fill(pt,aod->GetFiducialArea()); // remember to change method name
3856 if(fFillPileUpHistograms)
3858 if(GetReader()->IsPileUpFromSPD()) { fhENoIsoPileUp[0] ->Fill(energy) ; fhPtNoIsoPileUp[0]->Fill(pt) ; }
3859 if(GetReader()->IsPileUpFromEMCal()) { fhENoIsoPileUp[1] ->Fill(energy) ; fhPtNoIsoPileUp[1]->Fill(pt) ; }
3860 if(GetReader()->IsPileUpFromSPDOrEMCal()) { fhENoIsoPileUp[2] ->Fill(energy) ; fhPtNoIsoPileUp[2]->Fill(pt) ; }
3861 if(GetReader()->IsPileUpFromSPDAndEMCal()) { fhENoIsoPileUp[3] ->Fill(energy) ; fhPtNoIsoPileUp[3]->Fill(pt) ; }
3862 if(GetReader()->IsPileUpFromSPDAndNotEMCal()) { fhENoIsoPileUp[4] ->Fill(energy) ; fhPtNoIsoPileUp[4]->Fill(pt) ; }
3863 if(GetReader()->IsPileUpFromEMCalAndNotSPD()) { fhENoIsoPileUp[5] ->Fill(energy) ; fhPtNoIsoPileUp[5]->Fill(pt) ; }
3864 if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) { fhENoIsoPileUp[6] ->Fill(energy) ; fhPtNoIsoPileUp[6]->Fill(pt) ; }
3871 //______________________________________________________
3872 void AliAnaParticleIsolation::FillAcceptanceHistograms()
3874 // Fill acceptance histograms if MC data is available
3875 // Only when particle is in the calorimeter. Rethink if CTS is used.
3877 if(GetDebug() > 0) printf("AliAnaParticleIsolation::FillAcceptanceHistograms() - Start \n");
3879 //Double_t photonY = -100 ;
3880 Double_t photonE = -1 ;
3881 Double_t photonPt = -1 ;
3882 Double_t photonPhi = 100 ;
3883 Double_t photonEta = -1 ;
3891 TParticle * primStack = 0;
3892 AliAODMCParticle * primAOD = 0;
3895 // Calorimeter cluster merging angle
3896 // angle smaller than 3 cells 6 cm (0.014) in EMCal, 2.2 cm in PHOS (0.014*(2.2/6))
3897 Float_t overlapAngle = 0;
3898 Float_t minECalo = 0;
3899 if (fCalorimeter=="EMCAL")
3901 overlapAngle = 3*0.014 ;
3902 minECalo = GetReader()->GetEMCALEMin();
3904 else if (fCalorimeter=="PHOS" )
3906 overlapAngle = 3*0.00382;
3907 minECalo = GetReader()->GetPHOSEMin();
3910 // Get the ESD MC particles container
3911 AliStack * stack = 0;
3912 if( GetReader()->ReadStack() )
3914 stack = GetMCStack();
3916 nprim = stack->GetNtrack();
3919 // Get the AOD MC particles container
3920 TClonesArray * mcparticles = 0;
3921 if( GetReader()->ReadAODMCParticles() )
3923 mcparticles = GetReader()->GetAODMCParticles();
3924 if( !mcparticles ) return;
3925 nprim = mcparticles->GetEntriesFast();
3928 for(Int_t i=0 ; i < nprim; i++)
3930 if(GetReader()->AcceptOnlyHIJINGLabels() && !GetReader()->IsHIJINGLabel(i)) continue ;
3932 if(GetReader()->ReadStack())
3934 primStack = stack->Particle(i) ;
3937 printf("AliAnaParticleIsolation::FillAcceptanceHistograms() - ESD primaries pointer not available!!\n");
3941 pdg = primStack->GetPdgCode();
3942 status = primStack->GetStatusCode();
3944 if(primStack->Energy() == TMath::Abs(primStack->Pz())) continue ; //Protection against floating point exception
3946 //printf("i %d, %s %d %s %d \n",i, stack->Particle(i)->GetName(), stack->Particle(i)->GetPdgCode(),
3947 // prim->GetName(), prim->GetPdgCode());
3949 //photonY = 0.5*TMath::Log((prim->Energy()-prim->Pz())/(prim->Energy()+prim->Pz())) ;
3952 primStack->Momentum(lv);
3957 primAOD = (AliAODMCParticle *) mcparticles->At(i);
3960 printf("AliAnaParticleIsolation::FillAcceptanceHistograms() - AOD primaries pointer not available!!\n");
3964 pdg = primAOD->GetPdgCode();
3965 status = primAOD->GetStatus();
3967 if(primAOD->E() == TMath::Abs(primAOD->Pz())) continue ; //Protection against floating point exception
3969 //photonY = 0.5*TMath::Log((prim->Energy()-prim->Pz())/(prim->Energy()+prim->Pz())) ;
3972 lv.SetPxPyPzE(primAOD->Px(),primAOD->Py(),primAOD->Pz(),primAOD->E());
3975 // Select only photons in the final state
3976 if(pdg != 22 && pdg!=111) continue ;
3978 // Consider only final state particles, but this depends on generator,
3979 // status 1 is the usual one, in case of not being ok, leave the possibility
3980 // to not consider this.
3981 if(pdg == 22 && status != 1 && GetMCAnalysisUtils()->GetMCGenerator()!="" ) continue ;
3983 // If too small or too large pt, skip, same cut as for data analysis
3984 photonPt = lv.Pt () ;
3986 if(photonPt < GetMinPt() || photonPt > GetMaxPt() ) continue ;
3989 photonEta = lv.Eta() ;
3990 photonPhi = lv.Phi() ;
3992 if(photonPhi < 0) photonPhi+=TMath::TwoPi();
3994 // Check if photons hit the Calorimeter acceptance
3995 if(IsRealCaloAcceptanceOn() && fIsoDetector!="CTS") // defined on base class
3997 if(GetReader()->ReadStack() &&
3998 !GetCaloUtils()->IsMCParticleInCalorimeterAcceptance(fIsoDetector, primStack)) continue ;
3999 if(GetReader()->ReadAODMCParticles() &&
4000 !GetCaloUtils()->IsMCParticleInCalorimeterAcceptance(fIsoDetector, primAOD )) continue ;
4003 // Check same fidutial borders as in data analysis on top of real acceptance if real was requested.
4004 if(!GetFiducialCut()->IsInFiducialCut(lv,fIsoDetector)) continue ;
4006 // Get tag of this particle photon from fragmentation, decay, prompt ...
4007 // Set the origin of the photon.
4008 tag = GetMCAnalysisUtils()->CheckOrigin(i,GetReader(),fIsoDetector);
4010 if(pdg == 22 && !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
4012 // A conversion photon from a hadron, skip this kind of photon
4013 // printf("AliAnaPhoton::FillAcceptanceHistograms() - not a photon, weird!\n ");
4014 // GetMCAnalysisUtils()->PrintMCTag(tag);
4019 // Check the origin of the photon or if it is a pi0, assing a tag
4020 Int_t pi0d1Label = -1, pi0d2Label = -1;
4021 Bool_t overlapPi0 = kTRUE;
4024 mcIndex = kmcPrimPi0;
4025 //printf("check pi0\n");
4026 // Get the labels of the decay particles, remove them from isolation cone
4027 // Get also the opening angle and check if decays likely overlap
4028 Bool_t okpi0 = kFALSE;
4029 Int_t ndaugh = GetMCAnalysisUtils()->GetNDaughters(i,GetReader(), okpi0);
4030 //printf("OK pi0 %d, ndaugh %d\n",okpi0,ndaugh);
4031 Int_t d1Pdg = 0, d1Status = 0; Bool_t ok1 = kFALSE;
4032 Int_t d2Pdg = 0, d2Status = 0; Bool_t ok2 = kFALSE;
4033 TLorentzVector daugh1mom, daugh2mom;
4035 if ( ndaugh > 0 ) daugh1mom = GetMCAnalysisUtils()->GetDaughter(0,i,GetReader(),d1Pdg, d1Status,ok1, pi0d1Label);
4036 if ( ndaugh > 1 ) daugh2mom = GetMCAnalysisUtils()->GetDaughter(1,i,GetReader(),d2Pdg, d2Status,ok2, pi0d2Label);
4038 //printf("pi0 daug %d: a) %d, b) %d, c) %d\n", ndaugh,pi0d1Label,pi0d2Label);
4039 //if ( ndaugh !=2 ) printf("PDG: %d, %d, %d\n",d1Pdg,d2Pdg);
4041 // Select decays in 2 photons
4042 if( ndaugh!=2 || (d2Pdg != d1Pdg && d1Pdg!=22)) okpi0 = kFALSE;
4044 // Check overlap of decays
4045 if( okpi0 && fMakePrimaryPi0DecayStudy)
4047 Float_t d12Angle = daugh1mom.Angle(daugh2mom.Vect());
4048 if(d12Angle > overlapAngle) overlapPi0 = kFALSE;
4049 //printf(" -- d12 angle %2.3f, angle limit %2.3f, overlap %d\n",d12Angle,overlapAngle,overlapPi0);
4052 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt) )
4054 mcIndex = kmcPrimPrompt;
4056 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation) )
4058 mcIndex = kmcPrimFrag ;
4060 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR))
4062 mcIndex = kmcPrimISR;
4064 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))
4066 mcIndex = kmcPrimPi0Decay;
4068 else if( (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) ||
4069 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay)))
4071 mcIndex = kmcPrimOtherDecay;
4075 // Other decay but from non final state particle
4076 mcIndex = kmcPrimOtherDecay;
4079 // ////////////////////ISO MC/////////////////////////
4080 Double_t sumPtInCone = 0; Double_t dR=0. ;
4081 TParticle * mcisopStack = 0;
4082 AliAODMCParticle * mcisopAOD = 0;
4083 TLorentzVector mcisoLV;
4084 Int_t partInConeStatus = -1, partInConeMother = -1;
4085 Double_t partInConePt = 0, partInConeE = 0, partInConeEta = 0, partInConePhi = 0;
4086 Int_t partInConeCharge = 0, npart = 0;
4087 for(Int_t ip = 0; ip < nprim ; ip++)
4091 if( pdg==111 && ( ip == pi0d1Label || ip == pi0d2Label ) )
4093 //printf("Do not count pi0 decays in cone when isolating pi0 \n");
4097 if( GetReader()->ReadStack() )
4099 mcisopStack = static_cast<TParticle*>(stack->Particle(ip));
4100 if( !mcisopStack ) continue;
4101 partInConeStatus = mcisopStack->GetStatusCode();
4103 // Consider only final state particles, but this depends on generator,
4104 // status 1 is the usual one, in case of not being ok, leave the possibility
4105 // to not consider this.
4106 if( partInConeStatus != 1 && GetMCAnalysisUtils()->GetMCGenerator()!="" ) continue ;
4108 partInConeMother = mcisopStack->GetMother(0);
4109 partInConePt = mcisopStack->Pt();
4110 partInConeE = mcisopStack->Energy();
4111 partInConeEta = mcisopStack->Eta();
4112 partInConePhi = mcisopStack->Phi();
4113 partInConeCharge = TMath::Abs((Int_t) TDatabasePDG::Instance()->GetParticle(mcisopStack->GetPdgCode())->Charge());
4114 mcisopStack->Momentum(mcisoLV);
4118 mcisopAOD = (AliAODMCParticle *) mcparticles->At(ip);
4119 if( !mcisopAOD ) continue;
4121 partInConeStatus = mcisopAOD->GetStatus();
4122 // Consider only final state particles, but this depends on generator,
4123 // status 1 is the usual one, in case of not being ok, leave the possibility
4124 // to not consider this.
4125 if( partInConeStatus != 1 && GetMCAnalysisUtils()->GetMCGenerator()!="" ) continue ;
4127 partInConeMother = mcisopAOD->GetMother();
4128 partInConePt = mcisopAOD->Pt();
4129 partInConeE = mcisopAOD->E();
4130 partInConeEta = mcisopAOD->Eta();
4131 partInConePhi = mcisopAOD->Phi();
4132 partInConeCharge = TMath::Abs(mcisopAOD->Charge());
4133 mcisoLV.SetPxPyPzE(mcisopAOD->Px(),mcisopAOD->Py(),mcisopAOD->Pz(),mcisopAOD->E());
4136 if( partInConeMother == i ) continue;
4139 // Apply acceptance and energy/pt cut for particles in cone
4140 if(fSelectPrimariesInCone)
4142 if( partInConeCharge > 0) // charged pT cut and acceptance
4144 if( GetIsolationCut()->GetParticleTypeInCone() == AliIsolationCut::kOnlyNeutral ) continue;
4146 if( partInConePt < GetReader()->GetCTSPtMin () ) continue;
4148 if(!GetReader()->GetFiducialCut()->IsInFiducialCut(mcisoLV,"CTS")) continue ;
4150 else // neutrals E cut and acceptance
4152 if( GetIsolationCut()->GetParticleTypeInCone() == AliIsolationCut::kOnlyCharged ) continue;
4154 if( partInConeE <= minECalo ) continue;
4156 if(!GetReader()->GetFiducialCut()->IsInFiducialCut(mcisoLV,fCalorimeter)) continue ;
4158 if(IsRealCaloAcceptanceOn()) // defined on base class
4160 if(GetReader()->ReadStack() &&
4161 !GetCaloUtils()->IsMCParticleInCalorimeterAcceptance(fCalorimeter, mcisopStack)) continue ;
4162 if(GetReader()->ReadAODMCParticles() &&
4163 !GetCaloUtils()->IsMCParticleInCalorimeterAcceptance(fCalorimeter, mcisopAOD )) continue ;
4169 dR = GetIsolationCut()->Radius(photonEta, photonPhi, partInConeEta, partInConePhi);
4171 if(dR > GetIsolationCut()->GetConeSize())
4174 sumPtInCone += partInConePt;
4175 if(partInConePt > GetIsolationCut()->GetPtThreshold() &&
4176 partInConePt < GetIsolationCut()->GetPtThresholdMax()) npart++;
4179 ///////END ISO MC/////////////////////////
4181 // Fill the histograms, only those in the defined calorimeter acceptance
4183 fhEtaPrimMC[kmcPrimPhoton]->Fill(photonPt , photonEta) ;
4184 fhPhiPrimMC[kmcPrimPhoton]->Fill(photonPt , photonPhi) ;
4185 fhEPrimMC [kmcPrimPhoton]->Fill(photonE) ;
4186 fhPtPrimMC [kmcPrimPhoton]->Fill(photonPt) ;
4188 fhEtaPrimMC[mcIndex]->Fill(photonPt , photonEta) ;
4189 fhPhiPrimMC[mcIndex]->Fill(photonPt , photonPhi) ;
4190 fhEPrimMC [mcIndex]->Fill(photonE ) ;
4191 fhPtPrimMC [mcIndex]->Fill(photonPt) ;
4193 // In case the photon is a decay from pi0,
4194 // study how the decay kinematics affects the isolation
4195 TLorentzVector pi0mom, daugh1mom, daugh2mom;
4197 Bool_t okpi0 = 0, ok1 = 0, ok2 = 0;
4198 Int_t pi0label = -1, d1Label = -1, d2Label = -1;
4199 Bool_t d2Acc = kTRUE, overlap = kTRUE;
4201 Float_t dRdaugh2 = 0, d12Angle = 0;
4202 if(mcIndex == kmcPrimPi0Decay && fMakePrimaryPi0DecayStudy)
4204 pi0mom = GetMCAnalysisUtils()->GetMotherWithPDG(i,111,GetReader(),okpi0, pi0label);
4207 ndaugh = GetMCAnalysisUtils()->GetNDaughters(pi0label,GetReader(), okpi0);
4210 Int_t d1Pdg = 0, d1Status = 0;
4211 daugh1mom = GetMCAnalysisUtils()->GetDaughter(0,pi0label,GetReader(),d1Pdg, d1Status,ok1, d1Label);
4212 Int_t d2Pdg = 0, d2Status = 0;
4213 daugh2mom = GetMCAnalysisUtils()->GetDaughter(1,pi0label,GetReader(),d2Pdg, d2Status,ok2, d2Label);
4214 if(d2Pdg != d1Pdg && d1Pdg!=22) okpi0 = kFALSE;
4216 // Check the momentum and location of second daughter
4219 // assign current trigger to first daughter
4222 Int_t tmpLabel = d2Label;
4225 TLorentzVector tmpLV = daugh2mom;
4226 daugh2mom = daugh1mom;
4230 // Check if photons hit the Calorimeter acceptance
4231 if(IsRealCaloAcceptanceOn() && fIsoDetector!="CTS") // defined on base class
4232 d2Acc = GetCaloUtils()->IsMCParticleInCalorimeterAcceptance(fIsoDetector,daugh2mom,d2AbsId) ;
4234 //printf("D2 (eta %2.2f,phi %2.2f)in real calo %d, with absId %d\n",
4235 // daugh2mom.Eta(), daugh2mom.Phi()*TMath::RadToDeg(),d2Acc,d2AbsId);
4237 // Check same fidutial borders as in data analysis on top of real acceptance if real was requested.
4238 if(d2Acc) d2Acc = GetReader()->GetFiducialCut()->IsInFiducialCut(daugh2mom,fIsoDetector);
4239 //printf("D2 fidcut %d\n",d2Acc);
4241 Float_t phiDaugh2 = daugh2mom.Phi();
4242 if(phiDaugh2 < 0) phiDaugh2+=TMath::TwoPi();
4243 dRdaugh2 = GetIsolationCut()->Radius(photonEta, photonPhi, daugh2mom.Eta(),phiDaugh2);
4245 // Opening angle, check if pairs will likely overlap
4246 d12Angle = daugh1mom.Angle(daugh2mom.Vect());
4247 if(d12Angle > overlapAngle) overlap = kFALSE;
4253 //printf("Check mother of label %d: mom label %d, okmom %d ndaugh %d, daugh label1 %d, label2 %d, ok1 %d, ok2 %d, R %2.3f, opening angle %2.3f, overlap %d\n",
4254 // i, pi0label,okpi0,ndaugh,d1Label,d2Label,ok1,ok2, dRdaugh2, d12Angle, overlap);
4256 // Second decay out of cone
4257 if(dRdaugh2 > GetIsolationCut()->GetConeSize())
4258 fhPtPrimMCPi0DecayPairOutOfCone->Fill(photonPt);
4260 // Second decay out of acceptance
4261 if(!ok2 || !d2Acc || daugh2mom.E() <= minECalo)
4262 fhPtPrimMCPi0DecayPairOutOfAcceptance->Fill(photonPt);
4264 // Not Overlapped decay
4265 if(!overlap) fhPtPrimMCPi0DecayPairNoOverlap->Fill(photonPt);
4267 // Second decay pt smaller than threshold
4268 if(d2Acc && dRdaugh2 < GetIsolationCut()->GetConeSize() &&
4269 daugh2mom.E() < GetIsolationCut()->GetPtThreshold())
4271 fhPtPrimMCPi0DecayPairAcceptInConeLowPt->Fill(photonPt);
4274 fhPtPrimMCPi0DecayPairAcceptInConeLowPtNoOverlap->Fill(photonPt);
4275 if(daugh2mom.E() > minECalo)fhPtPrimMCPi0DecayPairAcceptInConeLowPtNoOverlapCaloE->Fill(photonPt);
4280 if(mcIndex == kmcPrimPi0 && fMakePrimaryPi0DecayStudy && overlapPi0)
4281 fhPtPrimMCPi0Overlap->Fill(photonPt);
4284 Bool_t isolated = kFALSE;
4285 if(GetIsolationCut()->GetICMethod() == AliIsolationCut::kSumPtIC &&
4286 (sumPtInCone < GetIsolationCut()->GetSumPtThreshold() ||
4287 sumPtInCone > GetIsolationCut()->GetSumPtThresholdMax()))
4290 if(GetIsolationCut()->GetICMethod() == AliIsolationCut::kPtThresIC &&
4296 fhPtPrimMCiso [mcIndex] ->Fill(photonPt) ;
4297 fhPtPrimMCiso [kmcPrimPhoton]->Fill(photonPt) ;
4298 if(mcIndex == kmcPrimPi0Decay && fMakePrimaryPi0DecayStudy)
4300 // Not Overlapped decay
4301 if(!overlap) fhPtPrimMCPi0DecayIsoPairNoOverlap->Fill(photonPt);
4303 // Second decay out of cone
4304 if(dRdaugh2 > GetIsolationCut()->GetConeSize())
4305 fhPtPrimMCPi0DecayIsoPairOutOfCone->Fill(photonPt);
4307 // Second decay out of acceptance
4308 if(!ok2 || !d2Acc || daugh2mom.E() <= minECalo)
4309 fhPtPrimMCPi0DecayIsoPairOutOfAcceptance->Fill(photonPt);
4311 // Second decay pt smaller than threshold
4312 if(d2Acc && dRdaugh2 < GetIsolationCut()->GetConeSize() &&
4313 daugh2mom.E() < GetIsolationCut()->GetPtThreshold())
4315 fhPtPrimMCPi0DecayIsoPairAcceptInConeLowPt->Fill(photonPt);
4318 fhPtPrimMCPi0DecayIsoPairAcceptInConeLowPtNoOverlap->Fill(photonPt);
4319 if(daugh2mom.E() > minECalo) fhPtPrimMCPi0DecayIsoPairAcceptInConeLowPtNoOverlapCaloE->Fill(photonPt);
4324 if(mcIndex == kmcPrimPi0 && fMakePrimaryPi0DecayStudy && overlapPi0)
4325 fhPtPrimMCPi0IsoOverlap->Fill(photonPt);
4328 }//loop on primaries
4330 if(GetDebug() > 0) printf("AliAnaParticleIsolation::FillAcceptanceHistograms() - End \n");
4335 //_____________________________________________________________________________________
4336 void AliAnaParticleIsolation::MakeSeveralICAnalysis(AliAODPWG4ParticleCorrelation* ph,
4340 //Isolation Cut Analysis for both methods and different pt cuts and cones
4341 Float_t ptC = ph->Pt();
4342 Float_t etaC = ph->Eta();
4343 Float_t phiC = ph->Phi();
4344 Int_t tag = ph->GetTag();
4347 if(fFillTaggedDecayHistograms)
4349 decayTag = ph->GetBtag(); // temporary
4350 if(decayTag < 0) decayTag = 0; // temporary
4354 printf("AliAnaParticleIsolation::MakeSeveralICAnalysis() - Isolate pT %2.2f, decay tag %d\n",ptC, decayTag);
4356 //Keep original setting used when filling AODs, reset at end of analysis
4357 Float_t ptthresorg = GetIsolationCut()->GetPtThreshold();
4358 Float_t ptfracorg = GetIsolationCut()->GetPtFraction();
4359 Float_t ptsumcorg = GetIsolationCut()->GetSumPtThreshold();
4360 Float_t rorg = GetIsolationCut()->GetConeSize();
4362 Float_t coneptsum = 0, coneptlead = 0;
4363 Int_t n [10][10];//[fNCones][fNPtThresFrac];
4364 Int_t nfrac[10][10];//[fNCones][fNPtThresFrac];
4365 Bool_t isolated = kFALSE;
4367 // Fill hist with all particles before isolation criteria
4368 fhENoIso ->Fill(ph->E());
4369 fhPtNoIso ->Fill(ptC);
4370 fhEtaPhiNoIso->Fill(etaC,phiC);
4374 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
4375 fhPtNoIsoMC[kmcPhoton]->Fill(ptC);
4377 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCDecayPairLost) && mcIndex==kmcPi0Decay )
4378 fhPtNoIsoMC[kmcPi0DecayLostPair]->Fill(ptC);
4380 fhPtNoIsoMC[mcIndex]->Fill(ptC);
4383 // Candidates tagged as decay in another analysis (AliAnaPi0EbE)
4384 if(fFillTaggedDecayHistograms && decayTag > 0)
4386 for(Int_t ibit = 0; ibit < fNDecayBits; ibit++)
4388 if(GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[ibit]))
4390 fhPtDecayNoIso[ibit] ->Fill(ptC);
4391 fhEtaPhiDecayNoIso[ibit]->Fill(etaC,phiC);
4395 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
4396 fhPtDecayNoIsoMC[ibit][kmcPhoton]->Fill(ptC);
4398 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCDecayPairLost) && mcIndex==kmcPi0Decay )
4399 fhPtDecayNoIsoMC[ibit][kmcPi0DecayLostPair]->Fill(ptC);
4401 fhPtDecayNoIsoMC[ibit][mcIndex]->Fill(ptC);
4405 } // decay histograms
4407 //Get vertex for photon momentum calculation
4408 Double_t vertex[] = {0,0,0} ; //vertex ;
4409 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
4410 GetReader()->GetVertex(vertex);
4412 //Loop on cone sizes
4413 for(Int_t icone = 0; icone<fNCones; icone++)
4415 //Recover reference arrays with clusters and tracks
4416 TObjArray * refclusters = ph->GetObjArray(GetAODObjArrayName()+"Clusters");
4417 TObjArray * reftracks = ph->GetObjArray(GetAODObjArrayName()+"Tracks");
4419 //If too small or too large pt, skip
4420 if(ptC < GetMinPt() || ptC > GetMaxPt() ) continue ;
4422 //In case a more strict IC is needed in the produced AOD
4424 isolated = kFALSE; coneptsum = 0; coneptlead = 0;
4426 GetIsolationCut()->SetSumPtThreshold(100);
4427 GetIsolationCut()->SetPtThreshold(100);
4428 GetIsolationCut()->SetPtFraction(100);
4429 GetIsolationCut()->SetConeSize(fConeSizes[icone]);
4431 // retreive pt tracks to fill histo vs. pt leading
4432 //Fill pt distribution of particles in cone
4433 //fhPtLeadingPt(),fhPerpSumPtLeadingPt(),fhPerpPtLeadingPt(),
4435 // Tracks in perpendicular cones
4436 Double_t sumptPerp = 0. ;
4437 TObjArray * trackList = GetCTSTracks() ;
4438 for(Int_t itrack=0; itrack < trackList->GetEntriesFast(); itrack++)
4440 AliVTrack* track = (AliVTrack *) trackList->At(itrack);
4441 //fill the histograms at forward range
4444 printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Track not available?");
4448 Double_t dPhi = phiC - track->Phi() + TMath::PiOver2();
4449 Double_t dEta = etaC - track->Eta();
4450 Double_t arg = dPhi*dPhi + dEta*dEta;
4451 if(TMath::Sqrt(arg) < fConeSizes[icone])
4453 fhPerpPtLeadingPt[icone]->Fill(ptC,TMath::Sqrt(track->Px()*track->Px()+track->Py()*track->Py()));
4454 sumptPerp+=track->Pt();
4457 dPhi = phiC - track->Phi() - TMath::PiOver2();
4458 arg = dPhi*dPhi + dEta*dEta;
4459 if(TMath::Sqrt(arg) < fConeSizes[icone])
4461 fhPerpPtLeadingPt[icone]->Fill(ptC,TMath::Sqrt(track->Px()*track->Px()+track->Py()*track->Py()));
4462 sumptPerp+=track->Pt();
4466 fhPerpSumPtLeadingPt[icone]->Fill(ptC,sumptPerp);
4468 // Tracks in isolation cone, pT distribution and sum
4469 if(reftracks && GetIsolationCut()->GetParticleTypeInCone()!= AliIsolationCut::kOnlyNeutral)
4471 for(Int_t itrack=0; itrack < reftracks->GetEntriesFast(); itrack++)
4473 AliVTrack* track = (AliVTrack *) reftracks->At(itrack);
4475 Float_t rad = GetIsolationCut()->Radius(etaC, phiC, track->Eta(), track->Phi());
4477 if(rad > fConeSizes[icone]) continue ;
4479 fhPtLeadingPt[icone]->Fill(ptC, track->Pt());
4480 coneptsum += track->Pt();
4484 // Clusters in isolation cone, pT distribution and sum
4485 if(refclusters && GetIsolationCut()->GetParticleTypeInCone()!= AliIsolationCut::kOnlyCharged)
4487 TLorentzVector mom ;
4488 for(Int_t icalo=0; icalo < refclusters->GetEntriesFast(); icalo++)
4490 AliVCluster* calo = (AliVCluster *) refclusters->At(icalo);
4492 calo->GetMomentum(mom,vertex) ;//Assume that come from vertex in straight line
4494 Float_t rad = GetIsolationCut()->Radius(etaC, phiC, mom.Eta(), mom.Phi());
4496 if(rad > fConeSizes[icone]) continue ;
4498 fhPtLeadingPt[icone]->Fill(ptC, mom.Pt());
4499 coneptsum += mom.Pt();
4503 fhSumPtLeadingPt[icone]->Fill(ptC,coneptsum);
4507 if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
4508 fhSumPtLeadingPtMC[kmcPhoton][icone]->Fill(ptC,coneptsum) ;
4510 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCDecayPairLost) && mcIndex==kmcPi0Decay )
4511 fhSumPtLeadingPtMC[kmcPi0DecayLostPair][icone]->Fill(ptC,coneptsum) ;
4513 fhSumPtLeadingPtMC[mcIndex][icone]->Fill(ptC,coneptsum) ;
4518 //Loop on pt thresholds
4519 for(Int_t ipt = 0; ipt < fNPtThresFrac ; ipt++)
4522 nfrac[icone][ipt]=0;
4523 GetIsolationCut()->SetPtThreshold(fPtThresholds[ipt]);
4524 GetIsolationCut()->SetPtFraction(fPtFractions[ipt]) ;
4525 GetIsolationCut()->SetSumPtThreshold(fSumPtThresholds[ipt]);
4527 GetIsolationCut()->MakeIsolationCut(reftracks, refclusters,
4528 GetReader(), GetCaloPID(),
4530 n[icone][ipt],nfrac[icone][ipt],
4531 coneptsum, coneptlead, isolated);
4533 // Normal pT threshold cut
4537 printf(" AliAnaParticleIsolation::MakeSeveralICAnalysis() - cone size %1.1f, ptThres %1.1f, sumptThresh %1.1f\n",
4538 fConeSizes[icone],fPtThresholds[ipt],fSumPtThresholds[ipt]);
4539 printf("\t n %d, nfrac %d, coneptsum %2.2f\n",
4540 n[icone][ipt],nfrac[icone][ipt],coneptsum);
4542 printf(" AliAnaParticleIsolation::MakeSeveralICAnalysis() - pt %1.1f, eta %1.1f, phi %1.1f\n",ptC, etaC, phiC);
4545 if(n[icone][ipt] == 0)
4548 printf(" AliAnaParticleIsolation::MakeSeveralICAnalysis() - filling pt threshold loop\n");
4550 fhPtThresIsolated [icone][ipt]->Fill(ptC);
4551 fhEtaPhiPtThresIso[icone][ipt]->Fill(etaC,phiC);
4553 if( fFillTaggedDecayHistograms && decayTag > 0 && fNDecayBits > 0)
4555 if(GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[0]))
4557 fhPtPtThresDecayIso [icone][ipt]->Fill(ptC);
4558 fhEtaPhiPtThresDecayIso[icone][ipt]->Fill(etaC,phiC);
4564 if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) )
4565 fhPtThresIsolatedMC[kmcPhoton][icone][ipt]->Fill(ptC) ;
4567 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCDecayPairLost) && mcIndex==kmcPi0Decay )
4568 fhPtThresIsolatedMC[kmcPi0DecayLostPair][icone][ipt]->Fill(ptC) ;
4570 fhPtThresIsolatedMC[mcIndex][icone][ipt]->Fill(ptC) ;
4575 // pt in cone fraction
4576 if(nfrac[icone][ipt] == 0)
4579 printf(" AliAnaParticleIsolation::MakeSeveralICAnalysis() - filling frac loop\n");
4581 fhPtFracIsolated [icone][ipt]->Fill(ptC);
4582 fhEtaPhiPtFracIso[icone][ipt]->Fill(etaC,phiC);
4584 if( fFillTaggedDecayHistograms && decayTag > 0 && fNDecayBits > 0)
4586 if(GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[0]))
4588 fhPtPtFracDecayIso [icone][ipt]->Fill(ptC);
4589 fhEtaPhiPtFracDecayIso[icone][ipt]->Fill(etaC,phiC);
4595 if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
4596 fhPtFracIsolatedMC[kmcPhoton][icone][ipt]->Fill(ptC) ;
4598 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCDecayPairLost) && mcIndex==kmcPi0Decay )
4599 fhPtFracIsolatedMC[kmcPi0DecayLostPair][icone][ipt]->Fill(ptC) ;
4601 fhPtFracIsolatedMC[mcIndex][icone][ipt]->Fill(ptC) ;
4606 printf(" AliAnaParticleIsolation::MakeSeveralICAnalysis() - checking IC method : %i\n",GetIsolationCut()->GetICMethod());
4608 //Pt threshold on pt cand/ sum in cone histograms
4609 if(coneptsum < fSumPtThresholds[ipt])
4612 printf(" AliAnaParticleIsolation::MakeSeveralICAnalysis() - filling sum loop\n");
4614 fhSumPtIsolated [icone][ipt]->Fill(ptC) ;
4615 fhEtaPhiPtSumIso[icone][ipt]->Fill(etaC, phiC) ;
4617 if( fFillTaggedDecayHistograms && decayTag > 0 && fNDecayBits > 0)
4619 if(GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[0]))
4621 fhPtPtSumDecayIso[icone][ipt]->Fill(ptC);
4622 fhEtaPhiPtSumDecayIso[icone][ipt]->Fill(etaC, phiC) ;
4628 if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
4629 fhSumPtIsolatedMC[kmcPhoton][icone][ipt]->Fill(ptC) ;
4631 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCDecayPairLost) && mcIndex==kmcPi0Decay )
4632 fhSumPtIsolatedMC[kmcPi0DecayLostPair][icone][ipt]->Fill(ptC) ;
4634 fhSumPtIsolatedMC[mcIndex][icone][ipt]->Fill(ptC) ;
4638 // pt sum pt frac method
4639 // if( ((fPtFractions[ipt]*ptC < fSumPtThresholds[ipt]) && (coneptsum < fSumPtThresholds[ipt])) || ((fPtFractions[ipt]*ptC > fSumPtThresholds[ipt]) && (coneptsum < fPtFractions[ipt]*ptC)) )
4641 if(coneptsum < fPtFractions[ipt]*ptC)
4644 printf(" AliAnaParticleIsolation::MakeSeveralICAnalysis() - filling PtFrac PtSum loop\n");
4646 fhPtFracPtSumIso [icone][ipt]->Fill(ptC) ;
4647 fhEtaPhiFracPtSumIso[icone][ipt]->Fill(etaC,phiC) ;
4649 if( fFillTaggedDecayHistograms && decayTag > 0 && fNDecayBits > 0)
4651 if(GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[0]))
4653 fhPtFracPtSumDecayIso [icone][ipt]->Fill(ptC);
4654 fhEtaPhiFracPtSumDecayIso[icone][ipt]->Fill(etaC,phiC);
4660 Float_t cellDensity = GetIsolationCut()->GetCellDensity( ph, GetReader());
4661 if(coneptsum < fSumPtThresholds[ipt]*cellDensity)
4664 printf(" AliAnaParticleIsolation::MakeSeveralICAnalysis() - filling density loop\n");
4666 fhPtSumDensityIso [icone][ipt]->Fill(ptC) ;
4667 fhEtaPhiSumDensityIso[icone][ipt]->Fill(etaC,phiC) ;
4669 if( fFillTaggedDecayHistograms && decayTag > 0 && fNDecayBits > 0)
4671 if(GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[0]))
4673 fhPtSumDensityDecayIso [icone][ipt]->Fill(ptC);
4674 fhEtaPhiSumDensityDecayIso[icone][ipt]->Fill(etaC, phiC);
4683 //Reset original parameters for AOD analysis
4684 GetIsolationCut()->SetPtThreshold(ptthresorg);
4685 GetIsolationCut()->SetPtFraction(ptfracorg);
4686 GetIsolationCut()->SetSumPtThreshold(ptsumcorg);
4687 GetIsolationCut()->SetConeSize(rorg);
4691 //_____________________________________________________________
4692 void AliAnaParticleIsolation::Print(const Option_t * opt) const
4695 //Print some relevant parameters set for the analysis
4699 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
4700 AliAnaCaloTrackCorrBaseClass::Print(" ");
4702 printf("ReMake Isolation = %d \n", fReMakeIC) ;
4703 printf("Make Several Isolation = %d \n", fMakeSeveralIC) ;
4704 printf("Calorimeter for isolation = %s \n", fCalorimeter.Data()) ;
4705 printf("Detector for candidate isolation = %s \n", fIsoDetector.Data()) ;
4709 printf("N Cone Sizes = %d\n", fNCones) ;
4710 printf("Cone Sizes = \n") ;
4711 for(Int_t i = 0; i < fNCones; i++)
4712 printf(" %1.2f;", fConeSizes[i]) ;
4715 printf("N pT thresholds/fractions = %d\n", fNPtThresFrac) ;
4716 printf(" pT thresholds = \n") ;
4717 for(Int_t i = 0; i < fNPtThresFrac; i++)
4718 printf(" %2.2f;", fPtThresholds[i]) ;
4722 printf(" pT fractions = \n") ;
4723 for(Int_t i = 0; i < fNPtThresFrac; i++)
4724 printf(" %2.2f;", fPtFractions[i]) ;
4728 printf("sum pT thresholds = \n") ;
4729 for(Int_t i = 0; i < fNPtThresFrac; i++)
4730 printf(" %2.2f;", fSumPtThresholds[i]) ;