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 anywhere in PHOS.//
23 //---------------------------------------------------------------------------//
28 #include "TLorentzVector.h"
30 #include "TRefArray.h"
34 #include "AliAnalysisTaskPHOSPi0CalibSelection.h"
35 #include "AliAODEvent.h"
36 #include "AliESDEvent.h"
37 #include "AliESDCaloCluster.h"
38 #include "AliESDCaloCells.h"
39 #include "AliPHOSAodCluster.h"
40 #include "AliPHOSPIDv1.h"
41 #include "AliCentrality.h"
43 ClassImp(AliAnalysisTaskPHOSPi0CalibSelection)
45 //__________________________________________________
46 AliAnalysisTaskPHOSPi0CalibSelection::AliAnalysisTaskPHOSPi0CalibSelection() :
47 AliAnalysisTaskSE(),fOutputContainer(0x0),fPhosGeo(0x0),fCalibData(0x0),fHmgg(0x0),
48 fEmin(0.), fLogWeight(4.5), fCopyAOD(kFALSE)
50 //Default constructor.
52 for(Int_t iMod=0; iMod<5; iMod++) {
53 for(Int_t iX=0; iX<64; iX++) {
54 for(Int_t iZ=0; iZ<56; iZ++) {
55 fHmpi0[iMod][iX][iZ]=0;
62 //__________________________________________________
63 AliAnalysisTaskPHOSPi0CalibSelection::AliAnalysisTaskPHOSPi0CalibSelection(const char* name) :
64 AliAnalysisTaskSE(name),fOutputContainer(0x0),fPhosGeo(0x0),fCalibData(0x0),fHmgg(0x0),
65 fEmin(0.), fLogWeight(4.5), fCopyAOD(kFALSE)
67 //Named constructor which should be used.
69 DefineOutput(1,TList::Class());
71 for(Int_t iMod=0; iMod<5; iMod++) {
72 for(Int_t iX=0; iX<64; iX++) {
73 for(Int_t iZ=0; iZ<56; iZ++) {
74 fHmpi0[iMod][iX][iZ]=0;
81 //__________________________________________________
82 AliAnalysisTaskPHOSPi0CalibSelection::~AliAnalysisTaskPHOSPi0CalibSelection()
87 fOutputContainer->Delete() ;
88 delete fOutputContainer ;
91 if(fCalibData) delete fCalibData;
92 if(fPhosGeo) delete fPhosGeo;
96 //__________________________________________________
97 void AliAnalysisTaskPHOSPi0CalibSelection::CreateAODFromAOD()
99 // Copy AOD header, vertex, CaloClusters and CaloCells to output AOD
101 AliAODEvent* aod = dynamic_cast<AliAODEvent*>(InputEvent());
104 printf("AliAnalysisTaskPHOSPi0CalibSelection::CreateAODFromAOD() - This event does not contain AODs?");
108 // set arrays and pointers
114 for (Int_t i = 0; i < 6; i++) covVtx[i] = 0.;
117 AliAODHeader*headerin = aod->GetHeader();
118 AliAODHeader* header = AODEvent()->GetHeader();
119 header->SetRunNumber(headerin->GetRunNumber());
120 header->SetBunchCrossNumber(headerin->GetBunchCrossNumber());
121 header->SetOrbitNumber(headerin->GetOrbitNumber());
122 header->SetPeriodNumber(headerin->GetPeriodNumber());
123 header->SetEventType(headerin->GetEventType());
124 header->SetMuonMagFieldScale(headerin->GetMuonMagFieldScale());
126 if(headerin->GetCentrality()){
127 header->SetCentrality(new AliCentrality(*(headerin->GetCentralityP())));
129 header->SetCentrality(0);
132 header->SetTriggerMask(headerin->GetTriggerMask());
133 header->SetTriggerCluster(headerin->GetTriggerCluster());
134 header->SetMagneticField(headerin->GetMagneticField());
135 header->SetZDCN1Energy(headerin->GetZDCN1Energy());
136 header->SetZDCP1Energy(headerin->GetZDCP1Energy());
137 header->SetZDCN2Energy(headerin->GetZDCN2Energy());
138 header->SetZDCP2Energy(headerin->GetZDCP2Energy());
139 header->SetZDCEMEnergy(headerin->GetZDCEMEnergy(0),headerin->GetZDCEMEnergy(1));
140 Float_t diamxy[2]={aod->GetDiamondX(),aod->GetDiamondY()};
141 Float_t diamcov[3]; aod->GetDiamondCovXY(diamcov);
142 header->SetDiamond(diamxy,diamcov);
145 Int_t nVertices = 1 ;/* = prim. vtx*/;
146 Int_t nCaloClus = aod->GetNumberOfCaloClusters();
148 AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
150 // Access to the AOD container of vertices
151 TClonesArray &vertices = *(AODEvent()->GetVertices());
154 // Add primary vertex. The primary tracks will be defined
155 // after the loops on the composite objects (V0, cascades, kinks)
156 const AliAODVertex *vtx = aod->GetPrimaryVertex();
158 vtx->GetXYZ(pos); // position
159 vtx->GetCovMatrix(covVtx); //covariance matrix
161 AliAODVertex * primary = new(vertices[jVertices++])
162 AliAODVertex(pos, covVtx, vtx->GetChi2perNDF(), NULL, -1, AliAODVertex::kPrimary);
163 primary->SetName(vtx->GetName());
164 primary->SetTitle(vtx->GetTitle());
166 // Access to the AOD container of clusters
167 TClonesArray &caloClusters = *(AODEvent()->GetCaloClusters());
170 for (Int_t iClust=0; iClust<nCaloClus; ++iClust) {
172 AliAODCaloCluster * cluster = aod->GetCaloCluster(iClust);
174 //Check if it is a PHOS cluster
175 if(!cluster->IsPHOS()) continue ;
177 Int_t id = cluster->GetID();
178 Float_t energy = cluster->E();
179 cluster->GetPosition(posF);
180 Char_t ttype = cluster->GetType();
182 AliAODCaloCluster *caloCluster = new(caloClusters[jClusters++])
183 AliAODCaloCluster(id,
191 caloCluster->SetCaloCluster(cluster->GetDistanceToBadChannel(),
192 cluster->GetDispersion(),
193 cluster->GetM20(), cluster->GetM02(),
194 cluster->GetEmcCpvDistance(),
195 cluster->GetNExMax(),cluster->GetTOF()) ;
197 caloCluster->SetPIDFromESD(cluster->GetPID());
198 caloCluster->SetNCells(cluster->GetNCells());
199 caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
200 caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
204 caloClusters.Expand(jClusters); // resize TObjArray to 'remove' slots for pseudo clusters
205 // end of loop on calo clusters
207 // fill PHOS cell info
208 if (aod->GetPHOSCells()) { // protection against missing AOD information
209 AliAODCaloCells &aodinPHcells = *(aod->GetPHOSCells());
210 Int_t nPHcell = aodinPHcells.GetNumberOfCells() ;
212 AliAODCaloCells &aodPHcells = *(AODEvent()->GetPHOSCells());
213 aodPHcells.CreateContainer(nPHcell);
214 aodPHcells.SetType(AliVCaloCells::kPHOSCell);
215 for (Int_t iCell = 0; iCell < nPHcell; iCell++) {
216 aodPHcells.SetCell(iCell,aodinPHcells.GetCellNumber(iCell),aodinPHcells.GetAmplitude(iCell));
224 //__________________________________________________
225 void AliAnalysisTaskPHOSPi0CalibSelection::CreateAODFromESD()
228 // Copy Header, Vertex, CaloClusters and CaloCells from ESDs to AODs
230 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
233 printf("AliAnalysisTaskPHOSPi0CalibSelection::CreateAODFromESD() - This event does not contain AODs?");
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(0); // 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->IsPHOS()) 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::kPHOSNeutral);
315 caloCluster->SetCaloCluster(cluster->GetDistanceToBadChannel(),
316 cluster->GetDispersion(),
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());
324 caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
327 caloClusters.Expand(jClusters); // resize TObjArray to 'remove' slots for pseudo clusters
328 // end of loop on calo clusters
330 // fill PHOS cell info
331 if( esd->GetPHOSCells()) { // protection against missing ESD information
332 AliESDCaloCells &esdPHcells = *(esd->GetPHOSCells());
333 Int_t nPHcell = esdPHcells.GetNumberOfCells() ;
335 AliAODCaloCells &aodPHcells = *(AODEvent()->GetPHOSCells());
336 aodPHcells.CreateContainer(nPHcell);
337 aodPHcells.SetType(AliVCaloCells::kPHOSCell);
338 for (Int_t iCell = 0; iCell < nPHcell; iCell++) {
339 aodPHcells.SetCell(iCell,esdPHcells.GetCellNumber(iCell),esdPHcells.GetAmplitude(iCell));
346 //__________________________________________________
347 void AliAnalysisTaskPHOSPi0CalibSelection::UserCreateOutputObjects()
349 //Create output container
350 fOutputContainer = new TList();
351 const Int_t buffersize = 255;
352 char hname[buffersize], htitl[buffersize];
354 for(Int_t iMod=0; iMod<5; iMod++) {
355 for(Int_t iX=0; iX<64; iX++) {
356 for(Int_t iZ=0; iZ<56; iZ++) {
357 snprintf(hname,buffersize,"%d_%d_%d",iMod,iX,iZ);
358 snprintf(htitl,buffersize,"Two-gamma inv. mass for mod %d, cell (%d,%d)",iMod,iX,iZ);
359 fHmpi0[iMod][iX][iZ] = new TH1F(hname,htitl,100,0.,300.);
360 fOutputContainer->Add(fHmpi0[iMod][iX][iZ]);
365 fHmgg = new TH1F("hmgg","2-cluster invariant mass",100,0.,300.);
366 fOutputContainer->Add(fHmgg);
368 fCalibData = new AliPHOSCalibData();
369 fPhosGeo = AliPHOSGeometry::GetInstance("IHEP") ;
373 //__________________________________________________
374 void AliAnalysisTaskPHOSPi0CalibSelection::UserExec(Option_t* /* option */)
376 //Analysis per event.
377 if(DebugLevel() > 1) printf("AliAnalysisTaskPHOSPi0CalibSelection <<< Event %d >>>\n",(Int_t)Entry());
379 AliAODEvent* aod = 0x0;
380 if(!strcmp(InputEvent()->GetName(),"AliAODEvent")) {
382 aod = dynamic_cast<AliAODEvent*>(InputEvent());
384 printf("AliAnalysisTaskEMCALPi0CalibSelection::UserExec() - This event does not contain AODs?");
388 // Create new AOD with only CaloClusters and CaloCells
389 if(fCopyAOD) CreateAODFromAOD();
391 else if(!strcmp(InputEvent()->GetName(),"AliESDEvent")) {
395 printf("AliAnalysisTaskPHOSPi0CalibSelection::UserExec() - This event does not contain AODs?");
399 // Create AOD with CaloClusters and use it as input.
400 // If filtering task is already executed, this is not needed.
401 if(fCopyAOD) CreateAODFromESD();
404 printf("AliAnalysisTaskPHOSPi0CalibSelection: Unknown event type, STOP!\n");
408 Double_t v[] = {aod->GetVertex(0)->GetX(),aod->GetVertex(0)->GetY(),aod->GetVertex(0)->GetZ()}; //to check!!
409 //aod->GetVertex()->GetXYZ(v) ;
410 TVector3 vtx(v); //Check
412 if(DebugLevel() > 1) printf("AliAnalysisTaskPHOSPi0CalibSelection Vertex: (%.3f,%.3f,%.3f)\n",vtx.X(),vtx.Y(),vtx.Z());
414 Int_t runNum = aod->GetRunNumber();
415 if(DebugLevel() > 1) printf("Run number: %d\n",runNum);
417 //Get the matrix with geometry information
418 //Still not implemented in AOD, just a workaround to be able to work at least with ESDs
419 if(!strcmp(InputEvent()->GetName(),"AliAODEvent")) {
421 printf("AliAnalysisTaskPHOSPi0CalibSelection Use ideal geometry, values geometry matrix not kept in AODs.\n");
424 if(DebugLevel() > 1) printf("AliAnalysisTaskPHOSPi0CalibSelection Load Misaligned matrices. \n");
425 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent()) ;
427 printf("AliAnalysisTaskPHOSPi0CalibSelection::UserExec() - This event does not contain ESDs?");
430 for(Int_t mod=0; mod<5; mod++){
431 if(esd->GetPHOSMatrix(mod))
432 fPhosGeo->SetMisalMatrix(esd->GetPHOSMatrix(mod),mod) ;
436 if(DebugLevel() > 1) printf("AliAnalysisTaskPHOSPi0CalibSelection Will use fLogWeight %.3f .\n",fLogWeight);
445 TRefArray * caloClustersArr = new TRefArray();
446 aod->GetPHOSClusters(caloClustersArr);
448 const Int_t kNumberOfPhosClusters = caloClustersArr->GetEntries() ;
449 if(DebugLevel() > 1) printf("AliAnalysisTaskPHOSPi0CalibSelection CaloClusters: %d\n", kNumberOfPhosClusters);
452 AliAODCaloCells *phsCells = aod->GetPHOSCells();
454 // loop over PHOS clusters
455 for(Int_t iClu=0; iClu<kNumberOfPhosClusters; iClu++) {
457 AliAODCaloCluster *c1 = (AliAODCaloCluster *) caloClustersArr->At(iClu);
458 if(!c1->IsPHOS()) continue; // EMCAL cluster!
460 Float_t e1i = c1->E(); // cluster energy before correction
461 if(e1i<fEmin) continue;
463 AliPHOSAodCluster clu1(*c1);
464 clu1.Recalibrate(fCalibData, phsCells);
465 clu1.EvalAll(fLogWeight,vtx);
466 clu1.EnergyCorrection() ;
468 clu1.GetMomentum(p1,v);
470 MaxEnergyCellPos(phsCells,&clu1,maxId);
471 fPhosGeo->AbsToRelNumbering(maxId, relid);
473 mod1 = relid[0]-1; // module
474 Int_t iX = relid[2]-1; // cluster X-coord
475 Int_t iZ = relid[3]-1; // cluster Z-coord
476 Float_t e1ii = clu1.E(); // cluster energy after correction
478 for (Int_t jClu=iClu; jClu<kNumberOfPhosClusters; jClu++) {
479 AliAODCaloCluster *c2 = (AliAODCaloCluster *) caloClustersArr->At(jClu);
480 if(!c2->IsPHOS()) continue; // EMCAL cluster!
481 if(c2->IsEqual(c1)) continue;
483 Float_t e2i = c2->E();
484 if(e2i<fEmin) continue;
486 AliPHOSAodCluster clu2(*c2);
487 clu2.Recalibrate(fCalibData, phsCells);
488 clu2.EvalAll(fLogWeight,vtx);
489 clu2.EnergyCorrection() ;
491 clu2.GetMomentum(p2,v);
492 Float_t e2ii = clu2.E();
496 fHmgg->Fill(p12.M()*1000);
497 fHmpi0[mod1][iX][iZ]->Fill(p12.M()*1000);
499 if(DebugLevel() > 1) printf("AliAnalysisTaskPHOSPi0CalibSelection Mass in (mod%d,%d,%d): %.3f GeV E1_i=%f E1_ii=%f E2_i=%f E2_ii=%f\n",
500 mod1,iX,iZ,p12.M(),e1i,e1ii,e2i,e2ii);
503 } // end of loop over PHOS clusters
505 delete caloClustersArr;
506 PostData(1,fOutputContainer);
509 //__________________________________________________
510 void AliAnalysisTaskPHOSPi0CalibSelection::MaxEnergyCellPos(AliAODCaloCells *cells, AliAODCaloCluster* clu, Int_t& maxId)
512 //For a given CaloCluster calculates the absId of the cell
513 //with maximum energy deposit.
515 Double_t eMax = -111;
517 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
518 Int_t cellAbsId = clu->GetCellAbsId(iDig);
519 Double_t eCell = cells->GetCellAmplitude(cellAbsId)*clu->GetCellAmplitudeFraction(iDig);
528 //__________________________________________________
529 void AliAnalysisTaskPHOSPi0CalibSelection::SetCalibCorrections(AliPHOSCalibData* cdata)
531 //Set new correction factors (~1) to calibration coefficients, delete previous.
533 if(fCalibData) delete fCalibData;