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