1)AliCaloPID: Posibility to recalculate PID bayesian in EMCAL
[u/mrichter/AliRoot.git] / PWG4 / CaloCalib / AliAnalysisTaskEMCALPi0CalibSelection.cxx
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"
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"
45
46 ClassImp(AliAnalysisTaskEMCALPi0CalibSelection)
47
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)
52 {
53   //Default constructor.
54
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;
59       }
60     } 
61   }
62                 
63 }
64
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)
69 {
70   //Named constructor which should be used.
71   
72   DefineOutput(1,TList::Class());
73   
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;
78       }
79     } 
80   }
81
82 }
83
84 //__________________________________________________
85 AliAnalysisTaskEMCALPi0CalibSelection::~AliAnalysisTaskEMCALPi0CalibSelection()
86 {
87   //Destructor.
88   
89   if(fOutputContainer){
90     fOutputContainer->Delete() ; 
91     delete fOutputContainer ;
92   }
93         
94   if(fCalibData)  delete fCalibData;
95   if(fEMCALGeo)   delete fEMCALGeo;
96         
97 }
98
99
100 //__________________________________________________
101 void AliAnalysisTaskEMCALPi0CalibSelection::CreateAODFromAOD()
102 {
103   // Copy AOD header, vertex, CaloClusters and CaloCells to output AOD
104   
105   AliAODEvent* aod = dynamic_cast<AliAODEvent*>(InputEvent());
106   
107   // set arrays and pointers
108   Float_t posF[3];
109   Double_t pos[3];
110   
111   Double_t covVtx[6];
112   
113   for (Int_t i = 0; i < 6; i++)  covVtx[i] = 0.;
114   
115   // Update the header
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()); 
125   
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);
137   //
138   //
139   Int_t nVertices = 1 ;/* = prim. vtx*/;
140   Int_t nCaloClus = aod->GetNCaloClusters();
141   
142   AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
143   
144   // Access to the AOD container of vertices
145   TClonesArray &vertices = *(AODEvent()->GetVertices());
146   Int_t jVertices=0;
147   
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();
151   
152   vtx->GetXYZ(pos); // position
153   vtx->GetCovMatrix(covVtx); //covariance matrix
154   
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());
159   
160   // Access to the AOD container of clusters
161   TClonesArray &caloClusters = *(AODEvent()->GetCaloClusters());
162   Int_t jClusters=0;
163   
164   for (Int_t iClust=0; iClust<nCaloClus; ++iClust) {
165     
166     AliAODCaloCluster * cluster = aod->GetCaloCluster(iClust);
167     
168     //Check if it is a EMCAL cluster
169     if(!cluster->IsEMCALCluster())  continue ;
170     
171     Int_t id       = cluster->GetID();
172     Float_t energy = cluster->E();
173     cluster->GetPosition(posF);
174     Char_t ttype   = cluster->GetType(); 
175     
176     AliAODCaloCluster *caloCluster = new(caloClusters[jClusters++]) 
177       AliAODCaloCluster(id,
178                         0,
179                         0x0,
180                         energy,
181                         posF,
182                         NULL,
183                         ttype);
184     
185     caloCluster->SetCaloCluster(cluster->GetDistToBadChannel(),
186                                 cluster->GetDispersion(),
187                                 cluster->GetM20(), cluster->GetM02(),
188                                 cluster->GetEmcCpvDistance(),  
189                                 cluster->GetNExMax(),cluster->GetTOF()) ;
190     
191     caloCluster->SetPIDFromESD(cluster->PID());
192     caloCluster->SetNCells(cluster->GetNCells());
193     caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
194         if(!fOldData){
195                 caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
196         }
197         else{
198                 Double_t *fraction = cluster->GetCellsAmplitudeFraction();
199                 for(Int_t i = 0; i < cluster->GetNCells() ; i++) fraction[i] = 1;
200                 caloCluster->SetCellsAmplitudeFraction(fraction);
201         }
202     
203   } 
204
205   caloClusters.Expand(jClusters); // resize TObjArray to 'remove' slots for pseudo clusters      
206   // end of loop on calo clusters
207   
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() ;
212     
213     AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
214     aodEMcells.CreateContainer(nEMcell);
215     aodEMcells.SetType(AliAODCaloCells::kEMCAL);
216         
217         Double_t calibFactor = 1;
218         if(fOldData) calibFactor = 0.0153;
219
220     for (Int_t iCell = 0; iCell < nEMcell; iCell++) {      
221       aodEMcells.SetCell(iCell,aodinEMcells.GetCellNumber(iCell),aodinEMcells.GetAmplitude(iCell)*calibFactor);
222     }
223     aodEMcells.Sort();
224           
225   }
226   
227 }
228
229 //__________________________________________________
230 void AliAnalysisTaskEMCALPi0CalibSelection::CreateAODFromESD()
231 {
232   
233   // Copy Header, Vertex, CaloClusters and CaloCells from ESDs to AODs
234   
235   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
236   
237   // set arrays and pointers
238   Float_t posF[3];
239   Double_t pos[3];
240   
241   Double_t covVtx[6];
242   
243   for (Int_t i = 0; i < 6; i++)  covVtx[i] = 0.;
244   
245   // Update the header
246   
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
255   
256   
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);
268   //
269   //
270   Int_t nVertices = 1 ;/* = prim. vtx*/;
271   Int_t nCaloClus = esd->GetNumberOfCaloClusters();
272   
273   AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
274   
275   // Access to the AOD container of vertices
276   TClonesArray &vertices = *(AODEvent()->GetVertices());
277   Int_t jVertices=0;
278   
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();
282   
283   vtx->GetXYZ(pos); // position
284   vtx->GetCovMatrix(covVtx); //covariance matrix
285   
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());
290   
291   // Access to the AOD container of clusters
292   TClonesArray &caloClusters = *(AODEvent()->GetCaloClusters());
293   Int_t jClusters=0;
294   
295   for (Int_t iClust=0; iClust<nCaloClus; ++iClust) {
296     
297     AliESDCaloCluster * cluster = esd->GetCaloCluster(iClust);
298     
299     //Check which calorimeter information we want to keep.
300     if(!cluster->IsEMCAL())  continue ;
301     
302     Int_t id       = cluster->GetID();
303     Float_t energy = cluster->E();
304     cluster->GetPosition(posF);
305     
306     AliAODCaloCluster *caloCluster = new(caloClusters[jClusters++]) 
307       AliAODCaloCluster(id,
308                         0,
309                         0x0,
310                         energy,
311                         posF,
312                         NULL,
313                         AliAODCluster::kEMCALClusterv1);
314     
315     caloCluster->SetCaloCluster(cluster->GetDistanceToBadChannel(),
316                                 cluster->GetClusterDisp(),
317                                 cluster->GetM20(), cluster->GetM02(),
318                                 cluster->GetEmcCpvDistance(),  
319                                 cluster->GetNExMax(),cluster->GetTOF()) ;
320     
321     caloCluster->SetPIDFromESD(cluster->GetPid());
322     caloCluster->SetNCells(cluster->GetNCells());
323     caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
324         if(!fOldData){
325                   caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
326         }
327         else{
328                   Double_t *fraction = cluster->GetCellsAmplitudeFraction();
329                   for(Int_t i = 0; i < cluster->GetNCells() ; i++) fraction[i] = 1;
330                   caloCluster->SetCellsAmplitudeFraction(fraction);
331           }    
332   } 
333   caloClusters.Expand(jClusters); // resize TObjArray to 'remove' slots for pseudo clusters      
334   // end of loop on calo clusters
335   
336   // fill EMCAL cell info
337
338   if( esd->GetEMCALCells()) { // protection against missing ESD information
339     AliESDCaloCells &esdEMcells = *(esd->GetEMCALCells());
340         Int_t nEMcell = esdEMcells.GetNumberOfCells() ;
341     
342     AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
343     aodEMcells.CreateContainer(nEMcell);
344     aodEMcells.SetType(AliAODCaloCells::kEMCAL);  
345           
346         Double_t calibFactor = 1;
347         if(fOldData) calibFactor = 0.0153;
348    
349         for (Int_t iCell = 0; iCell < nEMcell; iCell++) {      
350       aodEMcells.SetCell(iCell,esdEMcells.GetCellNumber(iCell),esdEMcells.GetAmplitude(iCell)*calibFactor);
351     }
352     aodEMcells.Sort();
353           
354   }
355
356 }
357
358 //__________________________________________________
359 void AliAnalysisTaskEMCALPi0CalibSelection::UserCreateOutputObjects()
360 {
361   //Create output container
362   fOutputContainer = new TList();
363   
364   char hname[128], htitl[128];
365   
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]);
373       }
374     }
375   }
376
377   fHmgg = new TH1F("hmgg","2-cluster invariant mass",100,0.,300.);
378   fOutputContainer->Add(fHmgg);
379         
380   fCalibData = new AliEMCALCalibData();
381         
382   if(!fCalibData)
383     AliFatal("Calibration parameters not found in CDB!");
384   
385   printf("Get Geometry : %s\n", fEMCALGeoName.Data());
386   fEMCALGeo =  AliEMCALGeometry::GetInstance(fEMCALGeoName) ;   
387
388 }
389
390 //__________________________________________________
391 void AliAnalysisTaskEMCALPi0CalibSelection::UserExec(Option_t* /* option */)
392 {
393   //Analysis per event.
394   if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection <<< Event %d >>>\n",(Int_t)Entry());
395   
396   AliAODEvent* aod = 0x0;
397   if(!strcmp(InputEvent()->GetName(),"AliAODEvent")) {
398     //Input are ESDs
399     aod = dynamic_cast<AliAODEvent*>(InputEvent());
400     // Create new AOD with only CaloClusters and CaloCells
401     if(fCopyAOD) CreateAODFromAOD();
402   }
403   else  if(!strcmp(InputEvent()->GetName(),"AliESDEvent")) {
404     //Input are ESDs
405     aod = AODEvent();
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();
409   }
410   else {
411     printf("AliAnalysisTaskEMCALPi0CalibSelection: Unknown event type, STOP!\n");
412     abort();
413   }     
414   
415   Double_t v[] = {aod->GetVertex(0)->GetX(),aod->GetVertex(0)->GetY(),aod->GetVertex(0)->GetZ()}; //to check!!
416   //aod->GetVertex()->GetXYZ(v) ;
417   //TVector3 vtx(v); 
418   
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]);
421   
422   Int_t runNum = aod->GetRunNumber();
423   if(DebugLevel() > 1) printf("Run number: %d\n",runNum);
424   
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")) {
428     if(DebugLevel() > 1) 
429       printf("AliAnalysisTaskEMCALPi0CalibSelection Use ideal geometry, values geometry matrix not kept in AODs.\n");
430   }
431   else{ 
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) ;
436     }
437   }
438   
439   if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Will use fLogWeight %.3f .\n",fLogWeight);
440   
441   //AliEMCALPID pid;
442   
443   Int_t maxId;
444   Int_t iSupMod = -1;
445   Int_t iTower  = -1;
446   Int_t iIphi   = -1;
447   Int_t iIeta   = -1;
448   Int_t iphi   = -1;
449   Int_t ieta   = -1;
450   
451   TLorentzVector p1;
452   TLorentzVector p2;
453   TLorentzVector p12;
454   
455   
456   TRefArray * caloClustersArr  = new TRefArray();
457   aod->GetEMCALClusters(caloClustersArr);
458   
459   const Int_t kNumberOfEMCALClusters   = caloClustersArr->GetEntries() ;
460   if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection CaloClusters: %d\n", kNumberOfEMCALClusters);
461   
462   // EMCAL cells
463   AliAODCaloCells *emCells = aod->GetEMCALCells();
464         
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);
472           }
473           
474           //Change fraction
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);        
481           }
482   }
483         
484   // loop over EMCAL clusters
485   for(Int_t iClu=0; iClu<kNumberOfEMCALClusters; iClu++) {
486     
487     AliAODCaloCluster *c1 = (AliAODCaloCluster *) caloClustersArr->At(iClu);
488     if(!c1->IsEMCALCluster()) continue; // EMCAL cluster!
489                   
490         Float_t e1i = c1->E();   // cluster energy before correction   
491         if(e1i<fEmin) continue;
492
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]);
498         }
499           
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) ;
505           
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]);
511         }
512           
513         clu1.GetMomentum(p1,v);
514
515         MaxEnergyCellPos(emCells,&clu1,maxId);
516     
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);    
522     
523     Float_t e1ii = clu1.E(); // cluster energy after correction
524     
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;
529       
530       Float_t e2i = c2->E();
531       if(e2i<fEmin) continue;
532       
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());
539
540       clu2.GetMomentum(p2,v);
541       Float_t e2ii = clu2.E();
542       
543       p12 = p1+p2;
544       Float_t invmass = p12.M()*1000; 
545
546       fHmgg->Fill(invmass); 
547
548       //printf("iSM %d, ieta %d, iphi %d,  mass %f  \n",iSupMod, ieta, iphi, invmass);
549       if(invmass < 300)
550           fHmpi0[iSupMod][ieta][iphi]->Fill(invmass);
551       
552       
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);
555     }
556     
557   } // end of loop over EMCAL clusters
558   
559   delete caloClustersArr;
560   PostData(1,fOutputContainer);
561 }
562
563 //__________________________________________________
564 void AliAnalysisTaskEMCALPi0CalibSelection::MaxEnergyCellPos(AliAODCaloCells* const cells, AliAODCaloCluster* const clu, Int_t& maxId)
565 {
566   //For a given CaloCluster calculates the absId of the cell 
567   //with maximum energy deposit.
568   
569   Double_t eMax = -111;
570   
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);
574     if(eCell>eMax)  { 
575       eMax = eCell; 
576       maxId = cellAbsId;
577     }
578   }
579
580 }
581
582 //__________________________________________________
583 void AliAnalysisTaskEMCALPi0CalibSelection::SetCalibCorrections(AliEMCALCalibData* const cdata)
584 {
585   //Set new correction factors (~1) to calibration coefficients, delete previous.
586
587    if(fCalibData) delete fCalibData;
588    fCalibData = cdata;
589         
590 }
591