]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGGA/CaloTrackCorrelations/AliAnaPhoton.cxx
Added Print method
[u/mrichter/AliRoot.git] / PWGGA / CaloTrackCorrelations / AliAnaPhoton.cxx
CommitLineData
a3aebfff 1 /**************************************************************************
1c5acb87 2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes hereby granted *
cadbb0f3 9 * without fee, provided that the above copyright notice appears in all *
1c5acb87 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 **************************************************************************/
1c5acb87 15
16//_________________________________________________________________________
17//
18// Class for the photon identification.
19// Clusters from calorimeters are identified as photons
20// and kept in the AOD. Few histograms produced.
6175da48 21// Produces input for other analysis classes like AliAnaPi0,
22// AliAnaParticleHadronCorrelation ...
1c5acb87 23//
24// -- Author: Gustavo Conesa (LNF-INFN)
25//////////////////////////////////////////////////////////////////////////////
26
27
28// --- ROOT system ---
29#include <TH2F.h>
2244659d 30#include <TH3D.h>
477d6cee 31#include <TClonesArray.h>
0c1383b5 32#include <TObjString.h>
123fc3bd 33#include "TParticle.h"
6175da48 34#include "TDatabasePDG.h"
1c5acb87 35
36// --- Analysis system ---
37#include "AliAnaPhoton.h"
38#include "AliCaloTrackReader.h"
123fc3bd 39#include "AliStack.h"
1c5acb87 40#include "AliCaloPID.h"
6639984f 41#include "AliMCAnalysisUtils.h"
ff45398a 42#include "AliFiducialCut.h"
0ae57829 43#include "AliVCluster.h"
591cc579 44#include "AliAODMCParticle.h"
c8fe2783 45#include "AliMixedEvent.h"
fc195fd0 46#include "AliAODEvent.h"
c8fe2783 47
c5693f62 48// --- Detectors ---
49#include "AliPHOSGeoUtils.h"
50#include "AliEMCALGeometry.h"
1c5acb87 51
52ClassImp(AliAnaPhoton)
53
5812a064 54//____________________________
521636d2 55AliAnaPhoton::AliAnaPhoton() :
09273901 56 AliAnaCaloTrackCorrBaseClass(), fCalorimeter(""),
521636d2 57 fMinDist(0.), fMinDist2(0.), fMinDist3(0.),
09273901 58 fRejectTrackMatch(0), fFillTMHisto(kFALSE),
59 fTimeCutMin(-10000), fTimeCutMax(10000),
60 fNCellsCut(0), fFillSSHistograms(kFALSE),
f66d95af 61 fNOriginHistograms(8), fNPrimaryHistograms(4),
521636d2 62
c4a7d28a 63 // Histograms
5c46c992 64 fhNCellsE(0), fhCellsE(0), // Control histograms
65 fhMaxCellDiffClusterE(0), fhTimeE(0), // Control histograms
f66d95af 66 fhEPhoton(0), fhPtPhoton(0),
521636d2 67 fhPhiPhoton(0), fhEtaPhoton(0),
68 fhEtaPhiPhoton(0), fhEtaPhi05Photon(0),
69
4c8f7c2e 70 // Shower shape histograms
521636d2 71 fhDispE(0), fhLam0E(0), fhLam1E(0),
521636d2 72 fhDispETRD(0), fhLam0ETRD(0), fhLam1ETRD(0),
b5dbb99b 73 fhDispETM(0), fhLam0ETM(0), fhLam1ETM(0),
74 fhDispETMTRD(0), fhLam0ETMTRD(0), fhLam1ETMTRD(0),
521636d2 75
76 fhNCellsLam0LowE(0), fhNCellsLam1LowE(0), fhNCellsDispLowE(0),
77 fhNCellsLam0HighE(0), fhNCellsLam1HighE(0), fhNCellsDispHighE(0),
521636d2 78
79 fhEtaLam0LowE(0), fhPhiLam0LowE(0),
80 fhEtaLam0HighE(0), fhPhiLam0HighE(0),
81 fhLam0DispLowE(0), fhLam0DispHighE(0),
82 fhLam1Lam0LowE(0), fhLam1Lam0HighE(0),
83 fhDispLam1LowE(0), fhDispLam1HighE(0),
34c16486 84 fhDispEtaE(0), fhDispPhiE(0),
85 fhSumEtaE(0), fhSumPhiE(0), fhSumEtaPhiE(0),
86 fhDispEtaPhiDiffE(0), fhSphericityE(0),
87 fhDispSumEtaDiffE(0), fhDispSumPhiDiffE(0),
521636d2 88
4c8f7c2e 89 // MC histograms
8d6b7f60 90 fhMCPhotonELambda0NoOverlap(0), fhMCPhotonELambda0TwoOverlap(0), fhMCPhotonELambda0NOverlap(0),
b5dbb99b 91 // Embedding
7c65ad18 92 fhEmbeddedSignalFractionEnergy(0),
8d6b7f60 93 fhEmbedPhotonELambda0FullSignal(0), fhEmbedPhotonELambda0MostlySignal(0),
94 fhEmbedPhotonELambda0MostlyBkg(0), fhEmbedPhotonELambda0FullBkg(0),
95 fhEmbedPi0ELambda0FullSignal(0), fhEmbedPi0ELambda0MostlySignal(0),
4bfeae64 96 fhEmbedPi0ELambda0MostlyBkg(0), fhEmbedPi0ELambda0FullBkg(0)
b5dbb99b 97
4bfeae64 98 {
1c5acb87 99 //default ctor
100
4bfeae64 101 for(Int_t i = 0; i < 14; i++)
102 {
4c8f7c2e 103 fhMCPt [i] = 0;
104 fhMCE [i] = 0;
105 fhMCPhi [i] = 0;
106 fhMCEta [i] = 0;
107 fhMCDeltaE [i] = 0;
108 fhMCDeltaPt[i] = 0;
4c8f7c2e 109 fhMC2E [i] = 0;
110 fhMC2Pt [i] = 0;
521636d2 111 }
112
4bfeae64 113 for(Int_t i = 0; i < 7; i++)
114 {
3d5d5078 115 fhPtPrimMC [i] = 0;
116 fhEPrimMC [i] = 0;
117 fhPhiPrimMC[i] = 0;
118 fhYPrimMC [i] = 0;
119
120 fhPtPrimMCAcc [i] = 0;
121 fhEPrimMCAcc [i] = 0;
122 fhPhiPrimMCAcc[i] = 0;
123 fhYPrimMCAcc [i] = 0;
124 }
125
4bfeae64 126 for(Int_t i = 0; i < 6; i++)
127 {
f66d95af 128 fhMCELambda0 [i] = 0;
129 fhMCELambda1 [i] = 0;
130 fhMCEDispersion [i] = 0;
f66d95af 131 fhMCNCellsE [i] = 0;
132 fhMCMaxCellDiffClusterE[i] = 0;
bfdcf7fb 133 fhLambda0DispEta[i] = 0;
134 fhLambda0DispPhi[i] = 0;
135
f66d95af 136 fhMCLambda0vsClusterMaxCellDiffE0[i] = 0;
137 fhMCLambda0vsClusterMaxCellDiffE2[i] = 0;
138 fhMCLambda0vsClusterMaxCellDiffE6[i] = 0;
139 fhMCNCellsvsClusterMaxCellDiffE0 [i] = 0;
140 fhMCNCellsvsClusterMaxCellDiffE2 [i] = 0;
141 fhMCNCellsvsClusterMaxCellDiffE6 [i] = 0;
34c16486 142
143 fhMCEDispEta [i] = 0;
144 fhMCEDispPhi [i] = 0;
145 fhMCESumEtaPhi [i] = 0;
146 fhMCEDispEtaPhiDiff[i] = 0;
147 fhMCESphericity [i] = 0;
521636d2 148 }
149
34c16486 150 for(Int_t i = 0; i < 5; i++)
151 {
152 fhClusterCuts[i] = 0;
bfdcf7fb 153 fhDispEtaDispPhi[i] = 0;
34c16486 154 for(Int_t j = 0; j < 6; j++)
bfdcf7fb 155 {
156 fhMCDispEtaDispPhi[i][j] = 0;
157 fhMCLambda0DispEta [i][j] = 0;
158 fhMCLambda0DispPhi [i][j] = 0;
159 }
34c16486 160 }
fc195fd0 161
4bfeae64 162 // Track matching residuals
163 for(Int_t i = 0; i < 2; i++)
164 {
165 fhTrackMatchedDEta[i] = 0; fhTrackMatchedDPhi[i] = 0; fhTrackMatchedDEtaDPhi[i] = 0;
166 fhTrackMatchedDEtaTRD[i] = 0; fhTrackMatchedDPhiTRD[i] = 0;
167 fhTrackMatchedDEtaMCOverlap[i] = 0; fhTrackMatchedDPhiMCOverlap[i] = 0;
168 fhTrackMatchedDEtaMCNoOverlap[i] = 0; fhTrackMatchedDPhiMCNoOverlap[i] = 0;
169 fhTrackMatchedDEtaMCConversion[i] = 0; fhTrackMatchedDPhiMCConversion[i] = 0;
170 fhTrackMatchedMCParticle[i] = 0; fhTrackMatchedMCParticle[i] = 0;
171 fhdEdx[i] = 0; fhEOverP[i] = 0;
172 fhEOverPTRD[i] = 0;
173 }
174
1c5acb87 175 //Initialize parameters
176 InitParameters();
177
1c5acb87 178}
179
5812a064 180//__________________________________________________________________________
c4a7d28a 181Bool_t AliAnaPhoton::ClusterSelected(AliVCluster* calo, TLorentzVector mom)
182{
183 //Select clusters if they pass different cuts
184 if(GetDebug() > 2)
185 printf("AliAnaPhoton::ClusterSelected() Current Event %d; Before selection : E %2.2f, pT %2.2f, Ecl %2.2f, phi %2.2f, eta %2.2f\n",
186 GetReader()->GetEventNumber(),
fc195fd0 187 calo->E(), mom.Pt(),calo->E(),mom.Phi()*TMath::RadToDeg(),mom.Eta());
188
189 fhClusterCuts[1]->Fill(calo->E());
c4a7d28a 190
191 //.......................................
192 //If too small or big energy, skip it
fc195fd0 193 if(calo->E() < GetMinEnergy() || calo->E() > GetMaxEnergy() ) return kFALSE ;
09273901 194
c4a7d28a 195 if(GetDebug() > 2) printf("\t Cluster %d Pass E Cut \n",calo->GetID());
09273901 196
fc195fd0 197 fhClusterCuts[2]->Fill(calo->E());
198
c4a7d28a 199 //.......................................
200 // TOF cut, BE CAREFUL WITH THIS CUT
201 Double_t tof = calo->GetTOF()*1e9;
202 if(tof < fTimeCutMin || tof > fTimeCutMax) return kFALSE;
09273901 203
c4a7d28a 204 if(GetDebug() > 2) printf("\t Cluster %d Pass Time Cut \n",calo->GetID());
09273901 205
fc195fd0 206 fhClusterCuts[3]->Fill(calo->E());
207
c4a7d28a 208 //.......................................
209 if(calo->GetNCells() <= fNCellsCut && GetReader()->GetDataType() != AliCaloTrackReader::kMC) return kFALSE;
09273901 210
c4a7d28a 211 if(GetDebug() > 2) printf("\t Cluster %d Pass NCell Cut \n",calo->GetID());
09273901 212
fc195fd0 213 fhClusterCuts[4]->Fill(calo->E());
214
c4a7d28a 215 //.......................................
216 //Check acceptance selection
34c16486 217 if(IsFiducialCutOn())
218 {
c4a7d28a 219 Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
220 if(! in ) return kFALSE ;
221 }
09273901 222
c4a7d28a 223 if(GetDebug() > 2) printf("Fiducial cut passed \n");
09273901 224
fc195fd0 225 fhClusterCuts[5]->Fill(calo->E());
226
c4a7d28a 227 //.......................................
228 //Skip matched clusters with tracks
09273901 229
4bfeae64 230 // Fill matching residual histograms before PID cuts
231 if(fFillTMHisto) FillTrackMatchingResidualHistograms(calo,0);
09273901 232
34c16486 233 if(fRejectTrackMatch)
234 {
235 if(IsTrackMatched(calo,GetReader()->GetInputEvent()))
236 {
c4a7d28a 237 if(GetDebug() > 2) printf("\t Reject track-matched clusters\n");
238 return kFALSE ;
239 }
240 else
241 if(GetDebug() > 2) printf(" Track-matching cut passed \n");
242 }// reject matched clusters
09273901 243
fc195fd0 244 fhClusterCuts[6]->Fill(calo->E());
245
c4a7d28a 246 //.......................................
247 //Check Distance to Bad channel, set bit.
248 Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
249 if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
34c16486 250 if(distBad < fMinDist)
251 {//In bad channel (PHOS cristal size 2.2x2.2 cm), EMCAL ( cell units )
c4a7d28a 252 return kFALSE ;
253 }
254 else if(GetDebug() > 2) printf("\t Bad channel cut passed %4.2f > %2.2f \n",distBad, fMinDist);
fc195fd0 255
256 fhClusterCuts[7]->Fill(calo->E());
09273901 257
c4a7d28a 258 if(GetDebug() > 0)
259 printf("AliAnaPhoton::ClusterSelected() Current Event %d; After selection : E %2.2f, pT %2.2f, Ecl %2.2f, phi %2.2f, eta %2.2f\n",
260 GetReader()->GetEventNumber(),
fc195fd0 261 calo->E(), mom.Pt(),calo->E(),mom.Phi()*TMath::RadToDeg(),mom.Eta());
c4a7d28a 262
263 //All checks passed, cluster selected
264 return kTRUE;
265
266}
267
34c16486 268//___________________________________________
269void AliAnaPhoton::FillAcceptanceHistograms()
270{
3d5d5078 271 //Fill acceptance histograms if MC data is available
272
34c16486 273 Double_t photonY = -100 ;
274 Double_t photonE = -1 ;
275 Double_t photonPt = -1 ;
276 Double_t photonPhi = 100 ;
277 Double_t photonEta = -1 ;
278
279 Int_t pdg = 0 ;
280 Int_t tag = 0 ;
281 Int_t mcIndex = 0 ;
282 Bool_t inacceptance = kFALSE;
283
284 if(GetReader()->ReadStack())
285 {
3d5d5078 286 AliStack * stack = GetMCStack();
34c16486 287 if(stack)
288 {
289 for(Int_t i=0 ; i<stack->GetNtrack(); i++)
290 {
3d5d5078 291 TParticle * prim = stack->Particle(i) ;
34c16486 292 pdg = prim->GetPdgCode();
3d5d5078 293 //printf("i %d, %s %d %s %d \n",i, stack->Particle(i)->GetName(), stack->Particle(i)->GetPdgCode(),
294 // prim->GetName(), prim->GetPdgCode());
295
34c16486 296 if(pdg == 22)
297 {
3d5d5078 298 // Get tag of this particle photon from fragmentation, decay, prompt ...
34c16486 299 tag = GetMCAnalysisUtils()->CheckOrigin(i,GetReader(), 0);
300 if(!GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
301 {
3d5d5078 302 //A conversion photon from a hadron, skip this kind of photon
34c16486 303 // printf("AliAnaPhoton::FillAcceptanceHistograms() - not a photon, weird!\n ");
304 // GetMCAnalysisUtils()->PrintMCTag(tag);
3d5d5078 305
306 return;
307 }
308
309 //Get photon kinematics
310 if(prim->Energy() == TMath::Abs(prim->Pz())) continue ; //Protection against floating point exception
311
34c16486 312 photonY = 0.5*TMath::Log((prim->Energy()-prim->Pz())/(prim->Energy()+prim->Pz())) ;
313 photonE = prim->Energy() ;
314 photonPt = prim->Pt() ;
315 photonPhi = TMath::RadToDeg()*prim->Phi() ;
3d5d5078 316 if(photonPhi < 0) photonPhi+=TMath::TwoPi();
34c16486 317 photonEta = prim->Eta() ;
3d5d5078 318
319 //Check if photons hit the Calorimeter
320 TLorentzVector lv;
321 prim->Momentum(lv);
34c16486 322 inacceptance = kFALSE;
323 if (fCalorimeter == "PHOS")
324 {
325 if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet())
326 {
3d5d5078 327 Int_t mod ;
328 Double_t x,z ;
329 if(GetPHOSGeometry()->ImpactOnEmc(prim,mod,z,x))
330 inacceptance = kTRUE;
331 if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
332 }
34c16486 333 else
334 {
3d5d5078 335 if(GetFiducialCut()->IsInFiducialCut(lv,fCalorimeter))
336 inacceptance = kTRUE ;
337 if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
338 }
339 }
34c16486 340 else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet())
341 {
342 if(GetEMCALGeometry())
343 {
3d5d5078 344 Int_t absID=0;
345
346 GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(prim->Eta(),prim->Phi(),absID);
347
348 if( absID >= 0)
349 inacceptance = kTRUE;
350
351 // if(GetEMCALGeometry()->Impact(phot1) && GetEMCALGeometry()->Impact(phot2))
352 // inacceptance = kTRUE;
353 if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
354 }
34c16486 355 else
356 {
3d5d5078 357 if(GetFiducialCut()->IsInFiducialCut(lv,fCalorimeter))
358 inacceptance = kTRUE ;
359 if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
360 }
361 } //In EMCAL
362
363 //Fill histograms
c5693f62 364 fhYPrimMC[kmcPPhoton]->Fill(photonPt, photonY) ;
3d5d5078 365 if(TMath::Abs(photonY) < 1.0)
366 {
c5693f62 367 fhEPrimMC [kmcPPhoton]->Fill(photonE ) ;
368 fhPtPrimMC [kmcPPhoton]->Fill(photonPt) ;
369 fhPhiPrimMC[kmcPPhoton]->Fill(photonE , photonPhi) ;
34c16486 370 fhYPrimMC [kmcPPhoton]->Fill(photonE , photonEta) ;
3d5d5078 371 }
34c16486 372 if(inacceptance)
373 {
374 fhEPrimMCAcc [kmcPPhoton]->Fill(photonE ) ;
375 fhPtPrimMCAcc [kmcPPhoton]->Fill(photonPt) ;
c5693f62 376 fhPhiPrimMCAcc[kmcPPhoton]->Fill(photonE , photonPhi) ;
34c16486 377 fhYPrimMCAcc [kmcPPhoton]->Fill(photonE , photonY) ;
3d5d5078 378 }//Accepted
379
380 //Origin of photon
c5693f62 381 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt) && fhEPrimMC[kmcPPrompt])
3d5d5078 382 {
34c16486 383 mcIndex = kmcPPrompt;
3d5d5078 384 }
c5693f62 385 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation) && fhEPrimMC[kmcPFragmentation])
3d5d5078 386 {
34c16486 387 mcIndex = kmcPFragmentation ;
3d5d5078 388 }
c5693f62 389 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR) && fhEPrimMC[kmcPISR])
3d5d5078 390 {
34c16486 391 mcIndex = kmcPISR;
3d5d5078 392 }
c5693f62 393 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay)&& fhEPrimMC[kmcPPi0Decay])
3d5d5078 394 {
34c16486 395 mcIndex = kmcPPi0Decay;
3d5d5078 396 }
f586f4aa 397 else if( (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) ||
c5693f62 398 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay)) && fhEPrimMC[kmcPOtherDecay])
3d5d5078 399 {
34c16486 400 mcIndex = kmcPOtherDecay;
3d5d5078 401 }
c5693f62 402 else if(fhEPrimMC[kmcPOther])
3d5d5078 403 {
34c16486 404 mcIndex = kmcPOther;
3d5d5078 405 }//Other origin
34c16486 406
407 fhYPrimMC[mcIndex]->Fill(photonPt, photonY) ;
408 if(TMath::Abs(photonY) < 1.0)
409 {
410 fhEPrimMC [mcIndex]->Fill(photonE ) ;
411 fhPtPrimMC [mcIndex]->Fill(photonPt) ;
412 fhPhiPrimMC[mcIndex]->Fill(photonE , photonPhi) ;
413 fhYPrimMC [mcIndex]->Fill(photonE , photonEta) ;
414 }
415 if(inacceptance)
416 {
417 fhEPrimMCAcc [mcIndex]->Fill(photonE ) ;
418 fhPtPrimMCAcc [mcIndex]->Fill(photonPt) ;
419 fhPhiPrimMCAcc[mcIndex]->Fill(photonE , photonPhi) ;
420 fhYPrimMCAcc [mcIndex]->Fill(photonE , photonY) ;
421 }//Accepted
422
3d5d5078 423 }// Primary photon
424 }//loop on primaries
425 }//stack exists and data is MC
426 }//read stack
34c16486 427 else if(GetReader()->ReadAODMCParticles())
428 {
3d5d5078 429 TClonesArray * mcparticles = GetReader()->GetAODMCParticles(0);
34c16486 430 if(mcparticles)
431 {
3d5d5078 432 Int_t nprim = mcparticles->GetEntriesFast();
433
434 for(Int_t i=0; i < nprim; i++)
435 {
436 AliAODMCParticle * prim = (AliAODMCParticle *) mcparticles->At(i);
437
34c16486 438 pdg = prim->GetPdgCode();
3d5d5078 439
34c16486 440 if(pdg == 22)
441 {
3d5d5078 442 // Get tag of this particle photon from fragmentation, decay, prompt ...
34c16486 443 tag = GetMCAnalysisUtils()->CheckOrigin(i,GetReader(), 0);
444 if(!GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
445 {
3d5d5078 446 //A conversion photon from a hadron, skip this kind of photon
34c16486 447 // printf("AliAnaPhoton::FillAcceptanceHistograms() - not a photon, weird!\n ");
448 // GetMCAnalysisUtils()->PrintMCTag(tag);
3d5d5078 449
450 return;
451 }
452
453 //Get photon kinematics
454 if(prim->E() == TMath::Abs(prim->Pz())) continue ; //Protection against floating point exception
455
34c16486 456 photonY = 0.5*TMath::Log((prim->E()-prim->Pz())/(prim->E()+prim->Pz())) ;
457 photonE = prim->E() ;
458 photonPt = prim->Pt() ;
459 photonPhi = TMath::RadToDeg()*prim->Phi() ;
3d5d5078 460 if(photonPhi < 0) photonPhi+=TMath::TwoPi();
34c16486 461 photonEta = prim->Eta() ;
3d5d5078 462
463 //Check if photons hit the Calorimeter
464 TLorentzVector lv;
465 lv.SetPxPyPzE(prim->Px(),prim->Py(),prim->Pz(),prim->E());
34c16486 466 inacceptance = kFALSE;
467 if (fCalorimeter == "PHOS")
468 {
469 if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet())
470 {
3d5d5078 471 Int_t mod ;
472 Double_t x,z ;
473 Double_t vtx[]={prim->Xv(),prim->Yv(),prim->Zv()};
474 if(GetPHOSGeometry()->ImpactOnEmc(vtx, prim->Theta(),prim->Phi(),mod,z,x))
475 inacceptance = kTRUE;
476 if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
477 }
34c16486 478 else
479 {
3d5d5078 480 if(GetFiducialCut()->IsInFiducialCut(lv,fCalorimeter))
481 inacceptance = kTRUE ;
482 if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
483 }
484 }
34c16486 485 else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet())
486 {
487 if(GetEMCALGeometry())
488 {
3d5d5078 489 Int_t absID=0;
490
491 GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(prim->Eta(),prim->Phi(),absID);
492
493 if( absID >= 0)
494 inacceptance = kTRUE;
495
496 if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
497 }
34c16486 498 else
499 {
3d5d5078 500 if(GetFiducialCut()->IsInFiducialCut(lv,fCalorimeter))
501 inacceptance = kTRUE ;
502 if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
503 }
504 } //In EMCAL
505
506 //Fill histograms
507
c5693f62 508 fhYPrimMC[kmcPPhoton]->Fill(photonPt, photonY) ;
3d5d5078 509 if(TMath::Abs(photonY) < 1.0)
510 {
c5693f62 511 fhEPrimMC [kmcPPhoton]->Fill(photonE ) ;
512 fhPtPrimMC [kmcPPhoton]->Fill(photonPt) ;
513 fhPhiPrimMC[kmcPPhoton]->Fill(photonE , photonPhi) ;
514 fhYPrimMC[kmcPPhoton]->Fill(photonE , photonEta) ;
3d5d5078 515 }
34c16486 516
517 if(inacceptance)
518 {
c5693f62 519 fhEPrimMCAcc[kmcPPhoton] ->Fill(photonE ) ;
520 fhPtPrimMCAcc[kmcPPhoton] ->Fill(photonPt) ;
521 fhPhiPrimMCAcc[kmcPPhoton]->Fill(photonE , photonPhi) ;
522 fhYPrimMCAcc[kmcPPhoton] ->Fill(photonE , photonY) ;
3d5d5078 523 }//Accepted
524
3d5d5078 525 //Origin of photon
c5693f62 526 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt) && fhEPrimMC[kmcPPrompt])
3d5d5078 527 {
34c16486 528 mcIndex = kmcPPrompt;
3d5d5078 529 }
34c16486 530 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation) && fhEPrimMC[kmcPFragmentation])
3d5d5078 531 {
34c16486 532 mcIndex = kmcPFragmentation ;
3d5d5078 533 }
c5693f62 534 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR) && fhEPrimMC[kmcPISR])
3d5d5078 535 {
34c16486 536 mcIndex = kmcPISR;
3d5d5078 537 }
c5693f62 538 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay)&& fhEPrimMC[kmcPPi0Decay])
3d5d5078 539 {
34c16486 540 mcIndex = kmcPPi0Decay;
3d5d5078 541 }
34c16486 542 else if( (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) ||
543 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay)) && fhEPrimMC[kmcPOtherDecay])
3d5d5078 544 {
34c16486 545 mcIndex = kmcPOtherDecay;
3d5d5078 546 }
c5693f62 547 else if(fhEPrimMC[kmcPOther])
3d5d5078 548 {
34c16486 549 mcIndex = kmcPOther;
3d5d5078 550 }//Other origin
34c16486 551
552 fhYPrimMC[mcIndex]->Fill(photonPt, photonY) ;
553 if(TMath::Abs(photonY) < 1.0)
554 {
555 fhEPrimMC [mcIndex]->Fill(photonE ) ;
556 fhPtPrimMC [mcIndex]->Fill(photonPt) ;
557 fhPhiPrimMC[mcIndex]->Fill(photonE , photonPhi) ;
558 fhYPrimMC [mcIndex]->Fill(photonE , photonEta) ;
559 }
560 if(inacceptance)
561 {
562 fhEPrimMCAcc [mcIndex]->Fill(photonE ) ;
563 fhPtPrimMCAcc [mcIndex]->Fill(photonPt) ;
564 fhPhiPrimMCAcc[mcIndex]->Fill(photonE , photonPhi) ;
565 fhYPrimMCAcc [mcIndex]->Fill(photonE , photonY) ;
566 }//Accepted
567
3d5d5078 568 }// Primary photon
569 }//loop on primaries
570
c5693f62 571 }//kmc array exists and data is MC
3d5d5078 572 } // read AOD MC
573}
521636d2 574
34c16486 575//____________________________________________________________________________________
576void AliAnaPhoton::FillShowerShapeHistograms(AliVCluster* cluster, const Int_t mcTag)
577{
578 //Fill cluster Shower Shape histograms
521636d2 579
580 if(!fFillSSHistograms || GetMixedEvent()) return;
8d6b7f60 581
521636d2 582 Float_t energy = cluster->E();
583 Int_t ncells = cluster->GetNCells();
521636d2 584 Float_t lambda0 = cluster->GetM02();
585 Float_t lambda1 = cluster->GetM20();
586 Float_t disp = cluster->GetDispersion()*cluster->GetDispersion();
587
588 TLorentzVector mom;
34c16486 589 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
590 {
591 cluster->GetMomentum(mom,GetVertex(0)) ;
592 }//Assume that come from vertex in straight line
593 else
594 {
521636d2 595 Double_t vertex[]={0,0,0};
596 cluster->GetMomentum(mom,vertex) ;
597 }
598
599 Float_t eta = mom.Eta();
600 Float_t phi = mom.Phi();
601 if(phi < 0) phi+=TMath::TwoPi();
602
603 fhLam0E ->Fill(energy,lambda0);
604 fhLam1E ->Fill(energy,lambda1);
605 fhDispE ->Fill(energy,disp);
b5dbb99b 606
34c16486 607 if(fCalorimeter == "EMCAL" && GetModuleNumber(cluster) > 5)
608 {
521636d2 609 fhLam0ETRD->Fill(energy,lambda0);
610 fhLam1ETRD->Fill(energy,lambda1);
611 fhDispETRD->Fill(energy,disp);
521636d2 612 }
613
34c16486 614 Float_t l0 = 0., l1 = 0.;
615 Float_t dispp= 0., dEta = 0., dPhi = 0.;
616 Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
617 if(fCalorimeter == "EMCAL")
618 {
619 GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), GetReader()->GetInputEvent()->GetEMCALCells(), cluster,
620 l0, l1, dispp, dEta, dPhi, sEta, sPhi, sEtaPhi);
621 //printf("AliAnaPhoton::FillShowerShapeHistogram - l0 %2.6f, l1 %2.6f, disp %2.6f, dEta %2.6f, dPhi %2.6f, sEta %2.6f, sPhi %2.6f, sEtaPhi %2.6f \n",
622 // l0, l1, dispp, dEta, dPhi, sEta, sPhi, sEtaPhi );
623 //printf("AliAnaPhoton::FillShowerShapeHistogram - dispersion %f, dispersion eta+phi %f \n",
624 // disp, dPhi+dEta );
625 fhDispEtaE -> Fill(energy,dEta);
626 fhDispPhiE -> Fill(energy,dPhi);
627 fhSumEtaE -> Fill(energy,sEta);
628 fhSumPhiE -> Fill(energy,sPhi);
629 fhSumEtaPhiE -> Fill(energy,sEtaPhi);
630 fhDispEtaPhiDiffE -> Fill(energy,dPhi-dEta);
631 if(dEta+dPhi>0)fhSphericityE -> Fill(energy,(dPhi-dEta)/(dEta+dPhi));
632 if(dEta+sEta>0)fhDispSumEtaDiffE -> Fill(energy,(dEta-sEta)/((dEta+sEta)/2.));
633 if(dPhi+sPhi>0)fhDispSumPhiDiffE -> Fill(energy,(dPhi-sPhi)/((dPhi+sPhi)/2.));
634
bfdcf7fb 635 Int_t ebin = -1;
636 if (energy < 2 ) ebin = 0;
637 else if (energy < 4 ) ebin = 1;
638 else if (energy < 6 ) ebin = 2;
639 else if (energy < 10) ebin = 3;
640 else ebin = 4;
641
642 fhDispEtaDispPhi[ebin]->Fill(dEta ,dPhi);
643 fhLambda0DispEta[ebin]->Fill(lambda0,dEta);
644 fhLambda0DispPhi[ebin]->Fill(lambda0,dPhi);
645
34c16486 646 }
647
b5dbb99b 648 // if track-matching was of, check effect of track-matching residual cut
649
650 if(!fRejectTrackMatch)
651 {
652 Float_t dZ = cluster->GetTrackDz();
653 Float_t dR = cluster->GetTrackDx();
34c16486 654 if(cluster->IsEMCAL() && GetCaloUtils()->IsRecalculationOfClusterTrackMatchingOn())
655 {
b5dbb99b 656 dR = 2000., dZ = 2000.;
657 GetCaloUtils()->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
658 }
659
660 if(TMath::Abs(dZ) < 0.05 && TMath::Abs(dR) < 0.05)
661 {
662 fhLam0ETM ->Fill(energy,lambda0);
663 fhLam1ETM ->Fill(energy,lambda1);
664 fhDispETM ->Fill(energy,disp);
665
34c16486 666 if(fCalorimeter == "EMCAL" && GetModuleNumber(cluster) > 5)
667 {
b5dbb99b 668 fhLam0ETMTRD->Fill(energy,lambda0);
669 fhLam1ETMTRD->Fill(energy,lambda1);
670 fhDispETMTRD->Fill(energy,disp);
671 }
672 }
673 }// if track-matching was of, check effect of matching residual cut
674
34c16486 675 if(energy < 2)
676 {
521636d2 677 fhNCellsLam0LowE ->Fill(ncells,lambda0);
678 fhNCellsLam1LowE ->Fill(ncells,lambda1);
679 fhNCellsDispLowE ->Fill(ncells,disp);
521636d2 680
681 fhLam1Lam0LowE ->Fill(lambda1,lambda0);
682 fhLam0DispLowE ->Fill(lambda0,disp);
683 fhDispLam1LowE ->Fill(disp,lambda1);
684 fhEtaLam0LowE ->Fill(eta,lambda0);
685 fhPhiLam0LowE ->Fill(phi,lambda0);
521636d2 686 }
34c16486 687 else
688 {
521636d2 689 fhNCellsLam0HighE ->Fill(ncells,lambda0);
690 fhNCellsLam1HighE ->Fill(ncells,lambda1);
691 fhNCellsDispHighE ->Fill(ncells,disp);
7c65ad18 692
521636d2 693 fhLam1Lam0HighE ->Fill(lambda1,lambda0);
694 fhLam0DispHighE ->Fill(lambda0,disp);
695 fhDispLam1HighE ->Fill(disp,lambda1);
696 fhEtaLam0HighE ->Fill(eta, lambda0);
697 fhPhiLam0HighE ->Fill(phi, lambda0);
521636d2 698 }
699
34c16486 700 if(IsDataMC())
701 {
f66d95af 702 AliVCaloCells* cells = 0;
703 if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
704 else cells = GetPHOSCells();
3d5d5078 705
706 //Fill histograms to check shape of embedded clusters
707 Float_t fraction = 0;
34c16486 708 if(GetReader()->IsEmbeddedClusterSelectionOn())
709 {//Only working for EMCAL
3d5d5078 710 Float_t clusterE = 0; // recalculate in case corrections applied.
711 Float_t cellE = 0;
34c16486 712 for(Int_t icell = 0; icell < cluster->GetNCells(); icell++)
713 {
3d5d5078 714 cellE = cells->GetCellAmplitude(cluster->GetCellAbsId(icell));
715 clusterE+=cellE;
716 fraction+=cellE*cluster->GetCellAmplitudeFraction(icell);
717 }
718
719 //Fraction of total energy due to the embedded signal
720 fraction/=clusterE;
721
8d6b7f60 722 if(GetDebug() > 1 )
723 printf("AliAnaPhoton::FillShowerShapeHistogram() - Energy fraction of embedded signal %2.3f, Energy %2.3f\n",fraction, clusterE);
3d5d5078 724
725 fhEmbeddedSignalFractionEnergy->Fill(clusterE,fraction);
726
727 } // embedded fraction
728
f66d95af 729 // Get the fraction of the cluster energy that carries the cell with highest energy
34c16486 730 Int_t absID =-1 ;
f66d95af 731 Float_t maxCellFraction = 0.;
732
733 absID = GetCaloUtils()->GetMaxEnergyCell(cells, cluster,maxCellFraction);
734
735 // Check the origin and fill histograms
34c16486 736
737 Int_t mcIndex = -1;
738
739 if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton) &&
3d5d5078 740 !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion) &&
34c16486 741 !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0) &&
742 !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta))
743 {
744 mcIndex = kmcssPhoton ;
745
746 if(!GetReader()->IsEmbeddedClusterSelectionOn())
747 {
3d5d5078 748 //Check particle overlaps in cluster
749
8d6b7f60 750 // Compare the primary depositing more energy with the rest,
751 // if no photon/electron as comon ancestor (conversions), count as other particle
3d5d5078 752 Int_t ancPDG = 0, ancStatus = -1;
753 TLorentzVector momentum; TVector3 prodVertex;
754 Int_t ancLabel = 0;
755 Int_t noverlaps = 1;
34c16486 756 for (UInt_t ilab = 0; ilab < cluster->GetNLabels(); ilab++ )
757 {
758 ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(cluster->GetLabels()[0],cluster->GetLabels()[ilab],
759 GetReader(),ancPDG,ancStatus,momentum,prodVertex);
3d5d5078 760 if(ancPDG!=22 && TMath::Abs(ancPDG)!=11) noverlaps++;
761 }
8d6b7f60 762 //printf("N overlaps %d \n",noverlaps);
3d5d5078 763
34c16486 764 if(noverlaps == 1)
765 {
3d5d5078 766 fhMCPhotonELambda0NoOverlap ->Fill(energy, lambda0);
3d5d5078 767 }
34c16486 768 else if(noverlaps == 2)
769 {
3d5d5078 770 fhMCPhotonELambda0TwoOverlap ->Fill(energy, lambda0);
3d5d5078 771 }
34c16486 772 else if(noverlaps > 2)
773 {
3d5d5078 774 fhMCPhotonELambda0NOverlap ->Fill(energy, lambda0);
3d5d5078 775 }
34c16486 776 else
777 {
3d5d5078 778 printf("AliAnaPhoton::FillShowerShapeHistogram() - n overlaps = %d!!", noverlaps);
779 }
780 }//No embedding
781
782 //Fill histograms to check shape of embedded clusters
34c16486 783 if(GetReader()->IsEmbeddedClusterSelectionOn())
784 {
3d5d5078 785 if (fraction > 0.9)
786 {
787 fhEmbedPhotonELambda0FullSignal ->Fill(energy, lambda0);
3d5d5078 788 }
789 else if(fraction > 0.5)
790 {
791 fhEmbedPhotonELambda0MostlySignal ->Fill(energy, lambda0);
3d5d5078 792 }
793 else if(fraction > 0.1)
794 {
795 fhEmbedPhotonELambda0MostlyBkg ->Fill(energy, lambda0);
3d5d5078 796 }
797 else
798 {
799 fhEmbedPhotonELambda0FullBkg ->Fill(energy, lambda0);
3d5d5078 800 }
801 } // embedded
802
521636d2 803 }//photon no conversion
34c16486 804 else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron))
805 {
806 mcIndex = kmcssElectron ;
521636d2 807 }//electron
3d5d5078 808 else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton) &&
34c16486 809 GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion) )
810 {
811 mcIndex = kmcssConversion ;
521636d2 812 }//conversion photon
34c16486 813 else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0) )
814 {
815 mcIndex = kmcssPi0 ;
3d5d5078 816
817 //Fill histograms to check shape of embedded clusters
34c16486 818 if(GetReader()->IsEmbeddedClusterSelectionOn())
819 {
3d5d5078 820 if (fraction > 0.9)
821 {
822 fhEmbedPi0ELambda0FullSignal ->Fill(energy, lambda0);
3d5d5078 823 }
824 else if(fraction > 0.5)
825 {
826 fhEmbedPi0ELambda0MostlySignal ->Fill(energy, lambda0);
3d5d5078 827 }
828 else if(fraction > 0.1)
829 {
830 fhEmbedPi0ELambda0MostlyBkg ->Fill(energy, lambda0);
3d5d5078 831 }
832 else
833 {
834 fhEmbedPi0ELambda0FullBkg ->Fill(energy, lambda0);
3d5d5078 835 }
836 } // embedded
837
521636d2 838 }//pi0
34c16486 839 else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta) )
840 {
841 mcIndex = kmcssEta ;
3d5d5078 842 }//eta
34c16486 843 else
844 {
845 mcIndex = kmcssOther ;
521636d2 846 }//other particles
847
34c16486 848 fhMCELambda0 [mcIndex]->Fill(energy, lambda0);
849 fhMCELambda1 [mcIndex]->Fill(energy, lambda1);
850 fhMCEDispersion [mcIndex]->Fill(energy, disp);
851 fhMCNCellsE [mcIndex]->Fill(energy, ncells);
852 fhMCMaxCellDiffClusterE[mcIndex]->Fill(energy, maxCellFraction);
853
854 if (energy < 2.)
855 {
856 fhMCLambda0vsClusterMaxCellDiffE0[mcIndex]->Fill(lambda0, maxCellFraction);
857 fhMCNCellsvsClusterMaxCellDiffE0 [mcIndex]->Fill(ncells, maxCellFraction);
858 }
859 else if(energy < 6.)
860 {
861 fhMCLambda0vsClusterMaxCellDiffE2[mcIndex]->Fill(lambda0, maxCellFraction);
862 fhMCNCellsvsClusterMaxCellDiffE2 [mcIndex]->Fill(ncells, maxCellFraction);
863 }
864 else
865 {
866 fhMCLambda0vsClusterMaxCellDiffE6[mcIndex]->Fill(lambda0, maxCellFraction);
867 fhMCNCellsvsClusterMaxCellDiffE6 [mcIndex]->Fill(ncells, maxCellFraction);
868 }
869
870 if(fCalorimeter == "EMCAL")
871 {
872 fhMCEDispEta [mcIndex]-> Fill(energy,dEta);
873 fhMCEDispPhi [mcIndex]-> Fill(energy,dPhi);
874 fhMCESumEtaPhi [mcIndex]-> Fill(energy,sEtaPhi);
875 fhMCEDispEtaPhiDiff [mcIndex]-> Fill(energy,dPhi-dEta);
bfdcf7fb 876 if(dEta+dPhi>0)fhMCESphericity[mcIndex]-> Fill(energy,(dPhi-dEta)/(dEta+dPhi));
877
878 Int_t ebin = -1;
879 if (energy < 2 ) ebin = 0;
880 else if (energy < 4 ) ebin = 1;
881 else if (energy < 6 ) ebin = 2;
882 else if (energy < 10) ebin = 3;
883 else ebin = 4;
884
885 fhMCDispEtaDispPhi[ebin][mcIndex]->Fill(dEta ,dPhi);
886 fhMCLambda0DispEta[ebin][mcIndex]->Fill(lambda0,dEta);
887 fhMCLambda0DispPhi[ebin][mcIndex]->Fill(lambda0,dPhi);
34c16486 888 }
889
521636d2 890 }//MC data
891
892}
893
4bfeae64 894//__________________________________________________________________________
895void AliAnaPhoton::FillTrackMatchingResidualHistograms(AliVCluster* cluster,
896 const Int_t cut)
897{
898 // If selected, fill histograms with residuals of matched clusters, help to define track matching cut
899 // Residual filled for different cuts 0 (No cut), after 1 PID cut
900
901 Float_t dZ = cluster->GetTrackDz();
902 Float_t dR = cluster->GetTrackDx();
903
904 if(cluster->IsEMCAL() && GetCaloUtils()->IsRecalculationOfClusterTrackMatchingOn())
905 {
906 dR = 2000., dZ = 2000.;
907 GetCaloUtils()->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
908 }
909
910 if(fhTrackMatchedDEta[cut] && TMath::Abs(dR) < 999)
911 {
912 fhTrackMatchedDEta[cut]->Fill(cluster->E(),dZ);
913 fhTrackMatchedDPhi[cut]->Fill(cluster->E(),dR);
914
915 if(cluster->E() > 0.5) fhTrackMatchedDEtaDPhi[cut]->Fill(dZ,dR);
916
917 Int_t nSMod = GetModuleNumber(cluster);
918
919 if(fCalorimeter=="EMCAL" && nSMod > 5)
920 {
921 fhTrackMatchedDEtaTRD[cut]->Fill(cluster->E(),dZ);
922 fhTrackMatchedDPhiTRD[cut]->Fill(cluster->E(),dR);
923 }
924
925 // Check dEdx and E/p of matched clusters
926
927 if(TMath::Abs(dZ) < 0.05 && TMath::Abs(dR) < 0.05)
928 {
929
930 AliVTrack *track = GetCaloUtils()->GetMatchedTrack(cluster, GetReader()->GetInputEvent());
931
932 if(track)
933 {
934
935 Float_t dEdx = track->GetTPCsignal();
936 Float_t eOverp = cluster->E()/track->P();
937
938 fhdEdx[cut] ->Fill(cluster->E(), dEdx);
939 fhEOverP[cut]->Fill(cluster->E(), eOverp);
940
941 if(fCalorimeter=="EMCAL" && nSMod > 5)
942 fhEOverPTRD[cut]->Fill(cluster->E(), eOverp);
943
944
945 }
946 else
947 printf("AliAnaPhoton::FillTrackMatchingResidualHistograms() - Residual OK but (dR, dZ)= (%2.4f,%2.4f) no track associated WHAT? \n", dR,dZ);
948
949
950
951 if(IsDataMC())
952 {
953
954 Int_t tag = GetMCAnalysisUtils()->CheckOrigin(cluster->GetLabels(),cluster->GetNLabels(),GetReader(), 0);
955
956 if ( !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) )
957 {
958 if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) ||
959 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) ) fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 2.5 );
960 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) ) fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 0.5 );
961 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron) ) fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 1.5 );
962 else fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 3.5 );
963
964 // Check if several particles contributed to cluster and discard overlapped mesons
965 if(!GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) ||
34c16486 966 !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta))
967 {
4bfeae64 968 if(cluster->GetNLabels()==1)
969 {
970 fhTrackMatchedDEtaMCNoOverlap[cut]->Fill(cluster->E(),dZ);
971 fhTrackMatchedDPhiMCNoOverlap[cut]->Fill(cluster->E(),dR);
972 }
973 else
974 {
975 fhTrackMatchedDEtaMCOverlap[cut]->Fill(cluster->E(),dZ);
976 fhTrackMatchedDPhiMCOverlap[cut]->Fill(cluster->E(),dR);
977 }
978
979 }// Check overlaps
980
981 }
982 else
983 {
984 if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) ||
985 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) ) fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 6.5 );
986 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) ) fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 4.5 );
987 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron) ) fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 5.5 );
988 else fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 7.5 );
989
990 // Check if several particles contributed to cluster
991 if(!GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) ||
34c16486 992 !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta))
993 {
4bfeae64 994 fhTrackMatchedDEtaMCConversion[cut]->Fill(cluster->E(),dZ);
995 fhTrackMatchedDPhiMCConversion[cut]->Fill(cluster->E(),dR);
996
997 }// Check overlaps
998
999 }
1000
1001 } // MC
1002
1003 } // residuals window
1004
1005 } // Small residual
1006
1007}
1008
1009//___________________________________________
0c1383b5 1010TObjString * AliAnaPhoton::GetAnalysisCuts()
1011{
1012 //Save parameters used for analysis
1013 TString parList ; //this will be list of parameters used for this analysis.
5ae09196 1014 const Int_t buffersize = 255;
1015 char onePar[buffersize] ;
0c1383b5 1016
5ae09196 1017 snprintf(onePar,buffersize,"--- AliAnaPhoton ---\n") ;
0c1383b5 1018 parList+=onePar ;
5ae09196 1019 snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
0c1383b5 1020 parList+=onePar ;
5ae09196 1021 snprintf(onePar,buffersize,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
0c1383b5 1022 parList+=onePar ;
5ae09196 1023 snprintf(onePar,buffersize,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
0c1383b5 1024 parList+=onePar ;
5ae09196 1025 snprintf(onePar,buffersize,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
0c1383b5 1026 parList+=onePar ;
5ae09196 1027 snprintf(onePar,buffersize,"fRejectTrackMatch: %d\n",fRejectTrackMatch) ;
0c1383b5 1028 parList+=onePar ;
1029
1030 //Get parameters set in base class.
1031 parList += GetBaseParametersList() ;
1032
1033 //Get parameters set in PID class.
1034 parList += GetCaloPID()->GetPIDParametersList() ;
1035
1036 //Get parameters set in FiducialCut class (not available yet)
1037 //parlist += GetFidCut()->GetFidCutParametersList()
1038
1039 return new TObjString(parList) ;
1040}
1041
1c5acb87 1042//________________________________________________________________________
1043TList * AliAnaPhoton::GetCreateOutputObjects()
1044{
477d6cee 1045 // Create histograms to be saved in output file and
1046 // store them in outputContainer
1047 TList * outputContainer = new TList() ;
1048 outputContainer->SetName("PhotonHistos") ;
4a745797 1049
745913ae 1050 Int_t nptbins = GetHistogramRanges()->GetHistoPtBins(); Float_t ptmax = GetHistogramRanges()->GetHistoPtMax(); Float_t ptmin = GetHistogramRanges()->GetHistoPtMin();
1051 Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins(); Float_t phimax = GetHistogramRanges()->GetHistoPhiMax(); Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();
1052 Int_t netabins = GetHistogramRanges()->GetHistoEtaBins(); Float_t etamax = GetHistogramRanges()->GetHistoEtaMax(); Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
1053 Int_t ssbins = GetHistogramRanges()->GetHistoShowerShapeBins(); Float_t ssmax = GetHistogramRanges()->GetHistoShowerShapeMax(); Float_t ssmin = GetHistogramRanges()->GetHistoShowerShapeMin();
1054 Int_t nbins = GetHistogramRanges()->GetHistoNClusterCellBins(); Int_t nmax = GetHistogramRanges()->GetHistoNClusterCellMax(); Int_t nmin = GetHistogramRanges()->GetHistoNClusterCellMin();
1055 Int_t ntimebins= GetHistogramRanges()->GetHistoTimeBins(); Float_t timemax = GetHistogramRanges()->GetHistoTimeMax(); Float_t timemin = GetHistogramRanges()->GetHistoTimeMin();
521636d2 1056
09273901 1057 Int_t nresetabins = GetHistogramRanges()->GetHistoTrackResidualEtaBins();
1058 Float_t resetamax = GetHistogramRanges()->GetHistoTrackResidualEtaMax();
1059 Float_t resetamin = GetHistogramRanges()->GetHistoTrackResidualEtaMin();
1060 Int_t nresphibins = GetHistogramRanges()->GetHistoTrackResidualPhiBins();
1061 Float_t resphimax = GetHistogramRanges()->GetHistoTrackResidualPhiMax();
1062 Float_t resphimin = GetHistogramRanges()->GetHistoTrackResidualPhiMin();
1063
31ae6d59 1064 Int_t ndedxbins = GetHistogramRanges()->GetHistodEdxBins();
1065 Float_t dedxmax = GetHistogramRanges()->GetHistodEdxMax();
1066 Float_t dedxmin = GetHistogramRanges()->GetHistodEdxMin();
1067 Int_t nPoverEbins = GetHistogramRanges()->GetHistoPOverEBins();
1068 Float_t pOverEmax = GetHistogramRanges()->GetHistoPOverEMax();
1069 Float_t pOverEmin = GetHistogramRanges()->GetHistoPOverEMin();
09273901 1070
fc195fd0 1071 TString cut[] = {"Open","Reader","E","Time","NCells","Fidutial","Matching","Bad","PID"};
1072 for (Int_t i = 0; i < 9 ; i++)
1073 {
1074 fhClusterCuts[i] = new TH1F(Form("hCut_%d_%s", i, cut[i].Data()),
1075 Form("Number of clusters that pass cuts <= %d, %s", i, cut[i].Data()),
1076 nptbins,ptmin,ptmax);
1077 fhClusterCuts[i]->SetYTitle("dN/dE ");
1078 fhClusterCuts[i]->SetXTitle("E (GeV)");
1079 outputContainer->Add(fhClusterCuts[i]) ;
1080 }
1081
e1e62b89 1082 fhNCellsE = new TH2F ("hNCellsE","# of cells in cluster vs E of clusters", nptbins,ptmin,ptmax, nbins,nmin,nmax);
c4a7d28a 1083 fhNCellsE->SetXTitle("E (GeV)");
1084 fhNCellsE->SetYTitle("# of cells in cluster");
f15c25da 1085 outputContainer->Add(fhNCellsE);
1086
5c46c992 1087 fhCellsE = new TH2F ("hCellsE","energy of cells in cluster vs E of clusters", nptbins,ptmin,ptmax, nptbins*2,ptmin,ptmax);
1088 fhCellsE->SetXTitle("E_{cluster} (GeV)");
1089 fhCellsE->SetYTitle("E_{cell} (GeV)");
1090 outputContainer->Add(fhCellsE);
1091
f15c25da 1092 fhTimeE = new TH2F ("hTimeE","time of cluster vs E of clusters", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
1093 fhTimeE->SetXTitle("E (GeV)");
1094 fhTimeE->SetYTitle("time (ns)");
1095 outputContainer->Add(fhTimeE);
6175da48 1096
f66d95af 1097 fhMaxCellDiffClusterE = new TH2F ("hMaxCellDiffClusterE","energy vs difference of cluster energy - max cell energy / cluster energy, good clusters",
1098 nptbins,ptmin,ptmax, 500,0,1.);
1099 fhMaxCellDiffClusterE->SetXTitle("E_{cluster} (GeV) ");
1100 fhMaxCellDiffClusterE->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
1101 outputContainer->Add(fhMaxCellDiffClusterE);
1102
20218aea 1103 fhEPhoton = new TH1F("hEPhoton","Number of #gamma over calorimeter vs energy",nptbins,ptmin,ptmax);
1104 fhEPhoton->SetYTitle("N");
1105 fhEPhoton->SetXTitle("E_{#gamma}(GeV)");
1106 outputContainer->Add(fhEPhoton) ;
1107
1108 fhPtPhoton = new TH1F("hPtPhoton","Number of #gamma over calorimeter vs p_{T}",nptbins,ptmin,ptmax);
477d6cee 1109 fhPtPhoton->SetYTitle("N");
1110 fhPtPhoton->SetXTitle("p_{T #gamma}(GeV/c)");
1111 outputContainer->Add(fhPtPhoton) ;
1112
1113 fhPhiPhoton = new TH2F
20218aea 1114 ("hPhiPhoton","#phi_{#gamma} vs p_{T}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
6175da48 1115 fhPhiPhoton->SetYTitle("#phi (rad)");
477d6cee 1116 fhPhiPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
1117 outputContainer->Add(fhPhiPhoton) ;
1118
1119 fhEtaPhoton = new TH2F
20218aea 1120 ("hEtaPhoton","#eta_{#gamma} vs p_{T}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 1121 fhEtaPhoton->SetYTitle("#eta");
1122 fhEtaPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
1123 outputContainer->Add(fhEtaPhoton) ;
1124
6175da48 1125 fhEtaPhiPhoton = new TH2F
1126 ("hEtaPhiPhoton","#eta vs #phi",netabins,etamin,etamax,nphibins,phimin,phimax);
1127 fhEtaPhiPhoton->SetYTitle("#phi (rad)");
1128 fhEtaPhiPhoton->SetXTitle("#eta");
1129 outputContainer->Add(fhEtaPhiPhoton) ;
34c16486 1130 if(GetMinPt() < 0.5)
1131 {
20218aea 1132 fhEtaPhi05Photon = new TH2F
1133 ("hEtaPhi05Photon","#eta vs #phi, E > 0.5",netabins,etamin,etamax,nphibins,phimin,phimax);
1134 fhEtaPhi05Photon->SetYTitle("#phi (rad)");
1135 fhEtaPhi05Photon->SetXTitle("#eta");
1136 outputContainer->Add(fhEtaPhi05Photon) ;
1137 }
6175da48 1138
521636d2 1139 //Shower shape
34c16486 1140 if(fFillSSHistograms)
1141 {
521636d2 1142 fhLam0E = new TH2F ("hLam0E","#lambda_{0}^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1143 fhLam0E->SetYTitle("#lambda_{0}^{2}");
1144 fhLam0E->SetXTitle("E (GeV)");
1145 outputContainer->Add(fhLam0E);
1146
1147 fhLam1E = new TH2F ("hLam1E","#lambda_{1}^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1148 fhLam1E->SetYTitle("#lambda_{1}^{2}");
1149 fhLam1E->SetXTitle("E (GeV)");
1150 outputContainer->Add(fhLam1E);
1151
1152 fhDispE = new TH2F ("hDispE"," dispersion^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1153 fhDispE->SetYTitle("D^{2}");
1154 fhDispE->SetXTitle("E (GeV) ");
1155 outputContainer->Add(fhDispE);
b5dbb99b 1156
1157 if(!fRejectTrackMatch)
1158 {
1159 fhLam0ETM = new TH2F ("hLam0ETM","#lambda_{0}^{2} vs E, cut on track-matching residual |#Delta #eta| < 0.05, |#Delta #phi| < 0.05", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1160 fhLam0ETM->SetYTitle("#lambda_{0}^{2}");
1161 fhLam0ETM->SetXTitle("E (GeV)");
1162 outputContainer->Add(fhLam0ETM);
1163
1164 fhLam1ETM = new TH2F ("hLam1ETM","#lambda_{1}^{2} vs E, cut on track-matching residual |#Delta #eta| < 0.05, |#Delta #phi| < 0.05", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1165 fhLam1ETM->SetYTitle("#lambda_{1}^{2}");
1166 fhLam1ETM->SetXTitle("E (GeV)");
1167 outputContainer->Add(fhLam1ETM);
1168
1169 fhDispETM = new TH2F ("hDispETM"," dispersion^{2} vs E, cut on track-matching residual |#Delta #eta| < 0.05, |#Delta #phi| < 0.05", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1170 fhDispETM->SetYTitle("D^{2}");
1171 fhDispETM->SetXTitle("E (GeV) ");
1172 outputContainer->Add(fhDispETM);
1173 }
521636d2 1174
b5dbb99b 1175 if(fCalorimeter == "EMCAL")
1176 {
521636d2 1177 fhLam0ETRD = new TH2F ("hLam0ETRD","#lambda_{0}^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1178 fhLam0ETRD->SetYTitle("#lambda_{0}^{2}");
1179 fhLam0ETRD->SetXTitle("E (GeV)");
1180 outputContainer->Add(fhLam0ETRD);
1181
1182 fhLam1ETRD = new TH2F ("hLam1ETRD","#lambda_{1}^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1183 fhLam1ETRD->SetYTitle("#lambda_{1}^{2}");
1184 fhLam1ETRD->SetXTitle("E (GeV)");
1185 outputContainer->Add(fhLam1ETRD);
1186
1187 fhDispETRD = new TH2F ("hDispETRD"," dispersion^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1188 fhDispETRD->SetYTitle("Dispersion^{2}");
1189 fhDispETRD->SetXTitle("E (GeV) ");
b5dbb99b 1190 outputContainer->Add(fhDispETRD);
1191
1192 if(!fRejectTrackMatch)
1193 {
1194 fhLam0ETMTRD = new TH2F ("hLam0ETMTRD","#lambda_{0}^{2} vs E, EMCAL SM covered by TRD, cut on track-matching residual |#Delta #eta| < 0.05, |#Delta #phi| < 0.05", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1195 fhLam0ETMTRD->SetYTitle("#lambda_{0}^{2}");
1196 fhLam0ETMTRD->SetXTitle("E (GeV)");
1197 outputContainer->Add(fhLam0ETMTRD);
1198
1199 fhLam1ETMTRD = new TH2F ("hLam1ETMTRD","#lambda_{1}^{2} vs E, EMCAL SM covered by TRD, cut on track-matching residual |#Delta #eta| < 0.05, |#Delta #phi| < 0.05", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1200 fhLam1ETMTRD->SetYTitle("#lambda_{1}^{2}");
1201 fhLam1ETMTRD->SetXTitle("E (GeV)");
1202 outputContainer->Add(fhLam1ETMTRD);
1203
1204 fhDispETMTRD = new TH2F ("hDispETMTRD"," dispersion^{2} vs E, EMCAL SM covered by TRD, cut on track-matching residual |#Delta #eta| < 0.05, |#Delta #phi| < 0.05", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1205 fhDispETMTRD->SetYTitle("Dispersion^{2}");
1206 fhDispETMTRD->SetXTitle("E (GeV) ");
1207 outputContainer->Add(fhDispETMTRD);
1208 }
521636d2 1209 }
1210
d9105d92 1211 fhNCellsLam0LowE = new TH2F ("hNCellsLam0LowE","N_{cells} in cluster vs #lambda_{0}^{2}, E < 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
521636d2 1212 fhNCellsLam0LowE->SetXTitle("N_{cells}");
1213 fhNCellsLam0LowE->SetYTitle("#lambda_{0}^{2}");
1214 outputContainer->Add(fhNCellsLam0LowE);
1215
d9105d92 1216 fhNCellsLam0HighE = new TH2F ("hNCellsLam0HighE","N_{cells} in cluster vs #lambda_{0}^{2}, E > 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
521636d2 1217 fhNCellsLam0HighE->SetXTitle("N_{cells}");
1218 fhNCellsLam0HighE->SetYTitle("#lambda_{0}^{2}");
1219 outputContainer->Add(fhNCellsLam0HighE);
1220
d9105d92 1221 fhNCellsLam1LowE = new TH2F ("hNCellsLam1LowE","N_{cells} in cluster vs #lambda_{1}^{2}, E < 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
521636d2 1222 fhNCellsLam1LowE->SetXTitle("N_{cells}");
1223 fhNCellsLam1LowE->SetYTitle("#lambda_{0}^{2}");
1224 outputContainer->Add(fhNCellsLam1LowE);
1225
d9105d92 1226 fhNCellsLam1HighE = new TH2F ("hNCellsLam1HighE","N_{cells} in cluster vs #lambda_{1}^{2}, E > 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
521636d2 1227 fhNCellsLam1HighE->SetXTitle("N_{cells}");
1228 fhNCellsLam1HighE->SetYTitle("#lambda_{0}^{2}");
1229 outputContainer->Add(fhNCellsLam1HighE);
1230
d9105d92 1231 fhNCellsDispLowE = new TH2F ("hNCellsDispLowE","N_{cells} in cluster vs dispersion^{2}, E < 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
521636d2 1232 fhNCellsDispLowE->SetXTitle("N_{cells}");
1233 fhNCellsDispLowE->SetYTitle("D^{2}");
1234 outputContainer->Add(fhNCellsDispLowE);
1235
d9105d92 1236 fhNCellsDispHighE = new TH2F ("hNCellsDispHighE","N_{cells} in cluster vs dispersion^{2}, E < 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
521636d2 1237 fhNCellsDispHighE->SetXTitle("N_{cells}");
1238 fhNCellsDispHighE->SetYTitle("D^{2}");
1239 outputContainer->Add(fhNCellsDispHighE);
1240
521636d2 1241 fhEtaLam0LowE = new TH2F ("hEtaLam0LowE","#eta vs #lambda_{0}^{2}, E < 2 GeV", netabins,etamin,etamax, ssbins,ssmin,ssmax);
1242 fhEtaLam0LowE->SetYTitle("#lambda_{0}^{2}");
1243 fhEtaLam0LowE->SetXTitle("#eta");
1244 outputContainer->Add(fhEtaLam0LowE);
1245
1246 fhPhiLam0LowE = new TH2F ("hPhiLam0LowE","#phi vs #lambda_{0}^{2}, E < 2 GeV", nphibins,phimin,phimax, ssbins,ssmin,ssmax);
1247 fhPhiLam0LowE->SetYTitle("#lambda_{0}^{2}");
1248 fhPhiLam0LowE->SetXTitle("#phi");
1249 outputContainer->Add(fhPhiLam0LowE);
1250
1251 fhEtaLam0HighE = new TH2F ("hEtaLam0HighE","#eta vs #lambda_{0}^{2}, E > 2 GeV", netabins,etamin,etamax, ssbins,ssmin,ssmax);
1252 fhEtaLam0HighE->SetYTitle("#lambda_{0}^{2}");
1253 fhEtaLam0HighE->SetXTitle("#eta");
1254 outputContainer->Add(fhEtaLam0HighE);
1255
1256 fhPhiLam0HighE = new TH2F ("hPhiLam0HighE","#phi vs #lambda_{0}^{2}, E > 2 GeV", nphibins,phimin,phimax, ssbins,ssmin,ssmax);
1257 fhPhiLam0HighE->SetYTitle("#lambda_{0}^{2}");
1258 fhPhiLam0HighE->SetXTitle("#phi");
1259 outputContainer->Add(fhPhiLam0HighE);
1260
1261 fhLam1Lam0LowE = new TH2F ("hLam1Lam0LowE","#lambda_{0}^{2} vs #lambda_{1}^{2} in cluster of E < 2 GeV", ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
1262 fhLam1Lam0LowE->SetYTitle("#lambda_{0}^{2}");
1263 fhLam1Lam0LowE->SetXTitle("#lambda_{1}^{2}");
1264 outputContainer->Add(fhLam1Lam0LowE);
1265
1266 fhLam1Lam0HighE = new TH2F ("hLam1Lam0HighE","#lambda_{0}^{2} vs #lambda_{1}^{2} in cluster of E > 2 GeV", ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
1267 fhLam1Lam0HighE->SetYTitle("#lambda_{0}^{2}");
1268 fhLam1Lam0HighE->SetXTitle("#lambda_{1}^{2}");
1269 outputContainer->Add(fhLam1Lam0HighE);
1270
1271 fhLam0DispLowE = new TH2F ("hLam0DispLowE","#lambda_{0}^{2} vs dispersion^{2} in cluster of E < 2 GeV", ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
1272 fhLam0DispLowE->SetXTitle("#lambda_{0}^{2}");
1273 fhLam0DispLowE->SetYTitle("D^{2}");
1274 outputContainer->Add(fhLam0DispLowE);
1275
1276 fhLam0DispHighE = new TH2F ("hLam0DispHighE","#lambda_{0}^{2} vs dispersion^{2} in cluster of E > 2 GeV", ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
1277 fhLam0DispHighE->SetXTitle("#lambda_{0}^{2}");
1278 fhLam0DispHighE->SetYTitle("D^{2}");
1279 outputContainer->Add(fhLam0DispHighE);
1280
1281 fhDispLam1LowE = new TH2F ("hDispLam1LowE","Dispersion^{2} vs #lambda_{1}^{2} in cluster of E < 2 GeV", ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
1282 fhDispLam1LowE->SetXTitle("D^{2}");
1283 fhDispLam1LowE->SetYTitle("#lambda_{1}^{2}");
1284 outputContainer->Add(fhDispLam1LowE);
1285
1286 fhDispLam1HighE = new TH2F ("hDispLam1HighE","Dispersion^{2} vs #lambda_{1^{2}} in cluster of E > 2 GeV", ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
1287 fhDispLam1HighE->SetXTitle("D^{2}");
1288 fhDispLam1HighE->SetYTitle("#lambda_{1}^{2}");
1289 outputContainer->Add(fhDispLam1HighE);
1290
34c16486 1291 if(fCalorimeter == "EMCAL")
1292 {
1293 fhDispEtaE = new TH2F ("hDispEtaE","#sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs E", nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1294 fhDispEtaE->SetXTitle("E (GeV)");
1295 fhDispEtaE->SetYTitle("#sigma^{2}_{#eta #eta}");
1296 outputContainer->Add(fhDispEtaE);
1297
1298 fhDispPhiE = new TH2F ("hDispPhiE","#sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs E", nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1299 fhDispPhiE->SetXTitle("E (GeV)");
1300 fhDispPhiE->SetYTitle("#sigma^{2}_{#phi #phi}");
1301 outputContainer->Add(fhDispPhiE);
1302
1303 fhSumEtaE = new TH2F ("hSumEtaE","#sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i})^{2} / #Sigma w_{i} - <#eta>^{2} vs E", nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1304 fhSumEtaE->SetXTitle("E (GeV)");
1305 fhSumEtaE->SetYTitle("#sigma'^{2}_{#eta #eta}");
1306 outputContainer->Add(fhSumEtaE);
1307
1308 fhSumPhiE = new TH2F ("hSumPhiE","#sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i})^{2}/ #Sigma w_{i} - <#phi>^{2} vs E",
1309 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1310 fhSumPhiE->SetXTitle("E (GeV)");
1311 fhSumPhiE->SetYTitle("#sigma'^{2}_{#phi #phi}");
1312 outputContainer->Add(fhSumPhiE);
1313
1314 fhSumEtaPhiE = new TH2F ("hSumEtaPhiE","#sigma'^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs E",
1315 nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
1316 fhSumEtaPhiE->SetXTitle("E (GeV)");
1317 fhSumEtaPhiE->SetYTitle("#sigma'^{2}_{#eta #phi}");
1318 outputContainer->Add(fhSumEtaPhiE);
1319
1320 fhDispEtaPhiDiffE = new TH2F ("hDispEtaPhiDiffE","#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs E",
1321 nptbins,ptmin,ptmax,200, -10,10);
1322 fhDispEtaPhiDiffE->SetXTitle("E (GeV)");
1323 fhDispEtaPhiDiffE->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
1324 outputContainer->Add(fhDispEtaPhiDiffE);
1325
1326 fhSphericityE = new TH2F ("hSphericityE","(#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi}) vs E",
1327 nptbins,ptmin,ptmax, 200, -1,1);
1328 fhSphericityE->SetXTitle("E (GeV)");
1329 fhSphericityE->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
1330 outputContainer->Add(fhSphericityE);
1331
1332 fhDispSumEtaDiffE = new TH2F ("hDispSumEtaDiffE","#sigma^{2}_{#eta #eta} - #sigma'^{2}_{#eta #eta} / average vs E", nptbins,ptmin,ptmax, 200,-0.01,0.01);
1333 fhDispSumEtaDiffE->SetXTitle("E (GeV)");
1334 fhDispSumEtaDiffE->SetYTitle("#sigma^{2}_{#eta #eta} - #sigma'^{2}_{#eta #eta} / average");
1335 outputContainer->Add(fhDispSumEtaDiffE);
1336
1337 fhDispSumPhiDiffE = new TH2F ("hDispSumPhiDiffE","#sigma^{2}_{#phi #phi} - #sigma'^{2}_{#phi #phi} / average vs E", nptbins,ptmin,ptmax, 200,-0.01,0.01);
1338 fhDispSumPhiDiffE->SetXTitle("E (GeV)");
1339 fhDispSumPhiDiffE->SetYTitle("#sigma^{2}_{#phi #phi} - #sigma'^{2}_{#phi #phi} / average");
1340 outputContainer->Add(fhDispSumPhiDiffE);
1341
1342
1343 Int_t bin[] = {0,2,4,6,10,1000};
1344 for(Int_t i = 0; i < 5; i++)
1345 {
bfdcf7fb 1346 fhDispEtaDispPhi[i] = new TH2F (Form("hDispEtaDispPhi_EBin%d",i),Form("#sigma^{2}_{#phi #phi} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",bin[i],bin[i+1]),
34c16486 1347 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
bfdcf7fb 1348 fhDispEtaDispPhi[i]->SetXTitle("#sigma^{2}_{#eta #eta}");
1349 fhDispEtaDispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1350 outputContainer->Add(fhDispEtaDispPhi[i]);
1351
1352 fhLambda0DispEta[i] = new TH2F (Form("hLambda0DispEta_EBin%d",i),Form("#lambda^{2}_{0} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",bin[i],bin[i+1]),
1353 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1354 fhLambda0DispEta[i]->SetXTitle("#lambda^{2}_{0}");
1355 fhLambda0DispEta[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1356 outputContainer->Add(fhLambda0DispEta[i]);
1357
1358 fhLambda0DispPhi[i] = new TH2F (Form("hLambda0DispPhi_EBin%d",i),Form("#lambda^{2}_{0}} vs #sigma^{2}_{#phi #phi} for %d < E < %d GeV",bin[i],bin[i+1]),
1359 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1360 fhLambda0DispPhi[i]->SetXTitle("#lambda^{2}_{0}");
1361 fhLambda0DispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1362 outputContainer->Add(fhLambda0DispPhi[i]);
34c16486 1363 }
1364 }
521636d2 1365 } // Shower shape
1366
09273901 1367 // Track Matching
1368
b5dbb99b 1369 if(fFillTMHisto)
1370 {
4bfeae64 1371 fhTrackMatchedDEta[0] = new TH2F
1372 ("hTrackMatchedDEtaNoCut",
1373 "d#eta of cluster-track vs cluster energy, no photon cuts",
09273901 1374 nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
4bfeae64 1375 fhTrackMatchedDEta[0]->SetYTitle("d#eta");
1376 fhTrackMatchedDEta[0]->SetXTitle("E_{cluster} (GeV)");
09273901 1377
4bfeae64 1378 fhTrackMatchedDPhi[0] = new TH2F
1379 ("hTrackMatchedDPhiNoCut",
1380 "d#phi of cluster-track vs cluster energy, no photon cuts",
09273901 1381 nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
4bfeae64 1382 fhTrackMatchedDPhi[0]->SetYTitle("d#phi (rad)");
1383 fhTrackMatchedDPhi[0]->SetXTitle("E_{cluster} (GeV)");
09273901 1384
4bfeae64 1385 fhTrackMatchedDEtaDPhi[0] = new TH2F
1386 ("hTrackMatchedDEtaDPhiNoCut",
1387 "d#eta vs d#phi of cluster-track vs cluster energy, no photon cuts",
09273901 1388 nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
4bfeae64 1389 fhTrackMatchedDEtaDPhi[0]->SetYTitle("d#phi (rad)");
1390 fhTrackMatchedDEtaDPhi[0]->SetXTitle("d#eta");
b5dbb99b 1391
4bfeae64 1392 fhdEdx[0] = new TH2F ("hdEdxNoCut","matched track <dE/dx> vs cluster E, no photon cuts ",
1393 nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
1394 fhdEdx[0]->SetXTitle("E (GeV)");
1395 fhdEdx[0]->SetYTitle("<dE/dx>");
1396
1397 fhEOverP[0] = new TH2F ("hEOverPNoCut","matched track E/p vs cluster E, no photon cuts ",
1398 nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
1399 fhEOverP[0]->SetXTitle("E (GeV)");
1400 fhEOverP[0]->SetYTitle("E/p");
1401
1402 outputContainer->Add(fhTrackMatchedDEta[0]) ;
1403 outputContainer->Add(fhTrackMatchedDPhi[0]) ;
1404 outputContainer->Add(fhTrackMatchedDEtaDPhi[0]) ;
1405 outputContainer->Add(fhdEdx[0]);
1406 outputContainer->Add(fhEOverP[0]);
1407
1408 fhTrackMatchedDEta[1] = new TH2F
1409 ("hTrackMatchedDEta",
09273901 1410 "d#eta of cluster-track vs cluster energy, no photon cuts",
1411 nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
4bfeae64 1412 fhTrackMatchedDEta[1]->SetYTitle("d#eta");
1413 fhTrackMatchedDEta[1]->SetXTitle("E_{cluster} (GeV)");
09273901 1414
4bfeae64 1415 fhTrackMatchedDPhi[1] = new TH2F
1416 ("hTrackMatchedDPhi",
09273901 1417 "d#phi of cluster-track vs cluster energy, no photon cuts",
1418 nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
4bfeae64 1419 fhTrackMatchedDPhi[1]->SetYTitle("d#phi (rad)");
1420 fhTrackMatchedDPhi[1]->SetXTitle("E_{cluster} (GeV)");
09273901 1421
4bfeae64 1422 fhTrackMatchedDEtaDPhi[1] = new TH2F
1423 ("hTrackMatchedDEtaDPhi",
09273901 1424 "d#eta vs d#phi of cluster-track vs cluster energy, no photon cuts",
1425 nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
4bfeae64 1426 fhTrackMatchedDEtaDPhi[1]->SetYTitle("d#phi (rad)");
1427 fhTrackMatchedDEtaDPhi[1]->SetXTitle("d#eta");
09273901 1428
4bfeae64 1429 fhdEdx[1] = new TH2F ("hdEdx","matched track <dE/dx> vs cluster E ",
1430 nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
1431 fhdEdx[1]->SetXTitle("E (GeV)");
1432 fhdEdx[1]->SetYTitle("<dE/dx>");
1433
1434 fhEOverP[1] = new TH2F ("hEOverP","matched track E/p vs cluster E ",
1435 nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
1436 fhEOverP[1]->SetXTitle("E (GeV)");
1437 fhEOverP[1]->SetYTitle("E/p");
1438
1439 outputContainer->Add(fhTrackMatchedDEta[1]) ;
1440 outputContainer->Add(fhTrackMatchedDPhi[1]) ;
1441 outputContainer->Add(fhTrackMatchedDEtaDPhi[1]) ;
1442 outputContainer->Add(fhdEdx[1]);
1443 outputContainer->Add(fhEOverP[1]);
31ae6d59 1444
b5dbb99b 1445 if(fCalorimeter=="EMCAL")
1446 {
4bfeae64 1447 fhTrackMatchedDEtaTRD[0] = new TH2F
b5dbb99b 1448 ("hTrackMatchedDEtaTRDNoCut",
1449 "d#eta of cluster-track vs cluster energy, SM behind TRD, no photon cuts",
1450 nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
4bfeae64 1451 fhTrackMatchedDEtaTRD[0]->SetYTitle("d#eta");
1452 fhTrackMatchedDEtaTRD[0]->SetXTitle("E_{cluster} (GeV)");
b5dbb99b 1453
4bfeae64 1454 fhTrackMatchedDPhiTRD[0] = new TH2F
b5dbb99b 1455 ("hTrackMatchedDPhiTRDNoCut",
1456 "d#phi of cluster-track vs cluster energy, SM behing TRD, no photon cuts",
1457 nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
4bfeae64 1458 fhTrackMatchedDPhiTRD[0]->SetYTitle("d#phi (rad)");
1459 fhTrackMatchedDPhiTRD[0]->SetXTitle("E_{cluster} (GeV)");
b5dbb99b 1460
4bfeae64 1461 fhEOverPTRD[0] = new TH2F ("hEOverPTRDNoCut","matched track E/p vs cluster E, behind TRD, no photon cuts ",
1462 nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
1463 fhEOverPTRD[0]->SetXTitle("E (GeV)");
1464 fhEOverPTRD[0]->SetYTitle("E/p");
1465
1466 outputContainer->Add(fhTrackMatchedDEtaTRD[0]) ;
1467 outputContainer->Add(fhTrackMatchedDPhiTRD[0]) ;
1468 outputContainer->Add(fhEOverPTRD[0]);
b5dbb99b 1469
4bfeae64 1470 fhTrackMatchedDEtaTRD[1] = new TH2F
1471 ("hTrackMatchedDEtaTRD",
1472 "d#eta of cluster-track vs cluster energy, SM behind TRD",
1473 nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1474 fhTrackMatchedDEtaTRD[1]->SetYTitle("d#eta");
1475 fhTrackMatchedDEtaTRD[1]->SetXTitle("E_{cluster} (GeV)");
1476
1477 fhTrackMatchedDPhiTRD[1] = new TH2F
1478 ("hTrackMatchedDPhiTRD",
1479 "d#phi of cluster-track vs cluster energy, SM behing TRD",
1480 nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1481 fhTrackMatchedDPhiTRD[1]->SetYTitle("d#phi (rad)");
1482 fhTrackMatchedDPhiTRD[1]->SetXTitle("E_{cluster} (GeV)");
1483
1484 fhEOverPTRD[1] = new TH2F ("hEOverPTRD","matched track E/p vs cluster E, behind TRD ",
1485 nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
1486 fhEOverPTRD[1]->SetXTitle("E (GeV)");
1487 fhEOverPTRD[1]->SetYTitle("E/p");
1488
1489 outputContainer->Add(fhTrackMatchedDEtaTRD[1]) ;
1490 outputContainer->Add(fhTrackMatchedDPhiTRD[1]) ;
1491 outputContainer->Add(fhEOverPTRD[1]);
b5dbb99b 1492
1493 }
1494
31ae6d59 1495 if(IsDataMC())
1496 {
8d6b7f60 1497
4bfeae64 1498 fhTrackMatchedDEtaMCNoOverlap[0] = new TH2F
1499 ("hTrackMatchedDEtaMCNoOverlapNoCut",
1500 "d#eta of cluster-track vs cluster energy, no other MC particles overlap",
1501 nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1502 fhTrackMatchedDEtaMCNoOverlap[0]->SetYTitle("d#eta");
1503 fhTrackMatchedDEtaMCNoOverlap[0]->SetXTitle("E_{cluster} (GeV)");
1504
1505 fhTrackMatchedDPhiMCNoOverlap[0] = new TH2F
1506 ("hTrackMatchedDPhiMCNoOverlapNoCut",
1507 "d#phi of cluster-track vs cluster energy, no other MC particles overlap",
1508 nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1509 fhTrackMatchedDPhiMCNoOverlap[0]->SetYTitle("d#phi (rad)");
1510 fhTrackMatchedDPhiMCNoOverlap[0]->SetXTitle("E_{cluster} (GeV)");
1511
1512 outputContainer->Add(fhTrackMatchedDEtaMCNoOverlap[0]) ;
1513 outputContainer->Add(fhTrackMatchedDPhiMCNoOverlap[0]) ;
1514
1515 fhTrackMatchedDEtaMCNoOverlap[1] = new TH2F
8d6b7f60 1516 ("hTrackMatchedDEtaMCNoOverlap",
1517 "d#eta of cluster-track vs cluster energy, no other MC particles overlap",
1518 nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
4bfeae64 1519 fhTrackMatchedDEtaMCNoOverlap[1]->SetYTitle("d#eta");
1520 fhTrackMatchedDEtaMCNoOverlap[1]->SetXTitle("E_{cluster} (GeV)");
8d6b7f60 1521
4bfeae64 1522 fhTrackMatchedDPhiMCNoOverlap[1] = new TH2F
8d6b7f60 1523 ("hTrackMatchedDPhiMCNoOverlap",
1524 "d#phi of cluster-track vs cluster energy, no other MC particles overlap",
1525 nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
4bfeae64 1526 fhTrackMatchedDPhiMCNoOverlap[1]->SetYTitle("d#phi (rad)");
1527 fhTrackMatchedDPhiMCNoOverlap[1]->SetXTitle("E_{cluster} (GeV)");
1528
1529 outputContainer->Add(fhTrackMatchedDEtaMCNoOverlap[1]) ;
1530 outputContainer->Add(fhTrackMatchedDPhiMCNoOverlap[1]) ;
1531
1532 fhTrackMatchedDEtaMCOverlap[0] = new TH2F
1533 ("hTrackMatchedDEtaMCOverlapNoCut",
1534 "d#eta of cluster-track vs cluster energy, several MC particles overlap",
1535 nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1536 fhTrackMatchedDEtaMCOverlap[0]->SetYTitle("d#eta");
1537 fhTrackMatchedDEtaMCOverlap[0]->SetXTitle("E_{cluster} (GeV)");
1538
1539 fhTrackMatchedDPhiMCOverlap[0] = new TH2F
1540 ("hTrackMatchedDPhiMCOverlapNoCut",
1541 "d#phi of cluster-track vs cluster energy, several MC particles overlap",
1542 nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1543 fhTrackMatchedDPhiMCOverlap[0]->SetYTitle("d#phi (rad)");
1544 fhTrackMatchedDPhiMCOverlap[0]->SetXTitle("E_{cluster} (GeV)");
8d6b7f60 1545
4bfeae64 1546 outputContainer->Add(fhTrackMatchedDEtaMCOverlap[0]) ;
1547 outputContainer->Add(fhTrackMatchedDPhiMCOverlap[0]) ;
8d6b7f60 1548
4bfeae64 1549 fhTrackMatchedDEtaMCOverlap[1] = new TH2F
8d6b7f60 1550 ("hTrackMatchedDEtaMCOverlap",
1551 "d#eta of cluster-track vs cluster energy, several MC particles overlap",
1552 nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
4bfeae64 1553 fhTrackMatchedDEtaMCOverlap[1]->SetYTitle("d#eta");
1554 fhTrackMatchedDEtaMCOverlap[1]->SetXTitle("E_{cluster} (GeV)");
8d6b7f60 1555
4bfeae64 1556 fhTrackMatchedDPhiMCOverlap[1] = new TH2F
8d6b7f60 1557 ("hTrackMatchedDPhiMCOverlap",
1558 "d#phi of cluster-track vs cluster energy, several MC particles overlap",
1559 nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
4bfeae64 1560 fhTrackMatchedDPhiMCOverlap[1]->SetYTitle("d#phi (rad)");
1561 fhTrackMatchedDPhiMCOverlap[1]->SetXTitle("E_{cluster} (GeV)");
1562
1563 outputContainer->Add(fhTrackMatchedDEtaMCOverlap[1]) ;
1564 outputContainer->Add(fhTrackMatchedDPhiMCOverlap[1]) ;
1565
1566 fhTrackMatchedDEtaMCConversion[0] = new TH2F
1567 ("hTrackMatchedDEtaMCConversionNoCut",
1568 "d#eta of cluster-track vs cluster energy, no other MC particles overlap appart from conversions",
1569 nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1570 fhTrackMatchedDEtaMCConversion[0]->SetYTitle("d#eta");
1571 fhTrackMatchedDEtaMCConversion[0]->SetXTitle("E_{cluster} (GeV)");
8d6b7f60 1572
4bfeae64 1573 fhTrackMatchedDPhiMCConversion[0] = new TH2F
1574 ("hTrackMatchedDPhiMCConversionNoCut",
1575 "d#phi of cluster-track vs cluster energy, no other MC particles overlap appart from conversions",
1576 nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1577 fhTrackMatchedDPhiMCConversion[0]->SetYTitle("d#phi (rad)");
1578 fhTrackMatchedDPhiMCConversion[0]->SetXTitle("E_{cluster} (GeV)");
8d6b7f60 1579
4bfeae64 1580 outputContainer->Add(fhTrackMatchedDEtaMCConversion[0]) ;
1581 outputContainer->Add(fhTrackMatchedDPhiMCConversion[0]) ;
1582
8d6b7f60 1583
4bfeae64 1584 fhTrackMatchedDEtaMCConversion[1] = new TH2F
b5dbb99b 1585 ("hTrackMatchedDEtaMCConversion",
4bfeae64 1586 "d#eta of cluster-track vs cluster energy, no other MC particles overlap appart from conversions",
8d6b7f60 1587 nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
4bfeae64 1588 fhTrackMatchedDEtaMCConversion[1]->SetYTitle("d#eta");
1589 fhTrackMatchedDEtaMCConversion[1]->SetXTitle("E_{cluster} (GeV)");
8d6b7f60 1590
4bfeae64 1591 fhTrackMatchedDPhiMCConversion[1] = new TH2F
b5dbb99b 1592 ("hTrackMatchedDPhiMCConversion",
8d6b7f60 1593 "d#phi of cluster-track vs cluster energy, no other MC particles overlap appart from conversions",
1594 nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
4bfeae64 1595 fhTrackMatchedDPhiMCConversion[1]->SetYTitle("d#phi (rad)");
1596 fhTrackMatchedDPhiMCConversion[1]->SetXTitle("E_{cluster} (GeV)");
8d6b7f60 1597
4bfeae64 1598 outputContainer->Add(fhTrackMatchedDEtaMCConversion[1]) ;
1599 outputContainer->Add(fhTrackMatchedDPhiMCConversion[1]) ;
8d6b7f60 1600
31ae6d59 1601
4bfeae64 1602 fhTrackMatchedMCParticle[0] = new TH2F
31ae6d59 1603 ("hTrackMatchedMCParticleNoCut",
1604 "Origin of particle vs energy",
1605 nptbins,ptmin,ptmax,8,0,8);
4bfeae64 1606 fhTrackMatchedMCParticle[0]->SetXTitle("E (GeV)");
1607 //fhTrackMatchedMCParticle[0]->SetYTitle("Particle type");
1608
1609 fhTrackMatchedMCParticle[0]->GetYaxis()->SetBinLabel(1 ,"Photon");
1610 fhTrackMatchedMCParticle[0]->GetYaxis()->SetBinLabel(2 ,"Electron");
1611 fhTrackMatchedMCParticle[0]->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
1612 fhTrackMatchedMCParticle[0]->GetYaxis()->SetBinLabel(4 ,"Rest");
1613 fhTrackMatchedMCParticle[0]->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
1614 fhTrackMatchedMCParticle[0]->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
1615 fhTrackMatchedMCParticle[0]->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
1616 fhTrackMatchedMCParticle[0]->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
1617
1618 fhTrackMatchedMCParticle[1] = new TH2F
1619 ("hTrackMatchedMCParticle",
1620 "Origin of particle vs energy",
1621 nptbins,ptmin,ptmax,8,0,8);
1622 fhTrackMatchedMCParticle[1]->SetXTitle("E (GeV)");
1623 //fhTrackMatchedMCParticle[1]->SetYTitle("Particle type");
31ae6d59 1624
4bfeae64 1625 fhTrackMatchedMCParticle[1]->GetYaxis()->SetBinLabel(1 ,"Photon");
1626 fhTrackMatchedMCParticle[1]->GetYaxis()->SetBinLabel(2 ,"Electron");
1627 fhTrackMatchedMCParticle[1]->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
1628 fhTrackMatchedMCParticle[1]->GetYaxis()->SetBinLabel(4 ,"Rest");
1629 fhTrackMatchedMCParticle[1]->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
1630 fhTrackMatchedMCParticle[1]->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
1631 fhTrackMatchedMCParticle[1]->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
1632 fhTrackMatchedMCParticle[1]->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
31ae6d59 1633
4bfeae64 1634 outputContainer->Add(fhTrackMatchedMCParticle[0]);
1635 outputContainer->Add(fhTrackMatchedMCParticle[1]);
8d6b7f60 1636
31ae6d59 1637 }
09273901 1638 }
1639
34c16486 1640 if(IsDataMC())
1641 {
f66d95af 1642 TString ptype[] = { "#gamma", "#gamma_{#pi decay}","#gamma_{other decay}", "#pi^{0}","#eta",
1643 "e^{#pm}","#gamma->e^{#pm}","hadron?","Anti-N","Anti-P",
1644 "#gamma_{prompt}","#gamma_{fragmentation}","#gamma_{ISR}","String" } ;
3d5d5078 1645
f66d95af 1646 TString pname[] = { "Photon","PhotonPi0Decay","PhotonOtherDecay","Pi0","Eta","Electron",
1647 "Conversion", "Hadron", "AntiNeutron","AntiProton",
1648 "PhotonPrompt","PhotonFragmentation","PhotonISR","String" } ;
521636d2 1649
34c16486 1650 for(Int_t i = 0; i < fNOriginHistograms; i++)
1651 {
3d5d5078 1652 fhMCE[i] = new TH1F(Form("hE_MC%s",pname[i].Data()),
521636d2 1653 Form("cluster from %s : E ",ptype[i].Data()),
1654 nptbins,ptmin,ptmax);
3d5d5078 1655 fhMCE[i]->SetXTitle("E (GeV)");
1656 outputContainer->Add(fhMCE[i]) ;
521636d2 1657
4c8f7c2e 1658 fhMCPt[i] = new TH1F(Form("hPt_MC%s",pname[i].Data()),
521636d2 1659 Form("cluster from %s : p_{T} ",ptype[i].Data()),
1660 nptbins,ptmin,ptmax);
4c8f7c2e 1661 fhMCPt[i]->SetXTitle("p_{T} (GeV/c)");
1662 outputContainer->Add(fhMCPt[i]) ;
521636d2 1663
4c8f7c2e 1664 fhMCEta[i] = new TH2F(Form("hEta_MC%s",pname[i].Data()),
521636d2 1665 Form("cluster from %s : #eta ",ptype[i].Data()),
1666 nptbins,ptmin,ptmax,netabins,etamin,etamax);
4c8f7c2e 1667 fhMCEta[i]->SetYTitle("#eta");
1668 fhMCEta[i]->SetXTitle("E (GeV)");
1669 outputContainer->Add(fhMCEta[i]) ;
521636d2 1670
4c8f7c2e 1671 fhMCPhi[i] = new TH2F(Form("hPhi_MC%s",pname[i].Data()),
521636d2 1672 Form("cluster from %s : #phi ",ptype[i].Data()),
1673 nptbins,ptmin,ptmax,nphibins,phimin,phimax);
4c8f7c2e 1674 fhMCPhi[i]->SetYTitle("#phi (rad)");
1675 fhMCPhi[i]->SetXTitle("E (GeV)");
1676 outputContainer->Add(fhMCPhi[i]) ;
1677
1678
d9105d92 1679 fhMCDeltaE[i] = new TH2F (Form("hDeltaE_MC%s",pname[i].Data()),
1680 Form("MC - Reco E from %s",pname[i].Data()),
1681 nptbins,ptmin,ptmax, 200,-50,50);
4c8f7c2e 1682 fhMCDeltaE[i]->SetXTitle("#Delta E (GeV)");
1683 outputContainer->Add(fhMCDeltaE[i]);
1684
d9105d92 1685 fhMCDeltaPt[i] = new TH2F (Form("hDeltaPt_MC%s",pname[i].Data()),
1686 Form("MC - Reco p_{T} from %s",pname[i].Data()),
1687 nptbins,ptmin,ptmax, 200,-50,50);
4c8f7c2e 1688 fhMCDeltaPt[i]->SetXTitle("#Delta p_{T} (GeV/c)");
1689 outputContainer->Add(fhMCDeltaPt[i]);
d9105d92 1690
4c8f7c2e 1691 fhMC2E[i] = new TH2F (Form("h2E_MC%s",pname[i].Data()),
1692 Form("E distribution, reconstructed vs generated from %s",pname[i].Data()),
1693 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1694 fhMC2E[i]->SetXTitle("E_{rec} (GeV)");
1695 fhMC2E[i]->SetYTitle("E_{gen} (GeV)");
1696 outputContainer->Add(fhMC2E[i]);
1697
1698 fhMC2Pt[i] = new TH2F (Form("h2Pt_MC%s",pname[i].Data()),
1699 Form("p_T distribution, reconstructed vs generated from %s",pname[i].Data()),
1700 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1701 fhMC2Pt[i]->SetXTitle("p_{T,rec} (GeV/c)");
1702 fhMC2Pt[i]->SetYTitle("p_{T,gen} (GeV/c)");
1703 outputContainer->Add(fhMC2Pt[i]);
1704
521636d2 1705
1706 }
3d5d5078 1707
f66d95af 1708 TString pptype[] = { "#gamma", "#gamma_{#pi decay}","#gamma_{other decay}","hadron?",
1709 "#gamma_{prompt}","#gamma_{fragmentation}","#gamma_{ISR}"} ;
1710
1711 TString ppname[] = { "Photon","PhotonPi0Decay","PhotonOtherDecay","Hadron",
1712 "PhotonPrompt","PhotonFragmentation","PhotonISR"} ;
1713
34c16486 1714 for(Int_t i = 0; i < fNPrimaryHistograms; i++)
1715 {
f66d95af 1716 fhEPrimMC[i] = new TH1F(Form("hEPrim_MC%s",ppname[i].Data()),
1717 Form("primary photon %s : E ",pptype[i].Data()),
3d5d5078 1718 nptbins,ptmin,ptmax);
1719 fhEPrimMC[i]->SetXTitle("E (GeV)");
1720 outputContainer->Add(fhEPrimMC[i]) ;
1721
f66d95af 1722 fhPtPrimMC[i] = new TH1F(Form("hPtPrim_MC%s",ppname[i].Data()),
1723 Form("primary photon %s : p_{T} ",pptype[i].Data()),
3d5d5078 1724 nptbins,ptmin,ptmax);
1725 fhPtPrimMC[i]->SetXTitle("p_{T} (GeV/c)");
1726 outputContainer->Add(fhPtPrimMC[i]) ;
1727
f66d95af 1728 fhYPrimMC[i] = new TH2F(Form("hYPrim_MC%s",ppname[i].Data()),
1729 Form("primary photon %s : Rapidity ",pptype[i].Data()),
3d5d5078 1730 nptbins,ptmin,ptmax,800,-8,8);
1731 fhYPrimMC[i]->SetYTitle("Rapidity");
1732 fhYPrimMC[i]->SetXTitle("E (GeV)");
1733 outputContainer->Add(fhYPrimMC[i]) ;
1734
f66d95af 1735 fhPhiPrimMC[i] = new TH2F(Form("hPhiPrim_MC%s",ppname[i].Data()),
1736 Form("primary photon %s : #phi ",pptype[i].Data()),
3d5d5078 1737 nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1738 fhPhiPrimMC[i]->SetYTitle("#phi (rad)");
1739 fhPhiPrimMC[i]->SetXTitle("E (GeV)");
1740 outputContainer->Add(fhPhiPrimMC[i]) ;
1741
1742
f66d95af 1743 fhEPrimMCAcc[i] = new TH1F(Form("hEPrimAcc_MC%s",ppname[i].Data()),
1744 Form("primary photon %s in acceptance: E ",pptype[i].Data()),
3d5d5078 1745 nptbins,ptmin,ptmax);
1746 fhEPrimMCAcc[i]->SetXTitle("E (GeV)");
1747 outputContainer->Add(fhEPrimMCAcc[i]) ;
1748
f66d95af 1749 fhPtPrimMCAcc[i] = new TH1F(Form("hPtPrimAcc_MC%s",ppname[i].Data()),
1750 Form("primary photon %s in acceptance: p_{T} ",pptype[i].Data()),
3d5d5078 1751 nptbins,ptmin,ptmax);
1752 fhPtPrimMCAcc[i]->SetXTitle("p_{T} (GeV/c)");
1753 outputContainer->Add(fhPtPrimMCAcc[i]) ;
1754
f66d95af 1755 fhYPrimMCAcc[i] = new TH2F(Form("hYPrimAcc_MC%s",ppname[i].Data()),
1756 Form("primary photon %s in acceptance: Rapidity ",pptype[i].Data()),
3d5d5078 1757 nptbins,ptmin,ptmax,100,-1,1);
1758 fhYPrimMCAcc[i]->SetYTitle("Rapidity");
1759 fhYPrimMCAcc[i]->SetXTitle("E (GeV)");
1760 outputContainer->Add(fhYPrimMCAcc[i]) ;
1761
f66d95af 1762 fhPhiPrimMCAcc[i] = new TH2F(Form("hPhiPrimAcc_MC%s",ppname[i].Data()),
1763 Form("primary photon %s in acceptance: #phi ",pptype[i].Data()),
3d5d5078 1764 nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1765 fhPhiPrimMCAcc[i]->SetYTitle("#phi (rad)");
1766 fhPhiPrimMCAcc[i]->SetXTitle("E (GeV)");
1767 outputContainer->Add(fhPhiPrimMCAcc[i]) ;
1768
1769 }
1770
34c16486 1771 if(fFillSSHistograms)
1772 {
3d5d5078 1773 TString ptypess[] = { "#gamma","hadron?","#pi^{0}","#eta","#gamma->e^{#pm}","e^{#pm}"} ;
1774
1775 TString pnamess[] = { "Photon","Hadron","Pi0","Eta","Conversion","Electron"} ;
1776
34c16486 1777 for(Int_t i = 0; i < 6; i++)
1778 {
3d5d5078 1779 fhMCELambda0[i] = new TH2F(Form("hELambda0_MC%s",pnamess[i].Data()),
1780 Form("cluster from %s : E vs #lambda_{0}^{2}",ptypess[i].Data()),
521636d2 1781 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3d5d5078 1782 fhMCELambda0[i]->SetYTitle("#lambda_{0}^{2}");
1783 fhMCELambda0[i]->SetXTitle("E (GeV)");
1784 outputContainer->Add(fhMCELambda0[i]) ;
521636d2 1785
3d5d5078 1786 fhMCELambda1[i] = new TH2F(Form("hELambda1_MC%s",pnamess[i].Data()),
1787 Form("cluster from %s : E vs #lambda_{1}^{2}",ptypess[i].Data()),
521636d2 1788 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3d5d5078 1789 fhMCELambda1[i]->SetYTitle("#lambda_{1}^{2}");
1790 fhMCELambda1[i]->SetXTitle("E (GeV)");
1791 outputContainer->Add(fhMCELambda1[i]) ;
34c16486 1792
3d5d5078 1793 fhMCEDispersion[i] = new TH2F(Form("hEDispersion_MC%s",pnamess[i].Data()),
1794 Form("cluster from %s : E vs dispersion^{2}",ptypess[i].Data()),
521636d2 1795 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3d5d5078 1796 fhMCEDispersion[i]->SetYTitle("D^{2}");
1797 fhMCEDispersion[i]->SetXTitle("E (GeV)");
1798 outputContainer->Add(fhMCEDispersion[i]) ;
34c16486 1799
f66d95af 1800 fhMCNCellsE[i] = new TH2F (Form("hNCellsE_MC%s",pnamess[i].Data()),
34c16486 1801 Form("# of cells in cluster from %s vs E of clusters",ptypess[i].Data()),
f66d95af 1802 nptbins,ptmin,ptmax, nbins,nmin,nmax);
1803 fhMCNCellsE[i]->SetXTitle("E (GeV)");
1804 fhMCNCellsE[i]->SetYTitle("# of cells in cluster");
1805 outputContainer->Add(fhMCNCellsE[i]);
1806
1807 fhMCMaxCellDiffClusterE[i] = new TH2F (Form("hMaxCellDiffClusterE_MC%s",pnamess[i].Data()),
34c16486 1808 Form("energy vs difference of cluster energy from %s - max cell energy / cluster energy, good clusters",ptypess[i].Data()),
1809 nptbins,ptmin,ptmax, 500,0,1.);
f66d95af 1810 fhMCMaxCellDiffClusterE[i]->SetXTitle("E_{cluster} (GeV) ");
1811 fhMCMaxCellDiffClusterE[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
1812 outputContainer->Add(fhMCMaxCellDiffClusterE[i]);
1813
1814 fhMCLambda0vsClusterMaxCellDiffE0[i] = new TH2F(Form("hLambda0vsClusterMaxCellDiffE0_MC%s",pnamess[i].Data()),
34c16486 1815 Form("cluster from %s : #lambda^{2}_{0} vs fraction of energy carried by max cell, E < 2 GeV",ptypess[i].Data()),
1816 ssbins,ssmin,ssmax,500,0,1.);
f66d95af 1817 fhMCLambda0vsClusterMaxCellDiffE0[i]->SetXTitle("#lambda_{0}^{2}");
1818 fhMCLambda0vsClusterMaxCellDiffE0[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
1819 outputContainer->Add(fhMCLambda0vsClusterMaxCellDiffE0[i]) ;
1820
1821 fhMCLambda0vsClusterMaxCellDiffE2[i] = new TH2F(Form("hLambda0vsClusterMaxCellDiffE2_MC%s",pnamess[i].Data()),
1822 Form("cluster from %s : #lambda^{2}_{0} vs fraction of energy carried by max cell, 2< E < 6 GeV",ptypess[i].Data()),
1823 ssbins,ssmin,ssmax,500,0,1.);
1824 fhMCLambda0vsClusterMaxCellDiffE2[i]->SetXTitle("#lambda_{0}^{2}");
1825 fhMCLambda0vsClusterMaxCellDiffE2[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
1826 outputContainer->Add(fhMCLambda0vsClusterMaxCellDiffE2[i]) ;
1827
1828 fhMCLambda0vsClusterMaxCellDiffE6[i] = new TH2F(Form("hLambda0vsClusterMaxCellDiffE6_MC%s",pnamess[i].Data()),
1829 Form("cluster from %s : #lambda^{2}_{0} vs fraction of energy carried by max cell, E > 6 GeV",ptypess[i].Data()),
1830 ssbins,ssmin,ssmax,500,0,1.);
1831 fhMCLambda0vsClusterMaxCellDiffE6[i]->SetXTitle("#lambda_{0}^{2}");
1832 fhMCLambda0vsClusterMaxCellDiffE6[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
1833 outputContainer->Add(fhMCLambda0vsClusterMaxCellDiffE6[i]) ;
1834
1835 fhMCNCellsvsClusterMaxCellDiffE0[i] = new TH2F(Form("hNCellsvsClusterMaxCellDiffE0_MC%s",pnamess[i].Data()),
34c16486 1836 Form("cluster from %s : N cells in cluster vs fraction of energy carried by max cell, E < 2 GeV",ptypess[i].Data()),
1837 nbins/5,nmin,nmax/5,500,0,1.);
f66d95af 1838 fhMCNCellsvsClusterMaxCellDiffE0[i]->SetXTitle("N cells in cluster");
1839 fhMCNCellsvsClusterMaxCellDiffE0[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
1840 outputContainer->Add(fhMCNCellsvsClusterMaxCellDiffE0[i]) ;
1841
1842 fhMCNCellsvsClusterMaxCellDiffE2[i] = new TH2F(Form("hNCellsvsClusterMaxCellDiffE2_MC%s",pnamess[i].Data()),
34c16486 1843 Form("cluster from %s : N cells in cluster vs fraction of energy carried by max cell, 2< E < 6 GeV",ptypess[i].Data()),
1844 nbins/5,nmin,nmax/5,500,0,1.);
f66d95af 1845 fhMCNCellsvsClusterMaxCellDiffE2[i]->SetXTitle("N cells in cluster");
1846 fhMCNCellsvsClusterMaxCellDiffE2[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
1847 outputContainer->Add(fhMCNCellsvsClusterMaxCellDiffE2[i]) ;
1848
1849 fhMCNCellsvsClusterMaxCellDiffE6[i] = new TH2F(Form("hNCellsvsClusterMaxCellDiffE6_MC%s",pnamess[i].Data()),
34c16486 1850 Form("cluster from %s : N cells in cluster vs fraction of energy carried by max cell, E > 6 GeV",ptypess[i].Data()),
1851 nbins/5,nmin,nmax/5,500,0,1.);
f66d95af 1852 fhMCNCellsvsClusterMaxCellDiffE6[i]->SetXTitle("N cells in cluster");
1853 fhMCNCellsvsClusterMaxCellDiffE6[i]->SetYTitle("E (GeV)");
1854 outputContainer->Add(fhMCNCellsvsClusterMaxCellDiffE6[i]) ;
1855
34c16486 1856 if(fCalorimeter=="EMCAL")
1857 {
1858 fhMCEDispEta[i] = new TH2F (Form("hEDispEtaE_MC%s",pnamess[i].Data()),
1859 Form("cluster from %s : #sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs E",ptypess[i].Data()),
1860 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1861 fhMCEDispEta[i]->SetXTitle("E (GeV)");
1862 fhMCEDispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
1863 outputContainer->Add(fhMCEDispEta[i]);
1864
1865 fhMCEDispPhi[i] = new TH2F (Form("hEDispPhiE_MC%s",pnamess[i].Data()),
1866 Form("cluster from %s : #sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs E",ptypess[i].Data()),
1867 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1868 fhMCEDispPhi[i]->SetXTitle("E (GeV)");
1869 fhMCEDispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1870 outputContainer->Add(fhMCEDispPhi[i]);
1871
1872 fhMCESumEtaPhi[i] = new TH2F (Form("hESumEtaPhiE_MC%s",pnamess[i].Data()),
1873 Form("cluster from %s : #sigma'^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs E",ptypess[i].Data()),
1874 nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
1875 fhMCESumEtaPhi[i]->SetXTitle("E (GeV)");
1876 fhMCESumEtaPhi[i]->SetYTitle("#sigma'^{2}_{#eta #phi}");
1877 outputContainer->Add(fhMCESumEtaPhi[i]);
1878
1879 fhMCEDispEtaPhiDiff[i] = new TH2F (Form("hEDispEtaPhiDiffE_MC%s",pnamess[i].Data()),
1880 Form("cluster from %s : #sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs E",ptypess[i].Data()),
1881 nptbins,ptmin,ptmax,200,-10,10);
1882 fhMCEDispEtaPhiDiff[i]->SetXTitle("E (GeV)");
1883 fhMCEDispEtaPhiDiff[i]->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
1884 outputContainer->Add(fhMCEDispEtaPhiDiff[i]);
1885
1886 fhMCESphericity[i] = new TH2F (Form("hESphericity_MC%s",pnamess[i].Data()),
1887 Form("cluster from %s : (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi}) vs E",ptypess[i].Data()),
1888 nptbins,ptmin,ptmax, 200,-1,1);
1889 fhMCESphericity[i]->SetXTitle("E (GeV)");
1890 fhMCESphericity[i]->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
1891 outputContainer->Add(fhMCESphericity[i]);
1892
1893 Int_t bin[] = {0,2,4,6,10,1000};
1894 for(Int_t ie = 0; ie < 5; ie++)
1895 {
bfdcf7fb 1896 fhMCDispEtaDispPhi[ie][i] = new TH2F (Form("hMCDispEtaDispPhi_EBin%d_MC%s",ie,pnamess[i].Data()),
34c16486 1897 Form("cluster from %s : #sigma^{2}_{#phi #phi} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",pnamess[i].Data(),bin[ie],bin[ie+1]),
1898 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
bfdcf7fb 1899 fhMCDispEtaDispPhi[ie][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
1900 fhMCDispEtaDispPhi[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1901 outputContainer->Add(fhMCDispEtaDispPhi[ie][i]);
1902
1903 fhMCLambda0DispEta[ie][i] = new TH2F (Form("hMCLambda0DispEta_EBin%d_MC%s",ie,pname[i].Data()),
1904 Form("cluster from %s : #lambda^{2}_{0} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
1905 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1906 fhMCLambda0DispEta[ie][i]->SetXTitle("#lambda^{2}_{0}");
1907 fhMCLambda0DispEta[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1908 outputContainer->Add(fhMCLambda0DispEta[ie][i]);
1909
1910 fhMCLambda0DispPhi[ie][i] = new TH2F (Form("hMCLambda0DispPhi_EBin%d_MC%s",ie,pname[i].Data()),
1911 Form("cluster from %s :#lambda^{2}_{0} vs #sigma^{2}_{#phi #phi} for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
1912 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1913 fhMCLambda0DispPhi[ie][i]->SetXTitle("#lambda^{2}_{0}");
1914 fhMCLambda0DispPhi[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1915 outputContainer->Add(fhMCLambda0DispPhi[ie][i]);
34c16486 1916 }
1917
1918 }
1919 }// loop
3d5d5078 1920
1921 if(!GetReader()->IsEmbeddedClusterSelectionOn())
1922 {
1923 fhMCPhotonELambda0NoOverlap = new TH2F("hELambda0_MCPhoton_NoOverlap",
1924 "cluster from Photon : E vs #lambda_{0}^{2}",
1925 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1926 fhMCPhotonELambda0NoOverlap->SetYTitle("#lambda_{0}^{2}");
1927 fhMCPhotonELambda0NoOverlap->SetXTitle("E (GeV)");
1928 outputContainer->Add(fhMCPhotonELambda0NoOverlap) ;
1929
3d5d5078 1930 fhMCPhotonELambda0TwoOverlap = new TH2F("hELambda0_MCPhoton_TwoOverlap",
1931 "cluster from Photon : E vs #lambda_{0}^{2}",
1932 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1933 fhMCPhotonELambda0TwoOverlap->SetYTitle("#lambda_{0}^{2}");
1934 fhMCPhotonELambda0TwoOverlap->SetXTitle("E (GeV)");
1935 outputContainer->Add(fhMCPhotonELambda0TwoOverlap) ;
1936
3d5d5078 1937 fhMCPhotonELambda0NOverlap = new TH2F("hELambda0_MCPhoton_NOverlap",
1938 "cluster from Photon : E vs #lambda_{0}^{2}",
1939 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1940 fhMCPhotonELambda0NOverlap->SetYTitle("#lambda_{0}^{2}");
1941 fhMCPhotonELambda0NOverlap->SetXTitle("E (GeV)");
1942 outputContainer->Add(fhMCPhotonELambda0NOverlap) ;
521636d2 1943
3d5d5078 1944 } //No embedding
1945
1946 //Fill histograms to check shape of embedded clusters
1947 if(GetReader()->IsEmbeddedClusterSelectionOn())
1948 {
1949
1950 fhEmbeddedSignalFractionEnergy = new TH2F("hEmbeddedSignal_FractionEnergy",
34c16486 1951 "Energy Fraction of embedded signal versus cluster energy",
1952 nptbins,ptmin,ptmax,100,0.,1.);
3d5d5078 1953 fhEmbeddedSignalFractionEnergy->SetYTitle("Fraction");
1954 fhEmbeddedSignalFractionEnergy->SetXTitle("E (GeV)");
1955 outputContainer->Add(fhEmbeddedSignalFractionEnergy) ;
1956
1957 fhEmbedPhotonELambda0FullSignal = new TH2F("hELambda0_EmbedPhoton_FullSignal",
34c16486 1958 "cluster from Photon embedded with more than 90% energy in cluster : E vs #lambda_{0}^{2}",
1959 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3d5d5078 1960 fhEmbedPhotonELambda0FullSignal->SetYTitle("#lambda_{0}^{2}");
1961 fhEmbedPhotonELambda0FullSignal->SetXTitle("E (GeV)");
1962 outputContainer->Add(fhEmbedPhotonELambda0FullSignal) ;
34c16486 1963
3d5d5078 1964 fhEmbedPhotonELambda0MostlySignal = new TH2F("hELambda0_EmbedPhoton_MostlySignal",
34c16486 1965 "cluster from Photon embedded with 50% to 90% energy in cluster : E vs #lambda_{0}^{2}",
1966 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3d5d5078 1967 fhEmbedPhotonELambda0MostlySignal->SetYTitle("#lambda_{0}^{2}");
1968 fhEmbedPhotonELambda0MostlySignal->SetXTitle("E (GeV)");
1969 outputContainer->Add(fhEmbedPhotonELambda0MostlySignal) ;
1970
3d5d5078 1971 fhEmbedPhotonELambda0MostlyBkg = new TH2F("hELambda0_EmbedPhoton_MostlyBkg",
34c16486 1972 "cluster from Photon embedded with 10% to 50% energy in cluster : E vs #lambda_{0}^{2}",
1973 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3d5d5078 1974 fhEmbedPhotonELambda0MostlyBkg->SetYTitle("#lambda_{0}^{2}");
1975 fhEmbedPhotonELambda0MostlyBkg->SetXTitle("E (GeV)");
1976 outputContainer->Add(fhEmbedPhotonELambda0MostlyBkg) ;
34c16486 1977
3d5d5078 1978 fhEmbedPhotonELambda0FullBkg = new TH2F("hELambda0_EmbedPhoton_FullBkg",
34c16486 1979 "cluster from Photonm embedded with 0% to 10% energy in cluster : E vs #lambda_{0}^{2}",
1980 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3d5d5078 1981 fhEmbedPhotonELambda0FullBkg->SetYTitle("#lambda_{0}^{2}");
1982 fhEmbedPhotonELambda0FullBkg->SetXTitle("E (GeV)");
1983 outputContainer->Add(fhEmbedPhotonELambda0FullBkg) ;
1984
3d5d5078 1985 fhEmbedPi0ELambda0FullSignal = new TH2F("hELambda0_EmbedPi0_FullSignal",
34c16486 1986 "cluster from Pi0 embedded with more than 90% energy in cluster : E vs #lambda_{0}^{2}",
1987 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3d5d5078 1988 fhEmbedPi0ELambda0FullSignal->SetYTitle("#lambda_{0}^{2}");
1989 fhEmbedPi0ELambda0FullSignal->SetXTitle("E (GeV)");
1990 outputContainer->Add(fhEmbedPi0ELambda0FullSignal) ;
34c16486 1991
3d5d5078 1992 fhEmbedPi0ELambda0MostlySignal = new TH2F("hELambda0_EmbedPi0_MostlySignal",
34c16486 1993 "cluster from Pi0 embedded with 50% to 90% energy in cluster : E vs #lambda_{0}^{2}",
1994 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3d5d5078 1995 fhEmbedPi0ELambda0MostlySignal->SetYTitle("#lambda_{0}^{2}");
1996 fhEmbedPi0ELambda0MostlySignal->SetXTitle("E (GeV)");
1997 outputContainer->Add(fhEmbedPi0ELambda0MostlySignal) ;
1998
3d5d5078 1999 fhEmbedPi0ELambda0MostlyBkg = new TH2F("hELambda0_EmbedPi0_MostlyBkg",
34c16486 2000 "cluster from Pi0 embedded with 10% to 50% energy in cluster : E vs #lambda_{0}^{2}",
2001 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3d5d5078 2002 fhEmbedPi0ELambda0MostlyBkg->SetYTitle("#lambda_{0}^{2}");
2003 fhEmbedPi0ELambda0MostlyBkg->SetXTitle("E (GeV)");
2004 outputContainer->Add(fhEmbedPi0ELambda0MostlyBkg) ;
2005
3d5d5078 2006 fhEmbedPi0ELambda0FullBkg = new TH2F("hELambda0_EmbedPi0_FullBkg",
34c16486 2007 "cluster from Pi0 embedded with 0% to 10% energy in cluster : E vs #lambda_{0}^{2}",
2008 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3d5d5078 2009 fhEmbedPi0ELambda0FullBkg->SetYTitle("#lambda_{0}^{2}");
2010 fhEmbedPi0ELambda0FullBkg->SetXTitle("E (GeV)");
2011 outputContainer->Add(fhEmbedPi0ELambda0FullBkg) ;
34c16486 2012
3d5d5078 2013 }// embedded histograms
2014
521636d2 2015
2016 }// Fill SS MC histograms
2017
477d6cee 2018 }//Histos with MC
09273901 2019
477d6cee 2020 return outputContainer ;
2021
1c5acb87 2022}
2023
34c16486 2024//_______________________
6639984f 2025void AliAnaPhoton::Init()
2026{
2027
2028 //Init
2029 //Do some checks
34c16486 2030 if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD())
2031 {
591cc579 2032 printf("AliAnaPhoton::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
6639984f 2033 abort();
2034 }
34c16486 2035 else if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD())
2036 {
591cc579 2037 printf("AliAnaPhoton::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
6639984f 2038 abort();
2039 }
2040
49b5c49b 2041 if(GetReader()->GetDataType() == AliCaloTrackReader::kMC) GetCaloPID()->SwitchOnBayesian();
2042
6639984f 2043}
2044
1c5acb87 2045//____________________________________________________________________________
2046void AliAnaPhoton::InitParameters()
2047{
2048
2049 //Initialize the parameters of the analysis.
a3aebfff 2050 AddToHistogramsName("AnaPhoton_");
521636d2 2051
6175da48 2052 fCalorimeter = "EMCAL" ;
2053 fMinDist = 2.;
2054 fMinDist2 = 4.;
2055 fMinDist3 = 5.;
1e86c71e 2056
caa8a222 2057 fTimeCutMin =-1000000;
2058 fTimeCutMax = 1000000;
6175da48 2059 fNCellsCut = 0;
2ac125bf 2060
1e86c71e 2061 fRejectTrackMatch = kTRUE ;
1e86c71e 2062
1c5acb87 2063}
2064
2065//__________________________________________________________________
2066void AliAnaPhoton::MakeAnalysisFillAOD()
2067{
f8006433 2068 //Do photon analysis and fill aods
f37fa8d2 2069
6175da48 2070 //Get the vertex
5025c139 2071 Double_t v[3] = {0,0,0}; //vertex ;
2072 GetReader()->GetVertex(v);
f8006433 2073
f37fa8d2 2074 //Select the Calorimeter of the photon
c8fe2783 2075 TObjArray * pl = 0x0;
71e3889f 2076 AliVCaloCells* cells = 0;
2077 if (fCalorimeter == "PHOS" )
2078 {
2079 pl = GetPHOSClusters();
2080 cells = GetPHOSCells();
2081 }
477d6cee 2082 else if (fCalorimeter == "EMCAL")
71e3889f 2083 {
2084 pl = GetEMCALClusters();
2085 cells = GetEMCALCells();
2086 }
5ae09196 2087
34c16486 2088 if(!pl)
2089 {
5ae09196 2090 Info("MakeAnalysisFillAOD","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
2091 return;
2092 }
521636d2 2093
fc195fd0 2094 // Loop on raw clusters before filtering in the reader and fill control histogram
34c16486 2095 if((GetReader()->GetEMCALClusterListName()=="" && fCalorimeter=="EMCAL") || fCalorimeter=="PHOS")
2096 {
2097 for(Int_t iclus = 0; iclus < GetReader()->GetInputEvent()->GetNumberOfCaloClusters(); iclus++ )
2098 {
fc195fd0 2099 AliVCluster * clus = GetReader()->GetInputEvent()->GetCaloCluster(iclus);
2100 if (fCalorimeter == "PHOS" && clus->IsPHOS() && clus->E() > GetReader()->GetPHOSPtMin() ) fhClusterCuts[0]->Fill(clus->E());
2101 else if(fCalorimeter == "EMCAL" && clus->IsEMCAL() && clus->E() > GetReader()->GetEMCALPtMin()) fhClusterCuts[0]->Fill(clus->E());
2102 }
2103 }
34c16486 2104 else
2105 { // reclusterized
fc195fd0 2106 TClonesArray * clusterList = 0;
2107 if(GetReader()->GetOutputEvent())
4a9e1073 2108 clusterList = dynamic_cast<TClonesArray*> (GetReader()->GetOutputEvent()->FindListObject(GetReader()->GetEMCALClusterListName()));
34c16486 2109 if(clusterList)
2110 {
fc195fd0 2111 Int_t nclusters = clusterList->GetEntriesFast();
34c16486 2112 for (Int_t iclus = 0; iclus < nclusters; iclus++)
2113 {
fc195fd0 2114 AliVCluster * clus = dynamic_cast<AliVCluster*> (clusterList->At(iclus));
2115 if(clus)fhClusterCuts[0]->Fill(clus->E());
4a9e1073 2116 }
fc195fd0 2117 }
2118 }
fc195fd0 2119
6175da48 2120 //Init arrays, variables, get number of clusters
1e86c71e 2121 TLorentzVector mom, mom2 ;
2122 Int_t nCaloClusters = pl->GetEntriesFast();
20218aea 2123
6175da48 2124 if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillAOD() - input %s cluster entries %d\n", fCalorimeter.Data(), nCaloClusters);
521636d2 2125
6175da48 2126 //----------------------------------------------------
2127 // Fill AOD with PHOS/EMCAL AliAODPWG4Particle objects
2128 //----------------------------------------------------
2129 // Loop on clusters
34c16486 2130 for(Int_t icalo = 0; icalo < nCaloClusters; icalo++)
2131 {
0ae57829 2132 AliVCluster * calo = (AliVCluster*) (pl->At(icalo));
2133 //printf("calo %d, %f\n",icalo,calo->E());
521636d2 2134
f8006433 2135 //Get the index where the cluster comes, to retrieve the corresponding vertex
c8fe2783 2136 Int_t evtIndex = 0 ;
34c16486 2137 if (GetMixedEvent())
2138 {
c8fe2783 2139 evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
5025c139 2140 //Get the vertex and check it is not too large in z
96539743 2141 if(TMath::Abs(GetVertex(evtIndex)[2])> GetZvertexCut()) continue;
c8fe2783 2142 }
521636d2 2143
2144 //Cluster selection, not charged, with photon id and in fiducial cut
34c16486 2145 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
2146 {
f8006433 2147 calo->GetMomentum(mom,GetVertex(evtIndex)) ;}//Assume that come from vertex in straight line
34c16486 2148 else
2149 {
f8006433 2150 Double_t vertex[]={0,0,0};
2151 calo->GetMomentum(mom,vertex) ;
2152 }
c8fe2783 2153
6175da48 2154 //--------------------------------------
2155 // Cluster selection
2156 //--------------------------------------
c4a7d28a 2157 if(!ClusterSelected(calo,mom)) continue;
6175da48 2158
2159 //----------------------------
2160 //Create AOD for analysis
2161 //----------------------------
2162 AliAODPWG4Particle aodph = AliAODPWG4Particle(mom);
2163
2164 //...............................................
2165 //Set the indeces of the original caloclusters (MC, ID), and calorimeter
2166 Int_t label = calo->GetLabel();
2167 aodph.SetLabel(label);
6175da48 2168 aodph.SetCaloLabel(calo->GetID(),-1);
2169 aodph.SetDetector(fCalorimeter);
c4a7d28a 2170 //printf("Index %d, Id %d, iaod %d\n",icalo, calo->GetID(),GetOutputAODBranch()->GetEntriesFast());
521636d2 2171
6175da48 2172 //...............................................
2173 //Set bad channel distance bit
c4a7d28a 2174 Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
f37fa8d2 2175 if (distBad > fMinDist3) aodph.SetDistToBad(2) ;
477d6cee 2176 else if(distBad > fMinDist2) aodph.SetDistToBad(1) ;
f37fa8d2 2177 else aodph.SetDistToBad(0) ;
af7b3903 2178 //printf("DistBad %f Bit %d\n",distBad, aodph.DistToBad());
c8fe2783 2179
8d6b7f60 2180 //--------------------------------------------------------------------------------------
2181 // Play with the MC stack if available
2182 //--------------------------------------------------------------------------------------
2183
2184 //Check origin of the candidates
2185 Int_t tag = -1;
2186
34c16486 2187 if(IsDataMC())
2188 {
8d6b7f60 2189 tag = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader(), aodph.GetInputFileIndex());
2190 aodph.SetTag(tag);
2191
2192 if(GetDebug() > 0)
2193 printf("AliAnaPhoton::MakeAnalysisFillAOD() - Origin of candidate, bit map %d\n",aodph.GetTag());
2194 }//Work with stack also
2195
2196
521636d2 2197 //--------------------------------------------------------------------------------------
2198 //Fill some shower shape histograms before PID is applied
2199 //--------------------------------------------------------------------------------------
2200
8d6b7f60 2201 FillShowerShapeHistograms(calo,tag);
6175da48 2202
2203 //-------------------------------------
f37fa8d2 2204 //PID selection or bit setting
6175da48 2205 //-------------------------------------
49b5c49b 2206
6175da48 2207 //...............................................
2208 // Data, PID check on
3c1d9afb 2209 if(IsCaloPIDOn())
2210 {
49b5c49b 2211 // Get most probable PID, 2 options check bayesian PID weights or redo PID
2212 // By default, redo PID
09273901 2213
3c1d9afb 2214 aodph.SetIdentifiedParticleType(GetCaloPID()->GetIdentifiedParticleType(calo));
477d6cee 2215
21a4b1c0 2216 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetIdentifiedParticleType());
477d6cee 2217
f37fa8d2 2218 //If cluster does not pass pid, not photon, skip it.
21a4b1c0 2219 if(aodph.GetIdentifiedParticleType() != AliCaloPID::kPhoton) continue ;
477d6cee 2220
2221 }
8d6b7f60 2222
6175da48 2223 //...............................................
2224 // Data, PID check off
3c1d9afb 2225 else
2226 {
f37fa8d2 2227 //Set PID bits for later selection (AliAnaPi0 for example)
49b5c49b 2228 //GetIdentifiedParticleType already called in SetPIDBits.
2229
3c1d9afb 2230 GetCaloPID()->SetPIDBits(calo,&aodph, GetCaloUtils(),GetReader()->GetInputEvent());
49b5c49b 2231
a3aebfff 2232 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PID Bits set \n");
477d6cee 2233 }
2234
3c1d9afb 2235 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Photon selection cuts passed: pT %3.2f, pdg %d\n",
2236 aodph.Pt(), aodph.GetIdentifiedParticleType());
09273901 2237
fc195fd0 2238 fhClusterCuts[8]->Fill(calo->E());
09273901 2239
2240 // Matching after cuts
4bfeae64 2241 if(fFillTMHisto) FillTrackMatchingResidualHistograms(calo,1);
09273901 2242
5c46c992 2243 // Add number of local maxima to AOD, method name in AOD to be FIXED
5c46c992 2244
71e3889f 2245 aodph.SetFiducialArea(GetCaloUtils()->GetNumberOfLocalMaxima(calo, cells));
5c46c992 2246
5c46c992 2247
f37fa8d2 2248 //Add AOD with photon object to aod branch
477d6cee 2249 AddAODParticle(aodph);
2250
2251 }//loop
5812a064 2252
f37fa8d2 2253 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() End fill AODs, with %d entries \n",GetOutputAODBranch()->GetEntriesFast());
477d6cee 2254
1c5acb87 2255}
2256
2257//__________________________________________________________________
2258void AliAnaPhoton::MakeAnalysisFillHistograms()
2259{
6175da48 2260 //Fill histograms
f8006433 2261
6175da48 2262 //-------------------------------------------------------------------
577d9801 2263 // Access MC information in stack if requested, check that it exists.
521636d2 2264 AliStack * stack = 0x0;
2265 TParticle * primary = 0x0;
2266 TClonesArray * mcparticles = 0x0;
2267 AliAODMCParticle * aodprimary = 0x0;
2268
34c16486 2269 if(IsDataMC())
2270 {
2271 if(GetReader()->ReadStack())
2272 {
577d9801 2273 stack = GetMCStack() ;
34c16486 2274 if(!stack)
2275 {
3d5d5078 2276 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called? STOP\n");
2277 abort();
577d9801 2278 }
f8006433 2279
577d9801 2280 }
34c16486 2281 else if(GetReader()->ReadAODMCParticles())
2282 {
f8006433 2283
577d9801 2284 //Get the list of MC particles
521636d2 2285 mcparticles = GetReader()->GetAODMCParticles(0);
34c16486 2286 if(!mcparticles && GetDebug() > 0)
2287 {
3d5d5078 2288 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Standard MCParticles not available!\n");
577d9801 2289 }
577d9801 2290 }
2291 }// is data and MC
521636d2 2292
6175da48 2293
2294 // Get vertex
2244659d 2295 Double_t v[3] = {0,0,0}; //vertex ;
2296 GetReader()->GetVertex(v);
6175da48 2297 //fhVertex->Fill(v[0],v[1],v[2]);
2298 if(TMath::Abs(v[2]) > GetZvertexCut()) return ; // done elsewhere for Single Event analysis, but there for mixed event
2299
2300 //----------------------------------
577d9801 2301 //Loop on stored AOD photons
2302 Int_t naod = GetOutputAODBranch()->GetEntriesFast();
577d9801 2303 if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
521636d2 2304
3c1d9afb 2305 for(Int_t iaod = 0; iaod < naod ; iaod++)
2306 {
577d9801 2307 AliAODPWG4Particle* ph = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
2308 Int_t pdg = ph->GetIdentifiedParticleType();
521636d2 2309
577d9801 2310 if(GetDebug() > 3)
3c1d9afb 2311 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - PDG %d, MC TAG %d, Calorimeter %s\n",
2312 ph->GetIdentifiedParticleType(),ph->GetTag(), (ph->GetDetector()).Data()) ;
521636d2 2313
577d9801 2314 //If PID used, fill histos with photons in Calorimeter fCalorimeter
2315 if(IsCaloPIDOn() && pdg != AliCaloPID::kPhoton) continue;
2316 if(ph->GetDetector() != fCalorimeter) continue;
521636d2 2317
577d9801 2318 if(GetDebug() > 2)
2319 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - ID Photon: pt %f, phi %f, eta %f\n", ph->Pt(),ph->Phi(),ph->Eta()) ;
521636d2 2320
6175da48 2321 //................................
577d9801 2322 //Fill photon histograms
2323 Float_t ptcluster = ph->Pt();
2324 Float_t phicluster = ph->Phi();
2325 Float_t etacluster = ph->Eta();
2326 Float_t ecluster = ph->E();
521636d2 2327
20218aea 2328 fhEPhoton ->Fill(ecluster);
577d9801 2329 fhPtPhoton ->Fill(ptcluster);
2330 fhPhiPhoton ->Fill(ptcluster,phicluster);
2331 fhEtaPhoton ->Fill(ptcluster,etacluster);
521636d2 2332 if (ecluster > 0.5) fhEtaPhiPhoton ->Fill(etacluster, phicluster);
20218aea 2333 else if(GetMinPt() < 0.5) fhEtaPhi05Photon->Fill(etacluster, phicluster);
2334
5812a064 2335
2336 //Get original cluster, to recover some information
2337 Int_t absID = 0;
2338 Float_t maxCellFraction = 0;
2339 AliVCaloCells* cells = 0;
2340 TObjArray * clusters = 0;
34c16486 2341 if(fCalorimeter == "EMCAL")
2342 {
5812a064 2343 cells = GetEMCALCells();
2344 clusters = GetEMCALClusters();
2345 }
34c16486 2346 else
2347 {
5812a064 2348 cells = GetPHOSCells();
2349 clusters = GetPHOSClusters();
6175da48 2350 }
20218aea 2351
5812a064 2352 Int_t iclus = -1;
2353 AliVCluster *cluster = FindCluster(clusters,ph->GetCaloLabel(0),iclus);
5c46c992 2354 if(cluster)
2355 {
06f1b12a 2356 absID = GetCaloUtils()->GetMaxEnergyCell(cells, cluster,maxCellFraction);
2357
2358 // Control histograms
2359 fhMaxCellDiffClusterE->Fill(ph->E(),maxCellFraction);
2360 fhNCellsE ->Fill(ph->E(),cluster->GetNCells());
2361 fhTimeE ->Fill(ph->E(),cluster->GetTOF()*1.e9);
5c46c992 2362 if(cells)
2363 {
2364 for(Int_t icell = 0; icell < cluster->GetNCells(); icell++)
2365 fhCellsE->Fill(ph->E(),cells->GetCellAmplitude(cluster->GetCellsAbsId()[icell]));
2366 }
06f1b12a 2367 }
5812a064 2368
6175da48 2369 //.......................................
577d9801 2370 //Play with the MC data if available
34c16486 2371 if(IsDataMC())
2372 {
3d5d5078 2373 FillAcceptanceHistograms();
2374
4c8f7c2e 2375 //....................................................................
2376 // Access MC information in stack if requested, check that it exists.
2377 Int_t label =ph->GetLabel();
34c16486 2378 if(label < 0)
2379 {
4c8f7c2e 2380 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** bad label ***: label %d \n", label);
2381 continue;
2382 }
2383
2384 Float_t eprim = 0;
2385 Float_t ptprim = 0;
34c16486 2386 if(GetReader()->ReadStack())
2387 {
2388 if(label >= stack->GetNtrack())
2389 {
4c8f7c2e 2390 if(GetDebug() > 2) printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n", label, stack->GetNtrack());
2391 continue ;
2392 }
2393
2394 primary = stack->Particle(label);
34c16486 2395 if(!primary)
2396 {
4c8f7c2e 2397 printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** no primary ***: label %d \n", label);
2398 continue;
2399 }
34c16486 2400
4c8f7c2e 2401 eprim = primary->Energy();
2402 ptprim = primary->Pt();
2403
2404 }
34c16486 2405 else if(GetReader()->ReadAODMCParticles())
2406 {
4c8f7c2e 2407 //Check which is the input
34c16486 2408 if(ph->GetInputFileIndex() == 0)
2409 {
4c8f7c2e 2410 if(!mcparticles) continue;
34c16486 2411 if(label >= mcparticles->GetEntriesFast())
2412 {
4c8f7c2e 2413 if(GetDebug() > 2) printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n",
2414 label, mcparticles->GetEntriesFast());
2415 continue ;
2416 }
2417 //Get the particle
2418 aodprimary = (AliAODMCParticle*) mcparticles->At(label);
2419
2420 }
2421
34c16486 2422 if(!aodprimary)
2423 {
4c8f7c2e 2424 printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** no primary ***: label %d \n", label);
2425 continue;
2426 }
2427
2428 eprim = aodprimary->E();
2429 ptprim = aodprimary->Pt();
2430
2431 }
2432
577d9801 2433 Int_t tag =ph->GetTag();
521636d2 2434
c5693f62 2435 if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) && fhMCE[kmcPhoton])
3d5d5078 2436 {
c5693f62 2437 fhMCE [kmcPhoton] ->Fill(ecluster);
2438 fhMCPt [kmcPhoton] ->Fill(ptcluster);
2439 fhMCPhi[kmcPhoton] ->Fill(ecluster,phicluster);
2440 fhMCEta[kmcPhoton] ->Fill(ecluster,etacluster);
4c8f7c2e 2441
c5693f62 2442 fhMC2E[kmcPhoton] ->Fill(ecluster, eprim);
2443 fhMC2Pt[kmcPhoton] ->Fill(ptcluster, ptprim);
2444 fhMCDeltaE[kmcPhoton] ->Fill(ecluster,eprim-ecluster);
2445 fhMCDeltaPt[kmcPhoton]->Fill(ptcluster,ptprim-ptcluster);
3d5d5078 2446
c5693f62 2447 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) && fhMCE[kmcConversion])
3d5d5078 2448 {
c5693f62 2449 fhMCE [kmcConversion] ->Fill(ecluster);
2450 fhMCPt [kmcConversion] ->Fill(ptcluster);
2451 fhMCPhi[kmcConversion] ->Fill(ecluster,phicluster);
2452 fhMCEta[kmcConversion] ->Fill(ecluster,etacluster);
4c8f7c2e 2453
c5693f62 2454 fhMC2E[kmcConversion] ->Fill(ecluster, eprim);
2455 fhMC2Pt[kmcConversion] ->Fill(ptcluster, ptprim);
2456 fhMCDeltaE[kmcConversion] ->Fill(ecluster,eprim-ecluster);
2457 fhMCDeltaPt[kmcConversion]->Fill(ptcluster,ptprim-ptcluster);
3d5d5078 2458 }
2459
34c16486 2460 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt) && fhMCE[kmcPrompt])
2461 {
c5693f62 2462 fhMCE [kmcPrompt] ->Fill(ecluster);
2463 fhMCPt [kmcPrompt] ->Fill(ptcluster);
2464 fhMCPhi[kmcPrompt] ->Fill(ecluster,phicluster);
2465 fhMCEta[kmcPrompt] ->Fill(ecluster,etacluster);
4c8f7c2e 2466
c5693f62 2467 fhMC2E[kmcPrompt] ->Fill(ecluster, eprim);
2468 fhMC2Pt[kmcPrompt] ->Fill(ptcluster, ptprim);
2469 fhMCDeltaE[kmcPrompt] ->Fill(ecluster,eprim-ecluster);
2470 fhMCDeltaPt[kmcPrompt]->Fill(ptcluster,ptprim-ptcluster);
4c8f7c2e 2471
3d5d5078 2472 }
c5693f62 2473 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation)&& fhMCE[kmcFragmentation])
3d5d5078 2474 {
c5693f62 2475 fhMCE [kmcFragmentation] ->Fill(ecluster);
2476 fhMCPt [kmcFragmentation] ->Fill(ptcluster);
2477 fhMCPhi[kmcFragmentation] ->Fill(ecluster,phicluster);
2478 fhMCEta[kmcFragmentation] ->Fill(ecluster,etacluster);
4c8f7c2e 2479
c5693f62 2480 fhMC2E[kmcFragmentation] ->Fill(ecluster, eprim);
2481 fhMC2Pt[kmcFragmentation] ->Fill(ptcluster, ptprim);
2482 fhMCDeltaE[kmcFragmentation] ->Fill(ecluster,eprim-ecluster);
2483 fhMCDeltaPt[kmcFragmentation]->Fill(ptcluster,ptprim-ptcluster);
3d5d5078 2484
2485 }
c5693f62 2486 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR)&& fhMCE[kmcISR])
3d5d5078 2487 {
c5693f62 2488 fhMCE [kmcISR] ->Fill(ecluster);
2489 fhMCPt [kmcISR] ->Fill(ptcluster);
2490 fhMCPhi[kmcISR] ->Fill(ecluster,phicluster);
2491 fhMCEta[kmcISR] ->Fill(ecluster,etacluster);
4c8f7c2e 2492
c5693f62 2493 fhMC2E[kmcISR] ->Fill(ecluster, eprim);
2494 fhMC2Pt[kmcISR] ->Fill(ptcluster, ptprim);
2495 fhMCDeltaE[kmcISR] ->Fill(ecluster,eprim-ecluster);
2496 fhMCDeltaPt[kmcISR]->Fill(ptcluster,ptprim-ptcluster);
4c8f7c2e 2497
3d5d5078 2498 }
2499 else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay) &&
c5693f62 2500 !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE[kmcPi0Decay])
3d5d5078 2501 {
c5693f62 2502 fhMCE [kmcPi0Decay] ->Fill(ecluster);
2503 fhMCPt [kmcPi0Decay] ->Fill(ptcluster);
2504 fhMCPhi[kmcPi0Decay] ->Fill(ecluster,phicluster);
2505 fhMCEta[kmcPi0Decay] ->Fill(ecluster,etacluster);
4c8f7c2e 2506
c5693f62 2507 fhMC2E[kmcPi0Decay] ->Fill(ecluster, eprim);
2508 fhMC2Pt[kmcPi0Decay] ->Fill(ptcluster, ptprim);
2509 fhMCDeltaE[kmcPi0Decay] ->Fill(ecluster,eprim-ecluster);
2510 fhMCDeltaPt[kmcPi0Decay]->Fill(ptcluster,ptprim-ptcluster);
4c8f7c2e 2511
3d5d5078 2512 }
f586f4aa 2513 else if( (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) ||
c5693f62 2514 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay) ) && fhMCE[kmcOtherDecay])
3d5d5078 2515 {
c5693f62 2516 fhMCE [kmcOtherDecay] ->Fill(ecluster);
2517 fhMCPt [kmcOtherDecay] ->Fill(ptcluster);
2518 fhMCPhi[kmcOtherDecay] ->Fill(ecluster,phicluster);
2519 fhMCEta[kmcOtherDecay] ->Fill(ecluster,etacluster);
4c8f7c2e 2520
c5693f62 2521 fhMC2E[kmcOtherDecay] ->Fill(ecluster, eprim);
2522 fhMC2Pt[kmcOtherDecay] ->Fill(ptcluster, ptprim);
2523 fhMCDeltaE[kmcOtherDecay] ->Fill(ecluster,eprim-ecluster);
2524 fhMCDeltaPt[kmcOtherDecay]->Fill(ptcluster,ptprim-ptcluster);
4c8f7c2e 2525
3d5d5078 2526 }
c5693f62 2527 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE [kmcPi0])
3d5d5078 2528 {
c5693f62 2529 fhMCE [kmcPi0] ->Fill(ecluster);
2530 fhMCPt [kmcPi0] ->Fill(ptcluster);
2531 fhMCPhi[kmcPi0] ->Fill(ecluster,phicluster);
2532 fhMCEta[kmcPi0] ->Fill(ecluster,etacluster);
4c8f7c2e 2533
c5693f62 2534 fhMC2E[kmcPi0] ->Fill(ecluster, eprim);
2535 fhMC2Pt[kmcPi0] ->Fill(ptcluster, ptprim);
2536 fhMCDeltaE[kmcPi0] ->Fill(ecluster,eprim-ecluster);
2537 fhMCDeltaPt[kmcPi0]->Fill(ptcluster,ptprim-ptcluster);
4c8f7c2e 2538
f66d95af 2539 }
c5693f62 2540 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) && fhMCE[kmcEta])
f66d95af 2541 {
c5693f62 2542 fhMCE [kmcEta] ->Fill(ecluster);
2543 fhMCPt [kmcEta] ->Fill(ptcluster);
2544 fhMCPhi[kmcEta] ->Fill(ecluster,phicluster);
2545 fhMCEta[kmcEta] ->Fill(ecluster,etacluster);
4c8f7c2e 2546
c5693f62 2547 fhMC2E[kmcEta] ->Fill(ecluster, eprim);
2548 fhMC2Pt[kmcEta] ->Fill(ptcluster, ptprim);
2549 fhMCDeltaE[kmcEta] ->Fill(ecluster,eprim-ecluster);
2550 fhMCDeltaPt[kmcEta]->Fill(ptcluster,ptprim-ptcluster);
4c8f7c2e 2551
f66d95af 2552 }
3d5d5078 2553 }
c5693f62 2554 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron) && fhMCE[kmcAntiNeutron])
3d5d5078 2555 {
c5693f62 2556 fhMCE [kmcAntiNeutron] ->Fill(ecluster);
2557 fhMCPt [kmcAntiNeutron] ->Fill(ptcluster);
2558 fhMCPhi[kmcAntiNeutron] ->Fill(ecluster,phicluster);
2559 fhMCEta[kmcAntiNeutron] ->Fill(ecluster,etacluster);
3d5d5078 2560
c5693f62 2561 fhMC2E[kmcAntiNeutron] ->Fill(ecluster, eprim);
2562 fhMC2Pt[kmcAntiNeutron] ->Fill(ptcluster, ptprim);
2563 fhMCDeltaE[kmcAntiNeutron] ->Fill(ecluster,eprim-ecluster);
2564 fhMCDeltaPt[kmcAntiNeutron]->Fill(ptcluster,ptprim-ptcluster);
4c8f7c2e 2565
3d5d5078 2566 }
c5693f62 2567 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton) && fhMCE[kmcAntiProton])
3d5d5078 2568 {
c5693f62 2569 fhMCE [kmcAntiProton] ->Fill(ecluster);
2570 fhMCPt [kmcAntiProton] ->Fill(ptcluster);
2571 fhMCPhi[kmcAntiProton] ->Fill(ecluster,phicluster);
2572 fhMCEta[kmcAntiProton] ->Fill(ecluster,etacluster);
4c8f7c2e 2573
c5693f62 2574 fhMC2E[kmcAntiProton] ->Fill(ecluster, eprim);
2575 fhMC2Pt[kmcAntiProton] ->Fill(ptcluster, ptprim);
2576 fhMCDeltaE[kmcAntiProton] ->Fill(ecluster,eprim-ecluster);
2577 fhMCDeltaPt[kmcAntiProton]->Fill(ecluster,ptprim-ptcluster);
4c8f7c2e 2578
3d5d5078 2579 }
c5693f62 2580 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron) && fhMCE[kmcElectron])
3d5d5078 2581 {
c5693f62 2582 fhMCE [kmcElectron] ->Fill(ecluster);
2583 fhMCPt [kmcElectron] ->Fill(ptcluster);
2584 fhMCPhi[kmcElectron] ->Fill(ecluster,phicluster);
2585 fhMCEta[kmcElectron] ->Fill(ecluster,etacluster);
4c8f7c2e 2586
c5693f62 2587 fhMC2E[kmcElectron] ->Fill(ecluster, eprim);
2588 fhMC2Pt[kmcElectron] ->Fill(ptcluster, ptprim);
2589 fhMCDeltaE[kmcElectron] ->Fill(ecluster,eprim-ecluster);
2590 fhMCDeltaPt[kmcElectron]->Fill(ecluster,ptprim-ptcluster);
3d5d5078 2591 }
34c16486 2592 else if( fhMCE[kmcOther])
2593 {
c5693f62 2594 fhMCE [kmcOther] ->Fill(ecluster);
2595 fhMCPt [kmcOther] ->Fill(ptcluster);
2596 fhMCPhi[kmcOther] ->Fill(ecluster,phicluster);
2597 fhMCEta[kmcOther] ->Fill(ecluster,etacluster);
521636d2 2598
c5693f62 2599 fhMC2E[kmcOther] ->Fill(ecluster, eprim);
2600 fhMC2Pt[kmcOther] ->Fill(ptcluster, ptprim);
2601 fhMCDeltaE[kmcOther] ->Fill(ecluster,eprim-ecluster);
2602 fhMCDeltaPt[kmcOther]->Fill(ecluster,ptprim-ptcluster);
4c8f7c2e 2603
f8006433 2604 // printf(" AliAnaPhoton::MakeAnalysisFillHistograms() - Label %d, pT %2.3f Unknown, bits set: ",
2605 // ph->GetLabel(),ph->Pt());
2606 // for(Int_t i = 0; i < 20; i++) {
2607 // if(GetMCAnalysisUtils()->CheckTagBit(tag,i)) printf(" %d, ",i);
2608 // }
2609 // printf("\n");
2610
577d9801 2611 }
521636d2 2612
577d9801 2613 }//Histograms with MC
521636d2 2614
577d9801 2615 }// aod loop
521636d2 2616
1c5acb87 2617}
2618
2619
2620//__________________________________________________________________
2621void AliAnaPhoton::Print(const Option_t * opt) const
2622{
477d6cee 2623 //Print some relevant parameters set for the analysis
2624
2625 if(! opt)
2626 return;
2627
2628 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
745913ae 2629 AliAnaCaloTrackCorrBaseClass::Print(" ");
a3aebfff 2630
477d6cee 2631 printf("Calorimeter = %s\n", fCalorimeter.Data()) ;
2632 printf("Min Distance to Bad Channel = %2.1f\n",fMinDist);
2633 printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
2634 printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
a3aebfff 2635 printf("Reject clusters with a track matched = %d\n",fRejectTrackMatch);
4cf55759 2636 printf("Time Cut: %3.1f < TOF < %3.1f\n", fTimeCutMin, fTimeCutMax);
2ac125bf 2637 printf("Number of cells in cluster is > %d \n", fNCellsCut);
477d6cee 2638 printf(" \n") ;
1c5acb87 2639
2640}