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