1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: Boris Polishchuk *
5 * Adapted to AOD reading by Gustavo Conesa *
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 **************************************************************************/
16 //---------------------------------------------------------------------------//
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.//
23 //---------------------------------------------------------------------------//
26 //#include <Riostream.h>
28 #include "TLorentzVector.h"
29 //#include "TVector3.h"
30 #include "TRefArray.h"
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"
46 ClassImp(AliAnalysisTaskEMCALPi0CalibSelection)
48 //__________________________________________________
49 AliAnalysisTaskEMCALPi0CalibSelection::AliAnalysisTaskEMCALPi0CalibSelection() :
50 AliAnalysisTaskSE(),fEMCALGeo(0x0),fCalibData(0x0), fEmin(0.), fLogWeight(4.5), fCopyAOD(kFALSE),
51 fEMCALGeoName("EMCAL_COMPLETE"), fOldData(kFALSE),fOutputContainer(0x0),fHmgg(0x0)
53 //Default constructor.
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++) {
58 fHmpi0[iMod][iX][iZ]=0;
65 //__________________________________________________
66 AliAnalysisTaskEMCALPi0CalibSelection::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)
70 //Named constructor which should be used.
72 DefineOutput(1,TList::Class());
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++) {
77 fHmpi0[iMod][iX][iZ]=0;
84 //__________________________________________________
85 AliAnalysisTaskEMCALPi0CalibSelection::~AliAnalysisTaskEMCALPi0CalibSelection()
90 fOutputContainer->Delete() ;
91 delete fOutputContainer ;
94 if(fCalibData) delete fCalibData;
95 if(fEMCALGeo) delete fEMCALGeo;
100 //__________________________________________________
101 void AliAnalysisTaskEMCALPi0CalibSelection::CreateAODFromAOD()
103 // Copy AOD header, vertex, CaloClusters and CaloCells to output AOD
105 AliAODEvent* aod = dynamic_cast<AliAODEvent*>(InputEvent());
107 // set arrays and pointers
113 for (Int_t i = 0; i < 6; i++) covVtx[i] = 0.;
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());
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);
139 Int_t nVertices = 1 ;/* = prim. vtx*/;
140 Int_t nCaloClus = aod->GetNCaloClusters();
142 AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
144 // Access to the AOD container of vertices
145 TClonesArray &vertices = *(AODEvent()->GetVertices());
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();
152 vtx->GetXYZ(pos); // position
153 vtx->GetCovMatrix(covVtx); //covariance matrix
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());
160 // Access to the AOD container of clusters
161 TClonesArray &caloClusters = *(AODEvent()->GetCaloClusters());
164 for (Int_t iClust=0; iClust<nCaloClus; ++iClust) {
166 AliAODCaloCluster * cluster = aod->GetCaloCluster(iClust);
168 //Check if it is a EMCAL cluster
169 if(!cluster->IsEMCALCluster()) continue ;
171 Int_t id = cluster->GetID();
172 Float_t energy = cluster->E();
173 cluster->GetPosition(posF);
174 Char_t ttype = cluster->GetType();
176 AliAODCaloCluster *caloCluster = new(caloClusters[jClusters++])
177 AliAODCaloCluster(id,
185 caloCluster->SetCaloCluster(cluster->GetDistToBadChannel(),
186 cluster->GetDispersion(),
187 cluster->GetM20(), cluster->GetM02(),
188 cluster->GetEmcCpvDistance(),
189 cluster->GetNExMax(),cluster->GetTOF()) ;
191 caloCluster->SetPIDFromESD(cluster->PID());
192 caloCluster->SetNCells(cluster->GetNCells());
193 caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
195 caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
198 Double_t *fraction = cluster->GetCellsAmplitudeFraction();
199 for(Int_t i = 0; i < cluster->GetNCells() ; i++) fraction[i] = 1;
200 caloCluster->SetCellsAmplitudeFraction(fraction);
205 caloClusters.Expand(jClusters); // resize TObjArray to 'remove' slots for pseudo clusters
206 // end of loop on calo clusters
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() ;
213 AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
214 aodEMcells.CreateContainer(nEMcell);
215 aodEMcells.SetType(AliAODCaloCells::kEMCAL);
217 Double_t calibFactor = 1;
218 if(fOldData) calibFactor = 0.0153;
220 for (Int_t iCell = 0; iCell < nEMcell; iCell++) {
221 aodEMcells.SetCell(iCell,aodinEMcells.GetCellNumber(iCell),aodinEMcells.GetAmplitude(iCell)*calibFactor);
229 //__________________________________________________
230 void AliAnalysisTaskEMCALPi0CalibSelection::CreateAODFromESD()
233 // Copy Header, Vertex, CaloClusters and CaloCells from ESDs to AODs
235 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
237 // set arrays and pointers
243 for (Int_t i = 0; i < 6; i++) covVtx[i] = 0.;
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
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);
270 Int_t nVertices = 1 ;/* = prim. vtx*/;
271 Int_t nCaloClus = esd->GetNumberOfCaloClusters();
273 AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
275 // Access to the AOD container of vertices
276 TClonesArray &vertices = *(AODEvent()->GetVertices());
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();
283 vtx->GetXYZ(pos); // position
284 vtx->GetCovMatrix(covVtx); //covariance matrix
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());
291 // Access to the AOD container of clusters
292 TClonesArray &caloClusters = *(AODEvent()->GetCaloClusters());
295 for (Int_t iClust=0; iClust<nCaloClus; ++iClust) {
297 AliESDCaloCluster * cluster = esd->GetCaloCluster(iClust);
299 //Check which calorimeter information we want to keep.
300 if(!cluster->IsEMCAL()) continue ;
302 Int_t id = cluster->GetID();
303 Float_t energy = cluster->E();
304 cluster->GetPosition(posF);
306 AliAODCaloCluster *caloCluster = new(caloClusters[jClusters++])
307 AliAODCaloCluster(id,
313 AliAODCluster::kEMCALClusterv1);
315 caloCluster->SetCaloCluster(cluster->GetDistanceToBadChannel(),
316 cluster->GetClusterDisp(),
317 cluster->GetM20(), cluster->GetM02(),
318 cluster->GetEmcCpvDistance(),
319 cluster->GetNExMax(),cluster->GetTOF()) ;
321 caloCluster->SetPIDFromESD(cluster->GetPid());
322 caloCluster->SetNCells(cluster->GetNCells());
323 caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
325 caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
328 Double_t *fraction = cluster->GetCellsAmplitudeFraction();
329 for(Int_t i = 0; i < cluster->GetNCells() ; i++) fraction[i] = 1;
330 caloCluster->SetCellsAmplitudeFraction(fraction);
333 caloClusters.Expand(jClusters); // resize TObjArray to 'remove' slots for pseudo clusters
334 // end of loop on calo clusters
336 // fill EMCAL cell info
338 if( esd->GetEMCALCells()) { // protection against missing ESD information
339 AliESDCaloCells &esdEMcells = *(esd->GetEMCALCells());
340 Int_t nEMcell = esdEMcells.GetNumberOfCells() ;
342 AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
343 aodEMcells.CreateContainer(nEMcell);
344 aodEMcells.SetType(AliAODCaloCells::kEMCAL);
346 Double_t calibFactor = 1;
347 if(fOldData) calibFactor = 0.0153;
349 for (Int_t iCell = 0; iCell < nEMcell; iCell++) {
350 aodEMcells.SetCell(iCell,esdEMcells.GetCellNumber(iCell),esdEMcells.GetAmplitude(iCell)*calibFactor);
358 //__________________________________________________
359 void AliAnalysisTaskEMCALPi0CalibSelection::UserCreateOutputObjects()
361 //Create output container
362 fOutputContainer = new TList();
364 char hname[128], htitl[128];
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]);
377 fHmgg = new TH1F("hmgg","2-cluster invariant mass",100,0.,300.);
378 fOutputContainer->Add(fHmgg);
380 fCalibData = new AliEMCALCalibData();
383 AliFatal("Calibration parameters not found in CDB!");
385 printf("Get Geometry : %s\n", fEMCALGeoName.Data());
386 fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName) ;
390 //__________________________________________________
391 void AliAnalysisTaskEMCALPi0CalibSelection::UserExec(Option_t* /* option */)
393 //Analysis per event.
394 if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection <<< Event %d >>>\n",(Int_t)Entry());
396 AliAODEvent* aod = 0x0;
397 if(!strcmp(InputEvent()->GetName(),"AliAODEvent")) {
399 aod = dynamic_cast<AliAODEvent*>(InputEvent());
400 // Create new AOD with only CaloClusters and CaloCells
401 if(fCopyAOD) CreateAODFromAOD();
403 else if(!strcmp(InputEvent()->GetName(),"AliESDEvent")) {
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();
411 printf("AliAnalysisTaskEMCALPi0CalibSelection: Unknown event type, STOP!\n");
415 Double_t v[] = {aod->GetVertex(0)->GetX(),aod->GetVertex(0)->GetY(),aod->GetVertex(0)->GetZ()}; //to check!!
416 //aod->GetVertex()->GetXYZ(v) ;
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]);
422 Int_t runNum = aod->GetRunNumber();
423 if(DebugLevel() > 1) printf("Run number: %d\n",runNum);
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")) {
429 printf("AliAnalysisTaskEMCALPi0CalibSelection Use ideal geometry, values geometry matrix not kept in AODs.\n");
432 if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Load Misaligned matrices. \n");
433 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent()) ;
434 for(Int_t mod=0; mod < 12; mod++){
435 if(esd->GetEMCALMatrix(mod)) fEMCALGeo->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
439 if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Will use fLogWeight %.3f .\n",fLogWeight);
456 TRefArray * caloClustersArr = new TRefArray();
457 aod->GetEMCALClusters(caloClustersArr);
459 const Int_t kNumberOfEMCALClusters = caloClustersArr->GetEntries() ;
460 if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection CaloClusters: %d\n", kNumberOfEMCALClusters);
463 AliAODCaloCells *emCells = aod->GetEMCALCells();
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);
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);
484 // loop over EMCAL clusters
485 for(Int_t iClu=0; iClu<kNumberOfEMCALClusters; iClu++) {
487 AliAODCaloCluster *c1 = (AliAODCaloCluster *) caloClustersArr->At(iClu);
488 if(!c1->IsEMCALCluster()) continue; // EMCAL cluster!
490 Float_t e1i = c1->E(); // cluster energy before correction
491 if(e1i<fEmin) continue;
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]);
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) ;
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]);
513 clu1.GetMomentum(p1,v);
515 MaxEnergyCellPos(emCells,&clu1,maxId);
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);
523 Float_t e1ii = clu1.E(); // cluster energy after correction
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;
530 Float_t e2i = c2->E();
531 if(e2i<fEmin) continue;
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());
540 clu2.GetMomentum(p2,v);
541 Float_t e2ii = clu2.E();
544 Float_t invmass = p12.M()*1000;
546 fHmgg->Fill(invmass);
548 //printf("iSM %d, ieta %d, iphi %d, mass %f \n",iSupMod, ieta, iphi, invmass);
550 fHmpi0[iSupMod][ieta][iphi]->Fill(invmass);
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);
557 } // end of loop over EMCAL clusters
559 delete caloClustersArr;
560 PostData(1,fOutputContainer);
563 //__________________________________________________
564 void AliAnalysisTaskEMCALPi0CalibSelection::MaxEnergyCellPos(AliAODCaloCells* const cells, AliAODCaloCluster* const clu, Int_t& maxId)
566 //For a given CaloCluster calculates the absId of the cell
567 //with maximum energy deposit.
569 Double_t eMax = -111;
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);
582 //__________________________________________________
583 void AliAnalysisTaskEMCALPi0CalibSelection::SetCalibCorrections(AliEMCALCalibData* const cdata)
585 //Set new correction factors (~1) to calibration coefficients, delete previous.
587 if(fCalibData) delete fCalibData;