PWGGA/EMCALTasks/AliAnalysisTaskEMCALPi0CalibSelection.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3
   4  * Permission to use, copy, modify and distribute this software and its   *
   5  * documentation strictly for non-commercial purposes is hereby granted   *
   6  * without fee, provided that the above copyright notice appears in all   *
   7  * copies and that both the copyright notice and this permission notice   *
   8  * appear in the supporting documentation. The authors make no claims     *
   9  * about the suitability of this software for any purpose. It is          *
  10  * provided "as is" without express or implied warranty.                  *
  11  **************************************************************************/
  12
  13 // $Id$
  14
  15 //---------------------------------------------------------------------------//
  16 //                                                                           //
  17 // Fill histograms (one per cell) with two-cluster invariant mass            //
  18 // using calibration coefficients of the previous iteration.                 //
  19 // Histogram for a given cell is filled if the most energy of one cluster    //
  20 // is deposited in this cell and the other cluster could be anywherein EMCAL.//
  21 //                                                                           //
  22 //                                                                           //
  23 // Author: Boris Polishchuk                                                  //
  24 // Adapted to AOD reading by Gustavo Conesa                                  //
  25 //                                                                           //
  26 //                                                                           //
  27 //---------------------------------------------------------------------------//
  28
  29 // Root
  30 #include "TLorentzVector.h"
  31 #include "TRefArray.h"
  32 #include "TList.h"
  33 #include "TH1F.h"
  34 #include <TGeoManager.h>
  35
  36 // AliRoot
  37 #include "AliAnalysisTaskEMCALPi0CalibSelection.h"
  38 #include "AliAODEvent.h"
  39 #include "AliESDEvent.h"
  40 #include "AliEMCALGeometry.h"
  41 #include "AliVCluster.h"
  42 #include "AliVCaloCells.h"
  43 #include "AliEMCALRecoUtils.h"
  44 #include "AliOADBContainer.h"
  45
  46 ClassImp(AliAnalysisTaskEMCALPi0CalibSelection)
  47
  48
  49 //______________________________________________________________________________________________
  50 AliAnalysisTaskEMCALPi0CalibSelection::AliAnalysisTaskEMCALPi0CalibSelection(const char* name) :
  51 AliAnalysisTaskSE(name),
  52 fEMCALGeo(0x0),           fLoadMatrices(0),
  53 fEMCALGeoName("EMCAL_COMPLETE12SMV1"),
  54 fTriggerName("EMC"),
  55 fRecoUtils(new AliEMCALRecoUtils),
  56 fOADBFilePath(""),        fCorrectClusters(kFALSE),
  57 fCaloClustersArr(0x0),    fEMCALCells(0x0),
  58 fCuts(0x0),               fOutputContainer(0x0),
  59 fVertex(),                fFilteredInput(kFALSE),
  60 fEmin(0.5),               fEmax(15.),
  61 fL0min(0.01),             fL0max(0.5),
  62 fDTimeCut(100.),          fTimeMax(1000000),        fTimeMin(-1000000),
  63 fAsyCut(1.),              fMinNCells(2),            fGroupNCells(0),
  64 fLogWeight(4.5),          fSameSM(kFALSE),
  65 fNMaskCellColumns(11),    fMaskCellColumns(0x0),
  66 fInvMassCutMin(110.),     fInvMassCutMax(160.),
  67 //Histograms
  68 fNbins(300),
  69 fMinBin(0.),              fMaxBin(300.),
  70 fNTimeBins(1000),         fMinTimeBin(0.),          fMaxTimeBin(1000.),
  71 fHmgg(0x0),               fHmggDifferentSM(0x0),
  72 fHmggMaskFrame(0x0),      fHmggDifferentSMMaskFrame(0x0),
  73 fHOpeningAngle(0x0),      fHOpeningAngleDifferentSM(0x0),
  74 fHAsymmetry(0x0),         fHAsymmetryDifferentSM(0x0),
  75 fhNEvents(0x0),
  76 fhClusterTime(0x0),       fhClusterPairDiffTime(0x0)
  77 {
  78   //Named constructor which should be used.
  79
  80   for(Int_t iMod=0; iMod < AliEMCALGeoParams::fgkEMCALModules; iMod++) {
  81     for(Int_t iX=0; iX<24; iX++) {
  82       for(Int_t iZ=0; iZ<48; iZ++) {
  83         fHmpi0[iMod][iZ][iX]   = 0 ;
  84       }
  85     }
  86   }
  87
  88   fVertex[0]=fVertex[1]=fVertex[2]=-1000;
  89
  90   fHTpi0[0]= 0 ;
  91   fHTpi0[1]= 0 ;
  92   fHTpi0[2]= 0 ;
  93   fHTpi0[3]= 0 ;
  94
  95   fMaskCellColumns = new Int_t[fNMaskCellColumns];
  96   fMaskCellColumns[0] = 6 ;  fMaskCellColumns[1] = 7 ;  fMaskCellColumns[2] = 8 ;
  97   fMaskCellColumns[3] = 35;  fMaskCellColumns[4] = 36;  fMaskCellColumns[5] = 37;
  98   fMaskCellColumns[6] = 12+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[7] = 13+AliEMCALGeoParams::fgkEMCALCols;
  99   fMaskCellColumns[8] = 40+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[9] = 41+AliEMCALGeoParams::fgkEMCALCols;
 100   fMaskCellColumns[10]= 42+AliEMCALGeoParams::fgkEMCALCols;
 101
 102   for(Int_t iSMPair = 0; iSMPair < AliEMCALGeoParams::fgkEMCALModules/2; iSMPair++)
 103   {
 104     fHmggPairSameSectorSM[iSMPair]          = 0;
 105     fHmggPairSameSectorSMMaskFrame[iSMPair] = 0;
 106     fhClusterPairDiffTimeSameSector[iSMPair]= 0;
 107   }
 108
 109   for(Int_t iSMPair = 0; iSMPair < AliEMCALGeoParams::fgkEMCALModules-2; iSMPair++)
 110   {
 111     fHmggPairSameSideSM[iSMPair]            = 0;
 112     fHmggPairSameSideSMMaskFrame[iSMPair]   = 0;
 113     fhClusterPairDiffTimeSameSide[iSMPair]  = 0;
 114   }
 115
 116   for(Int_t iSM = 0; iSM < AliEMCALGeoParams::fgkEMCALModules; iSM++)
 117   {
 118     fHmggSM[iSM]                     = 0;
 119     fHmggSMMaskFrame[iSM]            = 0;
 120     fHOpeningAngleSM[iSM]            = 0;
 121     fHOpeningAnglePairSM[iSM]        = 0;
 122     fHAsymmetrySM[iSM]               = 0;
 123     fHAsymmetryPairSM[iSM]           = 0;
 124     fhTowerDecayPhotonHit[iSM]       = 0;
 125     fhTowerDecayPhotonEnergy[iSM]    = 0;
 126     fhTowerDecayPhotonAsymmetry[iSM] = 0;
 127     fhTowerDecayPhotonHitMaskFrame[iSM]= 0;
 128     fMatrix[iSM]                     = 0x0;
 129     fhClusterTimeSM[iSM]             = 0;
 130     fhClusterPairDiffTimeSameSM[iSM] = 0;
 131   }
 132
 133   DefineOutput(1, TList::Class());
 134   DefineOutput(2, TList::Class());  // will contain cuts or local params
 135
 136 }
 137
 138 //_____________________________________________________________________________
 139 AliAnalysisTaskEMCALPi0CalibSelection::~AliAnalysisTaskEMCALPi0CalibSelection()
 140 {
 141   //Destructor.
 142
 143   if(fOutputContainer)
 144   {
 145     fOutputContainer->Delete() ;
 146     delete fOutputContainer ;
 147   }
 148
 149   if(fEMCALGeo)         delete fEMCALGeo  ;
 150   if(fRecoUtils)        delete fRecoUtils ;
 151   if(fNMaskCellColumns) delete [] fMaskCellColumns;
 152
 153 }
 154
 155 //____________________________________________________________
 156 void  AliAnalysisTaskEMCALPi0CalibSelection::CorrectClusters()
 157 {
 158   // loop over EMCAL clusters
 159   //----------------------------------------------------------
 160   // First recalibrate and recalculate energy and position
 161
 162
 163   if(fCorrectClusters)
 164   {
 165
 166     if(fRecoUtils->GetParticleType()!=AliEMCALRecoUtils::kPhoton)
 167     {
 168       printf("Wrong particle type for cluster position recalculation! = %d\n", fRecoUtils->GetParticleType());
 169       abort();
 170     }
 171
 172     if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Will use fLogWeight %.3f .\n",fLogWeight);
 173
 174     Float_t pos[]={0,0,0};
 175
 176     for(Int_t iClu=0; iClu < fCaloClustersArr->GetEntriesFast(); iClu++)
 177     {
 178       AliVCluster *c1 = (AliVCluster *) fCaloClustersArr->At(iClu);
 179
 180       Float_t e1i = c1->E();   // cluster energy before correction
 181       if      (e1i < fEmin) continue;
 182       else if (e1i > fEmax) continue;
 183
 184       else if (c1->GetNCells() < fMinNCells)                   continue;
 185
 186       else if (c1->GetM02() < fL0min || c1->GetM02() > fL0max) continue;
 187
 188       if(fRecoUtils->ClusterContainsBadChannel(fEMCALGeo, c1->GetCellsAbsId(), c1->GetNCells())) continue;
 189
 190       if(DebugLevel() > 2)
 191       {
 192         printf("Std  : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",c1->GetID(),c1->E(),c1->GetDispersion(),c1->GetM02(),c1->GetM20());
 193         c1->GetPosition(pos);
 194         printf("Std  : i %d, x %f, y %f, z %f\n",c1->GetID(), pos[0], pos[1], pos[2]);
 195       }
 196
 197       //Correct cluster energy and position if requested, and not corrected previously, by default Off
 198       if(fRecoUtils->IsRecalibrationOn())
 199       {
 200         fRecoUtils->RecalibrateClusterEnergy(fEMCALGeo, c1, fEMCALCells);
 201         fRecoUtils->RecalculateClusterShowerShapeParameters(fEMCALGeo, fEMCALCells,c1);
 202         fRecoUtils->RecalculateClusterPID(c1);
 203       }
 204
 205       if(DebugLevel() > 2)
 206         printf("Energy: after recalibration %f; \n",c1->E());
 207
 208       // Recalculate cluster position
 209       fRecoUtils->RecalculateClusterPosition(fEMCALGeo, fEMCALCells,c1);
 210
 211       // Correct Non-Linearity
 212       c1->SetE(fRecoUtils->CorrectClusterEnergyLinearity(c1));
 213
 214       if(DebugLevel() > 2)
 215         printf("\t after linearity correction %f\n",c1->E());
 216
 217       //In case of MC analysis, to match resolution/calibration in real data
 218       c1->SetE(fRecoUtils->SmearClusterEnergy(c1));
 219
 220       if(DebugLevel() > 2)
 221         printf("\t after smearing %f\n",c1->E());
 222
 223       if(DebugLevel() > 2)
 224       {
 225         printf("Cor  : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",c1->GetID(),c1->E(),c1->GetDispersion(),c1->GetM02(),c1->GetM20());
 226         c1->GetPosition(pos);
 227         printf("Cor  : i %d, x %f, y %f, z %f\n",c1->GetID(), pos[0], pos[1], pos[2]);
 228       }
 229     }
 230   }
 231 }
 232
 233 //__________________________________________________________
 234 void AliAnalysisTaskEMCALPi0CalibSelection::FillHistograms()
 235 {
 236   // Now fill the invariant mass analysis with the corrected clusters, and other general histograms
 237
 238   Int_t absId1   = -1;
 239   Int_t iSupMod1 = -1;
 240   Int_t iphi1    = -1;
 241   Int_t ieta1    = -1;
 242   Int_t absId2   = -1;
 243   Int_t iSupMod2 = -1;
 244   Int_t iphi2    = -1;
 245   Int_t ieta2    = -1;
 246   Bool_t shared  = kFALSE;
 247
 248   TLorentzVector p1;
 249   TLorentzVector p2;
 250   TLorentzVector p12;
 251
 252   Float_t pos[]={0,0,0};
 253
 254   Int_t bc  = InputEvent()->GetBunchCrossNumber();
 255   Int_t nSM = (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules();
 256
 257   for(Int_t iClu=0; iClu<fCaloClustersArr->GetEntriesFast()-1; iClu++)
 258   {
 259     AliVCluster *c1 = (AliVCluster *) fCaloClustersArr->At(iClu);
 260
 261     if(fRecoUtils->ClusterContainsBadChannel(fEMCALGeo, c1->GetCellsAbsId(), c1->GetNCells())) continue;
 262
 263     Float_t e1i = c1->E();   // cluster energy before correction
 264
 265     if      (e1i < fEmin) continue;
 266     else if (e1i > fEmax) continue;
 267
 268     else if (!fRecoUtils->IsGoodCluster(c1,fEMCALGeo,fEMCALCells,bc)) continue;
 269
 270     else if (c1->GetNCells() < fMinNCells)                        continue;
 271
 272     else if (c1->GetM02() < fL0min || c1->GetM02() > fL0max)      continue;
 273
 274     if(DebugLevel() > 2)
 275     {
 276       printf("IMA  : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",c1->GetID(),e1i,c1->GetDispersion(),c1->GetM02(),c1->GetM20());
 277       c1->GetPosition(pos);
 278       printf("IMA  : i %d, x %f, y %f, z %f\n",c1->GetID(), pos[0], pos[1], pos[2]);
 279     }
 280
 281     fRecoUtils->GetMaxEnergyCell(fEMCALGeo, fEMCALCells,c1,absId1,iSupMod1,ieta1,iphi1,shared);
 282     c1->GetMomentum(p1,fVertex);
 283
 284     //Check if cluster is in fidutial region, not too close to borders
 285     Bool_t in1 = fRecoUtils->CheckCellFiducialRegion(fEMCALGeo, c1, fEMCALCells);
 286
 287     // Clusters not facing frame structures
 288     Bool_t mask1 = MaskFrameCluster(iSupMod1, ieta1);
 289     //if(mask1) printf("Reject eta %d SM %d\n",ieta1, iSupMod1);
 290
 291     Double_t time1 = c1->GetTOF()*1.e9;
 292
 293     if(time1 > fTimeMax || time1 < fTimeMin) continue;
 294
 295     fhClusterTime            ->Fill(c1->E(),time1);
 296     fhClusterTimeSM[iSupMod1]->Fill(c1->E(),time1);
 297
 298     // Combine cluster with other clusters and get the invariant mass
 299     for (Int_t jClu=iClu+1; jClu < fCaloClustersArr->GetEntriesFast(); jClu++)
 300     {
 301       AliAODCaloCluster *c2 = (AliAODCaloCluster *) fCaloClustersArr->At(jClu);
 302
 303       Float_t e2i = c2->E();
 304       if      (e2i < fEmin) continue;
 305       else if (e2i > fEmax) continue;
 306
 307       else if (!fRecoUtils->IsGoodCluster(c2,fEMCALGeo,fEMCALCells,bc))continue;
 308
 309       else if (c2->GetNCells() < fMinNCells)                       continue;
 310
 311       else if (c2->GetM02() < fL0min || c2->GetM02() > fL0max)     continue;
 312
 313
 314       fRecoUtils->GetMaxEnergyCell(fEMCALGeo, fEMCALCells,c2,absId2,iSupMod2,ieta2,iphi2,shared);
 315       c2->GetMomentum(p2,fVertex);
 316
 317       p12 = p1+p2;
 318       Float_t invmass = p12.M()*1000;
 319
 320       //Asimetry cut
 321       Float_t asym = TMath::Abs(p1.E()-p2.E())/(p1.E()+p2.E());
 322
 323       if(asym > fAsyCut) continue;
 324
 325       //Time cut
 326       Double_t time2 = c2->GetTOF()*1.e9;
 327
 328       if(time2 > fTimeMax || time2 < fTimeMin) continue;
 329
 330       fhClusterPairDiffTime->Fill(p12.E(),time1-time2);
 331       if(TMath::Abs(time1-time2) > fDTimeCut) continue;
 332
 333       if(invmass < fMaxBin && invmass > fMinBin )
 334       {
 335         //Check if cluster is in fidutial region, not too close to borders
 336         Bool_t in2 = fRecoUtils->CheckCellFiducialRegion(fEMCALGeo, c2, fEMCALCells);
 337
 338         // Clusters not facing frame structures
 339         Bool_t mask2 = MaskFrameCluster(iSupMod2, ieta2);
 340         //if(mask2) printf("Reject eta %d SM %d\n",ieta2, iSupMod2);
 341
 342         if(in1 && in2)
 343         {
 344           fHmgg->Fill(invmass,p12.Pt());
 345
 346           if(iSupMod1==iSupMod2)
 347           {
 348             fHmggSM[iSupMod1]->Fill(invmass,p12.Pt());
 349             fhClusterPairDiffTimeSameSM[iSupMod1]->Fill(p12.E(),time1-time2);
 350           }
 351           else
 352             fHmggDifferentSM ->Fill(invmass,p12.Pt());
 353
 354           // Same sector
 355           Int_t j=0;
 356           for(Int_t i = 0; i < nSM/2; i++)
 357           {
 358             j=2*i;
 359             if((iSupMod1==j && iSupMod2==j+1) || (iSupMod1==j+1 && iSupMod2==j))
 360             {
 361               fHmggPairSameSectorSM[i]->Fill(invmass,p12.Pt());
 362               fhClusterPairDiffTimeSameSector[i]->Fill(p12.E(),time1-time2);
 363             }
 364           }
 365
 366           // Same side
 367           for(Int_t i = 0; i < nSM-2; i++)
 368           {
 369             if((iSupMod1==i && iSupMod2==i+2) || (iSupMod1==i+2 && iSupMod2==i))
 370             {
 371               fHmggPairSameSideSM[i]->Fill(invmass,p12.Pt());
 372               fhClusterPairDiffTimeSameSide[i]->Fill(p12.E(),time1-time2);
 373             }
 374           }
 375
 376
 377           if(!mask1 && !mask2)
 378           {
 379             fHmggMaskFrame->Fill(invmass,p12.Pt());
 380
 381             if(iSupMod1==iSupMod2) fHmggSMMaskFrame[iSupMod1]->Fill(invmass,p12.Pt());
 382             else                   fHmggDifferentSMMaskFrame ->Fill(invmass,p12.Pt());
 383
 384             // Same sector
 385             j=0;
 386             for(Int_t i = 0; i < nSM/2; i++)
 387             {
 388               j=2*i;
 389               if((iSupMod1==j && iSupMod2==j+1) || (iSupMod1==j+1 && iSupMod2==j)) fHmggPairSameSectorSMMaskFrame[i]->Fill(invmass,p12.Pt());
 390             }
 391
 392             // Same side
 393             for(Int_t i = 0; i < nSM-2; i++)
 394             {
 395               if((iSupMod1==i && iSupMod2==i+2) || (iSupMod1==i+2 && iSupMod2==i)) fHmggPairSameSideSMMaskFrame[i]->Fill(invmass,p12.Pt());
 396             }
 397
 398           }// Pair not facing frame
 399
 400
 401           if(invmass > fInvMassCutMin && invmass < fInvMassCutMax) //restrict to clusters really close to pi0 peak
 402           {
 403
 404             // Check time of cells in both clusters, and fill time histogram
 405             for(Int_t icell = 0; icell < c1->GetNCells(); icell++)
 406             {
 407               Int_t absID = c1->GetCellAbsId(icell);
 408               fHTpi0[bc%4]->Fill(absID, fEMCALCells->GetCellTime(absID)*1.e9);
 409             }
 410
 411             for(Int_t icell = 0; icell < c2->GetNCells(); icell++)
 412             {
 413               Int_t absID = c2->GetCellAbsId(icell);
 414               fHTpi0[bc%4]->Fill(absID, fEMCALCells->GetCellTime(absID)*1.e9);
 415             }
 416
 417             //Opening angle of 2 photons
 418             Float_t opangle = p1.Angle(p2.Vect())*TMath::RadToDeg();
 419             //printf("*******>>>>>>>> In PEAK pt %f, angle %f \n",p12.Pt(),opangle);
 420
 421
 422             fHOpeningAngle ->Fill(opangle,p12.Pt());
 423             fHAsymmetry    ->Fill(asym,p12.Pt());
 424
 425             if(iSupMod1==iSupMod2)
 426             {
 427               fHOpeningAngleSM[iSupMod1] ->Fill(opangle,p12.Pt());
 428               fHAsymmetrySM[iSupMod1]    ->Fill(asym,p12.Pt());
 429             }
 430             else
 431             {
 432               fHOpeningAngleDifferentSM  ->Fill(opangle,p12.Pt());
 433               fHAsymmetryDifferentSM     ->Fill(asym,p12.Pt());
 434             }
 435
 436             if((iSupMod1==0 && iSupMod2==2) || (iSupMod1==2 && iSupMod2==0))
 437             {
 438               fHOpeningAnglePairSM[0] ->Fill(opangle,p12.Pt());
 439               fHAsymmetryPairSM[0]    ->Fill(asym,p12.Pt());
 440
 441             }
 442             if((iSupMod1==1 && iSupMod2==3) || (iSupMod1==3 && iSupMod2==1))
 443             {
 444               fHOpeningAnglePairSM[1] ->Fill(opangle,p12.Pt());
 445               fHAsymmetryPairSM[1]    ->Fill(asym,p12.Pt());
 446             }
 447
 448             if((iSupMod1==0 && iSupMod2==1) || (iSupMod1==1 && iSupMod2==0))
 449             {
 450               fHOpeningAnglePairSM[2] ->Fill(opangle,p12.Pt());
 451               fHAsymmetryPairSM[2]    ->Fill(asym,p12.Pt());
 452             }
 453             if((iSupMod1==2 && iSupMod2==3) || (iSupMod1==3 && iSupMod2==2))
 454             {
 455               fHOpeningAnglePairSM[3] ->Fill(opangle,p12.Pt());
 456               fHAsymmetryPairSM[3]    ->Fill(asym,p12.Pt());
 457             }
 458
 459           }// pair in 100 < mass < 160
 460
 461         }//in acceptance cuts
 462
 463         //In case of filling only channels with second cluster in same SM
 464         if(fSameSM && iSupMod1!=iSupMod2) continue;
 465
 466         if (fGroupNCells == 0)
 467         {
 468           fHmpi0[iSupMod1][ieta1][iphi1]->Fill(invmass);
 469           fHmpi0[iSupMod2][ieta2][iphi2]->Fill(invmass);
 470
 471           if(invmass > fInvMassCutMin && invmass < fInvMassCutMax)//restrict to clusters really close to pi0 peak
 472           {
 473             fhTowerDecayPhotonHit      [iSupMod1]->Fill(ieta1,iphi1);
 474             fhTowerDecayPhotonEnergy   [iSupMod1]->Fill(ieta1,iphi1,p1.E());
 475             fhTowerDecayPhotonAsymmetry[iSupMod1]->Fill(ieta1,iphi1,asym);
 476
 477             fhTowerDecayPhotonHit      [iSupMod2]->Fill(ieta2,iphi2);
 478             fhTowerDecayPhotonEnergy   [iSupMod2]->Fill(ieta2,iphi2,p2.E());
 479             fhTowerDecayPhotonAsymmetry[iSupMod2]->Fill(ieta2,iphi2,asym);
 480
 481             if(!mask1)fhTowerDecayPhotonHitMaskFrame[iSupMod1]->Fill(ieta1,iphi1);
 482             if(!mask2)fhTowerDecayPhotonHitMaskFrame[iSupMod2]->Fill(ieta2,iphi2);
 483
 484           }// pair in mass of pi0
 485         }
 486         else  {
 487           //printf("Regroup N %d, eta1 %d, phi1 %d, eta2 %d, phi2 %d \n",fGroupNCells, ieta1, iphi1, ieta2, iphi2);
 488           for (Int_t i = -fGroupNCells; i < fGroupNCells+1; i++)
 489           {
 490             for (Int_t j = -fGroupNCells; j < fGroupNCells+1; j++)
 491             {
 492               Int_t absId11 = fEMCALGeo->GetAbsCellIdFromCellIndexes(iSupMod1, iphi1+j, ieta1+i);
 493               Int_t absId22 = fEMCALGeo->GetAbsCellIdFromCellIndexes(iSupMod2, iphi2+j, ieta2+i);
 494               Bool_t ok1 = kFALSE;
 495               Bool_t ok2 = kFALSE;
 496               for(Int_t icell = 0; icell < c1->GetNCells(); icell++){
 497                 if(c1->GetCellsAbsId()[icell] == absId11) ok1=kTRUE;
 498               }
 499               for(Int_t icell = 0; icell < c2->GetNCells(); icell++){
 500                 if(c2->GetCellsAbsId()[icell] == absId22) ok2=kTRUE;
 501               }
 502
 503               if(ok1 && (ieta1+i >= 0) && (iphi1+j >= 0) && (ieta1+i < 48) && (iphi1+j < 24))
 504               {
 505                 fHmpi0[iSupMod1][ieta1+i][iphi1+j]->Fill(invmass);
 506               }
 507               if(ok2 && (ieta2+i >= 0) && (iphi2+j >= 0) && (ieta2+i < 48) && (iphi2+j < 24))
 508               {
 509                 fHmpi0[iSupMod2][ieta2+i][iphi2+j]->Fill(invmass);
 510               }
 511             }// j loop
 512           }//i loop
 513         }//group cells
 514
 515         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",
 516                                     iSupMod1,iphi1,ieta1,iSupMod2,iphi2,ieta2,p12.M(),e1i,c1->E(),e2i,c2->E());
 517       }
 518
 519     }
 520
 521   } // end of loop over EMCAL clusters
 522 }
 523
 524 //________________________________________________________________
 525 void AliAnalysisTaskEMCALPi0CalibSelection::InitGeometryMatrices()
 526 {
 527   // Init geometry and set the geometry matrix, for the first event, skip the rest
 528   // Also set once the run dependent calibrations
 529
 530
 531   Int_t runnumber = InputEvent()->GetRunNumber() ;
 532
 533   if(fLoadMatrices)
 534   {
 535     printf("AliAnalysisTaskEMCALPi0CalibSelection::InitGeometryMatrices() - Load user defined EMCAL geometry matrices\n");
 536
 537     // OADB if available
 538     AliOADBContainer emcGeoMat("AliEMCALgeo");
 539
 540     if(fOADBFilePath=="") fOADBFilePath = "$ALICE_ROOT/OADB/EMCAL" ;
 541
 542     emcGeoMat.InitFromFile(Form("%s/EMCALlocal2master.root",fOADBFilePath.Data()),"AliEMCALgeo");
 543
 544     TObjArray *matEMCAL=(TObjArray*)emcGeoMat.GetObject(runnumber,"EmcalMatrices");
 545
 546     for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
 547     {
 548
 549       if (!fMatrix[mod]) // Get it from OADB
 550       {
 551         if(fDebug > 1 )
 552           printf("AliAnalysisTaskEMCALPi0CalibSelection::InitGeometryMatrices() - EMCAL matrices SM %d, %p\n",
 553                  mod,((TGeoHMatrix*) matEMCAL->At(mod)));
 554         //((TGeoHMatrix*) matEMCAL->At(mod))->Print();
 555
 556         fMatrix[mod] = (TGeoHMatrix*) matEMCAL->At(mod) ;
 557       }
 558
 559       if(fMatrix[mod])
 560       {
 561         if(DebugLevel() > 1)
 562           fMatrix[mod]->Print();
 563
 564         fEMCALGeo->SetMisalMatrix(fMatrix[mod],mod) ;
 565       }
 566
 567     }//SM loop
 568   }//Load matrices
 569   else if(!gGeoManager)
 570   {
 571     printf("AliAnalysisTaskEMCALPi0CalibSelection::InitGeometryMatrices() - Get geo matrices from data");
 572     //Still not implemented in AOD, just a workaround to be able to work at least with ESDs
 573     if(!strcmp(InputEvent()->GetName(),"AliAODEvent"))
 574     {
 575       if(DebugLevel() > 1)
 576         Warning("UserExec","Use ideal geometry, values geometry matrix not kept in AODs.");
 577     }//AOD
 578     else
 579     {
 580       if(DebugLevel() > 1)
 581         printf("AliAnalysisTaskEMCALPi0CalibSelection::InitGeometryMatrices() - AliAnalysisTaskEMCALClusterize Load Misaligned matrices.");
 582
 583       for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
 584       {
 585         if(DebugLevel() > 1)
 586           InputEvent()->GetEMCALMatrix(mod)->Print();
 587
 588         if(InputEvent()->GetEMCALMatrix(mod)) fEMCALGeo->SetMisalMatrix(InputEvent()->GetEMCALMatrix(mod),mod) ;
 589
 590       }
 591
 592     }//ESD
 593   }//Load matrices from Data
 594
 595 }
 596
 597 //______________________________________________________________________
 598 void AliAnalysisTaskEMCALPi0CalibSelection::InitTemperatureCorrections()
 599 {
 600   // Apply run dependent calibration correction
 601
 602   if(!fRecoUtils->IsRunDepRecalibrationOn()) return;
 603
 604   AliOADBContainer *contRFTD=new AliOADBContainer("");
 605
 606   contRFTD->InitFromFile(Form("%s/EMCALTemperatureCorrCalib.root",fOADBFilePath.Data()),"AliEMCALRunDepTempCalibCorrections");
 607
 608   Int_t runnumber = InputEvent()->GetRunNumber() ;
 609
 610   TH1S *htd=(TH1S*)contRFTD->GetObject(runnumber);
 611
 612   //If it did not exist for this run, get closes one
 613   if (!htd)
 614   {
 615     AliWarning(Form("No TemperatureCorrCalib Objects for run: %d",runnumber));
 616     // let's get the closest runnumber instead then..
 617     Int_t lower = 0;
 618     Int_t ic = 0;
 619     Int_t maxEntry = contRFTD->GetNumberOfEntries();
 620
 621     while ( (ic < maxEntry) && (contRFTD->UpperLimit(ic) < runnumber) ) {
 622       lower = ic;
 623       ic++;
 624     }
 625
 626     Int_t closest = lower;
 627     if ( (ic<maxEntry) &&
 628         (contRFTD->LowerLimit(ic)-runnumber) < (runnumber - contRFTD->UpperLimit(lower)) ) {
 629       closest = ic;
 630     }
 631
 632     AliWarning(Form("TemperatureCorrCalib Objects found closest id %d from run: %d", closest, contRFTD->LowerLimit(closest)));
 633     htd = (TH1S*) contRFTD->GetObjectByIndex(closest);
 634   }
 635
 636   // Fill parameters
 637   if(htd)
 638   {
 639     printf("AliAnalysisTaskEMCALPi0CalibSelection::SetOADBParameters() - Recalibrate (Temperature) EMCAL \n");
 640
 641     Int_t nSM = fEMCALGeo->GetNumberOfSuperModules();
 642
 643     for (Int_t ism = 0; ism < nSM; ++ism)
 644     {
 645       for (Int_t icol = 0; icol < 48; ++icol)
 646       {
 647         for (Int_t irow = 0; irow < 24; ++irow)
 648         {
 649           Float_t factor = fRecoUtils->GetEMCALChannelRecalibrationFactor(ism,icol,irow);
 650
 651           Int_t absID = fEMCALGeo->GetAbsCellIdFromCellIndexes(ism, irow, icol); // original calibration factor
 652
 653           if(DebugLevel() > 3)
 654             printf(" ism %d, icol %d, irow %d,absID %d - Calib factor %1.5f - ",ism, icol, irow, absID, factor);
 655
 656           factor *= htd->GetBinContent(absID) / 10000. ; // correction dependent on T
 657
 658           fRecoUtils->SetEMCALChannelRecalibrationFactor(ism,icol,irow,factor);
 659
 660           if(DebugLevel() > 3)
 661             printf(" T  factor %1.5f - final factor %1.5f \n",htd->GetBinContent(absID) / 10000.,
 662                    fRecoUtils->GetEMCALChannelRecalibrationFactor(ism,icol,irow));
 663
 664         } // columns
 665       } // rows
 666     } // SM loop
 667   }else printf("AliAnalysisTaskEMCALPi0CalibSelection::SetOADBParameters() - Do NOT recalibrate EMCAL with T variations, no params TH1 \n");
 668
 669 }
 670
 671 //___________________________________________________________________
 672 void AliAnalysisTaskEMCALPi0CalibSelection::UserCreateOutputObjects()
 673 {
 674   //Create output container, init geometry
 675
 676   fEMCALGeo =  AliEMCALGeometry::GetInstance(fEMCALGeoName) ;
 677   Int_t nSM = (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules();
 678
 679   fOutputContainer = new TList();
 680   const Int_t buffersize = 255;
 681   char hname[buffersize], htitl[buffersize];
 682
 683   fhNEvents        = new TH1I("hNEvents", "Number of analyzed events"   , 1 , 0 , 1  ) ;
 684   fOutputContainer->Add(fhNEvents);
 685
 686   fHmgg = new TH2F("hmgg","2-cluster invariant mass",fNbins,fMinBin,fMaxBin,100,0,10);
 687   fHmgg->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
 688   fHmgg->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 689   fOutputContainer->Add(fHmgg);
 690
 691   fHmggDifferentSM = new TH2F("hmggDifferentSM","2-cluster invariant mass, different SM",fNbins,fMinBin,fMaxBin,100,0,10);
 692   fHmggDifferentSM->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
 693   fHmggDifferentSM->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 694   fOutputContainer->Add(fHmggDifferentSM);
 695
 696   fHOpeningAngle = new TH2F("hopang","2-cluster opening angle",100,0.,50.,100,0,10);
 697   fHOpeningAngle->SetXTitle("#alpha_{#gamma #gamma}");
 698   fHOpeningAngle->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 699   fOutputContainer->Add(fHOpeningAngle);
 700
 701   fHOpeningAngleDifferentSM = new TH2F("hopangDifferentSM","2-cluster opening angle, different SM",100,0,50.,100,0,10);
 702   fHOpeningAngleDifferentSM->SetXTitle("#alpha_{#gamma #gamma}");
 703   fHOpeningAngleDifferentSM->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 704   fOutputContainer->Add(fHOpeningAngleDifferentSM);
 705
 706   fHAsymmetry = new TH2F("hasym","2-cluster opening angle",100,0.,1.,100,0,10);
 707   fHAsymmetry->SetXTitle("a");
 708   fHAsymmetry->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 709   fOutputContainer->Add(fHAsymmetry);
 710
 711   fHAsymmetryDifferentSM = new TH2F("hasymDifferentSM","2-cluster opening angle, different SM",100,0,1.,100,0,10);
 712   fHAsymmetryDifferentSM->SetXTitle("a");
 713   fHAsymmetryDifferentSM->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 714   fOutputContainer->Add(fHAsymmetryDifferentSM);
 715
 716
 717   //TString pairname[] = {"A side (0-2)", "C side (1-3)","Row 0 (0-1)", "Row 1 (2-3)"};
 718
 719   fHmggMaskFrame = new TH2F("hmggMaskFrame","2-cluster invariant mass, frame masked",fNbins,fMinBin,fMaxBin,100,0,10);
 720   fHmggMaskFrame->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
 721   fHmggMaskFrame->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 722   fOutputContainer->Add(fHmggMaskFrame);
 723
 724   fHmggDifferentSMMaskFrame = new TH2F("hmggDifferentSMMaskFrame","2-cluster invariant mass, different SM, frame masked",
 725                                        fNbins,fMinBin,fMaxBin,100,0,10);
 726   fHmggDifferentSMMaskFrame->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
 727   fHmggDifferentSMMaskFrame->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 728   fOutputContainer->Add(fHmggDifferentSMMaskFrame);
 729
 730
 731   for(Int_t iSM = 0; iSM < nSM; iSM++)
 732   {
 733     snprintf(hname, buffersize, "hmgg_SM%d",iSM);
 734     snprintf(htitl, buffersize, "Two-gamma inv. mass for super mod %d",iSM);
 735     fHmggSM[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
 736     fHmggSM[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
 737     fHmggSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 738     fOutputContainer->Add(fHmggSM[iSM]);
 739
 740     snprintf(hname, buffersize, "hmgg_SM%d_MaskFrame",iSM);
 741     snprintf(htitl, buffersize, "Two-gamma inv. mass for super mod %d",iSM);
 742     fHmggSMMaskFrame[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
 743     fHmggSMMaskFrame[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
 744     fHmggSMMaskFrame[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 745     fOutputContainer->Add(fHmggSMMaskFrame[iSM]);
 746
 747
 748     if(iSM < nSM/2)
 749     {
 750       snprintf(hname,buffersize, "hmgg_PairSameSectorSM%d",iSM);
 751       snprintf(htitl,buffersize, "Two-gamma inv. mass for SM pair Sector: %d",iSM);
 752       fHmggPairSameSectorSM[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
 753       fHmggPairSameSectorSM[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
 754       fHmggPairSameSectorSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 755       fOutputContainer->Add(fHmggPairSameSectorSM[iSM]);
 756
 757       snprintf(hname,buffersize, "hmgg_PairSameSectorSM%d_MaskFrame",iSM);
 758       snprintf(htitl,buffersize, "Two-gamma inv. mass for SM pair Sector: %d",iSM);
 759       fHmggPairSameSectorSMMaskFrame[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
 760       fHmggPairSameSectorSMMaskFrame[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
 761       fHmggPairSameSectorSMMaskFrame[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 762       fOutputContainer->Add(fHmggPairSameSectorSMMaskFrame[iSM]);
 763
 764       fhClusterPairDiffTimeSameSector[iSM] = new TH2F(Form("hClusterPairDiffTimeSameSector%d",iSM),
 765                                                       Form("cluster pair time difference vs E, Sector %d",iSM),
 766                                                       100,0,10, 200,-100,100);
 767       fhClusterPairDiffTimeSameSector[iSM]->SetXTitle("E_{pair} (GeV)");
 768       fhClusterPairDiffTimeSameSector[iSM]->SetYTitle("#Delta t (ns)");
 769       fOutputContainer->Add(fhClusterPairDiffTimeSameSector[iSM]);
 770
 771
 772     }
 773
 774     if(iSM < nSM-2)
 775     {
 776       snprintf(hname,buffersize, "hmgg_PairSameSideSM%d",iSM);
 777       snprintf(htitl,buffersize, "Two-gamma inv. mass for SM pair Sector: %d",iSM);
 778       fHmggPairSameSideSM[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
 779       fHmggPairSameSideSM[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
 780       fHmggPairSameSideSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 781       fOutputContainer->Add(fHmggPairSameSideSM[iSM]);
 782
 783       snprintf(hname,buffersize, "hmgg_PairSameSideSM%d_MaskFrame",iSM);
 784       snprintf(htitl,buffersize, "Two-gamma inv. mass for SM pair Sector: %d",iSM);
 785       fHmggPairSameSideSMMaskFrame[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
 786       fHmggPairSameSideSMMaskFrame[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
 787       fHmggPairSameSideSMMaskFrame[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 788       fOutputContainer->Add(fHmggPairSameSideSMMaskFrame[iSM]);
 789
 790       fhClusterPairDiffTimeSameSide[iSM] = new TH2F(Form("hClusterPairDiffTimeSameSide%d",iSM),
 791                                                     Form("cluster pair time difference vs E,  Side %d",iSM),
 792                                                     100,0,10, 200,-100,100);
 793       fhClusterPairDiffTimeSameSide[iSM]->SetXTitle("E_{pair} (GeV)");
 794       fhClusterPairDiffTimeSameSide[iSM]->SetYTitle("#Delta t (ns)");
 795       fOutputContainer->Add(fhClusterPairDiffTimeSameSide[iSM]);
 796
 797     }
 798
 799     snprintf(hname, buffersize, "hopang_SM%d",iSM);
 800     snprintf(htitl, buffersize, "Opening angle for super mod %d",iSM);
 801     fHOpeningAngleSM[iSM] = new TH2F(hname,htitl,100,0.,50.,100,0,10);
 802     fHOpeningAngleSM[iSM]->SetXTitle("#alpha_{#gamma #gamma} (deg)");
 803     fHOpeningAngleSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 804     fOutputContainer->Add(fHOpeningAngleSM[iSM]);
 805
 806     snprintf(hname,buffersize, "hopang_PairSM%d",iSM);
 807     snprintf(htitl,buffersize, "Opening angle for SM pair: %d",iSM);
 808     fHOpeningAnglePairSM[iSM] = new TH2F(hname,htitl,100,0.,50.,100,0,10);
 809     fHOpeningAnglePairSM[iSM]->SetXTitle("#alpha_{#gamma #gamma} (deg)");
 810     fHOpeningAnglePairSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 811     fOutputContainer->Add(fHOpeningAnglePairSM[iSM]);
 812
 813     snprintf(hname, buffersize, "hasym_SM%d",iSM);
 814     snprintf(htitl, buffersize, "Asymmetry for super mod %d",iSM);
 815     fHAsymmetrySM[iSM] = new TH2F(hname,htitl,100,0.,1.,100,0,10);
 816     fHAsymmetrySM[iSM]->SetXTitle("a");
 817     fHAsymmetrySM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 818     fOutputContainer->Add(fHAsymmetrySM[iSM]);
 819
 820     snprintf(hname,buffersize, "hasym_PairSM%d",iSM);
 821     snprintf(htitl,buffersize, "Asymmetry for SM pair: %d",iSM);
 822     fHAsymmetryPairSM[iSM] = new TH2F(hname,htitl,100,0.,1.,100,0,10);
 823     fHAsymmetryPairSM[iSM]->SetXTitle("a");
 824     fHAsymmetryPairSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 825     fOutputContainer->Add(fHAsymmetryPairSM[iSM]);
 826
 827     Int_t colmax = 48;
 828     Int_t rowmax = 24;
 829
 830     fhTowerDecayPhotonHit[iSM]  = new TH2F (Form("hTowerDecPhotonHit_Mod%d",iSM),
 831                                             Form("Entries in grid of cells in Module %d",iSM),
 832                                             colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5);
 833     fhTowerDecayPhotonHit[iSM]->SetYTitle("row (phi direction)");
 834     fhTowerDecayPhotonHit[iSM]->SetXTitle("column (eta direction)");
 835     fOutputContainer->Add(fhTowerDecayPhotonHit[iSM]);
 836
 837     fhTowerDecayPhotonEnergy[iSM]  = new TH2F (Form("hTowerDecPhotonEnergy_Mod%d",iSM),
 838                                                Form("Accumulated energy in grid of cells in Module %d",iSM),
 839                                                colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5);
 840     fhTowerDecayPhotonEnergy[iSM]->SetYTitle("row (phi direction)");
 841     fhTowerDecayPhotonEnergy[iSM]->SetXTitle("column (eta direction)");
 842     fOutputContainer->Add(fhTowerDecayPhotonEnergy[iSM]);
 843
 844     fhTowerDecayPhotonAsymmetry[iSM]  = new TH2F (Form("hTowerDecPhotonAsymmetry_Mod%d",iSM),
 845                                                   Form("Accumulated asymmetry in grid of cells in Module %d",iSM),
 846                                                   colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5);
 847     fhTowerDecayPhotonAsymmetry[iSM]->SetYTitle("row (phi direction)");
 848     fhTowerDecayPhotonAsymmetry[iSM]->SetXTitle("column (eta direction)");
 849     fOutputContainer->Add(fhTowerDecayPhotonAsymmetry[iSM]);
 850
 851     fhTowerDecayPhotonHitMaskFrame[iSM]  = new TH2F (Form("hTowerDecPhotonHit_Mod%d_MaskFrame",iSM),Form("Entries in grid of cells in Module %d",iSM),
 852                                                      colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5);
 853     fhTowerDecayPhotonHitMaskFrame[iSM]->SetYTitle("row (phi direction)");
 854     fhTowerDecayPhotonHitMaskFrame[iSM]->SetXTitle("column (eta direction)");
 855     fOutputContainer->Add(fhTowerDecayPhotonHitMaskFrame[iSM]);
 856
 857     fhClusterTimeSM[iSM] = new TH2F(Form("hClusterTime_SM%d",iSM),"cluster time vs E",100,0,10, 100,0,1000);
 858     fhClusterTimeSM[iSM]->SetXTitle("E (GeV)");
 859     fhClusterTimeSM[iSM]->SetYTitle("t (ns)");
 860     fOutputContainer->Add(fhClusterTimeSM[iSM]);
 861
 862     fhClusterPairDiffTimeSameSM[iSM] = new TH2F(Form("hClusterPairDiffTimeSameSM%d",iSM),
 863                                                 Form("cluster pair time difference vs E, SM %d",iSM),
 864                                                 100,0,10, 200,-100,100);
 865     fhClusterPairDiffTimeSameSM[iSM]->SetXTitle("E (GeV)");
 866     fhClusterPairDiffTimeSameSM[iSM]->SetYTitle("#Delta t (ns)");
 867     fOutputContainer->Add(fhClusterPairDiffTimeSameSM[iSM]);
 868
 869   }
 870
 871   Int_t nchannels = nSM*AliEMCALGeoParams::fgkEMCALRows*AliEMCALGeoParams::fgkEMCALCols;
 872   for(Int_t ibc = 0; ibc < 4; ibc++)
 873   {
 874     fHTpi0[ibc] = new TH2F(Form("hTime_BC%d",ibc),Form("Time of cell clusters under pi0 peak, bunch crossing %d",ibc),
 875                            nchannels,0,nchannels, fNTimeBins,fMinTimeBin,fMaxTimeBin);
 876     fOutputContainer->Add(fHTpi0[ibc]);
 877     fHTpi0[ibc]->SetYTitle("time (ns)");
 878     fHTpi0[ibc]->SetXTitle("abs. Id. ");
 879   }
 880
 881
 882   fhClusterTime = new TH2F("hClusterTime","cluster time vs E",100,0,10, 100,0,1000);
 883   fhClusterTime->SetXTitle("E (GeV)");
 884   fhClusterTime->SetYTitle("t (ns)");
 885   fOutputContainer->Add(fhClusterTime);
 886
 887   fhClusterPairDiffTime = new TH2F("hClusterPairDiffTime","cluster pair time difference vs E",100,0,10, 800,-400,400);
 888   fhClusterPairDiffTime->SetXTitle("E_{pair} (GeV)");
 889   fhClusterPairDiffTime->SetYTitle("#Delta t (ns)");
 890   fOutputContainer->Add(fhClusterPairDiffTime);
 891
 892   for(Int_t iMod=0; iMod < nSM; iMod++)
 893   {
 894     for(Int_t iRow=0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++)
 895     {
 896       for(Int_t iCol=0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++)
 897       {
 898         snprintf(hname,buffersize, "%d_%d_%d",iMod,iCol,iRow);
 899         snprintf(htitl,buffersize, "Two-gamma inv. mass for super mod %d, cell(col,row)=(%d,%d)",iMod,iCol,iRow);
 900         fHmpi0[iMod][iCol][iRow] = new TH1F(hname,htitl,fNbins,fMinBin,fMaxBin);
 901         fHmpi0[iMod][iCol][iRow]->SetXTitle("mass (MeV/c^{2})");
 902         fOutputContainer->Add(fHmpi0[iMod][iCol][iRow]);
 903       }
 904     }
 905   }
 906
 907   fOutputContainer->SetOwner(kTRUE);
 908
 909   PostData(1,fOutputContainer);
 910
 911   // cuts container, set in terminate but init and post here
 912
 913   fCuts = new TList();
 914
 915   fCuts ->SetOwner(kTRUE);
 916
 917   PostData(2, fCuts);
 918
 919 }
 920
 921 //______________________________________________________________________________________________________
 922 Bool_t AliAnalysisTaskEMCALPi0CalibSelection::MaskFrameCluster(const Int_t iSM,  const Int_t ieta) const
 923 {
 924   //Check if cell is in one of the regions where we have significant amount of material in front of EMCAL
 925
 926   Int_t icol = ieta;
 927   if(iSM%2) icol+=48; // Impair SM, shift index [0-47] to [48-96]
 928
 929   if (fNMaskCellColumns && fMaskCellColumns)
 930   {
 931     for (Int_t imask = 0; imask < fNMaskCellColumns; imask++)
 932     {
 933       if(icol==fMaskCellColumns[imask]) return kTRUE;
 934     }
 935   }
 936
 937   return kFALSE;
 938
 939 }
 940
 941 //__________________________________________________________________________
 942 void AliAnalysisTaskEMCALPi0CalibSelection::UserExec(Option_t* /* option */)
 943 {
 944   // Do analysis, first select the events, then correct the clusters if needed
 945   // and finally fill the histograms per channel after recalibration
 946
 947   //Event selection
 948   if(fTriggerName!="")
 949   {
 950     AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (InputEvent());
 951     AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (InputEvent());
 952
 953     TString triggerClass = "";
 954     if     (esdevent) triggerClass = esdevent->GetFiredTriggerClasses();
 955     else if(aodevent) triggerClass = aodevent->GetFiredTriggerClasses();
 956
 957     if(DebugLevel() > 1)
 958       printf("AliAnalysisTaskEMCALPi0CalibSelection::UserExec() - Event %d, FiredClass %s",
 959              (Int_t)Entry(),(((AliESDEvent*)InputEvent())->GetFiredTriggerClasses()).Data());
 960
 961     if(!triggerClass.Contains(fTriggerName))
 962     {
 963        if(DebugLevel() > 1) printf("Reject event! \n");
 964        return;
 965     }
 966     else
 967        if(DebugLevel() > 1) printf("Accept Event! \n");
 968   }
 969
 970   //Get the input event
 971   AliVEvent* event = 0;
 972   if(fFilteredInput) event = AODEvent();
 973   else               event = InputEvent();
 974
 975   if(!event)
 976   {
 977     printf("Input event not available!\n");
 978     return;
 979   }
 980
 981   if(DebugLevel() > 1)
 982     printf("AliAnalysisTaskEMCALPi0CalibSelection <<< %s: Event %d >>>\n",event->GetName(), (Int_t)Entry());
 983
 984   //Get the primary vertex
 985   event->GetPrimaryVertex()->GetXYZ(fVertex) ;
 986
 987   if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Vertex: (%.3f,%.3f,%.3f)\n",fVertex[0],fVertex[1],fVertex[2]);
 988
 989   //Int_t runNum = aod->GetRunNumber();
 990   //if(DebugLevel() > 1) printf("Run number: %d\n",runNum);
 991
 992   fhNEvents->Fill(0); //Count the events to be analyzed
 993
 994   // Acccess once the geometry matrix and temperature corrections
 995   if(fhNEvents->GetEntries()==1)
 996   {
 997     InitGeometryMatrices();
 998
 999     InitTemperatureCorrections();
1000   }
1001
1002   //Get the list of clusters and cells
1003   fEMCALCells       = event->GetEMCALCells();
1004
1005   fCaloClustersArr  = new TRefArray();
1006   event->GetEMCALClusters(fCaloClustersArr);
1007
1008   if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection - N CaloClusters: %d - N CaloCells %d \n",
1009                               fCaloClustersArr->GetEntriesFast(), fEMCALCells->GetNumberOfCells());
1010
1011   CorrectClusters(); // Non linearity, new calibration, T calibration
1012
1013   FillHistograms();
1014
1015   delete fCaloClustersArr;
1016
1017   PostData(1,fOutputContainer);
1018
1019 }
1020
1021 //_____________________________________________________
1022 void AliAnalysisTaskEMCALPi0CalibSelection::PrintInfo()
1023 {
1024   //Print settings
1025
1026   printf("Cluster cuts: %2.2f < E < %2.2f GeV; number of cells > %d; Assymetry < %1.2f, pair time diff < %2.2f, %2.2f < t < %2.2f ns\n",
1027          fEmin,fEmax, fMinNCells, fAsyCut, fDTimeCut,fTimeMin,fTimeMax) ;
1028
1029   printf("Group %d cells\n", fGroupNCells) ;
1030
1031   printf("Cluster maximal cell away from border at least %d cells\n", fRecoUtils->GetNumberOfCellsFromEMCALBorder()) ;
1032
1033   printf("Histograms: bins %d; energy range: %2.2f < E < %2.2f GeV\n",fNbins,fMinBin,fMaxBin) ;
1034
1035   printf("Switchs:\n \t Remove Bad Channels? %d; Use filtered input? %d;  Correct Clusters? %d, \n \t Mass per channel same SM clusters? %d\n",
1036          fRecoUtils->IsBadChannelsRemovalSwitchedOn(),fFilteredInput,fCorrectClusters, fSameSM) ;
1037
1038   printf("EMCAL Geometry name: < %s >, Load Matrices %d\n",fEMCALGeoName.Data(), fLoadMatrices) ;
1039   if(fLoadMatrices) {for(Int_t ism = 0; ism < AliEMCALGeoParams::fgkEMCALModules; ism++) if(fMatrix[ism]) fMatrix[ism]->Print() ; }
1040
1041 }
1042
1043 //____________________________________________________________________
1044 void AliAnalysisTaskEMCALPi0CalibSelection::Terminate(Option_t*)
1045 {
1046   // Create cuts/param objects and publish to slot
1047   const Int_t buffersize = 255;
1048   char onePar[buffersize] ;
1049
1050   snprintf(onePar,buffersize, "Custer cuts: %2.2f < E < %2.2f GeV; %2.2f < Lambda0_2 < %2.2f GeV; min number of cells %d; Assymetry cut %1.2f, time1-time2 < %2.2f; %2.2f < T < %2.2f ns; %3.1f < Mass < %3.1f",
1051            fEmin,fEmax, fL0min, fL0max, fMinNCells, fAsyCut, fDTimeCut, fTimeMin, fTimeMax, fInvMassCutMin, fInvMassCutMax) ;
1052   fCuts->Add(new TObjString(onePar));
1053   snprintf(onePar,buffersize, "Group %d cells;", fGroupNCells) ;
1054   fCuts->Add(new TObjString(onePar));
1055   snprintf(onePar,buffersize, "Cluster maximal cell away from border at least %d cells;", fRecoUtils->GetNumberOfCellsFromEMCALBorder()) ;
1056   fCuts->Add(new TObjString(onePar));
1057   snprintf(onePar,buffersize, "Histograms, Mass bins %d; energy range: %2.2f < E < %2.2f GeV;",fNbins,fMinBin,fMaxBin) ;
1058   fCuts->Add(new TObjString(onePar));
1059   snprintf(onePar,buffersize, "Histograms, Time bins %d; energy range: %2.2f < E < %2.2f GeV;",fNTimeBins,fMinTimeBin,fMaxTimeBin) ;
1060   fCuts->Add(new TObjString(onePar));
1061   snprintf(onePar,buffersize, "Switchs: Remove Bad Channels? %d; Use filtered input? %d;  Correct Clusters? %d, Mass per channel same SM clusters? %d ",
1062            fRecoUtils->IsBadChannelsRemovalSwitchedOn(),fFilteredInput,fCorrectClusters, fSameSM) ;
1063   fCuts->Add(new TObjString(onePar));
1064   snprintf(onePar,buffersize, "EMCAL Geometry name: < %s >, Load Matrices? %d",fEMCALGeoName.Data(),fLoadMatrices) ;
1065   fCuts->Add(new TObjString(onePar));
1066
1067   // Post Data
1068   PostData(2, fCuts);
1069
1070 }
1071