]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGGA/CaloTrackCorrelations/AliAnaElectron.cxx
remove obsolete lines for CaloPID setting, reduce trigger cut from 8 to 4 GeV
[u/mrichter/AliRoot.git] / PWGGA / CaloTrackCorrelations / AliAnaElectron.cxx
CommitLineData
d9105d92 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 **************************************************************************/
d9105d92 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
49ClassImp(AliAnaElectron)
50
78a28af3 51//________________________________
d9105d92 52AliAnaElectron::AliAnaElectron() :
745913ae 53 AliAnaCaloTrackCorrBaseClass(), fCalorimeter(""),
78a28af3 54 fMinDist(0.), fMinDist2(0.), fMinDist3(0.),
55 fTimeCutMin(-1), fTimeCutMax(999999),
56 fNCellsCut(0), fFillSSHistograms(kFALSE),
57 fFillWeightHistograms(kFALSE), fNOriginHistograms(8),
58 fdEdxMin(0.), fdEdxMax (200.),
59 fEOverPMin(0), fEOverPMax (2),
d9105d92 60 // Histograms
78a28af3 61 fhdEdxvsE(0), fhdEdxvsP(0),
62 fhEOverPvsE(0), fhEOverPvsP(0),
63 // Weight studies
64 fhECellClusterRatio(0), fhECellClusterLogRatio(0),
65 fhEMaxCellClusterRatio(0), fhEMaxCellClusterLogRatio(0),
d9105d92 66 // MC histograms
67 // Electron SS MC histograms
68 fhMCElectronELambda0NoOverlap(0),
69 fhMCElectronELambda0TwoOverlap(0), fhMCElectronELambda0NOverlap(0),
70
71 //Embedding
72 fhEmbeddedSignalFractionEnergy(0),
73 fhEmbedElectronELambda0FullSignal(0), fhEmbedElectronELambda0MostlySignal(0),
74 fhEmbedElectronELambda0MostlyBkg(0), fhEmbedElectronELambda0FullBkg(0)
75{
76 //default ctor
77 for(Int_t index = 0; index < 2; index++){
78
42d47cb7 79 fhNCellsE [index] = 0;
80 fhTimeE [index] = 0;
d9105d92 81 fhMaxCellDiffClusterE[index] = 0;
992b14a7 82 fhE [index] = 0;
83 fhPt [index] = 0;
d9105d92 84 fhPhi [index] = 0;
85 fhEta [index] = 0;
86 fhEtaPhi [index] = 0;
87 fhEtaPhi05[index] = 0;
88
89 // Shower shape histograms
90 fhDispE [index] = 0;
91 fhLam0E [index] = 0;
92 fhLam1E [index] = 0;
93 fhDispETRD[index] = 0;
94 fhLam0ETRD[index] = 0;
95 fhLam1ETRD[index] = 0;
96 fhNCellsLam0LowE [index] = 0;
97 fhNCellsLam0HighE[index] = 0;
98 fhEtaLam0LowE [index] = 0;
99 fhPhiLam0LowE [index] = 0;
100 fhEtaLam0HighE [index] = 0;
101 fhPhiLam0HighE [index] = 0;
102
103 for(Int_t i = 0; i < 10; i++){
104 fhMCPt [index][i] = 0;
105 fhMCE [index][i] = 0;
106 fhMCPhi [index][i] = 0;
107 fhMCEta [index][i] = 0;
108 fhMCDeltaE [index][i] = 0;
109 fhMC2E [index][i] = 0;
110 }
111
112 for(Int_t i = 0; i < 6; i++){
113 fhMCELambda0[index][i]= 0;
114 }
115
116 }
78a28af3 117
118 //Weight studies
1a72f6c5 119 for(Int_t i =0; i < 14; i++){
78a28af3 120 fhLambda0ForW0[i] = 0;
1a72f6c5 121 //fhLambda1ForW0[i] = 0;
78a28af3 122 }
123
d9105d92 124 //Initialize parameters
125 InitParameters();
126
127}
128
78a28af3 129//____________________________________________________________________________
d9105d92 130Bool_t AliAnaElectron::ClusterSelected(AliVCluster* calo, TLorentzVector mom)
131{
132 //Select clusters if they pass different cuts
133 if(GetDebug() > 2)
134 printf("AliAnaElectron::ClusterSelected() Current Event %d; Before selection : E %2.2f, pT %2.2f, Ecl %2.2f, phi %2.2f, eta %2.2f\n",
135 GetReader()->GetEventNumber(),
136 mom.E(), mom.Pt(),calo->E(),mom.Phi()*TMath::RadToDeg(),mom.Eta());
137
138 //.......................................
139 //If too small or big energy, skip it
140 if(mom.E() < GetMinEnergy() || mom.E() > GetMaxEnergy() ) return kFALSE ;
141 if(GetDebug() > 2) printf("\t Cluster %d Pass E Cut \n",calo->GetID());
142
143 //.......................................
144 // TOF cut, BE CAREFUL WITH THIS CUT
145 Double_t tof = calo->GetTOF()*1e9;
146 if(tof < fTimeCutMin || tof > fTimeCutMax) return kFALSE;
147 if(GetDebug() > 2) printf("\t Cluster %d Pass Time Cut \n",calo->GetID());
148
149 //.......................................
150 if(calo->GetNCells() <= fNCellsCut && GetReader()->GetDataType() != AliCaloTrackReader::kMC) return kFALSE;
151 if(GetDebug() > 2) printf("\t Cluster %d Pass NCell Cut \n",calo->GetID());
152
153 //.......................................
154 //Check acceptance selection
155 if(IsFiducialCutOn()){
156 Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
157 if(! in ) return kFALSE ;
158 }
159 if(GetDebug() > 2) printf("Fiducial cut passed \n");
160
161 //.......................................
162 //Skip not matched clusters with tracks
49b5c49b 163 if(!IsTrackMatched(calo, GetReader()->GetInputEvent())) {
d9105d92 164 if(GetDebug() > 2) printf("\t Reject non track-matched clusters\n");
165 return kFALSE ;
166 }
167 else if(GetDebug() > 2) printf(" Track-matching cut passed \n");
168
169 //.......................................
170 //Check Distance to Bad channel, set bit.
171 Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
172 if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
173 if(distBad < fMinDist) {//In bad channel (PHOS cristal size 2.2x2.2 cm), EMCAL ( cell units )
174 return kFALSE ;
175 }
176 else if(GetDebug() > 2) printf("\t Bad channel cut passed %4.2f > %2.2f \n",distBad, fMinDist);
177 //printf("Cluster %d Pass Bad Dist Cut \n",icalo);
178
179 if(GetDebug() > 0)
180 printf("AliAnaElectron::ClusterSelected() Current Event %d; After selection : E %2.2f, pT %2.2f, Ecl %2.2f, phi %2.2f, eta %2.2f\n",
181 GetReader()->GetEventNumber(),
182 mom.E(), mom.Pt(),calo->E(),mom.Phi()*TMath::RadToDeg(),mom.Eta());
183
184 //All checks passed, cluster selected
185 return kTRUE;
186
187}
188
78a28af3 189//__________________________________________________________________________________________________________
190void AliAnaElectron::FillShowerShapeHistograms(AliVCluster* cluster, const Int_t mcTag, const Int_t pidTag)
191{
d9105d92 192
193 //Fill cluster Shower Shape histograms
194
195 if(!fFillSSHistograms || GetMixedEvent()) return;
196
197 Int_t pidIndex = 0;// Electron
198 if (pidTag == AliCaloPID::kElectron) pidIndex = 0;
199 else if(pidTag == AliCaloPID::kChargedHadron) pidIndex = 1;
200 else return;
201
202 Float_t energy = cluster->E();
203 Int_t ncells = cluster->GetNCells();
204 Float_t lambda0 = cluster->GetM02();
205 Float_t lambda1 = cluster->GetM20();
206 Float_t disp = cluster->GetDispersion()*cluster->GetDispersion();
207
208 TLorentzVector mom;
209 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC){
210 cluster->GetMomentum(mom,GetVertex(0)) ;}//Assume that come from vertex in straight line
211 else{
212 Double_t vertex[]={0,0,0};
213 cluster->GetMomentum(mom,vertex) ;
214 }
215
216 Float_t eta = mom.Eta();
217 Float_t phi = mom.Phi();
218 if(phi < 0) phi+=TMath::TwoPi();
219
220 fhLam0E[pidIndex] ->Fill(energy,lambda0);
221 fhLam1E[pidIndex] ->Fill(energy,lambda1);
222 fhDispE[pidIndex] ->Fill(energy,disp);
223
224 if(fCalorimeter == "EMCAL" && GetModuleNumber(cluster) > 5){
225 fhLam0ETRD[pidIndex]->Fill(energy,lambda0);
226 fhLam1ETRD[pidIndex]->Fill(energy,lambda1);
227 fhDispETRD[pidIndex]->Fill(energy,disp);
228 }
229
230 if(energy < 2){
231 fhNCellsLam0LowE[pidIndex] ->Fill(ncells,lambda0);
232 fhEtaLam0LowE[pidIndex] ->Fill(eta, lambda0);
233 fhPhiLam0LowE[pidIndex] ->Fill(phi, lambda0);
234 }
235 else {
236 fhNCellsLam0HighE[pidIndex]->Fill(ncells,lambda0);
237 fhEtaLam0HighE[pidIndex] ->Fill(eta, lambda0);
238 fhPhiLam0HighE[pidIndex] ->Fill(phi, lambda0);
239 }
78a28af3 240
d9105d92 241 if(IsDataMC()){
242
243 AliVCaloCells* cells = 0;
244 if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
245 else cells = GetPHOSCells();
246
247 //Fill histograms to check shape of embedded clusters
248 Float_t fraction = 0;
249 if(GetReader()->IsEmbeddedClusterSelectionOn()){//Only working for EMCAL
250
251 Float_t clusterE = 0; // recalculate in case corrections applied.
252 Float_t cellE = 0;
253 for(Int_t icell = 0; icell < cluster->GetNCells(); icell++){
254 cellE = cells->GetCellAmplitude(cluster->GetCellAbsId(icell));
255 clusterE+=cellE;
256 fraction+=cellE*cluster->GetCellAmplitudeFraction(icell);
257 }
258
259 //Fraction of total energy due to the embedded signal
260 fraction/=clusterE;
261
262 if(GetDebug() > 1 ) printf("AliAnaElectron::FillShowerShapeHistogram() - Energy fraction of embedded signal %2.3f, Energy %2.3f\n",fraction, clusterE);
263
264 fhEmbeddedSignalFractionEnergy->Fill(clusterE,fraction);
265
266 } // embedded fraction
267
268 // Get the fraction of the cluster energy that carries the cell with highest energy
269 Int_t absID =-1 ;
270 Float_t maxCellFraction = 0.;
271
272 absID = GetCaloUtils()->GetMaxEnergyCell(cells, cluster,maxCellFraction);
273
274 // Check the origin and fill histograms
275 if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton) &&
276 !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion) &&
277 !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0) &&
278 !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta)){
c5693f62 279 fhMCELambda0[pidIndex][kmcssPhoton] ->Fill(energy, lambda0);
d9105d92 280
281
282 }//photon no conversion
283 else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron &&
284 !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion))){
c5693f62 285 fhMCELambda0[pidIndex][kmcssElectron] ->Fill(energy, lambda0);
d9105d92 286
287 if(!GetReader()->IsEmbeddedClusterSelectionOn()){
288 //Check particle overlaps in cluster
289
290 //Compare the primary depositing more energy with the rest, if no photon/electron as comon ancestor (conversions), count as other particle
291 Int_t ancPDG = 0, ancStatus = -1;
292 TLorentzVector momentum; TVector3 prodVertex;
293 Int_t ancLabel = 0;
294 Int_t noverlaps = 1;
295 for (UInt_t ilab = 0; ilab < cluster->GetNLabels(); ilab++ ) {
296 ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(cluster->GetLabels()[0],cluster->GetLabels()[ilab], GetReader(),ancPDG,ancStatus,momentum,prodVertex);
297 if(ancPDG!=22 && TMath::Abs(ancPDG)!=11) noverlaps++;
298 }
299
300 if(noverlaps == 1){
301 fhMCElectronELambda0NoOverlap ->Fill(energy, lambda0);
302 }
303 else if(noverlaps == 2){
304 fhMCElectronELambda0TwoOverlap ->Fill(energy, lambda0);
305 }
306 else if(noverlaps > 2){
307 fhMCElectronELambda0NOverlap ->Fill(energy, lambda0);
308 }
309 else {
310 printf("AliAnaElectron::FillShowerShapeHistogram() - n overlaps = %d!!", noverlaps);
311 }
312 }//No embedding
313
314 //Fill histograms to check shape of embedded clusters
315 if(GetReader()->IsEmbeddedClusterSelectionOn()){
316
317 if (fraction > 0.9)
318 {
319 fhEmbedElectronELambda0FullSignal ->Fill(energy, lambda0);
320 }
321 else if(fraction > 0.5)
322 {
323 fhEmbedElectronELambda0MostlySignal ->Fill(energy, lambda0);
324 }
325 else if(fraction > 0.1)
326 {
327 fhEmbedElectronELambda0MostlyBkg ->Fill(energy, lambda0);
328 }
329 else
330 {
331 fhEmbedElectronELambda0FullBkg ->Fill(energy, lambda0);
332 }
333 } // embedded
334 }//electron
335 else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron) &&
336 GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion) ){
c5693f62 337 fhMCELambda0[pidIndex][kmcssConversion] ->Fill(energy, lambda0);
d9105d92 338 }//conversion photon
339 else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0) ){
c5693f62 340 fhMCELambda0[pidIndex][kmcssPi0] ->Fill(energy, lambda0);
d9105d92 341 }//pi0
342 else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta) ){
c5693f62 343 fhMCELambda0[pidIndex][kmcssEta] ->Fill(energy, lambda0);
d9105d92 344
345 }//eta
346 else {
c5693f62 347 fhMCELambda0[pidIndex][kmcssOther] ->Fill(energy, lambda0);
d9105d92 348 }//other particles
349
350 }//MC data
351
352}
353
78a28af3 354//_____________________________________________
d9105d92 355TObjString * AliAnaElectron::GetAnalysisCuts()
356{
357 //Save parameters used for analysis
358 TString parList ; //this will be list of parameters used for this analysis.
359 const Int_t buffersize = 255;
360 char onePar[buffersize] ;
361
362 snprintf(onePar,buffersize,"--- AliAnaElectron ---\n") ;
363 parList+=onePar ;
364 snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
365 parList+=onePar ;
366 snprintf(onePar,buffersize," %2.2f < dEdx < %2.2f \n",fdEdxMin,fdEdxMax) ;
367 parList+=onePar ;
368 snprintf(onePar,buffersize," %2.2f < E/P < %2.2f \n",fEOverPMin, fEOverPMax) ;
369 parList+=onePar ;
370 snprintf(onePar,buffersize,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
371 parList+=onePar ;
372 snprintf(onePar,buffersize,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
373 parList+=onePar ;
374 snprintf(onePar,buffersize,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
375 parList+=onePar ;
376
377 //Get parameters set in base class.
378 parList += GetBaseParametersList() ;
379
380 //Get parameters set in PID class.
381 parList += GetCaloPID()->GetPIDParametersList() ;
382
383 //Get parameters set in FiducialCut class (not available yet)
384 //parlist += GetFidCut()->GetFidCutParametersList()
385
386 return new TObjString(parList) ;
387}
388
78a28af3 389//_______________________________________________
d9105d92 390TList * AliAnaElectron::GetCreateOutputObjects()
391{
392 // Create histograms to be saved in output file and
393 // store them in outputContainer
394 TList * outputContainer = new TList() ;
395 outputContainer->SetName("ElectronHistos") ;
396
745913ae 397 Int_t nptbins = GetHistogramRanges()->GetHistoPtBins(); Float_t ptmax = GetHistogramRanges()->GetHistoPtMax(); Float_t ptmin = GetHistogramRanges()->GetHistoPtMin();
398 Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins(); Float_t phimax = GetHistogramRanges()->GetHistoPhiMax(); Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();
399 Int_t netabins = GetHistogramRanges()->GetHistoEtaBins(); Float_t etamax = GetHistogramRanges()->GetHistoEtaMax(); Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
400 Int_t ssbins = GetHistogramRanges()->GetHistoShowerShapeBins(); Float_t ssmax = GetHistogramRanges()->GetHistoShowerShapeMax(); Float_t ssmin = GetHistogramRanges()->GetHistoShowerShapeMin();
401 Int_t nbins = GetHistogramRanges()->GetHistoNClusterCellBins(); Int_t nmax = GetHistogramRanges()->GetHistoNClusterCellMax(); Int_t nmin = GetHistogramRanges()->GetHistoNClusterCellMin();
402 Int_t ndedxbins = GetHistogramRanges()->GetHistodEdxBins(); Float_t dedxmax = GetHistogramRanges()->GetHistodEdxMax(); Float_t dedxmin = GetHistogramRanges()->GetHistodEdxMin();
403 Int_t nPoverEbins = GetHistogramRanges()->GetHistoPOverEBins(); Float_t pOverEmax = GetHistogramRanges()->GetHistoPOverEMax(); Float_t pOverEmin = GetHistogramRanges()->GetHistoPOverEMin();
404 Int_t tbins = GetHistogramRanges()->GetHistoTimeBins() ; Float_t tmax = GetHistogramRanges()->GetHistoTimeMax(); Float_t tmin = GetHistogramRanges()->GetHistoTimeMin();
d9105d92 405
406 fhdEdxvsE = new TH2F ("hdEdxvsE","matched track <dE/dx> vs cluster E ", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
407 fhdEdxvsE->SetXTitle("E (GeV)");
408 fhdEdxvsE->SetYTitle("<dE/dx>");
409 outputContainer->Add(fhdEdxvsE);
410
411 fhdEdxvsP = new TH2F ("hdEdxvsP","matched track <dE/dx> vs track P ", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
412 fhdEdxvsP->SetXTitle("P (GeV/c)");
413 fhdEdxvsP->SetYTitle("<dE/dx>");
414 outputContainer->Add(fhdEdxvsP);
415
416 fhEOverPvsE = new TH2F ("hEOverPvsE","matched track E/p vs cluster E ", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
417 fhEOverPvsE->SetXTitle("E (GeV)");
418 fhEOverPvsE->SetYTitle("E/p");
419 outputContainer->Add(fhEOverPvsE);
420
421 fhEOverPvsP = new TH2F ("hEOverPvsP","matched track E/p vs track P ", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
422 fhEOverPvsP->SetXTitle("P (GeV/c)");
423 fhEOverPvsP->SetYTitle("E/p");
424 outputContainer->Add(fhEOverPvsP);
425
426
427 TString pidParticle[] = {"Electron","ChargedHadron"} ;
428
1a72f6c5 429 if(fFillWeightHistograms){
430
431 fhECellClusterRatio = new TH2F ("hECellClusterRatio"," cell energy / cluster energy vs cluster energy, for selected electrons",
432 nptbins,ptmin,ptmax, 100,0,1.);
433 fhECellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
434 fhECellClusterRatio->SetYTitle("E_{cell i}/E_{cluster}");
435 outputContainer->Add(fhECellClusterRatio);
436
437 fhECellClusterLogRatio = new TH2F ("hECellClusterLogRatio"," Log(cell energy / cluster energy) vs cluster energy, for selected electrons",
438 nptbins,ptmin,ptmax, 100,-10,0);
439 fhECellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
440 fhECellClusterLogRatio->SetYTitle("Log (E_{max cell}/E_{cluster})");
441 outputContainer->Add(fhECellClusterLogRatio);
442
443 fhEMaxCellClusterRatio = new TH2F ("hEMaxCellClusterRatio"," max cell energy / cluster energy vs cluster energy, for selected electrons",
444 nptbins,ptmin,ptmax, 100,0,1.);
445 fhEMaxCellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
446 fhEMaxCellClusterRatio->SetYTitle("E_{max cell}/E_{cluster}");
447 outputContainer->Add(fhEMaxCellClusterRatio);
448
449 fhEMaxCellClusterLogRatio = new TH2F ("hEMaxCellClusterLogRatio"," Log(max cell energy / cluster energy) vs cluster energy, for selected electrons",
450 nptbins,ptmin,ptmax, 100,-10,0);
451 fhEMaxCellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
452 fhEMaxCellClusterLogRatio->SetYTitle("Log (E_{max cell}/E_{cluster})");
453 outputContainer->Add(fhEMaxCellClusterLogRatio);
454
455 for(Int_t iw = 0; iw < 14; iw++){
456 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),
457 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
458 fhLambda0ForW0[iw]->SetXTitle("E_{cluster}");
459 fhLambda0ForW0[iw]->SetYTitle("#lambda^{2}_{0}");
460 outputContainer->Add(fhLambda0ForW0[iw]);
461
462 // 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),
463 // nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
464 // fhLambda1ForW0[iw]->SetXTitle("E_{cluster}");
465 // fhLambda1ForW0[iw]->SetYTitle("#lambda^{2}_{1}");
466 // outputContainer->Add(fhLambda1ForW0[iw]);
467
468 }
469 }
470
d9105d92 471 for(Int_t pidIndex = 0; pidIndex < 2; pidIndex++){
472
473 //Shower shape
474 if(fFillSSHistograms){
475 fhLam0E[pidIndex] = new TH2F (Form("h%sLam0E",pidParticle[pidIndex].Data()),
476 Form("%s: #lambda_{0}^{2} vs E",pidParticle[pidIndex].Data()),
477 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
478 fhLam0E[pidIndex]->SetYTitle("#lambda_{0}^{2}");
479 fhLam0E[pidIndex]->SetXTitle("E (GeV)");
480 outputContainer->Add(fhLam0E[pidIndex]);
481
482 fhLam1E[pidIndex] = new TH2F (Form("h%sLam1E",pidParticle[pidIndex].Data()),
483 Form("%s: #lambda_{1}^{2} vs E",pidParticle[pidIndex].Data()),
484 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
485 fhLam1E[pidIndex]->SetYTitle("#lambda_{1}^{2}");
486 fhLam1E[pidIndex]->SetXTitle("E (GeV)");
487 outputContainer->Add(fhLam1E[pidIndex]);
488
489 fhDispE[pidIndex] = new TH2F (Form("h%sDispE",pidParticle[pidIndex].Data()),
490 Form("%s: dispersion^{2} vs E",pidParticle[pidIndex].Data()),
491 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
492 fhDispE[pidIndex]->SetYTitle("D^{2}");
493 fhDispE[pidIndex]->SetXTitle("E (GeV) ");
494 outputContainer->Add(fhDispE[pidIndex]);
495
496 if(fCalorimeter == "EMCAL"){
497 fhLam0ETRD[pidIndex] = new TH2F (Form("h%sLam0ETRD",pidParticle[pidIndex].Data()),
498 Form("%s: #lambda_{0}^{2} vs E, EMCAL SM covered by TRD",pidParticle[pidIndex].Data()),
499 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
500 fhLam0ETRD[pidIndex]->SetYTitle("#lambda_{0}^{2}");
501 fhLam0ETRD[pidIndex]->SetXTitle("E (GeV)");
502 outputContainer->Add(fhLam0ETRD[pidIndex]);
503
504 fhLam1ETRD[pidIndex] = new TH2F (Form("h%sLam1ETRD",pidParticle[pidIndex].Data()),
505 Form("%s: #lambda_{1}^{2} vs E, EMCAL SM covered by TRD",pidParticle[pidIndex].Data()),
506 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
507 fhLam1ETRD[pidIndex]->SetYTitle("#lambda_{1}^{2}");
508 fhLam1ETRD[pidIndex]->SetXTitle("E (GeV)");
509 outputContainer->Add(fhLam1ETRD[pidIndex]);
510
511 fhDispETRD[pidIndex] = new TH2F (Form("h%sDispETRD",pidParticle[pidIndex].Data()),
512 Form("%s: dispersion^{2} vs E, EMCAL SM covered by TRD",pidParticle[pidIndex].Data()),
513 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
514 fhDispETRD[pidIndex]->SetYTitle("Dispersion^{2}");
515 fhDispETRD[pidIndex]->SetXTitle("E (GeV) ");
516 outputContainer->Add(fhDispETRD[pidIndex]);
517 }
518
519 fhNCellsLam0LowE[pidIndex] = new TH2F (Form("h%sNCellsLam0LowE",pidParticle[pidIndex].Data()),
520 Form("%s: N_{cells} in cluster vs #lambda_{0}^{2}, E < 2 GeV",pidParticle[pidIndex].Data()),
521 nbins,nmin, nmax, ssbins,ssmin,ssmax);
522 fhNCellsLam0LowE[pidIndex]->SetXTitle("N_{cells}");
523 fhNCellsLam0LowE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
524 outputContainer->Add(fhNCellsLam0LowE[pidIndex]);
525
526 fhNCellsLam0HighE[pidIndex] = new TH2F (Form("h%sNCellsLam0HighE",pidParticle[pidIndex].Data()),
527 Form("%s: N_{cells} in cluster vs #lambda_{0}^{2}, E > 2 GeV",pidParticle[pidIndex].Data()),
528 nbins,nmin, nmax, ssbins,ssmin,ssmax);
529 fhNCellsLam0HighE[pidIndex]->SetXTitle("N_{cells}");
530 fhNCellsLam0HighE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
531 outputContainer->Add(fhNCellsLam0HighE[pidIndex]);
532
533
534 fhEtaLam0LowE[pidIndex] = new TH2F (Form("h%sEtaLam0LowE",pidParticle[pidIndex].Data()),
535 Form("%s: #eta vs #lambda_{0}^{2}, E < 2 GeV",pidParticle[pidIndex].Data()),
536 netabins,etamin,etamax, ssbins,ssmin,ssmax);
537 fhEtaLam0LowE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
538 fhEtaLam0LowE[pidIndex]->SetXTitle("#eta");
539 outputContainer->Add(fhEtaLam0LowE[pidIndex]);
540
541 fhPhiLam0LowE[pidIndex] = new TH2F (Form("h%sPhiLam0LowE",pidParticle[pidIndex].Data()),
542 Form("%s: #phi vs #lambda_{0}^{2}, E < 2 GeV",pidParticle[pidIndex].Data()),
543 nphibins,phimin,phimax, ssbins,ssmin,ssmax);
544 fhPhiLam0LowE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
545 fhPhiLam0LowE[pidIndex]->SetXTitle("#phi");
546 outputContainer->Add(fhPhiLam0LowE[pidIndex]);
547
548 fhEtaLam0HighE[pidIndex] = new TH2F (Form("h%sEtaLam0HighE",pidParticle[pidIndex].Data()),
549 Form("%s: #eta vs #lambda_{0}^{2}, E > 2 GeV",pidParticle[pidIndex].Data()),
550 netabins,etamin,etamax, ssbins,ssmin,ssmax);
551 fhEtaLam0HighE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
552 fhEtaLam0HighE[pidIndex]->SetXTitle("#eta");
553 outputContainer->Add(fhEtaLam0HighE[pidIndex]);
554
555 fhPhiLam0HighE[pidIndex] = new TH2F (Form("h%sPhiLam0HighE",pidParticle[pidIndex].Data()),
556 Form("%s: #phi vs #lambda_{0}^{2}, E > 2 GeV",pidParticle[pidIndex].Data()),
557 nphibins,phimin,phimax, ssbins,ssmin,ssmax);
558 fhPhiLam0HighE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
559 fhPhiLam0HighE[pidIndex]->SetXTitle("#phi");
560 outputContainer->Add(fhPhiLam0HighE[pidIndex]);
561
562 } // Shower shape
78a28af3 563
d9105d92 564 if(IsDataMC()){
565
566 if(fFillSSHistograms){
567
568 TString ptypess[] = { "#gamma","hadron?","#pi^{0}","#eta","#gamma->e^{#pm}","e^{#pm}"} ;
569
570 TString pnamess[] = { "Photon","Hadron","Pi0","Eta","Conversion","Electron"} ;
571
572 for(Int_t i = 0; i < 6; i++){
573
574 fhMCELambda0[pidIndex][i] = new TH2F(Form("h%sELambda0_MC%s",pidParticle[pidIndex].Data(),pnamess[i].Data()),
575 Form("%s like cluster from %s : E vs #lambda_{0}^{2}",pidParticle[pidIndex].Data(),ptypess[i].Data()),
576 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
577 fhMCELambda0[pidIndex][i]->SetYTitle("#lambda_{0}^{2}");
578 fhMCELambda0[pidIndex][i]->SetXTitle("E (GeV)");
579 outputContainer->Add(fhMCELambda0[pidIndex][i]) ;
580
581 }// loop
582 }
583 }
584
585 if(IsCaloPIDOn() && pidIndex > 0) continue;
586
587 fhNCellsE[pidIndex] = new TH2F (Form("h%sNCellsE",pidParticle[pidIndex].Data()),
588 Form("N cells in %s cluster vs E ",pidParticle[pidIndex].Data()),
589 nptbins,ptmin,ptmax, nbins,nmin,nmax);
590 fhNCellsE[pidIndex]->SetXTitle("E (GeV)");
591 fhNCellsE[pidIndex]->SetYTitle("# of cells in cluster");
592 outputContainer->Add(fhNCellsE[pidIndex]);
593
42d47cb7 594 fhTimeE[pidIndex] = new TH2F(Form("h%sTimeE",pidParticle[pidIndex].Data()),
595 Form("Time in %s cluster vs E ",pidParticle[pidIndex].Data())
596 ,nptbins,ptmin,ptmax, tbins,tmin,tmax);
597 fhTimeE[pidIndex]->SetXTitle("E (GeV)");
598 fhTimeE[pidIndex]->SetYTitle(" t (ns)");
599 outputContainer->Add(fhTimeE[pidIndex]);
600
d9105d92 601 fhMaxCellDiffClusterE[pidIndex] = new TH2F (Form("h%sMaxCellDiffClusterE",pidParticle[pidIndex].Data()),
602 Form("%s: energy vs difference of cluster energy - max cell energy / cluster energy, good clusters",pidParticle[pidIndex].Data()),
603 nptbins,ptmin,ptmax, 500,0,1.);
604 fhMaxCellDiffClusterE[pidIndex]->SetXTitle("E_{cluster} (GeV) ");
605 fhMaxCellDiffClusterE[pidIndex]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
606 outputContainer->Add(fhMaxCellDiffClusterE[pidIndex]);
607
608 fhE[pidIndex] = new TH1F(Form("h%sE",pidParticle[pidIndex].Data()),
609 Form("Number of %s over calorimeter vs energy",pidParticle[pidIndex].Data()),
610 nptbins,ptmin,ptmax);
611 fhE[pidIndex]->SetYTitle("N");
612 fhE[pidIndex]->SetXTitle("E_{#gamma}(GeV)");
613 outputContainer->Add(fhE[pidIndex]) ;
614
615 fhPt[pidIndex] = new TH1F(Form("h%sPtElectron",pidParticle[pidIndex].Data()),
616 Form("Number of %s over calorimeter vs p_{T}",pidParticle[pidIndex].Data()),
617 nptbins,ptmin,ptmax);
618 fhPt[pidIndex]->SetYTitle("N");
619 fhPt[pidIndex]->SetXTitle("p_{T #gamma}(GeV/c)");
620 outputContainer->Add(fhPt[pidIndex]) ;
621
622 fhPhi[pidIndex] = new TH2F(Form("h%sPhiElectron",pidParticle[pidIndex].Data()),
623 Form("%s: #phi vs p_{T}",pidParticle[pidIndex].Data()),
624 nptbins,ptmin,ptmax,nphibins,phimin,phimax);
625 fhPhi[pidIndex]->SetYTitle("#phi (rad)");
626 fhPhi[pidIndex]->SetXTitle("p_{T #gamma} (GeV/c)");
627 outputContainer->Add(fhPhi[pidIndex]) ;
628
629 fhEta[pidIndex] = new TH2F(Form("h%sEta",pidParticle[pidIndex].Data()),
630 Form("%s: #eta vs p_{T}",pidParticle[pidIndex].Data()),
631 nptbins,ptmin,ptmax,netabins,etamin,etamax);
632 fhEta[pidIndex]->SetYTitle("#eta");
633 fhEta[pidIndex]->SetXTitle("p_{T #gamma} (GeV/c)");
634 outputContainer->Add(fhEta[pidIndex]) ;
635
636 fhEtaPhi[pidIndex] = new TH2F(Form("h%sEtaPhi",pidParticle[pidIndex].Data()),
637 Form("%s: #eta vs #phi",pidParticle[pidIndex].Data()),
638 netabins,etamin,etamax,nphibins,phimin,phimax);
639 fhEtaPhi[pidIndex]->SetYTitle("#phi (rad)");
640 fhEtaPhi[pidIndex]->SetXTitle("#eta");
641 outputContainer->Add(fhEtaPhi[pidIndex]) ;
642 if(GetMinPt() < 0.5){
643 fhEtaPhi05[pidIndex] = new TH2F(Form("h%sEtaPhi05",pidParticle[pidIndex].Data()),
644 Form("%s: #eta vs #phi, E > 0.5",pidParticle[pidIndex].Data()),
645 netabins,etamin,etamax,nphibins,phimin,phimax);
646 fhEtaPhi05[pidIndex]->SetYTitle("#phi (rad)");
647 fhEtaPhi05[pidIndex]->SetXTitle("#eta");
648 outputContainer->Add(fhEtaPhi05[pidIndex]) ;
649 }
650
651
652 if(IsDataMC()){
653
654
655 TString ptype[] = { "#gamma", "#gamma_{#pi decay}","#gamma_{other decay}", "#pi^{0}","#eta",
656 "e^{#pm}","#gamma->e^{#pm}","hadron?","Anti-N","Anti-P" } ;
657
658 TString pname[] = { "Photon","PhotonPi0Decay","PhotonOtherDecay","Pi0","Eta","Electron",
659 "Conversion", "Hadron", "AntiNeutron","AntiProton" } ;
660
661 for(Int_t i = 0; i < fNOriginHistograms; i++){
662
663 fhMCE[pidIndex][i] = new TH1F(Form("h%sE_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
664 Form("%s like cluster from %s : E ",pidParticle[pidIndex].Data(),ptype[i].Data()),
665 nptbins,ptmin,ptmax);
666 fhMCE[pidIndex][i]->SetXTitle("E (GeV)");
667 outputContainer->Add(fhMCE[pidIndex][i]) ;
668
669 fhMCPt[pidIndex][i] = new TH1F(Form("h%sPt_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
670 Form("%s like cluster from %s : p_{T} ",pidParticle[pidIndex].Data(),ptype[i].Data()),
671 nptbins,ptmin,ptmax);
672 fhMCPt[pidIndex][i]->SetXTitle("p_{T} (GeV/c)");
673 outputContainer->Add(fhMCPt[pidIndex][i]) ;
674
675 fhMCEta[pidIndex][i] = new TH2F(Form("h%sEta_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
676 Form("%s like cluster from %s : #eta ",pidParticle[pidIndex].Data(),ptype[i].Data()),
677 nptbins,ptmin,ptmax,netabins,etamin,etamax);
678 fhMCEta[pidIndex][i]->SetYTitle("#eta");
679 fhMCEta[pidIndex][i]->SetXTitle("E (GeV)");
680 outputContainer->Add(fhMCEta[pidIndex][i]) ;
681
682 fhMCPhi[pidIndex][i] = new TH2F(Form("h%sPhi_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
683 Form("%s like cluster from %s : #phi ",pidParticle[pidIndex].Data(),ptype[i].Data()),
684 nptbins,ptmin,ptmax,nphibins,phimin,phimax);
685 fhMCPhi[pidIndex][i]->SetYTitle("#phi (rad)");
686 fhMCPhi[pidIndex][i]->SetXTitle("E (GeV)");
687 outputContainer->Add(fhMCPhi[pidIndex][i]) ;
688
689
690 fhMCDeltaE[pidIndex][i] = new TH2F (Form("h%sDeltaE_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
691 Form("%s like MC - Reco E from %s",pidParticle[pidIndex].Data(),pname[i].Data()),
692 nptbins,ptmin,ptmax, 200,-50,50);
693 fhMCDeltaE[pidIndex][i]->SetXTitle("#Delta E (GeV)");
694 outputContainer->Add(fhMCDeltaE[pidIndex][i]);
695
696 fhMC2E[pidIndex][i] = new TH2F (Form("h%s2E_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
697 Form("%s like E distribution, reconstructed vs generated from %s",pidParticle[pidIndex].Data(),pname[i].Data()),
698 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
699 fhMC2E[pidIndex][i]->SetXTitle("E_{rec} (GeV)");
700 fhMC2E[pidIndex][i]->SetYTitle("E_{gen} (GeV)");
701 outputContainer->Add(fhMC2E[pidIndex][i]);
702
703 }
704 } // MC
705
706 }// pid Index
707
708
709 if(fFillSSHistograms){
710
711 if(IsDataMC()){
712
713 if(!GetReader()->IsEmbeddedClusterSelectionOn())
714 {
715 fhMCElectronELambda0NoOverlap = new TH2F("hELambda0_MCElectron_NoOverlap",
716 "cluster from Electron : E vs #lambda_{0}^{2}",
717 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
718 fhMCElectronELambda0NoOverlap->SetYTitle("#lambda_{0}^{2}");
719 fhMCElectronELambda0NoOverlap->SetXTitle("E (GeV)");
720 outputContainer->Add(fhMCElectronELambda0NoOverlap) ;
721
722 fhMCElectronELambda0TwoOverlap = new TH2F("hELambda0_MCElectron_TwoOverlap",
723 "cluster from Electron : E vs #lambda_{0}^{2}",
724 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
725 fhMCElectronELambda0TwoOverlap->SetYTitle("#lambda_{0}^{2}");
726 fhMCElectronELambda0TwoOverlap->SetXTitle("E (GeV)");
727 outputContainer->Add(fhMCElectronELambda0TwoOverlap) ;
728
729 fhMCElectronELambda0NOverlap = new TH2F("hELambda0_MCElectron_NOverlap",
730 "cluster from Electron : E vs #lambda_{0}^{2}",
731 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
732 fhMCElectronELambda0NOverlap->SetYTitle("#lambda_{0}^{2}");
733 fhMCElectronELambda0NOverlap->SetXTitle("E (GeV)");
734 outputContainer->Add(fhMCElectronELambda0NOverlap) ;
735
736 } //No embedding
737
738 //Fill histograms to check shape of embedded clusters
739 if(GetReader()->IsEmbeddedClusterSelectionOn())
740 {
741
742 fhEmbeddedSignalFractionEnergy = new TH2F("hEmbeddedSignal_FractionEnergy",
743 "Energy Fraction of embedded signal versus cluster energy",
744 nptbins,ptmin,ptmax,100,0.,1.);
745 fhEmbeddedSignalFractionEnergy->SetYTitle("Fraction");
746 fhEmbeddedSignalFractionEnergy->SetXTitle("E (GeV)");
747 outputContainer->Add(fhEmbeddedSignalFractionEnergy) ;
748
749 fhEmbedElectronELambda0FullSignal = new TH2F("hELambda0_EmbedElectron_FullSignal",
750 "cluster from Electron embedded with more than 90% energy in cluster : E vs #lambda_{0}^{2}",
751 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
752 fhEmbedElectronELambda0FullSignal->SetYTitle("#lambda_{0}^{2}");
753 fhEmbedElectronELambda0FullSignal->SetXTitle("E (GeV)");
754 outputContainer->Add(fhEmbedElectronELambda0FullSignal) ;
755
756 fhEmbedElectronELambda0MostlySignal = new TH2F("hELambda0_EmbedElectron_MostlySignal",
757 "cluster from Electron embedded with 50% to 90% energy in cluster : E vs #lambda_{0}^{2}",
758 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
759 fhEmbedElectronELambda0MostlySignal->SetYTitle("#lambda_{0}^{2}");
760 fhEmbedElectronELambda0MostlySignal->SetXTitle("E (GeV)");
761 outputContainer->Add(fhEmbedElectronELambda0MostlySignal) ;
762
763 fhEmbedElectronELambda0MostlyBkg = new TH2F("hELambda0_EmbedElectron_MostlyBkg",
764 "cluster from Electron embedded with 10% to 50% energy in cluster : E vs #lambda_{0}^{2}",
765 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
766 fhEmbedElectronELambda0MostlyBkg->SetYTitle("#lambda_{0}^{2}");
767 fhEmbedElectronELambda0MostlyBkg->SetXTitle("E (GeV)");
768 outputContainer->Add(fhEmbedElectronELambda0MostlyBkg) ;
769
770 fhEmbedElectronELambda0FullBkg = new TH2F("hELambda0_EmbedElectron_FullBkg",
771 "cluster from Electronm embedded with 0% to 10% energy in cluster : E vs #lambda_{0}^{2}",
772 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
773 fhEmbedElectronELambda0FullBkg->SetYTitle("#lambda_{0}^{2}");
774 fhEmbedElectronELambda0FullBkg->SetXTitle("E (GeV)");
775 outputContainer->Add(fhEmbedElectronELambda0FullBkg) ;
776
777
778 }// embedded histograms
779
780 }//Histos with MC
781
782 }// Fill SS MC histograms
783
784
785 //Store calo PID histograms
786 if(IsCaloPIDOn()){
787 TList * caloPIDHistos = GetCaloPID()->GetCreateOutputObjects() ;
788 for(Int_t i = 0; i < caloPIDHistos->GetEntries(); i++) {
789 outputContainer->Add(caloPIDHistos->At(i)) ;
790 }
791 delete caloPIDHistos;
792 }
793 return outputContainer ;
794
795}
796
78a28af3 797//_________________________
d9105d92 798void AliAnaElectron::Init()
799{
800
801 //Init
802 //Do some checks
803 if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
804 printf("AliAnaElectron::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
805 abort();
806 }
807 else if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
808 printf("AliAnaElectron::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
809 abort();
810 }
811
812}
813
78a28af3 814//___________________________________
d9105d92 815void AliAnaElectron::InitParameters()
816{
817
818 //Initialize the parameters of the analysis.
819 AddToHistogramsName("AnaElectron_");
820
821 fCalorimeter = "EMCAL" ;
822 fMinDist = 2.;
823 fMinDist2 = 4.;
824 fMinDist3 = 5.;
825
826 fTimeCutMin = -1;
827 fTimeCutMax = 9999999;
828 fNCellsCut = 0;
829
830 fdEdxMin = 76.; // for LHC11a, but for LHC11c pass1 56.
831 fdEdxMax = 85.; // for LHC11a, but for LHC11c pass1 64.
832
833 fEOverPMin = 0.8; // for LHC11a, but for LHC11c pass1 0.9
834 fEOverPMax = 1.2; // for LHC11a and LHC11c pass1
835
836}
837
78a28af3 838//_________________________________________
d9105d92 839void AliAnaElectron::MakeAnalysisFillAOD()
840{
841 //Do photon analysis and fill aods
842
843 //Get the vertex
844 Double_t v[3] = {0,0,0}; //vertex ;
845 GetReader()->GetVertex(v);
846
847 //Select the Calorimeter of the photon
848 TObjArray * pl = 0x0;
849 if(fCalorimeter == "PHOS")
850 pl = GetPHOSClusters();
851 else if (fCalorimeter == "EMCAL")
852 pl = GetEMCALClusters();
853
854 if(!pl) {
855 Info("MakeAnalysisFillAOD","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
856 return;
857 }
858
859 //Init arrays, variables, get number of clusters
860 TLorentzVector mom, mom2 ;
861 Int_t nCaloClusters = pl->GetEntriesFast();
862 //List to be used in conversion analysis, to tag the cluster as candidate for conversion
863
864 if(GetDebug() > 0) printf("AliAnaElectron::MakeAnalysisFillAOD() - input %s cluster entries %d\n", fCalorimeter.Data(), nCaloClusters);
865
866 //----------------------------------------------------
867 // Fill AOD with PHOS/EMCAL AliAODPWG4Particle objects
868 //----------------------------------------------------
869 // Loop on clusters
870 for(Int_t icalo = 0; icalo < nCaloClusters; icalo++){
871
872 AliVCluster * calo = (AliVCluster*) (pl->At(icalo));
873 //printf("calo %d, %f\n",icalo,calo->E());
874
875 //Get the index where the cluster comes, to retrieve the corresponding vertex
876 Int_t evtIndex = 0 ;
877 if (GetMixedEvent()) {
878 evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
879 //Get the vertex and check it is not too large in z
880 if(TMath::Abs(GetVertex(evtIndex)[2])> GetZvertexCut()) continue;
881 }
882
883 //Cluster selection, not charged, with photon id and in fiducial cut
884 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC){
885 calo->GetMomentum(mom,GetVertex(evtIndex)) ;}//Assume that come from vertex in straight line
886 else{
887 Double_t vertex[]={0,0,0};
888 calo->GetMomentum(mom,vertex) ;
889 }
890
891 //--------------------------------------
892 // Cluster selection
893 //--------------------------------------
894 if(!ClusterSelected(calo,mom)) continue;
895
896 //----------------------------
897 //Create AOD for analysis
898 //----------------------------
899 AliAODPWG4Particle aodph = AliAODPWG4Particle(mom);
900
901 //...............................................
902 //Set the indeces of the original caloclusters (MC, ID), and calorimeter
903 Int_t label = calo->GetLabel();
904 aodph.SetLabel(label);
905 aodph.SetCaloLabel(calo->GetID(),-1);
906 aodph.SetDetector(fCalorimeter);
907 //printf("Index %d, Id %d, iaod %d\n",icalo, calo->GetID(),GetOutputAODBranch()->GetEntriesFast());
908
909 //...............................................
910 //Set bad channel distance bit
911 Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
912 if (distBad > fMinDist3) aodph.SetDistToBad(2) ;
913 else if(distBad > fMinDist2) aodph.SetDistToBad(1) ;
914 else aodph.SetDistToBad(0) ;
915 //printf("DistBad %f Bit %d\n",distBad, aodph.DistToBad());
916
917 //--------------------------------------------------------------------------------------
918 //Play with the MC stack if available
919 //--------------------------------------------------------------------------------------
920
921 //Check origin of the candidates
922 if(IsDataMC()){
923 aodph.SetTag(GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader(), aodph.GetInputFileIndex()));
924
925 if(GetDebug() > 0)
926 printf("AliAnaElectron::MakeAnalysisFillAOD() - Origin of candidate, bit map %d\n",aodph.GetTag());
927 }//Work with stack also
928
929
930 //-------------------------------------
78a28af3 931 //PID selection via dEdx
d9105d92 932 //-------------------------------------
78a28af3 933
4bfeae64 934 AliVTrack *track = GetCaloUtils()->GetMatchedTrack(calo, GetReader()->GetInputEvent());
935
d9105d92 936 if(!track) {
937 printf("AliAnaElectron::MakeAnalysisFillAOD() - Null track");
938 continue;
939 }
940
e416be7d 941 Float_t dEdx = track->GetTPCsignal();
d9105d92 942 fhdEdxvsE->Fill(calo->E(), dEdx);
943 fhdEdxvsP->Fill(track->P(),dEdx);
944
c5693f62 945 Int_t pid = AliCaloPID::kChargedHadron;
946
d9105d92 947 if( dEdx < fdEdxMax && dEdx > fdEdxMin) {
948
e416be7d 949 Float_t eOverp = calo->E()/track->P();
950 fhEOverPvsE->Fill(calo->E(), eOverp);
951 fhEOverPvsP->Fill(track->P(), eOverp);
d9105d92 952
953 if( eOverp < fEOverPMax && eOverp > fEOverPMin) {
954
955 pid = AliCaloPID::kElectron;
d9105d92 956 } //E/p
957
958 }// dEdx
959
960 aodph.SetIdentifiedParticleType(pid);
961
962 Int_t pidIndex = 0;// Electron
dbf54f1e 963 if(pid == AliCaloPID::kChargedHadron) pidIndex = 1;
d9105d92 964
78a28af3 965 //---------------------------------
d9105d92 966 //Fill some shower shape histograms
78a28af3 967 //---------------------------------
1a72f6c5 968
d9105d92 969 FillShowerShapeHistograms(calo,aodph.GetTag(),pid);
970
78a28af3 971 if(pid == AliCaloPID::kElectron)
972 WeightHistograms(calo);
973
974 //-----------------------------------------
d9105d92 975 //PID Shower Shape selection or bit setting
78a28af3 976 //-----------------------------------------
977
d9105d92 978 // Data, PID check on
3c1d9afb 979 if(IsCaloPIDOn())
980 {
49b5c49b 981 // Get most probable PID, 2 options check bayesian PID weights or redo PID
982 // By default, redo PID
d9105d92 983
3c1d9afb 984 if(GetCaloPID()->GetIdentifiedParticleType(calo)!=AliCaloPID::kPhoton) continue;
49b5c49b 985
d9105d92 986 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetIdentifiedParticleType());
987
988 }
49b5c49b 989
3c1d9afb 990 if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillAOD() - Photon selection cuts passed: pT %3.2f, pdg %d\n",
991 aodph.Pt(), aodph.GetIdentifiedParticleType());
d9105d92 992
d9105d92 993 //FIXME, this to MakeAnalysisFillHistograms ...
994 Int_t absID = 0;
995 Float_t maxCellFraction = 0;
996 AliVCaloCells* cells = 0;
997
998 if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
999 else cells = GetPHOSCells();
1000
1001 absID = GetCaloUtils()->GetMaxEnergyCell(cells, calo,maxCellFraction);
1002
1003 fhMaxCellDiffClusterE[pidIndex]->Fill(aodph.E(),maxCellFraction);
1004 fhNCellsE[pidIndex] ->Fill(aodph.E(),calo->GetNCells());
42d47cb7 1005 fhTimeE[pidIndex] ->Fill(aodph.E(),calo->GetTOF()*1.e9);
1006
d9105d92 1007 //Add AOD with photon object to aod branch
1008 AddAODParticle(aodph);
1009
1010 }//loop
1011
1012 if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillAOD() End fill AODs, with %d entries \n",GetOutputAODBranch()->GetEntriesFast());
1013
1014}
1015
78a28af3 1016//________________________________________________
d9105d92 1017void AliAnaElectron::MakeAnalysisFillHistograms()
1018{
1019 //Fill histograms
1020
1021 //-------------------------------------------------------------------
1022 // Access MC information in stack if requested, check that it exists.
1023 AliStack * stack = 0x0;
1024 TParticle * primary = 0x0;
1025 TClonesArray * mcparticles = 0x0;
1026 AliAODMCParticle * aodprimary = 0x0;
1027
1028 if(IsDataMC()){
1029
1030 if(GetReader()->ReadStack()){
1031 stack = GetMCStack() ;
1032 if(!stack) {
1033 printf("AliAnaElectron::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called? STOP\n");
1034 abort();
1035 }
1036
1037 }
1038 else if(GetReader()->ReadAODMCParticles()){
1039
1040 //Get the list of MC particles
1041 mcparticles = GetReader()->GetAODMCParticles(0);
1042 if(!mcparticles && GetDebug() > 0) {
1043 printf("AliAnaElectron::MakeAnalysisFillHistograms() - Standard MCParticles not available!\n");
1044 }
1045 }
1046 }// is data and MC
1047
1048
1049 // Get vertex
1050 Double_t v[3] = {0,0,0}; //vertex ;
1051 GetReader()->GetVertex(v);
1052 //fhVertex->Fill(v[0],v[1],v[2]);
1053 if(TMath::Abs(v[2]) > GetZvertexCut()) return ; // done elsewhere for Single Event analysis, but there for mixed event
1054
1055 //----------------------------------
1056 //Loop on stored AOD photons
1057 Int_t naod = GetOutputAODBranch()->GetEntriesFast();
1058 if(GetDebug() > 0) printf("AliAnaElectron::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
1059
3c1d9afb 1060 for(Int_t iaod = 0; iaod < naod ; iaod++)
1061 {
d9105d92 1062 AliAODPWG4Particle* ph = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
1063 Int_t pdg = ph->GetIdentifiedParticleType();
1064
1065 Int_t pidIndex = 0;// Electron
1066 if (pdg == AliCaloPID::kElectron) pidIndex = 0;
1067 else if(pdg == AliCaloPID::kChargedHadron) pidIndex = 1;
1068 else continue ;
1069
1070 if(ph->GetDetector() != fCalorimeter) continue;
1071
1072 if(GetDebug() > 2)
1073 printf("AliAnaElectron::MakeAnalysisFillHistograms() - ID Electron: pt %f, phi %f, eta %f\n", ph->Pt(),ph->Phi(),ph->Eta()) ;
1074
1075 //................................
1076 //Fill photon histograms
1077 Float_t ptcluster = ph->Pt();
1078 Float_t phicluster = ph->Phi();
1079 Float_t etacluster = ph->Eta();
1080 Float_t ecluster = ph->E();
1081
1082 fhE[pidIndex] ->Fill(ecluster);
1083 fhPt[pidIndex] ->Fill(ptcluster);
1084 fhPhi[pidIndex] ->Fill(ptcluster,phicluster);
1085 fhEta[pidIndex] ->Fill(ptcluster,etacluster);
1086 if (ecluster > 0.5) fhEtaPhi[pidIndex] ->Fill(etacluster, phicluster);
1087 else if(GetMinPt() < 0.5) fhEtaPhi05[pidIndex]->Fill(etacluster, phicluster);
1088
1089 //.......................................
1090 //Play with the MC data if available
1091 if(IsDataMC()){
1092
1093 //....................................................................
1094 // Access MC information in stack if requested, check that it exists.
1095 Int_t label =ph->GetLabel();
1096 if(label < 0) {
1097 if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillHistograms() *** bad label ***: label %d \n", label);
1098 continue;
1099 }
1100
1101 Float_t eprim = 0;
1102 Float_t ptprim = 0;
1103 if(GetReader()->ReadStack()){
1104
1105 if(label >= stack->GetNtrack()) {
1106 if(GetDebug() > 2) printf("AliAnaElectron::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n", label, stack->GetNtrack());
1107 continue ;
1108 }
1109
1110 primary = stack->Particle(label);
1111 if(!primary){
1112 printf("AliAnaElectron::MakeAnalysisFillHistograms() *** no primary ***: label %d \n", label);
1113 continue;
1114 }
1115
1116 eprim = primary->Energy();
1117 ptprim = primary->Pt();
1118
1119 }
1120 else if(GetReader()->ReadAODMCParticles()){
1121 //Check which is the input
1122 if(ph->GetInputFileIndex() == 0){
1123 if(!mcparticles) continue;
1124 if(label >= mcparticles->GetEntriesFast()) {
1125 if(GetDebug() > 2) printf("AliAnaElectron::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n",
1126 label, mcparticles->GetEntriesFast());
1127 continue ;
1128 }
1129 //Get the particle
1130 aodprimary = (AliAODMCParticle*) mcparticles->At(label);
1131
1132 }
1133
1134 if(!aodprimary){
1135 printf("AliAnaElectron::MakeAnalysisFillHistograms() *** no primary ***: label %d \n", label);
1136 continue;
1137 }
1138
1139 eprim = aodprimary->E();
1140 ptprim = aodprimary->Pt();
1141
1142 }
1143
1144 Int_t tag =ph->GetTag();
1145
c5693f62 1146 if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) && fhMCE[pidIndex][kmcPhoton])
d9105d92 1147 {
c5693f62 1148 fhMCE [pidIndex][kmcPhoton] ->Fill(ecluster);
1149 fhMCPt [pidIndex][kmcPhoton] ->Fill(ptcluster);
1150 fhMCPhi[pidIndex][kmcPhoton] ->Fill(ecluster,phicluster);
1151 fhMCEta[pidIndex][kmcPhoton] ->Fill(ecluster,etacluster);
d9105d92 1152
c5693f62 1153 fhMC2E[pidIndex][kmcPhoton] ->Fill(ecluster, eprim);
1154 fhMCDeltaE[pidIndex][kmcPhoton] ->Fill(ecluster,eprim-ecluster);
d9105d92 1155
c5693f62 1156 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) && fhMCE[pidIndex][kmcConversion])
d9105d92 1157 {
c5693f62 1158 fhMCE [pidIndex][kmcConversion] ->Fill(ecluster);
1159 fhMCPt [pidIndex][kmcConversion] ->Fill(ptcluster);
1160 fhMCPhi[pidIndex][kmcConversion] ->Fill(ecluster,phicluster);
1161 fhMCEta[pidIndex][kmcConversion] ->Fill(ecluster,etacluster);
d9105d92 1162
c5693f62 1163 fhMC2E[pidIndex][kmcConversion] ->Fill(ecluster, eprim);
1164 fhMCDeltaE[pidIndex][kmcConversion] ->Fill(ecluster,eprim-ecluster);
d9105d92 1165
1166 }
1167 else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay) &&
c5693f62 1168 !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE[pidIndex][kmcPi0Decay])
d9105d92 1169 {
c5693f62 1170 fhMCE [pidIndex][kmcPi0Decay] ->Fill(ecluster);
1171 fhMCPt [pidIndex][kmcPi0Decay] ->Fill(ptcluster);
1172 fhMCPhi[pidIndex][kmcPi0Decay] ->Fill(ecluster,phicluster);
1173 fhMCEta[pidIndex][kmcPi0Decay] ->Fill(ecluster,etacluster);
d9105d92 1174
c5693f62 1175 fhMC2E[pidIndex][kmcPi0Decay] ->Fill(ecluster, eprim);
1176 fhMCDeltaE[pidIndex][kmcPi0Decay] ->Fill(ecluster,eprim-ecluster);
d9105d92 1177 }
1178 else if( (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) ||
c5693f62 1179 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay) ) && fhMCE[pidIndex][kmcOtherDecay])
d9105d92 1180 {
c5693f62 1181 fhMCE [pidIndex][kmcOtherDecay] ->Fill(ecluster);
1182 fhMCPt [pidIndex][kmcOtherDecay] ->Fill(ptcluster);
1183 fhMCPhi[pidIndex][kmcOtherDecay] ->Fill(ecluster,phicluster);
1184 fhMCEta[pidIndex][kmcOtherDecay] ->Fill(ecluster,etacluster);
d9105d92 1185
c5693f62 1186 fhMC2E[pidIndex][kmcOtherDecay] ->Fill(ecluster, eprim);
1187 fhMCDeltaE[pidIndex][kmcOtherDecay] ->Fill(ecluster,eprim-ecluster);
d9105d92 1188 }
c5693f62 1189 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE [pidIndex][kmcPi0])
d9105d92 1190 {
c5693f62 1191 fhMCE [pidIndex][kmcPi0] ->Fill(ecluster);
1192 fhMCPt [pidIndex][kmcPi0] ->Fill(ptcluster);
1193 fhMCPhi[pidIndex][kmcPi0] ->Fill(ecluster,phicluster);
1194 fhMCEta[pidIndex][kmcPi0] ->Fill(ecluster,etacluster);
d9105d92 1195
c5693f62 1196 fhMC2E[pidIndex][kmcPi0] ->Fill(ecluster, eprim);
1197 fhMCDeltaE[pidIndex][kmcPi0] ->Fill(ecluster,eprim-ecluster);
d9105d92 1198
1199 }
c5693f62 1200 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) && fhMCE[pidIndex][kmcEta])
d9105d92 1201 {
c5693f62 1202 fhMCE [pidIndex][kmcEta] ->Fill(ecluster);
1203 fhMCPt [pidIndex][kmcEta] ->Fill(ptcluster);
1204 fhMCPhi[pidIndex][kmcEta] ->Fill(ecluster,phicluster);
1205 fhMCEta[pidIndex][kmcEta] ->Fill(ecluster,etacluster);
d9105d92 1206
c5693f62 1207 fhMC2E[pidIndex][kmcEta] ->Fill(ecluster, eprim);
1208 fhMCDeltaE[pidIndex][kmcEta] ->Fill(ecluster,eprim-ecluster);
d9105d92 1209
1210 }
1211 }
c5693f62 1212 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron) && fhMCE[pidIndex][kmcAntiNeutron])
d9105d92 1213 {
c5693f62 1214 fhMCE [pidIndex][kmcAntiNeutron] ->Fill(ecluster);
1215 fhMCPt [pidIndex][kmcAntiNeutron] ->Fill(ptcluster);
1216 fhMCPhi[pidIndex][kmcAntiNeutron] ->Fill(ecluster,phicluster);
1217 fhMCEta[pidIndex][kmcAntiNeutron] ->Fill(ecluster,etacluster);
d9105d92 1218
c5693f62 1219 fhMC2E[pidIndex][kmcAntiNeutron] ->Fill(ecluster, eprim);
1220 fhMCDeltaE[pidIndex][kmcAntiNeutron] ->Fill(ecluster,eprim-ecluster);
d9105d92 1221
1222 }
c5693f62 1223 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton) && fhMCE[pidIndex][kmcAntiProton])
d9105d92 1224 {
c5693f62 1225 fhMCE [pidIndex][kmcAntiProton] ->Fill(ecluster);
1226 fhMCPt [pidIndex][kmcAntiProton] ->Fill(ptcluster);
1227 fhMCPhi[pidIndex][kmcAntiProton] ->Fill(ecluster,phicluster);
1228 fhMCEta[pidIndex][kmcAntiProton] ->Fill(ecluster,etacluster);
d9105d92 1229
c5693f62 1230 fhMC2E[pidIndex][kmcAntiProton] ->Fill(ecluster, eprim);
1231 fhMCDeltaE[pidIndex][kmcAntiProton] ->Fill(ecluster,eprim-ecluster);
d9105d92 1232
1233 }
c5693f62 1234 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron) && fhMCE[pidIndex][kmcElectron])
d9105d92 1235 {
c5693f62 1236 fhMCE [pidIndex][kmcElectron] ->Fill(ecluster);
1237 fhMCPt [pidIndex][kmcElectron] ->Fill(ptcluster);
1238 fhMCPhi[pidIndex][kmcElectron] ->Fill(ecluster,phicluster);
1239 fhMCEta[pidIndex][kmcElectron] ->Fill(ecluster,etacluster);
d9105d92 1240
c5693f62 1241 fhMC2E[pidIndex][kmcElectron] ->Fill(ecluster, eprim);
1242 fhMCDeltaE[pidIndex][kmcElectron] ->Fill(ecluster,eprim-ecluster);
d9105d92 1243
1244 }
c5693f62 1245 else if( fhMCE[pidIndex][kmcOther]){
1246 fhMCE [pidIndex][kmcOther] ->Fill(ecluster);
1247 fhMCPt [pidIndex][kmcOther] ->Fill(ptcluster);
1248 fhMCPhi[pidIndex][kmcOther] ->Fill(ecluster,phicluster);
1249 fhMCEta[pidIndex][kmcOther] ->Fill(ecluster,etacluster);
d9105d92 1250
c5693f62 1251 fhMC2E[pidIndex][kmcOther] ->Fill(ecluster, eprim);
1252 fhMCDeltaE[pidIndex][kmcOther] ->Fill(ecluster,eprim-ecluster);
d9105d92 1253
1254 }
1255
1256 }//Histograms with MC
1257
1258 }// aod loop
1259
1260}
1261
78a28af3 1262//____________________________________________________
d9105d92 1263void AliAnaElectron::Print(const Option_t * opt) const
1264{
1265 //Print some relevant parameters set for the analysis
1266
1267 if(! opt)
1268 return;
1269
1270 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
745913ae 1271 AliAnaCaloTrackCorrBaseClass::Print(" ");
d9105d92 1272
1273 printf("Calorimeter = %s\n", fCalorimeter.Data()) ;
1274 printf(" %2.2f < dEdx < %2.2f \n",fdEdxMin,fdEdxMax) ;
1275 printf(" %2.2f < E/P < %2.2f \n",fEOverPMin,fEOverPMax) ;
1276 printf("Min Distance to Bad Channel = %2.1f\n",fMinDist);
1277 printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
1278 printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
1279 printf("Time Cut: %3.1f < TOF < %3.1f\n", fTimeCutMin, fTimeCutMax);
1280 printf("Number of cells in cluster is > %d \n", fNCellsCut);
1281 printf(" \n") ;
1282
1283}
78a28af3 1284
78a28af3 1285//______________________________________________________
1286void AliAnaElectron::WeightHistograms(AliVCluster *clus)
1287{
1288 // Calculate weights and fill histograms
1289
1290 if(!fFillWeightHistograms || GetMixedEvent()) return;
1291
1292 AliVCaloCells* cells = 0;
1293 if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
1294 else cells = GetPHOSCells();
1295
1296 // First recalculate energy in case non linearity was applied
1297 Float_t energy = 0;
1298 Float_t ampMax = 0;
1299 for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++) {
1300
1301 Int_t id = clus->GetCellsAbsId()[ipos];
1302
1303 //Recalibrate cell energy if needed
1304 Float_t amp = cells->GetCellAmplitude(id);
dbba06ca 1305 GetCaloUtils()->RecalibrateCellAmplitude(amp,fCalorimeter, id);
78a28af3 1306
1307 energy += amp;
1308
1309 if(amp> ampMax)
1310 ampMax = amp;
1311
1312 } // energy loop
1313
1314 if(energy <=0 ) {
1315 printf("AliAnaCalorimeterQA::WeightHistograms()- Wrong calculated energy %f\n",energy);
1316 return;
1317 }
1318
1a72f6c5 1319 //printf("AliAnaElectron::WeightHistograms() - energy %f, ampmax %f, rat %f, lograt %f\n",energy,ampMax,ampMax/energy,TMath::Log(ampMax/energy));
78a28af3 1320 fhEMaxCellClusterRatio ->Fill(energy,ampMax/energy);
1321 fhEMaxCellClusterLogRatio->Fill(energy,TMath::Log(ampMax/energy));
1322
1323 //Get the ratio and log ratio to all cells in cluster
1324 for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++) {
1325 Int_t id = clus->GetCellsAbsId()[ipos];
1326
1327 //Recalibrate cell energy if needed
1328 Float_t amp = cells->GetCellAmplitude(id);
dbba06ca 1329 GetCaloUtils()->RecalibrateCellAmplitude(amp, fCalorimeter, id);
78a28af3 1330
1331 //printf("energy %f, amp %f, rat %f, lograt %f\n",energy,amp,amp/energy,TMath::Log(amp/energy));
1332 fhECellClusterRatio ->Fill(energy,amp/energy);
1333 fhECellClusterLogRatio->Fill(energy,TMath::Log(amp/energy));
1334 }
1335
1336 //Recalculate shower shape for different W0
1337 if(fCalorimeter=="EMCAL"){
1338
1339 Float_t l0org = clus->GetM02();
1340 Float_t l1org = clus->GetM20();
1341 Float_t dorg = clus->GetDispersion();
1342
1a72f6c5 1343 for(Int_t iw = 0; iw < 14; iw++){
1344
1345 GetCaloUtils()->GetEMCALRecoUtils()->SetW0(1+iw*0.5);
78a28af3 1346 GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), cells, clus);
1347
1348 fhLambda0ForW0[iw]->Fill(energy,clus->GetM02());
1a72f6c5 1349 //fhLambda1ForW0[iw]->Fill(energy,clus->GetM20());
78a28af3 1350
1351 //printf("\t w %1.1f, l0 %f, l1 %f,\n",3+iw*0.5,clus->GetM02(),clus->GetM20());
1352
1353 } // w0 loop
1354
1355 // Set the original values back
1356 clus->SetM02(l0org);
1357 clus->SetM20(l1org);
1358 clus->SetDispersion(dorg);
1359
1360 }// EMCAL
1361}
1362
1363