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