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