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 "TRefArray.h"
32 #include <TGeoManager.h>
35 #include "AliAnalysisTaskEMCALPi0CalibSelection.h"
36 #include "AliAODEvent.h"
37 #include "AliESDEvent.h"
38 #include "AliEMCALGeometry.h"
39 #include "AliVCluster.h"
40 #include "AliVCaloCells.h"
41 #include "AliEMCALRecoUtils.h"
42 //#include "AliEMCALAodCluster.h"
43 //#include "AliEMCALCalibData.h"
45 ClassImp(AliAnalysisTaskEMCALPi0CalibSelection)
48 //__________________________________________________
49 AliAnalysisTaskEMCALPi0CalibSelection::AliAnalysisTaskEMCALPi0CalibSelection(const char* name) :
50 AliAnalysisTaskSE(name),fEMCALGeo(0x0),//fCalibData(0x0),
51 fEmin(0.5), fEmax(15.), fAsyCut(1.),fMinNCells(2), fGroupNCells(0),
52 fLogWeight(4.5), fSameSM(kFALSE), fOldAOD(kFALSE), fFilteredInput(kFALSE),
53 fCorrectClusters(kFALSE), fEMCALGeoName("EMCAL_FIRSTYEARV1"),
54 fRecoUtils(new AliEMCALRecoUtils),
55 fNbins(300), fMinBin(0.), fMaxBin(300.),fOutputContainer(0x0),
56 fHmgg(0x0), fHmggDifferentSM(0x0),
57 fHOpeningAngle(0x0), fHOpeningAngleDifferentSM(0x0),
58 fHIncidentAngle(0x0), fHIncidentAngleDifferentSM(0x0),
59 fHAsymmetry(0x0), fHAsymmetryDifferentSM(0x0),
60 fhNEvents(0x0),fCuts(0x0),fLoadMatrices(0)
62 //Named constructor which should be used.
64 for(Int_t iMod=0; iMod < AliEMCALGeoParams::fgkEMCALModules; iMod++) {
65 for(Int_t iX=0; iX<24; iX++) {
66 for(Int_t iZ=0; iZ<48; iZ++) {
67 fHmpi0[iMod][iZ][iX]=0;
72 for(Int_t iSM = 0; iSM < AliEMCALGeoParams::fgkEMCALModules; iSM++) {
75 fHOpeningAngleSM[iSM] = 0;
76 fHOpeningAnglePairSM[iSM] = 0;
77 fHAsymmetrySM[iSM] = 0;
78 fHAsymmetryPairSM[iSM] = 0;
79 fHIncidentAngleSM[iSM] = 0;
80 fHIncidentAnglePairSM[iSM] = 0;
81 fhTowerDecayPhotonHit[iSM] = 0;
82 fhTowerDecayPhotonEnergy[iSM] = 0;
83 fhTowerDecayPhotonAsymmetry[iSM] = 0;
87 DefineOutput(1, TList::Class());
88 DefineOutput(2, TList::Class()); // will contain cuts or local params
92 //__________________________________________________
93 AliAnalysisTaskEMCALPi0CalibSelection::~AliAnalysisTaskEMCALPi0CalibSelection()
98 fOutputContainer->Delete() ;
99 delete fOutputContainer ;
102 //if(fCalibData) delete fCalibData;
103 if(fEMCALGeo) delete fEMCALGeo ;
104 if(fRecoUtils) delete fRecoUtils ;
108 //_____________________________________________________
109 void AliAnalysisTaskEMCALPi0CalibSelection::LocalInit()
111 // Local Initialization
113 // Create cuts/param objects and publish to slot
114 const Int_t buffersize = 255;
115 char onePar[buffersize] ;
118 snprintf(onePar,buffersize, "Custer cuts: %2.2f < E < %2.2f GeV; min number of cells %d; Assymetry cut %1.2f", fEmin,fEmax, fMinNCells, fAsyCut) ;
119 fCuts->Add(new TObjString(onePar));
120 snprintf(onePar,buffersize, "Group %d cells;", fGroupNCells) ;
121 fCuts->Add(new TObjString(onePar));
122 snprintf(onePar,buffersize, "Cluster maximal cell away from border at least %d cells;", fRecoUtils->GetNumberOfCellsFromEMCALBorder()) ;
123 fCuts->Add(new TObjString(onePar));
124 snprintf(onePar,buffersize, "Histograms: bins %d; energy range: %2.2f < E < %2.2f GeV;",fNbins,fMinBin,fMaxBin) ;
125 fCuts->Add(new TObjString(onePar));
126 snprintf(onePar,buffersize, "Switchs: Remove Bad Channels? %d; Use filtered input? %d; Correct Clusters? %d, Analyze Old AODs? %d, Mass per channel same SM clusters? %d ",
127 fRecoUtils->IsBadChannelsRemovalSwitchedOn(),fFilteredInput,fCorrectClusters, fOldAOD, fSameSM) ;
128 fCuts->Add(new TObjString(onePar));
129 snprintf(onePar,buffersize, "EMCAL Geometry name: < %s >, Load Matrices? %d",fEMCALGeoName.Data(),fLoadMatrices) ;
130 fCuts->Add(new TObjString(onePar));
137 //_________________________________________________________________
138 Int_t AliAnalysisTaskEMCALPi0CalibSelection::GetEMCALClusters(AliVEvent * event, TRefArray *clusters) const
140 // fills the provided TRefArray with all found emcal clusters
144 Bool_t first = kTRUE;
145 for (Int_t i = 0; i < event->GetNumberOfCaloClusters(); i++) {
146 if ( (cl = event->GetCaloCluster(i)) ) {
147 if (IsEMCALCluster(cl)){
149 new (clusters) TRefArray(TProcessID::GetProcessWithUID(cl));
153 //printf("IsEMCal cluster %d, E %2.3f Size: %d \n",i,cl->E(),clusters->GetEntriesFast());
157 return clusters->GetEntriesFast();
161 //____________________________________________________________________________
162 Bool_t AliAnalysisTaskEMCALPi0CalibSelection::IsEMCALCluster(AliVCluster* cluster) const {
163 // Check if it is a cluster from EMCAL. For old AODs cluster type has
164 // different number and need to patch here
168 if (cluster->GetType() == 2) return kTRUE;
173 return cluster->IsEMCAL();
179 //__________________________________________________
180 void AliAnalysisTaskEMCALPi0CalibSelection::UserCreateOutputObjects()
182 //Create output container, init geometry
184 fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName) ;
185 Int_t nSM = (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules();
187 fOutputContainer = new TList();
188 const Int_t buffersize = 255;
189 char hname[buffersize], htitl[buffersize];
191 for(Int_t iMod=0; iMod < nSM; iMod++) {
192 for(Int_t iRow=0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
193 for(Int_t iCol=0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
194 snprintf(hname,buffersize, "%d_%d_%d",iMod,iCol,iRow);
195 snprintf(htitl,buffersize, "Two-gamma inv. mass for super mod %d, cell(col,row)=(%d,%d)",iMod,iCol,iRow);
196 fHmpi0[iMod][iCol][iRow] = new TH1F(hname,htitl,fNbins,fMinBin,fMaxBin);
197 fOutputContainer->Add(fHmpi0[iMod][iCol][iRow]);
202 fHmgg = new TH2F("hmgg","2-cluster invariant mass",fNbins,fMinBin,fMaxBin,100,0,10);
203 fHmgg->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
204 fHmgg->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
205 fOutputContainer->Add(fHmgg);
207 fHmggDifferentSM = new TH2F("hmggDifferentSM","2-cluster invariant mass, different SM",fNbins,fMinBin,fMaxBin,100,0,10);
208 fHmggDifferentSM->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
209 fHmggDifferentSM->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
210 fOutputContainer->Add(fHmggDifferentSM);
212 fHOpeningAngle = new TH2F("hopang","2-cluster opening angle",100,0.,50.,100,0,10);
213 fHOpeningAngle->SetXTitle("#alpha_{#gamma #gamma}");
214 fHOpeningAngle->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
215 fOutputContainer->Add(fHOpeningAngle);
217 fHOpeningAngleDifferentSM = new TH2F("hopangDifferentSM","2-cluster opening angle, different SM",100,0,50.,100,0,10);
218 fHOpeningAngleDifferentSM->SetXTitle("#alpha_{#gamma #gamma}");
219 fHOpeningAngleDifferentSM->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
220 fOutputContainer->Add(fHOpeningAngleDifferentSM);
222 fHIncidentAngle = new TH2F("hinang","#gamma incident angle in SM",100,0.,20.,100,0,10);
223 fHIncidentAngle->SetXTitle("#alpha_{#gamma SM center}");
224 fHIncidentAngle->SetYTitle("p_{T #gamma} (GeV/c)");
225 fOutputContainer->Add(fHIncidentAngle);
227 fHIncidentAngleDifferentSM = new TH2F("hinangDifferentSM","#gamma incident angle in SM, different SM pair",100,0,20.,100,0,10);
228 fHIncidentAngleDifferentSM->SetXTitle("#alpha_{#gamma - SM center}");
229 fHIncidentAngleDifferentSM->SetYTitle("p_{T #gamma} (GeV/c)");
230 fOutputContainer->Add(fHIncidentAngleDifferentSM);
232 fHAsymmetry = new TH2F("hasym","2-cluster opening angle",100,0.,1.,100,0,10);
233 fHAsymmetry->SetXTitle("a");
234 fHAsymmetry->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
235 fOutputContainer->Add(fHAsymmetry);
237 fHAsymmetryDifferentSM = new TH2F("hasymDifferentSM","2-cluster opening angle, different SM",100,0,1.,100,0,10);
238 fHAsymmetryDifferentSM->SetXTitle("a");
239 fHAsymmetryDifferentSM->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
240 fOutputContainer->Add(fHAsymmetryDifferentSM);
243 TString pairname[] = {"A side (0-2)", "C side (1-3)","Row 0 (0-1)", "Row 1 (2-3)"};
245 for(Int_t iSM = 0; iSM < nSM; iSM++) {
247 snprintf(hname, buffersize, "hmgg_SM%d",iSM);
248 snprintf(htitl, buffersize, "Two-gamma inv. mass for super mod %d",iSM);
249 fHmggSM[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
250 fHmggSM[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
251 fHmggSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
252 fOutputContainer->Add(fHmggSM[iSM]);
254 snprintf(hname,buffersize, "hmgg_PairSM%d",iSM);
255 snprintf(htitl,buffersize, "Two-gamma inv. mass for SM pair: %s",pairname[iSM].Data());
256 fHmggPairSM[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
257 fHmggPairSM[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
258 fHmggPairSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
259 fOutputContainer->Add(fHmggPairSM[iSM]);
262 snprintf(hname, buffersize, "hopang_SM%d",iSM);
263 snprintf(htitl, buffersize, "Opening angle for super mod %d",iSM);
264 fHOpeningAngleSM[iSM] = new TH2F(hname,htitl,100,0.,50.,100,0,10);
265 fHOpeningAngleSM[iSM]->SetXTitle("#alpha_{#gamma #gamma} (deg)");
266 fHOpeningAngleSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
267 fOutputContainer->Add(fHOpeningAngleSM[iSM]);
269 snprintf(hname,buffersize, "hopang_PairSM%d",iSM);
270 snprintf(htitl,buffersize, "Opening angle for SM pair: %s",pairname[iSM].Data());
271 fHOpeningAnglePairSM[iSM] = new TH2F(hname,htitl,100,0.,50.,100,0,10);
272 fHOpeningAnglePairSM[iSM]->SetXTitle("#alpha_{#gamma #gamma} (deg)");
273 fHOpeningAnglePairSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
274 fOutputContainer->Add(fHOpeningAnglePairSM[iSM]);
276 snprintf(hname, buffersize, "hinang_SM%d",iSM);
277 snprintf(htitl, buffersize, "Incident angle for super mod %d",iSM);
278 fHIncidentAngleSM[iSM] = new TH2F(hname,htitl,100,0.,20.,100,0,10);
279 fHIncidentAngleSM[iSM]->SetXTitle("#alpha_{#gamma - SM center} (deg)");
280 fHIncidentAngleSM[iSM]->SetYTitle("p_{T #gamma} (GeV/c)");
281 fOutputContainer->Add(fHIncidentAngleSM[iSM]);
283 snprintf(hname,buffersize, "hinang_PairSM%d",iSM);
284 snprintf(htitl,buffersize, "Incident angle for SM pair: %s",pairname[iSM].Data());
285 fHIncidentAnglePairSM[iSM] = new TH2F(hname,htitl,100,0.,20.,100,0,10);
286 fHIncidentAnglePairSM[iSM]->SetXTitle("#alpha_{#gamma - SM center} (deg)");
287 fHIncidentAnglePairSM[iSM]->SetYTitle("p_{T #gamma} (GeV/c)");
288 fOutputContainer->Add(fHIncidentAnglePairSM[iSM]);
290 snprintf(hname, buffersize, "hasym_SM%d",iSM);
291 snprintf(htitl, buffersize, "asymmetry for super mod %d",iSM);
292 fHAsymmetrySM[iSM] = new TH2F(hname,htitl,100,0.,1.,100,0,10);
293 fHAsymmetrySM[iSM]->SetXTitle("a");
294 fHAsymmetrySM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
295 fOutputContainer->Add(fHAsymmetrySM[iSM]);
297 snprintf(hname,buffersize, "hasym_PairSM%d",iSM);
298 snprintf(htitl,buffersize, "Asymmetry for SM pair: %s",pairname[iSM].Data());
299 fHAsymmetryPairSM[iSM] = new TH2F(hname,htitl,100,0.,1.,100,0,10);
300 fHAsymmetryPairSM[iSM]->SetXTitle("a");
301 fHAsymmetryPairSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
302 fOutputContainer->Add(fHAsymmetryPairSM[iSM]);
308 fhTowerDecayPhotonHit[iSM] = new TH2F (Form("hTowerDecPhotonHit_Mod%d",iSM),Form("Entries in grid of cells in Module %d",iSM),
309 colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5);
310 fhTowerDecayPhotonHit[iSM]->SetYTitle("row (phi direction)");
311 fhTowerDecayPhotonHit[iSM]->SetXTitle("column (eta direction)");
312 fOutputContainer->Add(fhTowerDecayPhotonHit[iSM]);
314 fhTowerDecayPhotonEnergy[iSM] = new TH2F (Form("hTowerDecPhotonEnergy_Mod%d",iSM),Form("Accumulated energy in grid of cells in Module %d",iSM),
315 colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5);
316 fhTowerDecayPhotonEnergy[iSM]->SetYTitle("row (phi direction)");
317 fhTowerDecayPhotonEnergy[iSM]->SetXTitle("column (eta direction)");
318 fOutputContainer->Add(fhTowerDecayPhotonEnergy[iSM]);
320 fhTowerDecayPhotonAsymmetry[iSM] = new TH2F (Form("hTowerDecPhotonAsymmetry_Mod%d",iSM),Form("Accumulated asymmetry in grid of cells in Module %d",iSM),
321 colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5);
322 fhTowerDecayPhotonAsymmetry[iSM]->SetYTitle("row (phi direction)");
323 fhTowerDecayPhotonAsymmetry[iSM]->SetXTitle("column (eta direction)");
324 fOutputContainer->Add(fhTowerDecayPhotonAsymmetry[iSM]);
328 fhNEvents = new TH1I("hNEvents", "Number of analyzed events" , 1 , 0 , 1 ) ;
329 fOutputContainer->Add(fhNEvents);
331 fOutputContainer->SetOwner(kTRUE);
333 // fCalibData = new AliEMCALCalibData();
335 PostData(1,fOutputContainer);
339 //__________________________________________________
340 void AliAnalysisTaskEMCALPi0CalibSelection::UserExec(Option_t* /* option */)
342 //Analysis per event.
344 if(fRecoUtils->GetParticleType()!=AliEMCALRecoUtils::kPhoton){
345 printf("Wrong particle type for cluster position recalculation! = %d\n", fRecoUtils->GetParticleType());
349 fhNEvents->Fill(0); //Event analyzed
351 //Get the input event
352 AliVEvent* event = 0;
353 if(fFilteredInput) event = AODEvent();
354 else event = InputEvent();
357 printf("Input event not available!\n");
362 printf("AliAnalysisTaskEMCALPi0CalibSelection <<< %s: Event %d >>>\n",event->GetName(), (Int_t)Entry());
365 //Get the primary vertex
367 event->GetPrimaryVertex()->GetXYZ(v) ;
369 if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Vertex: (%.3f,%.3f,%.3f)\n",v[0],v[1],v[2]);
371 //Int_t runNum = aod->GetRunNumber();
372 //if(DebugLevel() > 1) printf("Run number: %d\n",runNum);
374 //Get the matrix with geometry information
375 if(fhNEvents->GetEntries()==1){
377 printf("AliAnalysisTaskEMCALPi0CalibSelection::UserExec() - Load user defined geometry matrices\n");
378 for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
381 fMatrix[mod]->Print();
382 fEMCALGeo->SetMisalMatrix(fMatrix[mod],mod) ;
386 else if(!gGeoManager){
387 printf("AliAnalysisTaskEMCALPi0CalibSelection::UserExec() - Get geo matrices from data\n");
388 //Still not implemented in AOD, just a workaround to be able to work at least with ESDs
389 if(!strcmp(event->GetName(),"AliAODEvent")) {
391 printf("AliAnalysisTaskEMCALPi0CalibSelection Use ideal geometry, values geometry matrix not kept in AODs.\n");
394 if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Load Misaligned matrices. \n");
395 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event) ;
397 printf("AliAnalysisTaskEMCALPi0CalibSelection::UserExec() - This event does not contain ESDs?");
400 for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
401 //if(DebugLevel() > 1)
402 esd->GetEMCALMatrix(mod)->Print();
403 if(esd->GetEMCALMatrix(mod)) fEMCALGeo->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
406 }//Load matrices from Data
409 if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Will use fLogWeight %.3f .\n",fLogWeight);
418 Bool_t shared = kFALSE;
424 //Get the list of clusters
425 TRefArray * caloClustersArr = new TRefArray();
426 if(!fOldAOD) event->GetEMCALClusters(caloClustersArr);
427 else GetEMCALClusters(event,caloClustersArr);
428 const Int_t kNumberOfEMCALClusters = caloClustersArr->GetEntries() ;
429 if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection - N CaloClusters: %d \n", kNumberOfEMCALClusters);
432 AliVCaloCells *emCells = event->GetEMCALCells();
434 // loop over EMCAL clusters
435 //----------------------------------------------------------
436 // First recalibrate and recalculate energy and position
437 Float_t pos[]={0,0,0};
438 if(fCorrectClusters){
439 for(Int_t iClu=0; iClu<kNumberOfEMCALClusters-1; iClu++) {
440 AliVCluster *c1 = (AliVCluster *) caloClustersArr->At(iClu);
442 if(fRecoUtils->ClusterContainsBadChannel(fEMCALGeo, c1->GetCellsAbsId(), c1->GetNCells())) continue;
446 printf("Std : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",c1->GetID(),c1->E(),c1->GetDispersion(),c1->GetM02(),c1->GetM20());
447 c1->GetPosition(pos);
448 printf("Std : i %d, x %f, y %f, z %f\n",c1->GetID(), pos[0], pos[1], pos[2]);
451 //Correct cluster energy and position if requested, and not corrected previously, by default Off
452 if(fRecoUtils->IsRecalibrationOn()) {
453 fRecoUtils->RecalibrateClusterEnergy(fEMCALGeo, c1, emCells);
454 fRecoUtils->RecalculateClusterShowerShapeParameters(fEMCALGeo, emCells,c1);
455 fRecoUtils->RecalculateClusterPID(c1);
458 printf("Energy: after recalibration %f; ",c1->E());
460 // Recalculate cluster position
461 fRecoUtils->RecalculateClusterPosition(fEMCALGeo, emCells,c1);
463 // Correct Non-Linearity
464 c1->SetE(fRecoUtils->CorrectClusterEnergyLinearity(c1));
466 printf("after linearity correction %f\n",c1->E());
470 printf("Cor : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",c1->GetID(),c1->E(),c1->GetDispersion(),c1->GetM02(),c1->GetM20());
471 c1->GetPosition(pos);
472 printf("Cor : i %d, x %f, y %f, z %f\n",c1->GetID(), pos[0], pos[1], pos[2]);
477 //----------------------------------------------------------
478 //Now the invariant mass analysis with the corrected clusters
479 for(Int_t iClu=0; iClu<kNumberOfEMCALClusters-1; iClu++) {
481 AliVCluster *c1 = (AliVCluster *) caloClustersArr->At(iClu);
482 if(fRecoUtils->ClusterContainsBadChannel(fEMCALGeo, c1->GetCellsAbsId(), c1->GetNCells())) continue;
484 Float_t e1i = c1->E(); // cluster energy before correction
485 if(e1i < fEmin) continue;
486 else if(e1i > fEmax) continue;
487 else if (c1->GetNCells() < fMinNCells) continue;
491 printf("IMA : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",c1->GetID(),e1i,c1->GetDispersion(),c1->GetM02(),c1->GetM20());
492 c1->GetPosition(pos);
493 printf("IMA : i %d, x %f, y %f, z %f\n",c1->GetID(), pos[0], pos[1], pos[2]);
496 //AliEMCALAodCluster newc1(*((AliAODCaloCluster*)c1));
497 //newc1.EvalAllFromRecoUtils(fEMCALGeo,fRecoUtils,emCells);
498 //printf("i %d, recal? %d\n",iClu,newc1.IsRecalibrated());
499 //clu1.Recalibrate(fCalibData, emCells, fEMCALGeoName);
501 //clu1.EvalAll(fLogWeight, fEMCALGeoName);
503 fRecoUtils->GetMaxEnergyCell(fEMCALGeo, emCells,c1,absId1,iSupMod1,ieta1,iphi1,shared);
504 c1->GetMomentum(p1,v);
505 //newc1.GetMomentum(p1,v);
507 // Combine cluster with other clusters and get the invariant mass
508 for (Int_t jClu=iClu+1; jClu<kNumberOfEMCALClusters; jClu++) {
509 AliAODCaloCluster *c2 = (AliAODCaloCluster *) caloClustersArr->At(jClu);
510 //if(c2->IsEqual(c1)) continue;
511 if(fRecoUtils->ClusterContainsBadChannel(fEMCALGeo, c2->GetCellsAbsId(), c2->GetNCells())) continue;
513 Float_t e2i = c2->E();
514 if(e2i < fEmin) continue;
515 else if (e2i > fEmax) continue;
516 else if (c2->GetNCells() < fMinNCells) continue;
518 //AliEMCALAodCluster newc2(*((AliAODCaloCluster*)c2));
519 //newc2.EvalAllFromRecoUtils(fEMCALGeo,fRecoUtils,emCells);
520 //printf("\t j %d, recal? %d\n",jClu,newc2.IsRecalibrated());
521 //clu2.Recalibrate(fCalibData, emCells,fEMCALGeoName);
523 //clu2.EvalAll(fLogWeight,fEMCALGeoName);
525 fRecoUtils->GetMaxEnergyCell(fEMCALGeo, emCells,c2,absId2,iSupMod2,ieta2,iphi2,shared);
526 c2->GetMomentum(p2,v);
527 //newc2.GetMomentum(p2,v);
529 Float_t invmass = p12.M()*1000;
530 //printf("*** mass %f\n",invmass);
531 Float_t asym = TMath::Abs(p1.E()-p2.E())/(p1.E()+p2.E());
532 //printf("asymmetry %f\n",asym);
534 if(asym > fAsyCut) continue;
536 if(invmass < fMaxBin && invmass > fMinBin){
538 //Check if cluster is in fidutial region, not too close to borders
539 Bool_t in1 = fRecoUtils->CheckCellFiducialRegion(fEMCALGeo, c1, emCells);
540 Bool_t in2 = fRecoUtils->CheckCellFiducialRegion(fEMCALGeo, c2, emCells);
544 fHmgg->Fill(invmass,p12.Pt());
546 if(iSupMod1==iSupMod2) fHmggSM[iSupMod1]->Fill(invmass,p12.Pt());
547 else fHmggDifferentSM ->Fill(invmass,p12.Pt());
549 if((iSupMod1==0 && iSupMod2==2) || (iSupMod1==2 && iSupMod2==0)) fHmggPairSM[0]->Fill(invmass,p12.Pt());
550 if((iSupMod1==1 && iSupMod2==3) || (iSupMod1==3 && iSupMod2==1)) fHmggPairSM[1]->Fill(invmass,p12.Pt());
551 if((iSupMod1==0 && iSupMod2==1) || (iSupMod1==1 && iSupMod2==0)) fHmggPairSM[2]->Fill(invmass,p12.Pt());
552 if((iSupMod1==2 && iSupMod2==3) || (iSupMod1==3 && iSupMod2==2)) fHmggPairSM[3]->Fill(invmass,p12.Pt());
554 if(invmass > 100. && invmass < 160.){//restrict to clusters really close to pi0 peak
556 //Opening angle of 2 photons
557 Float_t opangle = p1.Angle(p2.Vect())*TMath::RadToDeg();
558 //printf("*******>>>>>>>> In PEAK pt %f, angle %f \n",p12.Pt(),opangle);
560 //Incident angle of each photon
561 Float_t inangle1 =0., inangle2=0.;
563 Float_t posSM1cen[3]={0.,0.,0.};
564 Float_t depth = fRecoUtils->GetDepth(p1.Energy(),fRecoUtils->GetParticleType(),iSupMod1);
565 fEMCALGeo->RecalculateTowerPosition(11.5, 23.5, iSupMod1, depth, fRecoUtils->GetMisalTransShiftArray(),fRecoUtils->GetMisalRotShiftArray(),posSM1cen);
566 Float_t posSM2cen[3]={0.,0.,0.};
567 depth = fRecoUtils->GetDepth(p2.Energy(),fRecoUtils->GetParticleType(),iSupMod2);
568 fEMCALGeo->RecalculateTowerPosition(11.5, 23.5, iSupMod2, depth, fRecoUtils->GetMisalTransShiftArray(),fRecoUtils->GetMisalRotShiftArray(),posSM2cen);
569 //printf("SM1 %d pos (%2.3f,%2.3f,%2.3f) \n",iSupMod1,posSM1cen[0],posSM1cen[1],posSM1cen[2]);
570 //printf("SM2 %d pos (%2.3f,%2.3f,%2.3f) \n",iSupMod2,posSM2cen[0],posSM2cen[1],posSM2cen[2]);
572 TVector3 vecSM1cen(posSM1cen[0]-v[0],posSM1cen[1]-v[1],posSM1cen[2]-v[2]);
573 TVector3 vecSM2cen(posSM2cen[0]-v[0],posSM2cen[1]-v[1],posSM2cen[2]-v[2]);
574 inangle1 = p1.Angle(vecSM1cen)*TMath::RadToDeg();
575 inangle2 = p2.Angle(vecSM2cen)*TMath::RadToDeg();
576 //printf("Incident angle: cluster 1 %2.3f; cluster 2 %2.3f\n",inangle1,inangle2);
578 fHOpeningAngle ->Fill(opangle,p12.Pt());
579 fHIncidentAngle->Fill(inangle1,p1.Pt());
580 fHIncidentAngle->Fill(inangle2,p2.Pt());
581 fHAsymmetry ->Fill(asym,p12.Pt());
583 if(iSupMod1==iSupMod2) {
584 fHOpeningAngleSM[iSupMod1] ->Fill(opangle,p12.Pt());
585 fHIncidentAngleSM[iSupMod1]->Fill(inangle1,p1.Pt());
586 fHIncidentAngleSM[iSupMod1]->Fill(inangle2,p2.Pt());
587 fHAsymmetrySM[iSupMod1] ->Fill(asym,p12.Pt());
590 fHOpeningAngleDifferentSM ->Fill(opangle,p12.Pt());
591 fHIncidentAngleDifferentSM ->Fill(inangle1,p1.Pt());
592 fHIncidentAngleDifferentSM ->Fill(inangle2,p2.Pt());
593 fHAsymmetryDifferentSM ->Fill(asym,p12.Pt());
596 if((iSupMod1==0 && iSupMod2==2) || (iSupMod1==2 && iSupMod2==0)) {
597 fHOpeningAnglePairSM[0] ->Fill(opangle,p12.Pt());
598 fHIncidentAnglePairSM[0]->Fill(inangle1,p1.Pt());
599 fHIncidentAnglePairSM[0]->Fill(inangle2,p2.Pt());
600 fHAsymmetryPairSM[0] ->Fill(asym,p12.Pt());
603 if((iSupMod1==1 && iSupMod2==3) || (iSupMod1==3 && iSupMod2==1)) {
604 fHOpeningAnglePairSM[1] ->Fill(opangle,p12.Pt());
605 fHIncidentAnglePairSM[1]->Fill(inangle1,p1.Pt());
606 fHIncidentAnglePairSM[1]->Fill(inangle2,p2.Pt());
607 fHAsymmetryPairSM[1] ->Fill(asym,p12.Pt());
611 if((iSupMod1==0 && iSupMod2==1) || (iSupMod1==1 && iSupMod2==0)) {
612 fHOpeningAnglePairSM[2] ->Fill(opangle,p12.Pt());
613 fHIncidentAnglePairSM[2]->Fill(inangle1,p1.Pt());
614 fHIncidentAnglePairSM[2]->Fill(inangle2,p2.Pt());
615 fHAsymmetryPairSM[2] ->Fill(asym,p12.Pt());
619 if((iSupMod1==2 && iSupMod2==3) || (iSupMod1==3 && iSupMod2==2)) {
620 fHOpeningAnglePairSM[3] ->Fill(opangle,p12.Pt());
621 fHIncidentAnglePairSM[3]->Fill(inangle1,p1.Pt());
622 fHIncidentAnglePairSM[3]->Fill(inangle2,p2.Pt());
623 fHAsymmetryPairSM[3] ->Fill(asym,p12.Pt());
626 }// pair in 100 < mass < 160
628 }//in acceptance cuts
630 //In case of filling only channels with second cluster in same SM
631 if(fSameSM && iSupMod1!=iSupMod2) continue;
633 if (fGroupNCells == 0){
634 fHmpi0[iSupMod1][ieta1][iphi1]->Fill(invmass);
635 fHmpi0[iSupMod2][ieta2][iphi2]->Fill(invmass);
637 if(invmass > 100. && invmass < 160.){//restrict to clusters really close to pi0 peak
638 fhTowerDecayPhotonHit [iSupMod1]->Fill(ieta1,iphi1);
639 fhTowerDecayPhotonEnergy [iSupMod1]->Fill(ieta1,iphi1,p1.E());
640 fhTowerDecayPhotonAsymmetry[iSupMod1]->Fill(ieta1,iphi1,asym);
642 fhTowerDecayPhotonHit [iSupMod2]->Fill(ieta2,iphi2);
643 fhTowerDecayPhotonEnergy [iSupMod2]->Fill(ieta2,iphi2,p2.E());
644 fhTowerDecayPhotonAsymmetry[iSupMod2]->Fill(ieta2,iphi2,asym);
646 }// pair in mass of pi0
649 //printf("Regroup N %d, eta1 %d, phi1 %d, eta2 %d, phi2 %d \n",fGroupNCells, ieta1, iphi1, ieta2, iphi2);
650 for (Int_t i = -fGroupNCells; i < fGroupNCells+1; i++) {
651 for (Int_t j = -fGroupNCells; j < fGroupNCells+1; j++) {
652 //printf("\t i %d, j %d\n",i,j);
653 if((ieta1+i >= 0) && (iphi1+j >= 0) && (ieta1+i < 48) && (iphi1+j < 24)){
654 //printf("\t \t eta1+i %d, phi1+j %d\n", ieta1+i, iphi1+j);
655 fHmpi0[iSupMod1][ieta1+i][iphi1+j]->Fill(invmass);
657 if((ieta2+i >= 0) && (iphi2+j >= 0) && (ieta2+i < 48) && (iphi2+j < 24)){
658 //printf("\t \t eta2+i %d, phi2+j %d\n", ieta2+i, iphi2+j);
659 fHmpi0[iSupMod2][ieta2+i][iphi2+j]->Fill(invmass);
665 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",
666 iSupMod1,iphi1,ieta1,iSupMod2,iphi2,ieta2,p12.M(),e1i,c1->E(),e2i,c2->E());
671 } // end of loop over EMCAL clusters
673 delete caloClustersArr;
675 PostData(1,fOutputContainer);
679 //_____________________________________________________
680 void AliAnalysisTaskEMCALPi0CalibSelection::PrintInfo(){
683 printf("Cluster cuts: %2.2f < E < %2.2f GeV; min number of cells %d; Assymetry cut %1.2f\n", fEmin,fEmax, fMinNCells, fAsyCut) ;
684 printf("Group %d cells\n", fGroupNCells) ;
685 printf("Cluster maximal cell away from border at least %d cells\n", fRecoUtils->GetNumberOfCellsFromEMCALBorder()) ;
686 printf("Histograms: bins %d; energy range: %2.2f < E < %2.2f GeV\n",fNbins,fMinBin,fMaxBin) ;
687 printf("Switchs:\n \t Remove Bad Channels? %d; Use filtered input? %d; Correct Clusters? %d, \n \t Analyze Old AODs? %d, Mass per channel same SM clusters? %d\n",
688 fRecoUtils->IsBadChannelsRemovalSwitchedOn(),fFilteredInput,fCorrectClusters, fOldAOD, fSameSM) ;
689 printf("EMCAL Geometry name: < %s >, Load Matrices %d\n",fEMCALGeoName.Data(), fLoadMatrices) ;
690 if(fLoadMatrices) {for(Int_t ism = 0; ism < AliEMCALGeoParams::fgkEMCALModules; ism++) fMatrix[ism]->Print();}
695 //__________________________________________________
696 //void AliAnalysisTaskEMCALPi0CalibSelection::SetCalibCorrections(AliEMCALCalibData* const cdata)
698 // //Set new correction factors (~1) to calibration coefficients, delete previous.
700 // if(fCalibData) delete fCalibData;
701 // fCalibData = cdata;