]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGGA/CaloTrackCorrelations/AliAnaCalorimeterQA.cxx
1. JoinTrees Index - Possible to specify the branch without prefix ()
[u/mrichter/AliRoot.git] / PWGGA / CaloTrackCorrelations / AliAnaCalorimeterQA.cxx
CommitLineData
9725fd2a 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 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 **************************************************************************/
9725fd2a 15
a6f26052 16//_________________________________________________________________________
17// Class to check results from simulations or reconstructed real data.
18// Fill few histograms and do some checking plots
19//
20//-- Author: Gustavo Conesa (INFN-LNF)
21//_________________________________________________________________________
9725fd2a 22
23
a6f26052 24// --- ROOT system ---
d55bb5e1 25#include <TObjArray.h>
26#include <TParticle.h>
27#include <TDatabasePDG.h>
28#include <TH3F.h>
0c1383b5 29#include <TObjString.h>
9725fd2a 30
a6f26052 31//---- AliRoot system ----
9725fd2a 32#include "AliAnaCalorimeterQA.h"
33#include "AliCaloTrackReader.h"
34#include "AliStack.h"
c8fe2783 35#include "AliVCaloCells.h"
ff45398a 36#include "AliFiducialCut.h"
c8fe2783 37#include "AliVCluster.h"
d55bb5e1 38#include "AliVTrack.h"
c8fe2783 39#include "AliVEvent.h"
902aa95c 40#include "AliVEventHandler.h"
902aa95c 41#include "AliAODMCParticle.h"
42#include "AliMCAnalysisUtils.h"
9725fd2a 43
c5693f62 44// --- Detectors ---
45#include "AliPHOSGeoUtils.h"
46#include "AliEMCALGeometry.h"
47
9725fd2a 48ClassImp(AliAnaCalorimeterQA)
c8fe2783 49
649b825d 50//________________________________________
c8fe2783 51AliAnaCalorimeterQA::AliAnaCalorimeterQA() :
765206a5 52AliAnaCaloTrackCorrBaseClass(), fCalorimeter(""),
649b825d 53
54//Switches
e6fec6f5 55fFillAllCellTimeHisto(kTRUE),
9e9f04cb 56fFillAllPosHisto(kFALSE), fFillAllPosHisto2(kTRUE),
57fFillAllTH12(kFALSE), fFillAllTH3(kTRUE),
35c71d5c 58fFillAllTMHisto(kTRUE), fFillAllPi0Histo(kTRUE),
649b825d 59fCorrelate(kTRUE), fStudyBadClusters(kFALSE),
f1538a5f 60fStudyClustersAsymmetry(kFALSE), fStudyExotic(kFALSE),
61fStudyWeight(kFALSE),
649b825d 62
63//Parameters and cuts
35c71d5c 64fNModules(12), fNRCU(2),
65fNMaxCols(48), fNMaxRows(24),
f15c25da 66fTimeCutMin(-10000), fTimeCutMax(10000),
9e9f04cb 67fEMCALCellAmpMin(0), fPHOSCellAmpMin(0),
649b825d 68
f1538a5f 69// Exotic
70fExoNECrossCuts(0), fExoECrossCuts(),
71fExoNDTimeCuts(0), fExoDTimeCuts(),
72
649b825d 73//Histograms
9e9f04cb 74fhE(0), fhPt(0),
75fhPhi(0), fhEta(0), fhEtaPhiE(0),
76fhECharged(0), fhPtCharged(0),
77fhPhiCharged(0), fhEtaCharged(0), fhEtaPhiECharged(0),
521636d2 78
79//Invariant mass
9e9f04cb 80fhIM(0 ), fhAsym(0),
a82b4462 81
82fhNCellsPerCluster(0), fhNCellsPerClusterNoCut(0), fhNClusters(0),
3129a79e 83
e1e62b89 84//Timing
9e9f04cb 85fhClusterTimeEnergy(0), fhCellTimeSpreadRespectToCellMax(0),
9e9f04cb 86fhCellIdCellLargeTimeSpread(0), fhClusterPairDiffTimeE(0),
9e9f04cb 87fhClusterMaxCellCloseCellRatio(0), fhClusterMaxCellCloseCellDiff(0),
88fhClusterMaxCellDiff(0), fhClusterMaxCellDiffNoCut(0),
a82b4462 89fhClusterMaxCellDiffAverageTime(0), fhClusterMaxCellDiffWeightedTime(0),
1a72f6c5 90fhClusterMaxCellECross(0),
649b825d 91fhLambda0(0), fhLambda1(0), fhDispersion(0),
715fd81f 92
649b825d 93//bad clusters
94fhBadClusterEnergy(0), fhBadClusterTimeEnergy(0),
95fhBadClusterPairDiffTimeE(0), fhBadCellTimeSpreadRespectToCellMax(0),
9e9f04cb 96fhBadClusterMaxCellCloseCellRatio(0), fhBadClusterMaxCellCloseCellDiff(0), fhBadClusterMaxCellDiff(0),
a82b4462 97fhBadClusterMaxCellDiffAverageTime(0), fhBadClusterMaxCellDiffWeightedTime(0),
1a72f6c5 98fhBadClusterMaxCellECross(0),
1a83b960 99fhBadClusterDeltaIEtaDeltaIPhiE0(0), fhBadClusterDeltaIEtaDeltaIPhiE2(0),
100fhBadClusterDeltaIEtaDeltaIPhiE6(0), fhBadClusterDeltaIA(0),
9e9f04cb 101
521636d2 102//Position
9e9f04cb 103fhRNCells(0), fhXNCells(0),
104fhYNCells(0), fhZNCells(0),
105fhRE(0), fhXE(0),
106fhYE(0), fhZE(0),
521636d2 107fhXYZ(0),
9e9f04cb 108fhRCellE(0), fhXCellE(0),
109fhYCellE(0), fhZCellE(0),
521636d2 110fhXYZCell(0),
9e9f04cb 111fhDeltaCellClusterRNCells(0), fhDeltaCellClusterXNCells(0),
112fhDeltaCellClusterYNCells(0), fhDeltaCellClusterZNCells(0),
113fhDeltaCellClusterRE(0), fhDeltaCellClusterXE(0),
114fhDeltaCellClusterYE(0), fhDeltaCellClusterZE(0),
649b825d 115
521636d2 116// Cells
9e9f04cb 117fhNCells(0), fhAmplitude(0),
118fhAmpId(0), fhEtaPhiAmp(0),
1a72f6c5 119fhTime(0), fhTimeVz(0),
120fhTimeId(0), fhTimeAmp(0),
121fhCellECross(0),
9e9f04cb 122fhCaloCorrNClusters(0), fhCaloCorrEClusters(0),
123fhCaloCorrNCells(0), fhCaloCorrECells(0),
124fhCaloV0SCorrNClusters(0), fhCaloV0SCorrEClusters(0),
125fhCaloV0SCorrNCells(0), fhCaloV0SCorrECells(0),
126fhCaloV0MCorrNClusters(0), fhCaloV0MCorrEClusters(0),
127fhCaloV0MCorrNCells(0), fhCaloV0MCorrECells(0),
128fhCaloTrackMCorrNClusters(0), fhCaloTrackMCorrEClusters(0),
129fhCaloTrackMCorrNCells(0), fhCaloTrackMCorrECells(0),
649b825d 130
521636d2 131//Super-Module dependent histgrams
649b825d 132fhEMod(0), fhAmpMod(0), fhTimeMod(0),
133fhNClustersMod(0), fhNCellsMod(0),
134fhNCellsPerClusterMod(0), fhNCellsPerClusterModNoCut(0),
135
136fhGridCells(0), fhGridCellsE(0), fhGridCellsTime(0),
137fhTimeAmpPerRCU(0), fhIMMod(0),
138
139// Weight studies
140fhECellClusterRatio(0), fhECellClusterLogRatio(0),
141fhEMaxCellClusterRatio(0), fhEMaxCellClusterLogRatio(0),
715fd81f 142
765206a5 143fhExoL0ECross(0), fhExoL1ECross(0),
144
715fd81f 145// MC and reco
649b825d 146fhRecoMCE(), fhRecoMCPhi(), fhRecoMCEta(),
147fhRecoMCDeltaE(), fhRecoMCRatioE(),
148fhRecoMCDeltaPhi(), fhRecoMCDeltaEta(),
715fd81f 149
521636d2 150// MC only
649b825d 151fhGenMCE(), fhGenMCEtaPhi(),
152fhGenMCAccE(), fhGenMCAccEtaPhi(),
35c71d5c 153
154//matched MC
649b825d 155fhEMVxyz(0), fhEMR(0),
156fhHaVxyz(0), fhHaR(0),
d55bb5e1 157fh1EOverP(0), fh2dR(0),
649b825d 158fh2EledEdx(0), fh2MatchdEdx(0),
d55bb5e1 159fhMCEle1EOverP(0), fhMCEle1dR(0), fhMCEle2MatchdEdx(0),
160fhMCChHad1EOverP(0), fhMCChHad1dR(0), fhMCChHad2MatchdEdx(0),
161fhMCNeutral1EOverP(0), fhMCNeutral1dR(0), fhMCNeutral2MatchdEdx(0), fh1EOverPR02(0),
162fhMCEle1EOverPR02(0), fhMCChHad1EOverPR02(0), fhMCNeutral1EOverPR02(0)
9725fd2a 163{
a6f26052 164 //Default Ctor
afabc52f 165
649b825d 166 //Weight studies
1a72f6c5 167 for(Int_t i =0; i < 14; i++){
649b825d 168 fhLambda0ForW0[i] = 0;
1a72f6c5 169 //fhLambda1ForW0[i] = 0;
649b825d 170
171 for(Int_t j = 0; j < 5; j++){
172 fhLambda0ForW0MC[i][j] = 0;
1a72f6c5 173 //fhLambda1ForW0MC[i][j] = 0;
649b825d 174 }
175
176 }
c8fe2783 177
649b825d 178 //Cluster size
179 fhDeltaIEtaDeltaIPhiE0[0] = 0 ; fhDeltaIEtaDeltaIPhiE2[0] = 0; fhDeltaIEtaDeltaIPhiE6[0] = 0;
180 fhDeltaIEtaDeltaIPhiE0[1] = 0 ; fhDeltaIEtaDeltaIPhiE2[1] = 0; fhDeltaIEtaDeltaIPhiE6[1] = 0;
181 fhDeltaIA[0] = 0 ; fhDeltaIAL0[0] = 0; fhDeltaIAL1[0] = 0;
182 fhDeltaIA[1] = 0 ; fhDeltaIAL0[1] = 0; fhDeltaIAL1[1] = 0;
183 fhDeltaIANCells[0] = 0 ; fhDeltaIANCells[1] = 0;
184 fhDeltaIAMC[0] = 0 ; fhDeltaIAMC[1] = 0;
185 fhDeltaIAMC[2] = 0 ; fhDeltaIAMC[3] = 0;
2302a644 186
f1538a5f 187 // Exotic
188 for (Int_t ie = 0; ie < 10 ; ie++)
189 {
190 fhExoDTime[ie] = 0;
191 for (Int_t idt = 0; idt < 5 ; idt++)
192 {
193 fhExoNCell [ie][idt] = 0;
194 fhExoL0 [ie][idt] = 0;
765206a5 195 fhExoL1 [ie][idt] = 0;
f1538a5f 196 fhExoECross [ie][idt] = 0;
197 fhExoTime [ie][idt] = 0;
765206a5 198 fhExoL0NCell [ie][idt] = 0;
199 fhExoL1NCell [ie][idt] = 0;
f1538a5f 200 }
201 }
202
649b825d 203 // MC
c8fe2783 204
649b825d 205 for(Int_t i = 0; i < 6; i++){
206
207 fhRecoMCE[i][0] = 0; fhRecoMCE[i][1] = 0;
208 fhRecoMCPhi[i][0] = 0; fhRecoMCPhi[i][1] = 0;
209 fhRecoMCEta[i][0] = 0; fhRecoMCEta[i][1] = 0;
210 fhRecoMCDeltaE[i][0] = 0; fhRecoMCDeltaE[i][1] = 0;
211 fhRecoMCRatioE[i][0] = 0; fhRecoMCRatioE[i][1] = 0;
212 fhRecoMCDeltaPhi[i][0] = 0; fhRecoMCDeltaPhi[i][1] = 0;
213 fhRecoMCDeltaEta[i][0] = 0; fhRecoMCDeltaEta[i][1] = 0;
214
215 }
2302a644 216
649b825d 217 //Initialize parameters
218 InitParameters();
219}
3f5990d6 220
1a83b960 221//____________________________________________________________________________________________________________________________
c5693f62 222void AliAnaCalorimeterQA::BadClusterHistograms(AliVCluster* clus, const TObjArray *caloClusters, AliVCaloCells * cells,
1a83b960 223 const Int_t absIdMax, const Double_t maxCellFraction, const Float_t eCrossFrac,
a82b4462 224 const Double_t tmax, Double_t timeAverages[2]
649b825d 225 )
226{
227 //Bad cluster histograms
b0114dba 228
229 // printf("AliAnaCalorimeterQA::BadClusterHistograms() - Event %d - Calorimeter %s \n \t E %f, n cells %d, max cell absId %d, maxCellFrac %f\n",
230 // GetReader()->GetEventNumber(), fCalorimeter.Data(),
231 // clus->E(),clus->GetNCells(),absIdMax,maxCellFraction);
1a72f6c5 232
649b825d 233 fhBadClusterEnergy ->Fill(clus->E());
234 Double_t tof = clus->GetTOF()*1.e9;
1a83b960 235 fhBadClusterTimeEnergy ->Fill(clus->E(),tof);
236 fhBadClusterMaxCellDiff ->Fill(clus->E(),maxCellFraction);
237 fhBadClusterMaxCellECross->Fill(clus->E(),eCrossFrac);
238
239 if(fStudyClustersAsymmetry) ClusterAsymmetryHistograms(clus,absIdMax,kFALSE);
1a72f6c5 240
a82b4462 241 //Clusters in event time differencem bad minus good
2302a644 242
649b825d 243 for(Int_t iclus2 = 0; iclus2 < caloClusters->GetEntriesFast(); iclus2++ ){
a6f26052 244
649b825d 245 AliVCluster* clus2 = (AliVCluster*)caloClusters->At(iclus2);
a6f26052 246
649b825d 247 if(clus->GetID()==clus2->GetID()) continue;
715fd81f 248
a82b4462 249 Float_t maxCellFraction2 = 0.;
250 Int_t absIdMax2 = GetCaloUtils()->GetMaxEnergyCell(cells, clus2,maxCellFraction2);
251 if(IsGoodCluster(absIdMax2,cells)){
252 Double_t tof2 = clus2->GetTOF()*1.e9;
649b825d 253 fhBadClusterPairDiffTimeE ->Fill(clus->E(), (tof-tof2));
649b825d 254 }
a82b4462 255
649b825d 256 } // loop
924e319f 257
649b825d 258 // Max cell compared to other cells in cluster
e6fec6f5 259 if(fFillAllCellTimeHisto)
2747966a 260 {
649b825d 261 fhBadClusterMaxCellDiffAverageTime ->Fill(clus->E(),tmax-timeAverages[0]);
649b825d 262 fhBadClusterMaxCellDiffWeightedTime ->Fill(clus->E(),tmax-timeAverages[1]);
649b825d 263 }
715fd81f 264
2747966a 265 for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
266 {
649b825d 267 Int_t absId = clus->GetCellsAbsId()[ipos];
2747966a 268 if(absId!=absIdMax && cells->GetCellAmplitude(absIdMax) > 0.01)
269 {
649b825d 270 Float_t frac = cells->GetCellAmplitude(absId)/cells->GetCellAmplitude(absIdMax);
271
272 fhBadClusterMaxCellCloseCellRatio->Fill(clus->E(),frac);
273 fhBadClusterMaxCellCloseCellDiff ->Fill(clus->E(),cells->GetCellAmplitude(absIdMax)-cells->GetCellAmplitude(absId));
274
e6fec6f5 275 if(fFillAllCellTimeHisto)
dbba06ca 276 {
649b825d 277 Double_t time = cells->GetCellTime(absId);
dbba06ca 278 GetCaloUtils()->RecalibrateCellTime(time, fCalorimeter, absId,GetReader()->GetInputEvent()->GetBunchCrossNumber());
279
649b825d 280 Float_t diff = (tmax-time*1e9);
281 fhBadCellTimeSpreadRespectToCellMax->Fill(clus->E(), diff);
282
e6fec6f5 283 }
649b825d 284 }// Not max
285 }//loop
715fd81f 286
649b825d 287}
715fd81f 288
dbf54f1e 289//______________________________________________________________________
290void AliAnaCalorimeterQA::CalculateAverageTime(AliVCluster *clus,
291 AliVCaloCells* cells,
292 Double_t timeAverages[2])
649b825d 293{
294 // Calculate time averages and weights
a95eac90 295
649b825d 296 // First recalculate energy in case non linearity was applied
297 Float_t energy = 0;
298 Float_t ampMax = 0, amp = 0;
299 Int_t absIdMax =-1;
2747966a 300 for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
301 {
649b825d 302 Int_t id = clus->GetCellsAbsId()[ipos];
e1e62b89 303
649b825d 304 //Recalibrate cell energy if needed
305 amp = cells->GetCellAmplitude(id);
dbba06ca 306 GetCaloUtils()->RecalibrateCellAmplitude(amp,fCalorimeter, id);
4c8f7c2e 307
649b825d 308 energy += amp;
e1e62b89 309
2747966a 310 if(amp> ampMax)
311 {
649b825d 312 ampMax = amp;
313 absIdMax = id;
314 }
4c8f7c2e 315
649b825d 316 } // energy loop
317
318 // Calculate average time of cells in cluster and weighted average
dbf54f1e 319 Double_t aTime = 0;
320 Double_t wTime = 0;
321 Float_t wTot = 0;
322 Double_t time = 0;
323 Int_t id =-1;
324 Double_t w = 0;
325 Int_t ncells = clus->GetNCells();
2747966a 326 for (Int_t ipos = 0; ipos < ncells; ipos++)
327 {
dbf54f1e 328 id = clus ->GetCellsAbsId()[ipos];
329 amp = cells->GetCellAmplitude(id);
330 time = cells->GetCellTime(id);
649b825d 331
332 //Recalibrate energy and time
dbba06ca 333 GetCaloUtils()->RecalibrateCellAmplitude(amp , fCalorimeter, id);
334 GetCaloUtils()->RecalibrateCellTime (time, fCalorimeter, id, GetReader()->GetInputEvent()->GetBunchCrossNumber());
335
dbf54f1e 336 w = GetCaloUtils()->GetEMCALRecoUtils()->GetCellWeight(cells->GetCellAmplitude(id),energy);
649b825d 337 aTime += time*1e9;
338 wTime += time*1e9 * w;
339 wTot += w;
35c71d5c 340
649b825d 341 }
f3138ecf 342
dbf54f1e 343 if(ncells > 0) aTime /= ncells;
344 else aTime = 0;
649b825d 345
dbf54f1e 346 if(wTot > 0) wTime /= wTot;
f3138ecf 347 else wTime = 0;
348
dbf54f1e 349 timeAverages[0] = aTime;
350 timeAverages[1] = wTime;
39de6caa 351
649b825d 352}
353
354//____________________________________________________________
355void AliAnaCalorimeterQA::CellHistograms(AliVCaloCells *cells)
356{
357 // Plot histograms related to cells only
9e9f04cb 358
649b825d 359 Int_t ncells = cells->GetNumberOfCells();
39de6caa 360
649b825d 361 if(GetDebug() > 0)
362 printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - %s cell entries %d\n", fCalorimeter.Data(), ncells );
363
364 //Init arrays and used variables
365 Int_t *nCellsInModule = new Int_t[fNModules];
366 for(Int_t imod = 0; imod < fNModules; imod++ ) nCellsInModule[imod] = 0;
a82b4462 367
649b825d 368 Int_t icol = -1;
369 Int_t irow = -1;
370 Int_t iRCU = -1;
371 Float_t amp = 0.;
372 Double_t time = 0.;
373 Int_t id = -1;
374 Float_t recalF = 1.;
a82b4462 375 Int_t bc = GetReader()->GetInputEvent()->GetBunchCrossNumber();
649b825d 376
377 for (Int_t iCell = 0; iCell < cells->GetNumberOfCells(); iCell++) {
378 if(GetDebug() > 2)
379 printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - Cell : amp %f, absId %d \n", cells->GetAmplitude(iCell), cells->GetCellNumber(iCell));
380 Int_t nModule = GetModuleNumberCellIndexes(cells->GetCellNumber(iCell),fCalorimeter, icol, irow, iRCU);
381 if(GetDebug() > 2)
382 printf("\t module %d, column %d, row %d \n", nModule,icol,irow);
383
f1538a5f 384 if(nModule < fNModules)
385 {
649b825d 386 //Check if the cell is a bad channel
387 if(GetCaloUtils()->IsBadChannelsRemovalSwitchedOn()){
ae2c2bc4 388 if(fCalorimeter=="EMCAL")
389 {
649b825d 390 if(GetCaloUtils()->GetEMCALChannelStatus(nModule,icol,irow)) continue;
391 }
ae2c2bc4 392 else
393 {
394 if(GetCaloUtils()->GetPHOSChannelStatus(nModule,icol,irow) ) continue;
649b825d 395 }
396 } // use bad channel map
397
398 amp = cells->GetAmplitude(iCell)*recalF;
399 time = cells->GetTime(iCell);
400 id = cells->GetCellNumber(iCell);
401
402 // Amplitude recalibration if set
dbba06ca 403 GetCaloUtils()->RecalibrateCellAmplitude(amp, fCalorimeter, id);
649b825d 404
405 // Time recalibration if set
dbba06ca 406 GetCaloUtils()->RecalibrateCellTime (time, fCalorimeter, id, GetReader()->GetInputEvent()->GetBunchCrossNumber());
649b825d 407
408 //Transform time to ns
409 time *= 1.0e9;
f15c25da 410
e6fec6f5 411 if(time < fTimeCutMin || time > fTimeCutMax)
a87e069d 412 {
822cc7aa 413 if(GetDebug() > 0 )
f1538a5f 414 printf("AliAnaCalorimeterQA - Remove cell with Time %f\n",time);
649b825d 415 continue;
e6fec6f5 416 }
06f1b12a 417
418 // Remove exotic cells, defined only for EMCAL
419 if(fCalorimeter=="EMCAL" &&
420 GetCaloUtils()->GetEMCALRecoUtils()->IsExoticCell(id, cells, bc)) continue;
421
649b825d 422
423 fhAmplitude->Fill(amp);
424 fhAmpId ->Fill(amp,id);
425 fhAmpMod ->Fill(amp,nModule);
426
427 if ((fCalorimeter=="EMCAL" && amp > fEMCALCellAmpMin) ||
2747966a 428 (fCalorimeter=="PHOS" && amp > fPHOSCellAmpMin ) )
429 {
430
431 //E cross for exotic cells
432 if(amp > 0.01) fhCellECross->Fill(amp,1-GetECross(id,cells)/amp);
e3300002 433
649b825d 434 nCellsInModule[nModule]++ ;
06f1b12a 435
649b825d 436 Int_t icols = icol;
437 Int_t irows = irow;
57d8227a 438
439 if(fCalorimeter=="EMCAL")
440 {
649b825d 441 icols = (nModule % 2) ? icol + fNMaxCols : icol;
57d8227a 442 if(nModule < 10 )
443 irows = irow + fNMaxRows * Int_t(nModule / 2);
444 else // 1/3 SM
445 irows = irow + (fNMaxRows / 3) * Int_t(nModule / 2);
649b825d 446 }
57d8227a 447 else
448 {
06f1b12a 449 irows = irow + fNMaxRows * nModule;
649b825d 450 }
06f1b12a 451
649b825d 452 fhGridCells ->Fill(icols,irows);
453 fhGridCellsE->Fill(icols,irows,amp);
454
e6fec6f5 455 if(fFillAllCellTimeHisto)
a87e069d 456 {
649b825d 457 //printf("%s: time %g\n",fCalorimeter.Data(), time);
1a72f6c5 458
459 Double_t v[3] = {0,0,0}; //vertex ;
460 GetReader()->GetVertex(v);
461 if(amp > 0.5) fhTimeVz ->Fill(TMath::Abs(v[2]),time);
462
649b825d 463 fhTime ->Fill(time);
464 fhTimeId ->Fill(time,id);
465 fhTimeAmp ->Fill(amp,time);
466 fhGridCellsTime->Fill(icols,irows,time);
467 fhTimeMod ->Fill(time,nModule);
468 fhTimeAmpPerRCU [nModule*fNRCU+iRCU]->Fill(amp, time);
469
470 }
471 }
472
473 //Get Eta-Phi position of Cell
474 if(fFillAllPosHisto)
475 {
476 if(fCalorimeter=="EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet()){
477 Float_t celleta = 0.;
478 Float_t cellphi = 0.;
479 GetEMCALGeometry()->EtaPhiFromIndex(id, celleta, cellphi);
480
481 fhEtaPhiAmp->Fill(celleta,cellphi,amp);
482 Double_t cellpos[] = {0, 0, 0};
483 GetEMCALGeometry()->GetGlobal(id, cellpos);
484 fhXCellE->Fill(cellpos[0],amp) ;
485 fhYCellE->Fill(cellpos[1],amp) ;
486 fhZCellE->Fill(cellpos[2],amp) ;
487 Float_t rcell = TMath::Sqrt(cellpos[0]*cellpos[0]+cellpos[1]*cellpos[1]);//+cellpos[2]*cellpos[2]);
488 fhRCellE->Fill(rcell,amp) ;
489 fhXYZCell->Fill(cellpos[0],cellpos[1],cellpos[2]) ;
490 }//EMCAL Cells
491 else if(fCalorimeter=="PHOS" && GetCaloUtils()->IsPHOSGeoMatrixSet()){
492 TVector3 xyz;
493 Int_t relId[4], module;
494 Float_t xCell, zCell;
495
496 GetPHOSGeometry()->AbsToRelNumbering(id,relId);
497 module = relId[0];
498 GetPHOSGeometry()->RelPosInModule(relId,xCell,zCell);
499 GetPHOSGeometry()->Local2Global(module,xCell,zCell,xyz);
500 Float_t rcell = TMath::Sqrt(xyz.X()*xyz.X()+xyz.Y()*xyz.Y());
501 fhXCellE ->Fill(xyz.X(),amp) ;
502 fhYCellE ->Fill(xyz.Y(),amp) ;
503 fhZCellE ->Fill(xyz.Z(),amp) ;
504 fhRCellE ->Fill(rcell ,amp) ;
505 fhXYZCell->Fill(xyz.X(),xyz.Y(),xyz.Z()) ;
506 }//PHOS cells
507 }//fill cell position histograms
508
509 if (fCalorimeter=="EMCAL" && amp > fEMCALCellAmpMin) ncells ++ ;
510 else if(fCalorimeter=="PHOS" && amp > fPHOSCellAmpMin) ncells ++ ;
511 //else
512 // printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - no %s CELLS passed the analysis cut\n",fCalorimeter.Data());
513 }//nmodules
514 }//cell loop
515
516 if(ncells > 0 )fhNCells->Fill(ncells) ; //fill the cells after the cut
517
518 //Number of cells per module
519 for(Int_t imod = 0; imod < fNModules; imod++ ) {
520
521 if(GetDebug() > 1)
522 printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - module %d calo %s cells %d\n", imod, fCalorimeter.Data(), nCellsInModule[imod]);
523
524 fhNCellsMod->Fill(nCellsInModule[imod],imod) ;
525
526 }
527
528 delete [] nCellsInModule;
529
530}
531
532//__________________________________________________________________________
533void AliAnaCalorimeterQA::CellInClusterPositionHistograms(AliVCluster* clus)
534{
535 // Fill histograms releated to cell position
536
537
538 Int_t nCaloCellsPerCluster = clus->GetNCells();
539 UShort_t * indexList = clus->GetCellsAbsId();
540 Float_t pos[3];
541 clus->GetPosition(pos);
542 Float_t clEnergy = clus->E();
543
544 //Loop on cluster cells
545 for (Int_t ipos = 0; ipos < nCaloCellsPerCluster; ipos++) {
546
547 // printf("Index %d\n",ipos);
548 Int_t absId = indexList[ipos];
549
550 //Get position of cell compare to cluster
551
552 if(fCalorimeter=="EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet()){
553
554 Double_t cellpos[] = {0, 0, 0};
555 GetEMCALGeometry()->GetGlobal(absId, cellpos);
556
557 fhDeltaCellClusterXNCells->Fill(pos[0]-cellpos[0],nCaloCellsPerCluster) ;
558 fhDeltaCellClusterYNCells->Fill(pos[1]-cellpos[1],nCaloCellsPerCluster) ;
559 fhDeltaCellClusterZNCells->Fill(pos[2]-cellpos[2],nCaloCellsPerCluster) ;
560
561 fhDeltaCellClusterXE->Fill(pos[0]-cellpos[0],clEnergy) ;
562 fhDeltaCellClusterYE->Fill(pos[1]-cellpos[1],clEnergy) ;
563 fhDeltaCellClusterZE->Fill(pos[2]-cellpos[2],clEnergy) ;
564
565 Float_t r = TMath::Sqrt(pos[0] *pos[0] + pos[1] * pos[1] );
566 Float_t rcell = TMath::Sqrt(cellpos[0]*cellpos[0] + cellpos[1]* cellpos[1]);
567
568 fhDeltaCellClusterRNCells->Fill(r-rcell, nCaloCellsPerCluster) ;
569 fhDeltaCellClusterRE ->Fill(r-rcell, clEnergy) ;
570
571 }//EMCAL and its matrices are available
572 else if(fCalorimeter=="PHOS" && GetCaloUtils()->IsPHOSGeoMatrixSet()){
573 TVector3 xyz;
574 Int_t relId[4], module;
575 Float_t xCell, zCell;
576
577 GetPHOSGeometry()->AbsToRelNumbering(absId,relId);
578 module = relId[0];
579 GetPHOSGeometry()->RelPosInModule(relId,xCell,zCell);
580 GetPHOSGeometry()->Local2Global(module,xCell,zCell,xyz);
581
582 fhDeltaCellClusterXNCells->Fill(pos[0]-xyz.X(),nCaloCellsPerCluster) ;
583 fhDeltaCellClusterYNCells->Fill(pos[1]-xyz.Y(),nCaloCellsPerCluster) ;
584 fhDeltaCellClusterZNCells->Fill(pos[2]-xyz.Z(),nCaloCellsPerCluster) ;
585
586 fhDeltaCellClusterXE->Fill(pos[0]-xyz.X(),clEnergy) ;
587 fhDeltaCellClusterYE->Fill(pos[1]-xyz.Y(),clEnergy) ;
588 fhDeltaCellClusterZE->Fill(pos[2]-xyz.Z(),clEnergy) ;
589
590 Float_t r = TMath::Sqrt(pos[0] * pos[0] + pos[1] * pos[1] );
591 Float_t rcell = TMath::Sqrt(xyz.X() * xyz.X() + xyz.Y() * xyz.Y());
592
593 fhDeltaCellClusterRNCells->Fill(r-rcell, nCaloCellsPerCluster) ;
594 fhDeltaCellClusterRE ->Fill(r-rcell, clEnergy) ;
595
596 }//PHOS and its matrices are available
597 }// cluster cell loop
598}
599
600//___________________________________________________________________________________________
1a83b960 601void AliAnaCalorimeterQA::ClusterAsymmetryHistograms(AliVCluster* clus, const Int_t absIdMax,
602 const Bool_t goodCluster)
649b825d 603{
604 // Study the shape of the cluster in cell units terms
605
606 //No use to study clusters with less than 4 cells
607 if(clus->GetNCells() <=3 ) return;
608
609 Int_t dIeta = 0;
610 Int_t dIphi = 0;
611
612 Int_t ietaMax=-1; Int_t iphiMax = 0; Int_t rcuMax = 0;
613 Int_t smMax = GetModuleNumberCellIndexes(absIdMax,fCalorimeter, ietaMax, iphiMax, rcuMax);
614
615 for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++) {
616
617 Int_t absId = clus->GetCellsAbsId()[ipos];
618
619 Int_t ieta=-1; Int_t iphi = 0; Int_t rcu = 0;
620 Int_t sm = GetModuleNumberCellIndexes(absId,fCalorimeter, ieta, iphi, rcu);
621
622 if(dIphi < TMath::Abs(iphi-iphiMax)) dIphi = TMath::Abs(iphi-iphiMax);
623
624 if(smMax==sm){
625 if(dIeta < TMath::Abs(ieta-ietaMax)) dIeta = TMath::Abs(ieta-ietaMax);
626 }
627 else {
628 Int_t ietaShift = ieta;
629 Int_t ietaMaxShift = ietaMax;
630 if (ieta > ietaMax) ietaMaxShift+=48;
631 else ietaShift +=48;
632 if(dIeta < TMath::Abs(ietaShift-ietaMaxShift)) dIeta = TMath::Abs(ietaShift-ietaMaxShift);
633 }
634
649b825d 635 }// fill cell-cluster histogram loop
636
649b825d 637
1a83b960 638 Float_t dIA = 1.*(dIphi-dIeta)/(dIeta+dIphi);
639
640 if(goodCluster)
641 {
649b825d 642
1a83b960 643 // Was cluster matched?
644 Bool_t matched = GetCaloPID()->IsTrackMatched(clus,GetCaloUtils(),GetReader()->GetInputEvent());
649b825d 645
1a83b960 646 if (clus->E() < 2 ) fhDeltaIEtaDeltaIPhiE0[matched]->Fill(dIeta,dIphi);
647 else if(clus->E() < 6 ) fhDeltaIEtaDeltaIPhiE2[matched]->Fill(dIeta,dIphi);
648 else fhDeltaIEtaDeltaIPhiE6[matched]->Fill(dIeta,dIphi);
649
650 fhDeltaIA[matched]->Fill(clus->E(),dIA);
651
652 if(clus->E() > 0.5){
653
654 fhDeltaIAL0[matched] ->Fill(clus->GetM02(),dIA);
655 fhDeltaIAL1[matched] ->Fill(clus->GetM20(),dIA);
656 fhDeltaIANCells[matched]->Fill(clus->GetNCells(),dIA);
657
649b825d 658 }
659
1a83b960 660 // Origin of clusters
661 Int_t nLabel = clus->GetNLabels();
662 Int_t* labels = clus->GetLabels();
663 if(IsDataMC()){
664 Int_t tag = GetMCAnalysisUtils()->CheckOrigin(labels,nLabel, GetReader(),0);
665 if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton) &&
666 !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) &&
667 !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta) &&
668 !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion) ){
669 fhDeltaIAMC[0]->Fill(clus->E(),dIA);//Pure Photon
670 }
671 else if ( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCElectron) &&
672 !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion) ){
673 fhDeltaIAMC[1]->Fill(clus->E(),dIA);//Pure electron
674 }
675 else if ( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton) &&
676 GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion) ){
677 fhDeltaIAMC[2]->Fill(clus->E(),dIA);//Converted cluster
678 }
679 else if(!GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton)){
680 fhDeltaIAMC[3]->Fill(clus->E(),dIA);//Hadrons
681 }
682
683 } // MC
684 } // good cluster
685 else
686 {
687 if (clus->E() < 2 ) fhBadClusterDeltaIEtaDeltaIPhiE0->Fill(dIeta,dIphi);
688 else if(clus->E() < 6 ) fhBadClusterDeltaIEtaDeltaIPhiE2->Fill(dIeta,dIphi);
689 else fhBadClusterDeltaIEtaDeltaIPhiE6->Fill(dIeta,dIphi);
690
691 fhBadClusterDeltaIA->Fill(clus->E(),dIA);
692
693 }
649b825d 694}
695
1a83b960 696//_________________________________________________________________________________________________________________________
c5693f62 697void AliAnaCalorimeterQA::ClusterHistograms(AliVCluster* clus,const TObjArray *caloClusters, AliVCaloCells * cells,
1a83b960 698 const Int_t absIdMax, const Double_t maxCellFraction, const Float_t eCrossFrac,
a82b4462 699 const Double_t tmax, Double_t timeAverages[2])
649b825d 700{
701 //Fill CaloCluster related histograms
702
1a83b960 703 Double_t tof = clus->GetTOF()*1.e9;
649b825d 704
1a83b960 705 fhLambda0 ->Fill(clus->E(),clus->GetM02());
706 fhLambda1 ->Fill(clus->E(),clus->GetM20());
707 fhDispersion ->Fill(clus->E(),clus->GetDispersion());
649b825d 708
a82b4462 709 fhClusterMaxCellDiff ->Fill(clus->E(),maxCellFraction);
710 fhClusterMaxCellECross->Fill(clus->E(),eCrossFrac);
711 fhClusterTimeEnergy ->Fill(clus->E(),tof);
649b825d 712
1a83b960 713 if(fStudyClustersAsymmetry) ClusterAsymmetryHistograms(clus,absIdMax,kTRUE);
714
649b825d 715 //Clusters in event time difference
716 for(Int_t iclus2 = 0; iclus2 < caloClusters->GetEntriesFast(); iclus2++ ){
717
718 AliVCluster* clus2 = (AliVCluster*) caloClusters->At(iclus2);
719
720 if(clus->GetID()==clus2->GetID()) continue;
721
a87e069d 722 if(clus->GetM02() > 0.01 && clus2->GetM02() > 0.01)
723 {
a82b4462 724 Double_t tof2 = clus2->GetTOF()*1.e9;
649b825d 725 fhClusterPairDiffTimeE ->Fill(clus->E(), tof-tof2);
726 }
727 }
728
1a83b960 729 Int_t nModule = GetModuleNumber(clus);
730 Int_t nCaloCellsPerCluster = clus->GetNCells();
731
649b825d 732 if(nCaloCellsPerCluster > 1){
733
734 // check time of cells respect to max energy cell
735
e6fec6f5 736 if(fFillAllCellTimeHisto)
a87e069d 737 {
649b825d 738 fhClusterMaxCellDiffAverageTime ->Fill(clus->E(),tmax-timeAverages[0]);
649b825d 739 fhClusterMaxCellDiffWeightedTime ->Fill(clus->E(),tmax-timeAverages[1]);
649b825d 740 }
741
dbba06ca 742 for (Int_t ipos = 0; ipos < nCaloCellsPerCluster; ipos++)
743 {
649b825d 744 Int_t absId = clus->GetCellsAbsId()[ipos];
2747966a 745 if(absId == absIdMax || cells->GetCellAmplitude(absIdMax) < 0.01) continue;
649b825d 746
747 Float_t frac = cells->GetCellAmplitude(absId)/cells->GetCellAmplitude(absIdMax);
748 fhClusterMaxCellCloseCellRatio->Fill(clus->E(),frac);
749 fhClusterMaxCellCloseCellDiff ->Fill(clus->E(),cells->GetCellAmplitude(absIdMax)-cells->GetCellAmplitude(absId));
750
e6fec6f5 751 if(fFillAllCellTimeHisto)
dbba06ca 752 {
649b825d 753 Double_t time = cells->GetCellTime(absId);
dbba06ca 754 GetCaloUtils()->RecalibrateCellTime(time, fCalorimeter, absId,GetReader()->GetInputEvent()->GetBunchCrossNumber());
649b825d 755
756 Float_t diff = (tmax-time*1.0e9);
757 fhCellTimeSpreadRespectToCellMax->Fill(clus->E(), diff);
758 if(TMath::Abs(TMath::Abs(diff) > 100) && clus->E() > 1 ) fhCellIdCellLargeTimeSpread->Fill(absId);
759 }
760
761 }// fill cell-cluster histogram loop
762
763 }//check time and energy of cells respect to max energy cell if cluster of more than 1 cell
764
765
766 // Get vertex for photon momentum calculation and event selection
767 Double_t v[3] = {0,0,0}; //vertex ;
1a83b960 768 //GetReader()->GetVertex(v); //
649b825d 769
770 TLorentzVector mom ;
1a83b960 771 clus->GetMomentum(mom,v);
649b825d 772
773 Float_t e = mom.E();
774 Float_t pt = mom.Pt();
775 Float_t eta = mom.Eta();
776 Float_t phi = mom.Phi();
777 if(phi < 0) phi +=TMath::TwoPi();
778
779 if(GetDebug() > 0) {
780 printf("AliAnaCalorimeterQA::ClusterHistograms() - cluster: E %2.3f, pT %2.3f, eta %2.3f, phi %2.3f \n",e,pt,eta,phi*TMath::RadToDeg());
781 }
782
783 fhE ->Fill(e);
784 if(nModule >=0 && nModule < fNModules) fhEMod->Fill(e,nModule);
785 if(fFillAllTH12){
786 fhPt ->Fill(pt);
787 fhPhi ->Fill(phi);
788 fhEta ->Fill(eta);
789 }
790
791 if(fFillAllTH3)
792 fhEtaPhiE->Fill(eta,phi,e);
793
794 //Cells per cluster
795 fhNCellsPerCluster ->Fill(e, nCaloCellsPerCluster);
a82b4462 796
649b825d 797 //Position
798 if(fFillAllPosHisto2){
799
800 Float_t pos[3] ;
801 clus->GetPosition(pos);
802
803 fhXE ->Fill(pos[0],e);
804 fhYE ->Fill(pos[1],e);
805 fhZE ->Fill(pos[2],e);
806 if(fFillAllTH3)
807 fhXYZ ->Fill(pos[0], pos[1],pos[2]);
808
809 fhXNCells->Fill(pos[0],nCaloCellsPerCluster);
810 fhYNCells->Fill(pos[1],nCaloCellsPerCluster);
811 fhZNCells->Fill(pos[2],nCaloCellsPerCluster);
812 Float_t rxyz = TMath::Sqrt(pos[0]*pos[0]+pos[1]*pos[1]);//+pos[2]*pos[2]);
813 fhRE ->Fill(rxyz,e);
814 fhRNCells->Fill(rxyz ,nCaloCellsPerCluster);
815 }
816
817 if(nModule >=0 && nModule < fNModules) fhNCellsPerClusterMod[nModule]->Fill(e, nCaloCellsPerCluster);
818
819}
820
f3138ecf 821//____________________________________________________________________________
822void AliAnaCalorimeterQA::ClusterLoopHistograms(const TObjArray *caloClusters,
823 AliVCaloCells* cells)
649b825d 824{
825 // Fill clusters related histograms
826
827 TLorentzVector mom ;
828 Int_t nLabel = 0 ;
829 Int_t *labels = 0x0;
830 Int_t nCaloClusters = caloClusters->GetEntriesFast() ;
831 Int_t nCaloClustersAccepted = 0 ;
832 Int_t nCaloCellsPerCluster = 0 ;
833 Bool_t matched = kFALSE;
834 Int_t nModule =-1 ;
835
836 // Get vertex for photon momentum calculation and event selection
837 Double_t v[3] = {0,0,0}; //vertex ;
1a83b960 838 //GetReader()->GetVertex(v);
649b825d 839
840 Int_t *nClustersInModule = new Int_t[fNModules];
841 for(Int_t imod = 0; imod < fNModules; imod++ ) nClustersInModule[imod] = 0;
842
843 if(GetDebug() > 0)
844 printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - In %s there are %d clusters \n", fCalorimeter.Data(), nCaloClusters);
845
846 // Loop over CaloClusters
847 for(Int_t iclus = 0; iclus < nCaloClusters; iclus++){
848
1a72f6c5 849 if(GetDebug() > 0)
850 printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - cluster: %d/%d, data %d \n",
f3138ecf 851 iclus+1,nCaloClusters,GetReader()->GetDataType());
649b825d 852
853 AliVCluster* clus = (AliVCluster*)caloClusters->At(iclus);
854
855 // Get the fraction of the cluster energy that carries the cell with highest energy and its absId
856 Float_t maxCellFraction = 0.;
857 Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(cells, clus,maxCellFraction);
858
859 //Cut on time of clusters
860 Double_t tof = clus->GetTOF()*1.e9;
e6fec6f5 861 if(tof < fTimeCutMin || tof > fTimeCutMax)
862 {
649b825d 863 if(GetDebug() > 0 )printf("AliAnaCalorimeterQA - Remove cluster with TOF %f\n",tof);
864 continue;
865 }
866
867 // Get cluster kinematics
1a83b960 868 clus->GetMomentum(mom,v);
649b825d 869
870 // Check only certain regions
871 Bool_t in = kTRUE;
872 if(IsFiducialCutOn()) in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
873 if(!in) continue;
874
1a83b960 875 // MC labels
649b825d 876 nLabel = clus->GetNLabels();
877 labels = clus->GetLabels();
878
1a83b960 879 // SuperModule number of cluster
880 nModule = GetModuleNumber(clus);
881
649b825d 882 // Cells per cluster
883 nCaloCellsPerCluster = clus->GetNCells();
884
885 // Cluster mathed with track?
49b5c49b 886 matched = GetCaloPID()->IsTrackMatched(clus,GetCaloUtils(), GetReader()->GetInputEvent());
649b825d 887
888 // Get some time averages
889 Double_t averTime[4] = {0.,0.,0.,0.};
890 CalculateAverageTime(clus, cells, averTime);
891
892 //Get time of max cell
893 Double_t tmax = cells->GetCellTime(absIdMax);
dbba06ca 894 GetCaloUtils()->RecalibrateCellTime(tmax, fCalorimeter, absIdMax,GetReader()->GetInputEvent()->GetBunchCrossNumber());
649b825d 895 tmax*=1.e9;
896
1a83b960 897 // Fill histograms related to single cluster
898
899
900 // Fill some histograms before applying the exotic cell / bad map cut
901 fhNCellsPerClusterNoCut ->Fill(clus->E(), nCaloCellsPerCluster);
902 if(nModule >=0 && nModule < fNModules) fhNCellsPerClusterModNoCut[nModule]->Fill(clus->E(), nCaloCellsPerCluster);
903
904 fhClusterMaxCellDiffNoCut->Fill(clus->E(),maxCellFraction);
905
906 Float_t ampMax = cells->GetCellAmplitude(absIdMax);
dbba06ca 907 GetCaloUtils()->RecalibrateCellAmplitude(ampMax,fCalorimeter, absIdMax);
1a83b960 908
f1538a5f 909 if(fStudyExotic) ExoticHistograms(absIdMax, ampMax, clus, cells);
910
649b825d 911 //Check bad clusters if requested and rejection was not on
a82b4462 912 Bool_t goodCluster = IsGoodCluster(absIdMax, cells);
f1538a5f 913
2747966a 914 Float_t eCrossFrac = 0;
915 if(ampMax > 0.01) eCrossFrac = 1-GetECross(absIdMax,cells)/ampMax;
649b825d 916
649b825d 917 if(!goodCluster)
1a83b960 918 {
649b825d 919 BadClusterHistograms(clus, caloClusters, cells, absIdMax,
1a83b960 920 maxCellFraction, eCrossFrac, tmax, averTime);
921 continue;
922 }
649b825d 923
924 ClusterHistograms(clus, caloClusters, cells, absIdMax,
1a83b960 925 maxCellFraction, eCrossFrac, tmax, averTime);
649b825d 926
927 nCaloClustersAccepted++;
49214ef9 928 nModule = GetModuleNumber(clus);
2747966a 929 if(nModule >=0 && nModule < fNModules)
930 {
649b825d 931 if (fCalorimeter=="EMCAL" && mom.E() > 2*fEMCALCellAmpMin) nClustersInModule[nModule]++;
932 else if(fCalorimeter=="PHOS" && mom.E() > 2*fPHOSCellAmpMin ) nClustersInModule[nModule]++;
933 }
1a83b960 934
649b825d 935 // Cluster weights
936 if(fStudyWeight) WeightHistograms(clus, cells);
937
938 // Cells in cluster position
939 if(fFillAllPosHisto) CellInClusterPositionHistograms(clus);
940
941 // Fill histograms related to single cluster, mc vs data
942 Int_t mcOK = kFALSE;
943 Int_t pdg = -1;
944 if(IsDataMC() && nLabel > 0 && labels)
945 mcOK = ClusterMCHistograms(mom, matched, labels, nLabel, pdg);
008693e5 946
649b825d 947 // Matched clusters with tracks, also do some MC comparison, needs input from ClusterMCHistograms
948 if( matched && fFillAllTMHisto)
949 ClusterMatchedWithTrackHistograms(clus,mom,mcOK,pdg);
950
951 // Invariant mass
d07278cf 952 // Try to reduce background with a mild shower shape cut and no more than 1 maxima
953 // in cluster and remove low energy clusters
954 if(fFillAllPi0Histo && nCaloClusters > 1 && nCaloCellsPerCluster > 1 &&
955 GetCaloUtils()->GetNumberOfLocalMaxima(clus,cells) == 1 &&
956 clus->GetM02() < 0.5 && clus->E() > 0.3)
a82b4462 957 InvariantMassHistograms(iclus, mom, nModule, caloClusters,cells);
649b825d 958
959 }//cluster loop
960
961 // Number of clusters histograms
962 if(nCaloClustersAccepted > 0) fhNClusters->Fill(nCaloClustersAccepted);
963
964 // Number of clusters per module
2747966a 965 for(Int_t imod = 0; imod < fNModules; imod++ )
966 {
649b825d 967 if(GetDebug() > 1)
968 printf("AliAnaCalorimeterQA::ClusterLoopHistograms() - module %d calo %s clusters %d\n", imod, fCalorimeter.Data(), nClustersInModule[imod]);
969 fhNClustersMod->Fill(nClustersInModule[imod],imod);
970 }
971
972 delete [] nClustersInModule;
973
974}
975
d07278cf 976//______________________________________________________________________________________________________
649b825d 977Bool_t AliAnaCalorimeterQA::ClusterMCHistograms(const TLorentzVector mom, const Bool_t matched,
978 const Int_t * labels, const Int_t nLabels, Int_t & pdg )
979{
980
981 //Fill histograms only possible when simulation
982
2747966a 983 if(!labels || nLabels<=0)
984 {
42d47cb7 985 if(GetDebug() > 1) printf("AliAnaCalorimeterQA::ClusterMCHistograms() - Strange, labels array %p, n labels %d \n", labels,nLabels);
986 return kFALSE;
987 }
988
2747966a 989 if(GetDebug() > 1)
990 {
42d47cb7 991 printf("AliAnaCalorimeterQA::ClusterMCHistograms() - Primaries: nlabels %d\n",nLabels);
649b825d 992 }
993
994 Float_t e = mom.E();
995 Float_t eta = mom.Eta();
996 Float_t phi = mom.Phi();
997 if(phi < 0) phi +=TMath::TwoPi();
998
999 AliAODMCParticle * aodprimary = 0x0;
1000 TParticle * primary = 0x0;
1001
1002 //Play with the MC stack if available
1003 Int_t label = labels[0];
1004
2747966a 1005 if(label < 0)
1006 {
008693e5 1007 if(GetDebug() >= 0) printf("AliAnaCalorimeterQA::ClusterMCHistograms() *** bad label ***: label %d \n", label);
649b825d 1008 return kFALSE;
1009 }
1010
1011 Int_t pdg0 =-1;Int_t status = -1; Int_t iMother = -1; Int_t iParent = -1;
1012 Float_t vxMC= 0; Float_t vyMC = 0;
1013 Float_t eMC = 0; Float_t ptMC= 0; Float_t phiMC =0; Float_t etaMC = 0;
1014 Int_t charge = 0;
1015
1016 //Check the origin.
1017 Int_t tag = GetMCAnalysisUtils()->CheckOrigin(labels,nLabels, GetReader(),0);
1018
d07278cf 1019 if ( GetReader()->ReadStack() &&
1020 !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCUnknown))
1021 { //it MC stack and known tag
649b825d 1022
d07278cf 1023 if( label >= GetMCStack()->GetNtrack())
1024 {
008693e5 1025 if(GetDebug() >= 0) printf("AliAnaCalorimeterQA::ClusterMCHistograms() *** large label ***: label %d, n tracks %d \n", label, GetMCStack()->GetNtrack());
649b825d 1026 return kFALSE;
1027 }
1028
1029 primary = GetMCStack()->Particle(label);
1030 iMother = label;
1031 pdg0 = TMath::Abs(primary->GetPdgCode());
1032 pdg = pdg0;
1033 status = primary->GetStatusCode();
1034 vxMC = primary->Vx();
1035 vyMC = primary->Vy();
1036 iParent = primary->GetFirstMother();
1037
d07278cf 1038 if(GetDebug() > 1 )
1039 {
008693e5 1040 printf("AliAnaCalorimeterQA::ClusterMCHistograms() - Cluster most contributing mother: \n");
649b825d 1041 printf("\t Mother label %d, pdg %d, %s, status %d, parent %d \n",iMother, pdg0, primary->GetName(),status, iParent);
1042 }
1043
1044 //Get final particle, no conversion products
d07278cf 1045 if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion))
1046 {
649b825d 1047 //Get the parent
1048 primary = GetMCStack()->Particle(iParent);
1049 pdg = TMath::Abs(primary->GetPdgCode());
008693e5 1050
1051 if(GetDebug() > 1 ) printf("AliAnaCalorimeterQA::ClusterMCHistograms() - Converted cluster!. Find before conversion: \n");
1052
d07278cf 1053 while((pdg == 22 || pdg == 11) && status != 1)
1054 {
008693e5 1055 Int_t iMotherOrg = iMother;
649b825d 1056 iMother = iParent;
1057 primary = GetMCStack()->Particle(iMother);
1058 status = primary->GetStatusCode();
649b825d 1059 pdg = TMath::Abs(primary->GetPdgCode());
008693e5 1060 iParent = primary->GetFirstMother();
1061
1062 // If gone too back and non stable, assign the decay photon/electron
1063 // there are other possible decays, ignore them for the moment
1064 if(pdg==111 || pdg==221)
1065 {
1066 primary = GetMCStack()->Particle(iMotherOrg);
1067 break;
1068 }
1069
1070 if( iParent < 0 )
1071 {
1072 iParent = iMother;
008693e5 1073 break;
1074 }
1075
649b825d 1076 if(GetDebug() > 1 )printf("\t pdg %d, index %d, %s, status %d \n",pdg, iMother, primary->GetName(),status);
1077 }
008693e5 1078
d07278cf 1079 if(GetDebug() > 1 )
1080 {
649b825d 1081 printf("AliAnaCalorimeterQA::ClusterHistograms() - Converted Cluster mother before conversion: \n");
1082 printf("\t Mother label %d, pdg %d, %s, status %d, parent %d \n",iMother, pdg, primary->GetName(), status, iParent);
1083 }
1084
1085 }
1086
1087 //Overlapped pi0 (or eta, there will be very few), get the meson
1088 if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) ||
d07278cf 1089 GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta))
1090 {
649b825d 1091 if(GetDebug() > 1 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped Meson decay!, Find it: \n");
008693e5 1092
d07278cf 1093 while(pdg != 111 && pdg != 221)
008693e5 1094 {
1095 //printf("iMother %d, pdg %d, iParent %d, pdg %d\n",iMother,pdg,iParent,GetMCStack()->Particle(iParent)->GetPdgCode());
649b825d 1096 iMother = iParent;
1097 primary = GetMCStack()->Particle(iMother);
1098 status = primary->GetStatusCode();
649b825d 1099 pdg = TMath::Abs(primary->GetPdgCode());
008693e5 1100 iParent = primary->GetFirstMother();
1101
1102 if( iParent < 0 )break;
1103
649b825d 1104 if(GetDebug() > 1 ) printf("\t pdg %d, %s, index %d\n",pdg, primary->GetName(),iMother);
008693e5 1105
d07278cf 1106 if(iMother==-1)
1107 {
649b825d 1108 printf("AliAnaCalorimeterQA::ClusterHistograms() - Tagged as Overlapped photon but meson not found, why?\n");
1109 //break;
1110 }
1111 }
008693e5 1112
649b825d 1113 if(GetDebug() > 2 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped %s decay, label %d \n",
1114 primary->GetName(),iMother);
1115 }
1116
1117 eMC = primary->Energy();
1118 ptMC = primary->Pt();
1119 phiMC = primary->Phi();
1120 etaMC = primary->Eta();
1121 pdg = TMath::Abs(primary->GetPdgCode());
1122 charge = (Int_t) TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
1123
1124 }
d07278cf 1125 else if( GetReader()->ReadAODMCParticles() &&
1126 !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCUnknown))
1127 {//it MC AOD and known tag
649b825d 1128 //Get the list of MC particles
1129 if(!GetReader()->GetAODMCParticles(0))
1130 AliFatal("MCParticles not available!");
1131
1132 aodprimary = (AliAODMCParticle*) (GetReader()->GetAODMCParticles(0))->At(label);
1133 iMother = label;
1134 pdg0 = TMath::Abs(aodprimary->GetPdgCode());
1135 pdg = pdg0;
1136 status = aodprimary->IsPrimary();
1137 vxMC = aodprimary->Xv();
1138 vyMC = aodprimary->Yv();
1139 iParent = aodprimary->GetMother();
1140
d07278cf 1141 if(GetDebug() > 1 )
1142 {
649b825d 1143 printf("AliAnaCalorimeterQA::ClusterHistograms() - Cluster most contributing mother: \n");
1144 printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d, parent %d \n",
1145 iMother, pdg0, aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary(), iParent);
1146 }
1147
1148 //Get final particle, no conversion products
008693e5 1149 if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion))
1150 {
649b825d 1151 if(GetDebug() > 1 )
1152 printf("AliAnaCalorimeterQA::ClusterHistograms() - Converted cluster!. Find before conversion: \n");
1153 //Get the parent
1154 aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles(0))->At(iParent);
1155 pdg = TMath::Abs(aodprimary->GetPdgCode());
d07278cf 1156 while ((pdg == 22 || pdg == 11) && !aodprimary->IsPhysicalPrimary())
1157 {
008693e5 1158 Int_t iMotherOrg = iMother;
649b825d 1159 iMother = iParent;
1160 aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles(0))->At(iMother);
1161 status = aodprimary->IsPrimary();
1162 iParent = aodprimary->GetMother();
1163 pdg = TMath::Abs(aodprimary->GetPdgCode());
008693e5 1164
1165 // If gone too back and non stable, assign the decay photon/electron
1166 // there are other possible decays, ignore them for the moment
1167 if(pdg==111 || pdg==221)
1168 {
1169 aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles(0))->At(iMotherOrg);
1170 break;
1171 }
1172
1173 if(iParent < 0 )
1174 {
1175 iParent = iMother;
1176 break;
1177 }
1178
649b825d 1179 if(GetDebug() > 1 )
1180 printf("\t pdg %d, index %d, Primary? %d, Physical Primary? %d \n",
1181 pdg, iMother, aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary());
1182 }
1183
d07278cf 1184 if(GetDebug() > 1 )
1185 {
649b825d 1186 printf("AliAnaCalorimeterQA::ClusterHistograms() - Converted Cluster mother before conversion: \n");
1187 printf("\t Mother label %d, pdg %d, parent %d, Primary? %d, Physical Primary? %d \n",
1188 iMother, pdg, iParent, aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary());
1189 }
1190
1191 }
1192
1193 //Overlapped pi0 (or eta, there will be very few), get the meson
1194 if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) ||
2747966a 1195 GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta))
1196 {
649b825d 1197 if(GetDebug() > 1 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped Meson decay!, Find it: PDG %d, mom %d \n",pdg, iMother);
008693e5 1198 while(pdg != 111 && pdg != 221)
1199 {
649b825d 1200 iMother = iParent;
1201 aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles(0))->At(iMother);
1202 status = aodprimary->IsPrimary();
1203 iParent = aodprimary->GetMother();
1204 pdg = TMath::Abs(aodprimary->GetPdgCode());
008693e5 1205
1206 if(iParent < 0 ) break;
649b825d 1207
1208 if(GetDebug() > 1 ) printf("\t pdg %d, index %d\n",pdg, iMother);
1209
2747966a 1210 if(iMother==-1)
1211 {
649b825d 1212 printf("AliAnaCalorimeterQA::ClusterHistograms() - Tagged as Overlapped photon but meson not found, why?\n");
1213 //break;
1214 }
1215 }
1216
1217 if(GetDebug() > 2 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped %s decay, label %d \n",
1218 aodprimary->GetName(),iMother);
1219 }
1220
1221 status = aodprimary->IsPrimary();
1222 eMC = aodprimary->E();
1223 ptMC = aodprimary->Pt();
1224 phiMC = aodprimary->Phi();
1225 etaMC = aodprimary->Eta();
1226 pdg = TMath::Abs(aodprimary->GetPdgCode());
1227 charge = aodprimary->Charge();
1228
1229 }
1230
1231 //Float_t vz = primary->Vz();
1232 Float_t rVMC = TMath::Sqrt(vxMC*vxMC + vyMC*vyMC);
2747966a 1233 if((pdg == 22 || TMath::Abs(pdg)==11) && status!=1)
1234 {
649b825d 1235 fhEMVxyz ->Fill(vxMC,vyMC);//,vz);
1236 fhEMR ->Fill(e,rVMC);
1237 }
1238
1239 //printf("reco e %f, pt %f, phi %f, eta %f \n", e, pt, phi, eta);
1240 //printf("prim e %f, pt %f, phi %f, eta %f \n", eMC,ptMC,phiMC ,etaMC );
1241 //printf("vertex: vx %f, vy %f, vz %f, r %f \n", vxMC, vyMC, vz, r);
1242
1243 //Overlapped pi0 (or eta, there will be very few)
2747966a 1244 if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0))
1245 {
c5693f62 1246 fhRecoMCE [kmcPi0][matched] ->Fill(e,eMC);
1247 if(e > 0.5 && eMC > 0.5) fhRecoMCEta[kmcPi0][(matched)]->Fill(eta,etaMC);
1248 if(e > 0.5 && eMC > 0.5) fhRecoMCPhi[kmcPi0][(matched)]->Fill(phi,phiMC);
1249 if(eMC > 0) fhRecoMCRatioE [kmcPi0][(matched)]->Fill(e,e/eMC);
1250 fhRecoMCDeltaE [kmcPi0][(matched)]->Fill(e,eMC-e);
1251 fhRecoMCDeltaPhi[kmcPi0][(matched)]->Fill(e,phiMC-phi);
1252 fhRecoMCDeltaEta[kmcPi0][(matched)]->Fill(e,etaMC-eta);
649b825d 1253 }//Overlapped pizero decay
2747966a 1254 else if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta))
1255 {
c5693f62 1256 fhRecoMCE [kmcEta][(matched)] ->Fill(e,eMC);
1257 if(e > 0.5 && eMC > 0.5) fhRecoMCEta[kmcEta][(matched)]->Fill(eta,etaMC);
1258 if(e > 0.5 && eMC > 0.5) fhRecoMCPhi[kmcEta][(matched)]->Fill(phi,phiMC);
1259 if(eMC > 0) fhRecoMCRatioE [kmcEta][(matched)]->Fill(e,e/eMC);
1260 fhRecoMCDeltaE [kmcEta][(matched)]->Fill(e,eMC-e);
1261 fhRecoMCDeltaPhi[kmcEta][(matched)]->Fill(e,phiMC-phi);
1262 fhRecoMCDeltaEta[kmcEta][(matched)]->Fill(e,etaMC-eta);
649b825d 1263 }//Overlapped eta decay
008693e5 1264 else if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton) &&
1265 !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion))
2747966a 1266 {
c5693f62 1267 fhRecoMCE [kmcPhoton][(matched)] ->Fill(e,eMC);
1268 if(e > 0.5 && eMC > 0.5) fhRecoMCEta[kmcPhoton][(matched)]->Fill(eta,etaMC);
1269 if(e > 0.5 && eMC > 0.5) fhRecoMCPhi[kmcPhoton][(matched)]->Fill(phi,phiMC);
1270 if(eMC > 0) fhRecoMCRatioE [kmcPhoton][(matched)]->Fill(e,e/eMC);
008693e5 1271
c5693f62 1272 fhRecoMCDeltaE [kmcPhoton][(matched)]->Fill(e,eMC-e);
1273 fhRecoMCDeltaPhi[kmcPhoton][(matched)]->Fill(e,phiMC-phi);
1274 fhRecoMCDeltaEta[kmcPhoton][(matched)]->Fill(e,etaMC-eta);
649b825d 1275 }//photon
008693e5 1276 else if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCElectron) &&
1277 !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion))
2747966a 1278 {
c5693f62 1279 fhRecoMCE [kmcElectron][(matched)] ->Fill(e,eMC);
1280 if(e > 0.5 && eMC > 0.5) fhRecoMCEta[kmcElectron][(matched)]->Fill(eta,etaMC);
1281 if(e > 0.5 && eMC > 0.5) fhRecoMCPhi[kmcElectron][(matched)]->Fill(phi,phiMC);
1282 if(eMC > 0) fhRecoMCRatioE [kmcElectron][(matched)]->Fill(e,e/eMC);
1283 fhRecoMCDeltaE [kmcElectron][(matched)]->Fill(e,eMC-e);
1284 fhRecoMCDeltaPhi[kmcElectron][(matched)]->Fill(e,phiMC-phi);
1285 fhRecoMCDeltaEta[kmcElectron][(matched)]->Fill(e,etaMC-eta);
649b825d 1286 fhEMVxyz ->Fill(vxMC,vyMC);//,vz);
1287 fhEMR ->Fill(e,rVMC);
1288 }
2747966a 1289 else if(charge == 0)
1290 {
c5693f62 1291 fhRecoMCE [kmcNeHadron][(matched)] ->Fill(e,eMC);
1292 if(e > 0.5 && eMC > 0.5) fhRecoMCEta[kmcNeHadron][(matched)]->Fill(eta,etaMC);
1293 if(e > 0.5 && eMC > 0.5) fhRecoMCPhi[kmcNeHadron][(matched)]->Fill(phi,phiMC);
1294 if(eMC > 0) fhRecoMCRatioE [kmcNeHadron][(matched)]->Fill(e,e/eMC);
1295 fhRecoMCDeltaE [kmcNeHadron][(matched)]->Fill(e,eMC-e);
1296 fhRecoMCDeltaPhi[kmcNeHadron][(matched)]->Fill(e,phiMC-phi);
1297 fhRecoMCDeltaEta[kmcNeHadron][(matched)]->Fill(e,etaMC-eta);
649b825d 1298 fhHaVxyz ->Fill(vxMC,vyMC);//,vz);
1299 fhHaR ->Fill(e,rVMC);
1300 }
2747966a 1301 else if(charge!=0)
1302 {
c5693f62 1303 fhRecoMCE [kmcChHadron][(matched)] ->Fill(e,eMC);
1304 if(e > 0.5 && eMC > 0.5) fhRecoMCEta[kmcChHadron][(matched)]->Fill(eta,etaMC);
1305 if(e > 0.5 && eMC > 0.5) fhRecoMCPhi[kmcChHadron][(matched)]->Fill(phi,phiMC);
1306 if(eMC > 0) fhRecoMCRatioE [kmcChHadron][(matched)]->Fill(e,e/eMC);
1307 fhRecoMCDeltaE [kmcChHadron][(matched)]->Fill(e,eMC-e);
1308 fhRecoMCDeltaPhi[kmcChHadron][(matched)]->Fill(e,phiMC-phi);
1309 fhRecoMCDeltaEta[kmcChHadron][(matched)]->Fill(e,etaMC-eta);
649b825d 1310 fhHaVxyz ->Fill(vxMC,vyMC);//,vz);
1311 fhHaR ->Fill(e,rVMC);
1312 }
1313
1314 if(primary || aodprimary) return kTRUE ;
1315 else return kFALSE;
1316
1317}
1318
1319//________________________________________________________________________________________________
1320void AliAnaCalorimeterQA::ClusterMatchedWithTrackHistograms(AliVCluster *clus, TLorentzVector mom,
1321 const Bool_t okPrimary, const Int_t pdg)
1322{
1323 //Histograms for clusters matched with tracks
42d47cb7 1324
649b825d 1325 Float_t e = mom.E();
1326 Float_t pt = mom.Pt();
1327 Float_t eta = mom.Eta();
1328 Float_t phi = mom.Phi();
1329 if(phi < 0) phi +=TMath::TwoPi();
1330
2747966a 1331 if(fFillAllTH12)
1332 {
649b825d 1333 fhECharged ->Fill(e);
1334 fhPtCharged ->Fill(pt);
1335 fhPhiCharged ->Fill(phi);
1336 fhEtaCharged ->Fill(eta);
1337 }
a82b4462 1338
42d47cb7 1339 //Study the track and matched cluster if track exists.
1340
4bfeae64 1341 AliVTrack *track = GetCaloUtils()->GetMatchedTrack(clus, GetReader()->GetInputEvent());
649b825d 1342
1343 if(!track) return ;
649b825d 1344
a87e069d 1345 Double_t tpt = track->Pt();
1346 Double_t tmom = track->P();
1347 Double_t dedx = track->GetTPCsignal();
1348 Int_t nITS = track->GetNcls(0);
1349 Int_t nTPC = track->GetNcls(1);
1350
1351 // Residuals
1352 Float_t deta = clus->GetTrackDz();
1353 Float_t dphi = clus->GetTrackDx();
1354 Double_t dR = TMath::Sqrt(dphi*dphi + deta*deta);
1355
1356 Double_t eOverP = e/tmom;
1357
d55bb5e1 1358 fh1EOverP->Fill(tpt, eOverP);
1359 if(dR < 0.02) fh1EOverPR02->Fill(tpt,eOverP);
a87e069d 1360
1361 fh2dR->Fill(e,dR);
1362 fh2MatchdEdx->Fill(tmom,dedx);
1363
1364 if(IsDataMC() && okPrimary)
1365 {
1366 Double_t charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
649b825d 1367
a87e069d 1368 if(TMath::Abs(pdg) == 11)
1369 {
d55bb5e1 1370 fhMCEle1EOverP->Fill(tpt,eOverP);
a87e069d 1371 fhMCEle1dR->Fill(dR);
1372 fhMCEle2MatchdEdx->Fill(tmom,dedx);
d55bb5e1 1373 if(dR < 0.02) fhMCEle1EOverPR02->Fill(tpt,eOverP);
a87e069d 1374 }
1375 else if(charge!=0)
1376 {
d55bb5e1 1377 fhMCChHad1EOverP->Fill(tpt,eOverP);
a87e069d 1378 fhMCChHad1dR->Fill(dR);
1379 fhMCChHad2MatchdEdx->Fill(tmom,dedx);
d55bb5e1 1380 if(dR < 0.02) fhMCChHad1EOverPR02->Fill(tpt,eOverP);
a87e069d 1381 }
1382 else if(charge == 0)
1383 {
d55bb5e1 1384 fhMCNeutral1EOverP->Fill(tpt,eOverP);
a87e069d 1385 fhMCNeutral1dR->Fill(dR);
1386 fhMCNeutral2MatchdEdx->Fill(tmom,dedx);
d55bb5e1 1387 if(dR < 0.02) fhMCNeutral1EOverPR02->Fill(tpt,eOverP);
649b825d 1388 }
a87e069d 1389 }//DataMC
1390
1391 if(dR < 0.02 && eOverP > 0.6 && eOverP< 1.2
1392 && clus->GetNCells() > 1 && nITS > 3 && nTPC > 20)
1393 {
1394 fh2EledEdx->Fill(tmom,dedx);
649b825d 1395 }
a87e069d 1396
649b825d 1397}
1398
1399//___________________________________
1400void AliAnaCalorimeterQA::Correlate()
1401{
1402 // Correlate information from PHOS and EMCAL and with V0 and track multiplicity
1403
1404 //Clusters
1405 TObjArray * caloClustersEMCAL = GetEMCALClusters();
1406 TObjArray * caloClustersPHOS = GetPHOSClusters();
1407
1408 Int_t nclEMCAL = caloClustersEMCAL->GetEntriesFast();
1409 Int_t nclPHOS = caloClustersPHOS ->GetEntriesFast();
1410
1411 Float_t sumClusterEnergyEMCAL = 0;
1412 Float_t sumClusterEnergyPHOS = 0;
1413 Int_t iclus = 0;
1414 for(iclus = 0 ; iclus < caloClustersEMCAL->GetEntriesFast() ; iclus++)
1415 sumClusterEnergyEMCAL += ((AliVCluster*)caloClustersEMCAL->At(iclus))->E();
1416 for(iclus = 0 ; iclus < caloClustersPHOS->GetEntriesFast(); iclus++)
1417 sumClusterEnergyPHOS += ((AliVCluster*)caloClustersPHOS->At(iclus))->E();
1418
1419
1420 //Cells
1421
1422 AliVCaloCells * cellsEMCAL = GetEMCALCells();
1423 AliVCaloCells * cellsPHOS = GetPHOSCells();
1424
1425 Int_t ncellsEMCAL = cellsEMCAL->GetNumberOfCells();
1426 Int_t ncellsPHOS = cellsPHOS ->GetNumberOfCells();
1427
1428 Float_t sumCellEnergyEMCAL = 0;
1429 Float_t sumCellEnergyPHOS = 0;
1430 Int_t icell = 0;
1431 for(icell = 0 ; icell < cellsEMCAL->GetNumberOfCells() ; icell++)
1432 sumCellEnergyEMCAL += cellsEMCAL->GetAmplitude(icell);
1433 for(icell = 0 ; icell < cellsPHOS->GetNumberOfCells(); icell++)
1434 sumCellEnergyPHOS += cellsPHOS->GetAmplitude(icell);
1435
1436
1437 //Fill Histograms
1438 fhCaloCorrNClusters->Fill(nclEMCAL,nclPHOS);
1439 fhCaloCorrEClusters->Fill(sumClusterEnergyEMCAL,sumClusterEnergyPHOS);
1440 fhCaloCorrNCells ->Fill(ncellsEMCAL,ncellsPHOS);
1441 fhCaloCorrECells ->Fill(sumCellEnergyEMCAL,sumCellEnergyPHOS);
1442
1443 Int_t v0S = GetV0Signal(0)+GetV0Signal(1);
1444 Int_t v0M = GetV0Multiplicity(0)+GetV0Multiplicity(1);
1445 Int_t trM = GetTrackMultiplicity();
1446 if(fCalorimeter=="PHOS"){
1447 fhCaloV0MCorrNClusters ->Fill(v0M,nclPHOS);
1448 fhCaloV0MCorrEClusters ->Fill(v0M,sumClusterEnergyPHOS);
1449 fhCaloV0MCorrNCells ->Fill(v0M,ncellsPHOS);
1450 fhCaloV0MCorrECells ->Fill(v0M,sumCellEnergyPHOS);
1451
1452 fhCaloV0SCorrNClusters ->Fill(v0S,nclPHOS);
1453 fhCaloV0SCorrEClusters ->Fill(v0S,sumClusterEnergyPHOS);
1454 fhCaloV0SCorrNCells ->Fill(v0S,ncellsPHOS);
1455 fhCaloV0SCorrECells ->Fill(v0S,sumCellEnergyPHOS);
1456
1457 fhCaloTrackMCorrNClusters->Fill(trM,nclPHOS);
1458 fhCaloTrackMCorrEClusters->Fill(trM,sumClusterEnergyPHOS);
1459 fhCaloTrackMCorrNCells ->Fill(trM,ncellsPHOS);
1460 fhCaloTrackMCorrECells ->Fill(trM,sumCellEnergyPHOS);
1461 }
1462 else{
1463 fhCaloV0MCorrNClusters ->Fill(v0M,nclEMCAL);
1464 fhCaloV0MCorrEClusters ->Fill(v0M,sumClusterEnergyEMCAL);
1465 fhCaloV0MCorrNCells ->Fill(v0M,ncellsEMCAL);
1466 fhCaloV0MCorrECells ->Fill(v0M,sumCellEnergyEMCAL);
1467
1468 fhCaloV0SCorrNClusters ->Fill(v0S,nclEMCAL);
1469 fhCaloV0SCorrEClusters ->Fill(v0S,sumClusterEnergyEMCAL);
1470 fhCaloV0SCorrNCells ->Fill(v0S,ncellsEMCAL);
1471 fhCaloV0SCorrECells ->Fill(v0S,sumCellEnergyEMCAL);
1472
1473 fhCaloTrackMCorrNClusters->Fill(trM,nclEMCAL);
1474 fhCaloTrackMCorrEClusters->Fill(trM,sumClusterEnergyEMCAL);
1475 fhCaloTrackMCorrNCells ->Fill(trM,ncellsEMCAL);
1476 fhCaloTrackMCorrECells ->Fill(trM,sumCellEnergyEMCAL);
1477 }
1478
1479 if(GetDebug() > 0 )
1480 {
1481 printf("AliAnaCalorimeterQA::Correlate(): \n");
1482 printf("\t EMCAL: N cells %d, N clusters %d, summed E cells %f, summed E clusters %f \n",
1483 ncellsEMCAL,nclEMCAL, sumCellEnergyEMCAL,sumClusterEnergyEMCAL);
1484 printf("\t PHOS : N cells %d, N clusters %d, summed E cells %f, summed E clusters %f \n",
1485 ncellsPHOS,nclPHOS,sumCellEnergyPHOS,sumClusterEnergyPHOS);
1486 printf("\t V0 : Signal %d, Multiplicity %d, Track Multiplicity %d \n", v0S,v0M,trM);
1487 }
1488
1489}
1490
1491//__________________________________________________
1492TObjString * AliAnaCalorimeterQA::GetAnalysisCuts()
1493{
1494 //Save parameters used for analysis
1495 TString parList ; //this will be list of parameters used for this analysis.
1496 const Int_t buffersize = 255;
1497 char onePar[buffersize] ;
1498
1499 snprintf(onePar,buffersize,"--- AliAnaCalorimeterQA ---\n") ;
1500 parList+=onePar ;
1501 snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
1502 parList+=onePar ;
1503 snprintf(onePar,buffersize,"Time Cut : %2.2f < T < %2.2f ns \n",fTimeCutMin, fTimeCutMax) ;
1504 parList+=onePar ;
1505 snprintf(onePar,buffersize,"PHOS Cell Amplitude > %2.2f GeV, EMCAL Cell Amplitude > %2.2f GeV \n",fPHOSCellAmpMin, fEMCALCellAmpMin) ;
1506 parList+=onePar ;
1507 //Get parameters set in base class.
1508 //parList += GetBaseParametersList() ;
1509
1510 //Get parameters set in FiducialCut class (not available yet)
1511 //parlist += GetFidCut()->GetFidCutParametersList()
1512
1513 return new TObjString(parList) ;
1514}
1515
f1538a5f 1516//___________________________________________________________________________________
1517void AliAnaCalorimeterQA::ExoticHistograms(const Int_t absIdMax, const Float_t ampMax,
1518 AliVCluster *clus, AliVCaloCells* cells)
1519{
1520 // Calculate weights
1521
1522 if(ampMax < 0.01)
1523 {
1524 printf("AliAnaCalorimeterQA::ExoticHistograms()- Low amplitude energy %f\n",ampMax);
1525 return;
1526 }
1527
cc11121e 1528 Float_t l0 = clus->GetM02();
765206a5 1529 Float_t l1 = clus->GetM20();
cc11121e 1530 Float_t en = clus->E();
1531 Int_t nc = clus->GetNCells();
765206a5 1532 Double_t tmax = clus->GetTOF()*1.e9; // recalibrated elsewhere
1533
1534 Float_t eCrossFrac = 1-GetECross(absIdMax,cells, 10000000)/ampMax;
1535
1536 if(en > 5)
1537 {
1538 fhExoL0ECross->Fill(eCrossFrac,l0);
1539 fhExoL1ECross->Fill(eCrossFrac,l1);
1540 }
f1538a5f 1541
1542 for(Int_t ie = 0; ie < fExoNECrossCuts; ie++)
1543 {
1544 for(Int_t idt = 0; idt < fExoNDTimeCuts; idt++)
1545 {
765206a5 1546 eCrossFrac = 1-GetECross(absIdMax,cells, fExoDTimeCuts[idt])/ampMax;
f1538a5f 1547
1548 if(eCrossFrac > fExoECrossCuts[ie])
1549 {
1550 //Exotic
1551 fhExoL0 [ie][idt]->Fill(en,l0 );
765206a5 1552 fhExoL1 [ie][idt]->Fill(en,l1 );
1553 fhExoTime [ie][idt]->Fill(en,tmax);
1554
1555 if(en > 5)
1556 {
1557 fhExoL0NCell[ie][idt]->Fill(nc,l0);
1558 fhExoL1NCell[ie][idt]->Fill(nc,l1);
1559 }
f1538a5f 1560
1561 // Diff time, do for one cut in e cross
1562 if(ie == 0)
1563 {
1564 for (Int_t icell = 0; icell < clus->GetNCells(); icell++)
1565 {
1566 Int_t absId = clus->GetCellsAbsId()[icell];
1567 Double_t time = cells->GetCellTime(absId);
1568 GetCaloUtils()->RecalibrateCellTime(time, fCalorimeter, absId,GetReader()->GetInputEvent()->GetBunchCrossNumber());
1569
1570 Float_t diff = (tmax-time)*1e9;
1571 fhExoDTime[idt]->Fill(en, diff);
1572 }
1573 }
1574 }
1575 else
1576 {
1577 fhExoECross[ie][idt]->Fill(en,eCrossFrac);
cc11121e 1578 fhExoNCell [ie][idt]->Fill(en,nc);
f1538a5f 1579 }
1580 } // D time cut loop
1581 } // e cross cut loop
1582}
1583
649b825d 1584//____________________________________________________
1585TList * AliAnaCalorimeterQA::GetCreateOutputObjects()
1586{
1587 // Create histograms to be saved in output file and
1588 // store them in outputContainer
1589
1590 TList * outputContainer = new TList() ;
1591 outputContainer->SetName("QAHistos") ;
1592
1593 //Histograms
745913ae 1594 Int_t nptbins = GetHistogramRanges()->GetHistoPtBins(); Float_t ptmax = GetHistogramRanges()->GetHistoPtMax(); Float_t ptmin = GetHistogramRanges()->GetHistoPtMin();
1595 Int_t nfineptbins = GetHistogramRanges()->GetHistoFinePtBins(); Float_t ptfinemax = GetHistogramRanges()->GetHistoFinePtMax(); Float_t ptfinemin = GetHistogramRanges()->GetHistoFinePtMin();
1596 Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins(); Float_t phimax = GetHistogramRanges()->GetHistoPhiMax(); Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();
1597 Int_t netabins = GetHistogramRanges()->GetHistoEtaBins(); Float_t etamax = GetHistogramRanges()->GetHistoEtaMax(); Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
1598 Int_t nmassbins = GetHistogramRanges()->GetHistoMassBins(); Float_t massmax = GetHistogramRanges()->GetHistoMassMax(); Float_t massmin = GetHistogramRanges()->GetHistoMassMin();
1599 Int_t nasymbins = GetHistogramRanges()->GetHistoAsymmetryBins(); Float_t asymmax = GetHistogramRanges()->GetHistoAsymmetryMax(); Float_t asymmin = GetHistogramRanges()->GetHistoAsymmetryMin();
ec58c056 1600 Int_t nPoverEbins = GetHistogramRanges()->GetHistoPOverEBins(); Float_t eOverPmax = GetHistogramRanges()->GetHistoPOverEMax(); Float_t eOverPmin = GetHistogramRanges()->GetHistoPOverEMin();
745913ae 1601 Int_t ndedxbins = GetHistogramRanges()->GetHistodEdxBins(); Float_t dedxmax = GetHistogramRanges()->GetHistodEdxMax(); Float_t dedxmin = GetHistogramRanges()->GetHistodEdxMin();
1602 Int_t ndRbins = GetHistogramRanges()->GetHistodRBins(); Float_t dRmax = GetHistogramRanges()->GetHistodRMax(); Float_t dRmin = GetHistogramRanges()->GetHistodRMin();
1603 Int_t ntimebins = GetHistogramRanges()->GetHistoTimeBins(); Float_t timemax = GetHistogramRanges()->GetHistoTimeMax(); Float_t timemin = GetHistogramRanges()->GetHistoTimeMin();
1604 Int_t nclbins = GetHistogramRanges()->GetHistoNClustersBins(); Int_t nclmax = GetHistogramRanges()->GetHistoNClustersMax(); Int_t nclmin = GetHistogramRanges()->GetHistoNClustersMin();
1605 Int_t ncebins = GetHistogramRanges()->GetHistoNCellsBins(); Int_t ncemax = GetHistogramRanges()->GetHistoNCellsMax(); Int_t ncemin = GetHistogramRanges()->GetHistoNCellsMin();
1606 Int_t nceclbins = GetHistogramRanges()->GetHistoNClusterCellBins(); Int_t nceclmax = GetHistogramRanges()->GetHistoNClusterCellMax(); Int_t nceclmin = GetHistogramRanges()->GetHistoNClusterCellMin();
1607 Int_t nvdistbins = GetHistogramRanges()->GetHistoVertexDistBins(); Float_t vdistmax = GetHistogramRanges()->GetHistoVertexDistMax(); Float_t vdistmin = GetHistogramRanges()->GetHistoVertexDistMin();
1608 Int_t rbins = GetHistogramRanges()->GetHistoRBins(); Float_t rmax = GetHistogramRanges()->GetHistoRMax(); Float_t rmin = GetHistogramRanges()->GetHistoRMin();
1609 Int_t xbins = GetHistogramRanges()->GetHistoXBins(); Float_t xmax = GetHistogramRanges()->GetHistoXMax(); Float_t xmin = GetHistogramRanges()->GetHistoXMin();
1610 Int_t ybins = GetHistogramRanges()->GetHistoYBins(); Float_t ymax = GetHistogramRanges()->GetHistoYMax(); Float_t ymin = GetHistogramRanges()->GetHistoYMin();
1611 Int_t zbins = GetHistogramRanges()->GetHistoZBins(); Float_t zmax = GetHistogramRanges()->GetHistoZMax(); Float_t zmin = GetHistogramRanges()->GetHistoZMin();
1612 Int_t ssbins = GetHistogramRanges()->GetHistoShowerShapeBins(); Float_t ssmax = GetHistogramRanges()->GetHistoShowerShapeMax(); Float_t ssmin = GetHistogramRanges()->GetHistoShowerShapeMin();
1613 Int_t tdbins = GetHistogramRanges()->GetHistoDiffTimeBins() ; Float_t tdmax = GetHistogramRanges()->GetHistoDiffTimeMax(); Float_t tdmin = GetHistogramRanges()->GetHistoDiffTimeMin();
1614
1615 Int_t nv0sbins = GetHistogramRanges()->GetHistoV0SignalBins(); Int_t nv0smax = GetHistogramRanges()->GetHistoV0SignalMax(); Int_t nv0smin = GetHistogramRanges()->GetHistoV0SignalMin();
1616 Int_t nv0mbins = GetHistogramRanges()->GetHistoV0MultiplicityBins(); Int_t nv0mmax = GetHistogramRanges()->GetHistoV0MultiplicityMax(); Int_t nv0mmin = GetHistogramRanges()->GetHistoV0MultiplicityMin();
1617 Int_t ntrmbins = GetHistogramRanges()->GetHistoTrackMultiplicityBins(); Int_t ntrmmax = GetHistogramRanges()->GetHistoTrackMultiplicityMax(); Int_t ntrmmin = GetHistogramRanges()->GetHistoTrackMultiplicityMin();
649b825d 1618
1619 //EMCAL
1620 fNMaxCols = 48;
1621 fNMaxRows = 24;
1622 fNRCU = 2 ;
1623 //PHOS
f1538a5f 1624 if(fCalorimeter=="PHOS")
1625 {
649b825d 1626 fNMaxCols = 56;
1627 fNMaxRows = 64;
1628 fNRCU = 4 ;
1629 }
1630
1631 fhE = new TH1F ("hE","E reconstructed clusters ", nptbins*5,ptmin,ptmax*5);
1632 fhE->SetXTitle("E (GeV)");
1633 outputContainer->Add(fhE);
1634
f1538a5f 1635 if(fFillAllTH12)
1636 {
649b825d 1637 fhPt = new TH1F ("hPt","p_{T} reconstructed clusters", nptbins,ptmin,ptmax);
1638 fhPt->SetXTitle("p_{T} (GeV/c)");
1639 outputContainer->Add(fhPt);
1640
1641 fhPhi = new TH1F ("hPhi","#phi reconstructed clusters ",nphibins,phimin,phimax);
1642 fhPhi->SetXTitle("#phi (rad)");
1643 outputContainer->Add(fhPhi);
1644
1645 fhEta = new TH1F ("hEta","#eta reconstructed clusters ",netabins,etamin,etamax);
1646 fhEta->SetXTitle("#eta ");
1647 outputContainer->Add(fhEta);
1648 }
1649
f1538a5f 1650 if(fFillAllTH3)
1651 {
649b825d 1652 fhEtaPhiE = new TH3F ("hEtaPhiE","#eta vs #phi vs energy, reconstructed clusters",
1653 netabins,etamin,etamax,nphibins,phimin,phimax,nptbins,ptmin,ptmax);
1654 fhEtaPhiE->SetXTitle("#eta ");
1655 fhEtaPhiE->SetYTitle("#phi (rad)");
1656 fhEtaPhiE->SetZTitle("E (GeV) ");
1657 outputContainer->Add(fhEtaPhiE);
1658 }
1659
1660 fhClusterTimeEnergy = new TH2F ("hClusterTimeEnergy","energy vs TOF, reconstructed clusters",
1661 nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
1662 fhClusterTimeEnergy->SetXTitle("E (GeV) ");
1663 fhClusterTimeEnergy->SetYTitle("TOF (ns)");
1664 outputContainer->Add(fhClusterTimeEnergy);
1665
1666 fhClusterPairDiffTimeE = new TH2F("hClusterPairDiffTimeE","cluster pair time difference vs E, only good clusters",
1667 nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1668 fhClusterPairDiffTimeE->SetXTitle("E_{cluster} (GeV)");
1669 fhClusterPairDiffTimeE->SetYTitle("#Delta t (ns)");
1670 outputContainer->Add(fhClusterPairDiffTimeE);
1671
1672 fhLambda0 = new TH2F ("hLambda0","shower shape, #lambda^{2}_{0} vs E",
1673 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1674 fhLambda0->SetXTitle("E_{cluster}");
1675 fhLambda0->SetYTitle("#lambda^{2}_{0}");
1676 outputContainer->Add(fhLambda0);
1677
1678 fhLambda1 = new TH2F ("hLambda1","shower shape, #lambda^{2}_{1} vs E for bad cluster ",
1679 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1680 fhLambda1->SetXTitle("E_{cluster}");
1681 fhLambda1->SetYTitle("#lambda^{2}_{1}");
1682 outputContainer->Add(fhLambda1);
1683
1684 fhDispersion = new TH2F ("hDispersion","shower shape, Dispersion^{2} vs E for bad cluster ",
1685 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1686 fhDispersion->SetXTitle("E_{cluster}");
1687 fhDispersion->SetYTitle("Dispersion");
1688 outputContainer->Add(fhDispersion);
1689
1690 fhClusterMaxCellCloseCellRatio = new TH2F ("hClusterMaxCellCloseCellRatio","energy vs ratio of max cell / neighbour cell, reconstructed clusters",
1691 nptbins,ptmin,ptmax, 100,0,1.);
1692 fhClusterMaxCellCloseCellRatio->SetXTitle("E_{cluster} (GeV) ");
1693 fhClusterMaxCellCloseCellRatio->SetYTitle("E_{cell i}/E_{cell max}");
1694 outputContainer->Add(fhClusterMaxCellCloseCellRatio);
1695
1696 fhClusterMaxCellCloseCellDiff = new TH2F ("hClusterMaxCellCloseCellDiff","energy vs ratio of max cell / neighbour cell, reconstructed clusters",
1697 nptbins,ptmin,ptmax, 500,0,100.);
1698 fhClusterMaxCellCloseCellDiff->SetXTitle("E_{cluster} (GeV) ");
1699 fhClusterMaxCellCloseCellDiff->SetYTitle("E_{cell max}-E_{cell i} (GeV)");
1700 outputContainer->Add(fhClusterMaxCellCloseCellDiff);
1701
1702 fhClusterMaxCellDiff = new TH2F ("hClusterMaxCellDiff","energy vs difference of cluster energy - max cell energy / cluster energy, good clusters",
1703 nptbins,ptmin,ptmax, 500,0,1.);
1704 fhClusterMaxCellDiff->SetXTitle("E_{cluster} (GeV) ");
1705 fhClusterMaxCellDiff->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
1706 outputContainer->Add(fhClusterMaxCellDiff);
1707
1708 fhClusterMaxCellDiffNoCut = new TH2F ("hClusterMaxCellDiffNoCut","energy vs difference of cluster energy - max cell energy / cluster energy",
1709 nptbins,ptmin,ptmax, 500,0,1.);
1710 fhClusterMaxCellDiffNoCut->SetXTitle("E_{cluster} (GeV) ");
1711 fhClusterMaxCellDiffNoCut->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
1712 outputContainer->Add(fhClusterMaxCellDiffNoCut);
1713
1a72f6c5 1714 fhClusterMaxCellECross = new TH2F ("hClusterMaxCellECross","1 - Energy in cross around max energy cell / max energy cell vs cluster energy, good clusters",
1715 nptbins,ptmin,ptmax, 400,-1,1.);
1716 fhClusterMaxCellECross->SetXTitle("E_{cluster} (GeV) ");
1717 fhClusterMaxCellECross->SetYTitle("1- E_{cross}/E_{cell max}");
1718 outputContainer->Add(fhClusterMaxCellECross);
1719
f1538a5f 1720 fhNCellsPerClusterNoCut = new TH2F ("hNCellsPerClusterNoCut","# cells per cluster vs energy, no bad clusters cut",
1721 nptbins,ptmin,ptmax, nceclbins,nceclmin,nceclmax);
1722 fhNCellsPerClusterNoCut->SetXTitle("E (GeV)");
1723 fhNCellsPerClusterNoCut->SetYTitle("n cells");
1724 outputContainer->Add(fhNCellsPerClusterNoCut);
1725
1726 fhNCellsPerCluster = new TH2F ("hNCellsPerCluster","# cells per cluster vs energy",nptbins,ptmin,ptmax, nceclbins,nceclmin,nceclmax);
1727 fhNCellsPerCluster->SetXTitle("E (GeV)");
1728 fhNCellsPerCluster->SetYTitle("n cells");
1729 outputContainer->Add(fhNCellsPerCluster);
1730
1731 fhNClusters = new TH1F ("hNClusters","# clusters", nclbins,nclmin,nclmax);
1732 fhNClusters->SetXTitle("number of clusters");
1733 outputContainer->Add(fhNClusters);
1734
1a83b960 1735 if(fStudyBadClusters)
1736 {
649b825d 1737 fhBadClusterEnergy = new TH1F ("hBadClusterEnergy","Bad cluster energy", nptbins,ptmin,ptmax);
1738 fhBadClusterEnergy->SetXTitle("E_{cluster} (GeV) ");
1739 outputContainer->Add(fhBadClusterEnergy);
1740
1741 fhBadClusterMaxCellCloseCellRatio = new TH2F ("hBadClusterMaxCellCloseCellRatio","energy vs ratio of max cell / neighbour cell constributing cell, reconstructed bad clusters",
1742 nptbins,ptmin,ptmax, 100,0,1.);
1743 fhBadClusterMaxCellCloseCellRatio->SetXTitle("E_{cluster} (GeV) ");
1744 fhBadClusterMaxCellCloseCellRatio->SetYTitle("ratio");
1745 outputContainer->Add(fhBadClusterMaxCellCloseCellRatio);
1746
1747 fhBadClusterMaxCellCloseCellDiff = new TH2F ("hBadClusterMaxCellCloseCellDiff","energy vs ratio of max cell - neighbour cell constributing cell, reconstructed bad clusters",
1748 nptbins,ptmin,ptmax, 500,0,100);
1749 fhBadClusterMaxCellCloseCellDiff->SetXTitle("E_{cluster} (GeV) ");
1750 fhBadClusterMaxCellCloseCellDiff->SetYTitle("E_{cell max} - E_{cell i} (GeV)");
1751 outputContainer->Add(fhBadClusterMaxCellCloseCellDiff);
1752
1753 fhBadClusterMaxCellDiff = new TH2F ("hBadClusterMaxCellDiff","energy vs difference of cluster energy - max cell energy / cluster energy for bad clusters",
1754 nptbins,ptmin,ptmax, 500,0,1.);
1755 fhBadClusterMaxCellDiff->SetXTitle("E_{cluster} (GeV) ");
1756 fhBadClusterMaxCellDiff->SetYTitle("(E_{cluster} - E_{cell max}) / E_{cluster}");
1757 outputContainer->Add(fhBadClusterMaxCellDiff);
1758
1759 fhBadClusterTimeEnergy = new TH2F ("hBadClusterTimeEnergy","energy vs TOF of reconstructed bad clusters",
1760 nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
1761 fhBadClusterTimeEnergy->SetXTitle("E_{cluster} (GeV) ");
1762 fhBadClusterTimeEnergy->SetYTitle("TOF (ns)");
1763 outputContainer->Add(fhBadClusterTimeEnergy);
1764
1765 fhBadClusterPairDiffTimeE = new TH2F("hBadClusterPairDiffTimeE","cluster pair time difference (bad - good) vs E from bad cluster",nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1766 fhBadClusterPairDiffTimeE->SetXTitle("E_{bad cluster} (GeV)");
1767 fhBadClusterPairDiffTimeE->SetYTitle("#Delta t (ns)");
1768 outputContainer->Add(fhBadClusterPairDiffTimeE);
1769
1a72f6c5 1770 fhBadClusterMaxCellECross = new TH2F ("hBadClusterMaxCellECross","1 - Energy in cross around max energy cell / max energy cell vs cluster energy, bad clusters",
1771 nptbins,ptmin,ptmax, 400,-1,1.);
1772 fhBadClusterMaxCellECross->SetXTitle("E_{cluster} (GeV) ");
1773 fhBadClusterMaxCellECross->SetYTitle("1- E_{cross}/E_{cell max}");
1774 outputContainer->Add(fhBadClusterMaxCellECross);
1775
e6fec6f5 1776 if(fFillAllCellTimeHisto)
1777 {
649b825d 1778 fhBadCellTimeSpreadRespectToCellMax = new TH2F ("hBadCellTimeSpreadRespectToCellMax","t_{cell max}-t_{cell i} from bad cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1779 fhBadCellTimeSpreadRespectToCellMax->SetXTitle("E (GeV)");
1780 fhBadCellTimeSpreadRespectToCellMax->SetYTitle("#Delta t_{cell max - i} (ns)");
1781 outputContainer->Add(fhBadCellTimeSpreadRespectToCellMax);
1782
1783 fhBadClusterMaxCellDiffAverageTime = new TH2F ("hBadClusterMaxCellDiffAverageTime","t_{cell max}-t_{average} from bad cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1784 fhBadClusterMaxCellDiffAverageTime->SetXTitle("E (GeV)");
1785 fhBadClusterMaxCellDiffAverageTime->SetYTitle("#Delta t_{cell max - average} (ns)");
1786 outputContainer->Add(fhBadClusterMaxCellDiffAverageTime);
a82b4462 1787
649b825d 1788 fhBadClusterMaxCellDiffWeightedTime = new TH2F ("hBadClusterMaxCellDiffWeightedTime","t_{cell max}-t_{weighted} from bad cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1789 fhBadClusterMaxCellDiffWeightedTime->SetXTitle("E (GeV)");
1790 fhBadClusterMaxCellDiffWeightedTime->SetYTitle("#Delta t_{cell max - weighted} (ns)");
1791 outputContainer->Add(fhBadClusterMaxCellDiffWeightedTime);
1792
649b825d 1793 }
1794
1795 }
1796
f1538a5f 1797 if(fStudyExotic)
1798 {
765206a5 1799 fhExoL0ECross = new TH2F("hExoL0_ECross",
1800 "#lambda^{2}_{0} vs 1-E_{+}/E_{max} for E > 5 GeV",
e49010f3 1801 400,0,1,ssbins,ssmin,ssmax);
cdacd584 1802 fhExoL0ECross ->SetXTitle("1-E_{+}/E_{cell max}");
1803 fhExoL0ECross ->SetYTitle("#lambda^{2}_{0}");
765206a5 1804 outputContainer->Add(fhExoL0ECross) ;
1805
1806 fhExoL1ECross = new TH2F("hExoL1_ECross",
1807 "#lambda^{2}_{1} vs 1-E_{+}/E_{max} for E > 5 GeV",
e49010f3 1808 400,0,1,ssbins,ssmin,ssmax);
cdacd584 1809 fhExoL1ECross ->SetXTitle("1-E_{+}/E_{cell max}");
1810 fhExoL1ECross ->SetYTitle("#lambda^{2}_{1}");
765206a5 1811 outputContainer->Add(fhExoL1ECross) ;
1812
f1538a5f 1813 for(Int_t ie = 0; ie <fExoNECrossCuts; ie++)
1814 {
1815
1816 fhExoDTime[ie] = new TH2F(Form("hExoDTime_ECross%d",ie),
1817 Form("#Delta time = t_{max}-t_{cells} vs E_{cluster} for exotic, 1-E_{+}/E_{max} < %2.2f",fExoECrossCuts[ie]),
1818 nptbins,ptmin,ptmax,tdbins,tdmin,tdmax);
1819 fhExoDTime[ie] ->SetYTitle("#Delta t (ns)");
1820 fhExoDTime[ie] ->SetXTitle("E (GeV)");
1821 outputContainer->Add(fhExoDTime[ie]) ;
1822
1823 for(Int_t idt = 0; idt < fExoNDTimeCuts; idt++)
1824 {
1825 fhExoNCell[ie][idt] = new TH2F(Form("hExoNCell_ECross%d_DT%d",ie,idt),
1826 Form("N cells per cluster vs E cluster, 1-E_{+}/E_{max} < %2.2f, #Delta t < %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
1827 nptbins,ptmin,ptmax,nceclbins,nceclmin,nceclmax);
1828 fhExoNCell[ie][idt] ->SetYTitle("N cells");
1829 fhExoNCell[ie][idt] ->SetXTitle("E (GeV)");
1830 outputContainer->Add(fhExoNCell[ie][idt]) ;
1831
1832 fhExoL0 [ie][idt] = new TH2F(Form("hExoL0_ECross%d_DT%d",ie,idt),
765206a5 1833 Form("#lambda^{2}_{0} vs E cluster for exotic, 1-E_{+}/E_{max} < %2.2f, #Delta t = %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
f1538a5f 1834 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1835 fhExoL0 [ie][idt] ->SetYTitle("#lambda^{2}_{0}");
1836 fhExoL0 [ie][idt] ->SetXTitle("E (GeV)");
1837 outputContainer->Add(fhExoL0[ie][idt]) ;
765206a5 1838
1839 fhExoL1 [ie][idt] = new TH2F(Form("hExoL1_ECross%d_DT%d",ie,idt),
1840 Form("#lambda^{2}_{1} vs E cluster for exotic, 1-E_{+}/E_{max} < %2.2f, #Delta t = %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
1841 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1842 fhExoL1 [ie][idt] ->SetYTitle("#lambda^{2}_{1}");
1843 fhExoL1 [ie][idt] ->SetXTitle("E (GeV)");
1844 outputContainer->Add(fhExoL1[ie][idt]) ;
f1538a5f 1845
1846 fhExoECross[ie][idt] = new TH2F(Form("hExoECross_ECross%d_DT%d",ie,idt),
1847 Form("E cross for cells vs E cell, 1-E_{+}/E_{max} < %2.2f, #Delta t < %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
1848 nptbins,ptmin,ptmax,400,0,1);
1849 fhExoECross[ie][idt] ->SetYTitle("1-E_{+}/E_{cell max}");
1850 fhExoECross[ie][idt] ->SetXTitle("E_{cell} (GeV)");
1851 outputContainer->Add(fhExoECross[ie][idt]) ;
1852
1853 fhExoTime [ie][idt] = new TH2F(Form("hExoTime_ECross%d_DT%d",ie,idt),
1854 Form("Time of cluster (max cell) vs E cluster for exotic, 1-E_{+}/E_{max} < %2.2f, #Delta t = %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
1855 nptbins,ptmin,ptmax,ntimebins,timemin,timemax);
1856 fhExoTime [ie][idt] ->SetYTitle("time_{max} (ns)");
1857 fhExoTime [ie][idt] ->SetXTitle("E (GeV)");
1858 outputContainer->Add(fhExoTime[ie][idt]) ;
1859
765206a5 1860 fhExoL0NCell[ie][idt] = new TH2F(Form("hExoL0_NCell%d_DT%d",ie,idt),
1861 Form("#lambda^{2}_{0} vs N cells per clusters for E > 5 GeV, 1-E_{+}/E_{max} < %2.2f, #Delta t = %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
1862 nptbins,ptmin,ptmax,ntimebins,timemin,timemax);
1863 fhExoL0NCell[ie][idt] ->SetYTitle("N cells");
1864 fhExoL0NCell[ie][idt] ->SetXTitle("#lambda^{2}_{0}");
1865 outputContainer->Add(fhExoL0NCell[ie][idt]) ;
1866
1867 fhExoL1NCell[ie][idt] = new TH2F(Form("hExoL1_NCell%d_DT%d",ie,idt),
1868 Form("#lambda^{2}_{1} vs N cells per clusters for E > 5 GeV, 1-E_{+}/E_{max} < %2.2f, #Delta t = %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
1869 nptbins,ptmin,ptmax,ntimebins,timemin,timemax);
1870 fhExoL1NCell[ie][idt] ->SetYTitle("N cells");
1871 fhExoL1NCell[ie][idt] ->SetXTitle("#lambda^{2}_{1}");
1872 outputContainer->Add(fhExoL1NCell[ie][idt]) ;
1873
f1538a5f 1874 }
1875 }
1876 }
1877
649b825d 1878 // Cluster size in terms of cells
f1538a5f 1879 if(fStudyClustersAsymmetry)
1880 {
649b825d 1881 fhDeltaIEtaDeltaIPhiE0[0] = new TH2F ("hDeltaIEtaDeltaIPhiE0"," Cluster size in columns vs rows for E < 2 GeV, n cells > 3",
1882 50,0,50,50,0,50);
1883 fhDeltaIEtaDeltaIPhiE0[0]->SetXTitle("#Delta Column");
1884 fhDeltaIEtaDeltaIPhiE0[0]->SetYTitle("#Delta Row");
1885 outputContainer->Add(fhDeltaIEtaDeltaIPhiE0[0]);
1886
1887 fhDeltaIEtaDeltaIPhiE2[0] = new TH2F ("hDeltaIEtaDeltaIPhiE2"," Cluster size in columns vs rows for 2 <E < 6 GeV, n cells > 3",
1888 50,0,50,50,0,50);
1889 fhDeltaIEtaDeltaIPhiE2[0]->SetXTitle("#Delta Column");
1890 fhDeltaIEtaDeltaIPhiE2[0]->SetYTitle("#Delta Row");
1891 outputContainer->Add(fhDeltaIEtaDeltaIPhiE2[0]);
1892
1893 fhDeltaIEtaDeltaIPhiE6[0] = new TH2F ("hDeltaIEtaDeltaIPhiE6"," Cluster size in columns vs rows for E > 6 GeV, n cells > 3",
1894 50,0,50,50,0,50);
1895 fhDeltaIEtaDeltaIPhiE6[0]->SetXTitle("#Delta Column");
1896 fhDeltaIEtaDeltaIPhiE6[0]->SetYTitle("#Delta Row");
1897 outputContainer->Add(fhDeltaIEtaDeltaIPhiE6[0]);
1898
1899 fhDeltaIA[0] = new TH2F ("hDeltaIA"," Cluster *asymmetry* in cell units vs E",
1900 nptbins,ptmin,ptmax,21,-1.05,1.05);
1901 fhDeltaIA[0]->SetXTitle("E_{cluster}");
1902 fhDeltaIA[0]->SetYTitle("A_{cell in cluster}");
1903 outputContainer->Add(fhDeltaIA[0]);
1904
1905 fhDeltaIAL0[0] = new TH2F ("hDeltaIAL0"," Cluster *asymmetry* in cell units vs #lambda^{2}_{0}",
1906 ssbins,ssmin,ssmax,21,-1.05,1.05);
1907 fhDeltaIAL0[0]->SetXTitle("#lambda^{2}_{0}");
1908 fhDeltaIAL0[0]->SetYTitle("A_{cell in cluster}");
1909 outputContainer->Add(fhDeltaIAL0[0]);
1910
1911 fhDeltaIAL1[0] = new TH2F ("hDeltaIAL1"," Cluster *asymmetry* in cell units vs #lambda^{2}_{1}",
1912 ssbins,ssmin,ssmax,21,-1.05,1.05);
1913 fhDeltaIAL1[0]->SetXTitle("#lambda^{2}_{1}");
1914 fhDeltaIAL1[0]->SetYTitle("A_{cell in cluster}");
1915 outputContainer->Add(fhDeltaIAL1[0]);
1916
1917 fhDeltaIANCells[0] = new TH2F ("hDeltaIANCells"," Cluster *asymmetry* in cell units vs N cells in cluster",
1918 nceclbins,nceclmin,nceclmax,21,-1.05,1.05);
1919 fhDeltaIANCells[0]->SetXTitle("N_{cell in cluster}");
1920 fhDeltaIANCells[0]->SetYTitle("A_{cell in cluster}");
1921 outputContainer->Add(fhDeltaIANCells[0]);
1922
1923
1924 fhDeltaIEtaDeltaIPhiE0[1] = new TH2F ("hDeltaIEtaDeltaIPhiE0Charged"," Cluster size in columns vs rows for E < 2 GeV, n cells > 3, matched with track",
1925 50,0,50,50,0,50);
1926 fhDeltaIEtaDeltaIPhiE0[1]->SetXTitle("#Delta Column");
1927 fhDeltaIEtaDeltaIPhiE0[1]->SetYTitle("#Delta Row");
1928 outputContainer->Add(fhDeltaIEtaDeltaIPhiE0[1]);
1929
1930 fhDeltaIEtaDeltaIPhiE2[1] = new TH2F ("hDeltaIEtaDeltaIPhiE2Charged"," Cluster size in columns vs rows for 2 <E < 6 GeV, n cells > 3, matched with track",
1931 50,0,50,50,0,50);
1932 fhDeltaIEtaDeltaIPhiE2[1]->SetXTitle("#Delta Column");
1933 fhDeltaIEtaDeltaIPhiE2[1]->SetYTitle("#Delta Row");
1934 outputContainer->Add(fhDeltaIEtaDeltaIPhiE2[1]);
1935
1936 fhDeltaIEtaDeltaIPhiE6[1] = new TH2F ("hDeltaIEtaDeltaIPhiE6Charged"," Cluster size in columns vs rows for E > 6 GeV, n cells > 3, matched with track",
1937 50,0,50,50,0,50);
1938 fhDeltaIEtaDeltaIPhiE6[1]->SetXTitle("#Delta Column");
1939 fhDeltaIEtaDeltaIPhiE6[1]->SetYTitle("#Delta Row");
1940 outputContainer->Add(fhDeltaIEtaDeltaIPhiE6[1]);
1941
1942 fhDeltaIA[1] = new TH2F ("hDeltaIACharged"," Cluster *asymmetry* in cell units vs E, matched with track",
1943 nptbins,ptmin,ptmax,21,-1.05,1.05);
1944 fhDeltaIA[1]->SetXTitle("E_{cluster}");
1945 fhDeltaIA[1]->SetYTitle("A_{cell in cluster}");
1946 outputContainer->Add(fhDeltaIA[1]);
1947
1948 fhDeltaIAL0[1] = new TH2F ("hDeltaIAL0Charged"," Cluster *asymmetry* in cell units vs #lambda^{2}_{0}, matched with track",
1949 ssbins,ssmin,ssmax,21,-1.05,1.05);
1950 fhDeltaIAL0[1]->SetXTitle("#lambda^{2}_{0}");
1951 fhDeltaIAL0[1]->SetYTitle("A_{cell in cluster}");
1952 outputContainer->Add(fhDeltaIAL0[1]);
1953
1954 fhDeltaIAL1[1] = new TH2F ("hDeltaIAL1Charged"," Cluster *asymmetry* in cell units vs #lambda^{2}_{1}, matched with track",
1955 ssbins,ssmin,ssmax,21,-1.05,1.05);
1956 fhDeltaIAL1[1]->SetXTitle("#lambda^{2}_{1}");
1957 fhDeltaIAL1[1]->SetYTitle("A_{cell in cluster}");
1958 outputContainer->Add(fhDeltaIAL1[1]);
1959
1960 fhDeltaIANCells[1] = new TH2F ("hDeltaIANCellsCharged"," Cluster *asymmetry* in cell units vs N cells in cluster, matched with track",
1961 nceclbins,nceclmin,nceclmax,21,-1.05,1.05);
1962 fhDeltaIANCells[1]->SetXTitle("N_{cell in cluster}");
1963 fhDeltaIANCells[1]->SetYTitle("A_{cell in cluster}");
1964 outputContainer->Add(fhDeltaIANCells[1]);
1965
1966 if(IsDataMC()){
1967 TString particle[]={"Photon","Electron","Conversion","Hadron"};
1968 for (Int_t iPart = 0; iPart < 4; iPart++) {
1969
1970 fhDeltaIAMC[iPart] = new TH2F (Form("hDeltaIA_MC%s",particle[iPart].Data()),Form(" Cluster *asymmetry* in cell units vs E, from %s",particle[iPart].Data()),
1971 nptbins,ptmin,ptmax,21,-1.05,1.05);
1972 fhDeltaIAMC[iPart]->SetXTitle("E_{cluster}");
1973 fhDeltaIAMC[iPart]->SetYTitle("A_{cell in cluster}");
1974 outputContainer->Add(fhDeltaIAMC[iPart]);
1975 }
1976 }
f1538a5f 1977
1a83b960 1978 if(fStudyBadClusters)
1979 {
1980 fhBadClusterDeltaIEtaDeltaIPhiE0 = new TH2F ("hBadClusterDeltaIEtaDeltaIPhiE0"," Cluster size in columns vs rows for E < 2 GeV, n cells > 3",
1981 50,0,50,50,0,50);
1982 fhBadClusterDeltaIEtaDeltaIPhiE0->SetXTitle("#Delta Column");
1983 fhBadClusterDeltaIEtaDeltaIPhiE0->SetYTitle("#Delta Row");
1984 outputContainer->Add(fhBadClusterDeltaIEtaDeltaIPhiE0);
1985
1986 fhBadClusterDeltaIEtaDeltaIPhiE2 = new TH2F ("hBadClusterDeltaIEtaDeltaIPhiE2"," Cluster size in columns vs rows for 2 <E < 6 GeV, n cells > 3",
1987 50,0,50,50,0,50);
1988 fhBadClusterDeltaIEtaDeltaIPhiE2->SetXTitle("#Delta Column");
1989 fhBadClusterDeltaIEtaDeltaIPhiE2->SetYTitle("#Delta Row");
1990 outputContainer->Add(fhBadClusterDeltaIEtaDeltaIPhiE2);
1991
1992 fhBadClusterDeltaIEtaDeltaIPhiE6 = new TH2F ("hBadClusterDeltaIEtaDeltaIPhiE6"," Cluster size in columns vs rows for E > 6 GeV, n cells > 3",
1993 50,0,50,50,0,50);
1994 fhBadClusterDeltaIEtaDeltaIPhiE6->SetXTitle("#Delta Column");
1995 fhBadClusterDeltaIEtaDeltaIPhiE6->SetYTitle("#Delta Row");
1996 outputContainer->Add(fhBadClusterDeltaIEtaDeltaIPhiE6);
1997
1998 fhBadClusterDeltaIA = new TH2F ("hBadClusterDeltaIA"," Cluster *asymmetry* in cell units vs E",
1999 nptbins,ptmin,ptmax,21,-1.05,1.05);
2000 fhBadClusterDeltaIA->SetXTitle("E_{cluster}");
2001 fhBadClusterDeltaIA->SetYTitle("A_{cell in cluster}");
2002 outputContainer->Add(fhBadClusterDeltaIA);
2003 }
649b825d 2004 }
2005
f1538a5f 2006 if(fStudyWeight)
2007 {
649b825d 2008 fhECellClusterRatio = new TH2F ("hECellClusterRatio"," cell energy / cluster energy vs cluster energy",
2009 nptbins,ptmin,ptmax, 100,0,1.);
2010 fhECellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
2011 fhECellClusterRatio->SetYTitle("E_{cell i}/E_{cluster}");
2012 outputContainer->Add(fhECellClusterRatio);
2013
2014 fhECellClusterLogRatio = new TH2F ("hECellClusterLogRatio"," Log(cell energy / cluster energy) vs cluster energy",
1a72f6c5 2015 nptbins,ptmin,ptmax, 100,-10,0);
649b825d 2016 fhECellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
1a72f6c5 2017 fhECellClusterLogRatio->SetYTitle("Log(E_{cell i}/E_{cluster})");
649b825d 2018 outputContainer->Add(fhECellClusterLogRatio);
2019
2020 fhEMaxCellClusterRatio = new TH2F ("hEMaxCellClusterRatio"," max cell energy / cluster energy vs cluster energy",
2021 nptbins,ptmin,ptmax, 100,0,1.);
2022 fhEMaxCellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
2023 fhEMaxCellClusterRatio->SetYTitle("E_{max cell}/E_{cluster}");
2024 outputContainer->Add(fhEMaxCellClusterRatio);
2025
2026 fhEMaxCellClusterLogRatio = new TH2F ("hEMaxCellClusterLogRatio"," Log(max cell energy / cluster energy) vs cluster energy",
1a72f6c5 2027 nptbins,ptmin,ptmax, 100,-10,0);
649b825d 2028 fhEMaxCellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
1a72f6c5 2029 fhEMaxCellClusterLogRatio->SetYTitle("Log (E_{max cell}/E_{cluster})");
649b825d 2030 outputContainer->Add(fhEMaxCellClusterLogRatio);
2031
1a72f6c5 2032 for(Int_t iw = 0; iw < 14; iw++){
2033 fhLambda0ForW0[iw] = new TH2F (Form("hLambda0ForW0%d",iw),Form("shower shape, #lambda^{2}_{0} vs E, w0 = %1.1f",1+0.5*iw),
649b825d 2034 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2035 fhLambda0ForW0[iw]->SetXTitle("E_{cluster}");
2036 fhLambda0ForW0[iw]->SetYTitle("#lambda^{2}_{0}");
2037 outputContainer->Add(fhLambda0ForW0[iw]);
2038
1a72f6c5 2039// fhLambda1ForW0[iw] = new TH2F (Form("hLambda1ForW0%d",iw),Form("shower shape, #lambda^{2}_{1} vs E, w0 = %1.1f",1+0.5*iw),
2040// nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2041// fhLambda1ForW0[iw]->SetXTitle("E_{cluster}");
2042// fhLambda1ForW0[iw]->SetYTitle("#lambda^{2}_{1}");
2043// outputContainer->Add(fhLambda1ForW0[iw]);
649b825d 2044
2045 if(IsDataMC()){
2046 TString mcnames[] = {"Photon", "Electron","Conversion","Pi0","Hadron"};
2047 for(Int_t imc = 0; imc < 5; imc++){
2048 fhLambda0ForW0MC[iw][imc] = new TH2F (Form("hLambda0ForW0%d_MC%s",iw,mcnames[imc].Data()),
1a72f6c5 2049 Form("shower shape, #lambda^{2}_{0} vs E, w0 = %1.1f, for MC %s",1+0.5*iw,mcnames[imc].Data()),
649b825d 2050 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2051 fhLambda0ForW0MC[iw][imc]->SetXTitle("E_{cluster}");
2052 fhLambda0ForW0MC[iw][imc]->SetYTitle("#lambda^{2}_{0}");
2053 outputContainer->Add(fhLambda0ForW0MC[iw][imc]);
2054
1a72f6c5 2055// fhLambda1ForW0MC[iw][imc] = new TH2F (Form("hLambda1ForW0%d_MC%s",iw,mcnames[imc].Data()),
2056// Form("shower shape, #lambda^{2}_{1} vs E, w0 = %1.1f, for MC %s",1+0.5*iw,mcnames[imc].Data()),
2057// nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2058// fhLambda1ForW0MC[iw][imc]->SetXTitle("E_{cluster}");
2059// fhLambda1ForW0MC[iw][imc]->SetYTitle("#lambda^{2}_{1}");
2060// outputContainer->Add(fhLambda1ForW0MC[iw][imc]);
649b825d 2061 }
2062 }
f1538a5f 2063 }
649b825d 2064 }
2065
2066 //Track Matching
f1538a5f 2067 if(fFillAllTMHisto)
2068 {
2069 if(fFillAllTH12)
2070 {
649b825d 2071 fhECharged = new TH1F ("hECharged","E reconstructed clusters, matched with track", nptbins,ptmin,ptmax);
2072 fhECharged->SetXTitle("E (GeV)");
2073 outputContainer->Add(fhECharged);
2074
2075 fhPtCharged = new TH1F ("hPtCharged","p_{T} reconstructed clusters, matched with track", nptbins,ptmin,ptmax);
2076 fhPtCharged->SetXTitle("p_{T} (GeV/c)");
2077 outputContainer->Add(fhPtCharged);
2078
2079 fhPhiCharged = new TH1F ("hPhiCharged","#phi reconstructed clusters, matched with track",nphibins,phimin,phimax);
2080 fhPhiCharged->SetXTitle("#phi (rad)");
2081 outputContainer->Add(fhPhiCharged);
55c05f8c 2082
2083 fhEtaCharged = new TH1F ("hEtaCharged","#eta reconstructed clusters, matched with track",netabins,etamin,etamax);
2084 fhEtaCharged->SetXTitle("#eta ");
2085 outputContainer->Add(fhEtaCharged);
2086 }
3748ffb5 2087 if(fFillAllTH3){
2088 fhEtaPhiECharged = new TH3F ("hEtaPhiECharged","#eta vs #phi, reconstructed clusters, matched with track",
521636d2 2089 netabins,etamin,etamax,nphibins,phimin,phimax,nptbins,ptmin,ptmax);
3748ffb5 2090 fhEtaPhiECharged->SetXTitle("#eta ");
2091 fhEtaPhiECharged->SetYTitle("#phi ");
2092 fhEtaPhiECharged->SetZTitle("E (GeV) ");
2093 outputContainer->Add(fhEtaPhiECharged);
2094 }
55c05f8c 2095
ec58c056 2096 fh1EOverP = new TH2F("h1EOverP","TRACK matches E/p",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
d55bb5e1 2097 fh1EOverP->SetYTitle("E/p");
2098 fh1EOverP->SetXTitle("p_{T} (GeV/c)");
2099 outputContainer->Add(fh1EOverP);
3bfc4732 2100
a87e069d 2101 fh2dR = new TH2F("h2dR","TRACK matches dR",nptbins,ptmin,ptmax,ndRbins,dRmin,dRmax);
2102 fh2dR->SetXTitle("#Delta R (rad)");
2103 fh2dR->SetXTitle("E cluster (GeV)");
2104 outputContainer->Add(fh2dR) ;
3bfc4732 2105
2106 fh2MatchdEdx = new TH2F("h2MatchdEdx","dE/dx vs. p for all matches",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
2107 fh2MatchdEdx->SetXTitle("p (GeV/c)");
2108 fh2MatchdEdx->SetYTitle("<dE/dx>");
2109 outputContainer->Add(fh2MatchdEdx);
2110
2111 fh2EledEdx = new TH2F("h2EledEdx","dE/dx vs. p for electrons",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
2112 fh2EledEdx->SetXTitle("p (GeV/c)");
2113 fh2EledEdx->SetYTitle("<dE/dx>");
2114 outputContainer->Add(fh2EledEdx) ;
2115
ec58c056 2116 fh1EOverPR02 = new TH2F("h1EOverPR02","TRACK matches E/p, all",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
d55bb5e1 2117 fh1EOverPR02->SetYTitle("E/p");
2118 fh1EOverPR02->SetXTitle("p_{T} (GeV/c)");
2119 outputContainer->Add(fh1EOverPR02);
55c05f8c 2120 }
2121
f1538a5f 2122 if(fFillAllPi0Histo)
2123 {
9e9f04cb 2124 fhIM = new TH2F ("hIM","Cluster pairs Invariant mass vs reconstructed pair energy, ncell > 1",nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
55c05f8c 2125 fhIM->SetXTitle("p_{T, cluster pairs} (GeV) ");
2126 fhIM->SetYTitle("M_{cluster pairs} (GeV/c^{2})");
2127 outputContainer->Add(fhIM);
a6f26052 2128
55c05f8c 2129 fhAsym = new TH2F ("hAssym","Cluster pairs Asymmetry vs reconstructed pair energy",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax);
2130 fhAsym->SetXTitle("p_{T, cluster pairs} (GeV) ");
2131 fhAsym->SetYTitle("Asymmetry");
2132 outputContainer->Add(fhAsym);
a6f26052 2133 }
649b825d 2134
c8fe2783 2135
f1538a5f 2136 if(fFillAllPosHisto2)
2137 {
2138 if(fFillAllTH3)
2139 {
55c05f8c 2140 fhXYZ = new TH3F ("hXYZ","Cluster: x vs y vs z",xbins,xmin,xmax,ybins,ymin,ymax,zbins,zmin,zmax);
2141 fhXYZ->SetXTitle("x (cm)");
2142 fhXYZ->SetYTitle("y (cm)");
2143 fhXYZ->SetZTitle("z (cm) ");
2144 outputContainer->Add(fhXYZ);
2145 }
2146
35c71d5c 2147 fhXNCells = new TH2F ("hXNCells","Cluster X position vs N Cells per Cluster",xbins,xmin,xmax,nceclbins,nceclmin,nceclmax);
55c05f8c 2148 fhXNCells->SetXTitle("x (cm)");
2149 fhXNCells->SetYTitle("N cells per cluster");
2150 outputContainer->Add(fhXNCells);
2151
35c71d5c 2152 fhZNCells = new TH2F ("hZNCells","Cluster Z position vs N Cells per Cluster",zbins,zmin,zmax,nceclbins,nceclmin,nceclmax);
55c05f8c 2153 fhZNCells->SetXTitle("z (cm)");
2154 fhZNCells->SetYTitle("N cells per cluster");
2155 outputContainer->Add(fhZNCells);
2156
2157 fhXE = new TH2F ("hXE","Cluster X position vs cluster energy",xbins,xmin,xmax,nptbins,ptmin,ptmax);
2158 fhXE->SetXTitle("x (cm)");
2159 fhXE->SetYTitle("E (GeV)");
2160 outputContainer->Add(fhXE);
2161
2162 fhZE = new TH2F ("hZE","Cluster Z position vs cluster energy",zbins,zmin,zmax,nptbins,ptmin,ptmax);
2163 fhZE->SetXTitle("z (cm)");
2164 fhZE->SetYTitle("E (GeV)");
2165 outputContainer->Add(fhZE);
2166
35c71d5c 2167 fhRNCells = new TH2F ("hRNCells","Cluster R position vs N Cells per Cluster",rbins,rmin,rmax,nceclbins,nceclmin,nceclmax);
55c05f8c 2168 fhRNCells->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
2169 fhRNCells->SetYTitle("N cells per cluster");
2170 outputContainer->Add(fhRNCells);
2171
2172
35c71d5c 2173 fhYNCells = new TH2F ("hYNCells","Cluster Y position vs N Cells per Cluster",ybins,ymin,ymax,nceclbins,nceclmin,nceclmax);
55c05f8c 2174 fhYNCells->SetXTitle("y (cm)");
2175 fhYNCells->SetYTitle("N cells per cluster");
2176 outputContainer->Add(fhYNCells);
2177
2178 fhRE = new TH2F ("hRE","Cluster R position vs cluster energy",rbins,rmin,rmax,nptbins,ptmin,ptmax);
2179 fhRE->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
2180 fhRE->SetYTitle("E (GeV)");
2181 outputContainer->Add(fhRE);
2182
2183 fhYE = new TH2F ("hYE","Cluster Y position vs cluster energy",ybins,ymin,ymax,nptbins,ptmin,ptmax);
2184 fhYE->SetXTitle("y (cm)");
2185 fhYE->SetYTitle("E (GeV)");
2186 outputContainer->Add(fhYE);
2187 }
f1538a5f 2188
2189 if(fFillAllPosHisto)
2190 {
a6f26052 2191 fhRCellE = new TH2F ("hRCellE","Cell R position vs cell energy",rbins,rmin,rmax,nptbins,ptmin,ptmax);
2192 fhRCellE->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
2193 fhRCellE->SetYTitle("E (GeV)");
2194 outputContainer->Add(fhRCellE);
2195
2196 fhXCellE = new TH2F ("hXCellE","Cell X position vs cell energy",xbins,xmin,xmax,nptbins,ptmin,ptmax);
2197 fhXCellE->SetXTitle("x (cm)");
2198 fhXCellE->SetYTitle("E (GeV)");
2199 outputContainer->Add(fhXCellE);
2200
2201 fhYCellE = new TH2F ("hYCellE","Cell Y position vs cell energy",ybins,ymin,ymax,nptbins,ptmin,ptmax);
2202 fhYCellE->SetXTitle("y (cm)");
2203 fhYCellE->SetYTitle("E (GeV)");
2204 outputContainer->Add(fhYCellE);
2205
2206 fhZCellE = new TH2F ("hZCellE","Cell Z position vs cell energy",zbins,zmin,zmax,nptbins,ptmin,ptmax);
2207 fhZCellE->SetXTitle("z (cm)");
2208 fhZCellE->SetYTitle("E (GeV)");
2209 outputContainer->Add(fhZCellE);
2210
2211 fhXYZCell = new TH3F ("hXYZCell","Cell : x vs y vs z",xbins,xmin,xmax,ybins,ymin,ymax,zbins,zmin,zmax);
2212 fhXYZCell->SetXTitle("x (cm)");
2213 fhXYZCell->SetYTitle("y (cm)");
2214 fhXYZCell->SetZTitle("z (cm)");
2215 outputContainer->Add(fhXYZCell);
2216
2217
2218 Float_t dx = TMath::Abs(xmin)+TMath::Abs(xmax);
2219 Float_t dy = TMath::Abs(ymin)+TMath::Abs(ymax);
2220 Float_t dz = TMath::Abs(zmin)+TMath::Abs(zmax);
2221 Float_t dr = TMath::Abs(rmin)+TMath::Abs(rmax);
2222
35c71d5c 2223 fhDeltaCellClusterRNCells = new TH2F ("hDeltaCellClusterRNCells","Cluster-Cell R position vs N Cells per Cluster",rbins*2,-dr,dr,nceclbins,nceclmin,nceclmax);
a6f26052 2224 fhDeltaCellClusterRNCells->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
2225 fhDeltaCellClusterRNCells->SetYTitle("N cells per cluster");
2226 outputContainer->Add(fhDeltaCellClusterRNCells);
2227
35c71d5c 2228 fhDeltaCellClusterXNCells = new TH2F ("hDeltaCellClusterXNCells","Cluster-Cell X position vs N Cells per Cluster",xbins*2,-dx,dx,nceclbins,nceclmin,nceclmax);
a6f26052 2229 fhDeltaCellClusterXNCells->SetXTitle("x (cm)");
2230 fhDeltaCellClusterXNCells->SetYTitle("N cells per cluster");
2231 outputContainer->Add(fhDeltaCellClusterXNCells);
2232
35c71d5c 2233 fhDeltaCellClusterYNCells = new TH2F ("hDeltaCellClusterYNCells","Cluster-Cell Y position vs N Cells per Cluster",ybins*2,-dy,dy,nceclbins,nceclmin,nceclmax);
a6f26052 2234 fhDeltaCellClusterYNCells->SetXTitle("y (cm)");
2235 fhDeltaCellClusterYNCells->SetYTitle("N cells per cluster");
2236 outputContainer->Add(fhDeltaCellClusterYNCells);
2237
35c71d5c 2238 fhDeltaCellClusterZNCells = new TH2F ("hDeltaCellClusterZNCells","Cluster-Cell Z position vs N Cells per Cluster",zbins*2,-dz,dz,nceclbins,nceclmin,nceclmax);
a6f26052 2239 fhDeltaCellClusterZNCells->SetXTitle("z (cm)");
2240 fhDeltaCellClusterZNCells->SetYTitle("N cells per cluster");
2241 outputContainer->Add(fhDeltaCellClusterZNCells);
2242
2243 fhDeltaCellClusterRE = new TH2F ("hDeltaCellClusterRE","Cluster-Cell R position vs cluster energy",rbins*2,-dr,dr,nptbins,ptmin,ptmax);
2244 fhDeltaCellClusterRE->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
2245 fhDeltaCellClusterRE->SetYTitle("E (GeV)");
2246 outputContainer->Add(fhDeltaCellClusterRE);
2247
2248 fhDeltaCellClusterXE = new TH2F ("hDeltaCellClusterXE","Cluster-Cell X position vs cluster energy",xbins*2,-dx,dx,nptbins,ptmin,ptmax);
2249 fhDeltaCellClusterXE->SetXTitle("x (cm)");
2250 fhDeltaCellClusterXE->SetYTitle("E (GeV)");
2251 outputContainer->Add(fhDeltaCellClusterXE);
2252
2253 fhDeltaCellClusterYE = new TH2F ("hDeltaCellClusterYE","Cluster-Cell Y position vs cluster energy",ybins*2,-dy,dy,nptbins,ptmin,ptmax);
2254 fhDeltaCellClusterYE->SetXTitle("y (cm)");
2255 fhDeltaCellClusterYE->SetYTitle("E (GeV)");
2256 outputContainer->Add(fhDeltaCellClusterYE);
2257
2258 fhDeltaCellClusterZE = new TH2F ("hDeltaCellClusterZE","Cluster-Cell Z position vs cluster energy",zbins*2,-dz,dz,nptbins,ptmin,ptmax);
2259 fhDeltaCellClusterZE->SetXTitle("z (cm)");
2260 fhDeltaCellClusterZE->SetYTitle("E (GeV)");
2261 outputContainer->Add(fhDeltaCellClusterZE);
2262
2263 fhEtaPhiAmp = new TH3F ("hEtaPhiAmp","Cell #eta vs cell #phi vs cell energy",netabins,etamin,etamax,nphibins,phimin,phimax,nptbins,ptmin,ptmax);
2264 fhEtaPhiAmp->SetXTitle("#eta ");
2265 fhEtaPhiAmp->SetYTitle("#phi (rad)");
2266 fhEtaPhiAmp->SetZTitle("E (GeV) ");
2267 outputContainer->Add(fhEtaPhiAmp);
2268
2302a644 2269 }
c8fe2783 2270
a6f26052 2271 //Calo cells
e3300002 2272 fhNCells = new TH1F ("hNCells","# cells", ncebins,ncemin+0.5,ncemax);
2302a644 2273 fhNCells->SetXTitle("n cells");
2274 outputContainer->Add(fhNCells);
2275
2276 fhAmplitude = new TH1F ("hAmplitude","Cell Energy", nptbins*2,ptmin,ptmax);
2277 fhAmplitude->SetXTitle("Cell Energy (GeV)");
2278 outputContainer->Add(fhAmplitude);
2279
35c71d5c 2280 fhAmpId = new TH2F ("hAmpId","Cell Energy", nfineptbins,ptfinemin,ptfinemax,fNMaxRows*fNMaxCols*fNModules,0,fNMaxRows*fNMaxCols*fNModules);
2302a644 2281 fhAmpId->SetXTitle("Cell Energy (GeV)");
2282 outputContainer->Add(fhAmpId);
2283
e6fec6f5 2284 if(fFillAllCellTimeHisto)
2285 {
649b825d 2286 fhCellTimeSpreadRespectToCellMax = new TH2F ("hCellTimeSpreadRespectToCellMax","t_{cell max}-t_{cell i} per cluster", nptbins,ptmin,ptmax,tdbins,tdmin,tdmax);
924e319f 2287 fhCellTimeSpreadRespectToCellMax->SetXTitle("E (GeV)");
2288 fhCellTimeSpreadRespectToCellMax->SetYTitle("#Delta t_{cell max-i} (ns)");
2302a644 2289 outputContainer->Add(fhCellTimeSpreadRespectToCellMax);
2290
649b825d 2291 fhClusterMaxCellDiffAverageTime = new TH2F ("hClusterMaxCellDiffAverageTime","t_{cell max}-t_{average} per cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
9e9f04cb 2292 fhClusterMaxCellDiffAverageTime->SetXTitle("E (GeV)");
2293 fhClusterMaxCellDiffAverageTime->SetYTitle("#Delta t_{cell max - average} (ns)");
2294 outputContainer->Add(fhClusterMaxCellDiffAverageTime);
a82b4462 2295
649b825d 2296 fhClusterMaxCellDiffWeightedTime = new TH2F ("hClusterMaxCellDiffWeightedTime","t_{cell max}-t_{weighted} per cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
2297 fhClusterMaxCellDiffWeightedTime->SetXTitle("E (GeV)");
2298 fhClusterMaxCellDiffWeightedTime->SetYTitle("#Delta t_{cell max - weighted} (ns)");
2299 outputContainer->Add(fhClusterMaxCellDiffWeightedTime);
35c71d5c 2300
924e319f 2301 fhCellIdCellLargeTimeSpread= new TH1F ("hCellIdCellLargeTimeSpread","Cells with time 100 ns larger than cell max in cluster ",
35c71d5c 2302 fNMaxCols*fNMaxRows*fNModules,0,fNMaxCols*fNMaxRows*fNModules);
2302a644 2303 fhCellIdCellLargeTimeSpread->SetXTitle("Absolute Cell Id");
2304 outputContainer->Add(fhCellIdCellLargeTimeSpread);
521636d2 2305
2302a644 2306 fhTime = new TH1F ("hTime","Cell Time",ntimebins,timemin,timemax);
2307 fhTime->SetXTitle("Cell Time (ns)");
2308 outputContainer->Add(fhTime);
2309
1a72f6c5 2310 fhTimeVz = new TH2F ("hTimeVz","Cell Time vs vertex, amplitude > 0.5 GeV",100, 0, 50,ntimebins,timemin,timemax);
2311 fhTimeVz->SetXTitle("|v_{z}| (cm)");
2312 fhTimeVz->SetYTitle("Cell Time (ns)");
2313 outputContainer->Add(fhTimeVz);
2314
35c71d5c 2315 fhTimeId = new TH2F ("hTimeId","Cell Time vs Absolute Id",
2316 ntimebins,timemin,timemax,fNMaxRows*fNMaxCols*fNModules,0,fNMaxRows*fNMaxCols*fNModules);
2302a644 2317 fhTimeId->SetXTitle("Cell Time (ns)");
2318 fhTimeId->SetYTitle("Cell Absolute Id");
2319 outputContainer->Add(fhTimeId);
2320
2321 fhTimeAmp = new TH2F ("hTimeAmp","Cell Time vs Cell Energy",nptbins*2,ptmin,ptmax,ntimebins,timemin,timemax);
2322 fhTimeAmp->SetYTitle("Cell Time (ns)");
2323 fhTimeAmp->SetXTitle("Cell Energy (GeV)");
2324 outputContainer->Add(fhTimeAmp);
2325
2302a644 2326 }
06f1b12a 2327
2328 fhCellECross = new TH2F ("hCellECross","1 - Energy in cross around cell / cell energy",
2329 nptbins,ptmin,ptmax, 400,-1,1.);
2330 fhCellECross->SetXTitle("E_{cell} (GeV) ");
2331 fhCellECross->SetYTitle("1- E_{cross}/E_{cell}");
2332 outputContainer->Add(fhCellECross);
2333
1a72f6c5 2334
f1538a5f 2335 if(fCorrelate)
2336 {
798a9b04 2337 //PHOS vs EMCAL
35c71d5c 2338 fhCaloCorrNClusters = new TH2F ("hCaloCorrNClusters","# clusters in EMCAL vs PHOS", nclbins,nclmin,nclmax,nclbins,nclmin,nclmax);
2302a644 2339 fhCaloCorrNClusters->SetXTitle("number of clusters in EMCAL");
2340 fhCaloCorrNClusters->SetYTitle("number of clusters in PHOS");
2341 outputContainer->Add(fhCaloCorrNClusters);
2342
798a9b04 2343 fhCaloCorrEClusters = new TH2F ("hCaloCorrEClusters","summed energy of clusters in EMCAL vs PHOS", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2302a644 2344 fhCaloCorrEClusters->SetXTitle("#Sigma E of clusters in EMCAL (GeV)");
2345 fhCaloCorrEClusters->SetYTitle("#Sigma E of clusters in PHOS (GeV)");
2346 outputContainer->Add(fhCaloCorrEClusters);
2347
35c71d5c 2348 fhCaloCorrNCells = new TH2F ("hCaloCorrNCells","# Cells in EMCAL vs PHOS", ncebins,ncemin,ncemax, ncebins,ncemin,ncemax);
649b825d 2349 fhCaloCorrNCells->SetXTitle("number of Cells in EMCAL");
2350 fhCaloCorrNCells->SetYTitle("number of Cells in PHOS");
2351 outputContainer->Add(fhCaloCorrNCells);
2302a644 2352
649b825d 2353 fhCaloCorrECells = new TH2F ("hCaloCorrECells","summed energy of Cells in EMCAL vs PHOS", nptbins*2,ptmin,ptmax*2,nptbins*2,ptmin,ptmax*2);
2354 fhCaloCorrECells->SetXTitle("#Sigma E of Cells in EMCAL (GeV)");
2355 fhCaloCorrECells->SetYTitle("#Sigma E of Cells in PHOS (GeV)");
2356 outputContainer->Add(fhCaloCorrECells);
2302a644 2357
649b825d 2358 //Calorimeter VS V0 signal
2359 fhCaloV0SCorrNClusters = new TH2F ("hCaloV0SNClusters",Form("# clusters in %s vs V0 signal",fCalorimeter.Data()), nv0sbins,nv0smin,nv0smax,nclbins,nclmin,nclmax);
2360 fhCaloV0SCorrNClusters->SetXTitle("V0 signal");
2361 fhCaloV0SCorrNClusters->SetYTitle(Form("number of clusters in %s",fCalorimeter.Data()));
2362 outputContainer->Add(fhCaloV0SCorrNClusters);
2302a644 2363
649b825d 2364 fhCaloV0SCorrEClusters = new TH2F ("hCaloV0SEClusters",Form("summed energy of clusters in %s vs V0 signal",fCalorimeter.Data()), nv0sbins,nv0smin,nv0smax,nptbins,ptmin,ptmax);
2365 fhCaloV0SCorrEClusters->SetXTitle("V0 signal");
2366 fhCaloV0SCorrEClusters->SetYTitle(Form("#Sigma E of clusters in %s (GeV)",fCalorimeter.Data()));
2367 outputContainer->Add(fhCaloV0SCorrEClusters);
2302a644 2368
649b825d 2369 fhCaloV0SCorrNCells = new TH2F ("hCaloV0SNCells",Form("# Cells in %s vs V0 signal",fCalorimeter.Data()), nv0sbins,nv0smin,nv0smax, ncebins,ncemin,ncemax);
2370 fhCaloV0SCorrNCells->SetXTitle("V0 signal");
2371 fhCaloV0SCorrNCells->SetYTitle(Form("number of Cells in %s",fCalorimeter.Data()));
2372 outputContainer->Add(fhCaloV0SCorrNCells);
3bfc4732 2373
649b825d 2374 fhCaloV0SCorrECells = new TH2F ("hCaloV0SECells",Form("summed energy of Cells in %s vs V0 signal",fCalorimeter.Data()), nv0sbins,nv0smin,nv0smax,nptbins,ptmin,ptmax);
2375 fhCaloV0SCorrECells->SetXTitle("V0 signal");
2376 fhCaloV0SCorrECells->SetYTitle(Form("#Sigma E of Cells in %s (GeV)",fCalorimeter.Data()));
2377 outputContainer->Add(fhCaloV0SCorrECells);
3bfc4732 2378
649b825d 2379 //Calorimeter VS V0 multiplicity
2380 fhCaloV0MCorrNClusters = new TH2F ("hCaloV0MNClusters",Form("# clusters in %s vs V0 signal",fCalorimeter.Data()), nv0mbins,nv0mmin,nv0mmax,nclbins,nclmin,nclmax);
2381 fhCaloV0MCorrNClusters->SetXTitle("V0 signal");
2382 fhCaloV0MCorrNClusters->SetYTitle(Form("number of clusters in %s",fCalorimeter.Data()));
2383 outputContainer->Add(fhCaloV0MCorrNClusters);
3bfc4732 2384
649b825d 2385 fhCaloV0MCorrEClusters = new TH2F ("hCaloV0MEClusters",Form("summed energy of clusters in %s vs V0 signal",fCalorimeter.Data()), nv0mbins,nv0mmin,nv0mmax,nptbins,ptmin,ptmax);
2386 fhCaloV0MCorrEClusters->SetXTitle("V0 signal");
2387 fhCaloV0MCorrEClusters->SetYTitle(Form("#Sigma E of clusters in %s (GeV)",fCalorimeter.Data()));
2388 outputContainer->Add(fhCaloV0MCorrEClusters);
3bfc4732 2389
649b825d 2390 fhCaloV0MCorrNCells = new TH2F ("hCaloV0MNCells",Form("# Cells in %s vs V0 signal",fCalorimeter.Data()), nv0mbins,nv0mmin,nv0mmax, ncebins,ncemin,ncemax);
2391 fhCaloV0MCorrNCells->SetXTitle("V0 signal");
2392 fhCaloV0MCorrNCells->SetYTitle(Form("number of Cells in %s",fCalorimeter.Data()));
2393 outputContainer->Add(fhCaloV0MCorrNCells);
3bfc4732 2394
649b825d 2395 fhCaloV0MCorrECells = new TH2F ("hCaloV0MECells",Form("summed energy of Cells in %s vs V0 signal",fCalorimeter.Data()), nv0mbins,nv0mmin,nv0mmax,nptbins,ptmin,ptmax);
2396 fhCaloV0MCorrECells->SetXTitle("V0 signal");
2397 fhCaloV0MCorrECells->SetYTitle(Form("#Sigma E of Cells in %s (GeV)",fCalorimeter.Data()));
2398 outputContainer->Add(fhCaloV0MCorrECells);
3bfc4732 2399
649b825d 2400 //Calorimeter VS Track multiplicity
2401 fhCaloTrackMCorrNClusters = new TH2F ("hCaloTrackMNClusters",Form("# clusters in %s vs # tracks",fCalorimeter.Data()), ntrmbins,ntrmmin,ntrmmax,nclbins,nclmin,nclmax);
2402 fhCaloTrackMCorrNClusters->SetXTitle("# tracks");
2403 fhCaloTrackMCorrNClusters->SetYTitle(Form("number of clusters in %s",fCalorimeter.Data()));
2404 outputContainer->Add(fhCaloTrackMCorrNClusters);
3bfc4732 2405
649b825d 2406 fhCaloTrackMCorrEClusters = new TH2F ("hCaloTrackMEClusters",Form("summed energy of clusters in %s vs # tracks",fCalorimeter.Data()), ntrmbins,ntrmmin,ntrmmax,nptbins,ptmin,ptmax);
2407 fhCaloTrackMCorrEClusters->SetXTitle("# tracks");
2408 fhCaloTrackMCorrEClusters->SetYTitle(Form("#Sigma E of clusters in %s (GeV)",fCalorimeter.Data()));
2409 outputContainer->Add(fhCaloTrackMCorrEClusters);
3bfc4732 2410
649b825d 2411 fhCaloTrackMCorrNCells = new TH2F ("hCaloTrackMNCells",Form("# Cells in %s vs # tracks",fCalorimeter.Data()), ntrmbins,ntrmmin,ntrmmax, ncebins,ncemin,ncemax);
2412 fhCaloTrackMCorrNCells->SetXTitle("# tracks");
2413 fhCaloTrackMCorrNCells->SetYTitle(Form("number of Cells in %s",fCalorimeter.Data()));
2414 outputContainer->Add(fhCaloTrackMCorrNCells);
3bfc4732 2415
649b825d 2416 fhCaloTrackMCorrECells = new TH2F ("hCaloTrackMECells",Form("summed energy of Cells in %s vs # tracks",fCalorimeter.Data()), ntrmbins,ntrmmin,ntrmmax,nptbins,ptmin,ptmax);
2417 fhCaloTrackMCorrECells->SetXTitle("# tracks");
2418 fhCaloTrackMCorrECells->SetYTitle(Form("#Sigma E of Cells in %s (GeV)",fCalorimeter.Data()));
2419 outputContainer->Add(fhCaloTrackMCorrECells);
3bfc4732 2420
3bfc4732 2421
649b825d 2422 }//correlate calorimeters
c8fe2783 2423
649b825d 2424 //Module histograms
9725fd2a 2425
649b825d 2426 fhEMod = new TH2F ("hE_Mod","Cluster reconstructed Energy in each present Module",nptbins,ptmin,ptmax,fNModules,0,fNModules);
2427 fhEMod->SetXTitle("E (GeV)");
2428 fhEMod->SetYTitle("Module");
2429 outputContainer->Add(fhEMod);
c8fe2783 2430
649b825d 2431 fhAmpMod = new TH2F ("hAmp_Mod","Cell energy in each present Module",nptbins,ptmin,ptmax,fNModules,0,fNModules);
2432 fhAmpMod->SetXTitle("E (GeV)");
2433 fhAmpMod->SetYTitle("Module");
2434 outputContainer->Add(fhAmpMod);
3bfc4732 2435
e6fec6f5 2436 if(fFillAllCellTimeHisto)
2437 {
0fb69ade 2438 fhTimeMod = new TH2F ("hTime_Mod","Cell time in each present Module",ntimebins,timemin,timemax,fNModules,0,fNModules);
2439 fhTimeMod->SetXTitle("t (ns)");
2440 fhTimeMod->SetYTitle("Module");
2441 outputContainer->Add(fhTimeMod);
2442 }
3bfc4732 2443
e3300002 2444 fhNClustersMod = new TH2F ("hNClusters_Mod","# clusters vs Module", nclbins,nclmin+0.5,nclmax,fNModules,0,fNModules);
649b825d 2445 fhNClustersMod->SetXTitle("number of clusters");
2446 fhNClustersMod->SetYTitle("Module");
2447 outputContainer->Add(fhNClustersMod);
2302a644 2448
e3300002 2449 fhNCellsMod = new TH2F ("hNCells_Mod","# cells vs Module", ncebins,ncemin+0.5,ncemax,fNModules,0,fNModules);
649b825d 2450 fhNCellsMod->SetXTitle("n cells");
2451 fhNCellsMod->SetYTitle("Module");
2452 outputContainer->Add(fhNCellsMod);
17708df9 2453
649b825d 2454 Int_t colmaxs = fNMaxCols;
2455 Int_t rowmaxs = fNMaxRows;
e6fec6f5 2456 if(fCalorimeter=="EMCAL")
2457 {
649b825d 2458 colmaxs=2*fNMaxCols;
2459 rowmaxs=Int_t(fNModules/2)*fNMaxRows;
2460 }
e6fec6f5 2461 else
2462 {
649b825d 2463 rowmaxs=fNModules*fNMaxRows;
2302a644 2464 }
3bfc4732 2465
649b825d 2466 fhGridCells = new TH2F ("hGridCells",Form("Entries in grid of cells"),
2467 colmaxs+2,-1.5,colmaxs+0.5, rowmaxs+2,-1.5,rowmaxs+0.5);
2468 fhGridCells->SetYTitle("row (phi direction)");
2469 fhGridCells->SetXTitle("column (eta direction)");
2470 outputContainer->Add(fhGridCells);
3bfc4732 2471
649b825d 2472 fhGridCellsE = new TH2F ("hGridCellsE","Accumulated energy in grid of cells",
2473 colmaxs+2,-1.5,colmaxs+0.5, rowmaxs+2,-1.5,rowmaxs+0.5);
2474 fhGridCellsE->SetYTitle("row (phi direction)");
2475 fhGridCellsE->SetXTitle("column (eta direction)");
2476 outputContainer->Add(fhGridCellsE);
3bfc4732 2477
e6fec6f5 2478 if(fFillAllCellTimeHisto)
2479 {
0fb69ade 2480 fhGridCellsTime = new TH2F ("hGridCellsTime","Accumulated time in grid of cells",
2481 colmaxs+2,-1.5,colmaxs+0.5, rowmaxs+2,-1.5,rowmaxs+0.5);
2482 fhGridCellsTime->SetYTitle("row (phi direction)");
2483 fhGridCellsTime->SetXTitle("column (eta direction)");
2484 outputContainer->Add(fhGridCellsTime);
2485 }
3bfc4732 2486
e6fec6f5 2487 fhNCellsPerClusterMod = new TH2F*[fNModules];
649b825d 2488 fhNCellsPerClusterModNoCut = new TH2F*[fNModules];
e6fec6f5 2489 fhIMMod = new TH2F*[fNModules];
2490 if(fFillAllCellTimeHisto) fhTimeAmpPerRCU = new TH2F*[fNModules*fNRCU];
0fb69ade 2491
e6fec6f5 2492 for(Int_t imod = 0; imod < fNModules; imod++)
2493 {
649b825d 2494 fhNCellsPerClusterMod[imod] = new TH2F (Form("hNCellsPerCluster_Mod%d",imod),
2495 Form("# cells per cluster vs cluster energy in Module %d",imod),
2496 nptbins,ptmin,ptmax, nceclbins,nceclmin,nceclmax);
2497 fhNCellsPerClusterMod[imod]->SetXTitle("E (GeV)");
2498 fhNCellsPerClusterMod[imod]->SetYTitle("n cells");
2499 outputContainer->Add(fhNCellsPerClusterMod[imod]);
3bfc4732 2500
649b825d 2501 fhNCellsPerClusterModNoCut[imod] = new TH2F (Form("hNCellsPerClusterNoCut_Mod%d",imod),
2502 Form("# cells per cluster vs cluster energy in Module %d, no cut",imod),
2503 nptbins,ptmin,ptmax, nceclbins,nceclmin,nceclmax);
2504 fhNCellsPerClusterModNoCut[imod]->SetXTitle("E (GeV)");
2505 fhNCellsPerClusterModNoCut[imod]->SetYTitle("n cells");
2506 outputContainer->Add(fhNCellsPerClusterModNoCut[imod]);
3bfc4732 2507
e6fec6f5 2508 if(fFillAllCellTimeHisto)
2509 {
2510 for(Int_t ircu = 0; ircu < fNRCU; ircu++)
2511 {
649b825d 2512 fhTimeAmpPerRCU[imod*fNRCU+ircu] = new TH2F (Form("hTimeAmp_Mod%d_RCU%d",imod,ircu),
2513 Form("Cell Energy vs Cell Time in Module %d, RCU %d ",imod,ircu),
2514 nptbins,ptmin,ptmax,ntimebins,timemin,timemax);
2515 fhTimeAmpPerRCU[imod*fNRCU+ircu]->SetXTitle("E (GeV)");
2516 fhTimeAmpPerRCU[imod*fNRCU+ircu]->SetYTitle("time (ns)");
2517 outputContainer->Add(fhTimeAmpPerRCU[imod*fNRCU+ircu]);
3bfc4732 2518
649b825d 2519 }
3bfc4732 2520 }
2521
e6fec6f5 2522 if(fFillAllPi0Histo)
2523 {
649b825d 2524 fhIMMod[imod] = new TH2F (Form("hIM_Mod%d",imod),
2525 Form("Cluster pairs Invariant mass vs reconstructed pair energy in Module %d, n cell > 1",imod),
2526 nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
2527 fhIMMod[imod]->SetXTitle("p_{T, cluster pairs} (GeV) ");
2528 fhIMMod[imod]->SetYTitle("M_{cluster pairs} (GeV/c^{2})");
2529 outputContainer->Add(fhIMMod[imod]);
3bfc4732 2530
649b825d 2531 }
2532 }
2533
a82b4462 2534 // Monte Carlo Histograms
649b825d 2535
2536 TString particleName[] = { "Photon", "Pi0", "Eta", "Electron", "NeutralHadron", "ChargedHadron" };
2537
f1538a5f 2538 if(IsDataMC())
2539 {
2540 for(Int_t iPart = 0; iPart < 6; iPart++)
2541 {
2542 for(Int_t iCh = 0; iCh < 2; iCh++)
2543 {
649b825d 2544 fhRecoMCRatioE[iPart][iCh] = new TH2F (Form("hRecoMCRatioE_%s_Match%d",particleName[iPart].Data(),iCh),
3fa37dff 2545 Form("Reconstructed/Generated E, %s, Matched %d",particleName[iPart].Data(),iCh),
649b825d 2546 nptbins, ptmin, ptmax, 200,0,2);
3fa37dff 2547 fhRecoMCRatioE[iPart][iCh]->SetYTitle("E_{reconstructed}/E_{generated}");
2548 fhRecoMCRatioE[iPart][iCh]->SetXTitle("E_{reconstructed} (GeV)");
649b825d 2549 outputContainer->Add(fhRecoMCRatioE[iPart][iCh]);
3bfc4732 2550
3bfc4732 2551
649b825d 2552 fhRecoMCDeltaE[iPart][iCh] = new TH2F (Form("hRecoMCDeltaE_%s_Match%d",particleName[iPart].Data(),iCh),
3fa37dff 2553 Form("Generated - Reconstructed E, %s, Matched %d",particleName[iPart].Data(),iCh),
649b825d 2554 nptbins, ptmin, ptmax, nptbins*2,-ptmax,ptmax);
3fa37dff 2555 fhRecoMCDeltaE[iPart][iCh]->SetYTitle("#Delta E (GeV)");
2556 fhRecoMCDeltaE[iPart][iCh]->SetXTitle("E_{reconstructed} (GeV)");
649b825d 2557 outputContainer->Add(fhRecoMCDeltaE[iPart][iCh]);
3bfc4732 2558
649b825d 2559 fhRecoMCDeltaPhi[iPart][iCh] = new TH2F (Form("hRecoMCDeltaPhi_%s_Match%d",particleName[iPart].Data(),iCh),
3fa37dff 2560 Form("Generated - Reconstructed #phi, %s, Matched %d",particleName[iPart].Data(),iCh),
649b825d 2561 nptbins, ptmin, ptmax, nphibins*2,-phimax,phimax);
3fa37dff 2562 fhRecoMCDeltaPhi[iPart][iCh]->SetYTitle("#Delta #phi (rad)");
2563 fhRecoMCDeltaPhi[iPart][iCh]->SetXTitle("E_{reconstructed} (GeV)");
649b825d 2564 outputContainer->Add(fhRecoMCDeltaPhi[iPart][iCh]);
3bfc4732 2565
649b825d 2566 fhRecoMCDeltaEta[iPart][iCh] = new TH2F (Form("hRecoMCDeltaEta_%s_Match%d",particleName[iPart].Data(),iCh),
3fa37dff 2567 Form("Generated - Reconstructed #eta, %s, Matched %d",particleName[iPart].Data(),iCh),
649b825d 2568 nptbins, ptmin, ptmax,netabins*2,-etamax,etamax);
3fa37dff 2569 fhRecoMCDeltaEta[iPart][iCh]->SetYTitle("#Delta #eta ");
2570 fhRecoMCDeltaEta[iPart][iCh]->SetXTitle("E_{reconstructed} (GeV)");
649b825d 2571 outputContainer->Add(fhRecoMCDeltaEta[iPart][iCh]);
e1e62b89 2572
649b825d 2573 fhRecoMCE[iPart][iCh] = new TH2F (Form("hRecoMCE_%s_Match%d",particleName[iPart].Data(),iCh),
2574 Form("E distribution, reconstructed vs generated, %s, Matched %d",particleName[iPart].Data(),iCh),
2575 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2576 fhRecoMCE[iPart][iCh]->SetXTitle("E_{rec} (GeV)");
2577 fhRecoMCE[iPart][iCh]->SetYTitle("E_{gen} (GeV)");
2578 outputContainer->Add(fhRecoMCE[iPart][iCh]);
3bfc4732 2579
649b825d 2580 fhRecoMCPhi[iPart][iCh] = new TH2F (Form("hRecoMCPhi_%s_Match%d",particleName[iPart].Data(),iCh),
2581 Form("#phi distribution, reconstructed vs generated, %s, Matched %d",particleName[iPart].Data(),iCh),
2582 nphibins,phimin,phimax, nphibins,phimin,phimax);
3fa37dff 2583 fhRecoMCPhi[iPart][iCh]->SetXTitle("#phi_{reconstructed} (rad)");
2584 fhRecoMCPhi[iPart][iCh]->SetYTitle("#phi_{generated} (rad)");
649b825d 2585 outputContainer->Add(fhRecoMCPhi[iPart][iCh]);
e1e62b89 2586
649b825d 2587 fhRecoMCEta[iPart][iCh] = new TH2F (Form("hRecoMCEta_%s_Match%d",particleName[iPart].Data(),iCh),
2588 Form("#eta distribution, reconstructed vs generated, %s, Matched %d",particleName[iPart].Data(),iCh),
2589 netabins,etamin,etamax,netabins,etamin,etamax);
3fa37dff 2590 fhRecoMCEta[iPart][iCh]->SetXTitle("#eta_{reconstructed} ");
2591 fhRecoMCEta[iPart][iCh]->SetYTitle("#eta_{generated} ");
649b825d 2592 outputContainer->Add(fhRecoMCEta[iPart][iCh]);
3bfc4732 2593 }
649b825d 2594 }
c8fe2783 2595
649b825d 2596 //Pure MC
f1538a5f 2597 for(Int_t iPart = 0; iPart < 4; iPart++)
2598 {
649b825d 2599 fhGenMCE[iPart] = new TH1F(Form("hGenMCE_%s",particleName[iPart].Data()) ,
2600 Form("p_{T} of generated %s",particleName[iPart].Data()),
2601 nptbins,ptmin,ptmax);
2602 fhGenMCEtaPhi[iPart] = new TH2F(Form("hGenMCEtaPhi_%s",particleName[iPart].Data()),
2603 Form("Y vs #phi of generated %s",particleName[iPart].Data()),
2604 netabins,etamin,etamax,nphibins,phimin,phimax);
2605
2606 fhGenMCE[iPart] ->SetXTitle("p_{T} (GeV/c)");
2607 fhGenMCEtaPhi[iPart]->SetXTitle("#eta");
2608 fhGenMCEtaPhi[iPart]->SetYTitle("#phi (rad)");
3bfc4732 2609
649b825d 2610 outputContainer->Add(fhGenMCE[iPart]);
2611 outputContainer->Add(fhGenMCEtaPhi[iPart]);
521636d2 2612
521636d2 2613
649b825d 2614 fhGenMCAccE[iPart] = new TH1F(Form("hGenMCAccE_%s",particleName[iPart].Data()) ,
2615 Form("p_{T} of generated %s",particleName[iPart].Data()),
2616 nptbins,ptmin,ptmax);
2617 fhGenMCAccEtaPhi[iPart] = new TH2F(Form("hGenMCAccEtaPhi_%s",particleName[iPart].Data()),
2618 Form("Y vs #phi of generated %s",particleName[iPart].Data()),
2619 netabins,etamin,etamax,nphibins,phimin,phimax);
2620
2621 fhGenMCAccE[iPart] ->SetXTitle("p_{T} (GeV/c)");
2622 fhGenMCAccEtaPhi[iPart]->SetXTitle("#eta");
2623 fhGenMCAccEtaPhi[iPart]->SetYTitle("#phi (rad)");
2624
2625 outputContainer->Add(fhGenMCAccE[iPart]);
2626 outputContainer->Add(fhGenMCAccEtaPhi[iPart]);
2627
2628 }
f5036bcb 2629
649b825d 2630 //Vertex of generated particles
2631
2632 fhEMVxyz = new TH2F ("hEMVxyz","Production vertex of reconstructed ElectroMagnetic particles",nvdistbins,vdistmin,vdistmax,nvdistbins,vdistmin,vdistmax);//,100,0,500);
2633 fhEMVxyz->SetXTitle("v_{x}");
2634 fhEMVxyz->SetYTitle("v_{y}");
2635 //fhEMVxyz->SetZTitle("v_{z}");
2636 outputContainer->Add(fhEMVxyz);
2637
2638 fhHaVxyz = new TH2F ("hHaVxyz","Production vertex of reconstructed hadrons",nvdistbins,vdistmin,vdistmax,nvdistbins,vdistmin,vdistmax);//,100,0,500);
2639 fhHaVxyz->SetXTitle("v_{x}");
2640 fhHaVxyz->SetYTitle("v_{y}");
2641 //fhHaVxyz->SetZTitle("v_{z}");
2642 outputContainer->Add(fhHaVxyz);
2643
2644 fhEMR = new TH2F ("hEMR","Distance to production vertex of reconstructed ElectroMagnetic particles vs E rec",nptbins,ptmin,ptmax,nvdistbins,vdistmin,vdistmax);
2645 fhEMR->SetXTitle("E (GeV)");
2646 fhEMR->SetYTitle("TMath::Sqrt(v_{x}^{2}+v_{y}^{2})");
2647 outputContainer->Add(fhEMR);
2648
2649 fhHaR = new TH2F ("hHaR","Distance to production vertex of reconstructed Hadrons vs E rec",nptbins,ptmin,ptmax,nvdistbins,vdistmin,vdistmax);
2650 fhHaR->SetXTitle("E (GeV)");
2651 fhHaR->SetYTitle("TMath::Sqrt(v_{x}^{2}+v_{y}^{2})");
2652 outputContainer->Add(fhHaR);
2653
2654
2655 //Track Matching
2656
ec58c056 2657 fhMCEle1EOverP = new TH2F("hMCEle1EOverP","TRACK matches E/p, MC electrons",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
d55bb5e1 2658 fhMCEle1EOverP->SetYTitle("E/p");
2659 fhMCEle1EOverP->SetXTitle("p_{T} (GeV/c)");
2660 outputContainer->Add(fhMCEle1EOverP);
649b825d 2661
2662 fhMCEle1dR = new TH1F("hMCEle1dR","TRACK matches dR, MC electrons",ndRbins,dRmin,dRmax);
2663 fhMCEle1dR->SetXTitle("#Delta R (rad)");
2664 outputContainer->Add(fhMCEle1dR) ;
2665
2666 fhMCEle2MatchdEdx = new TH2F("hMCEle2MatchdEdx","dE/dx vs. p for all matches, MC electrons",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
2667 fhMCEle2MatchdEdx->SetXTitle("p (GeV/c)");
2668 fhMCEle2MatchdEdx->SetYTitle("<dE/dx>");
2669 outputContainer->Add(fhMCEle2MatchdEdx);
2670
ec58c056 2671 fhMCChHad1EOverP = new TH2F("hMCChHad1EOverP","TRACK matches E/p, MC charged hadrons",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
d55bb5e1 2672 fhMCChHad1EOverP->SetYTitle("E/p");
2673 fhMCChHad1EOverP->SetXTitle("p_{T} (GeV/c)");
2674 outputContainer->Add(fhMCChHad1EOverP);
649b825d 2675
2676 fhMCChHad1dR = new TH1F("hMCChHad1dR","TRACK matches dR, MC charged hadrons",ndRbins,dRmin,dRmax);
2677 fhMCChHad1dR->SetXTitle("#Delta R (rad)");
2678 outputContainer->Add(fhMCChHad1dR) ;
2679
2680 fhMCChHad2MatchdEdx = new TH2F("hMCChHad2MatchdEdx","dE/dx vs. p for all matches, MC charged hadrons",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
2681 fhMCChHad2MatchdEdx->SetXTitle("p (GeV/c)");
2682 fhMCChHad2MatchdEdx->SetYTitle("<dE/dx>");
2683 outputContainer->Add(fhMCChHad2MatchdEdx);
2684
ec58c056 2685 fhMCNeutral1EOverP = new TH2F("hMCNeutral1EOverP","TRACK matches E/p, MC neutrals",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
d55bb5e1 2686 fhMCNeutral1EOverP->SetYTitle("E/p");
2687 fhMCNeutral1EOverP->SetXTitle("p_{T} (GeV/c)");
2688 outputContainer->Add(fhMCNeutral1EOverP);
649b825d 2689
2690 fhMCNeutral1dR = new TH1F("hMCNeutral1dR","TRACK matches dR, MC neutrals",ndRbins,dRmin,dRmax);
2691 fhMCNeutral1dR->SetXTitle("#Delta R (rad)");
2692 outputContainer->Add(fhMCNeutral1dR) ;
2693
2694 fhMCNeutral2MatchdEdx = new TH2F("hMCNeutral2MatchdEdx","dE/dx vs. p for all matches, MC neutrals",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
2695 fhMCNeutral2MatchdEdx->SetXTitle("p (GeV/c)");
2696 fhMCNeutral2MatchdEdx->SetYTitle("<dE/dx>");
2697 outputContainer->Add(fhMCNeutral2MatchdEdx);
2698
ec58c056 2699 fhMCEle1EOverPR02 = new TH2F("hMCEle1EOverPR02","TRACK matches E/p, MC electrons",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
d55bb5e1 2700 fhMCEle1EOverPR02->SetYTitle("E/p");
2701 fhMCEle1EOverPR02->SetXTitle("p_{T} (GeV/c)");
2702 outputContainer->Add(fhMCEle1EOverPR02);
649b825d 2703
ec58c056 2704 fhMCChHad1EOverPR02 = new TH2F("hMCChHad1EOverPR02","TRACK matches E/p, MC charged hadrons",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
d55bb5e1 2705 fhMCChHad1EOverPR02->SetYTitle("E/p");
2706 fhMCChHad1EOverPR02->SetXTitle("p_{T} (GeV/c)");
2707 outputContainer->Add(fhMCChHad1EOverPR02);
f5036bcb 2708
ec58c056 2709 fhMCNeutral1EOverPR02 = new TH2F("hMCNeutral1EOverPR02","TRACK matches E/p, MC neutrals",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
d55bb5e1 2710 fhMCNeutral1EOverPR02->SetYTitle("E/p");
2711 fhMCNeutral1EOverPR02->SetXTitle("p_{T} (GeV/c)");
2712 outputContainer->Add(fhMCNeutral1EOverPR02);
521636d2 2713 }
f5036bcb 2714
649b825d 2715 // for(Int_t i = 0; i < outputContainer->GetEntries() ; i++)
2716 // printf("i=%d, name= %s\n",i,outputContainer->At(i)->GetName());
521636d2 2717
649b825d 2718 return outputContainer;
902aa95c 2719}
2720
f1538a5f 2721//__________________________________________________________________________________________________
2722Float_t AliAnaCalorimeterQA::GetECross(const Int_t absID, AliVCaloCells* cells, const Float_t dtcut)
1a72f6c5 2723{
2724 // Get energy in cross axis around maximum cell, for EMCAL only
2725
06f1b12a 2726 Int_t icol =-1, irow=-1,iRCU = -1;
2727 Int_t imod = GetModuleNumberCellIndexes(absID, fCalorimeter, icol, irow, iRCU);
e3300002 2728
57d8227a 2729 if(fCalorimeter=="EMCAL")
2730 {
06f1b12a 2731 //Get close cells index, energy and time, not in corners
e3300002 2732
2733 Int_t absID1 = -1;
2734 Int_t absID2 = -1;
2735
2736 if( irow < AliEMCALGeoParams::fgkEMCALRows-1) absID1 = GetCaloUtils()->GetEMCALGeometry()->GetAbsCellIdFromCellIndexes(imod, irow+1, icol);
2737 if( irow > 0 ) absID2 = GetCaloUtils()->GetEMCALGeometry()->GetAbsCellIdFromCellIndexes(imod, irow-1, icol);
2747966a 2738
2739 // In case of cell in eta = 0 border, depending on SM shift the cross cell index
2740 Int_t absID3 = -1;
2741 Int_t absID4 = -1;
2742
2743 if ( icol == AliEMCALGeoParams::fgkEMCALCols - 1 && !(imod%2) )
2744 {
e3300002 2745 absID3 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod+1, irow, 0);
2746 absID4 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod , irow, icol-1);
2747966a 2747 }
2748 else if( icol == 0 && imod%2 )
2749 {
e3300002 2750 absID3 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod , irow, icol+1);
2751 absID4 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod-1, irow, AliEMCALGeoParams::fgkEMCALCols-1);
2747966a 2752 }
2753 else
2754 {
e3300002 2755 if( icol < AliEMCALGeoParams::fgkEMCALCols-1 )
2756 absID3 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod, irow, icol+1);
2757 if( icol > 0 )
2758 absID4 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod, irow, icol-1);
2747966a 2759 }
06f1b12a 2760
2761 //Recalibrate cell energy if needed
2762 //Float_t ecell = cells->GetCellAmplitude(absID);
dbba06ca 2763 //GetCaloUtils()->RecalibrateCellAmplitude(ecell,fCalorimeter, absID);
06f1b12a 2764 Double_t tcell = cells->GetCellTime(absID);
dbba06ca 2765 GetCaloUtils()->RecalibrateCellTime(tcell, fCalorimeter, absID,GetReader()->GetInputEvent()->GetBunchCrossNumber());
06f1b12a 2766
2767 Float_t ecell1 = 0, ecell2 = 0, ecell3 = 0, ecell4 = 0;
2768 Double_t tcell1 = 0, tcell2 = 0, tcell3 = 0, tcell4 = 0;
2769
2747966a 2770 if(absID1 >0 )
2771 {
06f1b12a 2772 ecell1 = cells->GetCellAmplitude(absID1);
dbba06ca 2773 GetCaloUtils()->RecalibrateCellAmplitude(ecell1, fCalorimeter, absID1);
06f1b12a 2774 tcell1 = cells->GetCellTime(absID1);
dbba06ca 2775 GetCaloUtils()->RecalibrateCellTime (tcell1, fCalorimeter, absID1,GetReader()->GetInputEvent()->GetBunchCrossNumber());
06f1b12a 2776 }
2747966a 2777 if(absID2 >0 )
2778 {
06f1b12a 2779 ecell2 = cells->GetCellAmplitude(absID2);
dbba06ca 2780 GetCaloUtils()->RecalibrateCellAmplitude(ecell2, fCalorimeter, absID2);
06f1b12a 2781 tcell2 = cells->GetCellTime(absID2);
dbba06ca 2782 GetCaloUtils()->RecalibrateCellTime (tcell2, fCalorimeter, absID2, GetReader()->GetInputEvent()->GetBunchCrossNumber());
06f1b12a 2783 }
2747966a 2784 if(absID3 >0 )
2785 {
06f1b12a 2786 ecell3 = cells->GetCellAmplitude(absID3);
dbba06ca 2787 GetCaloUtils()->RecalibrateCellAmplitude(ecell3, fCalorimeter, absID3);
06f1b12a 2788 tcell3 = cells->GetCellTime(absID3);
dbba06ca 2789 GetCaloUtils()->RecalibrateCellTime (tcell3, fCalorimeter, absID3, GetReader()->GetInputEvent()->GetBunchCrossNumber());
06f1b12a 2790 }
2747966a 2791 if(absID4 >0 )
2792 {
06f1b12a 2793 ecell4 = cells->GetCellAmplitude(absID4);
dbba06ca 2794 GetCaloUtils()->RecalibrateCellAmplitude(ecell4, fCalorimeter, absID4);
06f1b12a 2795 tcell4 = cells->GetCellTime(absID4);
dbba06ca 2796 GetCaloUtils()->RecalibrateCellTime (tcell4, fCalorimeter, absID4, GetReader()->GetInputEvent()->GetBunchCrossNumber());
06f1b12a 2797 }
f1538a5f 2798
2799 if(TMath::Abs(tcell-tcell1)*1.e9 > dtcut) ecell1 = 0 ;
2800 if(TMath::Abs(tcell-tcell2)*1.e9 > dtcut) ecell2 = 0 ;
2801 if(TMath::Abs(tcell-tcell3)*1.e9 > dtcut) ecell3 = 0 ;
2802 if(TMath::Abs(tcell-tcell4)*1.e9 > dtcut) ecell4 = 0 ;
06f1b12a 2803
2804 return ecell1+ecell2+ecell3+ecell4;
1a72f6c5 2805 }
57d8227a 2806 else //PHOS
2807 {
06f1b12a 2808
2809 Int_t absId1 = -1, absId2 = -1, absId3 = -1, absId4 = -1;
2810
2811 Int_t relId1[] = { imod+1, 0, irow+1, icol };
2812 Int_t relId2[] = { imod+1, 0, irow-1, icol };
2813 Int_t relId3[] = { imod+1, 0, irow , icol+1 };
2814 Int_t relId4[] = { imod+1, 0, irow , icol-1 };
2815
2816 GetCaloUtils()->GetPHOSGeometry()->RelToAbsNumbering(relId1, absId1);
2817 GetCaloUtils()->GetPHOSGeometry()->RelToAbsNumbering(relId2, absId2);
2818 GetCaloUtils()->GetPHOSGeometry()->RelToAbsNumbering(relId3, absId3);
2819 GetCaloUtils()->GetPHOSGeometry()->RelToAbsNumbering(relId4, absId4);
2820
2821 Float_t ecell1 = 0, ecell2 = 0, ecell3 = 0, ecell4 = 0;
2822
2823 if(absId1 > 0 ) ecell1 = cells->GetCellAmplitude(absId1);
2824 if(absId2 > 0 ) ecell2 = cells->GetCellAmplitude(absId2);
2825 if(absId3 > 0 ) ecell3 = cells->GetCellAmplitude(absId3);
2826 if(absId4 > 0 ) ecell4 = cells->GetCellAmplitude(absId4);
2827
2828 return ecell1+ecell2+ecell3+ecell4;
2829
1a72f6c5 2830 }
2831
1a72f6c5 2832}
2833
c5693f62 2834//__________________________________________________________________________________________________
649b825d 2835void AliAnaCalorimeterQA::InvariantMassHistograms(const Int_t iclus, const TLorentzVector mom,
c5693f62 2836 const Int_t nModule, const TObjArray* caloClusters,
a82b4462 2837 AliVCaloCells * cells)
649b825d 2838{
2839 // Fill Invariant mass histograms
c8fe2783 2840
649b825d 2841 if(GetDebug()>1) printf("AliAnaCalorimeterQA::InvariantMassHistograms() - Start \n");
3748ffb5 2842
649b825d 2843 //Get vertex for photon momentum calculation and event selection
2844 Double_t v[3] = {0,0,0}; //vertex ;
1a83b960 2845 //GetReader()->GetVertex(v);
a6f26052 2846
649b825d 2847 Int_t nModule2 = -1;
2848 TLorentzVector mom2 ;
2849 Int_t nCaloClusters = caloClusters->GetEntriesFast();
a6f26052 2850
d07278cf 2851 for(Int_t jclus = iclus + 1 ; jclus < nCaloClusters ; jclus++)
2852 {
649b825d 2853 AliVCluster* clus2 = (AliVCluster*)caloClusters->At(jclus);
a82b4462 2854
2855 Float_t maxCellFraction = 0.;
2856 Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(cells, clus2,maxCellFraction);
55c05f8c 2857
d07278cf 2858 // Try to rediuce background with a mild shower shape cut and no more than 1 maxima
2859 // in cluster and remove low energy clusters
2860 if( clus2->GetNCells() <= 1 || !IsGoodCluster(absIdMax,cells) ||
2861 GetCaloUtils()->GetNumberOfLocalMaxima(clus2,cells) > 1 ||
2862 clus2->GetM02() > 0.5 || clus2->E() < 0.3) continue;
c8fe2783 2863
649b825d 2864 //Get cluster kinematics
2865 clus2->GetMomentum(mom2,v);
c8fe2783 2866
649b825d 2867 //Check only certain regions
2868 Bool_t in2 = kTRUE;
2869 if(IsFiducialCutOn()) in2 = GetFiducialCut()->IsInFiducialCut(mom2,fCalorimeter) ;
2870 if(!in2) continue;
2302a644 2871
649b825d 2872 //Get module of cluster
2873 nModule2 = GetModuleNumber(clus2);
c8fe2783 2874
649b825d 2875 //Fill histograms
c8fe2783 2876
649b825d 2877 //All modules
2878 fhIM ->Fill((mom+mom2).Pt(),(mom+mom2).M());
49214ef9 2879
649b825d 2880 //Single module
d07278cf 2881 if(nModule == nModule2 && nModule >=0 && nModule < fNModules)
649b825d 2882 fhIMMod[nModule]->Fill((mom+mom2).Pt(),(mom+mom2).M());
d07278cf 2883
c8fe2783 2884
649b825d 2885 //Asymetry histograms
2886 fhAsym->Fill((mom+mom2).Pt(),TMath::Abs((mom.E()-mom2.E())/(mom.E()+mom2.E())));
2887
2888 }// 2nd cluster loop
2889
2890}
2891
2892//______________________________
2893void AliAnaCalorimeterQA::Init()
2894{
2895 //Check if the data or settings are ok
c8fe2783 2896
649b825d 2897 if(fCalorimeter != "PHOS" && fCalorimeter !="EMCAL")
2898 AliFatal(Form("Wrong calorimeter name <%s>", fCalorimeter.Data()));
521636d2 2899
b6002a83 2900 //if(GetReader()->GetDataType()== AliCaloTrackReader::kMC)
2901 // AliFatal("Analysis of reconstructed data, MC reader not aplicable");
649b825d 2902
2903}
3bfc4732 2904
649b825d 2905//________________________________________
2906void AliAnaCalorimeterQA::InitParameters()
2907{
2908 //Initialize the parameters of the analysis.
2909 AddToHistogramsName("AnaCaloQA_");
2910
2911 fCalorimeter = "EMCAL"; //or PHOS
2912 fNModules = 12; // set maximum to maximum number of EMCAL modules
2913 fNRCU = 2; // set maximum number of RCU in EMCAL per SM
e6fec6f5 2914 fTimeCutMin = -9999999;
2915 fTimeCutMax = 9999999;
649b825d 2916 fEMCALCellAmpMin = 0.2;
2917 fPHOSCellAmpMin = 0.2;
c8fe2783 2918
f1538a5f 2919 // Exotic studies
2920 fExoNECrossCuts = 10 ;
2921 fExoNDTimeCuts = 4 ;
2922
2923 fExoDTimeCuts [0] = 1.e4 ; fExoDTimeCuts [1] = 50.0 ; fExoDTimeCuts [2] = 25.0 ; fExoDTimeCuts [3] = 10.0 ;
2924 fExoECrossCuts[0] = 0.80 ; fExoECrossCuts[1] = 0.85 ; fExoECrossCuts[2] = 0.90 ; fExoECrossCuts[3] = 0.92 ; fExoECrossCuts[4] = 0.94 ;
2925 fExoECrossCuts[5] = 0.95 ; fExoECrossCuts[6] = 0.96 ; fExoECrossCuts[7] = 0.97 ; fExoECrossCuts[8] = 0.98 ; fExoECrossCuts[9] = 0.99 ;
2926
649b825d 2927}
c8fe2783 2928
a82b4462 2929//___________________________________________________________________________________
2930Bool_t AliAnaCalorimeterQA::IsGoodCluster(const Int_t absIdMax, AliVCaloCells* cells)
649b825d 2931{
2932 //Identify cluster as exotic or not
2933
06f1b12a 2934 if(!fStudyBadClusters) return kTRUE;
a82b4462 2935
2747966a 2936 if(fCalorimeter=="EMCAL")
2937 {
06f1b12a 2938 if(!GetCaloUtils()->GetEMCALRecoUtils()->IsRejectExoticCluster())
2747966a 2939 {
06f1b12a 2940 return !( GetCaloUtils()->GetEMCALRecoUtils()->IsExoticCell(absIdMax,cells,(GetReader()->GetInputEvent())->GetBunchCrossNumber()) );
2747966a 2941 }
2942 else
2943 {
06f1b12a 2944 return kTRUE;
2747966a 2945 }
649b825d 2946 }
06f1b12a 2947 else // PHOS
2948 {
1a83b960 2949 Float_t ampMax = cells->GetCellAmplitude(absIdMax);
dbba06ca 2950 GetCaloUtils()->RecalibrateCellAmplitude(ampMax, fCalorimeter, absIdMax);
2747966a 2951
2952 if(ampMax < 0.01) return kFALSE;
2953
1a83b960 2954 if(1-GetECross(absIdMax,cells)/ampMax > 0.95) return kFALSE;
2955 else return kTRUE;
06f1b12a 2956 }
2957
649b825d 2958}
17708df9 2959
a82b4462 2960//_________________________________________________________
649b825d 2961void AliAnaCalorimeterQA::Print(const Option_t * opt) const
2962{
2963 //Print some relevant parameters set for the analysis
2964 if(! opt)
2965 return;
521636d2 2966
649b825d 2967 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
745913ae 2968 AliAnaCaloTrackCorrBaseClass::Print(" ");
2302a644 2969
649b825d 2970 printf("Select Calorimeter %s \n",fCalorimeter.Data());
2971 printf("Time Cut: %3.1f < TOF < %3.1f\n", fTimeCutMin, fTimeCutMax);
2972 printf("EMCAL Min Amplitude : %2.1f GeV/c\n", fEMCALCellAmpMin) ;
2973 printf("PHOS Min Amplitude : %2.1f GeV/c\n", fPHOSCellAmpMin) ;
c8fe2783 2974
649b825d 2975}
2976
649b825d 2977//_____________________________________________________
2978void AliAnaCalorimeterQA::MakeAnalysisFillHistograms()
2979{
2980 //Fill Calorimeter QA histograms
c8fe2783 2981
649b825d 2982 //Play with the MC stack if available
2983 if(IsDataMC()) MCHistograms();
798a9b04 2984
649b825d 2985 //Get List with CaloClusters
2986 TObjArray * caloClusters = NULL;
2987 if (fCalorimeter == "PHOS") caloClusters = GetPHOSClusters();
2988 else if (fCalorimeter == "EMCAL") caloClusters = GetEMCALClusters();
2989 else
2990 AliFatal(Form("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - Wrong calorimeter name <%s>, END\n", fCalorimeter.Data()));
798a9b04 2991
649b825d 2992 // Do not do anything if there are no clusters
2993 if(caloClusters->GetEntriesFast() == 0) return;
521636d2 2994
649b825d 2995 //Get List with CaloCells
2996 AliVCaloCells * cells = 0x0;
2997 if(fCalorimeter == "PHOS") cells = GetPHOSCells();
2998 else cells = GetEMCALCells();
798a9b04 2999
649b825d 3000 if(!caloClusters || !cells) {
3001 AliFatal(Form("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - No CaloClusters or CaloCells available\n"));
3002 return; // trick coverity
c8fe2783 3003 }
649b825d 3004
1a72f6c5 3005 //printf("QA: N cells %d, N clusters %d \n",cells->GetNumberOfCells(),caloClusters->GetEntriesFast());
3006
649b825d 3007 // Correlate Calorimeters and V0 and track Multiplicity
3008 if(fCorrelate) Correlate();
3009
3010 // Clusters
3011 ClusterLoopHistograms(caloClusters,cells);
3012
3013 // Cells
3014 CellHistograms(cells);
3015
3016 if(GetDebug() > 0)
3017 printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - End \n");
3018
a0bb4dc0 3019}
3020
649b825d 3021//______________________________________
3022void AliAnaCalorimeterQA::MCHistograms()
3023{
3024 //Get the MC arrays and do some checks before filling MC histograms
9e9f04cb 3025
649b825d 3026 TLorentzVector mom ;
3027
3028 if(GetReader()->ReadStack()){
9e9f04cb 3029
649b825d 3030 if(!GetMCStack())
3031 AliFatal("Stack not available, is the MC handler called?\n");
9e9f04cb 3032
649b825d 3033 //Fill some pure MC histograms, only primaries.
2747966a 3034 for(Int_t i=0 ; i<GetMCStack()->GetNprimary(); i++)
3035 {//Only primary particles, for all MC transport put GetNtrack()
649b825d 3036 TParticle *primary = GetMCStack()->Particle(i) ;
9e9f04cb 3037
649b825d 3038 if (primary->GetStatusCode() > 11) continue; //Working for PYTHIA and simple generators, check for HERWIG
3039 primary->Momentum(mom);
3040 MCHistograms(mom,TMath::Abs(primary->GetPdgCode()));
3041 } //primary loop
3042 }
3043 else if(GetReader()->ReadAODMCParticles()){
3044
3045 if(!GetReader()->GetAODMCParticles(0))
3046 AliFatal("AODMCParticles not available!");
3047
3048 //Fill some pure MC histograms, only primaries.
2747966a 3049 for(Int_t i=0 ; i < (GetReader()->GetAODMCParticles(0))->GetEntriesFast(); i++)
3050 {
649b825d 3051 AliAODMCParticle *aodprimary = (AliAODMCParticle*) (GetReader()->GetAODMCParticles(0))->At(i) ;
9e9f04cb 3052
649b825d 3053 if (!aodprimary->IsPrimary()) continue; //accept all which is not MC transport generated. Don't know how to avoid partons
9e9f04cb 3054
649b825d 3055 mom.SetPxPyPzE(aodprimary->Px(),aodprimary->Py(),aodprimary->Pz(),aodprimary->E());
3056 MCHistograms(mom,TMath::Abs(aodprimary->GetPdgCode()));
3057 } //primary loop
3058
3059 }
3060}
17708df9 3061
649b825d 3062//_______________________________________________________________________________
3063void AliAnaCalorimeterQA::MCHistograms(const TLorentzVector mom, const Int_t pdg)
3064{
a6f26052 3065 //Fill pure monte carlo related histograms
4865e325 3066
2302a644 3067 Float_t eMC = mom.E();
2302a644 3068 Float_t phiMC = mom.Phi();
3069 if(phiMC < 0)
3070 phiMC += TMath::TwoPi();
3071 Float_t etaMC = mom.Eta();
3072
3073 if (TMath::Abs(etaMC) > 1) return;
3074
35c71d5c 3075 Bool_t in = kFALSE;
3076
3077 //Rough stimate of acceptance for pi0, Eta and electrons
2747966a 3078 if(fCalorimeter == "PHOS")
3079 {
35c71d5c 3080 if(GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter))
3081 in = kTRUE ;
3082 if(GetDebug() > 2) printf("AliAnaCalorimeterQA::MCHistograms() - In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),in);
3083
3084 }
2747966a 3085 else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet())
3086 {
3087 if(GetEMCALGeometry())
3088 {
35c71d5c 3089 Int_t absID=0;
3090 GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(mom.Eta(),mom.Phi(),absID);
3091
3092 if( absID >= 0)
3093 in = kTRUE;
3094
3095 if(GetDebug() > 2) printf("AliAnaCalorimeterQA::MCHistograms() - In %s Real acceptance? %d\n",fCalorimeter.Data(),in);
3096 }
2747966a 3097 else
3098 {
35c71d5c 3099 if(GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter))
3100 in = kTRUE ;
3101 if(GetDebug() > 2) printf("AliAnaCalorimeterQA::MCHistograms() - In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),in);
3102 }
3103 }
2302a644 3104
2747966a 3105 if (pdg==22)
3106 {
c5693f62 3107 fhGenMCE[kmcPhoton] ->Fill(eMC);
3108 if(eMC > 0.5) fhGenMCEtaPhi[kmcPhoton]->Fill(etaMC,phiMC);
2747966a 3109 if(in)
3110 {
c5693f62 3111 fhGenMCAccE[kmcPhoton] ->Fill(eMC);
3112 if(eMC > 0.5) fhGenMCAccEtaPhi[kmcPhoton]->Fill(etaMC,phiMC);
2302a644 3113 }
3114 }
2747966a 3115 else if (pdg==111)
3116 {
c5693f62 3117 fhGenMCE[kmcPi0] ->Fill(eMC);
3118 if(eMC > 0.5) fhGenMCEtaPhi[kmcPi0]->Fill(etaMC,phiMC);
2747966a 3119 if(in)
3120 {
c5693f62 3121 fhGenMCAccE[kmcPi0] ->Fill(eMC);
3122 if(eMC > 0.5) fhGenMCAccEtaPhi[kmcPi0]->Fill(etaMC,phiMC);
2302a644 3123 }
3124 }
2747966a 3125 else if (pdg==221)
3126 {
c5693f62 3127 fhGenMCE[kmcEta] ->Fill(eMC);
3128 if(eMC > 0.5) fhGenMCEtaPhi[kmcEta]->Fill(etaMC,phiMC);
2747966a 3129 if(in)
3130 {
c5693f62 3131 fhGenMCAccE[kmcEta] ->Fill(eMC);
3132 if(eMC > 0.5) fhGenMCAccEtaPhi[kmcEta]->Fill(etaMC,phiMC);
35c71d5c 3133 }
2302a644 3134 }
2747966a 3135 else if (TMath::Abs(pdg)==11)
3136 {
c5693f62 3137 fhGenMCE[kmcElectron] ->Fill(eMC);
3138 if(eMC > 0.5) fhGenMCEtaPhi[kmcElectron]->Fill(etaMC,phiMC);
2747966a 3139 if(in)
3140 {
c5693f62 3141 fhGenMCAccE[kmcElectron] ->Fill(eMC);
3142 if(eMC > 0.5) fhGenMCAccEtaPhi[kmcElectron]->Fill(etaMC,phiMC);
35c71d5c 3143 }
2302a644 3144 }
902aa95c 3145}
c8fe2783 3146
649b825d 3147//_________________________________________________________________________________
3148void AliAnaCalorimeterQA::WeightHistograms(AliVCluster *clus, AliVCaloCells* cells)
3149{
3150 // Calculate weights
3151
3152 // First recalculate energy in case non linearity was applied
3153 Float_t energy = 0;
3154 Float_t ampMax = 0;
2747966a 3155 for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
3156 {
649b825d 3157 Int_t id = clus->GetCellsAbsId()[ipos];
3158
3159 //Recalibrate cell energy if needed
3160 Float_t amp = cells->GetCellAmplitude(id);
dbba06ca 3161 GetCaloUtils()->RecalibrateCellAmplitude(amp, fCalorimeter, id);
649b825d 3162
3163 energy += amp;
3164
3165 if(amp> ampMax)
3166 ampMax = amp;
3167
3168 } // energy loop
3169
2747966a 3170 if(energy <=0 )
3171 {
649b825d 3172 printf("AliAnaCalorimeterQA::WeightHistograms()- Wrong calculated energy %f\n",energy);
3173 return;
3174 }
3175
3176 fhEMaxCellClusterRatio ->Fill(energy,ampMax/energy);
3177 fhEMaxCellClusterLogRatio->Fill(energy,TMath::Log(ampMax/energy));
3178
3179 //Get the ratio and log ratio to all cells in cluster
2747966a 3180 for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
3181 {
649b825d 3182 Int_t id = clus->GetCellsAbsId()[ipos];
3183
3184 //Recalibrate cell energy if needed
3185 Float_t amp = cells->GetCellAmplitude(id);
dbba06ca 3186 GetCaloUtils()->RecalibrateCellAmplitude(amp, fCalorimeter, id);
649b825d 3187
3188 fhECellClusterRatio ->Fill(energy,amp/energy);
3189 fhECellClusterLogRatio->Fill(energy,TMath::Log(amp/energy));
3190 }
3191
3192 //Recalculate shower shape for different W0
2747966a 3193 if(fCalorimeter=="EMCAL")
3194 {
649b825d 3195 Float_t l0org = clus->GetM02();
3196 Float_t l1org = clus->GetM20();
3197 Float_t dorg = clus->GetDispersion();
3198
1a72f6c5 3199 for(Int_t iw = 0; iw < 14; iw++){
3200 GetCaloUtils()->GetEMCALRecoUtils()->SetW0(1+iw*0.5);
649b825d 3201 GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), cells, clus);
3202
3203 fhLambda0ForW0[iw]->Fill(energy,clus->GetM02());
1a72f6c5 3204 //fhLambda1ForW0[iw]->Fill(energy,clus->GetM20());
649b825d 3205
3206 if(IsDataMC()){
3207
3208 Int_t tag = GetMCAnalysisUtils()->CheckOrigin(clus->GetLabels(),clus->GetNLabels(), GetReader(),0);
3209
3210 if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton) &&
3211 !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) &&
3212 !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta) &&
3213 !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion) ){
3214 fhLambda0ForW0MC[iw][0]->Fill(energy,clus->GetM02());
1a72f6c5 3215 //fhLambda1ForW0MC[iw][0]->Fill(energy,clus->GetM20());
649b825d 3216 } // Pure Photon
3217 else if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCElectron) &&
3218 !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion) ){
3219 fhLambda0ForW0MC[iw][1]->Fill(energy,clus->GetM02());
1a72f6c5 3220 //fhLambda1ForW0MC[iw][1]->Fill(energy,clus->GetM20());
649b825d 3221 } // Electron
3222 else if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion) ){
3223 fhLambda0ForW0MC[iw][2]->Fill(energy,clus->GetM02());
1a72f6c5 3224 //fhLambda1ForW0MC[iw][2]->Fill(energy,clus->GetM20());
649b825d 3225 } // Conversion
3226 else if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) ){
3227 fhLambda0ForW0MC[iw][3]->Fill(energy,clus->GetM02());
1a72f6c5 3228 //fhLambda1ForW0MC[iw][3]->Fill(energy,clus->GetM20());
649b825d 3229 }// Pi0
3230 else if(!GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta) &&
3231 !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton) ){
3232 fhLambda0ForW0MC[iw][4]->Fill(energy,clus->GetM02());
1a72f6c5 3233 //fhLambda1ForW0MC[iw][4]->Fill(energy,clus->GetM20());
649b825d 3234 }// Hadron
3235
3236 }// Is MC
3237 } // w0 loop
3238
3239 // Set the original values back
3240 clus->SetM02(l0org);
3241 clus->SetM20(l1org);
3242 clus->SetDispersion(dorg);
3243
3244 }// EMCAL
3245
3246}
3247
3248
3249