]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGGA/CaloTrackCorrelations/AliAnaCalorimeterQA.cxx
add dphi histograms for deta > 0.8
[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;
1073 printf("break\n");
1074 break;
1075 }
1076
649b825d 1077 if(GetDebug() > 1 )printf("\t pdg %d, index %d, %s, status %d \n",pdg, iMother, primary->GetName(),status);
1078 }
008693e5 1079
d07278cf 1080 if(GetDebug() > 1 )
1081 {
649b825d 1082 printf("AliAnaCalorimeterQA::ClusterHistograms() - Converted Cluster mother before conversion: \n");
1083 printf("\t Mother label %d, pdg %d, %s, status %d, parent %d \n",iMother, pdg, primary->GetName(), status, iParent);
1084 }
1085
1086 }
1087
1088 //Overlapped pi0 (or eta, there will be very few), get the meson
1089 if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) ||
d07278cf 1090 GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta))
1091 {
649b825d 1092 if(GetDebug() > 1 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped Meson decay!, Find it: \n");
008693e5 1093
d07278cf 1094 while(pdg != 111 && pdg != 221)
008693e5 1095 {
1096 //printf("iMother %d, pdg %d, iParent %d, pdg %d\n",iMother,pdg,iParent,GetMCStack()->Particle(iParent)->GetPdgCode());
649b825d 1097 iMother = iParent;
1098 primary = GetMCStack()->Particle(iMother);
1099 status = primary->GetStatusCode();
649b825d 1100 pdg = TMath::Abs(primary->GetPdgCode());
008693e5 1101 iParent = primary->GetFirstMother();
1102
1103 if( iParent < 0 )break;
1104
649b825d 1105 if(GetDebug() > 1 ) printf("\t pdg %d, %s, index %d\n",pdg, primary->GetName(),iMother);
008693e5 1106
d07278cf 1107 if(iMother==-1)
1108 {
649b825d 1109 printf("AliAnaCalorimeterQA::ClusterHistograms() - Tagged as Overlapped photon but meson not found, why?\n");
1110 //break;
1111 }
1112 }
008693e5 1113
649b825d 1114 if(GetDebug() > 2 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped %s decay, label %d \n",
1115 primary->GetName(),iMother);
1116 }
1117
1118 eMC = primary->Energy();
1119 ptMC = primary->Pt();
1120 phiMC = primary->Phi();
1121 etaMC = primary->Eta();
1122 pdg = TMath::Abs(primary->GetPdgCode());
1123 charge = (Int_t) TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
1124
1125 }
d07278cf 1126 else if( GetReader()->ReadAODMCParticles() &&
1127 !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCUnknown))
1128 {//it MC AOD and known tag
649b825d 1129 //Get the list of MC particles
1130 if(!GetReader()->GetAODMCParticles(0))
1131 AliFatal("MCParticles not available!");
1132
1133 aodprimary = (AliAODMCParticle*) (GetReader()->GetAODMCParticles(0))->At(label);
1134 iMother = label;
1135 pdg0 = TMath::Abs(aodprimary->GetPdgCode());
1136 pdg = pdg0;
1137 status = aodprimary->IsPrimary();
1138 vxMC = aodprimary->Xv();
1139 vyMC = aodprimary->Yv();
1140 iParent = aodprimary->GetMother();
1141
d07278cf 1142 if(GetDebug() > 1 )
1143 {
649b825d 1144 printf("AliAnaCalorimeterQA::ClusterHistograms() - Cluster most contributing mother: \n");
1145 printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d, parent %d \n",
1146 iMother, pdg0, aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary(), iParent);
1147 }
1148
1149 //Get final particle, no conversion products
008693e5 1150 if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion))
1151 {
649b825d 1152 if(GetDebug() > 1 )
1153 printf("AliAnaCalorimeterQA::ClusterHistograms() - Converted cluster!. Find before conversion: \n");
1154 //Get the parent
1155 aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles(0))->At(iParent);
1156 pdg = TMath::Abs(aodprimary->GetPdgCode());
d07278cf 1157 while ((pdg == 22 || pdg == 11) && !aodprimary->IsPhysicalPrimary())
1158 {
008693e5 1159 Int_t iMotherOrg = iMother;
649b825d 1160 iMother = iParent;
1161 aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles(0))->At(iMother);
1162 status = aodprimary->IsPrimary();
1163 iParent = aodprimary->GetMother();
1164 pdg = TMath::Abs(aodprimary->GetPdgCode());
008693e5 1165
1166 // If gone too back and non stable, assign the decay photon/electron
1167 // there are other possible decays, ignore them for the moment
1168 if(pdg==111 || pdg==221)
1169 {
1170 aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles(0))->At(iMotherOrg);
1171 break;
1172 }
1173
1174 if(iParent < 0 )
1175 {
1176 iParent = iMother;
1177 break;
1178 }
1179
649b825d 1180 if(GetDebug() > 1 )
1181 printf("\t pdg %d, index %d, Primary? %d, Physical Primary? %d \n",
1182 pdg, iMother, aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary());
1183 }
1184
d07278cf 1185 if(GetDebug() > 1 )
1186 {
649b825d 1187 printf("AliAnaCalorimeterQA::ClusterHistograms() - Converted Cluster mother before conversion: \n");
1188 printf("\t Mother label %d, pdg %d, parent %d, Primary? %d, Physical Primary? %d \n",
1189 iMother, pdg, iParent, aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary());
1190 }
1191
1192 }
1193
1194 //Overlapped pi0 (or eta, there will be very few), get the meson
1195 if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) ||
2747966a 1196 GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta))
1197 {
649b825d 1198 if(GetDebug() > 1 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped Meson decay!, Find it: PDG %d, mom %d \n",pdg, iMother);
008693e5 1199 while(pdg != 111 && pdg != 221)
1200 {
649b825d 1201 iMother = iParent;
1202 aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles(0))->At(iMother);
1203 status = aodprimary->IsPrimary();
1204 iParent = aodprimary->GetMother();
1205 pdg = TMath::Abs(aodprimary->GetPdgCode());
008693e5 1206
1207 if(iParent < 0 ) break;
649b825d 1208
1209 if(GetDebug() > 1 ) printf("\t pdg %d, index %d\n",pdg, iMother);
1210
2747966a 1211 if(iMother==-1)
1212 {
649b825d 1213 printf("AliAnaCalorimeterQA::ClusterHistograms() - Tagged as Overlapped photon but meson not found, why?\n");
1214 //break;
1215 }
1216 }
1217
1218 if(GetDebug() > 2 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped %s decay, label %d \n",
1219 aodprimary->GetName(),iMother);
1220 }
1221
1222 status = aodprimary->IsPrimary();
1223 eMC = aodprimary->E();
1224 ptMC = aodprimary->Pt();
1225 phiMC = aodprimary->Phi();
1226 etaMC = aodprimary->Eta();
1227 pdg = TMath::Abs(aodprimary->GetPdgCode());
1228 charge = aodprimary->Charge();
1229
1230 }
1231
1232 //Float_t vz = primary->Vz();
1233 Float_t rVMC = TMath::Sqrt(vxMC*vxMC + vyMC*vyMC);
2747966a 1234 if((pdg == 22 || TMath::Abs(pdg)==11) && status!=1)
1235 {
649b825d 1236 fhEMVxyz ->Fill(vxMC,vyMC);//,vz);
1237 fhEMR ->Fill(e,rVMC);
1238 }
1239
1240 //printf("reco e %f, pt %f, phi %f, eta %f \n", e, pt, phi, eta);
1241 //printf("prim e %f, pt %f, phi %f, eta %f \n", eMC,ptMC,phiMC ,etaMC );
1242 //printf("vertex: vx %f, vy %f, vz %f, r %f \n", vxMC, vyMC, vz, r);
1243
1244 //Overlapped pi0 (or eta, there will be very few)
2747966a 1245 if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0))
1246 {
c5693f62 1247 fhRecoMCE [kmcPi0][matched] ->Fill(e,eMC);
1248 if(e > 0.5 && eMC > 0.5) fhRecoMCEta[kmcPi0][(matched)]->Fill(eta,etaMC);
1249 if(e > 0.5 && eMC > 0.5) fhRecoMCPhi[kmcPi0][(matched)]->Fill(phi,phiMC);
1250 if(eMC > 0) fhRecoMCRatioE [kmcPi0][(matched)]->Fill(e,e/eMC);
1251 fhRecoMCDeltaE [kmcPi0][(matched)]->Fill(e,eMC-e);
1252 fhRecoMCDeltaPhi[kmcPi0][(matched)]->Fill(e,phiMC-phi);
1253 fhRecoMCDeltaEta[kmcPi0][(matched)]->Fill(e,etaMC-eta);
649b825d 1254 }//Overlapped pizero decay
2747966a 1255 else if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta))
1256 {
c5693f62 1257 fhRecoMCE [kmcEta][(matched)] ->Fill(e,eMC);
1258 if(e > 0.5 && eMC > 0.5) fhRecoMCEta[kmcEta][(matched)]->Fill(eta,etaMC);
1259 if(e > 0.5 && eMC > 0.5) fhRecoMCPhi[kmcEta][(matched)]->Fill(phi,phiMC);
1260 if(eMC > 0) fhRecoMCRatioE [kmcEta][(matched)]->Fill(e,e/eMC);
1261 fhRecoMCDeltaE [kmcEta][(matched)]->Fill(e,eMC-e);
1262 fhRecoMCDeltaPhi[kmcEta][(matched)]->Fill(e,phiMC-phi);
1263 fhRecoMCDeltaEta[kmcEta][(matched)]->Fill(e,etaMC-eta);
649b825d 1264 }//Overlapped eta decay
008693e5 1265 else if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton) &&
1266 !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion))
2747966a 1267 {
c5693f62 1268 fhRecoMCE [kmcPhoton][(matched)] ->Fill(e,eMC);
1269 if(e > 0.5 && eMC > 0.5) fhRecoMCEta[kmcPhoton][(matched)]->Fill(eta,etaMC);
1270 if(e > 0.5 && eMC > 0.5) fhRecoMCPhi[kmcPhoton][(matched)]->Fill(phi,phiMC);
1271 if(eMC > 0) fhRecoMCRatioE [kmcPhoton][(matched)]->Fill(e,e/eMC);
008693e5 1272
c5693f62 1273 fhRecoMCDeltaE [kmcPhoton][(matched)]->Fill(e,eMC-e);
1274 fhRecoMCDeltaPhi[kmcPhoton][(matched)]->Fill(e,phiMC-phi);
1275 fhRecoMCDeltaEta[kmcPhoton][(matched)]->Fill(e,etaMC-eta);
649b825d 1276 }//photon
008693e5 1277 else if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCElectron) &&
1278 !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion))
2747966a 1279 {
c5693f62 1280 fhRecoMCE [kmcElectron][(matched)] ->Fill(e,eMC);
1281 if(e > 0.5 && eMC > 0.5) fhRecoMCEta[kmcElectron][(matched)]->Fill(eta,etaMC);
1282 if(e > 0.5 && eMC > 0.5) fhRecoMCPhi[kmcElectron][(matched)]->Fill(phi,phiMC);
1283 if(eMC > 0) fhRecoMCRatioE [kmcElectron][(matched)]->Fill(e,e/eMC);
1284 fhRecoMCDeltaE [kmcElectron][(matched)]->Fill(e,eMC-e);
1285 fhRecoMCDeltaPhi[kmcElectron][(matched)]->Fill(e,phiMC-phi);
1286 fhRecoMCDeltaEta[kmcElectron][(matched)]->Fill(e,etaMC-eta);
649b825d 1287 fhEMVxyz ->Fill(vxMC,vyMC);//,vz);
1288 fhEMR ->Fill(e,rVMC);
1289 }
2747966a 1290 else if(charge == 0)
1291 {
c5693f62 1292 fhRecoMCE [kmcNeHadron][(matched)] ->Fill(e,eMC);
1293 if(e > 0.5 && eMC > 0.5) fhRecoMCEta[kmcNeHadron][(matched)]->Fill(eta,etaMC);
1294 if(e > 0.5 && eMC > 0.5) fhRecoMCPhi[kmcNeHadron][(matched)]->Fill(phi,phiMC);
1295 if(eMC > 0) fhRecoMCRatioE [kmcNeHadron][(matched)]->Fill(e,e/eMC);
1296 fhRecoMCDeltaE [kmcNeHadron][(matched)]->Fill(e,eMC-e);
1297 fhRecoMCDeltaPhi[kmcNeHadron][(matched)]->Fill(e,phiMC-phi);
1298 fhRecoMCDeltaEta[kmcNeHadron][(matched)]->Fill(e,etaMC-eta);
649b825d 1299 fhHaVxyz ->Fill(vxMC,vyMC);//,vz);
1300 fhHaR ->Fill(e,rVMC);
1301 }
2747966a 1302 else if(charge!=0)
1303 {
c5693f62 1304 fhRecoMCE [kmcChHadron][(matched)] ->Fill(e,eMC);
1305 if(e > 0.5 && eMC > 0.5) fhRecoMCEta[kmcChHadron][(matched)]->Fill(eta,etaMC);
1306 if(e > 0.5 && eMC > 0.5) fhRecoMCPhi[kmcChHadron][(matched)]->Fill(phi,phiMC);
1307 if(eMC > 0) fhRecoMCRatioE [kmcChHadron][(matched)]->Fill(e,e/eMC);
1308 fhRecoMCDeltaE [kmcChHadron][(matched)]->Fill(e,eMC-e);
1309 fhRecoMCDeltaPhi[kmcChHadron][(matched)]->Fill(e,phiMC-phi);
1310 fhRecoMCDeltaEta[kmcChHadron][(matched)]->Fill(e,etaMC-eta);
649b825d 1311 fhHaVxyz ->Fill(vxMC,vyMC);//,vz);
1312 fhHaR ->Fill(e,rVMC);
1313 }
1314
1315 if(primary || aodprimary) return kTRUE ;
1316 else return kFALSE;
1317
1318}
1319
1320//________________________________________________________________________________________________
1321void AliAnaCalorimeterQA::ClusterMatchedWithTrackHistograms(AliVCluster *clus, TLorentzVector mom,
1322 const Bool_t okPrimary, const Int_t pdg)
1323{
1324 //Histograms for clusters matched with tracks
42d47cb7 1325
649b825d 1326 Float_t e = mom.E();
1327 Float_t pt = mom.Pt();
1328 Float_t eta = mom.Eta();
1329 Float_t phi = mom.Phi();
1330 if(phi < 0) phi +=TMath::TwoPi();
1331
2747966a 1332 if(fFillAllTH12)
1333 {
649b825d 1334 fhECharged ->Fill(e);
1335 fhPtCharged ->Fill(pt);
1336 fhPhiCharged ->Fill(phi);
1337 fhEtaCharged ->Fill(eta);
1338 }
a82b4462 1339
42d47cb7 1340 //Study the track and matched cluster if track exists.
1341
4bfeae64 1342 AliVTrack *track = GetCaloUtils()->GetMatchedTrack(clus, GetReader()->GetInputEvent());
649b825d 1343
1344 if(!track) return ;
649b825d 1345
a87e069d 1346 Double_t tpt = track->Pt();
1347 Double_t tmom = track->P();
1348 Double_t dedx = track->GetTPCsignal();
1349 Int_t nITS = track->GetNcls(0);
1350 Int_t nTPC = track->GetNcls(1);
1351
1352 // Residuals
1353 Float_t deta = clus->GetTrackDz();
1354 Float_t dphi = clus->GetTrackDx();
1355 Double_t dR = TMath::Sqrt(dphi*dphi + deta*deta);
1356
1357 Double_t eOverP = e/tmom;
1358
d55bb5e1 1359 fh1EOverP->Fill(tpt, eOverP);
1360 if(dR < 0.02) fh1EOverPR02->Fill(tpt,eOverP);
a87e069d 1361
1362 fh2dR->Fill(e,dR);
1363 fh2MatchdEdx->Fill(tmom,dedx);
1364
1365 if(IsDataMC() && okPrimary)
1366 {
1367 Double_t charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
649b825d 1368
a87e069d 1369 if(TMath::Abs(pdg) == 11)
1370 {
d55bb5e1 1371 fhMCEle1EOverP->Fill(tpt,eOverP);
a87e069d 1372 fhMCEle1dR->Fill(dR);
1373 fhMCEle2MatchdEdx->Fill(tmom,dedx);
d55bb5e1 1374 if(dR < 0.02) fhMCEle1EOverPR02->Fill(tpt,eOverP);
a87e069d 1375 }
1376 else if(charge!=0)
1377 {
d55bb5e1 1378 fhMCChHad1EOverP->Fill(tpt,eOverP);
a87e069d 1379 fhMCChHad1dR->Fill(dR);
1380 fhMCChHad2MatchdEdx->Fill(tmom,dedx);
d55bb5e1 1381 if(dR < 0.02) fhMCChHad1EOverPR02->Fill(tpt,eOverP);
a87e069d 1382 }
1383 else if(charge == 0)
1384 {
d55bb5e1 1385 fhMCNeutral1EOverP->Fill(tpt,eOverP);
a87e069d 1386 fhMCNeutral1dR->Fill(dR);
1387 fhMCNeutral2MatchdEdx->Fill(tmom,dedx);
d55bb5e1 1388 if(dR < 0.02) fhMCNeutral1EOverPR02->Fill(tpt,eOverP);
649b825d 1389 }
a87e069d 1390 }//DataMC
1391
1392 if(dR < 0.02 && eOverP > 0.6 && eOverP< 1.2
1393 && clus->GetNCells() > 1 && nITS > 3 && nTPC > 20)
1394 {
1395 fh2EledEdx->Fill(tmom,dedx);
649b825d 1396 }
a87e069d 1397
649b825d 1398}
1399
1400//___________________________________
1401void AliAnaCalorimeterQA::Correlate()
1402{
1403 // Correlate information from PHOS and EMCAL and with V0 and track multiplicity
1404
1405 //Clusters
1406 TObjArray * caloClustersEMCAL = GetEMCALClusters();
1407 TObjArray * caloClustersPHOS = GetPHOSClusters();
1408
1409 Int_t nclEMCAL = caloClustersEMCAL->GetEntriesFast();
1410 Int_t nclPHOS = caloClustersPHOS ->GetEntriesFast();
1411
1412 Float_t sumClusterEnergyEMCAL = 0;
1413 Float_t sumClusterEnergyPHOS = 0;
1414 Int_t iclus = 0;
1415 for(iclus = 0 ; iclus < caloClustersEMCAL->GetEntriesFast() ; iclus++)
1416 sumClusterEnergyEMCAL += ((AliVCluster*)caloClustersEMCAL->At(iclus))->E();
1417 for(iclus = 0 ; iclus < caloClustersPHOS->GetEntriesFast(); iclus++)
1418 sumClusterEnergyPHOS += ((AliVCluster*)caloClustersPHOS->At(iclus))->E();
1419
1420
1421 //Cells
1422
1423 AliVCaloCells * cellsEMCAL = GetEMCALCells();
1424 AliVCaloCells * cellsPHOS = GetPHOSCells();
1425
1426 Int_t ncellsEMCAL = cellsEMCAL->GetNumberOfCells();
1427 Int_t ncellsPHOS = cellsPHOS ->GetNumberOfCells();
1428
1429 Float_t sumCellEnergyEMCAL = 0;
1430 Float_t sumCellEnergyPHOS = 0;
1431 Int_t icell = 0;
1432 for(icell = 0 ; icell < cellsEMCAL->GetNumberOfCells() ; icell++)
1433 sumCellEnergyEMCAL += cellsEMCAL->GetAmplitude(icell);
1434 for(icell = 0 ; icell < cellsPHOS->GetNumberOfCells(); icell++)
1435 sumCellEnergyPHOS += cellsPHOS->GetAmplitude(icell);
1436
1437
1438 //Fill Histograms
1439 fhCaloCorrNClusters->Fill(nclEMCAL,nclPHOS);
1440 fhCaloCorrEClusters->Fill(sumClusterEnergyEMCAL,sumClusterEnergyPHOS);
1441 fhCaloCorrNCells ->Fill(ncellsEMCAL,ncellsPHOS);
1442 fhCaloCorrECells ->Fill(sumCellEnergyEMCAL,sumCellEnergyPHOS);
1443
1444 Int_t v0S = GetV0Signal(0)+GetV0Signal(1);
1445 Int_t v0M = GetV0Multiplicity(0)+GetV0Multiplicity(1);
1446 Int_t trM = GetTrackMultiplicity();
1447 if(fCalorimeter=="PHOS"){
1448 fhCaloV0MCorrNClusters ->Fill(v0M,nclPHOS);
1449 fhCaloV0MCorrEClusters ->Fill(v0M,sumClusterEnergyPHOS);
1450 fhCaloV0MCorrNCells ->Fill(v0M,ncellsPHOS);
1451 fhCaloV0MCorrECells ->Fill(v0M,sumCellEnergyPHOS);
1452
1453 fhCaloV0SCorrNClusters ->Fill(v0S,nclPHOS);
1454 fhCaloV0SCorrEClusters ->Fill(v0S,sumClusterEnergyPHOS);
1455 fhCaloV0SCorrNCells ->Fill(v0S,ncellsPHOS);
1456 fhCaloV0SCorrECells ->Fill(v0S,sumCellEnergyPHOS);
1457
1458 fhCaloTrackMCorrNClusters->Fill(trM,nclPHOS);
1459 fhCaloTrackMCorrEClusters->Fill(trM,sumClusterEnergyPHOS);
1460 fhCaloTrackMCorrNCells ->Fill(trM,ncellsPHOS);
1461 fhCaloTrackMCorrECells ->Fill(trM,sumCellEnergyPHOS);
1462 }
1463 else{
1464 fhCaloV0MCorrNClusters ->Fill(v0M,nclEMCAL);
1465 fhCaloV0MCorrEClusters ->Fill(v0M,sumClusterEnergyEMCAL);
1466 fhCaloV0MCorrNCells ->Fill(v0M,ncellsEMCAL);
1467 fhCaloV0MCorrECells ->Fill(v0M,sumCellEnergyEMCAL);
1468
1469 fhCaloV0SCorrNClusters ->Fill(v0S,nclEMCAL);
1470 fhCaloV0SCorrEClusters ->Fill(v0S,sumClusterEnergyEMCAL);
1471 fhCaloV0SCorrNCells ->Fill(v0S,ncellsEMCAL);
1472 fhCaloV0SCorrECells ->Fill(v0S,sumCellEnergyEMCAL);
1473
1474 fhCaloTrackMCorrNClusters->Fill(trM,nclEMCAL);
1475 fhCaloTrackMCorrEClusters->Fill(trM,sumClusterEnergyEMCAL);
1476 fhCaloTrackMCorrNCells ->Fill(trM,ncellsEMCAL);
1477 fhCaloTrackMCorrECells ->Fill(trM,sumCellEnergyEMCAL);
1478 }
1479
1480 if(GetDebug() > 0 )
1481 {
1482 printf("AliAnaCalorimeterQA::Correlate(): \n");
1483 printf("\t EMCAL: N cells %d, N clusters %d, summed E cells %f, summed E clusters %f \n",
1484 ncellsEMCAL,nclEMCAL, sumCellEnergyEMCAL,sumClusterEnergyEMCAL);
1485 printf("\t PHOS : N cells %d, N clusters %d, summed E cells %f, summed E clusters %f \n",
1486 ncellsPHOS,nclPHOS,sumCellEnergyPHOS,sumClusterEnergyPHOS);
1487 printf("\t V0 : Signal %d, Multiplicity %d, Track Multiplicity %d \n", v0S,v0M,trM);
1488 }
1489
1490}
1491
1492//__________________________________________________
1493TObjString * AliAnaCalorimeterQA::GetAnalysisCuts()
1494{
1495 //Save parameters used for analysis
1496 TString parList ; //this will be list of parameters used for this analysis.
1497 const Int_t buffersize = 255;
1498 char onePar[buffersize] ;
1499
1500 snprintf(onePar,buffersize,"--- AliAnaCalorimeterQA ---\n") ;
1501 parList+=onePar ;
1502 snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
1503 parList+=onePar ;
1504 snprintf(onePar,buffersize,"Time Cut : %2.2f < T < %2.2f ns \n",fTimeCutMin, fTimeCutMax) ;
1505 parList+=onePar ;
1506 snprintf(onePar,buffersize,"PHOS Cell Amplitude > %2.2f GeV, EMCAL Cell Amplitude > %2.2f GeV \n",fPHOSCellAmpMin, fEMCALCellAmpMin) ;
1507 parList+=onePar ;
1508 //Get parameters set in base class.
1509 //parList += GetBaseParametersList() ;
1510
1511 //Get parameters set in FiducialCut class (not available yet)
1512 //parlist += GetFidCut()->GetFidCutParametersList()
1513
1514 return new TObjString(parList) ;
1515}
1516
f1538a5f 1517//___________________________________________________________________________________
1518void AliAnaCalorimeterQA::ExoticHistograms(const Int_t absIdMax, const Float_t ampMax,
1519 AliVCluster *clus, AliVCaloCells* cells)
1520{
1521 // Calculate weights
1522
1523 if(ampMax < 0.01)
1524 {
1525 printf("AliAnaCalorimeterQA::ExoticHistograms()- Low amplitude energy %f\n",ampMax);
1526 return;
1527 }
1528
cc11121e 1529 Float_t l0 = clus->GetM02();
765206a5 1530 Float_t l1 = clus->GetM20();
cc11121e 1531 Float_t en = clus->E();
1532 Int_t nc = clus->GetNCells();
765206a5 1533 Double_t tmax = clus->GetTOF()*1.e9; // recalibrated elsewhere
1534
1535 Float_t eCrossFrac = 1-GetECross(absIdMax,cells, 10000000)/ampMax;
1536
1537 if(en > 5)
1538 {
1539 fhExoL0ECross->Fill(eCrossFrac,l0);
1540 fhExoL1ECross->Fill(eCrossFrac,l1);
1541 }
f1538a5f 1542
1543 for(Int_t ie = 0; ie < fExoNECrossCuts; ie++)
1544 {
1545 for(Int_t idt = 0; idt < fExoNDTimeCuts; idt++)
1546 {
765206a5 1547 eCrossFrac = 1-GetECross(absIdMax,cells, fExoDTimeCuts[idt])/ampMax;
f1538a5f 1548
1549 if(eCrossFrac > fExoECrossCuts[ie])
1550 {
1551 //Exotic
1552 fhExoL0 [ie][idt]->Fill(en,l0 );
765206a5 1553 fhExoL1 [ie][idt]->Fill(en,l1 );
1554 fhExoTime [ie][idt]->Fill(en,tmax);
1555
1556 if(en > 5)
1557 {
1558 fhExoL0NCell[ie][idt]->Fill(nc,l0);
1559 fhExoL1NCell[ie][idt]->Fill(nc,l1);
1560 }
f1538a5f 1561
1562 // Diff time, do for one cut in e cross
1563 if(ie == 0)
1564 {
1565 for (Int_t icell = 0; icell < clus->GetNCells(); icell++)
1566 {
1567 Int_t absId = clus->GetCellsAbsId()[icell];
1568 Double_t time = cells->GetCellTime(absId);
1569 GetCaloUtils()->RecalibrateCellTime(time, fCalorimeter, absId,GetReader()->GetInputEvent()->GetBunchCrossNumber());
1570
1571 Float_t diff = (tmax-time)*1e9;
1572 fhExoDTime[idt]->Fill(en, diff);
1573 }
1574 }
1575 }
1576 else
1577 {
1578 fhExoECross[ie][idt]->Fill(en,eCrossFrac);
cc11121e 1579 fhExoNCell [ie][idt]->Fill(en,nc);
f1538a5f 1580 }
1581 } // D time cut loop
1582 } // e cross cut loop
1583}
1584
649b825d 1585//____________________________________________________
1586TList * AliAnaCalorimeterQA::GetCreateOutputObjects()
1587{
1588 // Create histograms to be saved in output file and
1589 // store them in outputContainer
1590
1591 TList * outputContainer = new TList() ;
1592 outputContainer->SetName("QAHistos") ;
1593
1594 //Histograms
745913ae 1595 Int_t nptbins = GetHistogramRanges()->GetHistoPtBins(); Float_t ptmax = GetHistogramRanges()->GetHistoPtMax(); Float_t ptmin = GetHistogramRanges()->GetHistoPtMin();
1596 Int_t nfineptbins = GetHistogramRanges()->GetHistoFinePtBins(); Float_t ptfinemax = GetHistogramRanges()->GetHistoFinePtMax(); Float_t ptfinemin = GetHistogramRanges()->GetHistoFinePtMin();
1597 Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins(); Float_t phimax = GetHistogramRanges()->GetHistoPhiMax(); Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();
1598 Int_t netabins = GetHistogramRanges()->GetHistoEtaBins(); Float_t etamax = GetHistogramRanges()->GetHistoEtaMax(); Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
1599 Int_t nmassbins = GetHistogramRanges()->GetHistoMassBins(); Float_t massmax = GetHistogramRanges()->GetHistoMassMax(); Float_t massmin = GetHistogramRanges()->GetHistoMassMin();
1600 Int_t nasymbins = GetHistogramRanges()->GetHistoAsymmetryBins(); Float_t asymmax = GetHistogramRanges()->GetHistoAsymmetryMax(); Float_t asymmin = GetHistogramRanges()->GetHistoAsymmetryMin();
d55bb5e1 1601 Int_t nPoverEbins = GetHistogramRanges()->GetHistoPOverEBins(); Float_t EOverPmax = GetHistogramRanges()->GetHistoPOverEMax(); Float_t EOverPmin = GetHistogramRanges()->GetHistoPOverEMin();
745913ae 1602 Int_t ndedxbins = GetHistogramRanges()->GetHistodEdxBins(); Float_t dedxmax = GetHistogramRanges()->GetHistodEdxMax(); Float_t dedxmin = GetHistogramRanges()->GetHistodEdxMin();
1603 Int_t ndRbins = GetHistogramRanges()->GetHistodRBins(); Float_t dRmax = GetHistogramRanges()->GetHistodRMax(); Float_t dRmin = GetHistogramRanges()->GetHistodRMin();
1604 Int_t ntimebins = GetHistogramRanges()->GetHistoTimeBins(); Float_t timemax = GetHistogramRanges()->GetHistoTimeMax(); Float_t timemin = GetHistogramRanges()->GetHistoTimeMin();
1605 Int_t nclbins = GetHistogramRanges()->GetHistoNClustersBins(); Int_t nclmax = GetHistogramRanges()->GetHistoNClustersMax(); Int_t nclmin = GetHistogramRanges()->GetHistoNClustersMin();
1606 Int_t ncebins = GetHistogramRanges()->GetHistoNCellsBins(); Int_t ncemax = GetHistogramRanges()->GetHistoNCellsMax(); Int_t ncemin = GetHistogramRanges()->GetHistoNCellsMin();
1607 Int_t nceclbins = GetHistogramRanges()->GetHistoNClusterCellBins(); Int_t nceclmax = GetHistogramRanges()->GetHistoNClusterCellMax(); Int_t nceclmin = GetHistogramRanges()->GetHistoNClusterCellMin();
1608 Int_t nvdistbins = GetHistogramRanges()->GetHistoVertexDistBins(); Float_t vdistmax = GetHistogramRanges()->GetHistoVertexDistMax(); Float_t vdistmin = GetHistogramRanges()->GetHistoVertexDistMin();
1609 Int_t rbins = GetHistogramRanges()->GetHistoRBins(); Float_t rmax = GetHistogramRanges()->GetHistoRMax(); Float_t rmin = GetHistogramRanges()->GetHistoRMin();
1610 Int_t xbins = GetHistogramRanges()->GetHistoXBins(); Float_t xmax = GetHistogramRanges()->GetHistoXMax(); Float_t xmin = GetHistogramRanges()->GetHistoXMin();
1611 Int_t ybins = GetHistogramRanges()->GetHistoYBins(); Float_t ymax = GetHistogramRanges()->GetHistoYMax(); Float_t ymin = GetHistogramRanges()->GetHistoYMin();
1612 Int_t zbins = GetHistogramRanges()->GetHistoZBins(); Float_t zmax = GetHistogramRanges()->GetHistoZMax(); Float_t zmin = GetHistogramRanges()->GetHistoZMin();
1613 Int_t ssbins = GetHistogramRanges()->GetHistoShowerShapeBins(); Float_t ssmax = GetHistogramRanges()->GetHistoShowerShapeMax(); Float_t ssmin = GetHistogramRanges()->GetHistoShowerShapeMin();
1614 Int_t tdbins = GetHistogramRanges()->GetHistoDiffTimeBins() ; Float_t tdmax = GetHistogramRanges()->GetHistoDiffTimeMax(); Float_t tdmin = GetHistogramRanges()->GetHistoDiffTimeMin();
1615
1616 Int_t nv0sbins = GetHistogramRanges()->GetHistoV0SignalBins(); Int_t nv0smax = GetHistogramRanges()->GetHistoV0SignalMax(); Int_t nv0smin = GetHistogramRanges()->GetHistoV0SignalMin();
1617 Int_t nv0mbins = GetHistogramRanges()->GetHistoV0MultiplicityBins(); Int_t nv0mmax = GetHistogramRanges()->GetHistoV0MultiplicityMax(); Int_t nv0mmin = GetHistogramRanges()->GetHistoV0MultiplicityMin();
1618 Int_t ntrmbins = GetHistogramRanges()->GetHistoTrackMultiplicityBins(); Int_t ntrmmax = GetHistogramRanges()->GetHistoTrackMultiplicityMax(); Int_t ntrmmin = GetHistogramRanges()->GetHistoTrackMultiplicityMin();
649b825d 1619
1620 //EMCAL
1621 fNMaxCols = 48;
1622 fNMaxRows = 24;
1623 fNRCU = 2 ;
1624 //PHOS
f1538a5f 1625 if(fCalorimeter=="PHOS")
1626 {
649b825d 1627 fNMaxCols = 56;
1628 fNMaxRows = 64;
1629 fNRCU = 4 ;
1630 }
1631
1632 fhE = new TH1F ("hE","E reconstructed clusters ", nptbins*5,ptmin,ptmax*5);
1633 fhE->SetXTitle("E (GeV)");
1634 outputContainer->Add(fhE);
1635
f1538a5f 1636 if(fFillAllTH12)
1637 {
649b825d 1638 fhPt = new TH1F ("hPt","p_{T} reconstructed clusters", nptbins,ptmin,ptmax);
1639 fhPt->SetXTitle("p_{T} (GeV/c)");
1640 outputContainer->Add(fhPt);
1641
1642 fhPhi = new TH1F ("hPhi","#phi reconstructed clusters ",nphibins,phimin,phimax);
1643 fhPhi->SetXTitle("#phi (rad)");
1644 outputContainer->Add(fhPhi);
1645
1646 fhEta = new TH1F ("hEta","#eta reconstructed clusters ",netabins,etamin,etamax);
1647 fhEta->SetXTitle("#eta ");
1648 outputContainer->Add(fhEta);
1649 }
1650
f1538a5f 1651 if(fFillAllTH3)
1652 {
649b825d 1653 fhEtaPhiE = new TH3F ("hEtaPhiE","#eta vs #phi vs energy, reconstructed clusters",
1654 netabins,etamin,etamax,nphibins,phimin,phimax,nptbins,ptmin,ptmax);
1655 fhEtaPhiE->SetXTitle("#eta ");
1656 fhEtaPhiE->SetYTitle("#phi (rad)");
1657 fhEtaPhiE->SetZTitle("E (GeV) ");
1658 outputContainer->Add(fhEtaPhiE);
1659 }
1660
1661 fhClusterTimeEnergy = new TH2F ("hClusterTimeEnergy","energy vs TOF, reconstructed clusters",
1662 nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
1663 fhClusterTimeEnergy->SetXTitle("E (GeV) ");
1664 fhClusterTimeEnergy->SetYTitle("TOF (ns)");
1665 outputContainer->Add(fhClusterTimeEnergy);
1666
1667 fhClusterPairDiffTimeE = new TH2F("hClusterPairDiffTimeE","cluster pair time difference vs E, only good clusters",
1668 nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1669 fhClusterPairDiffTimeE->SetXTitle("E_{cluster} (GeV)");
1670 fhClusterPairDiffTimeE->SetYTitle("#Delta t (ns)");
1671 outputContainer->Add(fhClusterPairDiffTimeE);
1672
1673 fhLambda0 = new TH2F ("hLambda0","shower shape, #lambda^{2}_{0} vs E",
1674 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1675 fhLambda0->SetXTitle("E_{cluster}");
1676 fhLambda0->SetYTitle("#lambda^{2}_{0}");
1677 outputContainer->Add(fhLambda0);
1678
1679 fhLambda1 = new TH2F ("hLambda1","shower shape, #lambda^{2}_{1} vs E for bad cluster ",
1680 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1681 fhLambda1->SetXTitle("E_{cluster}");
1682 fhLambda1->SetYTitle("#lambda^{2}_{1}");
1683 outputContainer->Add(fhLambda1);
1684
1685 fhDispersion = new TH2F ("hDispersion","shower shape, Dispersion^{2} vs E for bad cluster ",
1686 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1687 fhDispersion->SetXTitle("E_{cluster}");
1688 fhDispersion->SetYTitle("Dispersion");
1689 outputContainer->Add(fhDispersion);
1690
1691 fhClusterMaxCellCloseCellRatio = new TH2F ("hClusterMaxCellCloseCellRatio","energy vs ratio of max cell / neighbour cell, reconstructed clusters",
1692 nptbins,ptmin,ptmax, 100,0,1.);
1693 fhClusterMaxCellCloseCellRatio->SetXTitle("E_{cluster} (GeV) ");
1694 fhClusterMaxCellCloseCellRatio->SetYTitle("E_{cell i}/E_{cell max}");
1695 outputContainer->Add(fhClusterMaxCellCloseCellRatio);
1696
1697 fhClusterMaxCellCloseCellDiff = new TH2F ("hClusterMaxCellCloseCellDiff","energy vs ratio of max cell / neighbour cell, reconstructed clusters",
1698 nptbins,ptmin,ptmax, 500,0,100.);
1699 fhClusterMaxCellCloseCellDiff->SetXTitle("E_{cluster} (GeV) ");
1700 fhClusterMaxCellCloseCellDiff->SetYTitle("E_{cell max}-E_{cell i} (GeV)");
1701 outputContainer->Add(fhClusterMaxCellCloseCellDiff);
1702
1703 fhClusterMaxCellDiff = new TH2F ("hClusterMaxCellDiff","energy vs difference of cluster energy - max cell energy / cluster energy, good clusters",
1704 nptbins,ptmin,ptmax, 500,0,1.);
1705 fhClusterMaxCellDiff->SetXTitle("E_{cluster} (GeV) ");
1706 fhClusterMaxCellDiff->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
1707 outputContainer->Add(fhClusterMaxCellDiff);
1708
1709 fhClusterMaxCellDiffNoCut = new TH2F ("hClusterMaxCellDiffNoCut","energy vs difference of cluster energy - max cell energy / cluster energy",
1710 nptbins,ptmin,ptmax, 500,0,1.);
1711 fhClusterMaxCellDiffNoCut->SetXTitle("E_{cluster} (GeV) ");
1712 fhClusterMaxCellDiffNoCut->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
1713 outputContainer->Add(fhClusterMaxCellDiffNoCut);
1714
1a72f6c5 1715 fhClusterMaxCellECross = new TH2F ("hClusterMaxCellECross","1 - Energy in cross around max energy cell / max energy cell vs cluster energy, good clusters",
1716 nptbins,ptmin,ptmax, 400,-1,1.);
1717 fhClusterMaxCellECross->SetXTitle("E_{cluster} (GeV) ");
1718 fhClusterMaxCellECross->SetYTitle("1- E_{cross}/E_{cell max}");
1719 outputContainer->Add(fhClusterMaxCellECross);
1720
f1538a5f 1721 fhNCellsPerClusterNoCut = new TH2F ("hNCellsPerClusterNoCut","# cells per cluster vs energy, no bad clusters cut",
1722 nptbins,ptmin,ptmax, nceclbins,nceclmin,nceclmax);
1723 fhNCellsPerClusterNoCut->SetXTitle("E (GeV)");
1724 fhNCellsPerClusterNoCut->SetYTitle("n cells");
1725 outputContainer->Add(fhNCellsPerClusterNoCut);
1726
1727 fhNCellsPerCluster = new TH2F ("hNCellsPerCluster","# cells per cluster vs energy",nptbins,ptmin,ptmax, nceclbins,nceclmin,nceclmax);
1728 fhNCellsPerCluster->SetXTitle("E (GeV)");
1729 fhNCellsPerCluster->SetYTitle("n cells");
1730 outputContainer->Add(fhNCellsPerCluster);
1731
1732 fhNClusters = new TH1F ("hNClusters","# clusters", nclbins,nclmin,nclmax);
1733 fhNClusters->SetXTitle("number of clusters");
1734 outputContainer->Add(fhNClusters);
1735
1a83b960 1736 if(fStudyBadClusters)
1737 {
649b825d 1738 fhBadClusterEnergy = new TH1F ("hBadClusterEnergy","Bad cluster energy", nptbins,ptmin,ptmax);
1739 fhBadClusterEnergy->SetXTitle("E_{cluster} (GeV) ");
1740 outputContainer->Add(fhBadClusterEnergy);
1741
1742 fhBadClusterMaxCellCloseCellRatio = new TH2F ("hBadClusterMaxCellCloseCellRatio","energy vs ratio of max cell / neighbour cell constributing cell, reconstructed bad clusters",
1743 nptbins,ptmin,ptmax, 100,0,1.);
1744 fhBadClusterMaxCellCloseCellRatio->SetXTitle("E_{cluster} (GeV) ");
1745 fhBadClusterMaxCellCloseCellRatio->SetYTitle("ratio");
1746 outputContainer->Add(fhBadClusterMaxCellCloseCellRatio);
1747
1748 fhBadClusterMaxCellCloseCellDiff = new TH2F ("hBadClusterMaxCellCloseCellDiff","energy vs ratio of max cell - neighbour cell constributing cell, reconstructed bad clusters",
1749 nptbins,ptmin,ptmax, 500,0,100);
1750 fhBadClusterMaxCellCloseCellDiff->SetXTitle("E_{cluster} (GeV) ");
1751 fhBadClusterMaxCellCloseCellDiff->SetYTitle("E_{cell max} - E_{cell i} (GeV)");
1752 outputContainer->Add(fhBadClusterMaxCellCloseCellDiff);
1753
1754 fhBadClusterMaxCellDiff = new TH2F ("hBadClusterMaxCellDiff","energy vs difference of cluster energy - max cell energy / cluster energy for bad clusters",
1755 nptbins,ptmin,ptmax, 500,0,1.);
1756 fhBadClusterMaxCellDiff->SetXTitle("E_{cluster} (GeV) ");
1757 fhBadClusterMaxCellDiff->SetYTitle("(E_{cluster} - E_{cell max}) / E_{cluster}");
1758 outputContainer->Add(fhBadClusterMaxCellDiff);
1759
1760 fhBadClusterTimeEnergy = new TH2F ("hBadClusterTimeEnergy","energy vs TOF of reconstructed bad clusters",
1761 nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
1762 fhBadClusterTimeEnergy->SetXTitle("E_{cluster} (GeV) ");
1763 fhBadClusterTimeEnergy->SetYTitle("TOF (ns)");
1764 outputContainer->Add(fhBadClusterTimeEnergy);
1765
1766 fhBadClusterPairDiffTimeE = new TH2F("hBadClusterPairDiffTimeE","cluster pair time difference (bad - good) vs E from bad cluster",nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1767 fhBadClusterPairDiffTimeE->SetXTitle("E_{bad cluster} (GeV)");
1768 fhBadClusterPairDiffTimeE->SetYTitle("#Delta t (ns)");
1769 outputContainer->Add(fhBadClusterPairDiffTimeE);
1770
1a72f6c5 1771 fhBadClusterMaxCellECross = new TH2F ("hBadClusterMaxCellECross","1 - Energy in cross around max energy cell / max energy cell vs cluster energy, bad clusters",
1772 nptbins,ptmin,ptmax, 400,-1,1.);
1773 fhBadClusterMaxCellECross->SetXTitle("E_{cluster} (GeV) ");
1774 fhBadClusterMaxCellECross->SetYTitle("1- E_{cross}/E_{cell max}");
1775 outputContainer->Add(fhBadClusterMaxCellECross);
1776
e6fec6f5 1777 if(fFillAllCellTimeHisto)
1778 {
649b825d 1779 fhBadCellTimeSpreadRespectToCellMax = new TH2F ("hBadCellTimeSpreadRespectToCellMax","t_{cell max}-t_{cell i} from bad cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1780 fhBadCellTimeSpreadRespectToCellMax->SetXTitle("E (GeV)");
1781 fhBadCellTimeSpreadRespectToCellMax->SetYTitle("#Delta t_{cell max - i} (ns)");
1782 outputContainer->Add(fhBadCellTimeSpreadRespectToCellMax);
1783
1784 fhBadClusterMaxCellDiffAverageTime = new TH2F ("hBadClusterMaxCellDiffAverageTime","t_{cell max}-t_{average} from bad cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1785 fhBadClusterMaxCellDiffAverageTime->SetXTitle("E (GeV)");
1786 fhBadClusterMaxCellDiffAverageTime->SetYTitle("#Delta t_{cell max - average} (ns)");
1787 outputContainer->Add(fhBadClusterMaxCellDiffAverageTime);
a82b4462 1788
649b825d 1789 fhBadClusterMaxCellDiffWeightedTime = new TH2F ("hBadClusterMaxCellDiffWeightedTime","t_{cell max}-t_{weighted} from bad cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1790 fhBadClusterMaxCellDiffWeightedTime->SetXTitle("E (GeV)");
1791 fhBadClusterMaxCellDiffWeightedTime->SetYTitle("#Delta t_{cell max - weighted} (ns)");
1792 outputContainer->Add(fhBadClusterMaxCellDiffWeightedTime);
1793
649b825d 1794 }
1795
1796 }
1797
f1538a5f 1798 if(fStudyExotic)
1799 {
765206a5 1800 fhExoL0ECross = new TH2F("hExoL0_ECross",
1801 "#lambda^{2}_{0} vs 1-E_{+}/E_{max} for E > 5 GeV",
1802 nptbins,ptmin,ptmax,tdbins,tdmin,tdmax);
1803 fhExoL0ECross ->SetYTitle("1-E_{+}/E_{cell max}");
1804 fhExoL0ECross ->SetXTitle("#lambda^{2}_{0}");
1805 outputContainer->Add(fhExoL0ECross) ;
1806
1807 fhExoL1ECross = new TH2F("hExoL1_ECross",
1808 "#lambda^{2}_{1} vs 1-E_{+}/E_{max} for E > 5 GeV",
1809 nptbins,ptmin,ptmax,tdbins,tdmin,tdmax);
1810 fhExoL1ECross ->SetYTitle("1-E_{+}/E_{cell max}");
1811 fhExoL1ECross ->SetXTitle("#lambda^{2}_{1}");
1812 outputContainer->Add(fhExoL1ECross) ;
1813
f1538a5f 1814 for(Int_t ie = 0; ie <fExoNECrossCuts; ie++)
1815 {
1816
1817 fhExoDTime[ie] = new TH2F(Form("hExoDTime_ECross%d",ie),
1818 Form("#Delta time = t_{max}-t_{cells} vs E_{cluster} for exotic, 1-E_{+}/E_{max} < %2.2f",fExoECrossCuts[ie]),
1819 nptbins,ptmin,ptmax,tdbins,tdmin,tdmax);
1820 fhExoDTime[ie] ->SetYTitle("#Delta t (ns)");
1821 fhExoDTime[ie] ->SetXTitle("E (GeV)");
1822 outputContainer->Add(fhExoDTime[ie]) ;
1823
1824 for(Int_t idt = 0; idt < fExoNDTimeCuts; idt++)
1825 {
1826 fhExoNCell[ie][idt] = new TH2F(Form("hExoNCell_ECross%d_DT%d",ie,idt),
1827 Form("N cells per cluster vs E cluster, 1-E_{+}/E_{max} < %2.2f, #Delta t < %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
1828 nptbins,ptmin,ptmax,nceclbins,nceclmin,nceclmax);
1829 fhExoNCell[ie][idt] ->SetYTitle("N cells");
1830 fhExoNCell[ie][idt] ->SetXTitle("E (GeV)");
1831 outputContainer->Add(fhExoNCell[ie][idt]) ;
1832
1833 fhExoL0 [ie][idt] = new TH2F(Form("hExoL0_ECross%d_DT%d",ie,idt),
765206a5 1834 Form("#lambda^{2}_{0} vs E cluster for exotic, 1-E_{+}/E_{max} < %2.2f, #Delta t = %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
f1538a5f 1835 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1836 fhExoL0 [ie][idt] ->SetYTitle("#lambda^{2}_{0}");
1837 fhExoL0 [ie][idt] ->SetXTitle("E (GeV)");
1838 outputContainer->Add(fhExoL0[ie][idt]) ;
765206a5 1839
1840 fhExoL1 [ie][idt] = new TH2F(Form("hExoL1_ECross%d_DT%d",ie,idt),
1841 Form("#lambda^{2}_{1} vs E cluster for exotic, 1-E_{+}/E_{max} < %2.2f, #Delta t = %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
1842 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1843 fhExoL1 [ie][idt] ->SetYTitle("#lambda^{2}_{1}");
1844 fhExoL1 [ie][idt] ->SetXTitle("E (GeV)");
1845 outputContainer->Add(fhExoL1[ie][idt]) ;
f1538a5f 1846
1847 fhExoECross[ie][idt] = new TH2F(Form("hExoECross_ECross%d_DT%d",ie,idt),
1848 Form("E cross for cells vs E cell, 1-E_{+}/E_{max} < %2.2f, #Delta t < %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
1849 nptbins,ptmin,ptmax,400,0,1);
1850 fhExoECross[ie][idt] ->SetYTitle("1-E_{+}/E_{cell max}");
1851 fhExoECross[ie][idt] ->SetXTitle("E_{cell} (GeV)");
1852 outputContainer->Add(fhExoECross[ie][idt]) ;
1853
1854 fhExoTime [ie][idt] = new TH2F(Form("hExoTime_ECross%d_DT%d",ie,idt),
1855 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]),
1856 nptbins,ptmin,ptmax,ntimebins,timemin,timemax);
1857 fhExoTime [ie][idt] ->SetYTitle("time_{max} (ns)");
1858 fhExoTime [ie][idt] ->SetXTitle("E (GeV)");
1859 outputContainer->Add(fhExoTime[ie][idt]) ;
1860
765206a5 1861 fhExoL0NCell[ie][idt] = new TH2F(Form("hExoL0_NCell%d_DT%d",ie,idt),
1862 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]),
1863 nptbins,ptmin,ptmax,ntimebins,timemin,timemax);
1864 fhExoL0NCell[ie][idt] ->SetYTitle("N cells");
1865 fhExoL0NCell[ie][idt] ->SetXTitle("#lambda^{2}_{0}");
1866 outputContainer->Add(fhExoL0NCell[ie][idt]) ;
1867
1868 fhExoL1NCell[ie][idt] = new TH2F(Form("hExoL1_NCell%d_DT%d",ie,idt),
1869 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]),
1870 nptbins,ptmin,ptmax,ntimebins,timemin,timemax);
1871 fhExoL1NCell[ie][idt] ->SetYTitle("N cells");
1872 fhExoL1NCell[ie][idt] ->SetXTitle("#lambda^{2}_{1}");
1873 outputContainer->Add(fhExoL1NCell[ie][idt]) ;
1874
f1538a5f 1875 }
1876 }
1877 }
1878
649b825d 1879 // Cluster size in terms of cells
f1538a5f 1880 if(fStudyClustersAsymmetry)
1881 {
649b825d 1882 fhDeltaIEtaDeltaIPhiE0[0] = new TH2F ("hDeltaIEtaDeltaIPhiE0"," Cluster size in columns vs rows for E < 2 GeV, n cells > 3",
1883 50,0,50,50,0,50);
1884 fhDeltaIEtaDeltaIPhiE0[0]->SetXTitle("#Delta Column");
1885 fhDeltaIEtaDeltaIPhiE0[0]->SetYTitle("#Delta Row");
1886 outputContainer->Add(fhDeltaIEtaDeltaIPhiE0[0]);
1887
1888 fhDeltaIEtaDeltaIPhiE2[0] = new TH2F ("hDeltaIEtaDeltaIPhiE2"," Cluster size in columns vs rows for 2 <E < 6 GeV, n cells > 3",
1889 50,0,50,50,0,50);
1890 fhDeltaIEtaDeltaIPhiE2[0]->SetXTitle("#Delta Column");
1891 fhDeltaIEtaDeltaIPhiE2[0]->SetYTitle("#Delta Row");
1892 outputContainer->Add(fhDeltaIEtaDeltaIPhiE2[0]);
1893
1894 fhDeltaIEtaDeltaIPhiE6[0] = new TH2F ("hDeltaIEtaDeltaIPhiE6"," Cluster size in columns vs rows for E > 6 GeV, n cells > 3",
1895 50,0,50,50,0,50);
1896 fhDeltaIEtaDeltaIPhiE6[0]->SetXTitle("#Delta Column");
1897 fhDeltaIEtaDeltaIPhiE6[0]->SetYTitle("#Delta Row");
1898 outputContainer->Add(fhDeltaIEtaDeltaIPhiE6[0]);
1899
1900 fhDeltaIA[0] = new TH2F ("hDeltaIA"," Cluster *asymmetry* in cell units vs E",
1901 nptbins,ptmin,ptmax,21,-1.05,1.05);
1902 fhDeltaIA[0]->SetXTitle("E_{cluster}");
1903 fhDeltaIA[0]->SetYTitle("A_{cell in cluster}");
1904 outputContainer->Add(fhDeltaIA[0]);
1905
1906 fhDeltaIAL0[0] = new TH2F ("hDeltaIAL0"," Cluster *asymmetry* in cell units vs #lambda^{2}_{0}",
1907 ssbins,ssmin,ssmax,21,-1.05,1.05);
1908 fhDeltaIAL0[0]->SetXTitle("#lambda^{2}_{0}");
1909 fhDeltaIAL0[0]->SetYTitle("A_{cell in cluster}");
1910 outputContainer->Add(fhDeltaIAL0[0]);
1911
1912 fhDeltaIAL1[0] = new TH2F ("hDeltaIAL1"," Cluster *asymmetry* in cell units vs #lambda^{2}_{1}",
1913 ssbins,ssmin,ssmax,21,-1.05,1.05);
1914 fhDeltaIAL1[0]->SetXTitle("#lambda^{2}_{1}");
1915 fhDeltaIAL1[0]->SetYTitle("A_{cell in cluster}");
1916 outputContainer->Add(fhDeltaIAL1[0]);
1917
1918 fhDeltaIANCells[0] = new TH2F ("hDeltaIANCells"," Cluster *asymmetry* in cell units vs N cells in cluster",
1919 nceclbins,nceclmin,nceclmax,21,-1.05,1.05);
1920 fhDeltaIANCells[0]->SetXTitle("N_{cell in cluster}");
1921 fhDeltaIANCells[0]->SetYTitle("A_{cell in cluster}");
1922 outputContainer->Add(fhDeltaIANCells[0]);
1923
1924
1925 fhDeltaIEtaDeltaIPhiE0[1] = new TH2F ("hDeltaIEtaDeltaIPhiE0Charged"," Cluster size in columns vs rows for E < 2 GeV, n cells > 3, matched with track",
1926 50,0,50,50,0,50);
1927 fhDeltaIEtaDeltaIPhiE0[1]->SetXTitle("#Delta Column");
1928 fhDeltaIEtaDeltaIPhiE0[1]->SetYTitle("#Delta Row");
1929 outputContainer->Add(fhDeltaIEtaDeltaIPhiE0[1]);
1930
1931 fhDeltaIEtaDeltaIPhiE2[1] = new TH2F ("hDeltaIEtaDeltaIPhiE2Charged"," Cluster size in columns vs rows for 2 <E < 6 GeV, n cells > 3, matched with track",
1932 50,0,50,50,0,50);
1933 fhDeltaIEtaDeltaIPhiE2[1]->SetXTitle("#Delta Column");
1934 fhDeltaIEtaDeltaIPhiE2[1]->SetYTitle("#Delta Row");
1935 outputContainer->Add(fhDeltaIEtaDeltaIPhiE2[1]);
1936
1937 fhDeltaIEtaDeltaIPhiE6[1] = new TH2F ("hDeltaIEtaDeltaIPhiE6Charged"," Cluster size in columns vs rows for E > 6 GeV, n cells > 3, matched with track",
1938 50,0,50,50,0,50);
1939 fhDeltaIEtaDeltaIPhiE6[1]->SetXTitle("#Delta Column");
1940 fhDeltaIEtaDeltaIPhiE6[1]->SetYTitle("#Delta Row");
1941 outputContainer->Add(fhDeltaIEtaDeltaIPhiE6[1]);
1942
1943 fhDeltaIA[1] = new TH2F ("hDeltaIACharged"," Cluster *asymmetry* in cell units vs E, matched with track",
1944 nptbins,ptmin,ptmax,21,-1.05,1.05);
1945 fhDeltaIA[1]->SetXTitle("E_{cluster}");
1946 fhDeltaIA[1]->SetYTitle("A_{cell in cluster}");
1947 outputContainer->Add(fhDeltaIA[1]);
1948
1949 fhDeltaIAL0[1] = new TH2F ("hDeltaIAL0Charged"," Cluster *asymmetry* in cell units vs #lambda^{2}_{0}, matched with track",
1950 ssbins,ssmin,ssmax,21,-1.05,1.05);
1951 fhDeltaIAL0[1]->SetXTitle("#lambda^{2}_{0}");
1952 fhDeltaIAL0[1]->SetYTitle("A_{cell in cluster}");
1953 outputContainer->Add(fhDeltaIAL0[1]);
1954
1955 fhDeltaIAL1[1] = new TH2F ("hDeltaIAL1Charged"," Cluster *asymmetry* in cell units vs #lambda^{2}_{1}, matched with track",
1956 ssbins,ssmin,ssmax,21,-1.05,1.05);
1957 fhDeltaIAL1[1]->SetXTitle("#lambda^{2}_{1}");
1958 fhDeltaIAL1[1]->SetYTitle("A_{cell in cluster}");
1959 outputContainer->Add(fhDeltaIAL1[1]);
1960
1961 fhDeltaIANCells[1] = new TH2F ("hDeltaIANCellsCharged"," Cluster *asymmetry* in cell units vs N cells in cluster, matched with track",
1962 nceclbins,nceclmin,nceclmax,21,-1.05,1.05);
1963 fhDeltaIANCells[1]->SetXTitle("N_{cell in cluster}");
1964 fhDeltaIANCells[1]->SetYTitle("A_{cell in cluster}");
1965 outputContainer->Add(fhDeltaIANCells[1]);
1966
1967 if(IsDataMC()){
1968 TString particle[]={"Photon","Electron","Conversion","Hadron"};
1969 for (Int_t iPart = 0; iPart < 4; iPart++) {
1970
1971 fhDeltaIAMC[iPart] = new TH2F (Form("hDeltaIA_MC%s",particle[iPart].Data()),Form(" Cluster *asymmetry* in cell units vs E, from %s",particle[iPart].Data()),
1972 nptbins,ptmin,ptmax,21,-1.05,1.05);
1973 fhDeltaIAMC[iPart]->SetXTitle("E_{cluster}");
1974 fhDeltaIAMC[iPart]->SetYTitle("A_{cell in cluster}");
1975 outputContainer->Add(fhDeltaIAMC[iPart]);
1976 }
1977 }
f1538a5f 1978
1a83b960 1979 if(fStudyBadClusters)
1980 {
1981 fhBadClusterDeltaIEtaDeltaIPhiE0 = new TH2F ("hBadClusterDeltaIEtaDeltaIPhiE0"," Cluster size in columns vs rows for E < 2 GeV, n cells > 3",
1982 50,0,50,50,0,50);
1983 fhBadClusterDeltaIEtaDeltaIPhiE0->SetXTitle("#Delta Column");
1984 fhBadClusterDeltaIEtaDeltaIPhiE0->SetYTitle("#Delta Row");
1985 outputContainer->Add(fhBadClusterDeltaIEtaDeltaIPhiE0);
1986
1987 fhBadClusterDeltaIEtaDeltaIPhiE2 = new TH2F ("hBadClusterDeltaIEtaDeltaIPhiE2"," Cluster size in columns vs rows for 2 <E < 6 GeV, n cells > 3",
1988 50,0,50,50,0,50);
1989 fhBadClusterDeltaIEtaDeltaIPhiE2->SetXTitle("#Delta Column");
1990 fhBadClusterDeltaIEtaDeltaIPhiE2->SetYTitle("#Delta Row");
1991 outputContainer->Add(fhBadClusterDeltaIEtaDeltaIPhiE2);
1992
1993 fhBadClusterDeltaIEtaDeltaIPhiE6 = new TH2F ("hBadClusterDeltaIEtaDeltaIPhiE6"," Cluster size in columns vs rows for E > 6 GeV, n cells > 3",
1994 50,0,50,50,0,50);
1995 fhBadClusterDeltaIEtaDeltaIPhiE6->SetXTitle("#Delta Column");
1996 fhBadClusterDeltaIEtaDeltaIPhiE6->SetYTitle("#Delta Row");
1997 outputContainer->Add(fhBadClusterDeltaIEtaDeltaIPhiE6);
1998
1999 fhBadClusterDeltaIA = new TH2F ("hBadClusterDeltaIA"," Cluster *asymmetry* in cell units vs E",
2000 nptbins,ptmin,ptmax,21,-1.05,1.05);
2001 fhBadClusterDeltaIA->SetXTitle("E_{cluster}");
2002 fhBadClusterDeltaIA->SetYTitle("A_{cell in cluster}");
2003 outputContainer->Add(fhBadClusterDeltaIA);
2004 }
649b825d 2005 }
2006
f1538a5f 2007 if(fStudyWeight)
2008 {
649b825d 2009 fhECellClusterRatio = new TH2F ("hECellClusterRatio"," cell energy / cluster energy vs cluster energy",
2010 nptbins,ptmin,ptmax, 100,0,1.);
2011 fhECellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
2012 fhECellClusterRatio->SetYTitle("E_{cell i}/E_{cluster}");
2013 outputContainer->Add(fhECellClusterRatio);
2014
2015 fhECellClusterLogRatio = new TH2F ("hECellClusterLogRatio"," Log(cell energy / cluster energy) vs cluster energy",
1a72f6c5 2016 nptbins,ptmin,ptmax, 100,-10,0);
649b825d 2017 fhECellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
1a72f6c5 2018 fhECellClusterLogRatio->SetYTitle("Log(E_{cell i}/E_{cluster})");
649b825d 2019 outputContainer->Add(fhECellClusterLogRatio);
2020
2021 fhEMaxCellClusterRatio = new TH2F ("hEMaxCellClusterRatio"," max cell energy / cluster energy vs cluster energy",
2022 nptbins,ptmin,ptmax, 100,0,1.);
2023 fhEMaxCellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
2024 fhEMaxCellClusterRatio->SetYTitle("E_{max cell}/E_{cluster}");
2025 outputContainer->Add(fhEMaxCellClusterRatio);
2026
2027 fhEMaxCellClusterLogRatio = new TH2F ("hEMaxCellClusterLogRatio"," Log(max cell energy / cluster energy) vs cluster energy",
1a72f6c5 2028 nptbins,ptmin,ptmax, 100,-10,0);
649b825d 2029 fhEMaxCellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
1a72f6c5 2030 fhEMaxCellClusterLogRatio->SetYTitle("Log (E_{max cell}/E_{cluster})");
649b825d 2031 outputContainer->Add(fhEMaxCellClusterLogRatio);
2032
1a72f6c5 2033 for(Int_t iw = 0; iw < 14; iw++){
2034 fhLambda0ForW0[iw] = new TH2F (Form("hLambda0ForW0%d",iw),Form("shower shape, #lambda^{2}_{0} vs E, w0 = %1.1f",1+0.5*iw),
649b825d 2035 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2036 fhLambda0ForW0[iw]->SetXTitle("E_{cluster}");
2037 fhLambda0ForW0[iw]->SetYTitle("#lambda^{2}_{0}");
2038 outputContainer->Add(fhLambda0ForW0[iw]);
2039
1a72f6c5 2040// fhLambda1ForW0[iw] = new TH2F (Form("hLambda1ForW0%d",iw),Form("shower shape, #lambda^{2}_{1} vs E, w0 = %1.1f",1+0.5*iw),
2041// nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2042// fhLambda1ForW0[iw]->SetXTitle("E_{cluster}");
2043// fhLambda1ForW0[iw]->SetYTitle("#lambda^{2}_{1}");
2044// outputContainer->Add(fhLambda1ForW0[iw]);
649b825d 2045
2046 if(IsDataMC()){
2047 TString mcnames[] = {"Photon", "Electron","Conversion","Pi0","Hadron"};
2048 for(Int_t imc = 0; imc < 5; imc++){
2049 fhLambda0ForW0MC[iw][imc] = new TH2F (Form("hLambda0ForW0%d_MC%s",iw,mcnames[imc].Data()),
1a72f6c5 2050 Form("shower shape, #lambda^{2}_{0} vs E, w0 = %1.1f, for MC %s",1+0.5*iw,mcnames[imc].Data()),
649b825d 2051 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2052 fhLambda0ForW0MC[iw][imc]->SetXTitle("E_{cluster}");
2053 fhLambda0ForW0MC[iw][imc]->SetYTitle("#lambda^{2}_{0}");
2054 outputContainer->Add(fhLambda0ForW0MC[iw][imc]);
2055
1a72f6c5 2056// fhLambda1ForW0MC[iw][imc] = new TH2F (Form("hLambda1ForW0%d_MC%s",iw,mcnames[imc].Data()),
2057// Form("shower shape, #lambda^{2}_{1} vs E, w0 = %1.1f, for MC %s",1+0.5*iw,mcnames[imc].Data()),
2058// nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2059// fhLambda1ForW0MC[iw][imc]->SetXTitle("E_{cluster}");
2060// fhLambda1ForW0MC[iw][imc]->SetYTitle("#lambda^{2}_{1}");
2061// outputContainer->Add(fhLambda1ForW0MC[iw][imc]);
649b825d 2062 }
2063 }
f1538a5f 2064 }
649b825d 2065 }
2066
2067 //Track Matching
f1538a5f 2068 if(fFillAllTMHisto)
2069 {
2070 if(fFillAllTH12)
2071 {
649b825d 2072 fhECharged = new TH1F ("hECharged","E reconstructed clusters, matched with track", nptbins,ptmin,ptmax);
2073 fhECharged->SetXTitle("E (GeV)");
2074 outputContainer->Add(fhECharged);
2075
2076 fhPtCharged = new TH1F ("hPtCharged","p_{T} reconstructed clusters, matched with track", nptbins,ptmin,ptmax);
2077 fhPtCharged->SetXTitle("p_{T} (GeV/c)");
2078 outputContainer->Add(fhPtCharged);
2079
2080 fhPhiCharged = new TH1F ("hPhiCharged","#phi reconstructed clusters, matched with track",nphibins,phimin,phimax);
2081 fhPhiCharged->SetXTitle("#phi (rad)");
2082 outputContainer->Add(fhPhiCharged);
55c05f8c 2083
2084 fhEtaCharged = new TH1F ("hEtaCharged","#eta reconstructed clusters, matched with track",netabins,etamin,etamax);
2085 fhEtaCharged->SetXTitle("#eta ");
2086 outputContainer->Add(fhEtaCharged);
2087 }
3748ffb5 2088 if(fFillAllTH3){
2089 fhEtaPhiECharged = new TH3F ("hEtaPhiECharged","#eta vs #phi, reconstructed clusters, matched with track",
521636d2 2090 netabins,etamin,etamax,nphibins,phimin,phimax,nptbins,ptmin,ptmax);
3748ffb5 2091 fhEtaPhiECharged->SetXTitle("#eta ");
2092 fhEtaPhiECharged->SetYTitle("#phi ");
2093 fhEtaPhiECharged->SetZTitle("E (GeV) ");
2094 outputContainer->Add(fhEtaPhiECharged);
2095 }
55c05f8c 2096
d55bb5e1 2097 fh1EOverP = new TH2F("h1EOverP","TRACK matches E/p",nptbins,ptmin,ptmax, nPoverEbins,EOverPmin,EOverPmax);
2098 fh1EOverP->SetYTitle("E/p");
2099 fh1EOverP->SetXTitle("p_{T} (GeV/c)");
2100 outputContainer->Add(fh1EOverP);
3bfc4732 2101
a87e069d 2102 fh2dR = new TH2F("h2dR","TRACK matches dR",nptbins,ptmin,ptmax,ndRbins,dRmin,dRmax);
2103 fh2dR->SetXTitle("#Delta R (rad)");
2104 fh2dR->SetXTitle("E cluster (GeV)");
2105 outputContainer->Add(fh2dR) ;
3bfc4732 2106
2107 fh2MatchdEdx = new TH2F("h2MatchdEdx","dE/dx vs. p for all matches",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
2108 fh2MatchdEdx->SetXTitle("p (GeV/c)");
2109 fh2MatchdEdx->SetYTitle("<dE/dx>");
2110 outputContainer->Add(fh2MatchdEdx);
2111
2112 fh2EledEdx = new TH2F("h2EledEdx","dE/dx vs. p for electrons",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
2113 fh2EledEdx->SetXTitle("p (GeV/c)");
2114 fh2EledEdx->SetYTitle("<dE/dx>");
2115 outputContainer->Add(fh2EledEdx) ;
2116
d55bb5e1 2117 fh1EOverPR02 = new TH2F("h1EOverPR02","TRACK matches E/p, all",nptbins,ptmin,ptmax, nPoverEbins,EOverPmin,EOverPmax);
2118 fh1EOverPR02->SetYTitle("E/p");
2119 fh1EOverPR02->SetXTitle("p_{T} (GeV/c)");
2120 outputContainer->Add(fh1EOverPR02);
55c05f8c 2121 }
2122
f1538a5f 2123 if(fFillAllPi0Histo)
2124 {
9e9f04cb 2125 fhIM = new TH2F ("hIM","Cluster pairs Invariant mass vs reconstructed pair energy, ncell > 1",nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
55c05f8c 2126 fhIM->SetXTitle("p_{T, cluster pairs} (GeV) ");
2127 fhIM->SetYTitle("M_{cluster pairs} (GeV/c^{2})");
2128 outputContainer->Add(fhIM);
a6f26052 2129
55c05f8c 2130 fhAsym = new TH2F ("hAssym","Cluster pairs Asymmetry vs reconstructed pair energy",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax);
2131 fhAsym->SetXTitle("p_{T, cluster pairs} (GeV) ");
2132 fhAsym->SetYTitle("Asymmetry");
2133 outputContainer->Add(fhAsym);
a6f26052 2134 }
649b825d 2135
c8fe2783 2136
f1538a5f 2137 if(fFillAllPosHisto2)
2138 {
2139 if(fFillAllTH3)
2140 {
55c05f8c 2141 fhXYZ = new TH3F ("hXYZ","Cluster: x vs y vs z",xbins,xmin,xmax,ybins,ymin,ymax,zbins,zmin,zmax);
2142 fhXYZ->SetXTitle("x (cm)");
2143 fhXYZ->SetYTitle("y (cm)");
2144 fhXYZ->SetZTitle("z (cm) ");
2145 outputContainer->Add(fhXYZ);
2146 }
2147
35c71d5c 2148 fhXNCells = new TH2F ("hXNCells","Cluster X position vs N Cells per Cluster",xbins,xmin,xmax,nceclbins,nceclmin,nceclmax);
55c05f8c 2149 fhXNCells->SetXTitle("x (cm)");
2150 fhXNCells->SetYTitle("N cells per cluster");
2151 outputContainer->Add(fhXNCells);
2152
35c71d5c 2153 fhZNCells = new TH2F ("hZNCells","Cluster Z position vs N Cells per Cluster",zbins,zmin,zmax,nceclbins,nceclmin,nceclmax);
55c05f8c 2154 fhZNCells->SetXTitle("z (cm)");
2155 fhZNCells->SetYTitle("N cells per cluster");
2156 outputContainer->Add(fhZNCells);
2157
2158 fhXE = new TH2F ("hXE","Cluster X position vs cluster energy",xbins,xmin,xmax,nptbins,ptmin,ptmax);
2159 fhXE->SetXTitle("x (cm)");
2160 fhXE->SetYTitle("E (GeV)");
2161 outputContainer->Add(fhXE);
2162
2163 fhZE = new TH2F ("hZE","Cluster Z position vs cluster energy",zbins,zmin,zmax,nptbins,ptmin,ptmax);
2164 fhZE->SetXTitle("z (cm)");
2165 fhZE->SetYTitle("E (GeV)");
2166 outputContainer->Add(fhZE);
2167
35c71d5c 2168 fhRNCells = new TH2F ("hRNCells","Cluster R position vs N Cells per Cluster",rbins,rmin,rmax,nceclbins,nceclmin,nceclmax);
55c05f8c 2169 fhRNCells->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
2170 fhRNCells->SetYTitle("N cells per cluster");
2171 outputContainer->Add(fhRNCells);
2172
2173
35c71d5c 2174 fhYNCells = new TH2F ("hYNCells","Cluster Y position vs N Cells per Cluster",ybins,ymin,ymax,nceclbins,nceclmin,nceclmax);
55c05f8c 2175 fhYNCells->SetXTitle("y (cm)");
2176 fhYNCells->SetYTitle("N cells per cluster");
2177 outputContainer->Add(fhYNCells);
2178
2179 fhRE = new TH2F ("hRE","Cluster R position vs cluster energy",rbins,rmin,rmax,nptbins,ptmin,ptmax);
2180 fhRE->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
2181 fhRE->SetYTitle("E (GeV)");
2182 outputContainer->Add(fhRE);
2183
2184 fhYE = new TH2F ("hYE","Cluster Y position vs cluster energy",ybins,ymin,ymax,nptbins,ptmin,ptmax);
2185 fhYE->SetXTitle("y (cm)");
2186 fhYE->SetYTitle("E (GeV)");
2187 outputContainer->Add(fhYE);
2188 }
f1538a5f 2189
2190 if(fFillAllPosHisto)
2191 {
a6f26052 2192 fhRCellE = new TH2F ("hRCellE","Cell R position vs cell energy",rbins,rmin,rmax,nptbins,ptmin,ptmax);
2193 fhRCellE->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
2194 fhRCellE->SetYTitle("E (GeV)");
2195 outputContainer->Add(fhRCellE);
2196
2197 fhXCellE = new TH2F ("hXCellE","Cell X position vs cell energy",xbins,xmin,xmax,nptbins,ptmin,ptmax);
2198 fhXCellE->SetXTitle("x (cm)");
2199 fhXCellE->SetYTitle("E (GeV)");
2200 outputContainer->Add(fhXCellE);
2201
2202 fhYCellE = new TH2F ("hYCellE","Cell Y position vs cell energy",ybins,ymin,ymax,nptbins,ptmin,ptmax);
2203 fhYCellE->SetXTitle("y (cm)");
2204 fhYCellE->SetYTitle("E (GeV)");
2205 outputContainer->Add(fhYCellE);
2206
2207 fhZCellE = new TH2F ("hZCellE","Cell Z position vs cell energy",zbins,zmin,zmax,nptbins,ptmin,ptmax);
2208 fhZCellE->SetXTitle("z (cm)");
2209 fhZCellE->SetYTitle("E (GeV)");
2210 outputContainer->Add(fhZCellE);
2211
2212 fhXYZCell = new TH3F ("hXYZCell","Cell : x vs y vs z",xbins,xmin,xmax,ybins,ymin,ymax,zbins,zmin,zmax);
2213 fhXYZCell->SetXTitle("x (cm)");
2214 fhXYZCell->SetYTitle("y (cm)");
2215 fhXYZCell->SetZTitle("z (cm)");
2216 outputContainer->Add(fhXYZCell);
2217
2218
2219 Float_t dx = TMath::Abs(xmin)+TMath::Abs(xmax);
2220 Float_t dy = TMath::Abs(ymin)+TMath::Abs(ymax);
2221 Float_t dz = TMath::Abs(zmin)+TMath::Abs(zmax);
2222 Float_t dr = TMath::Abs(rmin)+TMath::Abs(rmax);
2223
35c71d5c 2224 fhDeltaCellClusterRNCells = new TH2F ("hDeltaCellClusterRNCells","Cluster-Cell R position vs N Cells per Cluster",rbins*2,-dr,dr,nceclbins,nceclmin,nceclmax);
a6f26052 2225 fhDeltaCellClusterRNCells->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
2226 fhDeltaCellClusterRNCells->SetYTitle("N cells per cluster");
2227 outputContainer->Add(fhDeltaCellClusterRNCells);
2228
35c71d5c 2229 fhDeltaCellClusterXNCells = new TH2F ("hDeltaCellClusterXNCells","Cluster-Cell X position vs N Cells per Cluster",xbins*2,-dx,dx,nceclbins,nceclmin,nceclmax);
a6f26052 2230 fhDeltaCellClusterXNCells->SetXTitle("x (cm)");
2231 fhDeltaCellClusterXNCells->SetYTitle("N cells per cluster");
2232 outputContainer->Add(fhDeltaCellClusterXNCells);
2233
35c71d5c 2234 fhDeltaCellClusterYNCells = new TH2F ("hDeltaCellClusterYNCells","Cluster-Cell Y position vs N Cells per Cluster",ybins*2,-dy,dy,nceclbins,nceclmin,nceclmax);
a6f26052 2235 fhDeltaCellClusterYNCells->SetXTitle("y (cm)");
2236 fhDeltaCellClusterYNCells->SetYTitle("N cells per cluster");
2237 outputContainer->Add(fhDeltaCellClusterYNCells);
2238
35c71d5c 2239 fhDeltaCellClusterZNCells = new TH2F ("hDeltaCellClusterZNCells","Cluster-Cell Z position vs N Cells per Cluster",zbins*2,-dz,dz,nceclbins,nceclmin,nceclmax);
a6f26052 2240 fhDeltaCellClusterZNCells->SetXTitle("z (cm)");
2241 fhDeltaCellClusterZNCells->SetYTitle("N cells per cluster");
2242 outputContainer->Add(fhDeltaCellClusterZNCells);
2243
2244 fhDeltaCellClusterRE = new TH2F ("hDeltaCellClusterRE","Cluster-Cell R position vs cluster energy",rbins*2,-dr,dr,nptbins,ptmin,ptmax);
2245 fhDeltaCellClusterRE->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
2246 fhDeltaCellClusterRE->SetYTitle("E (GeV)");
2247 outputContainer->Add(fhDeltaCellClusterRE);
2248
2249 fhDeltaCellClusterXE = new TH2F ("hDeltaCellClusterXE","Cluster-Cell X position vs cluster energy",xbins*2,-dx,dx,nptbins,ptmin,ptmax);
2250 fhDeltaCellClusterXE->SetXTitle("x (cm)");
2251 fhDeltaCellClusterXE->SetYTitle("E (GeV)");
2252 outputContainer->Add(fhDeltaCellClusterXE);
2253
2254 fhDeltaCellClusterYE = new TH2F ("hDeltaCellClusterYE","Cluster-Cell Y position vs cluster energy",ybins*2,-dy,dy,nptbins,ptmin,ptmax);
2255 fhDeltaCellClusterYE->SetXTitle("y (cm)");
2256 fhDeltaCellClusterYE->SetYTitle("E (GeV)");
2257 outputContainer->Add(fhDeltaCellClusterYE);
2258
2259 fhDeltaCellClusterZE = new TH2F ("hDeltaCellClusterZE","Cluster-Cell Z position vs cluster energy",zbins*2,-dz,dz,nptbins,ptmin,ptmax);
2260 fhDeltaCellClusterZE->SetXTitle("z (cm)");
2261 fhDeltaCellClusterZE->SetYTitle("E (GeV)");
2262 outputContainer->Add(fhDeltaCellClusterZE);
2263
2264 fhEtaPhiAmp = new TH3F ("hEtaPhiAmp","Cell #eta vs cell #phi vs cell energy",netabins,etamin,etamax,nphibins,phimin,phimax,nptbins,ptmin,ptmax);
2265 fhEtaPhiAmp->SetXTitle("#eta ");
2266 fhEtaPhiAmp->SetYTitle("#phi (rad)");
2267 fhEtaPhiAmp->SetZTitle("E (GeV) ");
2268 outputContainer->Add(fhEtaPhiAmp);
2269
2302a644 2270 }
c8fe2783 2271
a6f26052 2272 //Calo cells
e3300002 2273 fhNCells = new TH1F ("hNCells","# cells", ncebins,ncemin+0.5,ncemax);
2302a644 2274 fhNCells->SetXTitle("n cells");
2275 outputContainer->Add(fhNCells);
2276
2277 fhAmplitude = new TH1F ("hAmplitude","Cell Energy", nptbins*2,ptmin,ptmax);
2278 fhAmplitude->SetXTitle("Cell Energy (GeV)");
2279 outputContainer->Add(fhAmplitude);
2280
35c71d5c 2281 fhAmpId = new TH2F ("hAmpId","Cell Energy", nfineptbins,ptfinemin,ptfinemax,fNMaxRows*fNMaxCols*fNModules,0,fNMaxRows*fNMaxCols*fNModules);
2302a644 2282 fhAmpId->SetXTitle("Cell Energy (GeV)");
2283 outputContainer->Add(fhAmpId);
2284
e6fec6f5 2285 if(fFillAllCellTimeHisto)
2286 {
649b825d 2287 fhCellTimeSpreadRespectToCellMax = new TH2F ("hCellTimeSpreadRespectToCellMax","t_{cell max}-t_{cell i} per cluster", nptbins,ptmin,ptmax,tdbins,tdmin,tdmax);
924e319f 2288 fhCellTimeSpreadRespectToCellMax->SetXTitle("E (GeV)");
2289 fhCellTimeSpreadRespectToCellMax->SetYTitle("#Delta t_{cell max-i} (ns)");
2302a644 2290 outputContainer->Add(fhCellTimeSpreadRespectToCellMax);
2291
649b825d 2292 fhClusterMaxCellDiffAverageTime = new TH2F ("hClusterMaxCellDiffAverageTime","t_{cell max}-t_{average} per cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
9e9f04cb 2293 fhClusterMaxCellDiffAverageTime->SetXTitle("E (GeV)");
2294 fhClusterMaxCellDiffAverageTime->SetYTitle("#Delta t_{cell max - average} (ns)");
2295 outputContainer->Add(fhClusterMaxCellDiffAverageTime);
a82b4462 2296
649b825d 2297 fhClusterMaxCellDiffWeightedTime = new TH2F ("hClusterMaxCellDiffWeightedTime","t_{cell max}-t_{weighted} per cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
2298 fhClusterMaxCellDiffWeightedTime->SetXTitle("E (GeV)");
2299 fhClusterMaxCellDiffWeightedTime->SetYTitle("#Delta t_{cell max - weighted} (ns)");
2300 outputContainer->Add(fhClusterMaxCellDiffWeightedTime);
35c71d5c 2301
924e319f 2302 fhCellIdCellLargeTimeSpread= new TH1F ("hCellIdCellLargeTimeSpread","Cells with time 100 ns larger than cell max in cluster ",
35c71d5c 2303 fNMaxCols*fNMaxRows*fNModules,0,fNMaxCols*fNMaxRows*fNModules);
2302a644 2304 fhCellIdCellLargeTimeSpread->SetXTitle("Absolute Cell Id");
2305 outputContainer->Add(fhCellIdCellLargeTimeSpread);
521636d2 2306
2302a644 2307 fhTime = new TH1F ("hTime","Cell Time",ntimebins,timemin,timemax);
2308 fhTime->SetXTitle("Cell Time (ns)");
2309 outputContainer->Add(fhTime);
2310
1a72f6c5 2311 fhTimeVz = new TH2F ("hTimeVz","Cell Time vs vertex, amplitude > 0.5 GeV",100, 0, 50,ntimebins,timemin,timemax);
2312 fhTimeVz->SetXTitle("|v_{z}| (cm)");
2313 fhTimeVz->SetYTitle("Cell Time (ns)");
2314 outputContainer->Add(fhTimeVz);
2315
35c71d5c 2316 fhTimeId = new TH2F ("hTimeId","Cell Time vs Absolute Id",
2317 ntimebins,timemin,timemax,fNMaxRows*fNMaxCols*fNModules,0,fNMaxRows*fNMaxCols*fNModules);
2302a644 2318 fhTimeId->SetXTitle("Cell Time (ns)");
2319 fhTimeId->SetYTitle("Cell Absolute Id");
2320 outputContainer->Add(fhTimeId);
2321
2322 fhTimeAmp = new TH2F ("hTimeAmp","Cell Time vs Cell Energy",nptbins*2,ptmin,ptmax,ntimebins,timemin,timemax);
2323 fhTimeAmp->SetYTitle("Cell Time (ns)");
2324 fhTimeAmp->SetXTitle("Cell Energy (GeV)");
2325 outputContainer->Add(fhTimeAmp);
2326
2302a644 2327 }
06f1b12a 2328
2329 fhCellECross = new TH2F ("hCellECross","1 - Energy in cross around cell / cell energy",
2330 nptbins,ptmin,ptmax, 400,-1,1.);
2331 fhCellECross->SetXTitle("E_{cell} (GeV) ");
2332 fhCellECross->SetYTitle("1- E_{cross}/E_{cell}");
2333 outputContainer->Add(fhCellECross);
2334
1a72f6c5 2335
f1538a5f 2336 if(fCorrelate)
2337 {
798a9b04 2338 //PHOS vs EMCAL
35c71d5c 2339 fhCaloCorrNClusters = new TH2F ("hCaloCorrNClusters","# clusters in EMCAL vs PHOS", nclbins,nclmin,nclmax,nclbins,nclmin,nclmax);
2302a644 2340 fhCaloCorrNClusters->SetXTitle("number of clusters in EMCAL");
2341 fhCaloCorrNClusters->SetYTitle("number of clusters in PHOS");
2342 outputContainer->Add(fhCaloCorrNClusters);
2343
798a9b04 2344 fhCaloCorrEClusters = new TH2F ("hCaloCorrEClusters","summed energy of clusters in EMCAL vs PHOS", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2302a644 2345 fhCaloCorrEClusters->SetXTitle("#Sigma E of clusters in EMCAL (GeV)");
2346 fhCaloCorrEClusters->SetYTitle("#Sigma E of clusters in PHOS (GeV)");
2347 outputContainer->Add(fhCaloCorrEClusters);
2348
35c71d5c 2349 fhCaloCorrNCells = new TH2F ("hCaloCorrNCells","# Cells in EMCAL vs PHOS", ncebins,ncemin,ncemax, ncebins,ncemin,ncemax);
649b825d 2350 fhCaloCorrNCells->SetXTitle("number of Cells in EMCAL");
2351 fhCaloCorrNCells->SetYTitle("number of Cells in PHOS");
2352 outputContainer->Add(fhCaloCorrNCells);
2302a644 2353
649b825d 2354 fhCaloCorrECells = new TH2F ("hCaloCorrECells","summed energy of Cells in EMCAL vs PHOS", nptbins*2,ptmin,ptmax*2,nptbins*2,ptmin,ptmax*2);
2355 fhCaloCorrECells->SetXTitle("#Sigma E of Cells in EMCAL (GeV)");
2356 fhCaloCorrECells->SetYTitle("#Sigma E of Cells in PHOS (GeV)");
2357 outputContainer->Add(fhCaloCorrECells);
2302a644 2358
649b825d 2359 //Calorimeter VS V0 signal
2360 fhCaloV0SCorrNClusters = new TH2F ("hCaloV0SNClusters",Form("# clusters in %s vs V0 signal",fCalorimeter.Data()), nv0sbins,nv0smin,nv0smax,nclbins,nclmin,nclmax);
2361 fhCaloV0SCorrNClusters->SetXTitle("V0 signal");
2362 fhCaloV0SCorrNClusters->SetYTitle(Form("number of clusters in %s",fCalorimeter.Data()));
2363 outputContainer->Add(fhCaloV0SCorrNClusters);
2302a644 2364
649b825d 2365 fhCaloV0SCorrEClusters = new TH2F ("hCaloV0SEClusters",Form("summed energy of clusters in %s vs V0 signal",fCalorimeter.Data()), nv0sbins,nv0smin,nv0smax,nptbins,ptmin,ptmax);
2366 fhCaloV0SCorrEClusters->SetXTitle("V0 signal");
2367 fhCaloV0SCorrEClusters->SetYTitle(Form("#Sigma E of clusters in %s (GeV)",fCalorimeter.Data()));
2368 outputContainer->Add(fhCaloV0SCorrEClusters);
2302a644 2369
649b825d 2370 fhCaloV0SCorrNCells = new TH2F ("hCaloV0SNCells",Form("# Cells in %s vs V0 signal",fCalorimeter.Data()), nv0sbins,nv0smin,nv0smax, ncebins,ncemin,ncemax);
2371 fhCaloV0SCorrNCells->SetXTitle("V0 signal");
2372 fhCaloV0SCorrNCells->SetYTitle(Form("number of Cells in %s",fCalorimeter.Data()));
2373 outputContainer->Add(fhCaloV0SCorrNCells);
3bfc4732 2374
649b825d 2375 fhCaloV0SCorrECells = new TH2F ("hCaloV0SECells",Form("summed energy of Cells in %s vs V0 signal",fCalorimeter.Data()), nv0sbins,nv0smin,nv0smax,nptbins,ptmin,ptmax);
2376 fhCaloV0SCorrECells->SetXTitle("V0 signal");
2377 fhCaloV0SCorrECells->SetYTitle(Form("#Sigma E of Cells in %s (GeV)",fCalorimeter.Data()));
2378 outputContainer->Add(fhCaloV0SCorrECells);
3bfc4732 2379
649b825d 2380 //Calorimeter VS V0 multiplicity
2381 fhCaloV0MCorrNClusters = new TH2F ("hCaloV0MNClusters",Form("# clusters in %s vs V0 signal",fCalorimeter.Data()), nv0mbins,nv0mmin,nv0mmax,nclbins,nclmin,nclmax);
2382 fhCaloV0MCorrNClusters->SetXTitle("V0 signal");
2383 fhCaloV0MCorrNClusters->SetYTitle(Form("number of clusters in %s",fCalorimeter.Data()));
2384 outputContainer->Add(fhCaloV0MCorrNClusters);
3bfc4732 2385
649b825d 2386 fhCaloV0MCorrEClusters = new TH2F ("hCaloV0MEClusters",Form("summed energy of clusters in %s vs V0 signal",fCalorimeter.Data()), nv0mbins,nv0mmin,nv0mmax,nptbins,ptmin,ptmax);
2387 fhCaloV0MCorrEClusters->SetXTitle("V0 signal");
2388 fhCaloV0MCorrEClusters->SetYTitle(Form("#Sigma E of clusters in %s (GeV)",fCalorimeter.Data()));
2389 outputContainer->Add(fhCaloV0MCorrEClusters);
3bfc4732 2390
649b825d 2391 fhCaloV0MCorrNCells = new TH2F ("hCaloV0MNCells",Form("# Cells in %s vs V0 signal",fCalorimeter.Data()), nv0mbins,nv0mmin,nv0mmax, ncebins,ncemin,ncemax);
2392 fhCaloV0MCorrNCells->SetXTitle("V0 signal");
2393 fhCaloV0MCorrNCells->SetYTitle(Form("number of Cells in %s",fCalorimeter.Data()));
2394 outputContainer->Add(fhCaloV0MCorrNCells);
3bfc4732 2395
649b825d 2396 fhCaloV0MCorrECells = new TH2F ("hCaloV0MECells",Form("summed energy of Cells in %s vs V0 signal",fCalorimeter.Data()), nv0mbins,nv0mmin,nv0mmax,nptbins,ptmin,ptmax);
2397 fhCaloV0MCorrECells->SetXTitle("V0 signal");
2398 fhCaloV0MCorrECells->SetYTitle(Form("#Sigma E of Cells in %s (GeV)",fCalorimeter.Data()));
2399 outputContainer->Add(fhCaloV0MCorrECells);
3bfc4732 2400
649b825d 2401 //Calorimeter VS Track multiplicity
2402 fhCaloTrackMCorrNClusters = new TH2F ("hCaloTrackMNClusters",Form("# clusters in %s vs # tracks",fCalorimeter.Data()), ntrmbins,ntrmmin,ntrmmax,nclbins,nclmin,nclmax);
2403 fhCaloTrackMCorrNClusters->SetXTitle("# tracks");
2404 fhCaloTrackMCorrNClusters->SetYTitle(Form("number of clusters in %s",fCalorimeter.Data()));
2405 outputContainer->Add(fhCaloTrackMCorrNClusters);
3bfc4732 2406
649b825d 2407 fhCaloTrackMCorrEClusters = new TH2F ("hCaloTrackMEClusters",Form("summed energy of clusters in %s vs # tracks",fCalorimeter.Data()), ntrmbins,ntrmmin,ntrmmax,nptbins,ptmin,ptmax);
2408 fhCaloTrackMCorrEClusters->SetXTitle("# tracks");
2409 fhCaloTrackMCorrEClusters->SetYTitle(Form("#Sigma E of clusters in %s (GeV)",fCalorimeter.Data()));
2410 outputContainer->Add(fhCaloTrackMCorrEClusters);
3bfc4732 2411
649b825d 2412 fhCaloTrackMCorrNCells = new TH2F ("hCaloTrackMNCells",Form("# Cells in %s vs # tracks",fCalorimeter.Data()), ntrmbins,ntrmmin,ntrmmax, ncebins,ncemin,ncemax);
2413 fhCaloTrackMCorrNCells->SetXTitle("# tracks");
2414 fhCaloTrackMCorrNCells->SetYTitle(Form("number of Cells in %s",fCalorimeter.Data()));
2415 outputContainer->Add(fhCaloTrackMCorrNCells);
3bfc4732 2416
649b825d 2417 fhCaloTrackMCorrECells = new TH2F ("hCaloTrackMECells",Form("summed energy of Cells in %s vs # tracks",fCalorimeter.Data()), ntrmbins,ntrmmin,ntrmmax,nptbins,ptmin,ptmax);
2418 fhCaloTrackMCorrECells->SetXTitle("# tracks");
2419 fhCaloTrackMCorrECells->SetYTitle(Form("#Sigma E of Cells in %s (GeV)",fCalorimeter.Data()));
2420 outputContainer->Add(fhCaloTrackMCorrECells);
3bfc4732 2421
3bfc4732 2422
649b825d 2423 }//correlate calorimeters
c8fe2783 2424
649b825d 2425 //Module histograms
9725fd2a 2426
649b825d 2427 fhEMod = new TH2F ("hE_Mod","Cluster reconstructed Energy in each present Module",nptbins,ptmin,ptmax,fNModules,0,fNModules);
2428 fhEMod->SetXTitle("E (GeV)");
2429 fhEMod->SetYTitle("Module");
2430 outputContainer->Add(fhEMod);
c8fe2783 2431
649b825d 2432 fhAmpMod = new TH2F ("hAmp_Mod","Cell energy in each present Module",nptbins,ptmin,ptmax,fNModules,0,fNModules);
2433 fhAmpMod->SetXTitle("E (GeV)");
2434 fhAmpMod->SetYTitle("Module");
2435 outputContainer->Add(fhAmpMod);
3bfc4732 2436
e6fec6f5 2437 if(fFillAllCellTimeHisto)
2438 {
0fb69ade 2439 fhTimeMod = new TH2F ("hTime_Mod","Cell time in each present Module",ntimebins,timemin,timemax,fNModules,0,fNModules);
2440 fhTimeMod->SetXTitle("t (ns)");
2441 fhTimeMod->SetYTitle("Module");
2442 outputContainer->Add(fhTimeMod);
2443 }
3bfc4732 2444
e3300002 2445 fhNClustersMod = new TH2F ("hNClusters_Mod","# clusters vs Module", nclbins,nclmin+0.5,nclmax,fNModules,0,fNModules);
649b825d 2446 fhNClustersMod->SetXTitle("number of clusters");
2447 fhNClustersMod->SetYTitle("Module");
2448 outputContainer->Add(fhNClustersMod);
2302a644 2449
e3300002 2450 fhNCellsMod = new TH2F ("hNCells_Mod","# cells vs Module", ncebins,ncemin+0.5,ncemax,fNModules,0,fNModules);
649b825d 2451 fhNCellsMod->SetXTitle("n cells");
2452 fhNCellsMod->SetYTitle("Module");
2453 outputContainer->Add(fhNCellsMod);
17708df9 2454
649b825d 2455 Int_t colmaxs = fNMaxCols;
2456 Int_t rowmaxs = fNMaxRows;
e6fec6f5 2457 if(fCalorimeter=="EMCAL")
2458 {
649b825d 2459 colmaxs=2*fNMaxCols;
2460 rowmaxs=Int_t(fNModules/2)*fNMaxRows;
2461 }
e6fec6f5 2462 else
2463 {
649b825d 2464 rowmaxs=fNModules*fNMaxRows;
2302a644 2465 }
3bfc4732 2466
649b825d 2467 fhGridCells = new TH2F ("hGridCells",Form("Entries in grid of cells"),
2468 colmaxs+2,-1.5,colmaxs+0.5, rowmaxs+2,-1.5,rowmaxs+0.5);
2469 fhGridCells->SetYTitle("row (phi direction)");
2470 fhGridCells->SetXTitle("column (eta direction)");
2471 outputContainer->Add(fhGridCells);
3bfc4732 2472
649b825d 2473 fhGridCellsE = new TH2F ("hGridCellsE","Accumulated energy in grid of cells",
2474 colmaxs+2,-1.5,colmaxs+0.5, rowmaxs+2,-1.5,rowmaxs+0.5);
2475 fhGridCellsE->SetYTitle("row (phi direction)");
2476 fhGridCellsE->SetXTitle("column (eta direction)");
2477 outputContainer->Add(fhGridCellsE);
3bfc4732 2478
e6fec6f5 2479 if(fFillAllCellTimeHisto)
2480 {
0fb69ade 2481 fhGridCellsTime = new TH2F ("hGridCellsTime","Accumulated time in grid of cells",
2482 colmaxs+2,-1.5,colmaxs+0.5, rowmaxs+2,-1.5,rowmaxs+0.5);
2483 fhGridCellsTime->SetYTitle("row (phi direction)");
2484 fhGridCellsTime->SetXTitle("column (eta direction)");
2485 outputContainer->Add(fhGridCellsTime);
2486 }
3bfc4732 2487
e6fec6f5 2488 fhNCellsPerClusterMod = new TH2F*[fNModules];
649b825d 2489 fhNCellsPerClusterModNoCut = new TH2F*[fNModules];
e6fec6f5 2490 fhIMMod = new TH2F*[fNModules];
2491 if(fFillAllCellTimeHisto) fhTimeAmpPerRCU = new TH2F*[fNModules*fNRCU];
0fb69ade 2492
e6fec6f5 2493 for(Int_t imod = 0; imod < fNModules; imod++)
2494 {
649b825d 2495 fhNCellsPerClusterMod[imod] = new TH2F (Form("hNCellsPerCluster_Mod%d",imod),
2496 Form("# cells per cluster vs cluster energy in Module %d",imod),
2497 nptbins,ptmin,ptmax, nceclbins,nceclmin,nceclmax);
2498 fhNCellsPerClusterMod[imod]->SetXTitle("E (GeV)");
2499 fhNCellsPerClusterMod[imod]->SetYTitle("n cells");
2500 outputContainer->Add(fhNCellsPerClusterMod[imod]);
3bfc4732 2501
649b825d 2502 fhNCellsPerClusterModNoCut[imod] = new TH2F (Form("hNCellsPerClusterNoCut_Mod%d",imod),
2503 Form("# cells per cluster vs cluster energy in Module %d, no cut",imod),
2504 nptbins,ptmin,ptmax, nceclbins,nceclmin,nceclmax);
2505 fhNCellsPerClusterModNoCut[imod]->SetXTitle("E (GeV)");
2506 fhNCellsPerClusterModNoCut[imod]->SetYTitle("n cells");
2507 outputContainer->Add(fhNCellsPerClusterModNoCut[imod]);
3bfc4732 2508
e6fec6f5 2509 if(fFillAllCellTimeHisto)
2510 {
2511 for(Int_t ircu = 0; ircu < fNRCU; ircu++)
2512 {
649b825d 2513 fhTimeAmpPerRCU[imod*fNRCU+ircu] = new TH2F (Form("hTimeAmp_Mod%d_RCU%d",imod,ircu),
2514 Form("Cell Energy vs Cell Time in Module %d, RCU %d ",imod,ircu),
2515 nptbins,ptmin,ptmax,ntimebins,timemin,timemax);
2516 fhTimeAmpPerRCU[imod*fNRCU+ircu]->SetXTitle("E (GeV)");
2517 fhTimeAmpPerRCU[imod*fNRCU+ircu]->SetYTitle("time (ns)");
2518 outputContainer->Add(fhTimeAmpPerRCU[imod*fNRCU+ircu]);
3bfc4732 2519
649b825d 2520 }
3bfc4732 2521 }
2522
e6fec6f5 2523 if(fFillAllPi0Histo)
2524 {
649b825d 2525 fhIMMod[imod] = new TH2F (Form("hIM_Mod%d",imod),
2526 Form("Cluster pairs Invariant mass vs reconstructed pair energy in Module %d, n cell > 1",imod),
2527 nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
2528 fhIMMod[imod]->SetXTitle("p_{T, cluster pairs} (GeV) ");
2529 fhIMMod[imod]->SetYTitle("M_{cluster pairs} (GeV/c^{2})");
2530 outputContainer->Add(fhIMMod[imod]);
3bfc4732 2531
649b825d 2532 }
2533 }
2534
a82b4462 2535 // Monte Carlo Histograms
649b825d 2536
2537 TString particleName[] = { "Photon", "Pi0", "Eta", "Electron", "NeutralHadron", "ChargedHadron" };
2538
f1538a5f 2539 if(IsDataMC())
2540 {
2541 for(Int_t iPart = 0; iPart < 6; iPart++)
2542 {
2543 for(Int_t iCh = 0; iCh < 2; iCh++)
2544 {
649b825d 2545 fhRecoMCRatioE[iPart][iCh] = new TH2F (Form("hRecoMCRatioE_%s_Match%d",particleName[iPart].Data(),iCh),
2546 Form("Reco/Gen E, %s, Matched %d",particleName[iPart].Data(),iCh),
2547 nptbins, ptmin, ptmax, 200,0,2);
2548 fhRecoMCRatioE[iPart][iCh]->SetXTitle("E_{reconstructed}/E_{generated}");
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),
2553 Form("MC - Reco E, %s, Matched %d",particleName[iPart].Data(),iCh),
2554 nptbins, ptmin, ptmax, nptbins*2,-ptmax,ptmax);
2555 fhRecoMCDeltaE[iPart][iCh]->SetXTitle("#Delta E (GeV)");
2556 outputContainer->Add(fhRecoMCDeltaE[iPart][iCh]);
3bfc4732 2557
649b825d 2558 fhRecoMCDeltaPhi[iPart][iCh] = new TH2F (Form("hRecoMCDeltaPhi_%s_Match%d",particleName[iPart].Data(),iCh),
2559 Form("MC - Reco #phi, %s, Matched %d",particleName[iPart].Data(),iCh),
2560 nptbins, ptmin, ptmax, nphibins*2,-phimax,phimax);
2561 fhRecoMCDeltaPhi[iPart][iCh]->SetXTitle("#Delta #phi (rad)");
2562 outputContainer->Add(fhRecoMCDeltaPhi[iPart][iCh]);
3bfc4732 2563
649b825d 2564 fhRecoMCDeltaEta[iPart][iCh] = new TH2F (Form("hRecoMCDeltaEta_%s_Match%d",particleName[iPart].Data(),iCh),
2565 Form("MC- Reco #eta, %s, Matched %d",particleName[iPart].Data(),iCh),
2566 nptbins, ptmin, ptmax,netabins*2,-etamax,etamax);
2567 fhRecoMCDeltaEta[iPart][iCh]->SetXTitle("#Delta #eta ");
2568 outputContainer->Add(fhRecoMCDeltaEta[iPart][iCh]);
e1e62b89 2569
649b825d 2570 fhRecoMCE[iPart][iCh] = new TH2F (Form("hRecoMCE_%s_Match%d",particleName[iPart].Data(),iCh),
2571 Form("E distribution, reconstructed vs generated, %s, Matched %d",particleName[iPart].Data(),iCh),
2572 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2573 fhRecoMCE[iPart][iCh]->SetXTitle("E_{rec} (GeV)");
2574 fhRecoMCE[iPart][iCh]->SetYTitle("E_{gen} (GeV)");
2575 outputContainer->Add(fhRecoMCE[iPart][iCh]);
3bfc4732 2576
649b825d 2577 fhRecoMCPhi[iPart][iCh] = new TH2F (Form("hRecoMCPhi_%s_Match%d",particleName[iPart].Data(),iCh),
2578 Form("#phi distribution, reconstructed vs generated, %s, Matched %d",particleName[iPart].Data(),iCh),
2579 nphibins,phimin,phimax, nphibins,phimin,phimax);
2580 fhRecoMCPhi[iPart][iCh]->SetXTitle("#phi_{rec} (rad)");
2581 fhRecoMCPhi[iPart][iCh]->SetYTitle("#phi_{gen} (rad)");
2582 outputContainer->Add(fhRecoMCPhi[iPart][iCh]);
e1e62b89 2583
649b825d 2584 fhRecoMCEta[iPart][iCh] = new TH2F (Form("hRecoMCEta_%s_Match%d",particleName[iPart].Data(),iCh),
2585 Form("#eta distribution, reconstructed vs generated, %s, Matched %d",particleName[iPart].Data(),iCh),
2586 netabins,etamin,etamax,netabins,etamin,etamax);
2587 fhRecoMCEta[iPart][iCh]->SetXTitle("#eta_{rec} ");
2588 fhRecoMCEta[iPart][iCh]->SetYTitle("#eta_{gen} ");
2589 outputContainer->Add(fhRecoMCEta[iPart][iCh]);
3bfc4732 2590 }
649b825d 2591 }
c8fe2783 2592
649b825d 2593 //Pure MC
f1538a5f 2594 for(Int_t iPart = 0; iPart < 4; iPart++)
2595 {
649b825d 2596 fhGenMCE[iPart] = new TH1F(Form("hGenMCE_%s",particleName[iPart].Data()) ,
2597 Form("p_{T} of generated %s",particleName[iPart].Data()),
2598 nptbins,ptmin,ptmax);
2599 fhGenMCEtaPhi[iPart] = new TH2F(Form("hGenMCEtaPhi_%s",particleName[iPart].Data()),
2600 Form("Y vs #phi of generated %s",particleName[iPart].Data()),
2601 netabins,etamin,etamax,nphibins,phimin,phimax);
2602
2603 fhGenMCE[iPart] ->SetXTitle("p_{T} (GeV/c)");
2604 fhGenMCEtaPhi[iPart]->SetXTitle("#eta");
2605 fhGenMCEtaPhi[iPart]->SetYTitle("#phi (rad)");
3bfc4732 2606
649b825d 2607 outputContainer->Add(fhGenMCE[iPart]);
2608 outputContainer->Add(fhGenMCEtaPhi[iPart]);
521636d2 2609
521636d2 2610
649b825d 2611 fhGenMCAccE[iPart] = new TH1F(Form("hGenMCAccE_%s",particleName[iPart].Data()) ,
2612 Form("p_{T} of generated %s",particleName[iPart].Data()),
2613 nptbins,ptmin,ptmax);
2614 fhGenMCAccEtaPhi[iPart] = new TH2F(Form("hGenMCAccEtaPhi_%s",particleName[iPart].Data()),
2615 Form("Y vs #phi of generated %s",particleName[iPart].Data()),
2616 netabins,etamin,etamax,nphibins,phimin,phimax);
2617
2618 fhGenMCAccE[iPart] ->SetXTitle("p_{T} (GeV/c)");
2619 fhGenMCAccEtaPhi[iPart]->SetXTitle("#eta");
2620 fhGenMCAccEtaPhi[iPart]->SetYTitle("#phi (rad)");
2621
2622 outputContainer->Add(fhGenMCAccE[iPart]);
2623 outputContainer->Add(fhGenMCAccEtaPhi[iPart]);
2624
2625 }
f5036bcb 2626
649b825d 2627 //Vertex of generated particles
2628
2629 fhEMVxyz = new TH2F ("hEMVxyz","Production vertex of reconstructed ElectroMagnetic particles",nvdistbins,vdistmin,vdistmax,nvdistbins,vdistmin,vdistmax);//,100,0,500);
2630 fhEMVxyz->SetXTitle("v_{x}");
2631 fhEMVxyz->SetYTitle("v_{y}");
2632 //fhEMVxyz->SetZTitle("v_{z}");
2633 outputContainer->Add(fhEMVxyz);
2634
2635 fhHaVxyz = new TH2F ("hHaVxyz","Production vertex of reconstructed hadrons",nvdistbins,vdistmin,vdistmax,nvdistbins,vdistmin,vdistmax);//,100,0,500);
2636 fhHaVxyz->SetXTitle("v_{x}");
2637 fhHaVxyz->SetYTitle("v_{y}");
2638 //fhHaVxyz->SetZTitle("v_{z}");
2639 outputContainer->Add(fhHaVxyz);
2640
2641 fhEMR = new TH2F ("hEMR","Distance to production vertex of reconstructed ElectroMagnetic particles vs E rec",nptbins,ptmin,ptmax,nvdistbins,vdistmin,vdistmax);
2642 fhEMR->SetXTitle("E (GeV)");
2643 fhEMR->SetYTitle("TMath::Sqrt(v_{x}^{2}+v_{y}^{2})");
2644 outputContainer->Add(fhEMR);
2645
2646 fhHaR = new TH2F ("hHaR","Distance to production vertex of reconstructed Hadrons vs E rec",nptbins,ptmin,ptmax,nvdistbins,vdistmin,vdistmax);
2647 fhHaR->SetXTitle("E (GeV)");
2648 fhHaR->SetYTitle("TMath::Sqrt(v_{x}^{2}+v_{y}^{2})");
2649 outputContainer->Add(fhHaR);
2650
2651
2652 //Track Matching
2653
d55bb5e1 2654 fhMCEle1EOverP = new TH2F("hMCEle1EOverP","TRACK matches E/p, MC electrons",nptbins,ptmin,ptmax, nPoverEbins,EOverPmin,EOverPmax);
2655 fhMCEle1EOverP->SetYTitle("E/p");
2656 fhMCEle1EOverP->SetXTitle("p_{T} (GeV/c)");
2657 outputContainer->Add(fhMCEle1EOverP);
649b825d 2658
2659 fhMCEle1dR = new TH1F("hMCEle1dR","TRACK matches dR, MC electrons",ndRbins,dRmin,dRmax);
2660 fhMCEle1dR->SetXTitle("#Delta R (rad)");
2661 outputContainer->Add(fhMCEle1dR) ;
2662
2663 fhMCEle2MatchdEdx = new TH2F("hMCEle2MatchdEdx","dE/dx vs. p for all matches, MC electrons",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
2664 fhMCEle2MatchdEdx->SetXTitle("p (GeV/c)");
2665 fhMCEle2MatchdEdx->SetYTitle("<dE/dx>");
2666 outputContainer->Add(fhMCEle2MatchdEdx);
2667
d55bb5e1 2668 fhMCChHad1EOverP = new TH2F("hMCChHad1EOverP","TRACK matches E/p, MC charged hadrons",nptbins,ptmin,ptmax, nPoverEbins,EOverPmin,EOverPmax);
2669 fhMCChHad1EOverP->SetYTitle("E/p");
2670 fhMCChHad1EOverP->SetXTitle("p_{T} (GeV/c)");
2671 outputContainer->Add(fhMCChHad1EOverP);
649b825d 2672
2673 fhMCChHad1dR = new TH1F("hMCChHad1dR","TRACK matches dR, MC charged hadrons",ndRbins,dRmin,dRmax);
2674 fhMCChHad1dR->SetXTitle("#Delta R (rad)");
2675 outputContainer->Add(fhMCChHad1dR) ;
2676
2677 fhMCChHad2MatchdEdx = new TH2F("hMCChHad2MatchdEdx","dE/dx vs. p for all matches, MC charged hadrons",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
2678 fhMCChHad2MatchdEdx->SetXTitle("p (GeV/c)");
2679 fhMCChHad2MatchdEdx->SetYTitle("<dE/dx>");
2680 outputContainer->Add(fhMCChHad2MatchdEdx);
2681
d55bb5e1 2682 fhMCNeutral1EOverP = new TH2F("hMCNeutral1EOverP","TRACK matches E/p, MC neutrals",nptbins,ptmin,ptmax, nPoverEbins,EOverPmin,EOverPmax);
2683 fhMCNeutral1EOverP->SetYTitle("E/p");
2684 fhMCNeutral1EOverP->SetXTitle("p_{T} (GeV/c)");
2685 outputContainer->Add(fhMCNeutral1EOverP);
649b825d 2686
2687 fhMCNeutral1dR = new TH1F("hMCNeutral1dR","TRACK matches dR, MC neutrals",ndRbins,dRmin,dRmax);
2688 fhMCNeutral1dR->SetXTitle("#Delta R (rad)");
2689 outputContainer->Add(fhMCNeutral1dR) ;
2690
2691 fhMCNeutral2MatchdEdx = new TH2F("hMCNeutral2MatchdEdx","dE/dx vs. p for all matches, MC neutrals",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
2692 fhMCNeutral2MatchdEdx->SetXTitle("p (GeV/c)");
2693 fhMCNeutral2MatchdEdx->SetYTitle("<dE/dx>");
2694 outputContainer->Add(fhMCNeutral2MatchdEdx);
2695
d55bb5e1 2696 fhMCEle1EOverPR02 = new TH2F("hMCEle1EOverPR02","TRACK matches E/p, MC electrons",nptbins,ptmin,ptmax, nPoverEbins,EOverPmin,EOverPmax);
2697 fhMCEle1EOverPR02->SetYTitle("E/p");
2698 fhMCEle1EOverPR02->SetXTitle("p_{T} (GeV/c)");
2699 outputContainer->Add(fhMCEle1EOverPR02);
649b825d 2700
d55bb5e1 2701 fhMCChHad1EOverPR02 = new TH2F("hMCChHad1EOverPR02","TRACK matches E/p, MC charged hadrons",nptbins,ptmin,ptmax, nPoverEbins,EOverPmin,EOverPmax);
2702 fhMCChHad1EOverPR02->SetYTitle("E/p");
2703 fhMCChHad1EOverPR02->SetXTitle("p_{T} (GeV/c)");
2704 outputContainer->Add(fhMCChHad1EOverPR02);
f5036bcb 2705
d55bb5e1 2706 fhMCNeutral1EOverPR02 = new TH2F("hMCNeutral1EOverPR02","TRACK matches E/p, MC neutrals",nptbins,ptmin,ptmax, nPoverEbins,EOverPmin,EOverPmax);
2707 fhMCNeutral1EOverPR02->SetYTitle("E/p");
2708 fhMCNeutral1EOverPR02->SetXTitle("p_{T} (GeV/c)");
2709 outputContainer->Add(fhMCNeutral1EOverPR02);
521636d2 2710 }
f5036bcb 2711
649b825d 2712 // for(Int_t i = 0; i < outputContainer->GetEntries() ; i++)
2713 // printf("i=%d, name= %s\n",i,outputContainer->At(i)->GetName());
521636d2 2714
649b825d 2715 return outputContainer;
902aa95c 2716}
2717
f1538a5f 2718//__________________________________________________________________________________________________
2719Float_t AliAnaCalorimeterQA::GetECross(const Int_t absID, AliVCaloCells* cells, const Float_t dtcut)
1a72f6c5 2720{
2721 // Get energy in cross axis around maximum cell, for EMCAL only
2722
06f1b12a 2723 Int_t icol =-1, irow=-1,iRCU = -1;
2724 Int_t imod = GetModuleNumberCellIndexes(absID, fCalorimeter, icol, irow, iRCU);
e3300002 2725
57d8227a 2726 if(fCalorimeter=="EMCAL")
2727 {
06f1b12a 2728 //Get close cells index, energy and time, not in corners
e3300002 2729
2730 Int_t absID1 = -1;
2731 Int_t absID2 = -1;
2732
2733 if( irow < AliEMCALGeoParams::fgkEMCALRows-1) absID1 = GetCaloUtils()->GetEMCALGeometry()->GetAbsCellIdFromCellIndexes(imod, irow+1, icol);
2734 if( irow > 0 ) absID2 = GetCaloUtils()->GetEMCALGeometry()->GetAbsCellIdFromCellIndexes(imod, irow-1, icol);
2747966a 2735
2736 // In case of cell in eta = 0 border, depending on SM shift the cross cell index
2737 Int_t absID3 = -1;
2738 Int_t absID4 = -1;
2739
2740 if ( icol == AliEMCALGeoParams::fgkEMCALCols - 1 && !(imod%2) )
2741 {
e3300002 2742 absID3 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod+1, irow, 0);
2743 absID4 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod , irow, icol-1);
2747966a 2744 }
2745 else if( icol == 0 && imod%2 )
2746 {
e3300002 2747 absID3 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod , irow, icol+1);
2748 absID4 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod-1, irow, AliEMCALGeoParams::fgkEMCALCols-1);
2747966a 2749 }
2750 else
2751 {
e3300002 2752 if( icol < AliEMCALGeoParams::fgkEMCALCols-1 )
2753 absID3 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod, irow, icol+1);
2754 if( icol > 0 )
2755 absID4 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod, irow, icol-1);
2747966a 2756 }
06f1b12a 2757
2758 //Recalibrate cell energy if needed
2759 //Float_t ecell = cells->GetCellAmplitude(absID);
dbba06ca 2760 //GetCaloUtils()->RecalibrateCellAmplitude(ecell,fCalorimeter, absID);
06f1b12a 2761 Double_t tcell = cells->GetCellTime(absID);
dbba06ca 2762 GetCaloUtils()->RecalibrateCellTime(tcell, fCalorimeter, absID,GetReader()->GetInputEvent()->GetBunchCrossNumber());
06f1b12a 2763
2764 Float_t ecell1 = 0, ecell2 = 0, ecell3 = 0, ecell4 = 0;
2765 Double_t tcell1 = 0, tcell2 = 0, tcell3 = 0, tcell4 = 0;
2766
2747966a 2767 if(absID1 >0 )
2768 {
06f1b12a 2769 ecell1 = cells->GetCellAmplitude(absID1);
dbba06ca 2770 GetCaloUtils()->RecalibrateCellAmplitude(ecell1, fCalorimeter, absID1);
06f1b12a 2771 tcell1 = cells->GetCellTime(absID1);
dbba06ca 2772 GetCaloUtils()->RecalibrateCellTime (tcell1, fCalorimeter, absID1,GetReader()->GetInputEvent()->GetBunchCrossNumber());
06f1b12a 2773 }
2747966a 2774 if(absID2 >0 )
2775 {
06f1b12a 2776 ecell2 = cells->GetCellAmplitude(absID2);
dbba06ca 2777 GetCaloUtils()->RecalibrateCellAmplitude(ecell2, fCalorimeter, absID2);
06f1b12a 2778 tcell2 = cells->GetCellTime(absID2);
dbba06ca 2779 GetCaloUtils()->RecalibrateCellTime (tcell2, fCalorimeter, absID2, GetReader()->GetInputEvent()->GetBunchCrossNumber());
06f1b12a 2780 }
2747966a 2781 if(absID3 >0 )
2782 {
06f1b12a 2783 ecell3 = cells->GetCellAmplitude(absID3);
dbba06ca 2784 GetCaloUtils()->RecalibrateCellAmplitude(ecell3, fCalorimeter, absID3);
06f1b12a 2785 tcell3 = cells->GetCellTime(absID3);
dbba06ca 2786 GetCaloUtils()->RecalibrateCellTime (tcell3, fCalorimeter, absID3, GetReader()->GetInputEvent()->GetBunchCrossNumber());
06f1b12a 2787 }
2747966a 2788 if(absID4 >0 )
2789 {
06f1b12a 2790 ecell4 = cells->GetCellAmplitude(absID4);
dbba06ca 2791 GetCaloUtils()->RecalibrateCellAmplitude(ecell4, fCalorimeter, absID4);
06f1b12a 2792 tcell4 = cells->GetCellTime(absID4);
dbba06ca 2793 GetCaloUtils()->RecalibrateCellTime (tcell4, fCalorimeter, absID4, GetReader()->GetInputEvent()->GetBunchCrossNumber());
06f1b12a 2794 }
f1538a5f 2795
2796 if(TMath::Abs(tcell-tcell1)*1.e9 > dtcut) ecell1 = 0 ;
2797 if(TMath::Abs(tcell-tcell2)*1.e9 > dtcut) ecell2 = 0 ;
2798 if(TMath::Abs(tcell-tcell3)*1.e9 > dtcut) ecell3 = 0 ;
2799 if(TMath::Abs(tcell-tcell4)*1.e9 > dtcut) ecell4 = 0 ;
06f1b12a 2800
2801 return ecell1+ecell2+ecell3+ecell4;
1a72f6c5 2802 }
57d8227a 2803 else //PHOS
2804 {
06f1b12a 2805
2806 Int_t absId1 = -1, absId2 = -1, absId3 = -1, absId4 = -1;
2807
2808 Int_t relId1[] = { imod+1, 0, irow+1, icol };
2809 Int_t relId2[] = { imod+1, 0, irow-1, icol };
2810 Int_t relId3[] = { imod+1, 0, irow , icol+1 };
2811 Int_t relId4[] = { imod+1, 0, irow , icol-1 };
2812
2813 GetCaloUtils()->GetPHOSGeometry()->RelToAbsNumbering(relId1, absId1);
2814 GetCaloUtils()->GetPHOSGeometry()->RelToAbsNumbering(relId2, absId2);
2815 GetCaloUtils()->GetPHOSGeometry()->RelToAbsNumbering(relId3, absId3);
2816 GetCaloUtils()->GetPHOSGeometry()->RelToAbsNumbering(relId4, absId4);
2817
2818 Float_t ecell1 = 0, ecell2 = 0, ecell3 = 0, ecell4 = 0;
2819
2820 if(absId1 > 0 ) ecell1 = cells->GetCellAmplitude(absId1);
2821 if(absId2 > 0 ) ecell2 = cells->GetCellAmplitude(absId2);
2822 if(absId3 > 0 ) ecell3 = cells->GetCellAmplitude(absId3);
2823 if(absId4 > 0 ) ecell4 = cells->GetCellAmplitude(absId4);
2824
2825 return ecell1+ecell2+ecell3+ecell4;
2826
1a72f6c5 2827 }
2828
1a72f6c5 2829}
2830
c5693f62 2831//__________________________________________________________________________________________________
649b825d 2832void AliAnaCalorimeterQA::InvariantMassHistograms(const Int_t iclus, const TLorentzVector mom,
c5693f62 2833 const Int_t nModule, const TObjArray* caloClusters,
a82b4462 2834 AliVCaloCells * cells)
649b825d 2835{
2836 // Fill Invariant mass histograms
c8fe2783 2837
649b825d 2838 if(GetDebug()>1) printf("AliAnaCalorimeterQA::InvariantMassHistograms() - Start \n");
3748ffb5 2839
649b825d 2840 //Get vertex for photon momentum calculation and event selection
2841 Double_t v[3] = {0,0,0}; //vertex ;
1a83b960 2842 //GetReader()->GetVertex(v);
a6f26052 2843
649b825d 2844 Int_t nModule2 = -1;
2845 TLorentzVector mom2 ;
2846 Int_t nCaloClusters = caloClusters->GetEntriesFast();
a6f26052 2847
d07278cf 2848 for(Int_t jclus = iclus + 1 ; jclus < nCaloClusters ; jclus++)
2849 {
649b825d 2850 AliVCluster* clus2 = (AliVCluster*)caloClusters->At(jclus);
a82b4462 2851
2852 Float_t maxCellFraction = 0.;
2853 Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(cells, clus2,maxCellFraction);
55c05f8c 2854
d07278cf 2855 // Try to rediuce background with a mild shower shape cut and no more than 1 maxima
2856 // in cluster and remove low energy clusters
2857 if( clus2->GetNCells() <= 1 || !IsGoodCluster(absIdMax,cells) ||
2858 GetCaloUtils()->GetNumberOfLocalMaxima(clus2,cells) > 1 ||
2859 clus2->GetM02() > 0.5 || clus2->E() < 0.3) continue;
c8fe2783 2860
649b825d 2861 //Get cluster kinematics
2862 clus2->GetMomentum(mom2,v);
c8fe2783 2863
649b825d 2864 //Check only certain regions
2865 Bool_t in2 = kTRUE;
2866 if(IsFiducialCutOn()) in2 = GetFiducialCut()->IsInFiducialCut(mom2,fCalorimeter) ;
2867 if(!in2) continue;
2302a644 2868
649b825d 2869 //Get module of cluster
2870 nModule2 = GetModuleNumber(clus2);
c8fe2783 2871
649b825d 2872 //Fill histograms
c8fe2783 2873
649b825d 2874 //All modules
2875 fhIM ->Fill((mom+mom2).Pt(),(mom+mom2).M());
49214ef9 2876
649b825d 2877 //Single module
d07278cf 2878 if(nModule == nModule2 && nModule >=0 && nModule < fNModules)
649b825d 2879 fhIMMod[nModule]->Fill((mom+mom2).Pt(),(mom+mom2).M());
d07278cf 2880
c8fe2783 2881
649b825d 2882 //Asymetry histograms
2883 fhAsym->Fill((mom+mom2).Pt(),TMath::Abs((mom.E()-mom2.E())/(mom.E()+mom2.E())));
2884
2885 }// 2nd cluster loop
2886
2887}
2888
2889//______________________________
2890void AliAnaCalorimeterQA::Init()
2891{
2892 //Check if the data or settings are ok
c8fe2783 2893
649b825d 2894 if(fCalorimeter != "PHOS" && fCalorimeter !="EMCAL")
2895 AliFatal(Form("Wrong calorimeter name <%s>", fCalorimeter.Data()));
521636d2 2896
649b825d 2897 if(GetReader()->GetDataType()== AliCaloTrackReader::kMC)
2898 AliFatal("Analysis of reconstructed data, MC reader not aplicable");
2899
2900}
3bfc4732 2901
649b825d 2902//________________________________________
2903void AliAnaCalorimeterQA::InitParameters()
2904{
2905 //Initialize the parameters of the analysis.
2906 AddToHistogramsName("AnaCaloQA_");
2907
2908 fCalorimeter = "EMCAL"; //or PHOS
2909 fNModules = 12; // set maximum to maximum number of EMCAL modules
2910 fNRCU = 2; // set maximum number of RCU in EMCAL per SM
e6fec6f5 2911 fTimeCutMin = -9999999;
2912 fTimeCutMax = 9999999;
649b825d 2913 fEMCALCellAmpMin = 0.2;
2914 fPHOSCellAmpMin = 0.2;
c8fe2783 2915
f1538a5f 2916 // Exotic studies
2917 fExoNECrossCuts = 10 ;
2918 fExoNDTimeCuts = 4 ;
2919
2920 fExoDTimeCuts [0] = 1.e4 ; fExoDTimeCuts [1] = 50.0 ; fExoDTimeCuts [2] = 25.0 ; fExoDTimeCuts [3] = 10.0 ;
2921 fExoECrossCuts[0] = 0.80 ; fExoECrossCuts[1] = 0.85 ; fExoECrossCuts[2] = 0.90 ; fExoECrossCuts[3] = 0.92 ; fExoECrossCuts[4] = 0.94 ;
2922 fExoECrossCuts[5] = 0.95 ; fExoECrossCuts[6] = 0.96 ; fExoECrossCuts[7] = 0.97 ; fExoECrossCuts[8] = 0.98 ; fExoECrossCuts[9] = 0.99 ;
2923
649b825d 2924}
c8fe2783 2925
a82b4462 2926//___________________________________________________________________________________
2927Bool_t AliAnaCalorimeterQA::IsGoodCluster(const Int_t absIdMax, AliVCaloCells* cells)
649b825d 2928{
2929 //Identify cluster as exotic or not
2930
06f1b12a 2931 if(!fStudyBadClusters) return kTRUE;
a82b4462 2932
2747966a 2933 if(fCalorimeter=="EMCAL")
2934 {
06f1b12a 2935 if(!GetCaloUtils()->GetEMCALRecoUtils()->IsRejectExoticCluster())
2747966a 2936 {
06f1b12a 2937 return !( GetCaloUtils()->GetEMCALRecoUtils()->IsExoticCell(absIdMax,cells,(GetReader()->GetInputEvent())->GetBunchCrossNumber()) );
2747966a 2938 }
2939 else
2940 {
06f1b12a 2941 return kTRUE;
2747966a 2942 }
649b825d 2943 }
06f1b12a 2944 else // PHOS
2945 {
1a83b960 2946 Float_t ampMax = cells->GetCellAmplitude(absIdMax);
dbba06ca 2947 GetCaloUtils()->RecalibrateCellAmplitude(ampMax, fCalorimeter, absIdMax);
2747966a 2948
2949 if(ampMax < 0.01) return kFALSE;
2950
1a83b960 2951 if(1-GetECross(absIdMax,cells)/ampMax > 0.95) return kFALSE;
2952 else return kTRUE;
06f1b12a 2953 }
2954
649b825d 2955}
17708df9 2956
a82b4462 2957//_________________________________________________________
649b825d 2958void AliAnaCalorimeterQA::Print(const Option_t * opt) const
2959{
2960 //Print some relevant parameters set for the analysis
2961 if(! opt)
2962 return;
521636d2 2963
649b825d 2964 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
745913ae 2965 AliAnaCaloTrackCorrBaseClass::Print(" ");
2302a644 2966
649b825d 2967 printf("Select Calorimeter %s \n",fCalorimeter.Data());
2968 printf("Time Cut: %3.1f < TOF < %3.1f\n", fTimeCutMin, fTimeCutMax);
2969 printf("EMCAL Min Amplitude : %2.1f GeV/c\n", fEMCALCellAmpMin) ;
2970 printf("PHOS Min Amplitude : %2.1f GeV/c\n", fPHOSCellAmpMin) ;
c8fe2783 2971
649b825d 2972}
2973
649b825d 2974//_____________________________________________________
2975void AliAnaCalorimeterQA::MakeAnalysisFillHistograms()
2976{
2977 //Fill Calorimeter QA histograms
c8fe2783 2978
649b825d 2979 //Play with the MC stack if available
2980 if(IsDataMC()) MCHistograms();
798a9b04 2981
649b825d 2982 //Get List with CaloClusters
2983 TObjArray * caloClusters = NULL;
2984 if (fCalorimeter == "PHOS") caloClusters = GetPHOSClusters();
2985 else if (fCalorimeter == "EMCAL") caloClusters = GetEMCALClusters();
2986 else
2987 AliFatal(Form("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - Wrong calorimeter name <%s>, END\n", fCalorimeter.Data()));
798a9b04 2988
649b825d 2989 // Do not do anything if there are no clusters
2990 if(caloClusters->GetEntriesFast() == 0) return;
521636d2 2991
649b825d 2992 //Get List with CaloCells
2993 AliVCaloCells * cells = 0x0;
2994 if(fCalorimeter == "PHOS") cells = GetPHOSCells();
2995 else cells = GetEMCALCells();
798a9b04 2996
649b825d 2997 if(!caloClusters || !cells) {
2998 AliFatal(Form("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - No CaloClusters or CaloCells available\n"));
2999 return; // trick coverity
c8fe2783 3000 }
649b825d 3001
1a72f6c5 3002 //printf("QA: N cells %d, N clusters %d \n",cells->GetNumberOfCells(),caloClusters->GetEntriesFast());
3003
649b825d 3004 // Correlate Calorimeters and V0 and track Multiplicity
3005 if(fCorrelate) Correlate();
3006
3007 // Clusters
3008 ClusterLoopHistograms(caloClusters,cells);
3009
3010 // Cells
3011 CellHistograms(cells);
3012
3013 if(GetDebug() > 0)
3014 printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - End \n");
3015
a0bb4dc0 3016}
3017
649b825d 3018//______________________________________
3019void AliAnaCalorimeterQA::MCHistograms()
3020{
3021 //Get the MC arrays and do some checks before filling MC histograms
9e9f04cb 3022
649b825d 3023 TLorentzVector mom ;
3024
3025 if(GetReader()->ReadStack()){
9e9f04cb 3026
649b825d 3027 if(!GetMCStack())
3028 AliFatal("Stack not available, is the MC handler called?\n");
9e9f04cb 3029
649b825d 3030 //Fill some pure MC histograms, only primaries.
2747966a 3031 for(Int_t i=0 ; i<GetMCStack()->GetNprimary(); i++)
3032 {//Only primary particles, for all MC transport put GetNtrack()
649b825d 3033 TParticle *primary = GetMCStack()->Particle(i) ;
9e9f04cb 3034
649b825d 3035 if (primary->GetStatusCode() > 11) continue; //Working for PYTHIA and simple generators, check for HERWIG
3036 primary->Momentum(mom);
3037 MCHistograms(mom,TMath::Abs(primary->GetPdgCode()));
3038 } //primary loop
3039 }
3040 else if(GetReader()->ReadAODMCParticles()){
3041
3042 if(!GetReader()->GetAODMCParticles(0))
3043 AliFatal("AODMCParticles not available!");
3044
3045 //Fill some pure MC histograms, only primaries.
2747966a 3046 for(Int_t i=0 ; i < (GetReader()->GetAODMCParticles(0))->GetEntriesFast(); i++)
3047 {
649b825d 3048 AliAODMCParticle *aodprimary = (AliAODMCParticle*) (GetReader()->GetAODMCParticles(0))->At(i) ;
9e9f04cb 3049
649b825d 3050 if (!aodprimary->IsPrimary()) continue; //accept all which is not MC transport generated. Don't know how to avoid partons
9e9f04cb 3051
649b825d 3052 mom.SetPxPyPzE(aodprimary->Px(),aodprimary->Py(),aodprimary->Pz(),aodprimary->E());
3053 MCHistograms(mom,TMath::Abs(aodprimary->GetPdgCode()));
3054 } //primary loop
3055
3056 }
3057}
17708df9 3058
649b825d 3059//_______________________________________________________________________________
3060void AliAnaCalorimeterQA::MCHistograms(const TLorentzVector mom, const Int_t pdg)
3061{
a6f26052 3062 //Fill pure monte carlo related histograms
4865e325 3063
2302a644 3064 Float_t eMC = mom.E();
2302a644 3065 Float_t phiMC = mom.Phi();
3066 if(phiMC < 0)
3067 phiMC += TMath::TwoPi();
3068 Float_t etaMC = mom.Eta();
3069
3070 if (TMath::Abs(etaMC) > 1) return;
3071
35c71d5c 3072 Bool_t in = kFALSE;
3073
3074 //Rough stimate of acceptance for pi0, Eta and electrons
2747966a 3075 if(fCalorimeter == "PHOS")
3076 {
35c71d5c 3077 if(GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter))
3078 in = kTRUE ;
3079 if(GetDebug() > 2) printf("AliAnaCalorimeterQA::MCHistograms() - In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),in);
3080
3081 }
2747966a 3082 else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet())
3083 {
3084 if(GetEMCALGeometry())
3085 {
35c71d5c 3086 Int_t absID=0;
3087 GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(mom.Eta(),mom.Phi(),absID);
3088
3089 if( absID >= 0)
3090 in = kTRUE;
3091
3092 if(GetDebug() > 2) printf("AliAnaCalorimeterQA::MCHistograms() - In %s Real acceptance? %d\n",fCalorimeter.Data(),in);
3093 }
2747966a 3094 else
3095 {
35c71d5c 3096 if(GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter))
3097 in = kTRUE ;
3098 if(GetDebug() > 2) printf("AliAnaCalorimeterQA::MCHistograms() - In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),in);
3099 }
3100 }
2302a644 3101
2747966a 3102 if (pdg==22)
3103 {
c5693f62 3104 fhGenMCE[kmcPhoton] ->Fill(eMC);
3105 if(eMC > 0.5) fhGenMCEtaPhi[kmcPhoton]->Fill(etaMC,phiMC);
2747966a 3106 if(in)
3107 {
c5693f62 3108 fhGenMCAccE[kmcPhoton] ->Fill(eMC);
3109 if(eMC > 0.5) fhGenMCAccEtaPhi[kmcPhoton]->Fill(etaMC,phiMC);
2302a644 3110 }
3111 }
2747966a 3112 else if (pdg==111)
3113 {
c5693f62 3114 fhGenMCE[kmcPi0] ->Fill(eMC);
3115 if(eMC > 0.5) fhGenMCEtaPhi[kmcPi0]->Fill(etaMC,phiMC);
2747966a 3116 if(in)
3117 {
c5693f62 3118 fhGenMCAccE[kmcPi0] ->Fill(eMC);
3119 if(eMC > 0.5) fhGenMCAccEtaPhi[kmcPi0]->Fill(etaMC,phiMC);
2302a644 3120 }
3121 }
2747966a 3122 else if (pdg==221)
3123 {
c5693f62 3124 fhGenMCE[kmcEta] ->Fill(eMC);
3125 if(eMC > 0.5) fhGenMCEtaPhi[kmcEta]->Fill(etaMC,phiMC);
2747966a 3126 if(in)
3127 {
c5693f62 3128 fhGenMCAccE[kmcEta] ->Fill(eMC);
3129 if(eMC > 0.5) fhGenMCAccEtaPhi[kmcEta]->Fill(etaMC,phiMC);
35c71d5c 3130 }
2302a644 3131 }
2747966a 3132 else if (TMath::Abs(pdg)==11)
3133 {
c5693f62 3134 fhGenMCE[kmcElectron] ->Fill(eMC);
3135 if(eMC > 0.5) fhGenMCEtaPhi[kmcElectron]->Fill(etaMC,phiMC);
2747966a 3136 if(in)
3137 {
c5693f62 3138 fhGenMCAccE[kmcElectron] ->Fill(eMC);
3139 if(eMC > 0.5) fhGenMCAccEtaPhi[kmcElectron]->Fill(etaMC,phiMC);
35c71d5c 3140 }
2302a644 3141 }
902aa95c 3142}
c8fe2783 3143
649b825d 3144//_________________________________________________________________________________
3145void AliAnaCalorimeterQA::WeightHistograms(AliVCluster *clus, AliVCaloCells* cells)
3146{
3147 // Calculate weights
3148
3149 // First recalculate energy in case non linearity was applied
3150 Float_t energy = 0;
3151 Float_t ampMax = 0;
2747966a 3152 for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
3153 {
649b825d 3154 Int_t id = clus->GetCellsAbsId()[ipos];
3155
3156 //Recalibrate cell energy if needed
3157 Float_t amp = cells->GetCellAmplitude(id);
dbba06ca 3158 GetCaloUtils()->RecalibrateCellAmplitude(amp, fCalorimeter, id);
649b825d 3159
3160 energy += amp;
3161
3162 if(amp> ampMax)
3163 ampMax = amp;
3164
3165 } // energy loop
3166
2747966a 3167 if(energy <=0 )
3168 {
649b825d 3169 printf("AliAnaCalorimeterQA::WeightHistograms()- Wrong calculated energy %f\n",energy);
3170 return;
3171 }
3172
3173 fhEMaxCellClusterRatio ->Fill(energy,ampMax/energy);
3174 fhEMaxCellClusterLogRatio->Fill(energy,TMath::Log(ampMax/energy));
3175
3176 //Get the ratio and log ratio to all cells in cluster
2747966a 3177 for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
3178 {
649b825d 3179 Int_t id = clus->GetCellsAbsId()[ipos];
3180
3181 //Recalibrate cell energy if needed
3182 Float_t amp = cells->GetCellAmplitude(id);
dbba06ca 3183 GetCaloUtils()->RecalibrateCellAmplitude(amp, fCalorimeter, id);
649b825d 3184
3185 fhECellClusterRatio ->Fill(energy,amp/energy);
3186 fhECellClusterLogRatio->Fill(energy,TMath::Log(amp/energy));
3187 }
3188
3189 //Recalculate shower shape for different W0
2747966a 3190 if(fCalorimeter=="EMCAL")
3191 {
649b825d 3192 Float_t l0org = clus->GetM02();
3193 Float_t l1org = clus->GetM20();
3194 Float_t dorg = clus->GetDispersion();
3195
1a72f6c5 3196 for(Int_t iw = 0; iw < 14; iw++){
3197 GetCaloUtils()->GetEMCALRecoUtils()->SetW0(1+iw*0.5);
649b825d 3198 GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), cells, clus);
3199
3200 fhLambda0ForW0[iw]->Fill(energy,clus->GetM02());
1a72f6c5 3201 //fhLambda1ForW0[iw]->Fill(energy,clus->GetM20());
649b825d 3202
3203 if(IsDataMC()){
3204
3205 Int_t tag = GetMCAnalysisUtils()->CheckOrigin(clus->GetLabels(),clus->GetNLabels(), GetReader(),0);
3206
3207 if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton) &&
3208 !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) &&
3209 !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta) &&
3210 !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion) ){
3211 fhLambda0ForW0MC[iw][0]->Fill(energy,clus->GetM02());
1a72f6c5 3212 //fhLambda1ForW0MC[iw][0]->Fill(energy,clus->GetM20());
649b825d 3213 } // Pure Photon
3214 else if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCElectron) &&
3215 !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion) ){
3216 fhLambda0ForW0MC[iw][1]->Fill(energy,clus->GetM02());
1a72f6c5 3217 //fhLambda1ForW0MC[iw][1]->Fill(energy,clus->GetM20());
649b825d 3218 } // Electron
3219 else if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion) ){
3220 fhLambda0ForW0MC[iw][2]->Fill(energy,clus->GetM02());
1a72f6c5 3221 //fhLambda1ForW0MC[iw][2]->Fill(energy,clus->GetM20());
649b825d 3222 } // Conversion
3223 else if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) ){
3224 fhLambda0ForW0MC[iw][3]->Fill(energy,clus->GetM02());
1a72f6c5 3225 //fhLambda1ForW0MC[iw][3]->Fill(energy,clus->GetM20());
649b825d 3226 }// Pi0
3227 else if(!GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta) &&
3228 !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton) ){
3229 fhLambda0ForW0MC[iw][4]->Fill(energy,clus->GetM02());
1a72f6c5 3230 //fhLambda1ForW0MC[iw][4]->Fill(energy,clus->GetM20());
649b825d 3231 }// Hadron
3232
3233 }// Is MC
3234 } // w0 loop
3235
3236 // Set the original values back
3237 clus->SetM02(l0org);
3238 clus->SetM20(l1org);
3239 clus->SetDispersion(dorg);
3240
3241 }// EMCAL
3242
3243}
3244
3245
3246