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