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(htd)
 613   {
 614     printf("AliAnalysisTaskEMCALPi0CalibSelection::SetOADBParameters() - Recalibrate (Temperature) EMCAL \n");
 615
 616     Int_t nSM = fEMCALGeo->GetNumberOfSuperModules();
 617
 618     for (Int_t ism = 0; ism < nSM; ++ism)
 619     {
 620       for (Int_t icol = 0; icol < 48; ++icol)
 621       {
 622         for (Int_t irow = 0; irow < 24; ++irow)
 623         {
 624           Float_t factor = fRecoUtils->GetEMCALChannelRecalibrationFactor(ism,icol,irow);
 625
 626           Int_t absID = fEMCALGeo->GetAbsCellIdFromCellIndexes(ism, irow, icol); // original calibration factor
 627
 628           if(DebugLevel() > 3)
 629             printf(" ism %d, icol %d, irow %d,absID %d - Calib factor %1.5f - ",ism, icol, irow, absID, factor);
 630
 631           factor *= htd->GetBinContent(absID) / 10000. ; // correction dependent on T
 632
 633           fRecoUtils->SetEMCALChannelRecalibrationFactor(ism,icol,irow,factor);
 634
 635           if(DebugLevel() > 3)
 636             printf(" T  factor %1.5f - final factor %1.5f \n",htd->GetBinContent(absID) / 10000.,
 637                    fRecoUtils->GetEMCALChannelRecalibrationFactor(ism,icol,irow));
 638
 639         } // columns
 640       } // rows
 641     } // SM loop
 642   }else printf("AliAnalysisTaskEMCALPi0CalibSelection::SetOADBParameters() - Do NOT recalibrate EMCAL with T variations, no params TH1 \n");
 643
 644 }
 645
 646 //___________________________________________________________________
 647 void AliAnalysisTaskEMCALPi0CalibSelection::UserCreateOutputObjects()
 648 {
 649   //Create output container, init geometry
 650
 651   fEMCALGeo =  AliEMCALGeometry::GetInstance(fEMCALGeoName) ;
 652   Int_t nSM = (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules();
 653
 654   fOutputContainer = new TList();
 655   const Int_t buffersize = 255;
 656   char hname[buffersize], htitl[buffersize];
 657
 658   fhNEvents        = new TH1I("hNEvents", "Number of analyzed events"   , 1 , 0 , 1  ) ;
 659   fOutputContainer->Add(fhNEvents);
 660
 661   fHmgg = new TH2F("hmgg","2-cluster invariant mass",fNbins,fMinBin,fMaxBin,100,0,10);
 662   fHmgg->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
 663   fHmgg->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 664   fOutputContainer->Add(fHmgg);
 665
 666   fHmggDifferentSM = new TH2F("hmggDifferentSM","2-cluster invariant mass, different SM",fNbins,fMinBin,fMaxBin,100,0,10);
 667   fHmggDifferentSM->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
 668   fHmggDifferentSM->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 669   fOutputContainer->Add(fHmggDifferentSM);
 670
 671   fHOpeningAngle = new TH2F("hopang","2-cluster opening angle",100,0.,50.,100,0,10);
 672   fHOpeningAngle->SetXTitle("#alpha_{#gamma #gamma}");
 673   fHOpeningAngle->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 674   fOutputContainer->Add(fHOpeningAngle);
 675
 676   fHOpeningAngleDifferentSM = new TH2F("hopangDifferentSM","2-cluster opening angle, different SM",100,0,50.,100,0,10);
 677   fHOpeningAngleDifferentSM->SetXTitle("#alpha_{#gamma #gamma}");
 678   fHOpeningAngleDifferentSM->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 679   fOutputContainer->Add(fHOpeningAngleDifferentSM);
 680
 681   fHAsymmetry = new TH2F("hasym","2-cluster opening angle",100,0.,1.,100,0,10);
 682   fHAsymmetry->SetXTitle("a");
 683   fHAsymmetry->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 684   fOutputContainer->Add(fHAsymmetry);
 685
 686   fHAsymmetryDifferentSM = new TH2F("hasymDifferentSM","2-cluster opening angle, different SM",100,0,1.,100,0,10);
 687   fHAsymmetryDifferentSM->SetXTitle("a");
 688   fHAsymmetryDifferentSM->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 689   fOutputContainer->Add(fHAsymmetryDifferentSM);
 690
 691
 692   //TString pairname[] = {"A side (0-2)", "C side (1-3)","Row 0 (0-1)", "Row 1 (2-3)"};
 693
 694   fHmggMaskFrame = new TH2F("hmggMaskFrame","2-cluster invariant mass, frame masked",fNbins,fMinBin,fMaxBin,100,0,10);
 695   fHmggMaskFrame->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
 696   fHmggMaskFrame->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 697   fOutputContainer->Add(fHmggMaskFrame);
 698
 699   fHmggDifferentSMMaskFrame = new TH2F("hmggDifferentSMMaskFrame","2-cluster invariant mass, different SM, frame masked",
 700                                        fNbins,fMinBin,fMaxBin,100,0,10);
 701   fHmggDifferentSMMaskFrame->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
 702   fHmggDifferentSMMaskFrame->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 703   fOutputContainer->Add(fHmggDifferentSMMaskFrame);
 704
 705
 706   for(Int_t iSM = 0; iSM < nSM; iSM++)
 707   {
 708     snprintf(hname, buffersize, "hmgg_SM%d",iSM);
 709     snprintf(htitl, buffersize, "Two-gamma inv. mass for super mod %d",iSM);
 710     fHmggSM[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
 711     fHmggSM[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
 712     fHmggSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 713     fOutputContainer->Add(fHmggSM[iSM]);
 714
 715     snprintf(hname, buffersize, "hmgg_SM%d_MaskFrame",iSM);
 716     snprintf(htitl, buffersize, "Two-gamma inv. mass for super mod %d",iSM);
 717     fHmggSMMaskFrame[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
 718     fHmggSMMaskFrame[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
 719     fHmggSMMaskFrame[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 720     fOutputContainer->Add(fHmggSMMaskFrame[iSM]);
 721
 722
 723     if(iSM < nSM/2)
 724     {
 725       snprintf(hname,buffersize, "hmgg_PairSameSectorSM%d",iSM);
 726       snprintf(htitl,buffersize, "Two-gamma inv. mass for SM pair Sector: %d",iSM);
 727       fHmggPairSameSectorSM[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
 728       fHmggPairSameSectorSM[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
 729       fHmggPairSameSectorSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 730       fOutputContainer->Add(fHmggPairSameSectorSM[iSM]);
 731
 732       snprintf(hname,buffersize, "hmgg_PairSameSectorSM%d_MaskFrame",iSM);
 733       snprintf(htitl,buffersize, "Two-gamma inv. mass for SM pair Sector: %d",iSM);
 734       fHmggPairSameSectorSMMaskFrame[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
 735       fHmggPairSameSectorSMMaskFrame[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
 736       fHmggPairSameSectorSMMaskFrame[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 737       fOutputContainer->Add(fHmggPairSameSectorSMMaskFrame[iSM]);
 738
 739       fhClusterPairDiffTimeSameSector[iSM] = new TH2F(Form("hClusterPairDiffTimeSameSector%d",iSM),
 740                                                       Form("cluster pair time difference vs E, Sector %d",iSM),
 741                                                       100,0,10, 200,-100,100);
 742       fhClusterPairDiffTimeSameSector[iSM]->SetXTitle("E_{pair} (GeV)");
 743       fhClusterPairDiffTimeSameSector[iSM]->SetYTitle("#Delta t (ns)");
 744       fOutputContainer->Add(fhClusterPairDiffTimeSameSector[iSM]);
 745
 746
 747     }
 748
 749     if(iSM < nSM-2)
 750     {
 751       snprintf(hname,buffersize, "hmgg_PairSameSideSM%d",iSM);
 752       snprintf(htitl,buffersize, "Two-gamma inv. mass for SM pair Sector: %d",iSM);
 753       fHmggPairSameSideSM[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
 754       fHmggPairSameSideSM[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
 755       fHmggPairSameSideSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 756       fOutputContainer->Add(fHmggPairSameSideSM[iSM]);
 757
 758       snprintf(hname,buffersize, "hmgg_PairSameSideSM%d_MaskFrame",iSM);
 759       snprintf(htitl,buffersize, "Two-gamma inv. mass for SM pair Sector: %d",iSM);
 760       fHmggPairSameSideSMMaskFrame[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
 761       fHmggPairSameSideSMMaskFrame[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
 762       fHmggPairSameSideSMMaskFrame[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 763       fOutputContainer->Add(fHmggPairSameSideSMMaskFrame[iSM]);
 764
 765       fhClusterPairDiffTimeSameSide[iSM] = new TH2F(Form("hClusterPairDiffTimeSameSide%d",iSM),
 766                                                     Form("cluster pair time difference vs E,  Side %d",iSM),
 767                                                     100,0,10, 200,-100,100);
 768       fhClusterPairDiffTimeSameSide[iSM]->SetXTitle("E_{pair} (GeV)");
 769       fhClusterPairDiffTimeSameSide[iSM]->SetYTitle("#Delta t (ns)");
 770       fOutputContainer->Add(fhClusterPairDiffTimeSameSide[iSM]);
 771
 772     }
 773
 774     snprintf(hname, buffersize, "hopang_SM%d",iSM);
 775     snprintf(htitl, buffersize, "Opening angle for super mod %d",iSM);
 776     fHOpeningAngleSM[iSM] = new TH2F(hname,htitl,100,0.,50.,100,0,10);
 777     fHOpeningAngleSM[iSM]->SetXTitle("#alpha_{#gamma #gamma} (deg)");
 778     fHOpeningAngleSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 779     fOutputContainer->Add(fHOpeningAngleSM[iSM]);
 780
 781     snprintf(hname,buffersize, "hopang_PairSM%d",iSM);
 782     snprintf(htitl,buffersize, "Opening angle for SM pair: %d",iSM);
 783     fHOpeningAnglePairSM[iSM] = new TH2F(hname,htitl,100,0.,50.,100,0,10);
 784     fHOpeningAnglePairSM[iSM]->SetXTitle("#alpha_{#gamma #gamma} (deg)");
 785     fHOpeningAnglePairSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 786     fOutputContainer->Add(fHOpeningAnglePairSM[iSM]);
 787
 788     snprintf(hname, buffersize, "hasym_SM%d",iSM);
 789     snprintf(htitl, buffersize, "Asymmetry for super mod %d",iSM);
 790     fHAsymmetrySM[iSM] = new TH2F(hname,htitl,100,0.,1.,100,0,10);
 791     fHAsymmetrySM[iSM]->SetXTitle("a");
 792     fHAsymmetrySM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 793     fOutputContainer->Add(fHAsymmetrySM[iSM]);
 794
 795     snprintf(hname,buffersize, "hasym_PairSM%d",iSM);
 796     snprintf(htitl,buffersize, "Asymmetry for SM pair: %d",iSM);
 797     fHAsymmetryPairSM[iSM] = new TH2F(hname,htitl,100,0.,1.,100,0,10);
 798     fHAsymmetryPairSM[iSM]->SetXTitle("a");
 799     fHAsymmetryPairSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 800     fOutputContainer->Add(fHAsymmetryPairSM[iSM]);
 801
 802     Int_t colmax = 48;
 803     Int_t rowmax = 24;
 804
 805     fhTowerDecayPhotonHit[iSM]  = new TH2F (Form("hTowerDecPhotonHit_Mod%d",iSM),
 806                                             Form("Entries in grid of cells in Module %d",iSM),
 807                                             colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5);
 808     fhTowerDecayPhotonHit[iSM]->SetYTitle("row (phi direction)");
 809     fhTowerDecayPhotonHit[iSM]->SetXTitle("column (eta direction)");
 810     fOutputContainer->Add(fhTowerDecayPhotonHit[iSM]);
 811
 812     fhTowerDecayPhotonEnergy[iSM]  = new TH2F (Form("hTowerDecPhotonEnergy_Mod%d",iSM),
 813                                                Form("Accumulated energy in grid of cells in Module %d",iSM),
 814                                                colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5);
 815     fhTowerDecayPhotonEnergy[iSM]->SetYTitle("row (phi direction)");
 816     fhTowerDecayPhotonEnergy[iSM]->SetXTitle("column (eta direction)");
 817     fOutputContainer->Add(fhTowerDecayPhotonEnergy[iSM]);
 818
 819     fhTowerDecayPhotonAsymmetry[iSM]  = new TH2F (Form("hTowerDecPhotonAsymmetry_Mod%d",iSM),
 820                                                   Form("Accumulated asymmetry in grid of cells in Module %d",iSM),
 821                                                   colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5);
 822     fhTowerDecayPhotonAsymmetry[iSM]->SetYTitle("row (phi direction)");
 823     fhTowerDecayPhotonAsymmetry[iSM]->SetXTitle("column (eta direction)");
 824     fOutputContainer->Add(fhTowerDecayPhotonAsymmetry[iSM]);
 825
 826     fhTowerDecayPhotonHitMaskFrame[iSM]  = new TH2F (Form("hTowerDecPhotonHit_Mod%d_MaskFrame",iSM),Form("Entries in grid of cells in Module %d",iSM),
 827                                                      colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5);
 828     fhTowerDecayPhotonHitMaskFrame[iSM]->SetYTitle("row (phi direction)");
 829     fhTowerDecayPhotonHitMaskFrame[iSM]->SetXTitle("column (eta direction)");
 830     fOutputContainer->Add(fhTowerDecayPhotonHitMaskFrame[iSM]);
 831
 832     fhClusterTimeSM[iSM] = new TH2F(Form("hClusterTime_SM%d",iSM),"cluster time vs E",100,0,10, 100,0,1000);
 833     fhClusterTimeSM[iSM]->SetXTitle("E (GeV)");
 834     fhClusterTimeSM[iSM]->SetYTitle("t (ns)");
 835     fOutputContainer->Add(fhClusterTimeSM[iSM]);
 836
 837     fhClusterPairDiffTimeSameSM[iSM] = new TH2F(Form("hClusterPairDiffTimeSameSM%d",iSM),
 838                                                 Form("cluster pair time difference vs E, SM %d",iSM),
 839                                                 100,0,10, 200,-100,100);
 840     fhClusterPairDiffTimeSameSM[iSM]->SetXTitle("E (GeV)");
 841     fhClusterPairDiffTimeSameSM[iSM]->SetYTitle("#Delta t (ns)");
 842     fOutputContainer->Add(fhClusterPairDiffTimeSameSM[iSM]);
 843
 844   }
 845
 846   Int_t nchannels = nSM*AliEMCALGeoParams::fgkEMCALRows*AliEMCALGeoParams::fgkEMCALCols;
 847   for(Int_t ibc = 0; ibc < 4; ibc++)
 848   {
 849     fHTpi0[ibc] = new TH2F(Form("hTime_BC%d",ibc),Form("Time of cell clusters under pi0 peak, bunch crossing %d",ibc),
 850                            nchannels,0,nchannels, fNTimeBins,fMinTimeBin,fMaxTimeBin);
 851     fOutputContainer->Add(fHTpi0[ibc]);
 852     fHTpi0[ibc]->SetYTitle("time (ns)");
 853     fHTpi0[ibc]->SetXTitle("abs. Id. ");
 854   }
 855
 856
 857   fhClusterTime = new TH2F("hClusterTime","cluster time vs E",100,0,10, 100,0,1000);
 858   fhClusterTime->SetXTitle("E (GeV)");
 859   fhClusterTime->SetYTitle("t (ns)");
 860   fOutputContainer->Add(fhClusterTime);
 861
 862   fhClusterPairDiffTime = new TH2F("hClusterPairDiffTime","cluster pair time difference vs E",100,0,10, 800,-400,400);
 863   fhClusterPairDiffTime->SetXTitle("E_{pair} (GeV)");
 864   fhClusterPairDiffTime->SetYTitle("#Delta t (ns)");
 865   fOutputContainer->Add(fhClusterPairDiffTime);
 866
 867   for(Int_t iMod=0; iMod < nSM; iMod++)
 868   {
 869     for(Int_t iRow=0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++)
 870     {
 871       for(Int_t iCol=0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++)
 872       {
 873         snprintf(hname,buffersize, "%d_%d_%d",iMod,iCol,iRow);
 874         snprintf(htitl,buffersize, "Two-gamma inv. mass for super mod %d, cell(col,row)=(%d,%d)",iMod,iCol,iRow);
 875         fHmpi0[iMod][iCol][iRow] = new TH1F(hname,htitl,fNbins,fMinBin,fMaxBin);
 876         fHmpi0[iMod][iCol][iRow]->SetXTitle("mass (MeV/c^{2})");
 877         fOutputContainer->Add(fHmpi0[iMod][iCol][iRow]);
 878       }
 879     }
 880   }
 881
 882   fOutputContainer->SetOwner(kTRUE);
 883
 884   PostData(1,fOutputContainer);
 885
 886   // cuts container, set in terminate but init and post here
 887
 888   fCuts = new TList();
 889
 890   fCuts ->SetOwner(kTRUE);
 891
 892   PostData(2, fCuts);
 893
 894 }
 895
 896 //______________________________________________________________________________________________________
 897 Bool_t AliAnalysisTaskEMCALPi0CalibSelection::MaskFrameCluster(const Int_t iSM,  const Int_t ieta) const
 898 {
 899   //Check if cell is in one of the regions where we have significant amount of material in front of EMCAL
 900
 901   Int_t icol = ieta;
 902   if(iSM%2) icol+=48; // Impair SM, shift index [0-47] to [48-96]
 903
 904   if (fNMaskCellColumns && fMaskCellColumns)
 905   {
 906     for (Int_t imask = 0; imask < fNMaskCellColumns; imask++)
 907     {
 908       if(icol==fMaskCellColumns[imask]) return kTRUE;
 909     }
 910   }
 911
 912   return kFALSE;
 913
 914 }
 915
 916 //__________________________________________________________________________
 917 void AliAnalysisTaskEMCALPi0CalibSelection::UserExec(Option_t* /* option */)
 918 {
 919   // Do analysis, first select the events, then correct the clusters if needed
 920   // and finally fill the histograms per channel after recalibration
 921
 922   //Event selection
 923   if(fTriggerName!="")
 924   {
 925     AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (InputEvent());
 926     AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (InputEvent());
 927
 928     TString triggerClass = "";
 929     if     (esdevent) triggerClass = esdevent->GetFiredTriggerClasses();
 930     else if(aodevent) triggerClass = aodevent->GetFiredTriggerClasses();
 931
 932     if(triggerClass.Contains(fTriggerName))
 933     {
 934       //printf("Reject Event %d, FiredClass %s\n",(Int_t)Entry(),(((AliESDEvent*)InputEvent())->GetFiredTriggerClasses()).Data());
 935       return;
 936     }
 937   }
 938
 939   //Get the input event
 940   AliVEvent* event = 0;
 941   if(fFilteredInput) event = AODEvent();
 942   else               event = InputEvent();
 943
 944   if(!event)
 945   {
 946     printf("Input event not available!\n");
 947     return;
 948   }
 949
 950   if(DebugLevel() > 1)
 951     printf("AliAnalysisTaskEMCALPi0CalibSelection <<< %s: Event %d >>>\n",event->GetName(), (Int_t)Entry());
 952
 953   //Get the primary vertex
 954   event->GetPrimaryVertex()->GetXYZ(fVertex) ;
 955
 956   if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Vertex: (%.3f,%.3f,%.3f)\n",fVertex[0],fVertex[1],fVertex[2]);
 957
 958   //Int_t runNum = aod->GetRunNumber();
 959   //if(DebugLevel() > 1) printf("Run number: %d\n",runNum);
 960
 961   fhNEvents->Fill(0); //Count the events to be analyzed
 962
 963   // Acccess once the geometry matrix and temperature corrections
 964   if(fhNEvents->GetEntries()==1)
 965   {
 966     InitGeometryMatrices();
 967
 968     InitTemperatureCorrections();
 969   }
 970
 971   //Get the list of clusters and cells
 972   fEMCALCells       = event->GetEMCALCells();
 973
 974   fCaloClustersArr  = new TRefArray();
 975   event->GetEMCALClusters(fCaloClustersArr);
 976
 977   if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection - N CaloClusters: %d - N CaloCells %d \n",
 978                               fCaloClustersArr->GetEntriesFast(), fEMCALCells->GetNumberOfCells());
 979
 980   CorrectClusters(); // Non linearity, new calibration, T calibration
 981
 982   FillHistograms();
 983
 984   delete fCaloClustersArr;
 985
 986   PostData(1,fOutputContainer);
 987
 988 }
 989
 990 //_____________________________________________________
 991 void AliAnalysisTaskEMCALPi0CalibSelection::PrintInfo()
 992 {
 993   //Print settings
 994
 995   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",
 996          fEmin,fEmax, fMinNCells, fAsyCut, fDTimeCut,fTimeMin,fTimeMax) ;
 997
 998   printf("Group %d cells\n", fGroupNCells) ;
 999
1000   printf("Cluster maximal cell away from border at least %d cells\n", fRecoUtils->GetNumberOfCellsFromEMCALBorder()) ;
1001
1002   printf("Histograms: bins %d; energy range: %2.2f < E < %2.2f GeV\n",fNbins,fMinBin,fMaxBin) ;
1003
1004   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",
1005          fRecoUtils->IsBadChannelsRemovalSwitchedOn(),fFilteredInput,fCorrectClusters, fSameSM) ;
1006
1007   printf("EMCAL Geometry name: < %s >, Load Matrices %d\n",fEMCALGeoName.Data(), fLoadMatrices) ;
1008   if(fLoadMatrices) {for(Int_t ism = 0; ism < AliEMCALGeoParams::fgkEMCALModules; ism++) if(fMatrix[ism]) fMatrix[ism]->Print() ; }
1009
1010 }
1011
1012 //____________________________________________________________________
1013 void AliAnalysisTaskEMCALPi0CalibSelection::Terminate(Option_t*)
1014 {
1015   // Create cuts/param objects and publish to slot
1016   const Int_t buffersize = 255;
1017   char onePar[buffersize] ;
1018
1019   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",
1020            fEmin,fEmax, fL0min, fL0max, fMinNCells, fAsyCut, fDTimeCut, fTimeMin, fTimeMax, fInvMassCutMin, fInvMassCutMax) ;
1021   fCuts->Add(new TObjString(onePar));
1022   snprintf(onePar,buffersize, "Group %d cells;", fGroupNCells) ;
1023   fCuts->Add(new TObjString(onePar));
1024   snprintf(onePar,buffersize, "Cluster maximal cell away from border at least %d cells;", fRecoUtils->GetNumberOfCellsFromEMCALBorder()) ;
1025   fCuts->Add(new TObjString(onePar));
1026   snprintf(onePar,buffersize, "Histograms, Mass bins %d; energy range: %2.2f < E < %2.2f GeV;",fNbins,fMinBin,fMaxBin) ;
1027   fCuts->Add(new TObjString(onePar));
1028   snprintf(onePar,buffersize, "Histograms, Time bins %d; energy range: %2.2f < E < %2.2f GeV;",fNTimeBins,fMinTimeBin,fMaxTimeBin) ;
1029   fCuts->Add(new TObjString(onePar));
1030   snprintf(onePar,buffersize, "Switchs: Remove Bad Channels? %d; Use filtered input? %d;  Correct Clusters? %d, Mass per channel same SM clusters? %d ",
1031            fRecoUtils->IsBadChannelsRemovalSwitchedOn(),fFilteredInput,fCorrectClusters, fSameSM) ;
1032   fCuts->Add(new TObjString(onePar));
1033   snprintf(onePar,buffersize, "EMCAL Geometry name: < %s >, Load Matrices? %d",fEMCALGeoName.Data(),fLoadMatrices) ;
1034   fCuts->Add(new TObjString(onePar));
1035
1036   // Post Data
1037   PostData(2, fCuts);
1038
1039 }
1040