AliPhysics  b555aef (b555aef)
AliAnaElectron.cxx
Go to the documentation of this file.
1  /**************************************************************************
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 *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
16 // --- ROOT system ---
17 #include <TH2F.h>
18 #include <TH3D.h>
19 #include <TClonesArray.h>
20 #include <TObjString.h>
21 //#include <Riostream.h>
22 #include "TDatabasePDG.h"
23 #include "AliVTrack.h"
24 
25 // --- Analysis system ---
26 #include "AliAnaElectron.h"
27 #include "AliCaloTrackReader.h"
28 #include "AliMCEvent.h"
29 #include "AliCaloPID.h"
30 #include "AliMCAnalysisUtils.h"
31 #include "AliFiducialCut.h"
32 #include "AliVCluster.h"
33 #include "AliVParticle.h"
34 #include "AliMixedEvent.h"
35 
37 ClassImp(AliAnaElectron) ;
39 
40 //________________________________
42 //________________________________
45  fMinDist(0.), fMinDist2(0.), fMinDist3(0.),
46  fTimeCutMin(-1), fTimeCutMax(999999),
47  fNCellsCut(0), fNLMCutMin(-1), fNLMCutMax(10),
48  fFillSSHistograms(kFALSE), fFillOnlySimpleSSHisto(1),
49  fFillWeightHistograms(kFALSE), fNOriginHistograms(8),
50  fdEdxMin(0.), fdEdxMax (200.),
51  fEOverPMin(0), fEOverPMax (2),
52  fAODParticle(0),
53  fMomentum(), fMomentumMC(), fProdVertex(),
54  // Histograms
55  fhdEdxvsE(0), fhdEdxvsP(0),
56  fhEOverPvsE(0), fhEOverPvsP(0),
57  fhdEdxvsECutM02(0), fhdEdxvsPCutM02(0),
58  fhEOverPvsECutM02(0), fhEOverPvsPCutM02(0),
59  fhdEdxvsECutEOverP(0), fhdEdxvsPCutEOverP(0),
60  fhEOverPvsECutM02CutdEdx(0), fhEOverPvsPCutM02CutdEdx(0),
61  // Weight studies
62  fhECellClusterRatio(0), fhECellClusterLogRatio(0),
63  fhEMaxCellClusterRatio(0), fhEMaxCellClusterLogRatio(0),
64  // MC histograms
65  // Electron SS MC histograms
66  fhMCElectronELambda0NoOverlap(0),
67  fhMCElectronELambda0TwoOverlap(0), fhMCElectronELambda0NOverlap(0),
68 
69  //Embedding
70  fhEmbeddedSignalFractionEnergy(0),
71  fhEmbedElectronELambda0FullSignal(0), fhEmbedElectronELambda0MostlySignal(0),
72  fhEmbedElectronELambda0MostlyBkg(0), fhEmbedElectronELambda0FullBkg(0)
73 {
74  for(Int_t index = 0; index < 2; index++)
75  {
76  fhNCellsE [index] = 0;
77  fhNLME [index] = 0;
78  fhTimeE [index] = 0;
79  fhMaxCellDiffClusterE[index] = 0;
80  fhE [index] = 0;
81  fhPt [index] = 0;
82  fhPhi [index] = 0;
83  fhEta [index] = 0;
84  fhEtaPhi [index] = 0;
85  fhEtaPhi05[index] = 0;
86 
87  // Shower shape histograms
88  fhDispE [index] = 0;
89  fhLam0E [index] = 0;
90  fhLam1E [index] = 0;
91  fhDispETRD[index] = 0;
92  fhLam0ETRD[index] = 0;
93  fhLam1ETRD[index] = 0;
94  fhNCellsLam0LowE [index] = 0;
95  fhNCellsLam0HighE[index] = 0;
96  fhEtaLam0LowE [index] = 0;
97  fhPhiLam0LowE [index] = 0;
98  fhEtaLam0HighE [index] = 0;
99  fhPhiLam0HighE [index] = 0;
100 
101  fhDispEtaE [index] = 0;
102  fhDispPhiE [index] = 0;
103  fhSumEtaE [index] = 0;
104  fhSumPhiE [index] = 0;
105  fhSumEtaPhiE [index] = 0;
106  fhDispEtaPhiDiffE[index] = 0;
107  fhSphericityE [index] = 0;
108 
109  for(Int_t i = 0; i < 10; i++)
110  {
111  fhMCPt [index][i] = 0;
112  fhMCE [index][i] = 0;
113  fhMCPhi [index][i] = 0;
114  fhMCEta [index][i] = 0;
115  fhMCDeltaE [index][i] = 0;
116  fhMC2E [index][i] = 0;
117  fhMCdEdxvsE [i] = 0;
118  fhMCdEdxvsP [i] = 0;
119  fhMCEOverPvsE [i] = 0;
120  fhMCEOverPvsP [i] = 0;
121  }
122 
123  for(Int_t i = 0; i < 6; i++)
124  {
125  fhMCELambda0 [index][i] = 0;
126  fhMCEDispEta [index][i] = 0;
127  fhMCEDispPhi [index][i] = 0;
128  fhMCESumEtaPhi [index][i] = 0;
129  fhMCEDispEtaPhiDiff[index][i] = 0;
130  fhMCESphericity [index][i] = 0;
131  }
132 
133  for(Int_t i = 0; i < 5; i++)
134  {
135  fhDispEtaDispPhiEBin[index][i] = 0 ;
136  }
137  }
138 
139  // Weight studies
140  for(Int_t i =0; i < 14; i++)
141  {
142  fhLambda0ForW0[i] = 0;
143  //fhLambda1ForW0[i] = 0;
144  }
145 
146  // Initialize parameters
147  InitParameters();
148 }
149 
150 //____________________________________________________________________________
166 //____________________________________________________________________________
167 Bool_t AliAnaElectron::ClusterSelected(AliVCluster* calo, Int_t nMaxima)
168 {
169  AliDebug(1,Form("Current Event %d; Before selection : E %2.2f, pT %2.2f, Ecl %2.2f, phi %2.2f, eta %2.2f",
170  GetReader()->GetEventNumber(),fMomentum.E(),fMomentum.Pt(),calo->E(),fMomentum.Phi()*TMath::RadToDeg(),fMomentum.Eta()));
171 
172  //.......................................
173  //If too small or big energy, skip it
174  if(fMomentum.E() < GetMinEnergy() || fMomentum.E() > GetMaxEnergy() ) return kFALSE ;
175  AliDebug(2,Form("\t Cluster %d Pass E Cut",calo->GetID()));
176 
177  //.......................................
178  // TOF cut, BE CAREFUL WITH THIS CUT
179  Double_t tof = calo->GetTOF()*1e9;
180  if(tof < fTimeCutMin || tof > fTimeCutMax) return kFALSE;
181  AliDebug(2,Form("\t Cluster %d Pass Time Cut",calo->GetID()));
182 
183  //.......................................
184  if(calo->GetNCells() <= fNCellsCut && GetReader()->GetDataType() != AliCaloTrackReader::kMC) return kFALSE;
185  AliDebug(2,Form("\t Cluster %d Pass NCell Cut",calo->GetID()));
186 
187  //.......................................
188  //Check acceptance selection
189  if(IsFiducialCutOn())
190  {
192  if(! in ) return kFALSE ;
193  }
194  AliDebug(2,"\t Fiducial cut passed");
195 
196  //.......................................
197  //Skip not matched clusters with tracks
198  if(!IsTrackMatched(calo, GetReader()->GetInputEvent()))
199  {
200  AliDebug(1,"\t Reject non track-matched clusters");
201  return kFALSE ;
202  }
203  else AliDebug(2,"\t Track-matching cut passed");
204 
205  //...........................................
206  // skip clusters with too many maxima
207  if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax) return kFALSE ;
208  AliDebug(2,Form("\t Cluster %d pass NLM %d of out of range",calo->GetID(), nMaxima));
209 
210  //.......................................
211  //Check Distance to Bad channel, set bit.
212  Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
213  if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
214  if(distBad < fMinDist) {//In bad channel (PHOS cristal size 2.2x2.2 cm), EMCAL ( cell units )
215  return kFALSE ;
216  }
217  else AliDebug(2,Form("\t Bad channel cut passed %4.2f > %2.2f",distBad, fMinDist));
218  //printf("Cluster %d Pass Bad Dist Cut \n",icalo);
219 
220  AliDebug(1,Form("Current Event %d; After selection : E %2.2f, pT %2.2f, Ecl %2.2f, phi %2.2f, eta %2.2f",
222  fMomentum.E(), fMomentum.Pt(),calo->E(),fMomentum.Phi()*TMath::RadToDeg(),fMomentum.Eta()));
223 
224  //All checks passed, cluster selected
225  return kTRUE;
226 
227 }
228 
229 //______________________________________________________________________________________________
231 //______________________________________________________________________________________________
232 void AliAnaElectron::FillShowerShapeHistograms(AliVCluster* cluster, Int_t mcTag, Int_t pidTag)
233 {
234  if(!fFillSSHistograms || GetMixedEvent()) return;
235 
236  Int_t pidIndex = 0;// Electron
237  if (pidTag == AliCaloPID::kElectron) pidIndex = 0;
238  else if(pidTag == AliCaloPID::kChargedHadron) pidIndex = 1;
239  else return;
240 
241  Float_t energy = cluster->E();
242  Int_t ncells = cluster->GetNCells();
243  Float_t lambda0 = cluster->GetM02();
244  Float_t lambda1 = cluster->GetM20();
245  Float_t disp = cluster->GetDispersion()*cluster->GetDispersion();
246 
247  Float_t l0 = 0., l1 = 0.;
248  Float_t dispp= 0., dEta = 0., dPhi = 0.;
249  Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
250 
251  Float_t eta = fMomentum.Eta();
252  Float_t phi = fMomentum.Phi();
253  if(phi < 0) phi+=TMath::TwoPi();
254 
255  fhLam0E[pidIndex] ->Fill(energy, lambda0, GetEventWeight());
256  fhLam1E[pidIndex] ->Fill(energy, lambda1, GetEventWeight());
257  fhDispE[pidIndex] ->Fill(energy, disp , GetEventWeight());
258 
259  if(GetCalorimeter() == kEMCAL && GetFirstSMCoveredByTRD() >= 0 &&
261  {
262  fhLam0ETRD[pidIndex]->Fill(energy, lambda0, GetEventWeight());
263  fhLam1ETRD[pidIndex]->Fill(energy, lambda1, GetEventWeight());
264  fhDispETRD[pidIndex]->Fill(energy, disp , GetEventWeight());
265  }
266 
268  {
269  if(energy < 2)
270  {
271  fhNCellsLam0LowE[pidIndex] ->Fill(ncells, lambda0, GetEventWeight());
272  fhEtaLam0LowE[pidIndex] ->Fill(eta, lambda0, GetEventWeight());
273  fhPhiLam0LowE[pidIndex] ->Fill(phi, lambda0, GetEventWeight());
274  }
275  else
276  {
277  fhNCellsLam0HighE[pidIndex]->Fill(ncells, lambda0, GetEventWeight());
278  fhEtaLam0HighE[pidIndex] ->Fill(eta, lambda0, GetEventWeight());
279  fhPhiLam0HighE[pidIndex] ->Fill(phi, lambda0, GetEventWeight());
280  }
281 
282  if(GetCalorimeter() == kEMCAL)
283  {
284  GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), GetReader()->GetInputEvent()->GetEMCALCells(), cluster,
285  l0, l1, dispp, dEta, dPhi, sEta, sPhi, sEtaPhi);
286  fhDispEtaE [pidIndex]-> Fill(energy, dEta , GetEventWeight());
287  fhDispPhiE [pidIndex]-> Fill(energy, dPhi , GetEventWeight());
288  fhSumEtaE [pidIndex]-> Fill(energy, sEta , GetEventWeight());
289  fhSumPhiE [pidIndex]-> Fill(energy, sPhi , GetEventWeight());
290  fhSumEtaPhiE [pidIndex]-> Fill(energy, sEtaPhi , GetEventWeight());
291  fhDispEtaPhiDiffE [pidIndex]-> Fill(energy, dPhi-dEta, GetEventWeight());
292  if(dEta+dPhi>0)
293  fhSphericityE [pidIndex]-> Fill(energy, (dPhi-dEta)/(dEta+dPhi), GetEventWeight());
294 
295  if (energy < 2 ) fhDispEtaDispPhiEBin[pidIndex][0]->Fill(dEta, dPhi, GetEventWeight());
296  else if (energy < 4 ) fhDispEtaDispPhiEBin[pidIndex][1]->Fill(dEta, dPhi, GetEventWeight());
297  else if (energy < 6 ) fhDispEtaDispPhiEBin[pidIndex][2]->Fill(dEta, dPhi, GetEventWeight());
298  else if (energy < 10) fhDispEtaDispPhiEBin[pidIndex][3]->Fill(dEta, dPhi, GetEventWeight());
299  else fhDispEtaDispPhiEBin[pidIndex][4]->Fill(dEta, dPhi, GetEventWeight());
300  }
301  }
302 
303  if(IsDataMC())
304  {
305  AliVCaloCells* cells = 0;
306  if(GetCalorimeter() == kEMCAL) cells = GetEMCALCells();
307  else cells = GetPHOSCells();
308 
309  //Fill histograms to check shape of embedded clusters
310  Float_t fraction = 0;
311  if(GetReader()->IsEmbeddedClusterSelectionOn()){//Only working for EMCAL
312 
313  Float_t clusterE = 0; // recalculate in case corrections applied.
314  Float_t cellE = 0;
315  for(Int_t icell = 0; icell < cluster->GetNCells(); icell++){
316  cellE = cells->GetCellAmplitude(cluster->GetCellAbsId(icell));
317  clusterE+=cellE;
318  fraction+=cellE*cluster->GetCellAmplitudeFraction(icell);
319  }
320 
321  //Fraction of total energy due to the embedded signal
322  fraction/=clusterE;
323 
324  AliDebug(1,Form("Energy fraction of embedded signal %2.3f, Energy %2.3f",fraction, clusterE));
325 
326  fhEmbeddedSignalFractionEnergy->Fill(clusterE, fraction, GetEventWeight());
327  } // embedded fraction
328 
329  // Check the origin and fill histograms
330  Int_t index = -1;
331 
332  if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton) &&
333  !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion) &&
334  !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0) &&
335  !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta))
336  {
337  index = kmcssPhoton;
338  }//photon no conversion
339  else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron &&
341  {
342  index = kmcssElectron;
343 
344  if(!GetReader()->IsEmbeddedClusterSelectionOn())
345  {
346  //Check particle overlaps in cluster
347 
348  //Compare the primary depositing more energy with the rest, if no photon/electron as comon ancestor (conversions), count as other particle
349  Int_t ancPDG = 0, ancStatus = -1;
350  Int_t ancLabel = 0;
351  Int_t noverlaps = 1;
352  for (UInt_t ilab = 0; ilab < cluster->GetNLabels(); ilab++ )
353  {
354  ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(cluster->GetLabels()[0],cluster->GetLabels()[ilab], GetMC(),
355  ancPDG,ancStatus,fMomentumMC,fProdVertex);
356  if(ancPDG!=22 && TMath::Abs(ancPDG)!=11) noverlaps++;
357  }
358 
359  if(noverlaps == 1){
360  fhMCElectronELambda0NoOverlap ->Fill(energy, lambda0, GetEventWeight());
361  }
362  else if(noverlaps == 2){
363  fhMCElectronELambda0TwoOverlap ->Fill(energy, lambda0, GetEventWeight());
364  }
365  else if(noverlaps > 2){
366  fhMCElectronELambda0NOverlap ->Fill(energy, lambda0, GetEventWeight());
367  }
368  else
369  {
370  AliWarning(Form("N overlaps = %d for ancestor %d!!", noverlaps, ancLabel));
371  }
372  }//No embedding
373 
374  //Fill histograms to check shape of embedded clusters
375  if(GetReader()->IsEmbeddedClusterSelectionOn())
376  {
377  if (fraction > 0.9)
378  {
379  fhEmbedElectronELambda0FullSignal ->Fill(energy, lambda0, GetEventWeight());
380  }
381  else if(fraction > 0.5)
382  {
383  fhEmbedElectronELambda0MostlySignal ->Fill(energy, lambda0, GetEventWeight());
384  }
385  else if(fraction > 0.1)
386  {
387  fhEmbedElectronELambda0MostlyBkg ->Fill(energy, lambda0, GetEventWeight());
388  }
389  else
390  {
391  fhEmbedElectronELambda0FullBkg ->Fill(energy, lambda0, GetEventWeight());
392  }
393  } // embedded
394  } // electron
395  else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron) &&
397  {
398  index = kmcssConversion;
399  } // conversion photon
400  else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0) )
401  {
402  index = kmcssPi0;
403  } // pi0
404  else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta) )
405  {
406  index = kmcssEta;
407  } // eta
408  else
409  {
410  index = kmcssOther;
411  } // other particles
412 
413  fhMCELambda0[pidIndex][index] ->Fill(energy, lambda0, GetEventWeight());
414 
416  {
417  fhMCEDispEta [pidIndex][index]-> Fill(energy, dEta , GetEventWeight());
418  fhMCEDispPhi [pidIndex][index]-> Fill(energy, dPhi , GetEventWeight());
419  fhMCESumEtaPhi [pidIndex][index]-> Fill(energy, sEtaPhi , GetEventWeight());
420  fhMCEDispEtaPhiDiff [pidIndex][index]-> Fill(energy, dPhi-dEta, GetEventWeight());
421  if(dEta+dPhi>0)
422  fhMCESphericity [pidIndex][index]-> Fill(energy, (dPhi-dEta)/(dEta+dPhi), GetEventWeight());
423  }
424  } // MC data
425 }
426 
427 //_____________________________________________
429 //_____________________________________________
431 {
432  TString parList ; //this will be list of parameters used for this analysis.
433  const Int_t buffersize = 255;
434  char onePar[buffersize] ;
435 
436  snprintf(onePar,buffersize,"--- AliAnaElectron ---: ") ;
437  parList+=onePar ;
438  snprintf(onePar,buffersize,"Calorimeter: %s;",GetCalorimeterString().Data()) ;
439  parList+=onePar ;
440  snprintf(onePar,buffersize," %2.2f < dEdx < %2.2f;",fdEdxMin,fdEdxMax) ;
441  parList+=onePar ;
442  snprintf(onePar,buffersize," %2.2f < E/P < %2.2f;",fEOverPMin, fEOverPMax) ;
443  parList+=onePar ;
444  snprintf(onePar,buffersize,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster);",fMinDist) ;
445  parList+=onePar ;
446  snprintf(onePar,buffersize,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation);",fMinDist2) ;
447  parList+=onePar ;
448  snprintf(onePar,buffersize,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study);",fMinDist3) ;
449  parList+=onePar ;
450 
451  //Get parameters set in base class.
452  parList += GetBaseParametersList() ;
453 
454  //Get parameters set in PID class.
455  parList += GetCaloPID()->GetPIDParametersList() ;
456 
457  //Get parameters set in FiducialCut class (not available yet)
458  //parlist += GetFidCut()->GetFidCutParametersList()
459 
460  return new TObjString(parList) ;
461 }
462 
463 //_______________________________________________
464 // Create histograms to be saved in output file and
465 // store them in outputContainer.
466 //_______________________________________________
468 {
469  TList * outputContainer = new TList() ;
470  outputContainer->SetName("ElectronHistos") ;
471 
480 
481  // MC labels, titles, for originator particles
482  TString ptypess[] = { "#gamma","hadron?","#pi^{0}","#eta","#gamma->e^{#pm}","e^{#pm}"} ;
483  TString pnamess[] = { "Photon","Hadron" ,"Pi0" ,"Eta" ,"Conversion" ,"Electron"} ;
484  TString ptype[] = { "#gamma", "#gamma_{#pi decay}","#gamma_{other decay}", "#pi^{0}","#eta",
485  "e^{#pm}","#gamma->e^{#pm}","hadron?","Anti-N","Anti-P" } ;
486 
487  TString pname[] = { "Photon","PhotonPi0Decay","PhotonOtherDecay","Pi0","Eta","Electron",
488  "Conversion", "Hadron", "AntiNeutron","AntiProton" } ;
489 
490  fhdEdxvsE = new TH2F ("hdEdxvsE","matched track <dE/dx> vs cluster E ", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
491  fhdEdxvsE->SetXTitle("E (GeV)");
492  fhdEdxvsE->SetYTitle("<dE/dx>");
493  outputContainer->Add(fhdEdxvsE);
494 
495  fhdEdxvsP = new TH2F ("hdEdxvsP","matched track <dE/dx> vs track P ", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
496  fhdEdxvsP->SetXTitle("P (GeV/c)");
497  fhdEdxvsP->SetYTitle("<dE/dx>");
498  outputContainer->Add(fhdEdxvsP);
499 
500  fhEOverPvsE = new TH2F ("hEOverPvsE","matched track E/p vs cluster E ", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
501  fhEOverPvsE->SetXTitle("E (GeV)");
502  fhEOverPvsE->SetYTitle("E/p");
503  outputContainer->Add(fhEOverPvsE);
504 
505  fhEOverPvsP = new TH2F ("hEOverPvsP","matched track E/p vs track P ", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
506  fhEOverPvsP->SetXTitle("P (GeV/c)");
507  fhEOverPvsP->SetYTitle("E/p");
508  outputContainer->Add(fhEOverPvsP);
509 
510 
511  fhdEdxvsECutM02 = new TH2F ("hdEdxvsECutM02","matched track <dE/dx> vs cluster E, mild #lambda_{0}^{2} cut", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
512  fhdEdxvsECutM02->SetXTitle("E (GeV)");
513  fhdEdxvsECutM02->SetYTitle("<dE/dx>");
514  outputContainer->Add(fhdEdxvsECutM02);
515 
516  fhdEdxvsPCutM02 = new TH2F ("hdEdxvsPCutM02","matched track <dE/dx> vs track P, mild #lambda_{0}^{2} cut", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
517  fhdEdxvsPCutM02->SetXTitle("P (GeV/c)");
518  fhdEdxvsPCutM02->SetYTitle("<dE/dx>");
519  outputContainer->Add(fhdEdxvsPCutM02);
520 
521  fhEOverPvsECutM02 = new TH2F ("hEOverPvsECutM02","matched track E/p vs cluster E, mild #lambda_{0}^{2} cut", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
522  fhEOverPvsECutM02->SetXTitle("E (GeV)");
523  fhEOverPvsECutM02->SetYTitle("E/p");
524  outputContainer->Add(fhEOverPvsECutM02);
525 
526  fhEOverPvsPCutM02 = new TH2F ("hEOverPvsPCutM02","matched track E/p vs track P, mild #lambda_{0}^{2} cut", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
527  fhEOverPvsPCutM02->SetXTitle("P (GeV/c)");
528  fhEOverPvsPCutM02->SetYTitle("E/p");
529  outputContainer->Add(fhEOverPvsPCutM02);
530 
531 
532  fhdEdxvsECutEOverP = new TH2F ("hdEdxvsECutEOverP","matched track <dE/dx> vs cluster E, cut on E/p", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
533  fhdEdxvsECutEOverP->SetXTitle("E (GeV)");
534  fhdEdxvsECutEOverP->SetYTitle("<dE/dx>");
535  outputContainer->Add(fhdEdxvsECutEOverP);
536 
537  fhdEdxvsPCutEOverP = new TH2F ("hdEdxvsPCutEOverP","matched track <dE/dx> vs track P, cut on E/p", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
538  fhdEdxvsPCutEOverP->SetXTitle("P (GeV/c)");
539  fhdEdxvsPCutEOverP->SetYTitle("<dE/dx>");
540  outputContainer->Add(fhdEdxvsPCutEOverP);
541 
542  fhEOverPvsECutM02CutdEdx = new TH2F ("hEOverPvsECutM02CutdEdx","matched track E/p vs cluster E, dEdx cut, mild #lambda_{0}^{2} cut", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
543  fhEOverPvsECutM02CutdEdx->SetXTitle("E (GeV)");
544  fhEOverPvsECutM02CutdEdx->SetYTitle("E/p");
545  outputContainer->Add(fhEOverPvsECutM02CutdEdx);
546 
547  fhEOverPvsPCutM02CutdEdx = new TH2F ("hEOverPvsPCutM02CutdEdx","matched track E/p vs track P, dEdx cut, mild #lambda_{0}^{2} cut ", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
548  fhEOverPvsPCutM02CutdEdx->SetXTitle("P (GeV/c)");
549  fhEOverPvsPCutM02CutdEdx->SetYTitle("E/p");
550  outputContainer->Add(fhEOverPvsPCutM02CutdEdx);
551 
552  if(IsDataMC())
553  {
554  for(Int_t i = 0; i < fNOriginHistograms; i++)
555  {
556  fhMCdEdxvsE[i] = new TH2F(Form("hdEdxvsE_MC%s",pname[i].Data()),
557  Form("matched track <dE/dx> vs cluster E from %s : E ",ptype[i].Data()),
558  nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
559  fhMCdEdxvsE[i]->SetXTitle("E (GeV)");
560  fhMCdEdxvsE[i]->SetYTitle("<dE/dx>");
561  outputContainer->Add(fhMCdEdxvsE[i]) ;
562 
563  fhMCdEdxvsP[i] = new TH2F(Form("hdEdxvsP_MC%s",pname[i].Data()),
564  Form("matched track <dE/dx> vs track P from %s : E ",ptype[i].Data()),
565  nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
566  fhMCdEdxvsP[i]->SetXTitle("E (GeV)");
567  fhMCdEdxvsP[i]->SetYTitle("<dE/dx>");
568  outputContainer->Add(fhMCdEdxvsP[i]) ;
569 
570 
571  fhMCEOverPvsE[i] = new TH2F(Form("hEOverPvsE_MC%s",pname[i].Data()),
572  Form("matched track E/p vs cluster E from %s : E ",ptype[i].Data()),
573  nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
574  fhMCEOverPvsE[i]->SetXTitle("E (GeV)");
575  fhMCEOverPvsE[i]->SetYTitle("<dE/dx>");
576  outputContainer->Add(fhMCEOverPvsE[i]) ;
577 
578  fhMCEOverPvsP[i] = new TH2F(Form("hEOverPvsP_MC%s",pname[i].Data()),
579  Form("matched track E/pvs track P from %s : E ",ptype[i].Data()),
580  nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
581  fhMCEOverPvsP[i]->SetXTitle("E (GeV)");
582  fhMCEOverPvsP[i]->SetYTitle("<dE/dx>");
583  outputContainer->Add(fhMCEOverPvsP[i]) ;
584  }
585  }
586 
587  TString pidParticle[] = {"Electron","ChargedHadron"} ;
588 
590  {
591  fhECellClusterRatio = new TH2F ("hECellClusterRatio"," cell energy / cluster energy vs cluster energy, for selected electrons",
592  nptbins,ptmin,ptmax, 100,0,1.);
593  fhECellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
594  fhECellClusterRatio->SetYTitle("E_{cell i}/E_{cluster}");
595  outputContainer->Add(fhECellClusterRatio);
596 
597  fhECellClusterLogRatio = new TH2F ("hECellClusterLogRatio"," Log(cell energy / cluster energy) vs cluster energy, for selected electrons",
598  nptbins,ptmin,ptmax, 100,-10,0);
599  fhECellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
600  fhECellClusterLogRatio->SetYTitle("Log (E_{max cell}/E_{cluster})");
601  outputContainer->Add(fhECellClusterLogRatio);
602 
603  fhEMaxCellClusterRatio = new TH2F ("hEMaxCellClusterRatio"," max cell energy / cluster energy vs cluster energy, for selected electrons",
604  nptbins,ptmin,ptmax, 100,0,1.);
605  fhEMaxCellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
606  fhEMaxCellClusterRatio->SetYTitle("E_{max cell}/E_{cluster}");
607  outputContainer->Add(fhEMaxCellClusterRatio);
608 
609  fhEMaxCellClusterLogRatio = new TH2F ("hEMaxCellClusterLogRatio"," Log(max cell energy / cluster energy) vs cluster energy, for selected electrons",
610  nptbins,ptmin,ptmax, 100,-10,0);
611  fhEMaxCellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
612  fhEMaxCellClusterLogRatio->SetYTitle("Log (E_{max cell}/E_{cluster})");
613  outputContainer->Add(fhEMaxCellClusterLogRatio);
614 
615  for(Int_t iw = 0; iw < 14; iw++)
616  {
617  fhLambda0ForW0[iw] = new TH2F (Form("hLambda0ForW0%d",iw),Form("shower shape, #lambda^{2}_{0} vs E, w0 = %1.1f, for selected electrons",1+0.5*iw),
618  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
619  fhLambda0ForW0[iw]->SetXTitle("E_{cluster}");
620  fhLambda0ForW0[iw]->SetYTitle("#lambda^{2}_{0}");
621  outputContainer->Add(fhLambda0ForW0[iw]);
622 
623  // fhLambda1ForW0[iw] = new TH2F (Form("hLambda1ForW0%d",iw),Form("shower shape, #lambda^{2}_{1} vs E, w0 = %1.1f, for selected electrons",1+0.5*iw),
624  // nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
625  // fhLambda1ForW0[iw]->SetXTitle("E_{cluster}");
626  // fhLambda1ForW0[iw]->SetYTitle("#lambda^{2}_{1}");
627  // outputContainer->Add(fhLambda1ForW0[iw]);
628  }
629  }
630 
631  for(Int_t pidIndex = 0; pidIndex < 2; pidIndex++)
632  {
633  // Shower shape
635  {
636  fhLam0E[pidIndex] = new TH2F (Form("h%sLam0E",pidParticle[pidIndex].Data()),
637  Form("%s: #lambda_{0}^{2} vs E",pidParticle[pidIndex].Data()),
638  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
639  fhLam0E[pidIndex]->SetYTitle("#lambda_{0}^{2}");
640  fhLam0E[pidIndex]->SetXTitle("E (GeV)");
641  outputContainer->Add(fhLam0E[pidIndex]);
642 
643  fhLam1E[pidIndex] = new TH2F (Form("h%sLam1E",pidParticle[pidIndex].Data()),
644  Form("%s: #lambda_{1}^{2} vs E",pidParticle[pidIndex].Data()),
645  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
646  fhLam1E[pidIndex]->SetYTitle("#lambda_{1}^{2}");
647  fhLam1E[pidIndex]->SetXTitle("E (GeV)");
648  outputContainer->Add(fhLam1E[pidIndex]);
649 
650  fhDispE[pidIndex] = new TH2F (Form("h%sDispE",pidParticle[pidIndex].Data()),
651  Form("%s: dispersion^{2} vs E",pidParticle[pidIndex].Data()),
652  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
653  fhDispE[pidIndex]->SetYTitle("D^{2}");
654  fhDispE[pidIndex]->SetXTitle("E (GeV) ");
655  outputContainer->Add(fhDispE[pidIndex]);
656 
657  if(GetCalorimeter() == kEMCAL && GetFirstSMCoveredByTRD() >=0 )
658  {
659  fhLam0ETRD[pidIndex] = new TH2F (Form("h%sLam0ETRD",pidParticle[pidIndex].Data()),
660  Form("%s: #lambda_{0}^{2} vs E, EMCAL SM covered by TRD",pidParticle[pidIndex].Data()),
661  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
662  fhLam0ETRD[pidIndex]->SetYTitle("#lambda_{0}^{2}");
663  fhLam0ETRD[pidIndex]->SetXTitle("E (GeV)");
664  outputContainer->Add(fhLam0ETRD[pidIndex]);
665 
666  fhLam1ETRD[pidIndex] = new TH2F (Form("h%sLam1ETRD",pidParticle[pidIndex].Data()),
667  Form("%s: #lambda_{1}^{2} vs E, EMCAL SM covered by TRD",pidParticle[pidIndex].Data()),
668  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
669  fhLam1ETRD[pidIndex]->SetYTitle("#lambda_{1}^{2}");
670  fhLam1ETRD[pidIndex]->SetXTitle("E (GeV)");
671  outputContainer->Add(fhLam1ETRD[pidIndex]);
672 
673  fhDispETRD[pidIndex] = new TH2F (Form("h%sDispETRD",pidParticle[pidIndex].Data()),
674  Form("%s: dispersion^{2} vs E, EMCAL SM covered by TRD",pidParticle[pidIndex].Data()),
675  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
676  fhDispETRD[pidIndex]->SetYTitle("Dispersion^{2}");
677  fhDispETRD[pidIndex]->SetXTitle("E (GeV) ");
678  outputContainer->Add(fhDispETRD[pidIndex]);
679  }
680 
682  {
683  fhNCellsLam0LowE[pidIndex] = new TH2F (Form("h%sNCellsLam0LowE",pidParticle[pidIndex].Data()),
684  Form("%s: N_{cells} in cluster vs #lambda_{0}^{2}, E < 2 GeV",pidParticle[pidIndex].Data()),
685  nbins,nmin, nmax, ssbins,ssmin,ssmax);
686  fhNCellsLam0LowE[pidIndex]->SetXTitle("N_{cells}");
687  fhNCellsLam0LowE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
688  outputContainer->Add(fhNCellsLam0LowE[pidIndex]);
689 
690  fhNCellsLam0HighE[pidIndex] = new TH2F (Form("h%sNCellsLam0HighE",pidParticle[pidIndex].Data()),
691  Form("%s: N_{cells} in cluster vs #lambda_{0}^{2}, E > 2 GeV",pidParticle[pidIndex].Data()),
692  nbins,nmin, nmax, ssbins,ssmin,ssmax);
693  fhNCellsLam0HighE[pidIndex]->SetXTitle("N_{cells}");
694  fhNCellsLam0HighE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
695  outputContainer->Add(fhNCellsLam0HighE[pidIndex]);
696 
697 
698  fhEtaLam0LowE[pidIndex] = new TH2F (Form("h%sEtaLam0LowE",pidParticle[pidIndex].Data()),
699  Form("%s: #eta vs #lambda_{0}^{2}, E < 2 GeV",pidParticle[pidIndex].Data()),
700  netabins,etamin,etamax, ssbins,ssmin,ssmax);
701  fhEtaLam0LowE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
702  fhEtaLam0LowE[pidIndex]->SetXTitle("#eta");
703  outputContainer->Add(fhEtaLam0LowE[pidIndex]);
704 
705  fhPhiLam0LowE[pidIndex] = new TH2F (Form("h%sPhiLam0LowE",pidParticle[pidIndex].Data()),
706  Form("%s: #varphi vs #lambda_{0}^{2}, E < 2 GeV",pidParticle[pidIndex].Data()),
707  nphibins,phimin,phimax, ssbins,ssmin,ssmax);
708  fhPhiLam0LowE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
709  fhPhiLam0LowE[pidIndex]->SetXTitle("#varphi (rad)");
710  outputContainer->Add(fhPhiLam0LowE[pidIndex]);
711 
712  fhEtaLam0HighE[pidIndex] = new TH2F (Form("h%sEtaLam0HighE",pidParticle[pidIndex].Data()),
713  Form("%s: #eta vs #lambda_{0}^{2}, E > 2 GeV",pidParticle[pidIndex].Data()),
714  netabins,etamin,etamax, ssbins,ssmin,ssmax);
715  fhEtaLam0HighE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
716  fhEtaLam0HighE[pidIndex]->SetXTitle("#eta");
717  outputContainer->Add(fhEtaLam0HighE[pidIndex]);
718 
719  fhPhiLam0HighE[pidIndex] = new TH2F (Form("h%sPhiLam0HighE",pidParticle[pidIndex].Data()),
720  Form("%s: #varphi vs #lambda_{0}^{2}, E > 2 GeV",pidParticle[pidIndex].Data()),
721  nphibins,phimin,phimax, ssbins,ssmin,ssmax);
722  fhPhiLam0HighE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
723  fhPhiLam0HighE[pidIndex]->SetXTitle("#varphi (rad)");
724  outputContainer->Add(fhPhiLam0HighE[pidIndex]);
725 
726  if(GetCalorimeter() == kEMCAL)
727  {
728  fhDispEtaE[pidIndex] = new TH2F (Form("h%sDispEtaE",pidParticle[pidIndex].Data()),
729  Form("%s: #sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs E",pidParticle[pidIndex].Data()),
730  nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
731  fhDispEtaE[pidIndex]->SetXTitle("E (GeV)");
732  fhDispEtaE[pidIndex]->SetYTitle("#sigma^{2}_{#eta #eta}");
733  outputContainer->Add(fhDispEtaE[pidIndex]);
734 
735  fhDispPhiE[pidIndex] = new TH2F (Form("h%sDispPhiE",pidParticle[pidIndex].Data()),
736  Form("%s: #sigma^{2}_{#varphi #varphi} = #Sigma w_{i}(#varphi_{i} - <#varphi>)^{2} / #Sigma w_{i} vs E",pidParticle[pidIndex].Data()),
737  nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
738  fhDispPhiE[pidIndex]->SetXTitle("E (GeV)");
739  fhDispPhiE[pidIndex]->SetYTitle("#sigma^{2}_{#varphi #varphi}");
740  outputContainer->Add(fhDispPhiE[pidIndex]);
741 
742  fhSumEtaE[pidIndex] = new TH2F (Form("h%sSumEtaE",pidParticle[pidIndex].Data()),
743  Form("%s: #sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i})^{2} / #Sigma w_{i} - <#eta>^{2} vs E",pidParticle[pidIndex].Data()),
744  nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
745  fhSumEtaE[pidIndex]->SetXTitle("E (GeV)");
746  fhSumEtaE[pidIndex]->SetYTitle("#delta^{2}_{#eta #eta}");
747  outputContainer->Add(fhSumEtaE[pidIndex]);
748 
749  fhSumPhiE[pidIndex] = new TH2F (Form("h%sSumPhiE",pidParticle[pidIndex].Data()),
750  Form("%s: #sigma^{2}_{#varphi #varphi} = #Sigma w_{i}(#varphi_{i})^{2}/ #Sigma w_{i} - <#varphi>^{2} vs E",pidParticle[pidIndex].Data()),
751  nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
752  fhSumPhiE[pidIndex]->SetXTitle("E (GeV)");
753  fhSumPhiE[pidIndex]->SetYTitle("#delta^{2}_{#varphi #varphi}");
754  outputContainer->Add(fhSumPhiE[pidIndex]);
755 
756  fhSumEtaPhiE[pidIndex] = new TH2F (Form("h%sSumEtaPhiE",pidParticle[pidIndex].Data()),
757  Form("%s: #delta^{2}_{#eta #varphi} = #Sigma w_{i}(#varphi_{i} #eta_{i} ) / #Sigma w_{i} - <#varphi><#eta> vs E",pidParticle[pidIndex].Data()),
758  nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
759  fhSumEtaPhiE[pidIndex]->SetXTitle("E (GeV)");
760  fhSumEtaPhiE[pidIndex]->SetYTitle("#delta^{2}_{#eta #varphi}");
761  outputContainer->Add(fhSumEtaPhiE[pidIndex]);
762 
763  fhDispEtaPhiDiffE[pidIndex] = new TH2F (Form("h%sDispEtaPhiDiffE",pidParticle[pidIndex].Data()),
764  Form("%s: #sigma^{2}_{#varphi #varphi} - #sigma^{2}_{#eta #eta} vs E",pidParticle[pidIndex].Data()),
765  nptbins,ptmin,ptmax,200, -10,10);
766  fhDispEtaPhiDiffE[pidIndex]->SetXTitle("E (GeV)");
767  fhDispEtaPhiDiffE[pidIndex]->SetYTitle("#sigma^{2}_{#varphi #varphi}-#sigma^{2}_{#eta #eta}");
768  outputContainer->Add(fhDispEtaPhiDiffE[pidIndex]);
769 
770  fhSphericityE[pidIndex] = new TH2F (Form("h%sSphericityE",pidParticle[pidIndex].Data()),
771  Form("%s: (#sigma^{2}_{#varphi #varphi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#varphi #varphi}) vs E",pidParticle[pidIndex].Data()),
772  nptbins,ptmin,ptmax, 200, -1,1);
773  fhSphericityE[pidIndex]->SetXTitle("E (GeV)");
774  fhSphericityE[pidIndex]->SetYTitle("s = (#sigma^{2}_{#varphi #varphi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#varphi #varphi})");
775  outputContainer->Add(fhSphericityE[pidIndex]);
776 
777  Int_t bin[] = {0,2,4,6,10,1000};
778  for(Int_t i = 0; i < 5; i++)
779  {
780  fhDispEtaDispPhiEBin[pidIndex][i] = new TH2F (Form("h%sDispEtaDispPhi_EBin%d",pidParticle[pidIndex].Data(),i),
781  Form("%s: #sigma^{2}_{#varphi #varphi} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",pidParticle[pidIndex].Data(),bin[i],bin[i+1]),
782  ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
783  fhDispEtaDispPhiEBin[pidIndex][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
784  fhDispEtaDispPhiEBin[pidIndex][i]->SetYTitle("#sigma^{2}_{#varphi #varphi}");
785  outputContainer->Add(fhDispEtaDispPhiEBin[pidIndex][i]);
786  }
787  }
788  }
789  } // Shower shape
790 
791  if(IsDataMC())
792  {
794  {
795  for(Int_t i = 0; i < 6; i++)
796  {
797  fhMCELambda0[pidIndex][i] = new TH2F(Form("h%sELambda0_MC%s",pidParticle[pidIndex].Data(),pnamess[i].Data()),
798  Form("%s like cluster from %s : E vs #lambda_{0}^{2}",pidParticle[pidIndex].Data(),ptypess[i].Data()),
799  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
800  fhMCELambda0[pidIndex][i]->SetYTitle("#lambda_{0}^{2}");
801  fhMCELambda0[pidIndex][i]->SetXTitle("E (GeV)");
802  outputContainer->Add(fhMCELambda0[pidIndex][i]) ;
803 
805  {
806  fhMCEDispEta[pidIndex][i] = new TH2F (Form("h%sEDispEtaE_MC%s",pidParticle[pidIndex].Data(),pnamess[i].Data()),
807  Form("cluster from %s : %s like, #sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs E",ptypess[i].Data(),pidParticle[pidIndex].Data()),
808  nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
809  fhMCEDispEta[pidIndex][i]->SetXTitle("E (GeV)");
810  fhMCEDispEta[pidIndex][i]->SetYTitle("#sigma^{2}_{#eta #eta}");
811  outputContainer->Add(fhMCEDispEta[pidIndex][i]);
812 
813  fhMCEDispPhi[pidIndex][i] = new TH2F (Form("h%sEDispPhiE_MC%s",pidParticle[pidIndex].Data(),pnamess[i].Data()),
814  Form("cluster from %s : %s like, #sigma^{2}_{#varphi #varphi} = #Sigma w_{i}(#varphi_{i} - <#varphi>)^{2} / #Sigma w_{i} vs E",ptypess[i].Data(),pidParticle[pidIndex].Data()),
815  nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
816  fhMCEDispPhi[pidIndex][i]->SetXTitle("E (GeV)");
817  fhMCEDispPhi[pidIndex][i]->SetYTitle("#sigma^{2}_{#varphi #varphi}");
818  outputContainer->Add(fhMCEDispPhi[pidIndex][i]);
819 
820  fhMCESumEtaPhi[pidIndex][i] = new TH2F (Form("h%sESumEtaPhiE_MC%s",pidParticle[pidIndex].Data(),pnamess[i].Data()),
821  Form("cluster from %s : %s like, #delta^{2}_{#eta #varphi} = #Sigma w_{i}(#varphi_{i} #eta_{i} ) / #Sigma w_{i} - <#varphi><#eta> vs E",ptypess[i].Data(),pidParticle[pidIndex].Data()),
822  nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
823  fhMCESumEtaPhi[pidIndex][i]->SetXTitle("E (GeV)");
824  fhMCESumEtaPhi[pidIndex][i]->SetYTitle("#delta^{2}_{#eta #varphi}");
825  outputContainer->Add(fhMCESumEtaPhi[pidIndex][i]);
826 
827  fhMCEDispEtaPhiDiff[pidIndex][i] = new TH2F (Form("h%sEDispEtaPhiDiffE_MC%s",pidParticle[pidIndex].Data(),pnamess[i].Data()),
828  Form("cluster from %s : %s like, #sigma^{2}_{#varphi #varphi} - #sigma^{2}_{#eta #eta} vs E",ptypess[i].Data(),pidParticle[pidIndex].Data()),
829  nptbins,ptmin,ptmax,200,-10,10);
830  fhMCEDispEtaPhiDiff[pidIndex][i]->SetXTitle("E (GeV)");
831  fhMCEDispEtaPhiDiff[pidIndex][i]->SetYTitle("#sigma^{2}_{#varphi #varphi}-#sigma^{2}_{#eta #eta}");
832  outputContainer->Add(fhMCEDispEtaPhiDiff[pidIndex][i]);
833 
834  fhMCESphericity[pidIndex][i] = new TH2F (Form("h%sESphericity_MC%s",pidParticle[pidIndex].Data(),pnamess[i].Data()),
835  Form("cluster from %s : %s like, (#sigma^{2}_{#varphi #varphi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#varphi #varphi}) vs E",ptypess[i].Data(),pidParticle[pidIndex].Data()),
836  nptbins,ptmin,ptmax, 200,-1,1);
837  fhMCESphericity[pidIndex][i]->SetXTitle("E (GeV)");
838  fhMCESphericity[pidIndex][i]->SetYTitle("s = (#sigma^{2}_{#varphi #varphi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#varphi #varphi})");
839  outputContainer->Add(fhMCESphericity[pidIndex][i]);
840  }
841 
842  }// loop
843  }
844  }
845 
846  //if(IsCaloPIDOn() && pidIndex > 0) continue;
847 
848  fhNCellsE[pidIndex] = new TH2F (Form("h%sNCellsE",pidParticle[pidIndex].Data()),
849  Form("N cells in %s cluster vs E ",pidParticle[pidIndex].Data()),
850  nptbins,ptmin,ptmax, nbins,nmin,nmax);
851  fhNCellsE[pidIndex]->SetXTitle("E (GeV)");
852  fhNCellsE[pidIndex]->SetYTitle("# of cells in cluster");
853  outputContainer->Add(fhNCellsE[pidIndex]);
854 
855  fhNLME[pidIndex] = new TH2F (Form("h%sNLME",pidParticle[pidIndex].Data()),
856  Form("NLM in %s cluster vs E ",pidParticle[pidIndex].Data()),
857  nptbins,ptmin,ptmax, 10,0,10);
858  fhNLME[pidIndex]->SetXTitle("E (GeV)");
859  fhNLME[pidIndex]->SetYTitle("# of cells in cluster");
860  outputContainer->Add(fhNLME[pidIndex]);
861 
862  fhTimeE[pidIndex] = new TH2F(Form("h%sTimeE",pidParticle[pidIndex].Data()),
863  Form("Time in %s cluster vs E ",pidParticle[pidIndex].Data())
864  ,nptbins,ptmin,ptmax, tbins,tmin,tmax);
865  fhTimeE[pidIndex]->SetXTitle("E (GeV)");
866  fhTimeE[pidIndex]->SetYTitle(" t (ns)");
867  outputContainer->Add(fhTimeE[pidIndex]);
868 
869  fhMaxCellDiffClusterE[pidIndex] = new TH2F (Form("h%sMaxCellDiffClusterE",pidParticle[pidIndex].Data()),
870  Form("%s: energy vs difference of cluster energy - max cell energy / cluster energy, good clusters",pidParticle[pidIndex].Data()),
871  nptbins,ptmin,ptmax, 500,0,1.);
872  fhMaxCellDiffClusterE[pidIndex]->SetXTitle("E_{cluster} (GeV) ");
873  fhMaxCellDiffClusterE[pidIndex]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
874  outputContainer->Add(fhMaxCellDiffClusterE[pidIndex]);
875 
876  fhE[pidIndex] = new TH1F(Form("h%sE",pidParticle[pidIndex].Data()),
877  Form("Number of %s over calorimeter vs energy",pidParticle[pidIndex].Data()),
878  nptbins,ptmin,ptmax);
879  fhE[pidIndex]->SetYTitle("N");
880  fhE[pidIndex]->SetXTitle("E_{#gamma}(GeV)");
881  outputContainer->Add(fhE[pidIndex]) ;
882 
883  fhPt[pidIndex] = new TH1F(Form("h%sPtElectron",pidParticle[pidIndex].Data()),
884  Form("Number of %s over calorimeter vs p_{T}",pidParticle[pidIndex].Data()),
885  nptbins,ptmin,ptmax);
886  fhPt[pidIndex]->SetYTitle("N");
887  fhPt[pidIndex]->SetXTitle("p_{T #gamma}(GeV/c)");
888  outputContainer->Add(fhPt[pidIndex]) ;
889 
890  fhPhi[pidIndex] = new TH2F(Form("h%sPhiElectron",pidParticle[pidIndex].Data()),
891  Form("%s: #varphi vs p_{T}",pidParticle[pidIndex].Data()),
892  nptbins,ptmin,ptmax,nphibins,phimin,phimax);
893  fhPhi[pidIndex]->SetYTitle("#varphi (rad)");
894  fhPhi[pidIndex]->SetXTitle("p_{T #gamma} (GeV/c)");
895  outputContainer->Add(fhPhi[pidIndex]) ;
896 
897  fhEta[pidIndex] = new TH2F(Form("h%sEta",pidParticle[pidIndex].Data()),
898  Form("%s: #eta vs p_{T}",pidParticle[pidIndex].Data()),
899  nptbins,ptmin,ptmax,netabins,etamin,etamax);
900  fhEta[pidIndex]->SetYTitle("#eta");
901  fhEta[pidIndex]->SetXTitle("p_{T #gamma} (GeV/c)");
902  outputContainer->Add(fhEta[pidIndex]) ;
903 
904  fhEtaPhi[pidIndex] = new TH2F(Form("h%sEtaPhi",pidParticle[pidIndex].Data()),
905  Form("%s: #eta vs #varphi",pidParticle[pidIndex].Data()),
906  netabins,etamin,etamax,nphibins,phimin,phimax);
907  fhEtaPhi[pidIndex]->SetYTitle("#varphi (rad)");
908  fhEtaPhi[pidIndex]->SetXTitle("#eta");
909  outputContainer->Add(fhEtaPhi[pidIndex]) ;
910  if(GetMinPt() < 0.5)
911  {
912  fhEtaPhi05[pidIndex] = new TH2F(Form("h%sEtaPhi05",pidParticle[pidIndex].Data()),
913  Form("%s: #eta vs #varphi, E > 0.5",pidParticle[pidIndex].Data()),
914  netabins,etamin,etamax,nphibins,phimin,phimax);
915  fhEtaPhi05[pidIndex]->SetYTitle("#varphi (rad)");
916  fhEtaPhi05[pidIndex]->SetXTitle("#eta");
917  outputContainer->Add(fhEtaPhi05[pidIndex]) ;
918  }
919 
920 
921  if(IsDataMC())
922  {
923  for(Int_t i = 0; i < fNOriginHistograms; i++)
924  {
925  fhMCE[pidIndex][i] = new TH1F(Form("h%sE_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
926  Form("%s like cluster from %s : E ",pidParticle[pidIndex].Data(),ptype[i].Data()),
927  nptbins,ptmin,ptmax);
928  fhMCE[pidIndex][i]->SetXTitle("E (GeV)");
929  outputContainer->Add(fhMCE[pidIndex][i]) ;
930 
931  fhMCPt[pidIndex][i] = new TH1F(Form("h%sPt_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
932  Form("%s like cluster from %s : p_{T} ",pidParticle[pidIndex].Data(),ptype[i].Data()),
933  nptbins,ptmin,ptmax);
934  fhMCPt[pidIndex][i]->SetXTitle("p_{T} (GeV/c)");
935  outputContainer->Add(fhMCPt[pidIndex][i]) ;
936 
937  fhMCEta[pidIndex][i] = new TH2F(Form("h%sEta_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
938  Form("%s like cluster from %s : #eta ",pidParticle[pidIndex].Data(),ptype[i].Data()),
939  nptbins,ptmin,ptmax,netabins,etamin,etamax);
940  fhMCEta[pidIndex][i]->SetYTitle("#eta");
941  fhMCEta[pidIndex][i]->SetXTitle("E (GeV)");
942  outputContainer->Add(fhMCEta[pidIndex][i]) ;
943 
944  fhMCPhi[pidIndex][i] = new TH2F(Form("h%sPhi_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
945  Form("%s like cluster from %s : #varphi ",pidParticle[pidIndex].Data(),ptype[i].Data()),
946  nptbins,ptmin,ptmax,nphibins,phimin,phimax);
947  fhMCPhi[pidIndex][i]->SetYTitle("#varphi (rad)");
948  fhMCPhi[pidIndex][i]->SetXTitle("E (GeV)");
949  outputContainer->Add(fhMCPhi[pidIndex][i]) ;
950 
951 
952  fhMCDeltaE[pidIndex][i] = new TH2F (Form("h%sDeltaE_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
953  Form("%s like MC - Reco E from %s",pidParticle[pidIndex].Data(),pname[i].Data()),
954  nptbins,ptmin,ptmax, 200,-50,50);
955  fhMCDeltaE[pidIndex][i]->SetXTitle("#Delta E (GeV)");
956  outputContainer->Add(fhMCDeltaE[pidIndex][i]);
957 
958  fhMC2E[pidIndex][i] = new TH2F (Form("h%s2E_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
959  Form("%s like E distribution, reconstructed vs generated from %s",pidParticle[pidIndex].Data(),pname[i].Data()),
960  nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
961  fhMC2E[pidIndex][i]->SetXTitle("E_{rec} (GeV)");
962  fhMC2E[pidIndex][i]->SetYTitle("E_{gen} (GeV)");
963  outputContainer->Add(fhMC2E[pidIndex][i]);
964 
965  }
966  } // MC
967  }// pid Index
968 
969 
971  {
972  if(IsDataMC())
973  {
974  if(!GetReader()->IsEmbeddedClusterSelectionOn())
975  {
976  fhMCElectronELambda0NoOverlap = new TH2F("hELambda0_MCElectron_NoOverlap",
977  "cluster from Electron : E vs #lambda_{0}^{2}",
978  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
979  fhMCElectronELambda0NoOverlap->SetYTitle("#lambda_{0}^{2}");
980  fhMCElectronELambda0NoOverlap->SetXTitle("E (GeV)");
981  outputContainer->Add(fhMCElectronELambda0NoOverlap) ;
982 
983  fhMCElectronELambda0TwoOverlap = new TH2F("hELambda0_MCElectron_TwoOverlap",
984  "cluster from Electron : E vs #lambda_{0}^{2}",
985  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
986  fhMCElectronELambda0TwoOverlap->SetYTitle("#lambda_{0}^{2}");
987  fhMCElectronELambda0TwoOverlap->SetXTitle("E (GeV)");
988  outputContainer->Add(fhMCElectronELambda0TwoOverlap) ;
989 
990  fhMCElectronELambda0NOverlap = new TH2F("hELambda0_MCElectron_NOverlap",
991  "cluster from Electron : E vs #lambda_{0}^{2}",
992  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
993  fhMCElectronELambda0NOverlap->SetYTitle("#lambda_{0}^{2}");
994  fhMCElectronELambda0NOverlap->SetXTitle("E (GeV)");
995  outputContainer->Add(fhMCElectronELambda0NOverlap) ;
996  } // No embedding
997 
998  // Fill histograms to check shape of embedded clusters
999  if(GetReader()->IsEmbeddedClusterSelectionOn())
1000  {
1001  fhEmbeddedSignalFractionEnergy = new TH2F("hEmbeddedSignal_FractionEnergy",
1002  "Energy Fraction of embedded signal versus cluster energy",
1003  nptbins,ptmin,ptmax,100,0.,1.);
1004  fhEmbeddedSignalFractionEnergy->SetYTitle("Fraction");
1005  fhEmbeddedSignalFractionEnergy->SetXTitle("E (GeV)");
1006  outputContainer->Add(fhEmbeddedSignalFractionEnergy) ;
1007 
1008  fhEmbedElectronELambda0FullSignal = new TH2F("hELambda0_EmbedElectron_FullSignal",
1009  "cluster from Electron embedded with more than 90% energy in cluster : E vs #lambda_{0}^{2}",
1010  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1011  fhEmbedElectronELambda0FullSignal->SetYTitle("#lambda_{0}^{2}");
1012  fhEmbedElectronELambda0FullSignal->SetXTitle("E (GeV)");
1013  outputContainer->Add(fhEmbedElectronELambda0FullSignal) ;
1014 
1015  fhEmbedElectronELambda0MostlySignal = new TH2F("hELambda0_EmbedElectron_MostlySignal",
1016  "cluster from Electron embedded with 50% to 90% energy in cluster : E vs #lambda_{0}^{2}",
1017  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1018  fhEmbedElectronELambda0MostlySignal->SetYTitle("#lambda_{0}^{2}");
1019  fhEmbedElectronELambda0MostlySignal->SetXTitle("E (GeV)");
1020  outputContainer->Add(fhEmbedElectronELambda0MostlySignal) ;
1021 
1022  fhEmbedElectronELambda0MostlyBkg = new TH2F("hELambda0_EmbedElectron_MostlyBkg",
1023  "cluster from Electron embedded with 10% to 50% energy in cluster : E vs #lambda_{0}^{2}",
1024  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1025  fhEmbedElectronELambda0MostlyBkg->SetYTitle("#lambda_{0}^{2}");
1026  fhEmbedElectronELambda0MostlyBkg->SetXTitle("E (GeV)");
1027  outputContainer->Add(fhEmbedElectronELambda0MostlyBkg) ;
1028 
1029  fhEmbedElectronELambda0FullBkg = new TH2F("hELambda0_EmbedElectron_FullBkg",
1030  "cluster from Electronm embedded with 0% to 10% energy in cluster : E vs #lambda_{0}^{2}",
1031  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1032  fhEmbedElectronELambda0FullBkg->SetYTitle("#lambda_{0}^{2}");
1033  fhEmbedElectronELambda0FullBkg->SetXTitle("E (GeV)");
1034  outputContainer->Add(fhEmbedElectronELambda0FullBkg) ;
1035  } // Embedded histograms
1036  } // Histos with MC
1037  } // Fill SS MC histograms
1038 
1039  return outputContainer ;
1040 }
1041 
1042 //_________________________
1044 //_________________________
1046 {
1047  if ( GetCalorimeter() == kPHOS && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD() )
1048  AliFatal("STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!");
1049  else if ( GetCalorimeter() == kEMCAL && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD() )
1050  AliFatal("STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!");
1051 }
1052 
1053 //___________________________________
1055 //___________________________________
1057 {
1058  AddToHistogramsName("AnaElectron_");
1059 
1060  fMinDist = 2.;
1061  fMinDist2 = 4.;
1062  fMinDist3 = 5.;
1063 
1064  fTimeCutMin = -1;
1065  fTimeCutMax = 9999999;
1066  fNCellsCut = 0;
1067 
1068  fdEdxMin = 76.; // for LHC11a, but for LHC11c pass1 56.
1069  fdEdxMax = 85.; // for LHC11a, but for LHC11c pass1 64.
1070 
1071  fEOverPMin = 0.8; // for LHC11a, but for LHC11c pass1 0.9
1072  fEOverPMax = 1.2; // for LHC11a and LHC11c pass1
1073 }
1074 
1075 //_________________________________________
1077 //_________________________________________
1079 {
1080  // Get the vertex
1081  Double_t v[3] = {0,0,0}; //vertex ;
1082  GetReader()->GetVertex(v);
1083 
1084  // Select the Calorimeter of the photon
1085  TObjArray * pl = 0x0;
1086  if (GetCalorimeter() == kPHOS ) pl = GetPHOSClusters ();
1087  else if (GetCalorimeter() == kEMCAL) pl = GetEMCALClusters();
1088 
1089  if(!pl)
1090  {
1091  AliWarning(Form("TObjArray with %s clusters is NULL!",GetCalorimeterString().Data()));
1092  return;
1093  }
1094 
1095  //Init arrays, variables, get number of clusters
1096  Int_t nCaloClusters = pl->GetEntriesFast();
1097  //List to be used in conversion analysis, to tag the cluster as candidate for conversion
1098 
1099  AliDebug(1,Form("Input %s cluster entries %d", GetCalorimeterString().Data(), nCaloClusters));
1100 
1101  //----------------------------------------------------
1102  // Fill AOD with PHOS/EMCAL AliCaloTrackParticle objects
1103  //----------------------------------------------------
1104  // Loop on clusters
1105  for(Int_t icalo = 0; icalo < nCaloClusters; icalo++)
1106  {
1107  AliVCluster * calo = (AliVCluster*) (pl->At(icalo));
1108  //printf("calo %d, %f\n",icalo,calo->E());
1109 
1110  //Get the index where the cluster comes, to retrieve the corresponding vertex
1111  Int_t evtIndex = 0 ;
1112  if (GetMixedEvent())
1113  {
1114  evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
1115  //Get the vertex and check it is not too large in z
1116  if(TMath::Abs(GetVertex(evtIndex)[2])> GetZvertexCut()) continue;
1117  }
1118 
1119  //Cluster selection, not charged, with photon id and in fiducial cut
1121  {
1122  calo->GetMomentum(fMomentum,GetVertex(evtIndex)) ;
1123  }//Assume that come from vertex in straight line
1124  else
1125  {
1126  Double_t vertex[]={0,0,0};
1127  calo->GetMomentum(fMomentum,vertex) ;
1128  }
1129 
1130  //--------------------------------------
1131  // Cluster selection
1132  //--------------------------------------
1133  AliVCaloCells* cells = 0;
1134  if(GetCalorimeter() == kEMCAL) cells = GetEMCALCells();
1135  else cells = GetPHOSCells();
1136 
1137  Int_t nMaxima = GetCaloUtils()->GetNumberOfLocalMaxima(calo, cells); // NLM
1138  if(!ClusterSelected(calo,nMaxima)) continue;
1139 
1140  //-------------------------------------
1141  // PID selection via dE/dx
1142  //-------------------------------------
1143 
1144  AliVTrack *track = GetCaloUtils()->GetMatchedTrack(calo, GetReader()->GetInputEvent());
1145 
1146  if(!track)
1147  {
1148  AliWarning("Null track");
1149  continue;
1150  }
1151 
1152  //printf("track dedx %f, p %f, cluster E %f\n",track->GetTPCsignal(),track->P(),calo->E());
1153  Float_t dEdx = track->GetTPCsignal();
1154  Float_t eOverp = calo->E()/track->P();
1155 
1156  fhdEdxvsE->Fill(calo ->E(), dEdx, GetEventWeight());
1157  fhdEdxvsP->Fill(track->P(), dEdx, GetEventWeight());
1158 
1159  if( eOverp < fEOverPMax && eOverp > fEOverPMin)
1160  {
1161  fhdEdxvsECutEOverP ->Fill(calo ->E(), dEdx, GetEventWeight());
1162  fhdEdxvsPCutEOverP ->Fill(track->P(), dEdx, GetEventWeight());
1163  }
1164 
1165  // Apply a mild cut on the cluster SS and check the value of dEdX and EOverP
1166  Float_t m02 = calo->GetM02();
1167  if(m02 > 0.1 && m02 < 0.4)
1168  {
1169  fhdEdxvsECutM02 ->Fill(calo ->E(), dEdx , GetEventWeight());
1170  fhdEdxvsPCutM02 ->Fill(track->P(), dEdx , GetEventWeight());
1171  fhEOverPvsECutM02->Fill(calo ->E(), eOverp, GetEventWeight());
1172  fhEOverPvsPCutM02->Fill(track->P(), eOverp, GetEventWeight());
1173  }
1174 
1176 
1177  if( dEdx < fdEdxMax && dEdx > fdEdxMin)
1178  {
1179  fhEOverPvsE->Fill(calo ->E(), eOverp, GetEventWeight());
1180  fhEOverPvsP->Fill(track->P(), eOverp, GetEventWeight());
1181 
1182  if(m02 > 0.1 && m02 < 0.4)
1183  {
1184  fhEOverPvsECutM02CutdEdx->Fill(calo ->E(), eOverp, GetEventWeight());
1185  fhEOverPvsPCutM02CutdEdx->Fill(track->P(), eOverp, GetEventWeight());
1186  }
1187 
1188  if( eOverp < fEOverPMax && eOverp > fEOverPMin)
1189  {
1190  pid = AliCaloPID::kElectron;
1191  } // E/p
1192 
1193  }// dE/dx
1194 
1195  Int_t pidIndex = 0;// Electron
1196  if(pid == AliCaloPID::kChargedHadron) pidIndex = 1;
1197 
1198  //--------------------------------------------------------------------------------------
1199  // Play with the MC stack if available
1200  //--------------------------------------------------------------------------------------
1201 
1202  //Check origin of the candidates
1203  Int_t tag = -1 ;
1204  if(IsDataMC())
1205  {
1206  tag = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(), calo->GetNLabels(), GetMC());
1207 
1208  AliDebug(1,Form("Origin of candidate, bit map %d",tag));
1209 
1210  if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) && fhMCE[pidIndex][kmcPhoton])
1211  {
1212  fhMCdEdxvsE [kmcPhoton]->Fill(calo ->E(), dEdx , GetEventWeight());
1213  fhMCdEdxvsP [kmcPhoton]->Fill(track->P(), dEdx , GetEventWeight());
1214  fhMCEOverPvsE[kmcPhoton]->Fill(calo ->E(), eOverp, GetEventWeight());
1215  fhMCEOverPvsP[kmcPhoton]->Fill(track->P(), eOverp, GetEventWeight());
1216 
1218  {
1219  fhMCdEdxvsE [kmcConversion]->Fill(calo ->E(), dEdx , GetEventWeight());
1220  fhMCdEdxvsP [kmcConversion]->Fill(track->P(), dEdx , GetEventWeight());
1221  fhMCEOverPvsE[kmcConversion]->Fill(calo ->E(), eOverp, GetEventWeight());
1222  fhMCEOverPvsP[kmcConversion]->Fill(track->P(), eOverp, GetEventWeight());
1223  }
1226  {
1227  fhMCdEdxvsE [kmcPi0Decay]->Fill(calo ->E(), dEdx , GetEventWeight());
1228  fhMCdEdxvsP [kmcPi0Decay]->Fill(track->P(), dEdx , GetEventWeight());
1229  fhMCEOverPvsE[kmcPi0Decay]->Fill(calo ->E(), eOverp, GetEventWeight());
1230  fhMCEOverPvsP[kmcPi0Decay]->Fill(track->P(), eOverp, GetEventWeight());
1231  }
1232  else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE [pidIndex][kmcPi0])
1233  {
1234  fhMCdEdxvsE [kmcPi0]->Fill(calo ->E(), dEdx , GetEventWeight());
1235  fhMCdEdxvsP [kmcPi0]->Fill(track->P(), dEdx , GetEventWeight());
1236  fhMCEOverPvsE[kmcPi0]->Fill(calo ->E(), eOverp, GetEventWeight());
1237  fhMCEOverPvsP[kmcPi0]->Fill(track->P(), eOverp, GetEventWeight());
1238  }
1240  {
1241  fhMCdEdxvsE [kmcEta]->Fill(calo ->E(), dEdx , GetEventWeight());
1242  fhMCdEdxvsP [kmcEta]->Fill(track->P(), dEdx , GetEventWeight());
1243  fhMCEOverPvsE[kmcEta]->Fill(calo ->E(), eOverp, GetEventWeight());
1244  fhMCEOverPvsP[kmcEta]->Fill(track->P(), eOverp, GetEventWeight());
1245  }
1246  else if( fhMCE[pidIndex][kmcOtherDecay] )
1247  {
1248  fhMCdEdxvsE [kmcOtherDecay]->Fill(calo ->E(), dEdx , GetEventWeight());
1249  fhMCdEdxvsP [kmcOtherDecay]->Fill(track->P(), dEdx , GetEventWeight());
1250  fhMCEOverPvsE[kmcOtherDecay]->Fill(calo ->E(), eOverp, GetEventWeight());
1251  fhMCEOverPvsP[kmcOtherDecay]->Fill(track->P(), eOverp, GetEventWeight());
1252  }
1253  }
1255  {
1256  fhMCdEdxvsE [kmcAntiNeutron]->Fill(calo ->E(), dEdx , GetEventWeight());
1257  fhMCdEdxvsP [kmcAntiNeutron]->Fill(track->P(), dEdx , GetEventWeight());
1258  fhMCEOverPvsE[kmcAntiNeutron]->Fill(calo ->E(), eOverp, GetEventWeight());
1259  fhMCEOverPvsP[kmcAntiNeutron]->Fill(track->P(), eOverp, GetEventWeight());
1260  }
1262  {
1263  fhMCdEdxvsE [kmcAntiProton]->Fill(calo ->E(), dEdx , GetEventWeight());
1264  fhMCdEdxvsP [kmcAntiProton]->Fill(track->P(), dEdx , GetEventWeight());
1265  fhMCEOverPvsE[kmcAntiProton]->Fill(calo ->E(), eOverp, GetEventWeight());
1266  fhMCEOverPvsP[kmcAntiProton]->Fill(track->P(), eOverp, GetEventWeight());
1267  }
1269  {
1270  fhMCdEdxvsE [kmcElectron]->Fill(calo ->E(), dEdx , GetEventWeight());
1271  fhMCdEdxvsP [kmcElectron]->Fill(track->P(), dEdx , GetEventWeight());
1272  fhMCEOverPvsE[kmcElectron]->Fill(calo ->E(), eOverp, GetEventWeight());
1273  fhMCEOverPvsP[kmcElectron]->Fill(track->P(), eOverp, GetEventWeight());
1274  }
1275  else if( fhMCE[pidIndex][kmcOther])
1276  {
1277  fhMCdEdxvsE [kmcOther]->Fill(calo ->E(), dEdx , GetEventWeight());
1278  fhMCdEdxvsP [kmcOther]->Fill(track->P(), dEdx , GetEventWeight());
1279  fhMCEOverPvsE[kmcOther]->Fill(calo ->E(), eOverp, GetEventWeight());
1280  fhMCEOverPvsP[kmcOther]->Fill(track->P(), eOverp, GetEventWeight());
1281  }
1282  }// set MC tag and fill Histograms with MC
1283 
1284  //---------------------------------
1285  //Fill some shower shape histograms
1286  //---------------------------------
1287 
1288  FillShowerShapeHistograms(calo,tag,pid);
1289 
1290  if(pid == AliCaloPID::kElectron)
1291  WeightHistograms(calo);
1292 
1293  //-----------------------------------------
1294  // PID Shower Shape selection or bit setting
1295  //-----------------------------------------
1296 
1297  // Data, PID check on
1298  if(IsCaloPIDOn())
1299  {
1300  // Get most probable PID, 2 options check bayesian PID weights or redo PID
1301  // By default, redo PID
1302 
1303  if(GetCaloPID()->GetIdentifiedParticleType(calo)!=AliCaloPID::kPhoton)
1304  {
1306  continue;
1307 
1308  if(fAODParticle == 0 )
1310  }
1311 
1312  AliDebug(1,Form("PDG of identified particle %d",pid));
1313  }
1314 
1315  AliDebug(1,Form("Photon selection cuts passed: pT %3.2f, pdg %d",fMomentum.Pt(),pid));
1316 
1317  Float_t maxCellFraction = 0;
1318  Int_t absID = GetCaloUtils()->GetMaxEnergyCell(cells, calo,maxCellFraction);
1319  if ( absID >= 0 )fhMaxCellDiffClusterE[pidIndex]->Fill(fMomentum.E(), maxCellFraction, GetEventWeight());
1320 
1321  fhNCellsE[pidIndex] ->Fill(fMomentum.E(), calo->GetNCells() , GetEventWeight());
1322  fhNLME [pidIndex] ->Fill(fMomentum.E(), nMaxima , GetEventWeight());
1323  fhTimeE [pidIndex] ->Fill(fMomentum.E(), calo->GetTOF()*1.e9, GetEventWeight());
1324 
1325  //----------------------------
1326  // Create AOD for analysis
1327  //----------------------------
1328 
1329  //Add AOD with electron/hadron object to aod branch
1330  if ( pid == fAODParticle || fAODParticle == 0 )
1331  {
1333 
1334  //...............................................
1335  //Set the indeces of the original caloclusters (MC, ID), and calorimeter
1336  Int_t label = calo->GetLabel();
1337  aodpart.SetLabel(label);
1338  aodpart.SetCaloLabel (calo->GetID(),-1);
1339  aodpart.SetTrackLabel(GetReader()->GetTrackID(track),-1); // needed instead of track->GetID() since AOD needs some manipulations
1340 
1341  aodpart.SetDetectorTag(GetCalorimeter());
1342  //printf("Index %d, Id %d, iaod %d\n",icalo, calo->GetID(),GetOutputAODBranch()->GetEntriesFast());
1343 
1344  aodpart.SetM02(calo->GetM02());
1345  aodpart.SetM20(calo->GetM20());
1346  aodpart.SetNLM(nMaxima);
1347  aodpart.SetTime(calo->GetTOF()*1e9);
1348  aodpart.SetNCells(calo->GetNCells());
1349  Int_t nSM = GetModuleNumber(calo);
1350  aodpart.SetSModNumber(nSM);
1351 
1352  //...............................................
1353  //Set bad channel distance bit
1354  Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
1355  if (distBad > fMinDist3) aodpart.SetDistToBad(2) ;
1356  else if(distBad > fMinDist2) aodpart.SetDistToBad(1) ;
1357  else aodpart.SetDistToBad(0) ;
1358  //printf("DistBad %f Bit %d\n",distBad, aodpart.DistToBad());
1359 
1360  // MC tag
1361  aodpart.SetTag(tag);
1362 
1363  // PID tag
1364  aodpart.SetIdentifiedParticleType(pid);
1365 
1366  AddAODParticle(aodpart);
1367  }
1368 
1369  }//loop
1370 
1371  AliDebug(1,Form("End fill AODs, with %d entries",GetOutputAODBranch()->GetEntriesFast()));
1372 
1373 }
1374 
1375 //________________________________________________
1377 //________________________________________________
1379 {
1380  // Access MC information in stack if requested, check that it exists.
1381 
1382  AliVParticle * primary = 0x0;
1383 
1384  if( IsDataMC() && !GetMC() )
1385  {
1386  AliFatal("MCEvent not available! STOP");
1387  return;
1388  }
1389 
1390  // Get vertex
1391  Double_t v[3] = {0,0,0}; //vertex ;
1392  GetReader()->GetVertex(v);
1393  //fhVertex->Fill(v[0], v[1], v[2], GetEventWeight());
1394  if(TMath::Abs(v[2]) > GetZvertexCut()) return ; // done elsewhere for Single Event analysis, but there for mixed event
1395 
1396  //----------------------------------
1397  //Loop on stored AOD photons
1398  Int_t naod = GetOutputAODBranch()->GetEntriesFast();
1399  AliDebug(1,Form("AOD branch entries %d", naod));
1400 
1401  for(Int_t iaod = 0; iaod < naod ; iaod++)
1402  {
1405 
1406  Int_t pidIndex = 0;// Electron
1407  if (pdg == AliCaloPID::kElectron) pidIndex = 0;
1408  else if(pdg == AliCaloPID::kChargedHadron) pidIndex = 1;
1409  else continue ;
1410 
1411  if(((Int_t) ph->GetDetectorTag()) != GetCalorimeter()) continue;
1412 
1413  AliDebug(1,Form("ID Electron: pt %f, phi %f, eta %f", ph->Pt(),ph->Phi(),ph->Eta())) ;
1414 
1415  //................................
1416  //Fill photon histograms
1417  Float_t ptcluster = ph->Pt();
1418  Float_t phicluster = ph->Phi();
1419  Float_t etacluster = ph->Eta();
1420  Float_t ecluster = ph->E();
1421 
1422  fhE[pidIndex] ->Fill(ecluster, GetEventWeight());
1423  fhPt[pidIndex] ->Fill(ptcluster, GetEventWeight());
1424 
1425  fhPhi[pidIndex] ->Fill(ptcluster, phicluster, GetEventWeight());
1426  fhEta[pidIndex] ->Fill(ptcluster, etacluster, GetEventWeight());
1427 
1428  if (ecluster > 0.5) fhEtaPhi [pidIndex]->Fill(etacluster, phicluster, GetEventWeight());
1429  else if(GetMinPt() < 0.5) fhEtaPhi05[pidIndex]->Fill(etacluster, phicluster, GetEventWeight());
1430 
1431  //.......................................
1432  //Play with the MC data if available
1433  if(IsDataMC())
1434  {
1435  //....................................................................
1436  // Access MC information in stack if requested, check that it exists.
1437  Int_t label =ph->GetLabel();
1438  if(label < 0)
1439  {
1440  AliDebug(1,Form("*** bad label ***: label %d", label));
1441  continue;
1442  }
1443 
1444  Int_t nprim = GetMC()->GetNumberOfTracks();
1445  if ( label >= nprim )
1446  {
1447  AliDebug(1,Form("*** large label ***: label %d, n tracks %d", label, nprim));
1448  continue ;
1449  }
1450 
1451  Float_t eprim = 0;
1452  //Float_t ptprim = 0;
1453  primary = GetMC()->GetTrack(label);
1454 
1455  if(!primary)
1456  {
1457  AliWarning(Form("*** no primary ***: label %d", label));
1458  continue;
1459  }
1460 
1461  eprim = primary->E();
1462  //ptprim = aodprimary->Pt();
1463 
1464  Int_t tag =ph->GetTag();
1465 
1466  if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) && fhMCE[pidIndex][kmcPhoton])
1467  {
1468  fhMCE [pidIndex][kmcPhoton] ->Fill(ecluster , GetEventWeight());
1469  fhMCPt [pidIndex][kmcPhoton] ->Fill(ptcluster, GetEventWeight());
1470 
1471  fhMCPhi[pidIndex][kmcPhoton] ->Fill(ecluster, phicluster, GetEventWeight());
1472  fhMCEta[pidIndex][kmcPhoton] ->Fill(ecluster, etacluster, GetEventWeight());
1473 
1474  fhMC2E [pidIndex][kmcPhoton] ->Fill(ecluster, eprim , GetEventWeight());
1475  fhMCDeltaE[pidIndex][kmcPhoton] ->Fill(ecluster, eprim-ecluster, GetEventWeight());
1476 
1477  if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) && fhMCE[pidIndex][kmcConversion])
1478  {
1479  fhMCE [pidIndex][kmcConversion] ->Fill(ecluster , GetEventWeight());
1480  fhMCPt [pidIndex][kmcConversion] ->Fill(ptcluster, GetEventWeight());
1481 
1482  fhMCPhi[pidIndex][kmcConversion] ->Fill(ecluster, phicluster, GetEventWeight());
1483  fhMCEta[pidIndex][kmcConversion] ->Fill(ecluster, etacluster, GetEventWeight());
1484 
1485  fhMC2E[pidIndex][kmcConversion] ->Fill(ecluster, eprim , GetEventWeight());
1486  fhMCDeltaE[pidIndex][kmcConversion] ->Fill(ecluster, eprim-ecluster, GetEventWeight());
1487  }
1488  else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay) &&
1489  !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE[pidIndex][kmcPi0Decay])
1490  {
1491  fhMCE [pidIndex][kmcPi0Decay] ->Fill(ecluster , GetEventWeight());
1492  fhMCPt [pidIndex][kmcPi0Decay] ->Fill(ptcluster, GetEventWeight());
1493 
1494  fhMCPhi[pidIndex][kmcPi0Decay] ->Fill(ecluster, phicluster, GetEventWeight());
1495  fhMCEta[pidIndex][kmcPi0Decay] ->Fill(ecluster, etacluster, GetEventWeight());
1496 
1497  fhMC2E[pidIndex][kmcPi0Decay] ->Fill(ecluster, eprim , GetEventWeight());
1498  fhMCDeltaE[pidIndex][kmcPi0Decay] ->Fill(ecluster, eprim-ecluster, GetEventWeight());
1499  }
1500  else if( (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) ||
1501  GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay) ) && fhMCE[pidIndex][kmcOtherDecay])
1502  {
1503  fhMCE [pidIndex][kmcOtherDecay] ->Fill(ecluster , GetEventWeight());
1504  fhMCPt [pidIndex][kmcOtherDecay] ->Fill(ptcluster, GetEventWeight());
1505 
1506  fhMCPhi[pidIndex][kmcOtherDecay] ->Fill(ecluster, phicluster, GetEventWeight());
1507  fhMCEta[pidIndex][kmcOtherDecay] ->Fill(ecluster, etacluster, GetEventWeight());
1508 
1509  fhMC2E [pidIndex][kmcOtherDecay] ->Fill(ecluster, eprim , GetEventWeight());
1510  fhMCDeltaE[pidIndex][kmcOtherDecay] ->Fill(ecluster, eprim-ecluster, GetEventWeight());
1511  }
1512  else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE [pidIndex][kmcPi0])
1513  {
1514  fhMCE [pidIndex][kmcPi0] ->Fill(ecluster , GetEventWeight());
1515  fhMCPt [pidIndex][kmcPi0] ->Fill(ptcluster, GetEventWeight());
1516 
1517  fhMCPhi[pidIndex][kmcPi0] ->Fill(ecluster, phicluster, GetEventWeight());
1518  fhMCEta[pidIndex][kmcPi0] ->Fill(ecluster, etacluster, GetEventWeight());
1519 
1520  fhMC2E[pidIndex][kmcPi0] ->Fill(ecluster, eprim , GetEventWeight());
1521  fhMCDeltaE[pidIndex][kmcPi0] ->Fill(ecluster, eprim-ecluster, GetEventWeight());
1522  }
1523  else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) && fhMCE[pidIndex][kmcEta])
1524  {
1525  fhMCE [pidIndex][kmcEta] ->Fill(ecluster , GetEventWeight());
1526  fhMCPt [pidIndex][kmcEta] ->Fill(ptcluster, GetEventWeight());
1527 
1528  fhMCPhi[pidIndex][kmcEta] ->Fill(ecluster, phicluster, GetEventWeight());
1529  fhMCEta[pidIndex][kmcEta] ->Fill(ecluster, etacluster, GetEventWeight());
1530 
1531  fhMC2E [pidIndex][kmcEta] ->Fill(ecluster, eprim , GetEventWeight());
1532  fhMCDeltaE[pidIndex][kmcEta] ->Fill(ecluster, eprim-ecluster, GetEventWeight());
1533  }
1534  }
1535  else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron) && fhMCE[pidIndex][kmcAntiNeutron])
1536  {
1537  fhMCE [pidIndex][kmcAntiNeutron] ->Fill(ecluster , GetEventWeight());
1538  fhMCPt [pidIndex][kmcAntiNeutron] ->Fill(ptcluster, GetEventWeight());
1539 
1540  fhMCPhi[pidIndex][kmcAntiNeutron] ->Fill(ecluster, phicluster, GetEventWeight());
1541  fhMCEta[pidIndex][kmcAntiNeutron] ->Fill(ecluster, etacluster, GetEventWeight());
1542 
1543  fhMC2E[pidIndex][kmcAntiNeutron] ->Fill(ecluster, eprim , GetEventWeight());
1544  fhMCDeltaE[pidIndex][kmcAntiNeutron] ->Fill(ecluster, eprim-ecluster, GetEventWeight());
1545  }
1546  else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton) && fhMCE[pidIndex][kmcAntiProton])
1547  {
1548  fhMCE [pidIndex][kmcAntiProton] ->Fill(ecluster , GetEventWeight());
1549  fhMCPt [pidIndex][kmcAntiProton] ->Fill(ptcluster, GetEventWeight());
1550 
1551  fhMCPhi[pidIndex][kmcAntiProton] ->Fill(ecluster, phicluster, GetEventWeight());
1552  fhMCEta[pidIndex][kmcAntiProton] ->Fill(ecluster, etacluster, GetEventWeight());
1553 
1554  fhMC2E [pidIndex][kmcAntiProton] ->Fill(ecluster, eprim , GetEventWeight());
1555  fhMCDeltaE[pidIndex][kmcAntiProton] ->Fill(ecluster, eprim-ecluster, GetEventWeight());
1556  }
1557  else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron) && fhMCE[pidIndex][kmcElectron])
1558  {
1559  fhMCE [pidIndex][kmcElectron] ->Fill(ecluster , GetEventWeight());
1560  fhMCPt [pidIndex][kmcElectron] ->Fill(ptcluster, GetEventWeight());
1561 
1562  fhMCPhi[pidIndex][kmcElectron] ->Fill(ecluster, phicluster, GetEventWeight());
1563  fhMCEta[pidIndex][kmcElectron] ->Fill(ecluster, etacluster, GetEventWeight());
1564 
1565  fhMC2E[pidIndex][kmcElectron] ->Fill(ecluster, eprim , GetEventWeight());
1566  fhMCDeltaE[pidIndex][kmcElectron] ->Fill(ecluster, eprim-ecluster, GetEventWeight());
1567  }
1568  else if( fhMCE[pidIndex][kmcOther])
1569  {
1570  fhMCE [pidIndex][kmcOther] ->Fill(ecluster , GetEventWeight());
1571  fhMCPt [pidIndex][kmcOther] ->Fill(ptcluster, GetEventWeight());
1572 
1573  fhMCPhi[pidIndex][kmcOther] ->Fill(ecluster, phicluster, GetEventWeight());
1574  fhMCEta[pidIndex][kmcOther] ->Fill(ecluster, etacluster, GetEventWeight());
1575 
1576  fhMC2E [pidIndex][kmcOther] ->Fill(ecluster, eprim , GetEventWeight());
1577  fhMCDeltaE[pidIndex][kmcOther] ->Fill(ecluster, eprim-ecluster, GetEventWeight());
1578  }
1579  } // Histograms with MC
1580  } // aod loop
1581 }
1582 
1583 //____________________________________________________
1585 //____________________________________________________
1586 void AliAnaElectron::Print(const Option_t * opt) const
1587 {
1588  if(! opt)
1589  return;
1590 
1591  printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
1593 
1594  printf("Calorimeter = %s\n", GetCalorimeterString().Data()) ;
1595  printf(" %2.2f < dEdx < %2.2f \n",fdEdxMin,fdEdxMax) ;
1596  printf(" %2.2f < E/P < %2.2f \n",fEOverPMin,fEOverPMax) ;
1597  printf("Min Distance to Bad Channel = %2.1f\n",fMinDist);
1598  printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
1599  printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
1600  printf("Time Cut: %3.1f < TOF < %3.1f\n", fTimeCutMin, fTimeCutMax);
1601  printf("Number of cells in cluster is > %d \n", fNCellsCut);
1602  printf(" \n") ;
1603 }
1604 
1605 //______________________________________________________
1607 //______________________________________________________
1608 void AliAnaElectron::WeightHistograms(AliVCluster *clus)
1609 {
1610  if(!fFillWeightHistograms || GetMixedEvent()) return;
1611 
1612  AliVCaloCells* cells = 0;
1613  if(GetCalorimeter() == kEMCAL) cells = GetEMCALCells();
1614  else cells = GetPHOSCells();
1615 
1616  // First recalculate energy in case non linearity was applied
1617  Float_t energy = 0;
1618  Float_t ampMax = 0;
1619  for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++) {
1620 
1621  Int_t id = clus->GetCellsAbsId()[ipos];
1622 
1623  //Recalibrate cell energy if needed
1624  Float_t amp = cells->GetCellAmplitude(id);
1626 
1627  energy += amp;
1628 
1629  if(amp> ampMax)
1630  ampMax = amp;
1631  } // energy loop
1632 
1633  if ( energy <= 0 )
1634  {
1635  AliWarning(Form("Wrong calculated energy %f",energy));
1636  return;
1637  }
1638 
1639  //printf("AliAnaElectron::WeightHistograms() - energy %f, ampmax %f, rat %f, lograt %f\n",energy,ampMax,ampMax/energy,TMath::Log(ampMax/energy));
1640  fhEMaxCellClusterRatio ->Fill(energy, ampMax/energy , GetEventWeight());
1641  fhEMaxCellClusterLogRatio->Fill(energy, TMath::Log(ampMax/energy), GetEventWeight());
1642 
1643  // Get the ratio and log ratio to all cells in cluster
1644  for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
1645  {
1646  Int_t id = clus->GetCellsAbsId()[ipos];
1647 
1648  //Recalibrate cell energy if needed
1649  Float_t amp = cells->GetCellAmplitude(id);
1651 
1652  //printf("energy %f, amp %f, rat %f, lograt %f\n",energy,amp,amp/energy,TMath::Log(amp/energy));
1653  fhECellClusterRatio ->Fill(energy, amp/energy , GetEventWeight());
1654  fhECellClusterLogRatio->Fill(energy, TMath::Log(amp/energy), GetEventWeight());
1655  }
1656 
1657  // Recalculate shower shape for different W0
1658  if(GetCalorimeter()==kEMCAL)
1659  {
1660  Float_t l0org = clus->GetM02();
1661  Float_t l1org = clus->GetM20();
1662  Float_t dorg = clus->GetDispersion();
1663 
1664  for(Int_t iw = 0; iw < 14; iw++){
1665 
1666  GetCaloUtils()->GetEMCALRecoUtils()->SetW0(1+iw*0.5);
1667  GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), cells, clus);
1668 
1669  fhLambda0ForW0[iw]->Fill(energy, clus->GetM02(), GetEventWeight());
1670 // fhLambda1ForW0[iw]->Fill(energy, clus->GetM20(), GetEventWeight());
1671 
1672  //printf("\t w %1.1f, l0 %f, l1 %f,\n",3+iw*0.5,clus->GetM02(),clus->GetM20());
1673  } // w0 loop
1674 
1675  // Set the original values back
1676  clus->SetM02(l0org);
1677  clus->SetM20(l1org);
1678  clus->SetDispersion(dorg);
1679  } // EMCAL
1680 }
1681 
1682 
Float_t GetHistoPtMax() const
TH2F * fhDispEtaPhiDiffE[2]
! shower dispersion eta - phi
TH2F * fhEtaPhi05[2]
! Pseudorapidity vs Phi of identified electron for transerse momentum < 0.5
TH2F * fhMCESumEtaPhi[2][6]
! shower dispersion in eta vs phi direction from MC particle
Int_t pdg
Float_t fEOverPMax
Min E/p for electrons, after dEdx cut.
virtual Double_t Eta() const
TH2F * fhSumPhiE[2]
! shower dispersion in phi direction
void FillShowerShapeHistograms(AliVCluster *cluster, Int_t mcTag, Int_t pidTag)
Fill cluster Shower Shape histograms.
virtual void SetNLM(Int_t nlm)
TH2F * fhMCEDispPhi[2][6]
! shower dispersion in phi direction from MC particle
TH2F * fhLambda0ForW0[14]
! L0 for 7 defined w0= 3, 3.5 ... 6 for selected electrons
virtual Double_t Pt() const
Int_t GetHistoNClusterCellMin() const
Float_t GetHistoPtMin() const
double Double_t
Definition: External.C:58
Int_t GetHistoShowerShapeBins() const
Float_t GetHistodEdxMax() const
virtual void AddAODParticle(AliCaloTrackParticle part)
virtual void AddToHistogramsName(TString add)
virtual AliVCaloCells * GetEMCALCells() const
TH2F * fhNLME[2]
! Number of local maxima in cluster vs E
Definition: External.C:236
TH1F * fhMCPt[2][10]
! Number of identified electron vs cluster energy coming from MC particle
TH2F * fhEmbedElectronELambda0MostlyBkg
! Lambda0 vs E for embedded electrons with 50%<fraction<10%
AliEMCALRecoUtils * GetEMCALRecoUtils() const
TH2F * fhMCEDispEtaPhiDiff[2][6]
! shower dispersion in eta -phi direction from MC particle
virtual void GetVertex(Double_t vertex[3]) const
Int_t fNLMCutMax
Remove clusters/cells with number of local maxima larger than this value.
virtual Int_t GetLabel() const
TH2F * fhdEdxvsECutM02
! Matched track dEdx vs cluster E, mild M02 cut
TH2F * fhEOverPvsPCutM02
! Matched track E cluster over P track vs track P, after dEdx cut, mild M02 cut
virtual void SetNCells(Int_t nce)
TH2F * fhMCEta[2][10]
! eta of identified electron coming from MC particle
energy
Definition: HFPtSpectrum.C:44
virtual Float_t GetZvertexCut() const
Maximal number of events for mixin.
TH1F * fhMCE[2][10]
! Number of identified electron vs cluster energy coming from MC particle
TH2F * fhLam1ETRD[2]
! cluster lambda1 vs E, SM covered by TRD
TH2F * fhPhi[2]
! Azimuthal angle of identified electron vs transerse momentum
Bool_t fFillSSHistograms
Fill shower shape histograms.
virtual Bool_t IsTrackMatched(AliVCluster *cluster, AliVEvent *event)
TH2F * fhTimeE[2]
! E vs Time of selected cluster
Float_t fMinDist3
One more cut on distance used for acceptance-efficiency study.
TH2F * fhdEdxvsPCutEOverP
! Matched track dEdx vs track P, cut on EOverP
TH2F * fhMCdEdxvsE[10]
! Matched track dEdx vs cluster E, coming from MC particle
Int_t GetHistoPhiBins() const
TH2F * fhdEdxvsP
! Matched track dEdx vs track P
TH2F * fhMCEOverPvsP[10]
! Matched track E cluster over P track vs track P, after dEdx cut, coming from MC particle ...
TLorentzVector fMomentumMC
! mc particle momentum
virtual void SetDetectorTag(UInt_t d)
TH2F * fhEmbedElectronELambda0FullSignal
! Lambda0 vs E for embedded electrons with more than 90% of the cluster energy
TH2F * fhDispEtaE[2]
! shower dispersion in eta direction
Int_t GetHistoNClusterCellBins() const
virtual TClonesArray * GetOutputAODBranch() const
virtual void SetTime(Float_t tim)
TH2F * fhPhiLam0LowE[2]
! cluster phi vs lambda0, E<2
Selection of electron clusters in calorimeter.
TH2F * fhDispETRD[2]
! cluster dispersion vs E, SM covered by TRD
TH2F * fhMCESphericity[2][6]
! shower sphericity, eta vs phi from MC particle
TH2F * fhdEdxvsECutEOverP
! Matched track dEdx vs cluster E , cut on EOverP
virtual void GetVertex(Double_t v[3]) const
virtual UInt_t GetDetectorTag() const
TH2F * fhMCEDispEta[2][6]
! shower dispersion in eta direction from MC particle
Int_t fAODParticle
Select the type of particle to put in AODs for other analysis.
Int_t GetHistoPOverEBins() const
TH2F * fhECellClusterRatio
! E cell / e cluster vs e cluster for selected electrons
virtual void SetIdentifiedParticleType(Int_t pdg)
Float_t GetHistoPhiMin() const
virtual Double_t E() const
Int_t CheckOrigin(Int_t label, const AliMCEvent *mcevent)
void InitParameters()
Initialize the parameters of the analysis with default values.
TH2F * fhMCELambda0[2][6]
! E vs Lambda0 from MC particle
TH2F * fhECellClusterLogRatio
! log (E cell / E cluster) vs E cluster for selected electrons
TString GetPIDParametersList()
Put data member values in string to keep in output container.
Definition: AliCaloPID.cxx:943
TH2F * fhLam1E[2]
! cluster lambda1 vs E
TH2F * fhdEdxvsPCutM02
! Matched track dEdx vs track P, mild M02 cut
TH2F * fhNCellsLam0HighE[2]
! cluster N Cells vs lambda0, E>2
const Double_t etamin
Base class for CaloTrackCorr analysis algorithms.
virtual TString GetCalorimeterString() const
TObjString * GetAnalysisCuts()
Save parameters used for analysis.
TH2F * fhEOverPvsECutM02CutdEdx
! Matched track E cluster over P track vs cluster E, after dEdx cut and mild M02 cut ...
TH2F * fhSphericityE[2]
! shower sphericity in eta vs phi
Double_t fTimeCutMin
Remove clusters/cells with time smaller than this value, in ns.
Float_t GetHistodEdxMin() const
virtual AliFiducialCut * GetFiducialCut()
int Int_t
Definition: External.C:63
Container for input particle information on CaloTrackCorr package.
virtual AliHistogramRanges * GetHistogramRanges()
Float_t fEOverPMin
Max E/p for electrons, after dEdx cut.
TH2F * fhMCEOverPvsE[10]
! Matched track E cluster over P track vs cluster E, after dEdx cut, coming from MC particle ...
unsigned int UInt_t
Definition: External.C:33
TH2F * fhEOverPvsPCutM02CutdEdx
! Matched track E cluster over P track vs track P, after dEdx cut and mild M02 cut ...
float Float_t
Definition: External.C:68
Double_t fTimeCutMax
Remove clusters/cells with time larger than this value, in ns.
Int_t fNCellsCut
Accept for the analysis clusters with more than fNCellsCut cells.
TH2F * fhMCdEdxvsP[10]
! Matched track dEdx vs track P, coming from MC particle
TH2F * fhEmbeddedSignalFractionEnergy
! Fraction of electron energy of embedded signal vs cluster energy
const Double_t ptmax
void Init()
Init. Check if requested calorimeter is on, if not, abort.
TLorentzVector fMomentum
! cluster momentum
virtual AliEMCALGeometry * GetEMCALGeometry() const
TH2F * fhSumEtaPhiE[2]
! shower dispersion in eta and phi direction
TList * GetCreateOutputObjects()
TH2F * fhEmbedElectronELambda0MostlySignal
! Lambda0 vs E for embedded electrons with 90%<fraction<50%
Bool_t IsInFiducialCut(Float_t eta, Float_t phi, Int_t det) const
virtual void SetTrackLabel(Int_t a, Int_t b)
Float_t GetHistoShowerShapeMin() const
Float_t fdEdxMin
Max dEdx for electrons.
Int_t GetHistodEdxBins() const
virtual AliCalorimeterUtils * GetCaloUtils() const
Int_t GetHistoNClusterCellMax() const
Float_t fdEdxMax
Min dEdx for electrons.
const Double_t ptmin
TH2F * fhEtaLam0HighE[2]
! cluster eta vs lambda0, E>2
TH2F * fhdEdxvsE
! Matched track dEdx vs cluster E
void WeightHistograms(AliVCluster *clus)
Calculate weights and fill histograms.
virtual Double_t GetEventWeight() const
virtual void SetM02(Float_t m02)
TH2F * fhEtaPhi[2]
! Pseudorapidity vs Phi of identified electron for transerse momentum > 0.5
TH2F * fhEmbedElectronELambda0FullBkg
! Lambda0 vs E for embedded electrons with less than 10% of the cluster energy
Int_t fNOriginHistograms
Fill only NOriginHistograms of the 14 defined types.
Int_t GetNumberOfLocalMaxima(AliVCluster *cluster, AliVCaloCells *cells)
Find the number of local maxima in cluster.
TH2F * fhEMaxCellClusterLogRatio
! log (e max cell / e cluster) vs e cluster for selected electrons
TH2F * fhLam0ETRD[2]
! cluster lambda0 vs E, SM covered by TRD
Bool_t fFillWeightHistograms
Fill weigth histograms.
TH2F * fhPhiLam0HighE[2]
! cluster phi vs lambda0, E>2
TH2F * fhEOverPvsECutM02
! Matched track E cluster over P track vs cluster E, after dEdx cut, mild M02 cut ...
virtual TObjArray * GetPHOSClusters() const
Float_t GetHistoEtaMin() const
Float_t fMinDist2
Cuts on Minimal distance to study acceptance evaluation.
virtual void SetLabel(Int_t l)
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
virtual void SetTag(Int_t tag)
virtual void SetSModNumber(Int_t sm)
TH2F * fhMCDeltaE[2][10]
! MC-Reco E distribution coming from MC particle
TH2F * fhMCPhi[2][10]
! Phi of identified electron coming from MC particle
Bool_t fFillOnlySimpleSSHisto
Fill selected cluster histograms, selected SS histograms.
TH2F * fhMCElectronELambda0TwoOverlap
! E vs Lambda0 from MC electrons, 2 particles overlap
TH2F * fhEOverPvsE
! Matched track E cluster over P track vs cluster E, after dEdx cut
virtual void SetCaloLabel(Int_t a, Int_t b)
TH1F * fhPt[2]
! Number of identified electron vs transerse momentum
Float_t GetHistoEtaMax() const
TH2F * fhLam0E[2]
! cluster lambda0 vs E
TH2F * fhMCElectronELambda0NoOverlap
! E vs Lambda0 from MC electrons, no overlap
TH2F * fhEtaLam0LowE[2]
! cluster eta vs lambda0, E<2
TH2F * fhDispPhiE[2]
! shower dispersion in phi direction
Int_t GetHistoPtBins() const
TH2F * fhMC2E[2][10]
! E distribution, Reco vs MC coming from MC particle
AliVTrack * GetMatchedTrack(AliVCluster *cluster, AliVEvent *event, Int_t index=-1) const
virtual Int_t GetModuleNumber(AliCaloTrackParticle *part) const
Float_t GetHistoPOverEMax() const
Bool_t ClusterSelected(AliVCluster *cl, Int_t nMaxima)
void Print(const Option_t *opt) const
Print some relevant parameters set for the analysis.
const Double_t etamax
virtual Double_t Phi() const
void RecalibrateCellAmplitude(Float_t &amp, Int_t calo, Int_t absId) const
Recalculate cell energy if recalibration factor.
virtual AliMCAnalysisUtils * GetMCAnalysisUtils()
TH2F * fhEOverPvsP
! Matched track E cluster over P track vs track P, after dEdx cut
void MakeAnalysisFillAOD()
Do photon analysis selecting electron clusters (or charged non electron) and fill aods...
TH2F * fhEta[2]
! Pseudorapidity of identified electron vs transerse momentum
TH1F * fhE[2]
! Number of identified electron vs energy
virtual void Print(const Option_t *) const
Print some relevant parameters set for the analysis.
Int_t GetHistoTimeBins() const
Float_t GetHistoPOverEMin() const
TH2F * fhSumEtaE[2]
! shower dispersion in eta direction
Float_t fMinDist
Minimal distance to bad channel to accept cluster.
const char Option_t
Definition: External.C:48
Float_t GetHistoTimeMax() const
Int_t CheckCommonAncestor(Int_t index1, Int_t index2, const AliMCEvent *mcevent, Int_t &ancPDG, Int_t &ancStatus, TLorentzVector &momentum, TVector3 &prodVertex)
virtual Int_t GetDataType() const
Float_t GetHistoTimeMin() const
Float_t GetHistoShowerShapeMax() const
virtual void SetM20(Float_t m20)
const Int_t nbins
TH2F * fhNCellsE[2]
! Number of cells in cluster vs E
Float_t GetHistoPhiMax() const
virtual AliCaloTrackReader * GetReader() const
bool Bool_t
Definition: External.C:53
Int_t GetHistoEtaBins() const
TH2F * fhDispE[2]
! cluster dispersion vs E
virtual Int_t GetTag() const
virtual TObjArray * GetEMCALClusters() const
TH2F * fhMaxCellDiffClusterE[2]
! Fraction of energy carried by cell with maximum energy
TH2F * fhDispEtaDispPhiEBin[2][5]
! shower dispersion in eta direction vs phi direction for 5 E bins [0-2],[2-4],[4-6],[6-10],[> 10]
TH2F * fhNCellsLam0LowE[2]
! cluster N cells vs lambda0, E<2
virtual AliVCaloCells * GetPHOSCells() const
virtual void SetDistToBad(Int_t dist)
AliAnaElectron()
Default constructor. Initialize parameters.
Bool_t CheckTagBit(Int_t tag, UInt_t test) const
Int_t GetMaxEnergyCell(AliVCaloCells *cells, AliVCluster *clu, Float_t &fraction) const
For a given CaloCluster, it gets the absId of the cell with maximum energy deposit.
virtual Int_t GetIdentifiedParticleType() const
Int_t nptbins
TH2F * fhEMaxCellClusterRatio
! E max cell / E cluster vs E cluster for selected electrons
virtual AliMixedEvent * GetMixedEvent() const
void MakeAnalysisFillHistograms()
Fill histograms for selected clusters.
TH2F * fhMCElectronELambda0NOverlap
! E vs Lambda0 from MC electrons, N particles overlap
TVector3 fProdVertex
! mc particle production vertex
const Double_t phimin
Int_t GetFirstSMCoveredByTRD() const
Time cut in ns.