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