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