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