]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/CaloCalib/AliAnalysisTaskEMCALPi0CalibSelection.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Many improvements in EMCAL calibration class:
[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"
6eb2a715 43//#include "AliEMCALAodCluster.h"
44//#include "AliEMCALCalibData.h"
375cec9b 45
46ClassImp(AliAnalysisTaskEMCALPi0CalibSelection)
47
375cec9b 48
49//__________________________________________________
50AliAnalysisTaskEMCALPi0CalibSelection::AliAnalysisTaskEMCALPi0CalibSelection(const char* name) :
6eb2a715 51 AliAnalysisTaskSE(name),fEMCALGeo(0x0),//fCalibData(0x0),
52 fEmin(0.5), fEmax(15.), fMinNCells(2), fGroupNCells(0),
53 fLogWeight(4.5), fCopyAOD(kFALSE), fEMCALGeoName("EMCAL_FIRSTYEAR"),
54 fRemoveBadChannels(kFALSE),fEMCALBadChannelMap(0x0),
55 fRecalibration(kFALSE),fEMCALRecalibrationFactors(),
56 fNbins(300), fMinBin(0.), fMaxBin(300.),
57 fOutputContainer(0x0),fHmgg(0x0), fhNEvents(0x0),fCuts(0x0)
375cec9b 58{
59 //Named constructor which should be used.
60
bdd2a262 61 for(Int_t iMod=0; iMod < 12; iMod++) {
62 for(Int_t iX=0; iX<24; iX++) {
63 for(Int_t iZ=0; iZ<48; iZ++) {
375cec9b 64 fHmpi0[iMod][iX][iZ]=0;
65 }
66 }
67 }
6eb2a715 68
cf028690 69 DefineOutput(1, TList::Class());
6eb2a715 70 DefineOutput(2, TList::Class()); // will contain cuts or local params
375cec9b 71
72}
73
74//__________________________________________________
75AliAnalysisTaskEMCALPi0CalibSelection::~AliAnalysisTaskEMCALPi0CalibSelection()
76{
77 //Destructor.
78
79 if(fOutputContainer){
80 fOutputContainer->Delete() ;
81 delete fOutputContainer ;
82 }
83
6eb2a715 84 //if(fCalibData) delete fCalibData;
375cec9b 85 if(fEMCALGeo) delete fEMCALGeo;
86
6eb2a715 87 if(fEMCALBadChannelMap) {
88 fEMCALBadChannelMap->Clear();
89 delete fEMCALBadChannelMap;
90 }
91
92
93 if(fEMCALRecalibrationFactors) {
94 fEMCALRecalibrationFactors->Clear();
95 delete fEMCALRecalibrationFactors;
96 }
375cec9b 97}
98
6eb2a715 99//_____________________________________________________
100void AliAnalysisTaskEMCALPi0CalibSelection::LocalInit()
101{
102 // Local Initialization
103
104 // Create cuts/param objects and publish to slot
105
106 char onePar[255] ;
107 fCuts = new TList();
108
109 sprintf(onePar,"Custer cuts: %2.2f < E < %2.2f GeV; min number of cells %d;", fEmin,fEmax, fMinNCells) ;
110 fCuts->Add(new TObjString(onePar));
111 sprintf(onePar,"Group %d cells;", fGroupNCells) ;
112 fCuts->Add(new TObjString(onePar));
113 sprintf(onePar,"Histograms: bins %d; energy range: %2.2f < E < %2.2f GeV;",fNbins,fMinBin,fMaxBin) ;
114 fCuts->Add(new TObjString(onePar));
115 sprintf(onePar,"Switchs: Remove Bad Channels? %d; Copy AODs? %d; Recalibrate %d?",fRemoveBadChannels,fCopyAOD,fRecalibration) ;
116 fCuts->Add(new TObjString(onePar));
117 sprintf(onePar,"EMCAL Geometry name: < %s >",fEMCALGeoName.Data()) ;
118 fCuts->Add(new TObjString(onePar));
119
120 // Post Data
121 PostData(2, fCuts);
122
123}
375cec9b 124
125//__________________________________________________
126void AliAnalysisTaskEMCALPi0CalibSelection::CreateAODFromAOD()
127{
128 // Copy AOD header, vertex, CaloClusters and CaloCells to output AOD
129
130 AliAODEvent* aod = dynamic_cast<AliAODEvent*>(InputEvent());
131
132 // set arrays and pointers
133 Float_t posF[3];
134 Double_t pos[3];
135
136 Double_t covVtx[6];
137
138 for (Int_t i = 0; i < 6; i++) covVtx[i] = 0.;
139
140 // Update the header
141 AliAODHeader*headerin = aod->GetHeader();
142 AliAODHeader* header = AODEvent()->GetHeader();
143 header->SetRunNumber(headerin->GetRunNumber());
144 header->SetBunchCrossNumber(headerin->GetBunchCrossNumber());
145 header->SetOrbitNumber(headerin->GetOrbitNumber());
146 header->SetPeriodNumber(headerin->GetPeriodNumber());
147 header->SetEventType(headerin->GetEventType());
148 header->SetMuonMagFieldScale(headerin->GetMuonMagFieldScale());
149 header->SetCentrality(headerin->GetCentrality());
150
151 header->SetTriggerMask(headerin->GetTriggerMask());
152 header->SetTriggerCluster(headerin->GetTriggerCluster());
153 header->SetMagneticField(headerin->GetMagneticField());
154 header->SetZDCN1Energy(headerin->GetZDCN1Energy());
155 header->SetZDCP1Energy(headerin->GetZDCP1Energy());
156 header->SetZDCN2Energy(headerin->GetZDCN2Energy());
157 header->SetZDCP2Energy(headerin->GetZDCP2Energy());
158 header->SetZDCEMEnergy(headerin->GetZDCEMEnergy(0),headerin->GetZDCEMEnergy(1));
159 Float_t diamxy[2]={aod->GetDiamondX(),aod->GetDiamondY()};
160 Float_t diamcov[3]; aod->GetDiamondCovXY(diamcov);
161 header->SetDiamond(diamxy,diamcov);
162 //
163 //
164 Int_t nVertices = 1 ;/* = prim. vtx*/;
165 Int_t nCaloClus = aod->GetNCaloClusters();
166
167 AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
168
169 // Access to the AOD container of vertices
170 TClonesArray &vertices = *(AODEvent()->GetVertices());
171 Int_t jVertices=0;
172
173 // Add primary vertex. The primary tracks will be defined
174 // after the loops on the composite objects (V0, cascades, kinks)
175 const AliAODVertex *vtx = aod->GetPrimaryVertex();
176
177 vtx->GetXYZ(pos); // position
178 vtx->GetCovMatrix(covVtx); //covariance matrix
179
180 AliAODVertex * primary = new(vertices[jVertices++])
181 AliAODVertex(pos, covVtx, vtx->GetChi2perNDF(), NULL, -1, AliAODVertex::kPrimary);
182 primary->SetName(vtx->GetName());
183 primary->SetTitle(vtx->GetTitle());
184
185 // Access to the AOD container of clusters
186 TClonesArray &caloClusters = *(AODEvent()->GetCaloClusters());
187 Int_t jClusters=0;
188
189 for (Int_t iClust=0; iClust<nCaloClus; ++iClust) {
190
191 AliAODCaloCluster * cluster = aod->GetCaloCluster(iClust);
192
193 //Check if it is a EMCAL cluster
194 if(!cluster->IsEMCALCluster()) continue ;
195
6eb2a715 196 if(ClusterContainsBadChannel(cluster->GetCellsAbsId(), cluster->GetNCells())) continue;
197
198
375cec9b 199 Int_t id = cluster->GetID();
200 Float_t energy = cluster->E();
201 cluster->GetPosition(posF);
202 Char_t ttype = cluster->GetType();
375cec9b 203 AliAODCaloCluster *caloCluster = new(caloClusters[jClusters++])
204 AliAODCaloCluster(id,
205 0,
206 0x0,
207 energy,
208 posF,
209 NULL,
210 ttype);
6eb2a715 211
375cec9b 212 caloCluster->SetCaloCluster(cluster->GetDistToBadChannel(),
213 cluster->GetDispersion(),
214 cluster->GetM20(), cluster->GetM02(),
215 cluster->GetEmcCpvDistance(),
216 cluster->GetNExMax(),cluster->GetTOF()) ;
217
218 caloCluster->SetPIDFromESD(cluster->PID());
219 caloCluster->SetNCells(cluster->GetNCells());
220 caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
cf028690 221
6eb2a715 222 caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
375cec9b 223
224 }
225
6eb2a715 226 caloClusters.Expand(jClusters); // resize TObjArray
375cec9b 227 // end of loop on calo clusters
228
229 // fill EMCAL cell info
230 if (aod->GetEMCALCells()) { // protection against missing AOD information
231 AliAODCaloCells &aodinEMcells = *(aod->GetEMCALCells());
232 Int_t nEMcell = aodinEMcells.GetNumberOfCells() ;
233
234 AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
235 aodEMcells.CreateContainer(nEMcell);
236 aodEMcells.SetType(AliAODCaloCells::kEMCAL);
237
238 Double_t calibFactor = 1;
375cec9b 239 for (Int_t iCell = 0; iCell < nEMcell; iCell++) {
240 aodEMcells.SetCell(iCell,aodinEMcells.GetCellNumber(iCell),aodinEMcells.GetAmplitude(iCell)*calibFactor);
241 }
242 aodEMcells.Sort();
243
244 }
245
246}
247
248//__________________________________________________
249void AliAnalysisTaskEMCALPi0CalibSelection::CreateAODFromESD()
250{
251
252 // Copy Header, Vertex, CaloClusters and CaloCells from ESDs to AODs
253
254 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
255
256 // set arrays and pointers
257 Float_t posF[3];
258 Double_t pos[3];
259
260 Double_t covVtx[6];
261
262 for (Int_t i = 0; i < 6; i++) covVtx[i] = 0.;
263
264 // Update the header
265
266 AliAODHeader* header = AODEvent()->GetHeader();
267 header->SetRunNumber(esd->GetRunNumber());
268 header->SetBunchCrossNumber(esd->GetBunchCrossNumber());
269 header->SetOrbitNumber(esd->GetOrbitNumber());
270 header->SetPeriodNumber(esd->GetPeriodNumber());
271 header->SetEventType(esd->GetEventType());
272 header->SetMuonMagFieldScale(-999.); // FIXME
273 header->SetCentrality(-999.); // FIXME
274
275
276 header->SetTriggerMask(esd->GetTriggerMask());
277 header->SetTriggerCluster(esd->GetTriggerCluster());
278 header->SetMagneticField(esd->GetMagneticField());
279 header->SetZDCN1Energy(esd->GetZDCN1Energy());
280 header->SetZDCP1Energy(esd->GetZDCP1Energy());
281 header->SetZDCN2Energy(esd->GetZDCN2Energy());
282 header->SetZDCP2Energy(esd->GetZDCP2Energy());
283 header->SetZDCEMEnergy(esd->GetZDCEMEnergy(0),esd->GetZDCEMEnergy(1));
284 Float_t diamxy[2]={esd->GetDiamondX(),esd->GetDiamondY()};
285 Float_t diamcov[3]; esd->GetDiamondCovXY(diamcov);
286 header->SetDiamond(diamxy,diamcov);
287 //
288 //
289 Int_t nVertices = 1 ;/* = prim. vtx*/;
290 Int_t nCaloClus = esd->GetNumberOfCaloClusters();
291
292 AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
293
294 // Access to the AOD container of vertices
295 TClonesArray &vertices = *(AODEvent()->GetVertices());
296 Int_t jVertices=0;
297
298 // Add primary vertex. The primary tracks will be defined
299 // after the loops on the composite objects (V0, cascades, kinks)
300 const AliESDVertex *vtx = esd->GetPrimaryVertex();
301
302 vtx->GetXYZ(pos); // position
303 vtx->GetCovMatrix(covVtx); //covariance matrix
304
305 AliAODVertex * primary = new(vertices[jVertices++])
306 AliAODVertex(pos, covVtx, vtx->GetChi2toNDF(), NULL, -1, AliAODVertex::kPrimary);
307 primary->SetName(vtx->GetName());
308 primary->SetTitle(vtx->GetTitle());
309
310 // Access to the AOD container of clusters
311 TClonesArray &caloClusters = *(AODEvent()->GetCaloClusters());
312 Int_t jClusters=0;
313
314 for (Int_t iClust=0; iClust<nCaloClus; ++iClust) {
315
316 AliESDCaloCluster * cluster = esd->GetCaloCluster(iClust);
317
318 //Check which calorimeter information we want to keep.
319 if(!cluster->IsEMCAL()) continue ;
6eb2a715 320 if(ClusterContainsBadChannel(cluster->GetCellsAbsId(), cluster->GetNCells())) continue;
321
375cec9b 322 Int_t id = cluster->GetID();
323 Float_t energy = cluster->E();
324 cluster->GetPosition(posF);
325
326 AliAODCaloCluster *caloCluster = new(caloClusters[jClusters++])
327 AliAODCaloCluster(id,
328 0,
329 0x0,
330 energy,
331 posF,
332 NULL,
333 AliAODCluster::kEMCALClusterv1);
334
335 caloCluster->SetCaloCluster(cluster->GetDistanceToBadChannel(),
336 cluster->GetClusterDisp(),
337 cluster->GetM20(), cluster->GetM02(),
338 cluster->GetEmcCpvDistance(),
339 cluster->GetNExMax(),cluster->GetTOF()) ;
340
341 caloCluster->SetPIDFromESD(cluster->GetPid());
342 caloCluster->SetNCells(cluster->GetNCells());
343 caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
6eb2a715 344 caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
345
375cec9b 346 }
6eb2a715 347
348 caloClusters.Expand(jClusters); // resize TObjArray
375cec9b 349 // end of loop on calo clusters
350
351 // fill EMCAL cell info
352
353 if( esd->GetEMCALCells()) { // protection against missing ESD information
354 AliESDCaloCells &esdEMcells = *(esd->GetEMCALCells());
355 Int_t nEMcell = esdEMcells.GetNumberOfCells() ;
356
357 AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
358 aodEMcells.CreateContainer(nEMcell);
359 aodEMcells.SetType(AliAODCaloCells::kEMCAL);
360
6eb2a715 361 Double_t calibFactor = 1;
375cec9b 362 for (Int_t iCell = 0; iCell < nEMcell; iCell++) {
363 aodEMcells.SetCell(iCell,esdEMcells.GetCellNumber(iCell),esdEMcells.GetAmplitude(iCell)*calibFactor);
364 }
365 aodEMcells.Sort();
366
367 }
368
369}
370
371//__________________________________________________
372void AliAnalysisTaskEMCALPi0CalibSelection::UserCreateOutputObjects()
373{
cf028690 374 //Create output container, init geometry and calibration
375
cf028690 376 fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName) ;
377
375cec9b 378 fOutputContainer = new TList();
379
380 char hname[128], htitl[128];
381
6eb2a715 382 for(Int_t iMod=0; iMod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); iMod++) {
cf028690 383 for(Int_t iRow=0; iRow<AliEMCALGeoParams::fgkEMCALRows; iRow++) {
384 for(Int_t iCol=0; iCol<AliEMCALGeoParams::fgkEMCALCols; iCol++) {
385 sprintf(hname,"%d_%d_%d",iMod,iCol,iRow);
386 sprintf(htitl,"Two-gamma inv. mass for super mod %d, cell(col,row)=(%d,%d)",iMod,iCol,iRow);
6eb2a715 387 fHmpi0[iMod][iCol][iRow] = new TH1F(hname,htitl,fNbins,fMinBin,fMaxBin);
cf028690 388 fOutputContainer->Add(fHmpi0[iMod][iCol][iRow]);
375cec9b 389 }
390 }
391 }
392
6eb2a715 393 fHmgg = new TH1F("hmgg","2-cluster invariant mass",fNbins,fMinBin,fMaxBin);
375cec9b 394 fOutputContainer->Add(fHmgg);
6eb2a715 395
396 fhNEvents = new TH1I("hNEvents", "Number of analyzed events" , 1 , 0 , 1 ) ;
397 fOutputContainer->Add(fhNEvents);
398
399// fCalibData = new AliEMCALCalibData();
400
cf028690 401 PostData(1,fOutputContainer);
375cec9b 402
403}
404
405//__________________________________________________
406void AliAnalysisTaskEMCALPi0CalibSelection::UserExec(Option_t* /* option */)
407{
408 //Analysis per event.
409 if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection <<< Event %d >>>\n",(Int_t)Entry());
410
6eb2a715 411 fhNEvents->Fill(0); //Event analyzed
412
375cec9b 413 AliAODEvent* aod = 0x0;
6eb2a715 414 Bool_t kAOD = kFALSE;
415 if(!strcmp(InputEvent()->GetName(),"AliAODEvent")) kAOD=kTRUE;
416 Bool_t kESD = kFALSE;
417 if(!strcmp(InputEvent()->GetName(),"AliESDEvent")) kESD=kTRUE;
418
419 if(kAOD){
375cec9b 420 //Input are ESDs
421 aod = dynamic_cast<AliAODEvent*>(InputEvent());
422 // Create new AOD with only CaloClusters and CaloCells
423 if(fCopyAOD) CreateAODFromAOD();
424 }
6eb2a715 425 else if(kESD) {
375cec9b 426 //Input are ESDs
427 aod = AODEvent();
428 // Create AOD with CaloClusters and use it as input.
429 // If filtering task is already executed, this is not needed.
430 if(fCopyAOD) CreateAODFromESD();
431 }
432 else {
433 printf("AliAnalysisTaskEMCALPi0CalibSelection: Unknown event type, STOP!\n");
434 abort();
435 }
436
6eb2a715 437 Double_t v[3];// = {aod->GetVertex(0)->GetX(),aod->GetVertex(0)->GetY(),aod->GetVertex(0)->GetZ()}; //to check!!
cf028690 438 aod->GetPrimaryVertex()->GetXYZ(v) ;
375cec9b 439 //TVector3 vtx(v);
440
441 //if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Vertex: (%.3f,%.3f,%.3f)\n",vtx.X(),vtx.Y(),vtx.Z());
442 if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Vertex: (%.3f,%.3f,%.3f)\n",v[0],v[1],v[2]);
443
444 Int_t runNum = aod->GetRunNumber();
445 if(DebugLevel() > 1) printf("Run number: %d\n",runNum);
446
6eb2a715 447 //FIXME Not neede the matrices for the moment MEFIX
375cec9b 448 //Get the matrix with geometry information
449 //Still not implemented in AOD, just a workaround to be able to work at least with ESDs
6eb2a715 450// if(!strcmp(InputEvent()->GetName(),"AliAODEvent")) {
451// if(DebugLevel() > 1)
452// printf("AliAnalysisTaskEMCALPi0CalibSelection Use ideal geometry, values geometry matrix not kept in AODs.\n");
453// }
454// else{
455// if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Load Misaligned matrices. \n");
456// AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent()) ;
457// for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
458// if(esd->GetEMCALMatrix(mod)) fEMCALGeo->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
459// }
460// }
375cec9b 461
462 if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Will use fLogWeight %.3f .\n",fLogWeight);
cf028690 463
6eb2a715 464 Int_t iSupMod1 = -1;
465 Int_t iphi1 = -1;
466 Int_t ieta1 = -1;
467 Int_t iSupMod2 = -1;
468 Int_t iphi2 = -1;
469 Int_t ieta2 = -1;
470
375cec9b 471 TLorentzVector p1;
472 TLorentzVector p2;
473 TLorentzVector p12;
474
375cec9b 475 TRefArray * caloClustersArr = new TRefArray();
476 aod->GetEMCALClusters(caloClustersArr);
477
478 const Int_t kNumberOfEMCALClusters = caloClustersArr->GetEntries() ;
cf028690 479 if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection - N CaloClusters: %d \n", kNumberOfEMCALClusters);
375cec9b 480
481 // EMCAL cells
482 AliAODCaloCells *emCells = aod->GetEMCALCells();
6eb2a715 483
375cec9b 484 // loop over EMCAL clusters
485 for(Int_t iClu=0; iClu<kNumberOfEMCALClusters; iClu++) {
486
487 AliAODCaloCluster *c1 = (AliAODCaloCluster *) caloClustersArr->At(iClu);
6eb2a715 488 if(!fCopyAOD && kESD && ClusterContainsBadChannel(c1->GetCellsAbsId(), c1->GetNCells())) continue;
489
490 Float_t e1i = c1->E(); // cluster energy before correction
491 if(e1i < fEmin) continue;
492 else if(e1i > fEmax) continue;
493 else if (c1->GetNCells() < fMinNCells) continue;
494
495 if(DebugLevel() > 2)
496 {
497 printf("Std : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",iClu,e1i, c1->GetDispersion(),c1->GetM02(),c1->GetM20());
498 Float_t pos[]={0,0,0};
499 c1->GetPosition(pos);
500 printf("Std : i %d, x %f, y %f, z %f\n",iClu, pos[0], pos[1], pos[2]);
501 }
502
503 //AliEMCALAodCluster clu1(*c1);
504 //clu1.Recalibrate(fCalibData, emCells, fEMCALGeoName);
505 //clu1.EvalEnergy();
cf028690 506 //clu1.EvalAll(fLogWeight, fEMCALGeoName);
6eb2a715 507 if(IsRecalibrationOn()) {
508 Float_t energy = RecalibrateClusterEnergy(c1, emCells);
509 //clu1.SetE(energy);
510 c1->SetE(energy);
511 }
375cec9b 512
6eb2a715 513 //Float_t e1ii = clu1.E(); // cluster energy after correction
514 Float_t e1ii = c1->E(); // cluster energy after correction
375cec9b 515
6eb2a715 516 if(DebugLevel() > 2)
517 {
518 //printf("Recal: i %d, E %f, dispersion %f, m02 %f, m20 %f\n",iClu,e1ii, clu1.GetDispersion(),clu1.GetM02(),clu1.GetM20());
519 printf("Recal: i %d, E %f, dispersion %f, m02 %f, m20 %f\n",iClu,e1ii, c1->GetDispersion(),c1->GetM02(),c1->GetM20());
520 Float_t pos2[]={0,0,0};
521 //clu1.GetPosition(pos2);
522 c1->GetPosition(pos2);
523 printf("Recal: i %d, x %f, y %f, z %f\n",iClu, pos2[0], pos2[1], pos2[2]);
524 }
525
526 //clu1.GetMomentum(p1,v);
527 c1->GetMomentum(p1,v);
528
529 //Get tower with maximum energy and fill in the end the pi0 histogram for this cell.
530 //MaxEnergyCellPos(emCells,&clu1,iSupMod1,ieta1,iphi1);
531 MaxEnergyCellPos(emCells,c1,iSupMod1,ieta1,iphi1);
375cec9b 532
533 for (Int_t jClu=iClu; jClu<kNumberOfEMCALClusters; jClu++) {
534 AliAODCaloCluster *c2 = (AliAODCaloCluster *) caloClustersArr->At(jClu);
375cec9b 535 if(c2->IsEqual(c1)) continue;
6eb2a715 536 if(!fCopyAOD && kESD && ClusterContainsBadChannel(c2->GetCellsAbsId(), c2->GetNCells())) continue;
375cec9b 537
538 Float_t e2i = c2->E();
6eb2a715 539 if(e2i < fEmin) continue;
540 else if (e2i > fEmax) continue;
541 else if (c2->GetNCells() < fMinNCells) continue;
542
543 //AliEMCALAodCluster clu2(*c2);
544 //clu2.Recalibrate(fCalibData, emCells,fEMCALGeoName);
545 //clu2.EvalEnergy();
cf028690 546 //clu2.EvalAll(fLogWeight,fEMCALGeoName);
6eb2a715 547 if(IsRecalibrationOn()) {
548 Float_t energy = RecalibrateClusterEnergy(c2, emCells);
549 //clu2.SetE(energy);
550 c2->SetE(energy);
551 }
552
553 Float_t e2ii = c2->E();
554
555 //clu2.GetMomentum(p2,v);
556 c2->GetMomentum(p2,v);
557
558 //Get tower with maximum energy and fill in the end the pi0 histogram for this cell.
559 //MaxEnergyCellPos(emCells,&clu2,iSupMod2,ieta2,iphi2);
560 MaxEnergyCellPos(emCells,c2,iSupMod2,ieta2,iphi2);
375cec9b 561
562 p12 = p1+p2;
563 Float_t invmass = p12.M()*1000;
375cec9b 564
6eb2a715 565 if(invmass < fMaxBin && invmass > fMinBin){
566 fHmgg->Fill(invmass);
567
568 if (fGroupNCells == 0){
569 fHmpi0[iSupMod1][ieta1][iphi1]->Fill(invmass);
570 fHmpi0[iSupMod2][ieta2][iphi2]->Fill(invmass);
571 }
572 else {
573 //printf("Regroup N %d, eta1 %d, phi1 %d, eta2 %d, phi2 %d \n",fGroupNCells, ieta1, iphi1, ieta2, iphi2);
574 for (Int_t i = -fGroupNCells; i < fGroupNCells+1; i++) {
575 for (Int_t j = -fGroupNCells; j < fGroupNCells+1; j++) {
576 //printf("\t i %d, j %d\n",i,j);
577 if((ieta1+i >= 0) && (iphi1+j >= 0) && (ieta1+i < 48) && (iphi1+j < 24)){
578 //printf("\t \t eta1+i %d, phi1+j %d\n", ieta1+i, iphi1+j);
579 fHmpi0[iSupMod1][ieta1+i][iphi1+j]->Fill(invmass);
580 }
581 if((ieta2+i >= 0) && (iphi2+j >= 0) && (ieta2+i < 48) && (iphi2+j < 24)){
582 //printf("\t \t eta2+i %d, phi2+j %d\n", ieta2+i, iphi2+j);
583 fHmpi0[iSupMod2][ieta2+i][iphi2+j]->Fill(invmass);
584 }
585 }// j loop
586 }//i loop
587 }//group cells
588
589 if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Mass in (SM%d,%d,%d) and (SM%d,%d,%d): %.3f GeV E1_i=%f E1_ii=%f E2_i=%f E2_ii=%f\n",
590 iSupMod1,iphi1,ieta1,iSupMod2,iphi2,ieta2,p12.M(),e1i,e1ii,e2i,e2ii);
591 }
592
375cec9b 593 }
594
595 } // end of loop over EMCAL clusters
596
597 delete caloClustersArr;
6eb2a715 598
375cec9b 599 PostData(1,fOutputContainer);
6eb2a715 600
375cec9b 601}
602
603//__________________________________________________
6eb2a715 604void AliAnalysisTaskEMCALPi0CalibSelection::MaxEnergyCellPos(AliAODCaloCells* const cells, AliAODCaloCluster* const clu, Int_t& iSupMod, Int_t& ieta, Int_t& iphi)
375cec9b 605{
606 //For a given CaloCluster calculates the absId of the cell
607 //with maximum energy deposit.
608
6eb2a715 609 Double_t eMax = -1.;
610 Double_t eCell = -1.;
611 Float_t fraction = 1.;
612 Int_t cellAbsId = -1;
613 Float_t recalFactor = 1.;
614
615 Int_t maxId = -1;
616 Int_t iTower = -1;
617 Int_t iIphi = -1;
618 Int_t iIeta = -1;
619
375cec9b 620 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
6eb2a715 621 cellAbsId = clu->GetCellAbsId(iDig);
622 fraction = clu->GetCellAmplitudeFraction(iDig);
623 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
624 if(IsRecalibrationOn()) {
625 Int_t imodrc = -1, iphirc = -1, ietarc =-1;
626 Int_t iTowerrc = -1, iIphirc = -1, iIetarc =-1;
627 fEMCALGeo->GetCellIndex(cellAbsId,imodrc,iTowerrc,iIphirc,iIetarc);
628 fEMCALGeo->GetCellPhiEtaIndexInSModule(imodrc,iTowerrc,iIphirc, iIetarc,iphirc,ietarc);
629 recalFactor = GetEMCALChannelRecalibrationFactor(imodrc,ietarc,iphirc);
630 }
631 eCell = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
632 //printf("Max cell? cell %d, amplitude org %f, fraction %f, recalibration %f, amplitude new %f \n",cellAbsId, cells->GetCellAmplitude(cellAbsId), fraction, recalFactor, eCell) ;
633
634 if(eCell > eMax) {
635 eMax = eCell;
375cec9b 636 maxId = cellAbsId;
6eb2a715 637 //printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell);
375cec9b 638 }
639 }
6eb2a715 640
641 //Get from the absid the supermodule, tower and eta/phi numbers
642 fEMCALGeo->GetCellIndex(maxId,iSupMod,iTower,iIphi,iIeta);
643 //Gives SuperModule and Tower numbers
644 fEMCALGeo->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
645 iIphi, iIeta,iphi,ieta);
646
647 //printf("\t Max : cell %d, iSupMod %d, ieta %d, iphi %d \n",maxId,iSupMod, ieta,iphi);
648
375cec9b 649}
650
651//__________________________________________________
6eb2a715 652//void AliAnalysisTaskEMCALPi0CalibSelection::SetCalibCorrections(AliEMCALCalibData* const cdata)
653//{
654// //Set new correction factors (~1) to calibration coefficients, delete previous.
655//
656// if(fCalibData) delete fCalibData;
657// fCalibData = cdata;
658//
659//}
660
661//________________________________________________________________
662void AliAnalysisTaskEMCALPi0CalibSelection::InitEMCALBadChannelStatusMap(){
663 //Init EMCAL bad channels map
664 if(DebugLevel() > 0 )printf("AliAnalysisTaskEMCALPi0CalibSelection::InitEMCALBadChannelStatusMap()\n");
665 //In order to avoid rewriting the same histograms
666 Bool_t oldStatus = TH1::AddDirectoryStatus();
667 TH1::AddDirectory(kFALSE);
668
669 fEMCALBadChannelMap = new TObjArray(12);
670 //TH2F * hTemp = new TH2I("EMCALBadChannelMap","EMCAL SuperModule bad channel map", 48, 0, 48, 24, 0, 24);
671 for (int i = 0; i < 12; i++) {
672 fEMCALBadChannelMap->Add(new TH2I(Form("EMCALBadChannelMap_Mod%d",i),Form("EMCALBadChannelMap_Mod%d",i), 48, 0, 48, 24, 0, 24));
673 //fEMCALBadChannelMap->Add((TH2I*)hTemp->Clone(Form("EMCALBadChannelMap_Mod%d",i)));
674 }
675
676 //delete hTemp;
677
678 fEMCALBadChannelMap->SetOwner(kTRUE);
679 fEMCALBadChannelMap->Compress();
680
681 //In order to avoid rewriting the same histograms
682 TH1::AddDirectory(oldStatus);
683}
684
685
686//_________________________________________________________________________________________________________
687Bool_t AliAnalysisTaskEMCALPi0CalibSelection::ClusterContainsBadChannel(UShort_t* cellList, Int_t nCells){
688 // Check that in the cluster cells, there is no bad channel of those stored
689 // in fEMCALBadChannelMap or fPHOSBadChannelMap
690
691 if(!fRemoveBadChannels) return kFALSE;
692 if(!fEMCALBadChannelMap) return kFALSE;
693
694 Int_t icol = -1;
695 Int_t irow = -1;
696 Int_t imod = -1;
697 for(Int_t iCell = 0; iCell<nCells; iCell++){
698
699 //Get the column and row
700 Int_t iTower = -1, iIphi = -1, iIeta = -1;
701 fEMCALGeo->GetCellIndex(cellList[iCell],imod,iTower,iIphi,iIeta);
702 if(fEMCALBadChannelMap->GetEntries() <= imod) continue;
703 fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
704 if(GetEMCALChannelStatus(imod, icol, irow))return kTRUE;
705
706 }// cell cluster loop
707
708 return kFALSE;
709
710}
375cec9b 711
6eb2a715 712
713//________________________________________________________________
714void AliAnalysisTaskEMCALPi0CalibSelection::InitEMCALRecalibrationFactors(){
715 //Init EMCAL recalibration factors
716 if(DebugLevel() > 0 )printf("AliAnalysisTaskEMCALPi0CalibSelection::InitEMCALRecalibrationFactors()\n");
717 //In order to avoid rewriting the same histograms
718 Bool_t oldStatus = TH1::AddDirectoryStatus();
719 TH1::AddDirectory(kFALSE);
720
721 fEMCALRecalibrationFactors = new TObjArray(12);
722 for (int i = 0; i < 12; i++) fEMCALRecalibrationFactors->Add(new TH2F(Form("EMCALRecalFactors_SM%d",i),Form("EMCALRecalFactors_SM%d",i), 48, 0, 48, 24, 0, 24));
723 //Init the histograms with 1
724 for (Int_t sm = 0; sm < 12; sm++) {
725 for (Int_t i = 0; i < 48; i++) {
726 for (Int_t j = 0; j < 24; j++) {
727 SetEMCALChannelRecalibrationFactor(sm,i,j,1.);
728 }
729 }
730 }
731 fEMCALRecalibrationFactors->SetOwner(kTRUE);
732 fEMCALRecalibrationFactors->Compress();
733
734 //In order to avoid rewriting the same histograms
735 TH1::AddDirectory(oldStatus);
736}
737
738//________________________________________________________________
739Float_t AliAnalysisTaskEMCALPi0CalibSelection::RecalibrateClusterEnergy(AliAODCaloCluster * cluster, AliAODCaloCells * cells){
740 // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
741 // AOD case
742
743 if(!cells) {
744 printf("AliAnalysisTaskEMCALPi0CalibSelection::RecalibrateClusterEnergy(AOD) - Cells pointer does not exist, stop!");
745 abort();
746 }
747
748 //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
749 UShort_t * index = cluster->GetCellsAbsId() ;
750 Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
751 Int_t ncells = cluster->GetNCells();
752
753 //Initialize some used variables
754 Float_t energy = 0;
755 Int_t absId = -1;
756 Int_t icol = -1, irow = -1, imod=1;
757 Float_t factor = 1, frac = 0;
758
759 //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
760 for(Int_t icell = 0; icell < ncells; icell++){
761 absId = index[icell];
762 frac = fraction[icell];
763 if(frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off
764 Int_t iTower = -1, iIphi = -1, iIeta = -1;
765 fEMCALGeo->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
766 if(fEMCALRecalibrationFactors->GetEntries() <= imod) continue;
767 fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
768 factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
769 if(DebugLevel()>2)
770 printf("AliAnalysisTaskEMCALPi0CalibSelection::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n",
771 imod,icol,irow,frac,factor,cells->GetCellAmplitude(absId));
772
773 energy += cells->GetCellAmplitude(absId)*factor*frac;
774 }
775
776 if(DebugLevel()>1)
777 printf("AliAnalysisTaskEMCALPi0CalibSelection::RecalibrateClusterEnergy - Energy before %f, after %f\n",cluster->E(),energy);
778
779 return energy;
375cec9b 780
781}
782