1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is 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 **************************************************************************/
16 //_________________________________________________________________________
17 // Class for the analysis of high pT pi0 event by event
18 // Pi0/Eta identified by one of the following:
19 // -Invariant mass of 2 cluster in calorimeter
20 // -Shower shape analysis in calorimeter
21 // -Invariant mass of one cluster in calorimeter and one photon reconstructed in CTS
23 // -- Author: Gustavo Conesa (LNF-INFN) & Raphaelle Ichou (SUBATECH)
24 //////////////////////////////////////////////////////////////////////////////
27 // --- ROOT system ---
29 #include <TClonesArray.h>
30 #include <TObjString.h>
32 // --- Analysis system ---
33 #include "AliAnaPi0EbE.h"
34 #include "AliCaloTrackReader.h"
35 #include "AliIsolationCut.h"
36 #include "AliNeutralMesonSelection.h"
37 #include "AliCaloPID.h"
38 #include "AliMCAnalysisUtils.h"
40 #include "AliFiducialCut.h"
41 #include "TParticle.h"
42 #include "AliVCluster.h"
43 #include "AliESDEvent.h"
44 #include "AliAODEvent.h"
45 #include "AliAODMCParticle.h"
47 ClassImp(AliAnaPi0EbE)
49 //____________________________
50 AliAnaPi0EbE::AliAnaPi0EbE() :
51 AliAnaCaloTrackCorrBaseClass(),fAnaType(kIMCalo), fCalorimeter(""),
52 fMinDist(0.),fMinDist2(0.), fMinDist3(0.),
53 fNLMCutMin(-1), fNLMCutMax(10),
54 fUseSplitAsyCut(kFALSE),
55 fTimeCutMin(-10000), fTimeCutMax(10000),
56 fFillPileUpHistograms(0),
57 fFillWeightHistograms(kFALSE), fFillTMHisto(0),
58 fFillSelectClHisto(0), fFillOnlySimpleSSHisto(1),
59 fInputAODGammaConvName(""),
62 fhEEta(0), fhEPhi(0), fhEtaPhi(0),
63 fhPtReject(0), fhEReject(0),
64 fhEEtaReject(0), fhEPhiReject(0), fhEtaPhiReject(0),
65 fhMass(0), fhAsymmetry(0),
66 fhSelectedMass(0), fhSelectedAsymmetry(0),
67 fhPtDecay(0), fhEDecay(0),
68 // Shower shape histos
69 fhEDispersion(0), fhELambda0(0), fhELambda1(0),
70 fhELambda0NoTRD(0), fhELambda0FracMaxCellCut(0),
71 fhEFracMaxCell(0), fhEFracMaxCellNoTRD(0),
72 fhENCells(0), fhETime(0), fhEPairDiffTime(0),
73 fhDispEtaE(0), fhDispPhiE(0),
74 fhSumEtaE(0), fhSumPhiE(0), fhSumEtaPhiE(0),
75 fhDispEtaPhiDiffE(0), fhSphericityE(0),
80 fhMCEReject(), fhMCPtReject(),
81 fhMCPi0PtGenRecoFraction(0), fhMCEtaPtGenRecoFraction(0),
82 fhMCPi0DecayPt(0), fhMCPi0DecayPtFraction(0),
83 fhMCEtaDecayPt(0), fhMCEtaDecayPtFraction(0),
85 fhMassPairMCPi0(0), fhMassPairMCEta(0),
86 fhAnglePairMCPi0(0), fhAnglePairMCEta(0),
88 fhECellClusterRatio(0), fhECellClusterLogRatio(0),
89 fhEMaxCellClusterRatio(0), fhEMaxCellClusterLogRatio(0),
90 fhTrackMatchedDEta(0), fhTrackMatchedDPhi(0), fhTrackMatchedDEtaDPhi(0),
91 fhTrackMatchedMCParticle(0), fhdEdx(0),
92 fhEOverP(0), fhEOverPNoTRD(0),
93 // Number of local maxima in cluster
96 fhTimeENoCut(0), fhTimeESPD(0), fhTimeESPDMulti(0),
97 fhTimeNPileUpVertSPD(0), fhTimeNPileUpVertTrack(0),
98 fhTimeNPileUpVertContributors(0),
99 fhTimePileUpMainVertexZDistance(0), fhTimePileUpMainVertexZDiamond(0)
103 for(Int_t i = 0; i < 6; i++)
109 fhEMCLambda0 [i] = 0;
110 fhEMCLambda0NoTRD [i] = 0;
111 fhEMCLambda0FracMaxCellCut[i]= 0;
112 fhEMCFracMaxCell [i] = 0;
113 fhEMCLambda1 [i] = 0;
114 fhEMCDispersion [i] = 0;
116 fhMCEDispEta [i] = 0;
117 fhMCEDispPhi [i] = 0;
118 fhMCESumEtaPhi [i] = 0;
119 fhMCEDispEtaPhiDiff[i] = 0;
120 fhMCESphericity [i] = 0;
121 fhMCEAsymmetry [i] = 0;
123 for(Int_t j = 0; j < 7; j++)
125 fhMCLambda0DispEta [j][i] = 0;
126 fhMCLambda0DispPhi [j][i] = 0;
127 fhMCDispEtaDispPhi [j][i] = 0;
128 fhMCAsymmetryLambda0 [j][i] = 0;
129 fhMCAsymmetryDispEta [j][i] = 0;
130 fhMCAsymmetryDispPhi [j][i] = 0;
134 for(Int_t j = 0; j < 7; j++)
136 fhLambda0DispEta [j] = 0;
137 fhLambda0DispPhi [j] = 0;
138 fhDispEtaDispPhi [j] = 0;
139 fhAsymmetryLambda0 [j] = 0;
140 fhAsymmetryDispEta [j] = 0;
141 fhAsymmetryDispPhi [j] = 0;
143 fhPtPi0PileUp [j] = 0;
146 for(Int_t i = 0; i < 3; i++)
148 fhELambda0LocMax [i] = 0;
149 fhELambda1LocMax [i] = 0;
150 fhEDispersionLocMax [i] = 0;
151 fhEDispEtaLocMax [i] = 0;
152 fhEDispPhiLocMax [i] = 0;
153 fhESumEtaPhiLocMax [i] = 0;
154 fhEDispEtaPhiDiffLocMax[i] = 0;
155 fhESphericityLocMax [i] = 0;
156 fhEAsymmetryLocMax [i] = 0;
160 for(Int_t i =0; i < 14; i++){
161 fhLambda0ForW0[i] = 0;
162 //fhLambda1ForW0[i] = 0;
163 if(i<8)fhMassPairLocMax[i] = 0;
166 //Initialize parameters
171 //_______________________________________________________________________________
172 void AliAnaPi0EbE::FillPileUpHistograms(const Float_t energy, const Float_t time)
174 // Fill some histograms to understand pile-up
175 if(!fFillPileUpHistograms) return;
177 //printf("E %f, time %f\n",energy,time);
178 AliVEvent * event = GetReader()->GetInputEvent();
180 fhTimeENoCut->Fill(energy,time);
181 if(GetReader()->IsPileUpFromSPD()) fhTimeESPD ->Fill(energy,time);
182 if(event->IsPileupFromSPDInMultBins()) fhTimeESPDMulti->Fill(energy,time);
184 if(energy < 8) return; // Fill time figures for high energy clusters not too close to trigger threshold
186 AliESDEvent* esdEv = dynamic_cast<AliESDEvent*> (event);
187 AliAODEvent* aodEv = dynamic_cast<AliAODEvent*> (event);
189 // N pile up vertices
190 Int_t nVerticesSPD = -1;
191 Int_t nVerticesTracks = -1;
195 nVerticesSPD = esdEv->GetNumberOfPileupVerticesSPD();
196 nVerticesTracks = esdEv->GetNumberOfPileupVerticesTracks();
201 nVerticesSPD = aodEv->GetNumberOfPileupVerticesSPD();
202 nVerticesTracks = aodEv->GetNumberOfPileupVerticesTracks();
205 fhTimeNPileUpVertSPD ->Fill(time,nVerticesSPD);
206 fhTimeNPileUpVertTrack->Fill(time,nVerticesTracks);
208 //printf("Is SPD %d, Is SPD Multi %d, n spd %d, n track %d\n",
209 // GetReader()->IsPileUpFromSPD(),event->IsPileupFromSPDInMultBins(),nVerticesSPD,nVerticesTracks);
212 Float_t z1 = -1, z2 = -1;
214 for(Int_t iVert=0; iVert<nVerticesSPD;iVert++)
218 const AliESDVertex* pv=esdEv->GetPileupVertexSPD(iVert);
219 ncont=pv->GetNContributors();
220 z1 = esdEv->GetPrimaryVertexSPD()->GetZ();
222 diamZ = esdEv->GetDiamondZ();
226 AliAODVertex *pv=aodEv->GetVertex(iVert);
227 if(pv->GetType()!=AliAODVertex::kPileupSPD) continue;
228 ncont=pv->GetNContributors();
229 z1=aodEv->GetPrimaryVertexSPD()->GetZ();
231 diamZ = aodEv->GetDiamondZ();
234 Double_t distZ = TMath::Abs(z2-z1);
235 diamZ = TMath::Abs(z2-diamZ);
237 fhTimeNPileUpVertContributors ->Fill(time,ncont);
238 fhTimePileUpMainVertexZDistance->Fill(time,distZ);
239 fhTimePileUpMainVertexZDiamond ->Fill(time,diamZ);
245 //___________________________________________________________________________________________
246 void AliAnaPi0EbE::FillRejectedClusterHistograms(const TLorentzVector mom, const Int_t mctag)
248 // Fill histograms that do not pass the identification (SS case only)
250 Float_t ener = mom.E();
251 Float_t pt = mom.Pt();
252 Float_t phi = mom.Phi();
253 if(phi < 0) phi+=TMath::TwoPi();
254 Float_t eta = mom.Eta();
256 fhPtReject ->Fill(pt);
257 fhEReject ->Fill(ener);
259 fhEEtaReject ->Fill(ener,eta);
260 fhEPhiReject ->Fill(ener,phi);
261 fhEtaPhiReject ->Fill(eta,phi);
265 Int_t mcIndex = GetMCIndex(mctag);
266 fhMCEReject [mcIndex] ->Fill(ener);
267 fhMCPtReject [mcIndex] ->Fill(pt);
271 //_____________________________________________________________________________________
272 void AliAnaPi0EbE::FillSelectedClusterHistograms(AliVCluster* cluster,
277 // Fill shower shape, timing and other histograms for selected clusters from decay
279 Float_t e = cluster->E();
280 Float_t disp = cluster->GetDispersion()*cluster->GetDispersion();
281 Float_t l0 = cluster->GetM02();
282 Float_t l1 = cluster->GetM20();
283 Int_t nSM = GetModuleNumber(cluster);
286 if (e < 2 ) ebin = 0;
287 else if (e < 4 ) ebin = 1;
288 else if (e < 6 ) ebin = 2;
289 else if (e < 10) ebin = 3;
290 else if (e < 15) ebin = 4;
291 else if (e < 20) ebin = 5;
295 if (nMaxima==1) indexMax = 0 ;
296 else if(nMaxima==2) indexMax = 1 ;
300 AliVCaloCells * cell = 0x0;
301 if(fCalorimeter == "PHOS")
302 cell = GetPHOSCells();
304 cell = GetEMCALCells();
306 Float_t maxCellFraction = 0;
307 GetCaloUtils()->GetMaxEnergyCell(cell, cluster, maxCellFraction);
308 fhEFracMaxCell->Fill(e,maxCellFraction);
310 FillWeightHistograms(cluster);
312 fhEDispersion->Fill(e, disp);
313 fhELambda0 ->Fill(e, l0 );
314 fhELambda1 ->Fill(e, l1 );
316 Float_t ll0 = 0., ll1 = 0.;
317 Float_t dispp= 0., dEta = 0., dPhi = 0.;
318 Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
319 if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
321 GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), GetReader()->GetInputEvent()->GetEMCALCells(), cluster,
322 ll0, ll1, dispp, dEta, dPhi, sEta, sPhi, sEtaPhi);
324 fhDispEtaE -> Fill(e,dEta);
325 fhDispPhiE -> Fill(e,dPhi);
326 fhSumEtaE -> Fill(e,sEta);
327 fhSumPhiE -> Fill(e,sPhi);
328 fhSumEtaPhiE -> Fill(e,sEtaPhi);
329 fhDispEtaPhiDiffE -> Fill(e,dPhi-dEta);
330 if(dEta+dPhi>0)fhSphericityE -> Fill(e,(dPhi-dEta)/(dEta+dPhi));
332 fhDispEtaDispPhi[ebin]->Fill(dEta,dPhi);
333 fhLambda0DispEta[ebin]->Fill(l0 ,dEta);
334 fhLambda0DispPhi[ebin]->Fill(l0 ,dPhi);
336 if (fAnaType==kSSCalo)
338 // Asymmetry histograms
339 fhAsymmetryLambda0[ebin]->Fill(l0 ,asy);
340 fhAsymmetryDispEta[ebin]->Fill(dEta,asy);
341 fhAsymmetryDispPhi[ebin]->Fill(dPhi,asy);
345 fhNLocMax->Fill(e,nMaxima);
347 fhELambda0LocMax [indexMax]->Fill(e,l0);
348 fhELambda1LocMax [indexMax]->Fill(e,l1);
349 fhEDispersionLocMax[indexMax]->Fill(e,disp);
351 if(fCalorimeter=="EMCAL" && !fFillOnlySimpleSSHisto)
353 fhEDispEtaLocMax [indexMax]-> Fill(e,dEta);
354 fhEDispPhiLocMax [indexMax]-> Fill(e,dPhi);
355 fhESumEtaPhiLocMax [indexMax]-> Fill(e,sEtaPhi);
356 fhEDispEtaPhiDiffLocMax[indexMax]-> Fill(e,dPhi-dEta);
357 if(dEta+dPhi>0) fhESphericityLocMax[indexMax]->Fill(e,(dPhi-dEta)/(dEta+dPhi));
358 if(fAnaType==kSSCalo) fhEAsymmetryLocMax [indexMax]->Fill(e ,asy);
362 if(fCalorimeter=="EMCAL" && nSM < 6)
364 fhELambda0NoTRD->Fill(e, l0 );
365 fhEFracMaxCellNoTRD->Fill(e,maxCellFraction);
368 if(maxCellFraction < 0.5)
369 fhELambda0FracMaxCellCut->Fill(e, l0 );
371 fhETime ->Fill(e, cluster->GetTOF()*1.e9);
372 fhENCells->Fill(e, cluster->GetNCells());
374 // Fill Track matching control histograms
377 Float_t dZ = cluster->GetTrackDz();
378 Float_t dR = cluster->GetTrackDx();
380 if(cluster->IsEMCAL() && GetCaloUtils()->IsRecalculationOfClusterTrackMatchingOn())
382 dR = 2000., dZ = 2000.;
383 GetCaloUtils()->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
385 //printf("Pi0EbE: dPhi %f, dEta %f\n",dR,dZ);
387 if(fhTrackMatchedDEta && TMath::Abs(dR) < 999)
389 fhTrackMatchedDEta->Fill(e,dZ);
390 fhTrackMatchedDPhi->Fill(e,dR);
391 if(e > 0.5) fhTrackMatchedDEtaDPhi->Fill(dZ,dR);
394 // Check dEdx and E/p of matched clusters
396 if(TMath::Abs(dZ) < 0.05 && TMath::Abs(dR) < 0.05)
398 AliVTrack *track = GetCaloUtils()->GetMatchedTrack(cluster, GetReader()->GetInputEvent());
402 Float_t dEdx = track->GetTPCsignal();
403 fhdEdx->Fill(e, dEdx);
405 Float_t eOverp = e/track->P();
406 fhEOverP->Fill(e, eOverp);
408 if(fCalorimeter=="EMCAL" && nSM < 6) fhEOverPNoTRD->Fill(e, eOverp);
412 // printf("AliAnaPi0EbE::FillSelectedClusterHistograms() - Residual OK but (dR, dZ)= (%2.4f,%2.4f) no track associated WHAT? \n", dR,dZ);
418 if ( !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) )
420 if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) ||
421 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) ) fhTrackMatchedMCParticle->Fill(e, 2.5 );
422 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) ) fhTrackMatchedMCParticle->Fill(e, 0.5 );
423 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron) ) fhTrackMatchedMCParticle->Fill(e, 1.5 );
424 else fhTrackMatchedMCParticle->Fill(e, 3.5 );
429 if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) ||
430 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) ) fhTrackMatchedMCParticle->Fill(e, 6.5 );
431 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) ) fhTrackMatchedMCParticle->Fill(e, 4.5 );
432 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron) ) fhTrackMatchedMCParticle->Fill(e, 5.5 );
433 else fhTrackMatchedMCParticle->Fill(e, 7.5 );
437 }// Track matching histograms
441 Int_t mcIndex = GetMCIndex(tag);
443 fhEMCLambda0[mcIndex] ->Fill(e, l0);
444 fhEMCLambda1[mcIndex] ->Fill(e, l1);
445 fhEMCDispersion[mcIndex] ->Fill(e, disp);
446 fhEMCFracMaxCell[mcIndex]->Fill(e,maxCellFraction);
448 if(fCalorimeter=="EMCAL" && nSM < 6)
449 fhEMCLambda0NoTRD[mcIndex]->Fill(e, l0 );
451 if(maxCellFraction < 0.5)
452 fhEMCLambda0FracMaxCellCut[mcIndex]->Fill(e, l0 );
454 if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
456 fhMCEDispEta [mcIndex]-> Fill(e,dEta);
457 fhMCEDispPhi [mcIndex]-> Fill(e,dPhi);
458 fhMCESumEtaPhi [mcIndex]-> Fill(e,sEtaPhi);
459 fhMCEDispEtaPhiDiff [mcIndex]-> Fill(e,dPhi-dEta);
460 if(dEta+dPhi>0)fhMCESphericity[mcIndex]-> Fill(e,(dPhi-dEta)/(dEta+dPhi));
462 if (fAnaType==kSSCalo)
464 fhMCAsymmetryLambda0[ebin][mcIndex]->Fill(l0 ,asy);
465 fhMCAsymmetryDispEta[ebin][mcIndex]->Fill(dEta,asy);
466 fhMCAsymmetryDispPhi[ebin][mcIndex]->Fill(dPhi,asy);
469 fhMCDispEtaDispPhi[ebin][mcIndex]->Fill(dEta,dPhi);
470 fhMCLambda0DispEta[ebin][mcIndex]->Fill(l0 ,dEta);
471 fhMCLambda0DispPhi[ebin][mcIndex]->Fill(l0 ,dPhi);
479 //________________________________________________________
480 void AliAnaPi0EbE::FillWeightHistograms(AliVCluster *clus)
482 // Calculate weights and fill histograms
484 if(!fFillWeightHistograms || GetMixedEvent()) return;
486 AliVCaloCells* cells = 0;
487 if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
488 else cells = GetPHOSCells();
490 // First recalculate energy in case non linearity was applied
493 for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
496 Int_t id = clus->GetCellsAbsId()[ipos];
498 //Recalibrate cell energy if needed
499 Float_t amp = cells->GetCellAmplitude(id);
500 GetCaloUtils()->RecalibrateCellAmplitude(amp,fCalorimeter, id);
511 printf("AliAnaPi0EbE::WeightHistograms()- Wrong calculated energy %f\n",energy);
515 fhEMaxCellClusterRatio ->Fill(energy,ampMax/energy);
516 fhEMaxCellClusterLogRatio->Fill(energy,TMath::Log(ampMax/energy));
518 //Get the ratio and log ratio to all cells in cluster
519 for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
521 Int_t id = clus->GetCellsAbsId()[ipos];
523 //Recalibrate cell energy if needed
524 Float_t amp = cells->GetCellAmplitude(id);
525 GetCaloUtils()->RecalibrateCellAmplitude(amp,fCalorimeter, id);
527 fhECellClusterRatio ->Fill(energy,amp/energy);
528 fhECellClusterLogRatio->Fill(energy,TMath::Log(amp/energy));
531 //Recalculate shower shape for different W0
532 if(fCalorimeter=="EMCAL"){
534 Float_t l0org = clus->GetM02();
535 Float_t l1org = clus->GetM20();
536 Float_t dorg = clus->GetDispersion();
538 for(Int_t iw = 0; iw < 14; iw++)
540 GetCaloUtils()->GetEMCALRecoUtils()->SetW0(1+iw*0.5);
541 GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), cells, clus);
543 fhLambda0ForW0[iw]->Fill(energy,clus->GetM02());
544 //fhLambda1ForW0[iw]->Fill(energy,clus->GetM20());
548 // Set the original values back
551 clus->SetDispersion(dorg);
556 //__________________________________________
557 TObjString * AliAnaPi0EbE::GetAnalysisCuts()
559 //Save parameters used for analysis
560 TString parList ; //this will be list of parameters used for this analysis.
561 const Int_t buffersize = 255;
562 char onePar[buffersize] ;
564 snprintf(onePar,buffersize,"--- AliAnaPi0EbE ---\n") ;
566 snprintf(onePar,buffersize,"fAnaType=%d (Pi0 selection type) \n",fAnaType) ;
569 if(fAnaType == kSSCalo)
571 snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
573 snprintf(onePar,buffersize,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
575 snprintf(onePar,buffersize,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
577 snprintf(onePar,buffersize,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
581 //Get parameters set in base class.
582 parList += GetBaseParametersList() ;
584 //Get parameters set in PID class.
585 if(fAnaType == kSSCalo) parList += GetCaloPID()->GetPIDParametersList() ;
587 return new TObjString(parList) ;
590 //_____________________________________________
591 TList * AliAnaPi0EbE::GetCreateOutputObjects()
593 // Create histograms to be saved in output file and
594 // store them in outputContainer
595 TList * outputContainer = new TList() ;
596 outputContainer->SetName("Pi0EbEHistos") ;
598 Int_t nptbins = GetHistogramRanges()->GetHistoPtBins(); Float_t ptmax = GetHistogramRanges()->GetHistoPtMax(); Float_t ptmin = GetHistogramRanges()->GetHistoPtMin();
599 Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins(); Float_t phimax = GetHistogramRanges()->GetHistoPhiMax(); Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();
600 Int_t netabins = GetHistogramRanges()->GetHistoEtaBins(); Float_t etamax = GetHistogramRanges()->GetHistoEtaMax(); Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
601 Int_t ssbins = GetHistogramRanges()->GetHistoShowerShapeBins(); Float_t ssmax = GetHistogramRanges()->GetHistoShowerShapeMax(); Float_t ssmin = GetHistogramRanges()->GetHistoShowerShapeMin();
602 Int_t tdbins = GetHistogramRanges()->GetHistoDiffTimeBins() ; Float_t tdmax = GetHistogramRanges()->GetHistoDiffTimeMax(); Float_t tdmin = GetHistogramRanges()->GetHistoDiffTimeMin();
603 Int_t tbins = GetHistogramRanges()->GetHistoTimeBins() ; Float_t tmax = GetHistogramRanges()->GetHistoTimeMax(); Float_t tmin = GetHistogramRanges()->GetHistoTimeMin();
604 Int_t nbins = GetHistogramRanges()->GetHistoNClusterCellBins(); Int_t nmax = GetHistogramRanges()->GetHistoNClusterCellMax(); Int_t nmin = GetHistogramRanges()->GetHistoNClusterCellMin();
606 Int_t nmassbins = GetHistogramRanges()->GetHistoMassBins();
607 Float_t massmin = GetHistogramRanges()->GetHistoMassMin();
608 Float_t massmax = GetHistogramRanges()->GetHistoMassMax();
610 Int_t nresetabins = GetHistogramRanges()->GetHistoTrackResidualEtaBins();
611 Float_t resetamax = GetHistogramRanges()->GetHistoTrackResidualEtaMax();
612 Float_t resetamin = GetHistogramRanges()->GetHistoTrackResidualEtaMin();
613 Int_t nresphibins = GetHistogramRanges()->GetHistoTrackResidualPhiBins();
614 Float_t resphimax = GetHistogramRanges()->GetHistoTrackResidualPhiMax();
615 Float_t resphimin = GetHistogramRanges()->GetHistoTrackResidualPhiMin();
617 Int_t ndedxbins = GetHistogramRanges()->GetHistodEdxBins();
618 Float_t dedxmax = GetHistogramRanges()->GetHistodEdxMax();
619 Float_t dedxmin = GetHistogramRanges()->GetHistodEdxMin();
620 Int_t nPoverEbins = GetHistogramRanges()->GetHistoPOverEBins();
621 Float_t pOverEmax = GetHistogramRanges()->GetHistoPOverEMax();
622 Float_t pOverEmin = GetHistogramRanges()->GetHistoPOverEMin();
624 Int_t ntimebins= GetHistogramRanges()->GetHistoTimeBins();
625 Float_t timemax = GetHistogramRanges()->GetHistoTimeMax();
626 Float_t timemin = GetHistogramRanges()->GetHistoTimeMin();
628 TString nlm[] ={"1 Local Maxima","2 Local Maxima", "NLM > 2"};
629 TString ptype[] ={"#gamma","#gamma->e^{#pm}","#pi^{0}","#eta","e^{#pm}", "hadron"};
630 TString pname[] ={"Photon","Conversion", "Pi0", "Eta", "Electron","Hadron"};
631 Int_t bin[] = {0,2,4,6,10,15,20,100}; // energy bins
633 fhPt = new TH1F("hPt","Number of identified #pi^{0} (#eta) decay",nptbins,ptmin,ptmax);
634 fhPt->SetYTitle("N");
635 fhPt->SetXTitle("p_{T} (GeV/c)");
636 outputContainer->Add(fhPt) ;
638 fhE = new TH1F("hE","Number of identified #pi^{0} (#eta) decay pairs",nptbins,ptmin,ptmax);
640 fhE->SetXTitle("E (GeV)");
641 outputContainer->Add(fhE) ;
644 ("hEPhi","Selected #pi^{0} (#eta) pairs: E vs #phi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
645 fhEPhi->SetYTitle("#phi (rad)");
646 fhEPhi->SetXTitle("E (GeV)");
647 outputContainer->Add(fhEPhi) ;
650 ("hEEta","Selected #pi^{0} (#eta) pairs: E vs #eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
651 fhEEta->SetYTitle("#eta");
652 fhEEta->SetXTitle("E (GeV)");
653 outputContainer->Add(fhEEta) ;
656 ("hEtaPhi","Selected #pi^{0} (#eta) pairs: #eta vs #phi",netabins,etamin,etamax, nphibins,phimin,phimax);
657 fhEtaPhi->SetYTitle("#phi (rad)");
658 fhEtaPhi->SetXTitle("#eta");
659 outputContainer->Add(fhEtaPhi) ;
661 if(fAnaType == kSSCalo)
663 fhPtReject = new TH1F("hPtReject","Number of rejected as #pi^{0} (#eta) decay",nptbins,ptmin,ptmax);
664 fhPtReject->SetYTitle("N");
665 fhPtReject->SetXTitle("p_{T} (GeV/c)");
666 outputContainer->Add(fhPtReject) ;
668 fhEReject = new TH1F("hEReject","Number of rejected as #pi^{0} (#eta) decay pairs",nptbins,ptmin,ptmax);
669 fhEReject->SetYTitle("N");
670 fhEReject->SetXTitle("E (GeV)");
671 outputContainer->Add(fhEReject) ;
673 fhEPhiReject = new TH2F
674 ("hEPhiReject","Rejected #pi^{0} (#eta) cluster: E vs #phi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
675 fhEPhiReject->SetYTitle("#phi (rad)");
676 fhEPhiReject->SetXTitle("E (GeV)");
677 outputContainer->Add(fhEPhiReject) ;
679 fhEEtaReject = new TH2F
680 ("hEEtaReject","Rejected #pi^{0} (#eta) cluster: E vs #eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
681 fhEEtaReject->SetYTitle("#eta");
682 fhEEtaReject->SetXTitle("E (GeV)");
683 outputContainer->Add(fhEEtaReject) ;
685 fhEtaPhiReject = new TH2F
686 ("hEtaPhiReject","Rejected #pi^{0} (#eta) cluster: #eta vs #phi",netabins,etamin,etamax, nphibins,phimin,phimax);
687 fhEtaPhiReject->SetYTitle("#phi (rad)");
688 fhEtaPhiReject->SetXTitle("#eta");
689 outputContainer->Add(fhEtaPhiReject) ;
693 ("hMass","all pairs mass: E vs mass",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
694 fhMass->SetYTitle("mass (GeV/c^{2})");
695 fhMass->SetXTitle("E (GeV)");
696 outputContainer->Add(fhMass) ;
698 fhSelectedMass = new TH2F
699 ("hSelectedMass","Selected #pi^{0} (#eta) pairs mass: E vs mass",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
700 fhSelectedMass->SetYTitle("mass (GeV/c^{2})");
701 fhSelectedMass->SetXTitle("E (GeV)");
702 outputContainer->Add(fhSelectedMass) ;
704 if(fAnaType != kSSCalo)
706 fhPtDecay = new TH1F("hPtDecay","Number of identified #pi^{0} (#eta) decay photons",nptbins,ptmin,ptmax);
707 fhPtDecay->SetYTitle("N");
708 fhPtDecay->SetXTitle("p_{T} (GeV/c)");
709 outputContainer->Add(fhPtDecay) ;
711 fhEDecay = new TH1F("hEDecay","Number of identified #pi^{0} (#eta) decay photons",nptbins,ptmin,ptmax);
712 fhEDecay->SetYTitle("N");
713 fhEDecay->SetXTitle("E (GeV)");
714 outputContainer->Add(fhEDecay) ;
719 if( fFillSelectClHisto )
722 fhEDispersion = new TH2F
723 ("hEDispersion","Selected #pi^{0} (#eta) pairs: E vs dispersion",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
724 fhEDispersion->SetYTitle("D^{2}");
725 fhEDispersion->SetXTitle("E (GeV)");
726 outputContainer->Add(fhEDispersion) ;
728 fhELambda0 = new TH2F
729 ("hELambda0","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
730 fhELambda0->SetYTitle("#lambda_{0}^{2}");
731 fhELambda0->SetXTitle("E (GeV)");
732 outputContainer->Add(fhELambda0) ;
734 fhELambda1 = new TH2F
735 ("hELambda1","Selected #pi^{0} (#eta) pairs: E vs #lambda_{1}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
736 fhELambda1->SetYTitle("#lambda_{1}^{2}");
737 fhELambda1->SetXTitle("E (GeV)");
738 outputContainer->Add(fhELambda1) ;
740 fhELambda0FracMaxCellCut = new TH2F
741 ("hELambda0FracMaxCellCut","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, Max cell fraction of energy < 0.5",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
742 fhELambda0FracMaxCellCut->SetYTitle("#lambda_{0}^{2}");
743 fhELambda0FracMaxCellCut->SetXTitle("E (GeV)");
744 outputContainer->Add(fhELambda0FracMaxCellCut) ;
746 fhEFracMaxCell = new TH2F
747 ("hEFracMaxCell","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, Max cell fraction of energy",nptbins,ptmin,ptmax,100,0,1);
748 fhEFracMaxCell->SetYTitle("Fraction");
749 fhEFracMaxCell->SetXTitle("E (GeV)");
750 outputContainer->Add(fhEFracMaxCell) ;
752 if(fCalorimeter=="EMCAL")
754 fhELambda0NoTRD = new TH2F
755 ("hELambda0NoTRD","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, not behind TRD",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
756 fhELambda0NoTRD->SetYTitle("#lambda_{0}^{2}");
757 fhELambda0NoTRD->SetXTitle("E (GeV)");
758 outputContainer->Add(fhELambda0NoTRD) ;
760 fhEFracMaxCellNoTRD = new TH2F
761 ("hEFracMaxCellNoTRD","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, Max cell fraction of energy, not behind TRD",nptbins,ptmin,ptmax,100,0,1);
762 fhEFracMaxCellNoTRD->SetYTitle("Fraction");
763 fhEFracMaxCellNoTRD->SetXTitle("E (GeV)");
764 outputContainer->Add(fhEFracMaxCellNoTRD) ;
766 if(!fFillOnlySimpleSSHisto)
768 fhDispEtaE = new TH2F ("hDispEtaE","#sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs E", nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
769 fhDispEtaE->SetXTitle("E (GeV)");
770 fhDispEtaE->SetYTitle("#sigma^{2}_{#eta #eta}");
771 outputContainer->Add(fhDispEtaE);
773 fhDispPhiE = new TH2F ("hDispPhiE","#sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs E", nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
774 fhDispPhiE->SetXTitle("E (GeV)");
775 fhDispPhiE->SetYTitle("#sigma^{2}_{#phi #phi}");
776 outputContainer->Add(fhDispPhiE);
778 fhSumEtaE = new TH2F ("hSumEtaE","#sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i})^{2} / #Sigma w_{i} - <#eta>^{2} vs E", nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
779 fhSumEtaE->SetXTitle("E (GeV)");
780 fhSumEtaE->SetYTitle("#delta^{2}_{#eta #eta}");
781 outputContainer->Add(fhSumEtaE);
783 fhSumPhiE = new TH2F ("hSumPhiE","#sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i})^{2}/ #Sigma w_{i} - <#phi>^{2} vs E",
784 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
785 fhSumPhiE->SetXTitle("E (GeV)");
786 fhSumPhiE->SetYTitle("#delta^{2}_{#phi #phi}");
787 outputContainer->Add(fhSumPhiE);
789 fhSumEtaPhiE = new TH2F ("hSumEtaPhiE","#delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs E",
790 nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
791 fhSumEtaPhiE->SetXTitle("E (GeV)");
792 fhSumEtaPhiE->SetYTitle("#delta^{2}_{#eta #phi}");
793 outputContainer->Add(fhSumEtaPhiE);
795 fhDispEtaPhiDiffE = new TH2F ("hDispEtaPhiDiffE","#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs E",
796 nptbins,ptmin,ptmax,200, -10,10);
797 fhDispEtaPhiDiffE->SetXTitle("E (GeV)");
798 fhDispEtaPhiDiffE->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
799 outputContainer->Add(fhDispEtaPhiDiffE);
801 fhSphericityE = new TH2F ("hSphericityE","(#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi}) vs E",
802 nptbins,ptmin,ptmax, 200, -1,1);
803 fhSphericityE->SetXTitle("E (GeV)");
804 fhSphericityE->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
805 outputContainer->Add(fhSphericityE);
807 for(Int_t i = 0; i < 7; i++)
809 fhDispEtaDispPhi[i] = new TH2F (Form("hDispEtaDispPhi_EBin%d",i),Form("#sigma^{2}_{#phi #phi} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",bin[i],bin[i+1]),
810 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
811 fhDispEtaDispPhi[i]->SetXTitle("#sigma^{2}_{#eta #eta}");
812 fhDispEtaDispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
813 outputContainer->Add(fhDispEtaDispPhi[i]);
815 fhLambda0DispEta[i] = new TH2F (Form("hLambda0DispEta_EBin%d",i),Form("#lambda^{2}_{0} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",bin[i],bin[i+1]),
816 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
817 fhLambda0DispEta[i]->SetXTitle("#lambda^{2}_{0}");
818 fhLambda0DispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
819 outputContainer->Add(fhLambda0DispEta[i]);
821 fhLambda0DispPhi[i] = new TH2F (Form("hLambda0DispPhi_EBin%d",i),Form("#lambda^{2}_{0}} vs #sigma^{2}_{#phi #phi} for %d < E < %d GeV",bin[i],bin[i+1]),
822 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
823 fhLambda0DispPhi[i]->SetXTitle("#lambda^{2}_{0}");
824 fhLambda0DispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
825 outputContainer->Add(fhLambda0DispPhi[i]);
831 fhNLocMax = new TH2F("hNLocMax","Number of local maxima in cluster",
832 nptbins,ptmin,ptmax,10,0,10);
833 fhNLocMax ->SetYTitle("N maxima");
834 fhNLocMax ->SetXTitle("E (GeV)");
835 outputContainer->Add(fhNLocMax) ;
837 for (Int_t i = 0; i < 3; i++)
839 fhELambda0LocMax[i] = new TH2F(Form("hELambda0LocMax%d",i+1),
840 Form("Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, %s",nlm[i].Data()),
841 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
842 fhELambda0LocMax[i]->SetYTitle("#lambda_{0}^{2}");
843 fhELambda0LocMax[i]->SetXTitle("E (GeV)");
844 outputContainer->Add(fhELambda0LocMax[i]) ;
846 fhELambda1LocMax[i] = new TH2F(Form("hELambda1LocMax%d",i+1),
847 Form("Selected #pi^{0} (#eta) pairs: E vs #lambda_{1}, %s",nlm[i].Data()),
848 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
849 fhELambda1LocMax[i]->SetYTitle("#lambda_{1}^{2}");
850 fhELambda1LocMax[i]->SetXTitle("E (GeV)");
851 outputContainer->Add(fhELambda1LocMax[i]) ;
853 fhEDispersionLocMax[i] = new TH2F(Form("hEDispersionLocMax%d",i+1),
854 Form("Selected #pi^{0} (#eta) pairs: E vs dispersion^{2}, %s",nlm[i].Data()),
855 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
856 fhEDispersionLocMax[i]->SetYTitle("dispersion^{2}");
857 fhEDispersionLocMax[i]->SetXTitle("E (GeV)");
858 outputContainer->Add(fhEDispersionLocMax[i]) ;
860 if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
862 fhEDispEtaLocMax[i] = new TH2F(Form("hEDispEtaLocMax%d",i+1),
863 Form("Selected #pi^{0} (#eta) pairs: E vs #sigma_{#eta #eta}, %s",nlm[i].Data()),
864 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
865 fhEDispEtaLocMax[i]->SetYTitle("#sigma_{#eta #eta}");
866 fhEDispEtaLocMax[i]->SetXTitle("E (GeV)");
867 outputContainer->Add(fhEDispEtaLocMax[i]) ;
869 fhEDispPhiLocMax[i] = new TH2F(Form("hEDispPhiLocMax%d",i+1),
870 Form("Selected #pi^{0} (#eta) pairs: E vs #sigma_{#phi #phi}, %s",nlm[i].Data()),
871 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
872 fhEDispPhiLocMax[i]->SetYTitle("#sigma_{#phi #phi}");
873 fhEDispPhiLocMax[i]->SetXTitle("E (GeV)");
874 outputContainer->Add(fhEDispPhiLocMax[i]) ;
876 fhESumEtaPhiLocMax[i] = new TH2F(Form("hESumEtaPhiLocMax%d",i+1),
877 Form("Selected #pi^{0} (#eta) pairs: E vs #sigma_{#eta #phi}, %s",nlm[i].Data()),
878 nptbins,ptmin,ptmax,2*ssbins,-ssmax,ssmax);
879 fhESumEtaPhiLocMax[i]->SetYTitle("#sigma_{#eta #phi}");
880 fhESumEtaPhiLocMax[i]->SetXTitle("E (GeV)");
881 outputContainer->Add(fhESumEtaPhiLocMax[i]) ;
883 fhEDispEtaPhiDiffLocMax[i] = new TH2F(Form("hEDispEtaPhiDiffLocMax%d",i+1),
884 Form("Selected #pi^{0} (#eta) pairs: E vs #sigma_{#phi #phi} - #sigma_{#eta #eta}, %s",nlm[i].Data()),
885 nptbins,ptmin,ptmax,200, -10,10);
886 fhEDispEtaPhiDiffLocMax[i]->SetYTitle("#sigma_{#phi #phi} - #sigma_{#eta #eta}");
887 fhEDispEtaPhiDiffLocMax[i]->SetXTitle("E (GeV)");
888 outputContainer->Add(fhEDispEtaPhiDiffLocMax[i]) ;
890 fhESphericityLocMax[i] = new TH2F(Form("hESphericityLocMax%d",i+1),
891 Form("Selected #pi^{0} (#eta) pairs: E vs #sigma_{#phi #phi} - #sigma_{#eta #eta} / (#sigma_{#phi #phi} + #sigma_{#eta #eta}), %s",nlm[i].Data()),
892 nptbins,ptmin,ptmax,200, -1,1);
893 fhESphericityLocMax[i]->SetYTitle("#sigma_{#phi #phi} - #sigma_{#eta #eta} / (#sigma_{#phi #phi} + #sigma_{#eta #eta})");
894 fhESphericityLocMax[i]->SetXTitle("E (GeV)");
895 outputContainer->Add(fhESphericityLocMax[i]) ;
900 fhENCells = new TH2F ("hENCells","N cells in cluster vs E ", nptbins,ptmin,ptmax, nbins,nmin,nmax);
901 fhENCells->SetXTitle("E (GeV)");
902 fhENCells->SetYTitle("# of cells in cluster");
903 outputContainer->Add(fhENCells);
905 fhETime = new TH2F("hETime","cluster time vs pair E",nptbins,ptmin,ptmax, tbins,tmin,tmax);
906 fhETime->SetXTitle("E (GeV)");
907 fhETime->SetYTitle("t (ns)");
908 outputContainer->Add(fhETime);
912 if(fAnaType == kIMCalo)
914 fhEPairDiffTime = new TH2F("hEPairDiffTime","cluster pair time difference vs E",nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
915 fhEPairDiffTime->SetXTitle("E_{pair} (GeV)");
916 fhEPairDiffTime->SetYTitle("#Delta t (ns)");
917 outputContainer->Add(fhEPairDiffTime);
919 TString combiName [] = {"1LocMax","2LocMax","NLocMax","1LocMax2LocMax","1LocMaxNLocMax","2LocMaxNLocMax","1LocMaxSSBad","NLocMaxSSGood"};
920 TString combiTitle[] = {"1 Local Maxima in both clusters","2 Local Maxima in both clusters","more than 2 Local Maxima in both clusters",
921 "1 Local Maxima paired with 2 Local Maxima","1 Local Maxima paired with more than 2 Local Maxima",
922 "2 Local Maxima paired with more than 2 Local Maxima",
923 "1 Local Maxima paired with #lambda_{0}^{2}>0.3","N Local Maxima paired with 0.1<#lambda_{0}^{2}<0.3"};
925 for (Int_t i = 0; i < 8 ; i++)
928 if (fAnaType == kIMCaloTracks && i > 2 ) continue ;
930 fhMassPairLocMax[i] = new TH2F
931 (Form("MassPairLocMax%s",combiName[i].Data()),
932 Form("Mass for decay #gamma pair vs E_{pair}, origin #pi^{0}, %s", combiTitle[i].Data()),
933 nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
934 fhMassPairLocMax[i]->SetYTitle("Mass (MeV/c^{2})");
935 fhMassPairLocMax[i]->SetXTitle("E_{pair} (GeV)");
936 outputContainer->Add(fhMassPairLocMax[i]) ;
942 fhTrackMatchedDEta = new TH2F
943 ("hTrackMatchedDEta",
944 "d#eta of cluster-track vs cluster energy",
945 nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
946 fhTrackMatchedDEta->SetYTitle("d#eta");
947 fhTrackMatchedDEta->SetXTitle("E_{cluster} (GeV)");
949 fhTrackMatchedDPhi = new TH2F
950 ("hTrackMatchedDPhi",
951 "d#phi of cluster-track vs cluster energy",
952 nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
953 fhTrackMatchedDPhi->SetYTitle("d#phi (rad)");
954 fhTrackMatchedDPhi->SetXTitle("E_{cluster} (GeV)");
956 fhTrackMatchedDEtaDPhi = new TH2F
957 ("hTrackMatchedDEtaDPhi",
958 "d#eta vs d#phi of cluster-track vs cluster energy",
959 nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
960 fhTrackMatchedDEtaDPhi->SetYTitle("d#phi (rad)");
961 fhTrackMatchedDEtaDPhi->SetXTitle("d#eta");
963 outputContainer->Add(fhTrackMatchedDEta) ;
964 outputContainer->Add(fhTrackMatchedDPhi) ;
965 outputContainer->Add(fhTrackMatchedDEtaDPhi) ;
967 fhdEdx = new TH2F ("hdEdx","matched track <dE/dx> vs cluster E ", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
968 fhdEdx->SetXTitle("E (GeV)");
969 fhdEdx->SetYTitle("<dE/dx>");
970 outputContainer->Add(fhdEdx);
972 fhEOverP = new TH2F ("hEOverP","matched track E/p vs cluster E ", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
973 fhEOverP->SetXTitle("E (GeV)");
974 fhEOverP->SetYTitle("E/p");
975 outputContainer->Add(fhEOverP);
977 if(fCalorimeter=="EMCAL")
979 fhEOverPNoTRD = new TH2F ("hEOverPNoTRD","matched track E/p vs cluster E, SM not behind TRD ", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
980 fhEOverPNoTRD->SetXTitle("E (GeV)");
981 fhEOverPNoTRD->SetYTitle("E/p");
982 outputContainer->Add(fhEOverPNoTRD);
985 if(IsDataMC() && fFillTMHisto)
987 fhTrackMatchedMCParticle = new TH2F
988 ("hTrackMatchedMCParticle",
989 "Origin of particle vs energy",
990 nptbins,ptmin,ptmax,8,0,8);
991 fhTrackMatchedMCParticle->SetXTitle("E (GeV)");
992 //fhTrackMatchedMCParticle->SetYTitle("Particle type");
994 fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(1 ,"Photon");
995 fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(2 ,"Electron");
996 fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
997 fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(4 ,"Rest");
998 fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
999 fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
1000 fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
1001 fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
1003 outputContainer->Add(fhTrackMatchedMCParticle);
1007 if(fFillWeightHistograms)
1009 fhECellClusterRatio = new TH2F ("hECellClusterRatio"," cell energy / cluster energy vs cluster energy, for selected decay photons from neutral meson",
1010 nptbins,ptmin,ptmax, 100,0,1.);
1011 fhECellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
1012 fhECellClusterRatio->SetYTitle("E_{cell i}/E_{cluster}");
1013 outputContainer->Add(fhECellClusterRatio);
1015 fhECellClusterLogRatio = new TH2F ("hECellClusterLogRatio"," Log(cell energy / cluster energy) vs cluster energy, for selected decay photons from neutral meson",
1016 nptbins,ptmin,ptmax, 100,-10,0);
1017 fhECellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
1018 fhECellClusterLogRatio->SetYTitle("Log (E_{max cell}/E_{cluster})");
1019 outputContainer->Add(fhECellClusterLogRatio);
1021 fhEMaxCellClusterRatio = new TH2F ("hEMaxCellClusterRatio"," max cell energy / cluster energy vs cluster energy, for selected decay photons from neutral meson",
1022 nptbins,ptmin,ptmax, 100,0,1.);
1023 fhEMaxCellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
1024 fhEMaxCellClusterRatio->SetYTitle("E_{max cell}/E_{cluster}");
1025 outputContainer->Add(fhEMaxCellClusterRatio);
1027 fhEMaxCellClusterLogRatio = new TH2F ("hEMaxCellClusterLogRatio"," Log(max cell energy / cluster energy) vs cluster energy, for selected decay photons from neutral meson",
1028 nptbins,ptmin,ptmax, 100,-10,0);
1029 fhEMaxCellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
1030 fhEMaxCellClusterLogRatio->SetYTitle("Log (E_{max cell}/E_{cluster})");
1031 outputContainer->Add(fhEMaxCellClusterLogRatio);
1033 for(Int_t iw = 0; iw < 14; iw++)
1035 fhLambda0ForW0[iw] = new TH2F (Form("hLambda0ForW0%d",iw),Form("shower shape, #lambda^{2}_{0} vs E, w0 = %1.1f, for selected decay photons from neutral meson",1+0.5*iw),
1036 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1037 fhLambda0ForW0[iw]->SetXTitle("E_{cluster}");
1038 fhLambda0ForW0[iw]->SetYTitle("#lambda^{2}_{0}");
1039 outputContainer->Add(fhLambda0ForW0[iw]);
1041 // fhLambda1ForW0[iw] = new TH2F (Form("hLambda1ForW0%d",iw),Form("shower shape, #lambda^{2}_{1} vs E, w0 = %1.1f, for selected decay photons from neutral meson",0.5+0.5*iw),
1042 // nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1043 // fhLambda1ForW0[iw]->SetXTitle("E_{cluster}");
1044 // fhLambda1ForW0[iw]->SetYTitle("#lambda^{2}_{1}");
1045 // outputContainer->Add(fhLambda1ForW0[iw]);
1052 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC && fAnaType==kSSCalo)
1054 fhMCPi0PtGenRecoFraction = new TH2F("hMCPi0PtGenRecoFraction","Number of clusters from #pi^{0} (2 #gamma) identified as #pi^{0} (#eta), pT versus E primary #pi^{0} / E reco",
1055 nptbins,ptmin,ptmax,200,0,2);
1056 fhMCPi0PtGenRecoFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
1057 fhMCPi0PtGenRecoFraction->SetYTitle("E^{ #pi^{0} mother} / E^{rec}");
1058 outputContainer->Add(fhMCPi0PtGenRecoFraction) ;
1060 fhMCEtaPtGenRecoFraction = new TH2F("hMCEtaPtGenRecoFraction","Number of clusters from #eta (2 #gamma) identified as #pi^{0} (#eta),pT versus E primary #eta / E reco",
1061 nptbins,ptmin,ptmax,200,0,2);
1062 fhMCEtaPtGenRecoFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
1063 fhMCEtaPtGenRecoFraction->SetYTitle("E^{ #eta mother} / E^{rec}");
1064 outputContainer->Add(fhMCEtaPtGenRecoFraction) ;
1066 fhMCPi0DecayPt = new TH1F("hMCPi0DecayPt","Number of #gamma from #pi^{0} decay identified as #pi^{0} (#eta)",nptbins,ptmin,ptmax);
1067 fhMCPi0DecayPt->SetYTitle("N");
1068 fhMCPi0DecayPt->SetXTitle("p^{rec}_{T} (GeV/c)");
1069 outputContainer->Add(fhMCPi0DecayPt) ;
1071 fhMCPi0DecayPtFraction = new TH2F("hMCPi0DecayPtFraction","Number of #gamma from #pi^{0} decay identified as #pi^{0} (#eta), pT versus E primary #gamma / E primary #pi^{0}",
1072 nptbins,ptmin,ptmax,100,0,1);
1073 fhMCPi0DecayPtFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
1074 fhMCPi0DecayPtFraction->SetYTitle("E^{gen} / E^{gen-mother}");
1075 outputContainer->Add(fhMCPi0DecayPtFraction) ;
1077 fhMCEtaDecayPt = new TH1F("hMCEtaDecayPt","Number of #gamma from #eta decay identified as #pi^{0} (#eta)",nptbins,ptmin,ptmax);
1078 fhMCEtaDecayPt->SetYTitle("N");
1079 fhMCEtaDecayPt->SetXTitle("p^{rec}_{T} (GeV/c)");
1080 outputContainer->Add(fhMCEtaDecayPt) ;
1082 fhMCEtaDecayPtFraction = new TH2F("hMCEtaDecayPtFraction","Number of #gamma from #eta decay identified as #pi^{0} (#eta), pT versus E primary #gamma / E primary #eta",
1083 nptbins,ptmin,ptmax,100,0,1);
1084 fhMCEtaDecayPtFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
1085 fhMCEtaDecayPtFraction->SetYTitle("E^{gen} / E^{gen-mother}");
1086 outputContainer->Add(fhMCEtaDecayPtFraction) ;
1088 fhMCOtherDecayPt = new TH1F("hMCOtherDecayPt","Number of #gamma decay (not #eta or #pi^{0}) identified as #pi^{0} (#eta)",nptbins,ptmin,ptmax);
1089 fhMCOtherDecayPt->SetYTitle("N");
1090 fhMCOtherDecayPt->SetXTitle("p^{rec}_{T} (GeV/c)");
1091 outputContainer->Add(fhMCOtherDecayPt) ;
1095 if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) ||
1096 GetReader()->GetDataType() != AliCaloTrackReader::kMC)
1099 fhAnglePairMCPi0 = new TH2F
1101 "Angle between decay #gamma pair vs E_{pair}, origin #pi^{0}",nptbins,ptmin,ptmax,250,0,0.5);
1102 fhAnglePairMCPi0->SetYTitle("#alpha (rad)");
1103 fhAnglePairMCPi0->SetXTitle("E_{pair} (GeV)");
1104 outputContainer->Add(fhAnglePairMCPi0) ;
1106 if (fAnaType!= kSSCalo)
1108 fhAnglePairMCEta = new TH2F
1110 "Angle between decay #gamma pair vs E_{pair}, origin #eta",nptbins,ptmin,ptmax,250,0,0.5);
1111 fhAnglePairMCEta->SetYTitle("#alpha (rad)");
1112 fhAnglePairMCEta->SetXTitle("E_{pair} (GeV)");
1113 outputContainer->Add(fhAnglePairMCEta) ;
1115 fhMassPairMCPi0 = new TH2F
1117 "Mass for decay #gamma pair vs E_{pair}, origin #pi^{0}",nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1118 fhMassPairMCPi0->SetYTitle("Mass (MeV/c^{2})");
1119 fhMassPairMCPi0->SetXTitle("E_{pair} (GeV)");
1120 outputContainer->Add(fhMassPairMCPi0) ;
1122 fhMassPairMCEta = new TH2F
1124 "Mass for decay #gamma pair vs E_{pair}, origin #eta",nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1125 fhMassPairMCEta->SetYTitle("Mass (MeV/c^{2})");
1126 fhMassPairMCEta->SetXTitle("E_{pair} (GeV)");
1127 outputContainer->Add(fhMassPairMCEta) ;
1130 for(Int_t i = 0; i < 6; i++)
1134 (Form("hE_MC%s",pname[i].Data()),
1135 Form("Identified as #pi^{0} (#eta), cluster from %s",
1137 nptbins,ptmin,ptmax);
1138 fhMCE[i]->SetYTitle("N");
1139 fhMCE[i]->SetXTitle("E (GeV)");
1140 outputContainer->Add(fhMCE[i]) ;
1142 fhMCPt[i] = new TH1F
1143 (Form("hPt_MC%s",pname[i].Data()),
1144 Form("Identified as #pi^{0} (#eta), cluster from %s",
1146 nptbins,ptmin,ptmax);
1147 fhMCPt[i]->SetYTitle("N");
1148 fhMCPt[i]->SetXTitle("p_{T} (GeV/c)");
1149 outputContainer->Add(fhMCPt[i]) ;
1151 if(fAnaType == kSSCalo)
1153 fhMCEReject[i] = new TH1F
1154 (Form("hEReject_MC%s",pname[i].Data()),
1155 Form("Rejected as #pi^{0} (#eta), cluster from %s",
1157 nptbins,ptmin,ptmax);
1158 fhMCEReject[i]->SetYTitle("N");
1159 fhMCEReject[i]->SetXTitle("E (GeV)");
1160 outputContainer->Add(fhMCEReject[i]) ;
1162 fhMCPtReject[i] = new TH1F
1163 (Form("hPtReject_MC%s",pname[i].Data()),
1164 Form("Rejected as #pi^{0} (#eta), cluster from %s",
1166 nptbins,ptmin,ptmax);
1167 fhMCPtReject[i]->SetYTitle("N");
1168 fhMCPtReject[i]->SetXTitle("p_{T} (GeV/c)");
1169 outputContainer->Add(fhMCPtReject[i]) ;
1172 fhMCPhi[i] = new TH2F
1173 (Form("hPhi_MC%s",pname[i].Data()),
1174 Form("Identified as #pi^{0} (#eta), cluster from %s",ptype[i].Data()),
1175 nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1176 fhMCPhi[i]->SetYTitle("#phi");
1177 fhMCPhi[i]->SetXTitle("p_{T} (GeV/c)");
1178 outputContainer->Add(fhMCPhi[i]) ;
1180 fhMCEta[i] = new TH2F
1181 (Form("hEta_MC%s",pname[i].Data()),
1182 Form("Identified as #pi^{0} (#eta), cluster from %s",
1183 ptype[i].Data()),nptbins,ptmin,ptmax,netabins,etamin,etamax);
1184 fhMCEta[i]->SetYTitle("#eta");
1185 fhMCEta[i]->SetXTitle("p_{T} (GeV/c)");
1186 outputContainer->Add(fhMCEta[i]) ;
1189 if( fFillSelectClHisto )
1191 fhEMCLambda0[i] = new TH2F(Form("hELambda0_MC%s",pname[i].Data()),
1192 Form("Selected pair, cluster from %s : E vs #lambda_{0}^{2}",ptype[i].Data()),
1193 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1194 fhEMCLambda0[i]->SetYTitle("#lambda_{0}^{2}");
1195 fhEMCLambda0[i]->SetXTitle("E (GeV)");
1196 outputContainer->Add(fhEMCLambda0[i]) ;
1198 fhEMCLambda1[i] = new TH2F(Form("hELambda1_MC%s",pname[i].Data()),
1199 Form("Selected pair, cluster from %s : E vs #lambda_{1}^{2}",ptype[i].Data()),
1200 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1201 fhEMCLambda1[i]->SetYTitle("#lambda_{1}^{2}");
1202 fhEMCLambda1[i]->SetXTitle("E (GeV)");
1203 outputContainer->Add(fhEMCLambda1[i]) ;
1205 fhEMCDispersion[i] = new TH2F(Form("hEDispersion_MC%s",pname[i].Data()),
1206 Form("Selected pair, cluster from %s : E vs dispersion^{2}",ptype[i].Data()),
1207 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1208 fhEMCDispersion[i]->SetYTitle("D^{2}");
1209 fhEMCDispersion[i]->SetXTitle("E (GeV)");
1210 outputContainer->Add(fhEMCDispersion[i]) ;
1212 if(fCalorimeter=="EMCAL")
1214 fhEMCLambda0NoTRD[i] = new TH2F(Form("hELambda0NoTRD_MC%s",pname[i].Data()),
1215 Form("Selected pair, cluster from %s : E vs #lambda_{0}^{2}, NoTRD",ptype[i].Data()),
1216 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1217 fhEMCLambda0NoTRD[i]->SetYTitle("#lambda_{0}^{2}");
1218 fhEMCLambda0NoTRD[i]->SetXTitle("E (GeV)");
1219 outputContainer->Add(fhEMCLambda0NoTRD[i]) ;
1221 if(!fFillOnlySimpleSSHisto)
1223 fhMCEDispEta[i] = new TH2F (Form("hEDispEtaE_MC%s",pname[i].Data()),
1224 Form("cluster from %s : #sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs E",ptype[i].Data()),
1225 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1226 fhMCEDispEta[i]->SetXTitle("E (GeV)");
1227 fhMCEDispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
1228 outputContainer->Add(fhMCEDispEta[i]);
1230 fhMCEDispPhi[i] = new TH2F (Form("hEDispPhiE_MC%s",pname[i].Data()),
1231 Form("cluster from %s : #sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs E",ptype[i].Data()),
1232 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1233 fhMCEDispPhi[i]->SetXTitle("E (GeV)");
1234 fhMCEDispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1235 outputContainer->Add(fhMCEDispPhi[i]);
1237 fhMCESumEtaPhi[i] = new TH2F (Form("hESumEtaPhiE_MC%s",pname[i].Data()),
1238 Form("cluster from %s : #delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs E",ptype[i].Data()),
1239 nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
1240 fhMCESumEtaPhi[i]->SetXTitle("E (GeV)");
1241 fhMCESumEtaPhi[i]->SetYTitle("#delta^{2}_{#eta #phi}");
1242 outputContainer->Add(fhMCESumEtaPhi[i]);
1244 fhMCEDispEtaPhiDiff[i] = new TH2F (Form("hEDispEtaPhiDiffE_MC%s",pname[i].Data()),
1245 Form("cluster from %s : #sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs E",ptype[i].Data()),
1246 nptbins,ptmin,ptmax,200,-10,10);
1247 fhMCEDispEtaPhiDiff[i]->SetXTitle("E (GeV)");
1248 fhMCEDispEtaPhiDiff[i]->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
1249 outputContainer->Add(fhMCEDispEtaPhiDiff[i]);
1251 fhMCESphericity[i] = new TH2F (Form("hESphericity_MC%s",pname[i].Data()),
1252 Form("cluster from %s : (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi}) vs E",ptype[i].Data()),
1253 nptbins,ptmin,ptmax, 200,-1,1);
1254 fhMCESphericity[i]->SetXTitle("E (GeV)");
1255 fhMCESphericity[i]->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
1256 outputContainer->Add(fhMCESphericity[i]);
1258 for(Int_t ie = 0; ie < 7; ie++)
1260 fhMCDispEtaDispPhi[ie][i] = new TH2F (Form("hMCDispEtaDispPhi_EBin%d_MC%s",ie,pname[i].Data()),
1261 Form("cluster from %s : #sigma^{2}_{#phi #phi} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
1262 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1263 fhMCDispEtaDispPhi[ie][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
1264 fhMCDispEtaDispPhi[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1265 outputContainer->Add(fhMCDispEtaDispPhi[ie][i]);
1267 fhMCLambda0DispEta[ie][i] = new TH2F (Form("hMCLambda0DispEta_EBin%d_MC%s",ie,pname[i].Data()),
1268 Form("cluster from %s : #lambda^{2}_{0} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
1269 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1270 fhMCLambda0DispEta[ie][i]->SetXTitle("#lambda^{2}_{0}");
1271 fhMCLambda0DispEta[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1272 outputContainer->Add(fhMCLambda0DispEta[ie][i]);
1274 fhMCLambda0DispPhi[ie][i] = new TH2F (Form("hMCLambda0DispPhi_EBin%d_MC%s",ie,pname[i].Data()),
1275 Form("cluster from %s :#lambda^{2}_{0} vs #sigma^{2}_{#phi #phi} for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
1276 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1277 fhMCLambda0DispPhi[ie][i]->SetXTitle("#lambda^{2}_{0}");
1278 fhMCLambda0DispPhi[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1279 outputContainer->Add(fhMCLambda0DispPhi[ie][i]);
1285 fhEMCLambda0FracMaxCellCut[i] = new TH2F(Form("hELambda0FracMaxCellCut_MC%s",pname[i].Data()),
1286 Form("Selected pair, cluster from %s : E vs #lambda_{0}^{2}, Max cell fraction of energy < 0.5 ",ptype[i].Data()),
1287 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1288 fhEMCLambda0FracMaxCellCut[i]->SetYTitle("#lambda_{0}^{2}");
1289 fhEMCLambda0FracMaxCellCut[i]->SetXTitle("E (GeV)");
1290 outputContainer->Add(fhEMCLambda0FracMaxCellCut[i]) ;
1292 fhEMCFracMaxCell[i] = new TH2F(Form("hEFracMaxCell_MC%s",pname[i].Data()),
1293 Form("Selected pair, cluster from %s : E vs Max cell fraction of energy",ptype[i].Data()),
1294 nptbins,ptmin,ptmax,100,0,1);
1295 fhEMCFracMaxCell[i]->SetYTitle("Fraction");
1296 fhEMCFracMaxCell[i]->SetXTitle("E (GeV)");
1297 outputContainer->Add(fhEMCFracMaxCell[i]) ;
1300 } // shower shape histo
1305 if(fAnaType==kSSCalo)
1307 fhAsymmetry = new TH2F ("hAsymmetry","A = ( E1 - E2 ) / ( E1 + E2 ) vs E",
1308 nptbins,ptmin,ptmax, 200, -1,1);
1309 fhAsymmetry->SetXTitle("E (GeV)");
1310 fhAsymmetry->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1311 outputContainer->Add(fhAsymmetry);
1313 fhSelectedAsymmetry = new TH2F ("hSelectedAsymmetry","A = ( E1 - E2 ) / ( E1 + E2 ) vs E",
1314 nptbins,ptmin,ptmax, 200, -1,1);
1315 fhSelectedAsymmetry->SetXTitle("E (GeV)");
1316 fhSelectedAsymmetry->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1317 outputContainer->Add(fhSelectedAsymmetry);
1321 for(Int_t i = 0; i< 6; i++)
1323 fhMCEAsymmetry[i] = new TH2F (Form("hEAsymmetry_MC%s",pname[i].Data()),
1324 Form("cluster from %s : A = ( E1 - E2 ) / ( E1 + E2 ) vs E",ptype[i].Data()),
1325 nptbins,ptmin,ptmax, 200,-1,1);
1326 fhMCEAsymmetry[i]->SetXTitle("E (GeV)");
1327 fhMCEAsymmetry[i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1328 outputContainer->Add(fhMCEAsymmetry[i]);
1333 if(fAnaType==kSSCalo && fFillSelectClHisto && !fFillOnlySimpleSSHisto )
1337 for(Int_t i = 0; i< 3; i++)
1339 fhEAsymmetryLocMax[i] = new TH2F(Form("hEAsymmetryLocMax%d",i+1),
1340 Form("Selected #pi^{0} (#eta) pairs: E vs A = ( E1 - E2 ) / ( E1 + E2 ), %s",nlm[i].Data()),
1341 nptbins,ptmin,ptmax,200, -1,1);
1342 fhEAsymmetryLocMax[i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1343 fhEAsymmetryLocMax[i]->SetXTitle("E (GeV)");
1344 outputContainer->Add(fhEAsymmetryLocMax[i]) ;
1347 for(Int_t ie = 0; ie< 7; ie++)
1350 fhAsymmetryLambda0[ie] = new TH2F (Form("hAsymmetryLambda0_EBin%d",ie),
1351 Form("#lambda_{0}^{2} vs A for %d < E < %d GeV",bin[ie],bin[ie+1]),
1352 ssbins,ssmin,ssmax , 200,-1,1);
1353 fhAsymmetryLambda0[ie]->SetXTitle("#lambda_{0}^{2}");
1354 fhAsymmetryLambda0[ie]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1355 outputContainer->Add(fhAsymmetryLambda0[ie]);
1357 fhAsymmetryDispEta[ie] = new TH2F (Form("hAsymmetryDispEta_EBin%d",ie),
1358 Form("#sigma^{2}_{#eta #eta} vs A for %d < E < %d GeV",bin[ie],bin[ie+1]),
1359 ssbins,ssmin,ssmax , 200,-1,1);
1360 fhAsymmetryDispEta[ie]->SetXTitle("#sigma^{2}_{#eta #eta}");
1361 fhAsymmetryDispEta[ie]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1362 outputContainer->Add(fhAsymmetryDispEta[ie]);
1364 fhAsymmetryDispPhi[ie] = new TH2F (Form("hAsymmetryDispPhi_EBin%d",ie),
1365 Form("#sigma^{2}_{#phi #phi} vs A for %d < E < %d GeV",bin[ie],bin[ie+1]),
1366 ssbins,ssmin,ssmax , 200,-1,1);
1367 fhAsymmetryDispPhi[ie]->SetXTitle("#sigma^{2}_{#phi #phi}");
1368 fhAsymmetryDispPhi[ie]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1369 outputContainer->Add(fhAsymmetryDispPhi[ie]);
1375 for(Int_t i = 0; i< 6; i++)
1377 for(Int_t ie = 0; ie < 7; ie++)
1379 fhMCAsymmetryLambda0[ie][i] = new TH2F (Form("hMCAsymmetryLambda0_EBin%d_MC%s",ie,pname[i].Data()),
1380 Form("cluster from %s : #lambda_{0}^{2} vs A for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
1381 ssbins,ssmin,ssmax , 200,-1,1);
1382 fhMCAsymmetryLambda0[ie][i]->SetXTitle("#lambda_{0}^{2}");
1383 fhMCAsymmetryLambda0[ie][i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1384 outputContainer->Add(fhMCAsymmetryLambda0[ie][i]);
1386 fhMCAsymmetryDispEta[ie][i] = new TH2F (Form("hMCAsymmetryDispEta_EBin%d_MC%s",ie,pname[i].Data()),
1387 Form("cluster from %s : #sigma^{2}_{#eta #eta} vs A for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
1388 ssbins,ssmin,ssmax , 200,-1,1);
1389 fhMCAsymmetryDispEta[ie][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
1390 fhMCAsymmetryDispEta[ie][i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1391 outputContainer->Add(fhMCAsymmetryDispEta[ie][i]);
1393 fhMCAsymmetryDispPhi[ie][i] = new TH2F (Form("hMCAsymmetryDispPhi_EBin%d_MC%s",ie,pname[i].Data()),
1394 Form("cluster from %s : #sigma^{2}_{#phi #phi} vs A for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
1395 ssbins,ssmin,ssmax , 200,-1,1);
1396 fhMCAsymmetryDispPhi[ie][i]->SetXTitle("#sigma^{2}_{#phi #phi}");
1397 fhMCAsymmetryDispPhi[ie][i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1398 outputContainer->Add(fhMCAsymmetryDispPhi[ie][i]);
1404 if(fFillPileUpHistograms)
1407 TString pileUpName[] = {"SPD","EMCAL","SPDOrEMCAL","SPDAndEMCAL","SPDAndNotEMCAL","EMCALAndNotSPD","NotSPDAndNotEMCAL"} ;
1409 for(Int_t i = 0 ; i < 7 ; i++)
1411 fhPtPi0PileUp[i] = new TH1F(Form("hPtPi0PileUp%s",pileUpName[i].Data()),
1412 Form("Selected #pi^{0} (#eta) p_{T} distribution, %s Pile-Up event",pileUpName[i].Data()), nptbins,ptmin,ptmax);
1413 fhPtPi0PileUp[i]->SetXTitle("p_{T} (GeV/c)");
1414 outputContainer->Add(fhPtPi0PileUp[i]);
1417 fhTimeENoCut = new TH2F ("hTimeE_NoCut","time of cluster vs E of clusters, no cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
1418 fhTimeENoCut->SetXTitle("E (GeV)");
1419 fhTimeENoCut->SetYTitle("time (ns)");
1420 outputContainer->Add(fhTimeENoCut);
1422 fhTimeESPD = new TH2F ("hTimeE_SPD","time of cluster vs E of clusters, SPD cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
1423 fhTimeESPD->SetXTitle("E (GeV)");
1424 fhTimeESPD->SetYTitle("time (ns)");
1425 outputContainer->Add(fhTimeESPD);
1427 fhTimeESPDMulti = new TH2F ("hTimeE_SPDMulti","time of cluster vs E of clusters, SPD multi cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
1428 fhTimeESPDMulti->SetXTitle("E (GeV)");
1429 fhTimeESPDMulti->SetYTitle("time (ns)");
1430 outputContainer->Add(fhTimeESPDMulti);
1432 fhTimeNPileUpVertSPD = new TH2F ("hTime_NPileUpVertSPD","time of cluster vs N pile-up SPD vertex", ntimebins,timemin,timemax,50,0,50);
1433 fhTimeNPileUpVertSPD->SetYTitle("# vertex ");
1434 fhTimeNPileUpVertSPD->SetXTitle("time (ns)");
1435 outputContainer->Add(fhTimeNPileUpVertSPD);
1437 fhTimeNPileUpVertTrack = new TH2F ("hTime_NPileUpVertTracks","time of cluster vs N pile-up Tracks vertex", ntimebins,timemin,timemax, 50,0,50 );
1438 fhTimeNPileUpVertTrack->SetYTitle("# vertex ");
1439 fhTimeNPileUpVertTrack->SetXTitle("time (ns)");
1440 outputContainer->Add(fhTimeNPileUpVertTrack);
1442 fhTimeNPileUpVertContributors = new TH2F ("hTime_NPileUpVertContributors","time of cluster vs N constributors to pile-up SPD vertex", ntimebins,timemin,timemax,50,0,50);
1443 fhTimeNPileUpVertContributors->SetYTitle("# vertex ");
1444 fhTimeNPileUpVertContributors->SetXTitle("time (ns)");
1445 outputContainer->Add(fhTimeNPileUpVertContributors);
1447 fhTimePileUpMainVertexZDistance = new TH2F ("hTime_PileUpMainVertexZDistance","time of cluster vs distance in Z pile-up SPD vertex - main SPD vertex",ntimebins,timemin,timemax,100,0,50);
1448 fhTimePileUpMainVertexZDistance->SetYTitle("distance Z (cm) ");
1449 fhTimePileUpMainVertexZDistance->SetXTitle("time (ns)");
1450 outputContainer->Add(fhTimePileUpMainVertexZDistance);
1452 fhTimePileUpMainVertexZDiamond = new TH2F ("hTime_PileUpMainVertexZDiamond","time of cluster vs distance in Z pile-up SPD vertex - z diamond",ntimebins,timemin,timemax,100,0,50);
1453 fhTimePileUpMainVertexZDiamond->SetYTitle("diamond distance Z (cm) ");
1454 fhTimePileUpMainVertexZDiamond->SetXTitle("time (ns)");
1455 outputContainer->Add(fhTimePileUpMainVertexZDiamond);
1459 //Keep neutral meson selection histograms if requiered
1460 //Setting done in AliNeutralMesonSelection
1462 if(fAnaType!=kSSCalo && GetNeutralMesonSelection())
1464 TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;
1466 if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())
1467 for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) outputContainer->Add(nmsHistos->At(i)) ;
1472 return outputContainer ;
1476 //_____________________________________________
1477 Int_t AliAnaPi0EbE::GetMCIndex(const Int_t tag)
1480 // Assign mc index depending on MC bit set
1482 if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) )
1486 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) )
1490 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) &&
1491 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) )
1493 return kmcConversion ;
1494 }//conversion photon
1495 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) )
1498 }//photon no conversion
1499 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron))
1501 return kmcElectron ;
1510 //__________________________________________________________________
1511 void AliAnaPi0EbE::HasPairSameMCMother(AliAODPWG4Particle * photon1,
1512 AliAODPWG4Particle * photon2,
1513 Int_t & label, Int_t & tag)
1515 // Check the labels of pare in case mother was same pi0 or eta
1516 // Set the new AOD accordingly
1518 Int_t label1 = photon1->GetLabel();
1519 Int_t label2 = photon2->GetLabel();
1521 if(label1 < 0 || label2 < 0 ) return ;
1523 //Int_t tag1 = GetMCAnalysisUtils()->CheckOrigin(label1, GetReader(), photon1->GetInputFileIndex());
1524 //Int_t tag2 = GetMCAnalysisUtils()->CheckOrigin(label2, GetReader(), photon2->GetInputFileIndex());
1525 Int_t tag1 = photon1->GetTag();
1526 Int_t tag2 = photon2->GetTag();
1528 if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Origin of: photon1 %d; photon2 %d \n",tag1, tag2);
1529 if( (GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay) &&
1530 GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCPi0Decay) ) ||
1531 (GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCEtaDecay) &&
1532 GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCEtaDecay) )
1536 //Check if pi0/eta mother is the same
1537 if(GetReader()->ReadStack())
1541 TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree
1542 label1 = mother1->GetFirstMother();
1543 //mother1 = GetMCStack()->Particle(label1);//pi0
1547 TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree
1548 label2 = mother2->GetFirstMother();
1549 //mother2 = GetMCStack()->Particle(label2);//pi0
1552 else if(GetReader()->ReadAODMCParticles())
1553 {//&& (input > -1)){
1556 AliAODMCParticle * mother1 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(photon1->GetInputFileIndex()))->At(label1);//photon in kine tree
1557 label1 = mother1->GetMother();
1558 //mother1 = GetMCStack()->Particle(label1);//pi0
1562 AliAODMCParticle * mother2 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(photon2->GetInputFileIndex()))->At(label2);//photon in kine tree
1563 label2 = mother2->GetMother();
1564 //mother2 = GetMCStack()->Particle(label2);//pi0
1568 //printf("mother1 %d, mother2 %d\n",label1,label2);
1569 if( label1 == label2 && label1>=0 )
1574 TLorentzVector mom1 = *(photon1->Momentum());
1575 TLorentzVector mom2 = *(photon2->Momentum());
1577 Double_t angle = mom2.Angle(mom1.Vect());
1578 Double_t mass = (mom1+mom2).M();
1579 Double_t epair = (mom1+mom2).E();
1581 if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay))
1583 fhMassPairMCPi0 ->Fill(epair,mass);
1584 fhAnglePairMCPi0->Fill(epair,angle);
1585 GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCPi0);
1589 fhMassPairMCEta ->Fill(epair,mass);
1590 fhAnglePairMCEta->Fill(epair,angle);
1591 GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCEta);
1595 } // both from eta or pi0 decay
1599 //____________________________________________________________________________
1600 void AliAnaPi0EbE::Init()
1604 if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
1605 printf("AliAnaPi0EbE::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
1608 else if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
1609 printf("AliAnaPi0EbE::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
1615 //____________________________________________________________________________
1616 void AliAnaPi0EbE::InitParameters()
1618 //Initialize the parameters of the analysis.
1619 AddToHistogramsName("AnaPi0EbE_");
1621 fInputAODGammaConvName = "PhotonsCTS" ;
1622 fAnaType = kIMCalo ;
1623 fCalorimeter = "EMCAL" ;
1630 //__________________________________________________________________
1631 void AliAnaPi0EbE::MakeAnalysisFillAOD()
1633 //Do analysis and fill aods
1638 MakeInvMassInCalorimeter();
1642 MakeShowerShapeIdentification();
1646 MakeInvMassInCalorimeterAndCTS();
1652 //____________________________________________
1653 void AliAnaPi0EbE::MakeInvMassInCalorimeter()
1655 //Do analysis and fill aods
1656 //Search for the photon decay in calorimeters
1657 //Read photon list from AOD, produced in class AliAnaPhoton
1658 //Check if 2 photons have the mass of the pi0.
1660 TLorentzVector mom1;
1661 TLorentzVector mom2;
1662 TLorentzVector mom ;
1667 if(!GetInputAODBranch()){
1668 printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - No input calo photons in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
1672 //Get shower shape information of clusters
1673 TObjArray *clusters = 0;
1674 if (fCalorimeter=="EMCAL") clusters = GetEMCALClusters();
1675 else if(fCalorimeter=="PHOS") clusters = GetPHOSClusters() ;
1677 for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast()-1; iphoton++){
1678 AliAODPWG4Particle * photon1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
1680 //Vertex cut in case of mixed events
1681 Int_t evtIndex1 = 0 ;
1683 evtIndex1 = GetMixedEvent()->EventIndexForCaloCluster(photon1->GetCaloLabel(0)) ;
1684 if(TMath::Abs(GetVertex(evtIndex1)[2]) > GetZvertexCut()) continue ; //vertex cut
1685 mom1 = *(photon1->Momentum());
1687 //Get original cluster, to recover some information
1689 AliVCluster *cluster1 = FindCluster(clusters,photon1->GetCaloLabel(0),iclus);
1692 printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - First cluster not found\n");
1696 for(Int_t jphoton = iphoton+1; jphoton < GetInputAODBranch()->GetEntriesFast(); jphoton++)
1698 AliAODPWG4Particle * photon2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(jphoton));
1700 Int_t evtIndex2 = 0 ;
1702 evtIndex2 = GetMixedEvent()->EventIndexForCaloCluster(photon2->GetCaloLabel(0)) ;
1704 if(GetMixedEvent() && (evtIndex1 == evtIndex2))
1707 if(TMath::Abs(GetVertex(evtIndex2)[2]) > GetZvertexCut()) continue ; //vertex cut
1709 mom2 = *(photon2->Momentum());
1711 //Get original cluster, to recover some information
1713 AliVCluster *cluster2 = FindCluster(clusters,photon2->GetCaloLabel(0),iclus2,iclus+1);
1717 printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Second cluster not found\n");
1721 Float_t e1 = photon1->E();
1722 Float_t e2 = photon2->E();
1724 //Select clusters with good time window difference
1725 Float_t tof1 = cluster1->GetTOF()*1e9;;
1726 Float_t tof2 = cluster2->GetTOF()*1e9;;
1727 Double_t t12diff = tof1-tof2;
1728 fhEPairDiffTime->Fill(e1+e2, t12diff);
1729 if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
1731 //Play with the MC stack if available
1732 if(IsDataMC()) HasPairSameMCMother(photon1, photon2, label, tag) ;
1734 // Check the invariant mass for different selection on the local maxima
1735 // Name of AOD method TO BE FIXED
1736 Int_t nMaxima1 = photon1->GetFiducialArea();
1737 Int_t nMaxima2 = photon2->GetFiducialArea();
1739 Double_t mass = (mom1+mom2).M();
1740 Double_t epair = (mom1+mom2).E();
1742 if(nMaxima1==nMaxima2)
1744 if (nMaxima1==1) fhMassPairLocMax[0]->Fill(epair,mass);
1745 else if(nMaxima1==2) fhMassPairLocMax[1]->Fill(epair,mass);
1746 else fhMassPairLocMax[2]->Fill(epair,mass);
1748 else if(nMaxima1==1 || nMaxima2==1)
1750 if (nMaxima1==2 || nMaxima2==2) fhMassPairLocMax[3]->Fill(epair,mass);
1751 else fhMassPairLocMax[4]->Fill(epair,mass);
1754 fhMassPairLocMax[5]->Fill(epair,mass);
1756 // combinations with SS axis cut and NLM cut
1757 if(nMaxima1 == 1 && cluster2->GetM02() > 0.3) fhMassPairLocMax[6]->Fill(epair,mass);
1758 if(nMaxima2 == 1 && cluster1->GetM02() > 0.3) fhMassPairLocMax[6]->Fill(epair,mass);
1759 if(nMaxima1 > 1 && cluster2->GetM02() < 0.3 && cluster2->GetM02()> 0.1 ) fhMassPairLocMax[7]->Fill(epair,mass);
1760 if(nMaxima2 > 1 && cluster1->GetM02() < 0.3 && cluster1->GetM02()> 0.1 ) fhMassPairLocMax[7]->Fill(epair,mass);
1762 //Skip events with too few or too many NLM
1763 if((nMaxima1 < fNLMCutMin || nMaxima1 > fNLMCutMax) || (nMaxima2 < fNLMCutMin || nMaxima2 > fNLMCutMax)) continue ;
1765 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - NLM of out of range: cluster1 %d, cluster2 %d \n",nMaxima1, nMaxima2);
1768 fhMass->Fill(epair,(mom1+mom2).M());
1770 //Select good pair (good phi, pt cuts, aperture and invariant mass)
1771 if(GetNeutralMesonSelection()->SelectPair(mom1, mom2,fCalorimeter))
1774 printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Selected gamma pair: pt %f, phi %f, eta%f \n",(mom1+mom2).Pt(), (mom1+mom2).Phi()*180./3.1416, (mom1+mom2).Eta());
1776 //Fill some histograms about shower shape
1777 if(fFillSelectClHisto && clusters && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
1779 FillSelectedClusterHistograms(cluster1, nMaxima1, photon1->GetTag());
1780 FillSelectedClusterHistograms(cluster2, nMaxima2, photon2->GetTag());
1783 // Tag both photons as decay
1784 photon1->SetTagged(kTRUE);
1785 photon2->SetTagged(kTRUE);
1787 fhPtDecay->Fill(photon1->Pt());
1788 fhEDecay ->Fill(photon1->E() );
1790 fhPtDecay->Fill(photon2->Pt());
1791 fhEDecay ->Fill(photon2->E() );
1793 //Create AOD for analysis
1796 //Mass of selected pairs
1797 fhSelectedMass->Fill(epair,mom.M());
1799 // Fill histograms to undertand pile-up before other cuts applied
1800 // Remember to relax time cuts in the reader
1801 FillPileUpHistograms(mom.E(),((cluster1->GetTOF()+cluster2->GetTOF())*1e9) /2);
1803 AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
1805 pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
1806 pi0.SetDetector(photon1->GetDetector());
1809 pi0.SetLabel(label);
1812 //Set the indeces of the original caloclusters
1813 pi0.SetCaloLabel(photon1->GetCaloLabel(0), photon2->GetCaloLabel(0));
1814 //pi0.SetInputFileIndex(input);
1816 AddAODParticle(pi0);
1824 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - End fill AODs \n");
1828 //__________________________________________________
1829 void AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS()
1831 //Do analysis and fill aods
1832 //Search for the photon decay in calorimeters
1833 //Read photon list from AOD, produced in class AliAnaPhoton and AliGammaConversion
1834 //Check if 2 photons have the mass of the pi0.
1836 TLorentzVector mom1;
1837 TLorentzVector mom2;
1838 TLorentzVector mom ;
1843 // Check calorimeter input
1844 if(!GetInputAODBranch()){
1845 printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input calo photons in AOD branch with name < %s > , STOP\n",GetInputAODName().Data());
1849 // Get the array with conversion photons
1850 TClonesArray * inputAODGammaConv = (TClonesArray *) GetReader()->GetOutputEvent()->FindListObject(fInputAODGammaConvName);
1851 if(!inputAODGammaConv) {
1853 inputAODGammaConv = (TClonesArray *) GetReader()->GetInputEvent()->FindListObject(fInputAODGammaConvName);
1855 if(!inputAODGammaConv) {
1856 printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input gamma conversions in AOD branch with name < %s >\n",fInputAODGammaConvName.Data());
1862 //Get shower shape information of clusters
1863 TObjArray *clusters = 0;
1864 if (fCalorimeter=="EMCAL") clusters = GetEMCALClusters();
1865 else if(fCalorimeter=="PHOS") clusters = GetPHOSClusters() ;
1867 Int_t nCTS = inputAODGammaConv->GetEntriesFast();
1868 Int_t nCalo = GetInputAODBranch()->GetEntriesFast();
1869 if(nCTS<=0 || nCalo <=0)
1871 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - nCalo %d, nCTS %d, cannot loop\n",nCalo,nCTS);
1876 printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - Number of conversion photons %d\n",nCTS);
1878 // Do the loop, first calo, second CTS
1879 for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast(); iphoton++){
1880 AliAODPWG4Particle * photon1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
1881 mom1 = *(photon1->Momentum());
1883 //Get original cluster, to recover some information
1885 AliVCluster *cluster = FindCluster(clusters,photon1->GetCaloLabel(0),iclus);
1887 for(Int_t jphoton = 0; jphoton < nCTS; jphoton++){
1888 AliAODPWG4Particle * photon2 = (AliAODPWG4Particle*) (inputAODGammaConv->At(jphoton));
1890 evtIndex = GetMixedEvent()->EventIndexForCaloCluster(photon2->GetCaloLabel(0)) ;
1891 if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ; //vertex cut
1893 mom2 = *(photon2->Momentum());
1895 Double_t mass = (mom1+mom2).M();
1896 Double_t epair = (mom1+mom2).E();
1898 Int_t nMaxima = photon1->GetFiducialArea();
1899 if (nMaxima==1) fhMassPairLocMax[0]->Fill(epair,mass);
1900 else if(nMaxima==2) fhMassPairLocMax[1]->Fill(epair,mass);
1901 else fhMassPairLocMax[2]->Fill(epair,mass);
1903 if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax) continue ;
1904 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - NLM %d of out of range \n",nMaxima);
1906 //Play with the MC stack if available
1909 Int_t label2 = photon2->GetLabel();
1910 if(label2 >= 0 )photon2->SetTag(GetMCAnalysisUtils()->CheckOrigin(label2, GetReader(), photon2->GetInputFileIndex()));
1912 HasPairSameMCMother(photon1, photon2, label, tag) ;
1915 //Mass of selected pairs
1916 fhMass->Fill(epair,(mom1+mom2).M());
1918 //Select good pair (good phi, pt cuts, aperture and invariant mass)
1919 if(GetNeutralMesonSelection()->SelectPair(mom1, mom2,fCalorimeter))
1921 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - Selected gamma pair: pt %f, phi %f, eta%f\n",(mom1+mom2).Pt(), (mom1+mom2).Phi()*180./3.1416, (mom1+mom2).Eta());
1923 //Fill some histograms about shower shape
1924 if(fFillSelectClHisto && cluster && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
1926 FillSelectedClusterHistograms(cluster, nMaxima, photon1->GetTag());
1929 // Tag both photons as decay
1930 photon1->SetTagged(kTRUE);
1931 photon2->SetTagged(kTRUE);
1933 fhPtDecay->Fill(photon1->Pt());
1934 fhEDecay ->Fill(photon1->E() );
1936 //Create AOD for analysis
1940 //Mass of selected pairs
1941 fhSelectedMass->Fill(epair,mom.M());
1943 // Fill histograms to undertand pile-up before other cuts applied
1944 // Remember to relax time cuts in the reader
1945 if(cluster)FillPileUpHistograms(mom.E(),cluster->GetTOF()*1e9);
1947 AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
1949 pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
1950 pi0.SetDetector(photon1->GetDetector());
1953 pi0.SetLabel(label);
1956 //Set the indeces of the original tracks or caloclusters
1957 pi0.SetCaloLabel(photon1->GetCaloLabel(0), -1);
1958 pi0.SetTrackLabel(photon2->GetTrackLabel(0), photon2->GetTrackLabel(1));
1959 //pi0.SetInputFileIndex(input);
1961 AddAODParticle(pi0);
1968 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - End fill AODs \n");
1973 //_________________________________________________
1974 void AliAnaPi0EbE::MakeShowerShapeIdentification()
1976 //Search for pi0 in fCalorimeter with shower shape analysis
1978 TObjArray * pl = 0x0;
1979 AliVCaloCells * cells = 0x0;
1980 //Select the Calorimeter of the photon
1981 if (fCalorimeter == "PHOS" )
1983 pl = GetPHOSClusters();
1984 cells = GetPHOSCells();
1986 else if (fCalorimeter == "EMCAL")
1988 pl = GetEMCALClusters();
1989 cells = GetEMCALCells();
1994 Info("MakeShowerShapeIdentification","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
1998 TLorentzVector mom ;
1999 for(Int_t icalo = 0; icalo < pl->GetEntriesFast(); icalo++)
2001 AliVCluster * calo = (AliVCluster*) (pl->At(icalo));
2003 Int_t evtIndex = 0 ;
2004 if (GetMixedEvent())
2006 evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
2009 if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ; //vertex cut
2011 //Get Momentum vector,
2012 Double_t vertex[]={0,0,0};
2013 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
2015 calo->GetMomentum(mom,GetVertex(evtIndex)) ;
2016 }//Assume that come from vertex in straight line
2019 calo->GetMomentum(mom,vertex) ;
2022 //If too small or big pt, skip it
2023 if(mom.E() < GetMinEnergy() || mom.E() > GetMaxEnergy() ) continue ;
2025 //Check acceptance selection
2026 if(IsFiducialCutOn())
2028 Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
2029 if(! in ) continue ;
2033 printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: Min pt cut and fiducial cut passed: pt %3.2f, phi %2.2f, eta %1.2f\n",mom.Pt(),mom.Phi(),mom.Eta());
2035 //Check Distance to Bad channel, set bit.
2036 Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
2037 if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
2038 if(distBad < fMinDist) //In bad channel (PHOS cristal size 2.2x2.2 cm)
2041 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: Bad channel cut passed %4.2f\n",distBad);
2043 //.......................................
2044 // TOF cut, BE CAREFUL WITH THIS CUT
2045 Double_t tof = calo->GetTOF()*1e9;
2046 if(tof < fTimeCutMin || tof > fTimeCutMax) continue ;
2048 //Play with the MC stack if available
2049 //Check origin of the candidates
2053 tag = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader(),0);
2054 //GetMCAnalysisUtils()->CheckMultipleOrigin(calo->GetLabels(),calo->GetNLabels(), GetReader(), aodpi0.GetInputFileIndex(), tag);
2055 if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Origin of candidate %d\n",tag);
2058 //Skip matched clusters with tracks
2059 if(IsTrackMatched(calo, GetReader()->GetInputEvent()))
2061 FillRejectedClusterHistograms(mom,tag);
2067 //PID selection or bit setting
2069 Double_t mass = 0 , angle = 0;
2070 Double_t e1 = 0 , e2 = 0;
2071 Int_t idPartType = GetCaloPID()->GetIdentifiedParticleTypeFromClusterSplitting(calo,cells,GetCaloUtils(),
2072 GetVertex(evtIndex),nMaxima,
2075 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - PDG of identified particle %d\n",idPartType);
2078 //Skip events with too few or too many NLM
2079 if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax)
2081 FillRejectedClusterHistograms(mom,tag);
2086 printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - NLM %d accepted \n",nMaxima);
2088 //mass of all clusters
2089 fhMass->Fill(mom.E(),mass);
2091 // Asymmetry of all clusters
2093 if(e1+e2 > 0) asy = (e1-e2) / (e1+e2);
2094 fhAsymmetry->Fill(mom.E(),asy);
2098 Int_t mcIndex = GetMCIndex(tag);
2099 fhMCEAsymmetry[mcIndex]->Fill(mom.E(),asy);
2102 // If cluster does not pass pid, not pi0/eta, skip it.
2103 if (GetOutputAODName().Contains("Pi0") && idPartType != AliCaloPID::kPi0)
2105 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Cluster is not Pi0\n");
2106 FillRejectedClusterHistograms(mom,tag);
2110 else if(GetOutputAODName().Contains("Eta") && idPartType != AliCaloPID::kEta)
2112 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Cluster is not Eta\n");
2113 FillRejectedClusterHistograms(mom,tag);
2118 printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Pi0/Eta selection cuts passed: pT %3.2f, pdg %d\n",
2119 mom.Pt(), idPartType);
2121 fhSelectedAsymmetry->Fill(mom.E(),asy);
2123 if( fUseSplitAsyCut && GetCaloPID()->IsInPi0SplitAsymmetryRange(mom.E(),asy,nMaxima) )
2125 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Too large asymmetry\n");
2126 FillRejectedClusterHistograms(mom,tag);
2130 //Mass of selected pairs
2131 fhSelectedMass ->Fill(mom.E(),mass);
2133 //-----------------------
2134 //Create AOD for analysis
2136 AliAODPWG4Particle aodpi0 = AliAODPWG4Particle(mom);
2137 aodpi0.SetLabel(calo->GetLabel());
2139 //Set the indeces of the original caloclusters
2140 aodpi0.SetCaloLabel(calo->GetID(),-1);
2141 aodpi0.SetDetector(fCalorimeter);
2143 if (distBad > fMinDist3) aodpi0.SetDistToBad(2) ;
2144 else if(distBad > fMinDist2) aodpi0.SetDistToBad(1) ;
2145 else aodpi0.SetDistToBad(0) ;
2147 // Check if cluster is pi0 via cluster splitting
2148 aodpi0.SetIdentifiedParticleType(idPartType);
2150 // Add number of local maxima to AOD, method name in AOD to be FIXED
2151 aodpi0.SetFiducialArea(nMaxima);
2155 //Fill some histograms about shower shape
2156 if(fFillSelectClHisto && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
2158 FillSelectedClusterHistograms(calo, nMaxima, tag, asy);
2161 // Fill histograms to undertand pile-up before other cuts applied
2162 // Remember to relax time cuts in the reader
2163 FillPileUpHistograms(calo->E(),calo->GetTOF()*1e9);
2165 //Add AOD with pi0 object to aod branch
2166 AddAODParticle(aodpi0);
2170 if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - End fill AODs \n");
2173 //______________________________________________
2174 void AliAnaPi0EbE::MakeAnalysisFillHistograms()
2176 //Do analysis and fill histograms
2178 if(!GetOutputAODBranch())
2180 printf("AliAnaPi0EbE::MakeAnalysisFillHistograms() - No output pi0 in AOD branch with name < %s >,STOP \n",GetOutputAODName().Data());
2183 //Loop on stored AOD pi0
2184 Int_t naod = GetOutputAODBranch()->GetEntriesFast();
2185 if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
2187 for(Int_t iaod = 0; iaod < naod ; iaod++)
2190 AliAODPWG4Particle* pi0 = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
2191 Int_t pdg = pi0->GetIdentifiedParticleType();
2193 if(IsCaloPIDOn() && pdg != AliCaloPID::kPi0) continue;
2195 //Fill pi0 histograms
2196 Float_t ener = pi0->E();
2197 Float_t pt = pi0->Pt();
2198 Float_t phi = pi0->Phi();
2199 if(phi < 0) phi+=TMath::TwoPi();
2200 Float_t eta = pi0->Eta();
2205 fhEEta ->Fill(ener,eta);
2206 fhEPhi ->Fill(ener,phi);
2207 fhEtaPhi ->Fill(eta,phi);
2209 if(fFillPileUpHistograms)
2211 if(GetReader()->IsPileUpFromSPD()) fhPtPi0PileUp[0]->Fill(pt);
2212 if(GetReader()->IsPileUpFromEMCal()) fhPtPi0PileUp[1]->Fill(pt);
2213 if(GetReader()->IsPileUpFromSPDOrEMCal()) fhPtPi0PileUp[2]->Fill(pt);
2214 if(GetReader()->IsPileUpFromSPDAndEMCal()) fhPtPi0PileUp[3]->Fill(pt);
2215 if(GetReader()->IsPileUpFromSPDAndNotEMCal()) fhPtPi0PileUp[4]->Fill(pt);
2216 if(GetReader()->IsPileUpFromEMCalAndNotSPD()) fhPtPi0PileUp[5]->Fill(pt);
2217 if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) fhPtPi0PileUp[6]->Fill(pt);
2223 Int_t tag = pi0->GetTag();
2224 Int_t mcIndex = GetMCIndex(tag);
2226 fhMCE [mcIndex] ->Fill(ener);
2227 fhMCPt [mcIndex] ->Fill(pt);
2228 fhMCPhi[mcIndex] ->Fill(pt,phi);
2229 fhMCEta[mcIndex] ->Fill(pt,eta);
2231 if((mcIndex==kmcPhoton || mcIndex==kmcPi0 || mcIndex==kmcEta) && fAnaType==kSSCalo)
2233 Float_t efracMC = 0;
2234 Int_t label = pi0->GetLabel();
2237 TLorentzVector mom = GetMCAnalysisUtils()->GetMother(label,GetReader(),ok);
2240 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0))
2242 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,111,GetReader(),ok);
2243 if(grandmom.E() > 0 && ok)
2245 efracMC = grandmom.E()/ener;
2246 fhMCPi0PtGenRecoFraction ->Fill(pt,efracMC);
2249 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))
2251 fhMCPi0DecayPt->Fill(pt);
2252 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,111,GetReader(),ok);
2253 if(grandmom.E() > 0 && ok)
2255 efracMC = mom.E()/grandmom.E();
2256 fhMCPi0DecayPtFraction ->Fill(pt,efracMC);
2259 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta))
2261 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,221,GetReader(),ok);
2262 if(grandmom.E() > 0 && ok)
2264 efracMC = grandmom.E()/ener;
2265 fhMCEtaPtGenRecoFraction ->Fill(pt,efracMC);
2268 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay))
2270 fhMCEtaDecayPt->Fill(pt);
2271 TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,221,GetReader(),ok);
2272 if(grandmom.E() > 0 && ok)
2274 efracMC = mom.E()/grandmom.E();
2275 fhMCEtaDecayPtFraction ->Fill(pt,efracMC);
2278 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))
2280 fhMCOtherDecayPt->Fill(pt);
2285 }//Histograms with MC
2291 //__________________________________________________________________
2292 void AliAnaPi0EbE::Print(const Option_t * opt) const
2294 //Print some relevant parameters set for the analysis
2298 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
2299 AliAnaCaloTrackCorrBaseClass::Print("");
2300 printf("Analysis Type = %d \n", fAnaType) ;
2301 if(fAnaType == kSSCalo){
2302 printf("Calorimeter = %s\n", fCalorimeter.Data()) ;
2303 printf("Min Distance to Bad Channel = %2.1f\n",fMinDist);
2304 printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
2305 printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);