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 < 12; 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<4; 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) ;
186 fOutputContainer = new TList();
187 const Int_t buffersize = 255;
188 char hname[buffersize], htitl[buffersize];
190 for(Int_t iMod=0; iMod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); iMod++) {
191 for(Int_t iRow=0; iRow<AliEMCALGeoParams::fgkEMCALRows; iRow++) {
192 for(Int_t iCol=0; iCol<AliEMCALGeoParams::fgkEMCALCols; iCol++) {
193 snprintf(hname,buffersize, "%d_%d_%d",iMod,iCol,iRow);
194 snprintf(htitl,buffersize, "Two-gamma inv. mass for super mod %d, cell(col,row)=(%d,%d)",iMod,iCol,iRow);
195 fHmpi0[iMod][iCol][iRow] = new TH1F(hname,htitl,fNbins,fMinBin,fMaxBin);
196 fOutputContainer->Add(fHmpi0[iMod][iCol][iRow]);
201 fHmgg = new TH2F("hmgg","2-cluster invariant mass",fNbins,fMinBin,fMaxBin,100,0,10);
202 fHmgg->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
203 fHmgg->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
204 fOutputContainer->Add(fHmgg);
206 fHmggDifferentSM = new TH2F("hmggDifferentSM","2-cluster invariant mass, different SM",fNbins,fMinBin,fMaxBin,100,0,10);
207 fHmggDifferentSM->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
208 fHmggDifferentSM->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
209 fOutputContainer->Add(fHmggDifferentSM);
211 fHOpeningAngle = new TH2F("hopang","2-cluster opening angle",100,0.,50.,100,0,10);
212 fHOpeningAngle->SetXTitle("#alpha_{#gamma #gamma}");
213 fHOpeningAngle->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
214 fOutputContainer->Add(fHOpeningAngle);
216 fHOpeningAngleDifferentSM = new TH2F("hopangDifferentSM","2-cluster opening angle, different SM",100,0,50.,100,0,10);
217 fHOpeningAngleDifferentSM->SetXTitle("#alpha_{#gamma #gamma}");
218 fHOpeningAngleDifferentSM->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
219 fOutputContainer->Add(fHOpeningAngleDifferentSM);
221 fHIncidentAngle = new TH2F("hinang","#gamma incident angle in SM",100,0.,20.,100,0,10);
222 fHIncidentAngle->SetXTitle("#alpha_{#gamma SM center}");
223 fHIncidentAngle->SetYTitle("p_{T #gamma} (GeV/c)");
224 fOutputContainer->Add(fHIncidentAngle);
226 fHIncidentAngleDifferentSM = new TH2F("hinangDifferentSM","#gamma incident angle in SM, different SM pair",100,0,20.,100,0,10);
227 fHIncidentAngleDifferentSM->SetXTitle("#alpha_{#gamma - SM center}");
228 fHIncidentAngleDifferentSM->SetYTitle("p_{T #gamma} (GeV/c)");
229 fOutputContainer->Add(fHIncidentAngleDifferentSM);
231 fHAsymmetry = new TH2F("hasym","2-cluster opening angle",100,0.,1.,100,0,10);
232 fHAsymmetry->SetXTitle("a");
233 fHAsymmetry->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
234 fOutputContainer->Add(fHAsymmetry);
236 fHAsymmetryDifferentSM = new TH2F("hasymDifferentSM","2-cluster opening angle, different SM",100,0,1.,100,0,10);
237 fHAsymmetryDifferentSM->SetXTitle("a");
238 fHAsymmetryDifferentSM->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
239 fOutputContainer->Add(fHAsymmetryDifferentSM);
242 TString pairname[] = {"A side (0-2)", "C side (1-3)","Row 0 (0-1)", "Row 1 (2-3)"};
244 for(Int_t iSM=0; iSM<4; iSM++) {
246 snprintf(hname, buffersize, "hmgg_SM%d",iSM);
247 snprintf(htitl, buffersize, "Two-gamma inv. mass for super mod %d",iSM);
248 fHmggSM[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
249 fHmggSM[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
250 fHmggSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
251 fOutputContainer->Add(fHmggSM[iSM]);
253 snprintf(hname,buffersize, "hmgg_PairSM%d",iSM);
254 snprintf(htitl,buffersize, "Two-gamma inv. mass for SM pair: %s",pairname[iSM].Data());
255 fHmggPairSM[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
256 fHmggPairSM[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
257 fHmggPairSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
258 fOutputContainer->Add(fHmggPairSM[iSM]);
261 snprintf(hname, buffersize, "hopang_SM%d",iSM);
262 snprintf(htitl, buffersize, "Opening angle for super mod %d",iSM);
263 fHOpeningAngleSM[iSM] = new TH2F(hname,htitl,100,0.,50.,100,0,10);
264 fHOpeningAngleSM[iSM]->SetXTitle("#alpha_{#gamma #gamma} (deg)");
265 fHOpeningAngleSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
266 fOutputContainer->Add(fHOpeningAngleSM[iSM]);
268 snprintf(hname,buffersize, "hopang_PairSM%d",iSM);
269 snprintf(htitl,buffersize, "Opening angle for SM pair: %s",pairname[iSM].Data());
270 fHOpeningAnglePairSM[iSM] = new TH2F(hname,htitl,100,0.,50.,100,0,10);
271 fHOpeningAnglePairSM[iSM]->SetXTitle("#alpha_{#gamma #gamma} (deg)");
272 fHOpeningAnglePairSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
273 fOutputContainer->Add(fHOpeningAnglePairSM[iSM]);
275 snprintf(hname, buffersize, "hinang_SM%d",iSM);
276 snprintf(htitl, buffersize, "Incident angle for super mod %d",iSM);
277 fHIncidentAngleSM[iSM] = new TH2F(hname,htitl,100,0.,20.,100,0,10);
278 fHIncidentAngleSM[iSM]->SetXTitle("#alpha_{#gamma - SM center} (deg)");
279 fHIncidentAngleSM[iSM]->SetYTitle("p_{T #gamma} (GeV/c)");
280 fOutputContainer->Add(fHIncidentAngleSM[iSM]);
282 snprintf(hname,buffersize, "hinang_PairSM%d",iSM);
283 snprintf(htitl,buffersize, "Incident angle for SM pair: %s",pairname[iSM].Data());
284 fHIncidentAnglePairSM[iSM] = new TH2F(hname,htitl,100,0.,20.,100,0,10);
285 fHIncidentAnglePairSM[iSM]->SetXTitle("#alpha_{#gamma - SM center} (deg)");
286 fHIncidentAnglePairSM[iSM]->SetYTitle("p_{T #gamma} (GeV/c)");
287 fOutputContainer->Add(fHIncidentAnglePairSM[iSM]);
289 snprintf(hname, buffersize, "hasym_SM%d",iSM);
290 snprintf(htitl, buffersize, "asymmetry for super mod %d",iSM);
291 fHAsymmetrySM[iSM] = new TH2F(hname,htitl,100,0.,1.,100,0,10);
292 fHAsymmetrySM[iSM]->SetXTitle("a");
293 fHAsymmetrySM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
294 fOutputContainer->Add(fHAsymmetrySM[iSM]);
296 snprintf(hname,buffersize, "hasym_PairSM%d",iSM);
297 snprintf(htitl,buffersize, "Asymmetry for SM pair: %s",pairname[iSM].Data());
298 fHAsymmetryPairSM[iSM] = new TH2F(hname,htitl,100,0.,1.,100,0,10);
299 fHAsymmetryPairSM[iSM]->SetXTitle("a");
300 fHAsymmetryPairSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
301 fOutputContainer->Add(fHAsymmetryPairSM[iSM]);
307 fhTowerDecayPhotonHit[iSM] = new TH2F (Form("hTowerDecPhotonHit_Mod%d",iSM),Form("Entries in grid of cells in Module %d",iSM),
308 colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5);
309 fhTowerDecayPhotonHit[iSM]->SetYTitle("row (phi direction)");
310 fhTowerDecayPhotonHit[iSM]->SetXTitle("column (eta direction)");
311 fOutputContainer->Add(fhTowerDecayPhotonHit[iSM]);
313 fhTowerDecayPhotonEnergy[iSM] = new TH2F (Form("hTowerDecPhotonEnergy_Mod%d",iSM),Form("Accumulated energy in grid of cells in Module %d",iSM),
314 colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5);
315 fhTowerDecayPhotonEnergy[iSM]->SetYTitle("row (phi direction)");
316 fhTowerDecayPhotonEnergy[iSM]->SetXTitle("column (eta direction)");
317 fOutputContainer->Add(fhTowerDecayPhotonEnergy[iSM]);
319 fhTowerDecayPhotonAsymmetry[iSM] = new TH2F (Form("hTowerDecPhotonAsymmetry_Mod%d",iSM),Form("Accumulated asymmetry in grid of cells in Module %d",iSM),
320 colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5);
321 fhTowerDecayPhotonAsymmetry[iSM]->SetYTitle("row (phi direction)");
322 fhTowerDecayPhotonAsymmetry[iSM]->SetXTitle("column (eta direction)");
323 fOutputContainer->Add(fhTowerDecayPhotonAsymmetry[iSM]);
327 fhNEvents = new TH1I("hNEvents", "Number of analyzed events" , 1 , 0 , 1 ) ;
328 fOutputContainer->Add(fhNEvents);
330 fOutputContainer->SetOwner(kTRUE);
332 // fCalibData = new AliEMCALCalibData();
334 PostData(1,fOutputContainer);
338 //__________________________________________________
339 void AliAnalysisTaskEMCALPi0CalibSelection::UserExec(Option_t* /* option */)
341 //Analysis per event.
343 if(fRecoUtils->GetParticleType()!=AliEMCALRecoUtils::kPhoton){
344 printf("Wrong particle type for cluster position recalculation! = %d\n", fRecoUtils->GetParticleType());
348 fhNEvents->Fill(0); //Event analyzed
350 //Get the input event
351 AliVEvent* event = 0;
352 if(fFilteredInput) event = AODEvent();
353 else event = InputEvent();
356 printf("Input event not available!\n");
361 printf("AliAnalysisTaskEMCALPi0CalibSelection <<< %s: Event %d >>>\n",event->GetName(), (Int_t)Entry());
364 //Get the primary vertex
366 event->GetPrimaryVertex()->GetXYZ(v) ;
368 if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Vertex: (%.3f,%.3f,%.3f)\n",v[0],v[1],v[2]);
370 //Int_t runNum = aod->GetRunNumber();
371 //if(DebugLevel() > 1) printf("Run number: %d\n",runNum);
373 //Get the matrix with geometry information
374 if(fhNEvents->GetEntries()==1){
376 printf("AliAnalysisTaskEMCALPi0CalibSelection::UserExec() - Load user defined geometry matrices\n");
377 for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
380 fMatrix[mod]->Print();
381 fEMCALGeo->SetMisalMatrix(fMatrix[mod],mod) ;
385 else if(!gGeoManager){
386 printf("AliAnalysisTaskEMCALPi0CalibSelection::UserExec() - Get geo matrices from data\n");
387 //Still not implemented in AOD, just a workaround to be able to work at least with ESDs
388 if(!strcmp(event->GetName(),"AliAODEvent")) {
390 printf("AliAnalysisTaskEMCALPi0CalibSelection Use ideal geometry, values geometry matrix not kept in AODs.\n");
393 if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Load Misaligned matrices. \n");
394 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event) ;
396 printf("AliAnalysisTaskEMCALPi0CalibSelection::UserExec() - This event does not contain ESDs?");
399 for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
400 //if(DebugLevel() > 1)
401 esd->GetEMCALMatrix(mod)->Print();
402 if(esd->GetEMCALMatrix(mod)) fEMCALGeo->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
405 }//Load matrices from Data
408 if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Will use fLogWeight %.3f .\n",fLogWeight);
417 Bool_t shared = kFALSE;
423 //Get the list of clusters
424 TRefArray * caloClustersArr = new TRefArray();
425 if(!fOldAOD) event->GetEMCALClusters(caloClustersArr);
426 else GetEMCALClusters(event,caloClustersArr);
427 const Int_t kNumberOfEMCALClusters = caloClustersArr->GetEntries() ;
428 if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection - N CaloClusters: %d \n", kNumberOfEMCALClusters);
431 AliVCaloCells *emCells = event->GetEMCALCells();
433 // loop over EMCAL clusters
434 //----------------------------------------------------------
435 // First recalibrate and recalculate energy and position
436 Float_t pos[]={0,0,0};
437 if(fCorrectClusters){
438 for(Int_t iClu=0; iClu<kNumberOfEMCALClusters-1; iClu++) {
439 AliVCluster *c1 = (AliVCluster *) caloClustersArr->At(iClu);
441 if(fRecoUtils->ClusterContainsBadChannel(fEMCALGeo, c1->GetCellsAbsId(), c1->GetNCells())) continue;
445 printf("Std : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",c1->GetID(),c1->E(),c1->GetDispersion(),c1->GetM02(),c1->GetM20());
446 c1->GetPosition(pos);
447 printf("Std : i %d, x %f, y %f, z %f\n",c1->GetID(), pos[0], pos[1], pos[2]);
450 //Correct cluster energy and position if requested, and not corrected previously, by default Off
451 if(fRecoUtils->IsRecalibrationOn()) {
452 fRecoUtils->RecalibrateClusterEnergy(fEMCALGeo, c1, emCells);
453 fRecoUtils->RecalculateClusterShowerShapeParameters(fEMCALGeo, emCells,c1);
454 fRecoUtils->RecalculateClusterPID(c1);
457 printf("Energy: after recalibration %f; ",c1->E());
459 // Correct Non-Linearity
460 c1->SetE(fRecoUtils->CorrectClusterEnergyLinearity(c1));
463 printf("after linearity correction %f\n",c1->E());
464 // Recalculate cluster position
465 fRecoUtils->RecalculateClusterPosition(fEMCALGeo, emCells,c1);
468 printf("Cor : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",c1->GetID(),c1->E(),c1->GetDispersion(),c1->GetM02(),c1->GetM20());
469 c1->GetPosition(pos);
470 printf("Cor : i %d, x %f, y %f, z %f\n",c1->GetID(), pos[0], pos[1], pos[2]);
475 //----------------------------------------------------------
476 //Now the invariant mass analysis with the corrected clusters
477 for(Int_t iClu=0; iClu<kNumberOfEMCALClusters-1; iClu++) {
479 AliVCluster *c1 = (AliVCluster *) caloClustersArr->At(iClu);
480 if(fRecoUtils->ClusterContainsBadChannel(fEMCALGeo, c1->GetCellsAbsId(), c1->GetNCells())) continue;
482 Float_t e1i = c1->E(); // cluster energy before correction
483 if(e1i < fEmin) continue;
484 else if(e1i > fEmax) continue;
485 else if (c1->GetNCells() < fMinNCells) continue;
489 printf("IMA : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",c1->GetID(),e1i,c1->GetDispersion(),c1->GetM02(),c1->GetM20());
490 c1->GetPosition(pos);
491 printf("IMA : i %d, x %f, y %f, z %f\n",c1->GetID(), pos[0], pos[1], pos[2]);
494 //AliEMCALAodCluster newc1(*((AliAODCaloCluster*)c1));
495 //newc1.EvalAllFromRecoUtils(fEMCALGeo,fRecoUtils,emCells);
496 //printf("i %d, recal? %d\n",iClu,newc1.IsRecalibrated());
497 //clu1.Recalibrate(fCalibData, emCells, fEMCALGeoName);
499 //clu1.EvalAll(fLogWeight, fEMCALGeoName);
501 fRecoUtils->GetMaxEnergyCell(fEMCALGeo, emCells,c1,absId1,iSupMod1,ieta1,iphi1,shared);
502 c1->GetMomentum(p1,v);
503 //newc1.GetMomentum(p1,v);
505 // Combine cluster with other clusters and get the invariant mass
506 for (Int_t jClu=iClu+1; jClu<kNumberOfEMCALClusters; jClu++) {
507 AliAODCaloCluster *c2 = (AliAODCaloCluster *) caloClustersArr->At(jClu);
508 //if(c2->IsEqual(c1)) continue;
509 if(fRecoUtils->ClusterContainsBadChannel(fEMCALGeo, c2->GetCellsAbsId(), c2->GetNCells())) continue;
511 Float_t e2i = c2->E();
512 if(e2i < fEmin) continue;
513 else if (e2i > fEmax) continue;
514 else if (c2->GetNCells() < fMinNCells) continue;
516 //AliEMCALAodCluster newc2(*((AliAODCaloCluster*)c2));
517 //newc2.EvalAllFromRecoUtils(fEMCALGeo,fRecoUtils,emCells);
518 //printf("\t j %d, recal? %d\n",jClu,newc2.IsRecalibrated());
519 //clu2.Recalibrate(fCalibData, emCells,fEMCALGeoName);
521 //clu2.EvalAll(fLogWeight,fEMCALGeoName);
523 fRecoUtils->GetMaxEnergyCell(fEMCALGeo, emCells,c2,absId2,iSupMod2,ieta2,iphi2,shared);
524 c2->GetMomentum(p2,v);
525 //newc2.GetMomentum(p2,v);
527 Float_t invmass = p12.M()*1000;
528 //printf("*** mass %f\n",invmass);
529 Float_t asym = TMath::Abs(p1.E()-p2.E())/(p1.E()+p2.E());
530 //printf("asymmetry %f\n",asym);
532 if(asym > fAsyCut) continue;
534 if(invmass < fMaxBin && invmass > fMinBin){
536 //Check if cluster is in fidutial region, not too close to borders
537 Bool_t in1 = fRecoUtils->CheckCellFiducialRegion(fEMCALGeo, c1, emCells);
538 Bool_t in2 = fRecoUtils->CheckCellFiducialRegion(fEMCALGeo, c2, emCells);
542 fHmgg->Fill(invmass,p12.Pt());
544 if(iSupMod1==iSupMod2) fHmggSM[iSupMod1]->Fill(invmass,p12.Pt());
545 else fHmggDifferentSM ->Fill(invmass,p12.Pt());
547 if((iSupMod1==0 && iSupMod2==2) || (iSupMod1==2 && iSupMod2==0)) fHmggPairSM[0]->Fill(invmass,p12.Pt());
548 if((iSupMod1==1 && iSupMod2==3) || (iSupMod1==3 && iSupMod2==1)) fHmggPairSM[1]->Fill(invmass,p12.Pt());
549 if((iSupMod1==0 && iSupMod2==1) || (iSupMod1==1 && iSupMod2==0)) fHmggPairSM[2]->Fill(invmass,p12.Pt());
550 if((iSupMod1==2 && iSupMod2==3) || (iSupMod1==3 && iSupMod2==2)) fHmggPairSM[3]->Fill(invmass,p12.Pt());
552 if(invmass > 100. && invmass < 160.){//restrict to clusters really close to pi0 peak
554 //Opening angle of 2 photons
555 Float_t opangle = p1.Angle(p2.Vect())*TMath::RadToDeg();
556 //printf("*******>>>>>>>> In PEAK pt %f, angle %f \n",p12.Pt(),opangle);
558 //Incident angle of each photon
559 Float_t inangle1 =0., inangle2=0.;
561 Float_t posSM1cen[3]={0.,0.,0.};
562 Float_t depth = fRecoUtils->GetDepth(p1.Energy(),fRecoUtils->GetParticleType(),iSupMod1);
563 fEMCALGeo->RecalculateTowerPosition(11.5, 23.5, iSupMod1, depth, fRecoUtils->GetMisalTransShiftArray(),fRecoUtils->GetMisalRotShiftArray(),posSM1cen);
564 Float_t posSM2cen[3]={0.,0.,0.};
565 depth = fRecoUtils->GetDepth(p2.Energy(),fRecoUtils->GetParticleType(),iSupMod2);
566 fEMCALGeo->RecalculateTowerPosition(11.5, 23.5, iSupMod2, depth, fRecoUtils->GetMisalTransShiftArray(),fRecoUtils->GetMisalRotShiftArray(),posSM2cen);
567 //printf("SM1 %d pos (%2.3f,%2.3f,%2.3f) \n",iSupMod1,posSM1cen[0],posSM1cen[1],posSM1cen[2]);
568 //printf("SM2 %d pos (%2.3f,%2.3f,%2.3f) \n",iSupMod2,posSM2cen[0],posSM2cen[1],posSM2cen[2]);
570 TVector3 vecSM1cen(posSM1cen[0]-v[0],posSM1cen[1]-v[1],posSM1cen[2]-v[2]);
571 TVector3 vecSM2cen(posSM2cen[0]-v[0],posSM2cen[1]-v[1],posSM2cen[2]-v[2]);
572 inangle1 = p1.Angle(vecSM1cen)*TMath::RadToDeg();
573 inangle2 = p2.Angle(vecSM2cen)*TMath::RadToDeg();
574 //printf("Incident angle: cluster 1 %2.3f; cluster 2 %2.3f\n",inangle1,inangle2);
576 fHOpeningAngle ->Fill(opangle,p12.Pt());
577 fHIncidentAngle->Fill(inangle1,p1.Pt());
578 fHIncidentAngle->Fill(inangle2,p2.Pt());
579 fHAsymmetry ->Fill(asym,p12.Pt());
581 if(iSupMod1==iSupMod2) {
582 fHOpeningAngleSM[iSupMod1] ->Fill(opangle,p12.Pt());
583 fHIncidentAngleSM[iSupMod1]->Fill(inangle1,p1.Pt());
584 fHIncidentAngleSM[iSupMod1]->Fill(inangle2,p2.Pt());
585 fHAsymmetrySM[iSupMod1] ->Fill(asym,p12.Pt());
588 fHOpeningAngleDifferentSM ->Fill(opangle,p12.Pt());
589 fHIncidentAngleDifferentSM ->Fill(inangle1,p1.Pt());
590 fHIncidentAngleDifferentSM ->Fill(inangle2,p2.Pt());
591 fHAsymmetryDifferentSM ->Fill(asym,p12.Pt());
594 if((iSupMod1==0 && iSupMod2==2) || (iSupMod1==2 && iSupMod2==0)) {
595 fHOpeningAnglePairSM[0] ->Fill(opangle,p12.Pt());
596 fHIncidentAnglePairSM[0]->Fill(inangle1,p1.Pt());
597 fHIncidentAnglePairSM[0]->Fill(inangle2,p2.Pt());
598 fHAsymmetryPairSM[0] ->Fill(asym,p12.Pt());
601 if((iSupMod1==1 && iSupMod2==3) || (iSupMod1==3 && iSupMod2==1)) {
602 fHOpeningAnglePairSM[1] ->Fill(opangle,p12.Pt());
603 fHIncidentAnglePairSM[1]->Fill(inangle1,p1.Pt());
604 fHIncidentAnglePairSM[1]->Fill(inangle2,p2.Pt());
605 fHAsymmetryPairSM[1] ->Fill(asym,p12.Pt());
609 if((iSupMod1==0 && iSupMod2==1) || (iSupMod1==1 && iSupMod2==0)) {
610 fHOpeningAnglePairSM[2] ->Fill(opangle,p12.Pt());
611 fHIncidentAnglePairSM[2]->Fill(inangle1,p1.Pt());
612 fHIncidentAnglePairSM[2]->Fill(inangle2,p2.Pt());
613 fHAsymmetryPairSM[2] ->Fill(asym,p12.Pt());
617 if((iSupMod1==2 && iSupMod2==3) || (iSupMod1==3 && iSupMod2==2)) {
618 fHOpeningAnglePairSM[3] ->Fill(opangle,p12.Pt());
619 fHIncidentAnglePairSM[3]->Fill(inangle1,p1.Pt());
620 fHIncidentAnglePairSM[3]->Fill(inangle2,p2.Pt());
621 fHAsymmetryPairSM[3] ->Fill(asym,p12.Pt());
624 }// pair in 100 < mass < 160
626 }//in acceptance cuts
628 //In case of filling only channels with second cluster in same SM
629 if(fSameSM && iSupMod1!=iSupMod2) continue;
631 if (fGroupNCells == 0){
632 fHmpi0[iSupMod1][ieta1][iphi1]->Fill(invmass);
633 fHmpi0[iSupMod2][ieta2][iphi2]->Fill(invmass);
635 if(invmass > 100. && invmass < 160.){//restrict to clusters really close to pi0 peak
636 fhTowerDecayPhotonHit [iSupMod1]->Fill(ieta1,iphi1);
637 fhTowerDecayPhotonEnergy [iSupMod1]->Fill(ieta1,iphi1,p1.E());
638 fhTowerDecayPhotonAsymmetry[iSupMod1]->Fill(ieta1,iphi1,asym);
640 fhTowerDecayPhotonHit [iSupMod2]->Fill(ieta2,iphi2);
641 fhTowerDecayPhotonEnergy [iSupMod2]->Fill(ieta2,iphi2,p2.E());
642 fhTowerDecayPhotonAsymmetry[iSupMod2]->Fill(ieta2,iphi2,asym);
644 }// pair in mass of pi0
647 //printf("Regroup N %d, eta1 %d, phi1 %d, eta2 %d, phi2 %d \n",fGroupNCells, ieta1, iphi1, ieta2, iphi2);
648 for (Int_t i = -fGroupNCells; i < fGroupNCells+1; i++) {
649 for (Int_t j = -fGroupNCells; j < fGroupNCells+1; j++) {
650 //printf("\t i %d, j %d\n",i,j);
651 if((ieta1+i >= 0) && (iphi1+j >= 0) && (ieta1+i < 48) && (iphi1+j < 24)){
652 //printf("\t \t eta1+i %d, phi1+j %d\n", ieta1+i, iphi1+j);
653 fHmpi0[iSupMod1][ieta1+i][iphi1+j]->Fill(invmass);
655 if((ieta2+i >= 0) && (iphi2+j >= 0) && (ieta2+i < 48) && (iphi2+j < 24)){
656 //printf("\t \t eta2+i %d, phi2+j %d\n", ieta2+i, iphi2+j);
657 fHmpi0[iSupMod2][ieta2+i][iphi2+j]->Fill(invmass);
663 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",
664 iSupMod1,iphi1,ieta1,iSupMod2,iphi2,ieta2,p12.M(),e1i,c1->E(),e2i,c2->E());
669 } // end of loop over EMCAL clusters
671 delete caloClustersArr;
673 PostData(1,fOutputContainer);
677 //_____________________________________________________
678 void AliAnalysisTaskEMCALPi0CalibSelection::PrintInfo(){
681 printf("Cluster cuts: %2.2f < E < %2.2f GeV; min number of cells %d; Assymetry cut %1.2f\n", fEmin,fEmax, fMinNCells, fAsyCut) ;
682 printf("Group %d cells\n", fGroupNCells) ;
683 printf("Cluster maximal cell away from border at least %d cells\n", fRecoUtils->GetNumberOfCellsFromEMCALBorder()) ;
684 printf("Histograms: bins %d; energy range: %2.2f < E < %2.2f GeV\n",fNbins,fMinBin,fMaxBin) ;
685 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",
686 fRecoUtils->IsBadChannelsRemovalSwitchedOn(),fFilteredInput,fCorrectClusters, fOldAOD, fSameSM) ;
687 printf("EMCAL Geometry name: < %s >, Load Matrices %d\n",fEMCALGeoName.Data(), fLoadMatrices) ;
688 if(fLoadMatrices) {for(Int_t ism = 0; ism < 4; ism++) fMatrix[ism]->Print();}
693 //__________________________________________________
694 //void AliAnalysisTaskEMCALPi0CalibSelection::SetCalibCorrections(AliEMCALCalibData* const cdata)
696 // //Set new correction factors (~1) to calibration coefficients, delete previous.
698 // if(fCalibData) delete fCalibData;
699 // fCalibData = cdata;