1)AliCaloPID: Posibility to recalculate PID bayesian in EMCAL
[u/mrichter/AliRoot.git] / PWG4 / CaloCalib / AliAnalysisTaskEMCALPi0CalibSelection.cxx
CommitLineData
375cec9b 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: Boris Polishchuk *
5 * Adapted to AOD reading by Gustavo Conesa *
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 **************************************************************************/
15
16//---------------------------------------------------------------------------//
17// //
18// Fill histograms (one per cell) with two-cluster invariant mass //
19// using calibration coefficients of the previous iteration. //
20// Histogram for a given cell is filled if the most energy of one cluster //
21// is deposited in this cell and the other cluster could be anywherein EMCAL.//
22// //
23//---------------------------------------------------------------------------//
24
25//#include <cstdlib>
26//#include <Riostream.h>
27// Root
28#include "TLorentzVector.h"
29//#include "TVector3.h"
30#include "TRefArray.h"
31#include "TList.h"
32#include "TH1F.h"
33
34// AliRoot
35#include "AliAnalysisTaskEMCALPi0CalibSelection.h"
36#include "AliAODEvent.h"
37#include "AliESDEvent.h"
38#include "AliESDCaloCluster.h"
39#include "AliESDCaloCells.h"
40#include "AliEMCALGeometry.h"
41#include "AliAODCaloCluster.h"
42#include "AliAODCaloCells.h"
43#include "AliEMCALAodCluster.h"
44#include "AliEMCALCalibData.h"
45
46ClassImp(AliAnalysisTaskEMCALPi0CalibSelection)
47
48//__________________________________________________
49AliAnalysisTaskEMCALPi0CalibSelection::AliAnalysisTaskEMCALPi0CalibSelection() :
50AliAnalysisTaskSE(),fEMCALGeo(0x0),fCalibData(0x0), fEmin(0.), fLogWeight(4.5), fCopyAOD(kFALSE),
51fEMCALGeoName("EMCAL_COMPLETE"), fOldData(kFALSE),fOutputContainer(0x0),fHmgg(0x0)
52{
53 //Default constructor.
54
bdd2a262 55 for(Int_t iMod=0; iMod < 12; iMod++) {
56 for(Int_t iX=0; iX<24; iX++) {
57 for(Int_t iZ=0; iZ<48; iZ++) {
375cec9b 58 fHmpi0[iMod][iX][iZ]=0;
59 }
60 }
61 }
62
63}
64
65//__________________________________________________
66AliAnalysisTaskEMCALPi0CalibSelection::AliAnalysisTaskEMCALPi0CalibSelection(const char* name) :
67 AliAnalysisTaskSE(name),fEMCALGeo(0x0),fCalibData(0x0), fEmin(0.), fLogWeight(4.5), fCopyAOD(kFALSE),
68 fEMCALGeoName("EMCAL_COMPLETE"), fOldData(kFALSE), fOutputContainer(0x0),fHmgg(0x0)
69{
70 //Named constructor which should be used.
71
72 DefineOutput(1,TList::Class());
73
bdd2a262 74 for(Int_t iMod=0; iMod < 12; iMod++) {
75 for(Int_t iX=0; iX<24; iX++) {
76 for(Int_t iZ=0; iZ<48; iZ++) {
375cec9b 77 fHmpi0[iMod][iX][iZ]=0;
78 }
79 }
80 }
81
82}
83
84//__________________________________________________
85AliAnalysisTaskEMCALPi0CalibSelection::~AliAnalysisTaskEMCALPi0CalibSelection()
86{
87 //Destructor.
88
89 if(fOutputContainer){
90 fOutputContainer->Delete() ;
91 delete fOutputContainer ;
92 }
93
94 if(fCalibData) delete fCalibData;
95 if(fEMCALGeo) delete fEMCALGeo;
96
97}
98
99
100//__________________________________________________
101void AliAnalysisTaskEMCALPi0CalibSelection::CreateAODFromAOD()
102{
103 // Copy AOD header, vertex, CaloClusters and CaloCells to output AOD
104
105 AliAODEvent* aod = dynamic_cast<AliAODEvent*>(InputEvent());
106
107 // set arrays and pointers
108 Float_t posF[3];
109 Double_t pos[3];
110
111 Double_t covVtx[6];
112
113 for (Int_t i = 0; i < 6; i++) covVtx[i] = 0.;
114
115 // Update the header
116 AliAODHeader*headerin = aod->GetHeader();
117 AliAODHeader* header = AODEvent()->GetHeader();
118 header->SetRunNumber(headerin->GetRunNumber());
119 header->SetBunchCrossNumber(headerin->GetBunchCrossNumber());
120 header->SetOrbitNumber(headerin->GetOrbitNumber());
121 header->SetPeriodNumber(headerin->GetPeriodNumber());
122 header->SetEventType(headerin->GetEventType());
123 header->SetMuonMagFieldScale(headerin->GetMuonMagFieldScale());
124 header->SetCentrality(headerin->GetCentrality());
125
126 header->SetTriggerMask(headerin->GetTriggerMask());
127 header->SetTriggerCluster(headerin->GetTriggerCluster());
128 header->SetMagneticField(headerin->GetMagneticField());
129 header->SetZDCN1Energy(headerin->GetZDCN1Energy());
130 header->SetZDCP1Energy(headerin->GetZDCP1Energy());
131 header->SetZDCN2Energy(headerin->GetZDCN2Energy());
132 header->SetZDCP2Energy(headerin->GetZDCP2Energy());
133 header->SetZDCEMEnergy(headerin->GetZDCEMEnergy(0),headerin->GetZDCEMEnergy(1));
134 Float_t diamxy[2]={aod->GetDiamondX(),aod->GetDiamondY()};
135 Float_t diamcov[3]; aod->GetDiamondCovXY(diamcov);
136 header->SetDiamond(diamxy,diamcov);
137 //
138 //
139 Int_t nVertices = 1 ;/* = prim. vtx*/;
140 Int_t nCaloClus = aod->GetNCaloClusters();
141
142 AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
143
144 // Access to the AOD container of vertices
145 TClonesArray &vertices = *(AODEvent()->GetVertices());
146 Int_t jVertices=0;
147
148 // Add primary vertex. The primary tracks will be defined
149 // after the loops on the composite objects (V0, cascades, kinks)
150 const AliAODVertex *vtx = aod->GetPrimaryVertex();
151
152 vtx->GetXYZ(pos); // position
153 vtx->GetCovMatrix(covVtx); //covariance matrix
154
155 AliAODVertex * primary = new(vertices[jVertices++])
156 AliAODVertex(pos, covVtx, vtx->GetChi2perNDF(), NULL, -1, AliAODVertex::kPrimary);
157 primary->SetName(vtx->GetName());
158 primary->SetTitle(vtx->GetTitle());
159
160 // Access to the AOD container of clusters
161 TClonesArray &caloClusters = *(AODEvent()->GetCaloClusters());
162 Int_t jClusters=0;
163
164 for (Int_t iClust=0; iClust<nCaloClus; ++iClust) {
165
166 AliAODCaloCluster * cluster = aod->GetCaloCluster(iClust);
167
168 //Check if it is a EMCAL cluster
169 if(!cluster->IsEMCALCluster()) continue ;
170
171 Int_t id = cluster->GetID();
172 Float_t energy = cluster->E();
173 cluster->GetPosition(posF);
174 Char_t ttype = cluster->GetType();
175
176 AliAODCaloCluster *caloCluster = new(caloClusters[jClusters++])
177 AliAODCaloCluster(id,
178 0,
179 0x0,
180 energy,
181 posF,
182 NULL,
183 ttype);
184
185 caloCluster->SetCaloCluster(cluster->GetDistToBadChannel(),
186 cluster->GetDispersion(),
187 cluster->GetM20(), cluster->GetM02(),
188 cluster->GetEmcCpvDistance(),
189 cluster->GetNExMax(),cluster->GetTOF()) ;
190
191 caloCluster->SetPIDFromESD(cluster->PID());
192 caloCluster->SetNCells(cluster->GetNCells());
193 caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
194 if(!fOldData){
195 caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
196 }
197 else{
198 Double_t *fraction = cluster->GetCellsAmplitudeFraction();
199 for(Int_t i = 0; i < cluster->GetNCells() ; i++) fraction[i] = 1;
200 caloCluster->SetCellsAmplitudeFraction(fraction);
201 }
202
203 }
204
205 caloClusters.Expand(jClusters); // resize TObjArray to 'remove' slots for pseudo clusters
206 // end of loop on calo clusters
207
208 // fill EMCAL cell info
209 if (aod->GetEMCALCells()) { // protection against missing AOD information
210 AliAODCaloCells &aodinEMcells = *(aod->GetEMCALCells());
211 Int_t nEMcell = aodinEMcells.GetNumberOfCells() ;
212
213 AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
214 aodEMcells.CreateContainer(nEMcell);
215 aodEMcells.SetType(AliAODCaloCells::kEMCAL);
216
217 Double_t calibFactor = 1;
218 if(fOldData) calibFactor = 0.0153;
219
220 for (Int_t iCell = 0; iCell < nEMcell; iCell++) {
221 aodEMcells.SetCell(iCell,aodinEMcells.GetCellNumber(iCell),aodinEMcells.GetAmplitude(iCell)*calibFactor);
222 }
223 aodEMcells.Sort();
224
225 }
226
227}
228
229//__________________________________________________
230void AliAnalysisTaskEMCALPi0CalibSelection::CreateAODFromESD()
231{
232
233 // Copy Header, Vertex, CaloClusters and CaloCells from ESDs to AODs
234
235 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
236
237 // set arrays and pointers
238 Float_t posF[3];
239 Double_t pos[3];
240
241 Double_t covVtx[6];
242
243 for (Int_t i = 0; i < 6; i++) covVtx[i] = 0.;
244
245 // Update the header
246
247 AliAODHeader* header = AODEvent()->GetHeader();
248 header->SetRunNumber(esd->GetRunNumber());
249 header->SetBunchCrossNumber(esd->GetBunchCrossNumber());
250 header->SetOrbitNumber(esd->GetOrbitNumber());
251 header->SetPeriodNumber(esd->GetPeriodNumber());
252 header->SetEventType(esd->GetEventType());
253 header->SetMuonMagFieldScale(-999.); // FIXME
254 header->SetCentrality(-999.); // FIXME
255
256
257 header->SetTriggerMask(esd->GetTriggerMask());
258 header->SetTriggerCluster(esd->GetTriggerCluster());
259 header->SetMagneticField(esd->GetMagneticField());
260 header->SetZDCN1Energy(esd->GetZDCN1Energy());
261 header->SetZDCP1Energy(esd->GetZDCP1Energy());
262 header->SetZDCN2Energy(esd->GetZDCN2Energy());
263 header->SetZDCP2Energy(esd->GetZDCP2Energy());
264 header->SetZDCEMEnergy(esd->GetZDCEMEnergy(0),esd->GetZDCEMEnergy(1));
265 Float_t diamxy[2]={esd->GetDiamondX(),esd->GetDiamondY()};
266 Float_t diamcov[3]; esd->GetDiamondCovXY(diamcov);
267 header->SetDiamond(diamxy,diamcov);
268 //
269 //
270 Int_t nVertices = 1 ;/* = prim. vtx*/;
271 Int_t nCaloClus = esd->GetNumberOfCaloClusters();
272
273 AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
274
275 // Access to the AOD container of vertices
276 TClonesArray &vertices = *(AODEvent()->GetVertices());
277 Int_t jVertices=0;
278
279 // Add primary vertex. The primary tracks will be defined
280 // after the loops on the composite objects (V0, cascades, kinks)
281 const AliESDVertex *vtx = esd->GetPrimaryVertex();
282
283 vtx->GetXYZ(pos); // position
284 vtx->GetCovMatrix(covVtx); //covariance matrix
285
286 AliAODVertex * primary = new(vertices[jVertices++])
287 AliAODVertex(pos, covVtx, vtx->GetChi2toNDF(), NULL, -1, AliAODVertex::kPrimary);
288 primary->SetName(vtx->GetName());
289 primary->SetTitle(vtx->GetTitle());
290
291 // Access to the AOD container of clusters
292 TClonesArray &caloClusters = *(AODEvent()->GetCaloClusters());
293 Int_t jClusters=0;
294
295 for (Int_t iClust=0; iClust<nCaloClus; ++iClust) {
296
297 AliESDCaloCluster * cluster = esd->GetCaloCluster(iClust);
298
299 //Check which calorimeter information we want to keep.
300 if(!cluster->IsEMCAL()) continue ;
301
302 Int_t id = cluster->GetID();
303 Float_t energy = cluster->E();
304 cluster->GetPosition(posF);
305
306 AliAODCaloCluster *caloCluster = new(caloClusters[jClusters++])
307 AliAODCaloCluster(id,
308 0,
309 0x0,
310 energy,
311 posF,
312 NULL,
313 AliAODCluster::kEMCALClusterv1);
314
315 caloCluster->SetCaloCluster(cluster->GetDistanceToBadChannel(),
316 cluster->GetClusterDisp(),
317 cluster->GetM20(), cluster->GetM02(),
318 cluster->GetEmcCpvDistance(),
319 cluster->GetNExMax(),cluster->GetTOF()) ;
320
321 caloCluster->SetPIDFromESD(cluster->GetPid());
322 caloCluster->SetNCells(cluster->GetNCells());
323 caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
324 if(!fOldData){
325 caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
326 }
327 else{
328 Double_t *fraction = cluster->GetCellsAmplitudeFraction();
329 for(Int_t i = 0; i < cluster->GetNCells() ; i++) fraction[i] = 1;
330 caloCluster->SetCellsAmplitudeFraction(fraction);
331 }
332 }
333 caloClusters.Expand(jClusters); // resize TObjArray to 'remove' slots for pseudo clusters
334 // end of loop on calo clusters
335
336 // fill EMCAL cell info
337
338 if( esd->GetEMCALCells()) { // protection against missing ESD information
339 AliESDCaloCells &esdEMcells = *(esd->GetEMCALCells());
340 Int_t nEMcell = esdEMcells.GetNumberOfCells() ;
341
342 AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
343 aodEMcells.CreateContainer(nEMcell);
344 aodEMcells.SetType(AliAODCaloCells::kEMCAL);
345
346 Double_t calibFactor = 1;
347 if(fOldData) calibFactor = 0.0153;
348
349 for (Int_t iCell = 0; iCell < nEMcell; iCell++) {
350 aodEMcells.SetCell(iCell,esdEMcells.GetCellNumber(iCell),esdEMcells.GetAmplitude(iCell)*calibFactor);
351 }
352 aodEMcells.Sort();
353
354 }
355
356}
357
358//__________________________________________________
359void AliAnalysisTaskEMCALPi0CalibSelection::UserCreateOutputObjects()
360{
361 //Create output container
362 fOutputContainer = new TList();
363
364 char hname[128], htitl[128];
365
366 for(Int_t iMod=0; iMod < AliEMCALGeoParams::fgkEMCALModules; iMod++) {
367 for(Int_t iX=0; iX<AliEMCALGeoParams::fgkEMCALCols; iX++) {
368 for(Int_t iZ=0; iZ<AliEMCALGeoParams::fgkEMCALRows; iZ++) {
369 sprintf(hname,"%d_%d_%d",iMod,iX,iZ);
370 sprintf(htitl,"Two-gamma inv. mass for mod %d, cell (%d,%d)",iMod,iX,iZ);
371 fHmpi0[iMod][iX][iZ] = new TH1F(hname,htitl,100,0.,300.);
372 fOutputContainer->Add(fHmpi0[iMod][iX][iZ]);
373 }
374 }
375 }
376
377 fHmgg = new TH1F("hmgg","2-cluster invariant mass",100,0.,300.);
378 fOutputContainer->Add(fHmgg);
379
380 fCalibData = new AliEMCALCalibData();
381
382 if(!fCalibData)
383 AliFatal("Calibration parameters not found in CDB!");
384
385 printf("Get Geometry : %s\n", fEMCALGeoName.Data());
386 fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName) ;
387
388}
389
390//__________________________________________________
391void AliAnalysisTaskEMCALPi0CalibSelection::UserExec(Option_t* /* option */)
392{
393 //Analysis per event.
394 if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection <<< Event %d >>>\n",(Int_t)Entry());
395
396 AliAODEvent* aod = 0x0;
397 if(!strcmp(InputEvent()->GetName(),"AliAODEvent")) {
398 //Input are ESDs
399 aod = dynamic_cast<AliAODEvent*>(InputEvent());
400 // Create new AOD with only CaloClusters and CaloCells
401 if(fCopyAOD) CreateAODFromAOD();
402 }
403 else if(!strcmp(InputEvent()->GetName(),"AliESDEvent")) {
404 //Input are ESDs
405 aod = AODEvent();
406 // Create AOD with CaloClusters and use it as input.
407 // If filtering task is already executed, this is not needed.
408 if(fCopyAOD) CreateAODFromESD();
409 }
410 else {
411 printf("AliAnalysisTaskEMCALPi0CalibSelection: Unknown event type, STOP!\n");
412 abort();
413 }
414
415 Double_t v[] = {aod->GetVertex(0)->GetX(),aod->GetVertex(0)->GetY(),aod->GetVertex(0)->GetZ()}; //to check!!
416 //aod->GetVertex()->GetXYZ(v) ;
417 //TVector3 vtx(v);
418
419 //if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Vertex: (%.3f,%.3f,%.3f)\n",vtx.X(),vtx.Y(),vtx.Z());
420 if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Vertex: (%.3f,%.3f,%.3f)\n",v[0],v[1],v[2]);
421
422 Int_t runNum = aod->GetRunNumber();
423 if(DebugLevel() > 1) printf("Run number: %d\n",runNum);
424
425 //Get the matrix with geometry information
426 //Still not implemented in AOD, just a workaround to be able to work at least with ESDs
427 if(!strcmp(InputEvent()->GetName(),"AliAODEvent")) {
428 if(DebugLevel() > 1)
429 printf("AliAnalysisTaskEMCALPi0CalibSelection Use ideal geometry, values geometry matrix not kept in AODs.\n");
430 }
431 else{
432 if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Load Misaligned matrices. \n");
433 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent()) ;
bdd2a262 434 for(Int_t mod=0; mod < 12; mod++){
375cec9b 435 if(esd->GetEMCALMatrix(mod)) fEMCALGeo->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
436 }
437 }
438
439 if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Will use fLogWeight %.3f .\n",fLogWeight);
440
441 //AliEMCALPID pid;
442
443 Int_t maxId;
444 Int_t iSupMod = -1;
445 Int_t iTower = -1;
446 Int_t iIphi = -1;
447 Int_t iIeta = -1;
448 Int_t iphi = -1;
449 Int_t ieta = -1;
450
451 TLorentzVector p1;
452 TLorentzVector p2;
453 TLorentzVector p12;
454
455
456 TRefArray * caloClustersArr = new TRefArray();
457 aod->GetEMCALClusters(caloClustersArr);
458
459 const Int_t kNumberOfEMCALClusters = caloClustersArr->GetEntries() ;
460 if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection CaloClusters: %d\n", kNumberOfEMCALClusters);
461
462 // EMCAL cells
463 AliAODCaloCells *emCells = aod->GetEMCALCells();
464
465 // Check if is old data not modified for calibration and change the fraction factor and calibrate amplitudes
466 // Only for aliroot older than release 17 tag 3
467 if(fOldData && !fCopyAOD){
468 //Calibrate amplitudes
469 AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
470 for (Int_t iCell = 0; iCell < aodEMcells.GetNumberOfCells(); iCell++) {
471 aodEMcells.SetCell(iCell,aodEMcells.GetCellNumber(iCell),aodEMcells.GetAmplitude(iCell)*0.0153);
472 }
473
474 //Change fraction
475 for(Int_t iClu=0; iClu<kNumberOfEMCALClusters; iClu++) {
476 AliAODCaloCluster *cl = (AliAODCaloCluster *) caloClustersArr->At(iClu);
477 if(!cl->IsEMCALCluster()) continue; // EMCAL cluster!
478 Double_t *fraction = cl->GetCellsAmplitudeFraction();
479 for(Int_t i = 0; i < cl->GetNCells() ; i++) fraction[i] = 1;
480 cl->SetCellsAmplitudeFraction(fraction);
481 }
482 }
483
484 // loop over EMCAL clusters
485 for(Int_t iClu=0; iClu<kNumberOfEMCALClusters; iClu++) {
486
487 AliAODCaloCluster *c1 = (AliAODCaloCluster *) caloClustersArr->At(iClu);
488 if(!c1->IsEMCALCluster()) continue; // EMCAL cluster!
489
490 Float_t e1i = c1->E(); // cluster energy before correction
491 if(e1i<fEmin) continue;
492
493 if(DebugLevel() > 2){
494 printf("Std : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",iClu,e1i, c1->GetDispersion(),c1->GetM02(),c1->GetM20());
495 Double_t pos[]={0,0,0};
496 c1->GetPosition(pos);
497 printf("Std : i %d, x %f, y %f, z %f\n",iClu, pos[0], pos[1], pos[2]);
498 }
499
500 AliEMCALAodCluster clu1(*c1);
501 clu1.Recalibrate(fCalibData, emCells, fEMCALGeoName);
502 printf("recalibrated, now calculate the rest\n");
503 clu1.EvalAll(fLogWeight, fEMCALGeoName);
504 //clu1.EnergyCorrection(&pid) ;
505
506 if(DebugLevel() > 2){
507 printf("Recal: i %d, E %f, dispersion %f, m02 %f, m20 %f\n",iClu,clu1.E(), clu1.GetDispersion(),clu1.GetM02(),clu1.GetM20());
508 Double_t pos2[]={0,0,0};
509 clu1.GetPosition(pos2);
510 printf("Recal: i %d, x %f, y %f, z %f\n",iClu, pos2[0], pos2[1], pos2[2]);
511 }
512
513 clu1.GetMomentum(p1,v);
514
515 MaxEnergyCellPos(emCells,&clu1,maxId);
516
517 //Get from the absid the supermodule, tower and eta/phi numbers
518 fEMCALGeo->GetCellIndex(maxId,iSupMod,iTower,iIphi,iIeta);
519 //Gives SuperModule and Tower numbers
520 fEMCALGeo->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
521 iIphi, iIeta,iphi,ieta);
522
523 Float_t e1ii = clu1.E(); // cluster energy after correction
524
525 for (Int_t jClu=iClu; jClu<kNumberOfEMCALClusters; jClu++) {
526 AliAODCaloCluster *c2 = (AliAODCaloCluster *) caloClustersArr->At(jClu);
527 if(!c2->IsEMCALCluster()) continue; // EMCAL cluster!
528 if(c2->IsEqual(c1)) continue;
529
530 Float_t e2i = c2->E();
531 if(e2i<fEmin) continue;
532
533 AliEMCALAodCluster clu2(*c2);
534 //printf("i2 %d, E %f\n",iClu,e2i);
535 clu2.Recalibrate(fCalibData, emCells,fEMCALGeoName);
536 clu2.EvalAll(fLogWeight,fEMCALGeoName);
537 //clu2.EnergyCorrection(&pid) ;
538 // printf("i2 %d, Erc %f\n",iClu,clu2.E());
539
540 clu2.GetMomentum(p2,v);
541 Float_t e2ii = clu2.E();
542
543 p12 = p1+p2;
544 Float_t invmass = p12.M()*1000;
545
546 fHmgg->Fill(invmass);
547
548 //printf("iSM %d, ieta %d, iphi %d, mass %f \n",iSupMod, ieta, iphi, invmass);
549 if(invmass < 300)
550 fHmpi0[iSupMod][ieta][iphi]->Fill(invmass);
551
552
553 if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Mass in (SM%d,%d,%d): %.3f GeV E1_i=%f E1_ii=%f E2_i=%f E2_ii=%f\n",
554 iSupMod,iphi,ieta,p12.M(),e1i,e1ii,e2i,e2ii);
555 }
556
557 } // end of loop over EMCAL clusters
558
559 delete caloClustersArr;
560 PostData(1,fOutputContainer);
561}
562
563//__________________________________________________
bdd2a262 564void AliAnalysisTaskEMCALPi0CalibSelection::MaxEnergyCellPos(AliAODCaloCells* const cells, AliAODCaloCluster* const clu, Int_t& maxId)
375cec9b 565{
566 //For a given CaloCluster calculates the absId of the cell
567 //with maximum energy deposit.
568
569 Double_t eMax = -111;
570
571 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
572 Int_t cellAbsId = clu->GetCellAbsId(iDig);
573 Double_t eCell = cells->GetCellAmplitude(cellAbsId)*clu->GetCellAmplitudeFraction(iDig);
574 if(eCell>eMax) {
575 eMax = eCell;
576 maxId = cellAbsId;
577 }
578 }
579
580}
581
582//__________________________________________________
bdd2a262 583void AliAnalysisTaskEMCALPi0CalibSelection::SetCalibCorrections(AliEMCALCalibData* const cdata)
375cec9b 584{
585 //Set new correction factors (~1) to calibration coefficients, delete previous.
586
587 if(fCalibData) delete fCalibData;
588 fCalibData = cdata;
589
590}
591