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