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