]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/PartCorrDep/AliAnaElectron.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Update
[u/mrichter/AliRoot.git] / PWG4 / PartCorrDep / 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 **************************************************************************/
15/* $Id: AliAnaElectron.cxx 28688 2008-09-11 15:04:07Z gconesab $ */
16
17//_________________________________________________________________________
18//
19// Class for the photon identification.
20// Clusters from calorimeters are identified as photons
21// and kept in the AOD. Few histograms produced.
22// Copy of AliAnaPhoton just add electron id.
23//
24// -- Author: Gustavo Conesa (LPSC-IN2P3-CRNS)
25//////////////////////////////////////////////////////////////////////////////
26
27
28// --- ROOT system ---
29#include <TH2F.h>
30#include <TH3D.h>
31#include <TClonesArray.h>
32#include <TObjString.h>
33//#include <Riostream.h>
34#include "TParticle.h"
35#include "TDatabasePDG.h"
36#include "AliVTrack.h"
37
38// --- Analysis system ---
39#include "AliAnaElectron.h"
40#include "AliCaloTrackReader.h"
41#include "AliStack.h"
42#include "AliCaloPID.h"
43#include "AliMCAnalysisUtils.h"
44#include "AliFiducialCut.h"
45#include "AliVCluster.h"
46#include "AliAODMCParticle.h"
47#include "AliMixedEvent.h"
48
49
50ClassImp(AliAnaElectron)
51
78a28af3 52//________________________________
d9105d92 53AliAnaElectron::AliAnaElectron() :
78a28af3 54 AliAnaPartCorrBaseClass(), fCalorimeter(""),
55 fMinDist(0.), fMinDist2(0.), fMinDist3(0.),
56 fTimeCutMin(-1), fTimeCutMax(999999),
57 fNCellsCut(0), fFillSSHistograms(kFALSE),
58 fFillWeightHistograms(kFALSE), fNOriginHistograms(8),
59 fdEdxMin(0.), fdEdxMax (200.),
60 fEOverPMin(0), fEOverPMax (2),
d9105d92 61 // Histograms
78a28af3 62 fhdEdxvsE(0), fhdEdxvsP(0),
63 fhEOverPvsE(0), fhEOverPvsP(0),
64 // Weight studies
65 fhECellClusterRatio(0), fhECellClusterLogRatio(0),
66 fhEMaxCellClusterRatio(0), fhEMaxCellClusterLogRatio(0),
d9105d92 67 // MC histograms
68 // Electron SS MC histograms
69 fhMCElectronELambda0NoOverlap(0),
70 fhMCElectronELambda0TwoOverlap(0), fhMCElectronELambda0NOverlap(0),
71
72 //Embedding
73 fhEmbeddedSignalFractionEnergy(0),
74 fhEmbedElectronELambda0FullSignal(0), fhEmbedElectronELambda0MostlySignal(0),
75 fhEmbedElectronELambda0MostlyBkg(0), fhEmbedElectronELambda0FullBkg(0)
76{
77 //default ctor
78 for(Int_t index = 0; index < 2; index++){
79
42d47cb7 80 fhNCellsE [index] = 0;
81 fhTimeE [index] = 0;
d9105d92 82 fhMaxCellDiffClusterE[index] = 0;
83 fhE [index] = 0;
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
163 if(!IsTrackMatched(calo)) {
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)){
279 fhMCELambda0[pidIndex][mcssPhoton] ->Fill(energy, lambda0);
280
281
282 }//photon no conversion
283 else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron &&
284 !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion))){
285 fhMCELambda0[pidIndex][mcssElectron] ->Fill(energy, lambda0);
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) ){
337 fhMCELambda0[pidIndex][mcssConversion] ->Fill(energy, lambda0);
338 }//conversion photon
339 else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0) ){
340 fhMCELambda0[pidIndex][mcssPi0] ->Fill(energy, lambda0);
341 }//pi0
342 else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta) ){
343 fhMCELambda0[pidIndex][mcssEta] ->Fill(energy, lambda0);
344
345 }//eta
346 else {
347 fhMCELambda0[pidIndex][mcssOther] ->Fill(energy, lambda0);
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
397 Int_t nptbins = GetHistoPtBins(); Float_t ptmax = GetHistoPtMax(); Float_t ptmin = GetHistoPtMin();
398 Int_t nphibins = GetHistoPhiBins(); Float_t phimax = GetHistoPhiMax(); Float_t phimin = GetHistoPhiMin();
399 Int_t netabins = GetHistoEtaBins(); Float_t etamax = GetHistoEtaMax(); Float_t etamin = GetHistoEtaMin();
400 Int_t ssbins = GetHistoShowerShapeBins(); Float_t ssmax = GetHistoShowerShapeMax(); Float_t ssmin = GetHistoShowerShapeMin();
401 Int_t nbins = GetHistoNClusterCellBins(); Int_t nmax = GetHistoNClusterCellMax(); Int_t nmin = GetHistoNClusterCellMin();
402 Int_t ndedxbins = GetHistodEdxBins(); Float_t dedxmax = GetHistodEdxMax(); Float_t dedxmin = GetHistodEdxMin();
403 Int_t nPoverEbins = GetHistoPOverEBins(); Float_t pOverEmax = GetHistoPOverEMax(); Float_t pOverEmin = GetHistoPOverEMin();
42d47cb7 404 Int_t tbins = GetHistoTimeBins() ; Float_t tmax = GetHistoTimeMax(); Float_t tmin = 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
d9105d92 934 Int_t pid = AliCaloPID::kChargedHadron;
935 AliVTrack *track = 0;
936 if(!strcmp("AliESDCaloCluster",Form("%s",calo->ClassName()))){
937 Int_t iESDtrack = calo->GetTrackMatchedIndex();
938 if(iESDtrack<0) printf("AliAnaElectron::MakeAnalysisFillAOD() - Wrong track index\n");
939 AliVEvent * event = GetReader()->GetInputEvent();
940 track = dynamic_cast<AliVTrack*> (event->GetTrack(iESDtrack));
941 }
942 else {
943 track = dynamic_cast<AliVTrack*>(calo->GetTrackMatched(0));
944 }
945
946 if(!track) {
947 printf("AliAnaElectron::MakeAnalysisFillAOD() - Null track");
948 continue;
949 }
950
e416be7d 951 Float_t dEdx = track->GetTPCsignal();
d9105d92 952 fhdEdxvsE->Fill(calo->E(), dEdx);
953 fhdEdxvsP->Fill(track->P(),dEdx);
954
d9105d92 955 if( dEdx < fdEdxMax && dEdx > fdEdxMin) {
956
e416be7d 957 Float_t eOverp = calo->E()/track->P();
958 fhEOverPvsE->Fill(calo->E(), eOverp);
959 fhEOverPvsP->Fill(track->P(), eOverp);
d9105d92 960
961 if( eOverp < fEOverPMax && eOverp > fEOverPMin) {
962
963 pid = AliCaloPID::kElectron;
d9105d92 964 } //E/p
965
966 }// dEdx
967
968 aodph.SetIdentifiedParticleType(pid);
969
970 Int_t pidIndex = 0;// Electron
971 if (pid == AliCaloPID::kElectron) pidIndex = 0;
972 else if(pid == AliCaloPID::kChargedHadron) pidIndex = 1;
973 else continue ;
974
78a28af3 975 //---------------------------------
d9105d92 976 //Fill some shower shape histograms
78a28af3 977 //---------------------------------
1a72f6c5 978
d9105d92 979 FillShowerShapeHistograms(calo,aodph.GetTag(),pid);
980
78a28af3 981 if(pid == AliCaloPID::kElectron)
982 WeightHistograms(calo);
983
984 //-----------------------------------------
d9105d92 985 //PID Shower Shape selection or bit setting
78a28af3 986 //-----------------------------------------
987
d9105d92 988 // Data, PID check on
989 if(IsCaloPIDOn()){
990 //Get most probable PID, 2 options check PID weights
991 //or redo PID, recommended option for EMCal.
992 if(!IsCaloPIDRecalculationOn()){
993 if(GetCaloPID()->GetIdentifiedParticleType(fCalorimeter,calo->GetPID(),mom.E())!=AliCaloPID::kPhoton) continue;
994 }
995 else{
996 if(GetCaloPID()->GetIdentifiedParticleType(fCalorimeter,calo->GetPID(),mom.E())!=AliCaloPID::kPhoton) continue;
997 }
998
999 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetIdentifiedParticleType());
1000
1001 }
1002
1003 // //...............................................
1004 // // Data, PID check off
1005 // else{
1006 // //Set PID bits for later selection (AliAnaPi0 for example)
1007 // //GetPDG already called in SetPIDBits.
1008 // GetCaloPID()->SetPIDBits(fCalorimeter,calo,&aodph, GetCaloUtils());
1009 // if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillAOD() - PID Bits set \n");
1010 // }
1011
1012 if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillAOD() - Photon selection cuts passed: pT %3.2f, pdg %d\n",aodph.Pt(), aodph.GetIdentifiedParticleType());
1013
1014
1015 //FIXME, this to MakeAnalysisFillHistograms ...
1016 Int_t absID = 0;
1017 Float_t maxCellFraction = 0;
1018 AliVCaloCells* cells = 0;
1019
1020 if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
1021 else cells = GetPHOSCells();
1022
1023 absID = GetCaloUtils()->GetMaxEnergyCell(cells, calo,maxCellFraction);
1024
1025 fhMaxCellDiffClusterE[pidIndex]->Fill(aodph.E(),maxCellFraction);
1026 fhNCellsE[pidIndex] ->Fill(aodph.E(),calo->GetNCells());
42d47cb7 1027 fhTimeE[pidIndex] ->Fill(aodph.E(),calo->GetTOF()*1.e9);
1028
d9105d92 1029 //Add AOD with photon object to aod branch
1030 AddAODParticle(aodph);
1031
1032 }//loop
1033
1034 if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillAOD() End fill AODs, with %d entries \n",GetOutputAODBranch()->GetEntriesFast());
1035
1036}
1037
78a28af3 1038//________________________________________________
d9105d92 1039void AliAnaElectron::MakeAnalysisFillHistograms()
1040{
1041 //Fill histograms
1042
1043 //-------------------------------------------------------------------
1044 // Access MC information in stack if requested, check that it exists.
1045 AliStack * stack = 0x0;
1046 TParticle * primary = 0x0;
1047 TClonesArray * mcparticles = 0x0;
1048 AliAODMCParticle * aodprimary = 0x0;
1049
1050 if(IsDataMC()){
1051
1052 if(GetReader()->ReadStack()){
1053 stack = GetMCStack() ;
1054 if(!stack) {
1055 printf("AliAnaElectron::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called? STOP\n");
1056 abort();
1057 }
1058
1059 }
1060 else if(GetReader()->ReadAODMCParticles()){
1061
1062 //Get the list of MC particles
1063 mcparticles = GetReader()->GetAODMCParticles(0);
1064 if(!mcparticles && GetDebug() > 0) {
1065 printf("AliAnaElectron::MakeAnalysisFillHistograms() - Standard MCParticles not available!\n");
1066 }
1067 }
1068 }// is data and MC
1069
1070
1071 // Get vertex
1072 Double_t v[3] = {0,0,0}; //vertex ;
1073 GetReader()->GetVertex(v);
1074 //fhVertex->Fill(v[0],v[1],v[2]);
1075 if(TMath::Abs(v[2]) > GetZvertexCut()) return ; // done elsewhere for Single Event analysis, but there for mixed event
1076
1077 //----------------------------------
1078 //Loop on stored AOD photons
1079 Int_t naod = GetOutputAODBranch()->GetEntriesFast();
1080 if(GetDebug() > 0) printf("AliAnaElectron::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
1081
1082 for(Int_t iaod = 0; iaod < naod ; iaod++){
1083 AliAODPWG4Particle* ph = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
1084 Int_t pdg = ph->GetIdentifiedParticleType();
1085
1086 Int_t pidIndex = 0;// Electron
1087 if (pdg == AliCaloPID::kElectron) pidIndex = 0;
1088 else if(pdg == AliCaloPID::kChargedHadron) pidIndex = 1;
1089 else continue ;
1090
1091 if(ph->GetDetector() != fCalorimeter) continue;
1092
1093 if(GetDebug() > 2)
1094 printf("AliAnaElectron::MakeAnalysisFillHistograms() - ID Electron: pt %f, phi %f, eta %f\n", ph->Pt(),ph->Phi(),ph->Eta()) ;
1095
1096 //................................
1097 //Fill photon histograms
1098 Float_t ptcluster = ph->Pt();
1099 Float_t phicluster = ph->Phi();
1100 Float_t etacluster = ph->Eta();
1101 Float_t ecluster = ph->E();
1102
1103 fhE[pidIndex] ->Fill(ecluster);
1104 fhPt[pidIndex] ->Fill(ptcluster);
1105 fhPhi[pidIndex] ->Fill(ptcluster,phicluster);
1106 fhEta[pidIndex] ->Fill(ptcluster,etacluster);
1107 if (ecluster > 0.5) fhEtaPhi[pidIndex] ->Fill(etacluster, phicluster);
1108 else if(GetMinPt() < 0.5) fhEtaPhi05[pidIndex]->Fill(etacluster, phicluster);
1109
1110 //.......................................
1111 //Play with the MC data if available
1112 if(IsDataMC()){
1113
1114 //....................................................................
1115 // Access MC information in stack if requested, check that it exists.
1116 Int_t label =ph->GetLabel();
1117 if(label < 0) {
1118 if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillHistograms() *** bad label ***: label %d \n", label);
1119 continue;
1120 }
1121
1122 Float_t eprim = 0;
1123 Float_t ptprim = 0;
1124 if(GetReader()->ReadStack()){
1125
1126 if(label >= stack->GetNtrack()) {
1127 if(GetDebug() > 2) printf("AliAnaElectron::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n", label, stack->GetNtrack());
1128 continue ;
1129 }
1130
1131 primary = stack->Particle(label);
1132 if(!primary){
1133 printf("AliAnaElectron::MakeAnalysisFillHistograms() *** no primary ***: label %d \n", label);
1134 continue;
1135 }
1136
1137 eprim = primary->Energy();
1138 ptprim = primary->Pt();
1139
1140 }
1141 else if(GetReader()->ReadAODMCParticles()){
1142 //Check which is the input
1143 if(ph->GetInputFileIndex() == 0){
1144 if(!mcparticles) continue;
1145 if(label >= mcparticles->GetEntriesFast()) {
1146 if(GetDebug() > 2) printf("AliAnaElectron::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n",
1147 label, mcparticles->GetEntriesFast());
1148 continue ;
1149 }
1150 //Get the particle
1151 aodprimary = (AliAODMCParticle*) mcparticles->At(label);
1152
1153 }
1154
1155 if(!aodprimary){
1156 printf("AliAnaElectron::MakeAnalysisFillHistograms() *** no primary ***: label %d \n", label);
1157 continue;
1158 }
1159
1160 eprim = aodprimary->E();
1161 ptprim = aodprimary->Pt();
1162
1163 }
1164
1165 Int_t tag =ph->GetTag();
1166
1167 if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) && fhMCE[pidIndex][mcPhoton])
1168 {
1169 fhMCE [pidIndex][mcPhoton] ->Fill(ecluster);
1170 fhMCPt [pidIndex][mcPhoton] ->Fill(ptcluster);
1171 fhMCPhi[pidIndex][mcPhoton] ->Fill(ecluster,phicluster);
1172 fhMCEta[pidIndex][mcPhoton] ->Fill(ecluster,etacluster);
1173
1174 fhMC2E[pidIndex][mcPhoton] ->Fill(ecluster, eprim);
1175 fhMCDeltaE[pidIndex][mcPhoton] ->Fill(ecluster,eprim-ecluster);
1176
1177 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) && fhMCE[pidIndex][mcConversion])
1178 {
1179 fhMCE [pidIndex][mcConversion] ->Fill(ecluster);
1180 fhMCPt [pidIndex][mcConversion] ->Fill(ptcluster);
1181 fhMCPhi[pidIndex][mcConversion] ->Fill(ecluster,phicluster);
1182 fhMCEta[pidIndex][mcConversion] ->Fill(ecluster,etacluster);
1183
1184 fhMC2E[pidIndex][mcConversion] ->Fill(ecluster, eprim);
1185 fhMCDeltaE[pidIndex][mcConversion] ->Fill(ecluster,eprim-ecluster);
1186
1187 }
1188 else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay) &&
1189 !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE[pidIndex][mcPi0Decay])
1190 {
1191 fhMCE [pidIndex][mcPi0Decay] ->Fill(ecluster);
1192 fhMCPt [pidIndex][mcPi0Decay] ->Fill(ptcluster);
1193 fhMCPhi[pidIndex][mcPi0Decay] ->Fill(ecluster,phicluster);
1194 fhMCEta[pidIndex][mcPi0Decay] ->Fill(ecluster,etacluster);
1195
1196 fhMC2E[pidIndex][mcPi0Decay] ->Fill(ecluster, eprim);
1197 fhMCDeltaE[pidIndex][mcPi0Decay] ->Fill(ecluster,eprim-ecluster);
1198 }
1199 else if( (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) ||
1200 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay) ) && fhMCE[pidIndex][mcOtherDecay])
1201 {
1202 fhMCE [pidIndex][mcOtherDecay] ->Fill(ecluster);
1203 fhMCPt [pidIndex][mcOtherDecay] ->Fill(ptcluster);
1204 fhMCPhi[pidIndex][mcOtherDecay] ->Fill(ecluster,phicluster);
1205 fhMCEta[pidIndex][mcOtherDecay] ->Fill(ecluster,etacluster);
1206
1207 fhMC2E[pidIndex][mcOtherDecay] ->Fill(ecluster, eprim);
1208 fhMCDeltaE[pidIndex][mcOtherDecay] ->Fill(ecluster,eprim-ecluster);
1209 }
1210 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE [pidIndex][mcPi0])
1211 {
1212 fhMCE [pidIndex][mcPi0] ->Fill(ecluster);
1213 fhMCPt [pidIndex][mcPi0] ->Fill(ptcluster);
1214 fhMCPhi[pidIndex][mcPi0] ->Fill(ecluster,phicluster);
1215 fhMCEta[pidIndex][mcPi0] ->Fill(ecluster,etacluster);
1216
1217 fhMC2E[pidIndex][mcPi0] ->Fill(ecluster, eprim);
1218 fhMCDeltaE[pidIndex][mcPi0] ->Fill(ecluster,eprim-ecluster);
1219
1220 }
1221 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) && fhMCE[pidIndex][mcEta])
1222 {
1223 fhMCE [pidIndex][mcEta] ->Fill(ecluster);
1224 fhMCPt [pidIndex][mcEta] ->Fill(ptcluster);
1225 fhMCPhi[pidIndex][mcEta] ->Fill(ecluster,phicluster);
1226 fhMCEta[pidIndex][mcEta] ->Fill(ecluster,etacluster);
1227
1228 fhMC2E[pidIndex][mcEta] ->Fill(ecluster, eprim);
1229 fhMCDeltaE[pidIndex][mcEta] ->Fill(ecluster,eprim-ecluster);
1230
1231 }
1232 }
1233 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron) && fhMCE[pidIndex][mcAntiNeutron])
1234 {
1235 fhMCE [pidIndex][mcAntiNeutron] ->Fill(ecluster);
1236 fhMCPt [pidIndex][mcAntiNeutron] ->Fill(ptcluster);
1237 fhMCPhi[pidIndex][mcAntiNeutron] ->Fill(ecluster,phicluster);
1238 fhMCEta[pidIndex][mcAntiNeutron] ->Fill(ecluster,etacluster);
1239
1240 fhMC2E[pidIndex][mcAntiNeutron] ->Fill(ecluster, eprim);
1241 fhMCDeltaE[pidIndex][mcAntiNeutron] ->Fill(ecluster,eprim-ecluster);
1242
1243 }
1244 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton) && fhMCE[pidIndex][mcAntiProton])
1245 {
1246 fhMCE [pidIndex][mcAntiProton] ->Fill(ecluster);
1247 fhMCPt [pidIndex][mcAntiProton] ->Fill(ptcluster);
1248 fhMCPhi[pidIndex][mcAntiProton] ->Fill(ecluster,phicluster);
1249 fhMCEta[pidIndex][mcAntiProton] ->Fill(ecluster,etacluster);
1250
1251 fhMC2E[pidIndex][mcAntiProton] ->Fill(ecluster, eprim);
1252 fhMCDeltaE[pidIndex][mcAntiProton] ->Fill(ecluster,eprim-ecluster);
1253
1254 }
1255 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron) && fhMCE[pidIndex][mcElectron])
1256 {
1257 fhMCE [pidIndex][mcElectron] ->Fill(ecluster);
1258 fhMCPt [pidIndex][mcElectron] ->Fill(ptcluster);
1259 fhMCPhi[pidIndex][mcElectron] ->Fill(ecluster,phicluster);
1260 fhMCEta[pidIndex][mcElectron] ->Fill(ecluster,etacluster);
1261
1262 fhMC2E[pidIndex][mcElectron] ->Fill(ecluster, eprim);
1263 fhMCDeltaE[pidIndex][mcElectron] ->Fill(ecluster,eprim-ecluster);
1264
1265 }
1266 else if( fhMCE[pidIndex][mcOther]){
1267 fhMCE [pidIndex][mcOther] ->Fill(ecluster);
1268 fhMCPt [pidIndex][mcOther] ->Fill(ptcluster);
1269 fhMCPhi[pidIndex][mcOther] ->Fill(ecluster,phicluster);
1270 fhMCEta[pidIndex][mcOther] ->Fill(ecluster,etacluster);
1271
1272 fhMC2E[pidIndex][mcOther] ->Fill(ecluster, eprim);
1273 fhMCDeltaE[pidIndex][mcOther] ->Fill(ecluster,eprim-ecluster);
1274
1275 }
1276
1277 }//Histograms with MC
1278
1279 }// aod loop
1280
1281}
1282
78a28af3 1283//____________________________________________________
d9105d92 1284void AliAnaElectron::Print(const Option_t * opt) const
1285{
1286 //Print some relevant parameters set for the analysis
1287
1288 if(! opt)
1289 return;
1290
1291 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
1292 AliAnaPartCorrBaseClass::Print(" ");
1293
1294 printf("Calorimeter = %s\n", fCalorimeter.Data()) ;
1295 printf(" %2.2f < dEdx < %2.2f \n",fdEdxMin,fdEdxMax) ;
1296 printf(" %2.2f < E/P < %2.2f \n",fEOverPMin,fEOverPMax) ;
1297 printf("Min Distance to Bad Channel = %2.1f\n",fMinDist);
1298 printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
1299 printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
1300 printf("Time Cut: %3.1f < TOF < %3.1f\n", fTimeCutMin, fTimeCutMax);
1301 printf("Number of cells in cluster is > %d \n", fNCellsCut);
1302 printf(" \n") ;
1303
1304}
78a28af3 1305
1306//__________________________________________________________________________
1307void AliAnaElectron::RecalibrateCellAmplitude(Float_t & amp, const Int_t id)
1308{
1309 //Recaculate cell energy if recalibration factor
1310
1311 Int_t icol = -1; Int_t irow = -1; Int_t iRCU = -1;
1312 Int_t nModule = GetModuleNumberCellIndexes(id,fCalorimeter, icol, irow, iRCU);
1313
1314 if (GetCaloUtils()->IsRecalibrationOn()) {
1315 if(fCalorimeter == "PHOS") {
1316 amp *= GetCaloUtils()->GetPHOSChannelRecalibrationFactor(nModule,icol,irow);
1317 }
1318 else {
1319 amp *= GetCaloUtils()->GetEMCALChannelRecalibrationFactor(nModule,icol,irow);
1320 }
1321 }
1322}
1323
1324//______________________________________________________
1325void AliAnaElectron::WeightHistograms(AliVCluster *clus)
1326{
1327 // Calculate weights and fill histograms
1328
1329 if(!fFillWeightHistograms || GetMixedEvent()) return;
1330
1331 AliVCaloCells* cells = 0;
1332 if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
1333 else cells = GetPHOSCells();
1334
1335 // First recalculate energy in case non linearity was applied
1336 Float_t energy = 0;
1337 Float_t ampMax = 0;
1338 for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++) {
1339
1340 Int_t id = clus->GetCellsAbsId()[ipos];
1341
1342 //Recalibrate cell energy if needed
1343 Float_t amp = cells->GetCellAmplitude(id);
1344 RecalibrateCellAmplitude(amp,id);
1345
1346 energy += amp;
1347
1348 if(amp> ampMax)
1349 ampMax = amp;
1350
1351 } // energy loop
1352
1353 if(energy <=0 ) {
1354 printf("AliAnaCalorimeterQA::WeightHistograms()- Wrong calculated energy %f\n",energy);
1355 return;
1356 }
1357
1a72f6c5 1358 //printf("AliAnaElectron::WeightHistograms() - energy %f, ampmax %f, rat %f, lograt %f\n",energy,ampMax,ampMax/energy,TMath::Log(ampMax/energy));
78a28af3 1359 fhEMaxCellClusterRatio ->Fill(energy,ampMax/energy);
1360 fhEMaxCellClusterLogRatio->Fill(energy,TMath::Log(ampMax/energy));
1361
1362 //Get the ratio and log ratio to all cells in cluster
1363 for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++) {
1364 Int_t id = clus->GetCellsAbsId()[ipos];
1365
1366 //Recalibrate cell energy if needed
1367 Float_t amp = cells->GetCellAmplitude(id);
1368 RecalibrateCellAmplitude(amp,id);
1369
1370 //printf("energy %f, amp %f, rat %f, lograt %f\n",energy,amp,amp/energy,TMath::Log(amp/energy));
1371 fhECellClusterRatio ->Fill(energy,amp/energy);
1372 fhECellClusterLogRatio->Fill(energy,TMath::Log(amp/energy));
1373 }
1374
1375 //Recalculate shower shape for different W0
1376 if(fCalorimeter=="EMCAL"){
1377
1378 Float_t l0org = clus->GetM02();
1379 Float_t l1org = clus->GetM20();
1380 Float_t dorg = clus->GetDispersion();
1381
1a72f6c5 1382 for(Int_t iw = 0; iw < 14; iw++){
1383
1384 GetCaloUtils()->GetEMCALRecoUtils()->SetW0(1+iw*0.5);
78a28af3 1385 GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), cells, clus);
1386
1387 fhLambda0ForW0[iw]->Fill(energy,clus->GetM02());
1a72f6c5 1388 //fhLambda1ForW0[iw]->Fill(energy,clus->GetM20());
78a28af3 1389
1390 //printf("\t w %1.1f, l0 %f, l1 %f,\n",3+iw*0.5,clus->GetM02(),clus->GetM20());
1391
1392 } // w0 loop
1393
1394 // Set the original values back
1395 clus->SetM02(l0org);
1396 clus->SetM20(l1org);
1397 clus->SetDispersion(dorg);
1398
1399 }// EMCAL
1400}
1401
1402