PWG4/CaloCalib/AliAnalysisTaskEMCALPi0CalibSelection.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: Boris Polishchuk                                               *
   5  * Adapted to AOD reading by Gustavo Conesa                               *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 //---------------------------------------------------------------------------//
  17 //                                                                           //
  18 // Fill histograms (one per cell) with two-cluster invariant mass            //
  19 // using calibration coefficients of the previous iteration.                 //
  20 // Histogram for a given cell is filled if the most energy of one cluster    //
  21 // is deposited in this cell and the other cluster could be anywherein EMCAL.//
  22 //                                                                           //
  23 //---------------------------------------------------------------------------//
  24
  25 //#include <cstdlib>
  26 //#include <Riostream.h>
  27 // Root
  28 #include "TLorentzVector.h"
  29 #include "TRefArray.h"
  30 #include "TList.h"
  31 #include "TH1F.h"
  32 #include <TGeoManager.h>
  33
  34 // AliRoot
  35 #include "AliAnalysisTaskEMCALPi0CalibSelection.h"
  36 #include "AliAODEvent.h"
  37 #include "AliESDEvent.h"
  38 #include "AliEMCALGeometry.h"
  39 #include "AliVCluster.h"
  40 #include "AliVCaloCells.h"
  41 #include "AliEMCALRecoUtils.h"
  42
  43 ClassImp(AliAnalysisTaskEMCALPi0CalibSelection)
  44
  45
  46 //__________________________________________________
  47 AliAnalysisTaskEMCALPi0CalibSelection::AliAnalysisTaskEMCALPi0CalibSelection(const char* name) :
  48   AliAnalysisTaskSE(name),fEMCALGeo(0x0),
  49   fEmin(0.5),               fEmax(15.),      fDTimeCut(20.),
  50   fAsyCut(1.),              fMinNCells(2),   fGroupNCells(0),
  51   fLogWeight(4.5),          fSameSM(kFALSE), fFilteredInput(kFALSE),
  52   fCorrectClusters(kFALSE), fEMCALGeoName("EMCAL_COMPLETEV1"),
  53   fTriggerName("EMC"),      fRecoUtils(new AliEMCALRecoUtils),
  54   fCuts(0x0),               fLoadMatrices(0),
  55   fNMaskCellColumns(11),    fMaskCellColumns(0x0),
  56   //Histograms
  57   fNbins(300), fMinBin(0.), fMaxBin(300.),   fOutputContainer(0x0),
  58   fHmgg(0x0),               fHmggDifferentSM(0x0),
  59   fHmggMaskFrame(0x0),      fHmggDifferentSMMaskFrame(0x0),
  60   fHOpeningAngle(0x0),      fHOpeningAngleDifferentSM(0x0),
  61   fHIncidentAngle(0x0),     fHIncidentAngleDifferentSM(0x0),
  62   fHAsymmetry(0x0),         fHAsymmetryDifferentSM(0x0),
  63   fhNEvents(0x0),
  64   fhClusterTime(0x0),       fhClusterPairDiffTime(0x0)
  65 {
  66   //Named constructor which should be used.
  67
  68   for(Int_t iMod=0; iMod < AliEMCALGeoParams::fgkEMCALModules; iMod++) {
  69     for(Int_t iX=0; iX<24; iX++) {
  70       for(Int_t iZ=0; iZ<48; iZ++) {
  71         fHmpi0[iMod][iZ][iX]=0;
  72       }
  73     }
  74   }
  75
  76   fMaskCellColumns = new Int_t[fNMaskCellColumns];
  77   fMaskCellColumns[0] = 6 ;  fMaskCellColumns[1] = 7 ;  fMaskCellColumns[2] = 8 ;
  78   fMaskCellColumns[3] = 35;  fMaskCellColumns[4] = 36;  fMaskCellColumns[5] = 37;
  79   fMaskCellColumns[6] = 12+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[7] = 13+AliEMCALGeoParams::fgkEMCALCols;
  80   fMaskCellColumns[8] = 40+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[9] = 41+AliEMCALGeoParams::fgkEMCALCols;
  81   fMaskCellColumns[10]= 42+AliEMCALGeoParams::fgkEMCALCols;
  82
  83   for(Int_t iSMPair = 0; iSMPair < AliEMCALGeoParams::fgkEMCALModules/2; iSMPair++) {
  84     fHmggPairSameSectorSM[iSMPair]          = 0;
  85     fHmggPairSameSectorSMMaskFrame[iSMPair] = 0;
  86     fhClusterPairDiffTimeSameSector[iSMPair]= 0;
  87   }
  88   for(Int_t iSMPair = 0; iSMPair < AliEMCALGeoParams::fgkEMCALModules-2; iSMPair++){
  89     fHmggPairSameSideSM[iSMPair]            = 0;
  90     fHmggPairSameSideSMMaskFrame[iSMPair]   = 0;
  91     fhClusterPairDiffTimeSameSide[iSMPair]  = 0;
  92   }
  93
  94   for(Int_t iSM = 0; iSM < AliEMCALGeoParams::fgkEMCALModules; iSM++) {
  95     fHmggSM[iSM]                     = 0;
  96     fHmggSMMaskFrame[iSM]            = 0;
  97     fHOpeningAngleSM[iSM]            = 0;
  98     fHOpeningAnglePairSM[iSM]        = 0;
  99     fHAsymmetrySM[iSM]               = 0;
 100     fHAsymmetryPairSM[iSM]           = 0;
 101     fHIncidentAngleSM[iSM]           = 0;
 102     fHIncidentAnglePairSM[iSM]       = 0;
 103     fhTowerDecayPhotonHit[iSM]       = 0;
 104     fhTowerDecayPhotonEnergy[iSM]    = 0;
 105     fhTowerDecayPhotonAsymmetry[iSM] = 0;
 106     fhTowerDecayPhotonHitMaskFrame[iSM]= 0;
 107     fMatrix[iSM]                     = 0x0;
 108     fhClusterTimeSM[iSM]             = 0;
 109     fhClusterPairDiffTimeSameSM[iSM] = 0;
 110   }
 111
 112   DefineOutput(1, TList::Class());
 113   DefineOutput(2, TList::Class());  // will contain cuts or local params
 114
 115 }
 116
 117 //__________________________________________________
 118 AliAnalysisTaskEMCALPi0CalibSelection::~AliAnalysisTaskEMCALPi0CalibSelection()
 119 {
 120   //Destructor.
 121
 122   if(fOutputContainer){
 123     fOutputContainer->Delete() ;
 124     delete fOutputContainer ;
 125   }
 126
 127   if(fEMCALGeo)         delete fEMCALGeo  ;
 128   if(fRecoUtils)        delete fRecoUtils ;
 129   if(fNMaskCellColumns) delete [] fMaskCellColumns;
 130
 131 }
 132
 133 //_____________________________________________________
 134 void AliAnalysisTaskEMCALPi0CalibSelection::LocalInit()
 135 {
 136   // Local Initialization
 137
 138   // Create cuts/param objects and publish to slot
 139   const Int_t buffersize = 255;
 140   char onePar[buffersize] ;
 141   fCuts = new TList();
 142
 143   snprintf(onePar,buffersize, "Custer cuts: %2.2f < E < %2.2f GeV; min number of cells %d; Assymetry cut %1.2f", fEmin,fEmax, fMinNCells, fAsyCut) ;
 144   fCuts->Add(new TObjString(onePar));
 145   snprintf(onePar,buffersize, "Group %d cells;", fGroupNCells) ;
 146   fCuts->Add(new TObjString(onePar));
 147   snprintf(onePar,buffersize, "Cluster maximal cell away from border at least %d cells;", fRecoUtils->GetNumberOfCellsFromEMCALBorder()) ;
 148   fCuts->Add(new TObjString(onePar));
 149   snprintf(onePar,buffersize, "Histograms: bins %d; energy range: %2.2f < E < %2.2f GeV;",fNbins,fMinBin,fMaxBin) ;
 150   fCuts->Add(new TObjString(onePar));
 151   snprintf(onePar,buffersize, "Switchs: Remove Bad Channels? %d; Use filtered input? %d;  Correct Clusters? %d, Mass per channel same SM clusters? %d ",
 152            fRecoUtils->IsBadChannelsRemovalSwitchedOn(),fFilteredInput,fCorrectClusters, fSameSM) ;
 153   fCuts->Add(new TObjString(onePar));
 154   snprintf(onePar,buffersize, "EMCAL Geometry name: < %s >, Load Matrices? %d",fEMCALGeoName.Data(),fLoadMatrices) ;
 155   fCuts->Add(new TObjString(onePar));
 156
 157   fCuts ->SetOwner(kTRUE);
 158
 159   // Post Data
 160   PostData(2, fCuts);
 161
 162 }
 163
 164 //__________________________________________________
 165 void AliAnalysisTaskEMCALPi0CalibSelection::UserCreateOutputObjects()
 166 {
 167   //Create output container, init geometry
 168
 169   fEMCALGeo =  AliEMCALGeometry::GetInstance(fEMCALGeoName) ;
 170   Int_t nSM = (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules();
 171
 172   fOutputContainer = new TList();
 173   const Int_t buffersize = 255;
 174   char hname[buffersize], htitl[buffersize];
 175
 176   for(Int_t iMod=0; iMod < nSM; iMod++) {
 177     for(Int_t iRow=0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
 178       for(Int_t iCol=0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
 179         snprintf(hname,buffersize, "%d_%d_%d",iMod,iCol,iRow);
 180         snprintf(htitl,buffersize, "Two-gamma inv. mass for super mod %d, cell(col,row)=(%d,%d)",iMod,iCol,iRow);
 181         fHmpi0[iMod][iCol][iRow] = new TH1F(hname,htitl,fNbins,fMinBin,fMaxBin);
 182         fOutputContainer->Add(fHmpi0[iMod][iCol][iRow]);
 183       }
 184     }
 185   }
 186
 187   fHmgg = new TH2F("hmgg","2-cluster invariant mass",fNbins,fMinBin,fMaxBin,100,0,10);
 188   fHmgg->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
 189   fHmgg->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 190   fOutputContainer->Add(fHmgg);
 191
 192   fHmggDifferentSM = new TH2F("hmggDifferentSM","2-cluster invariant mass, different SM",fNbins,fMinBin,fMaxBin,100,0,10);
 193   fHmggDifferentSM->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
 194   fHmggDifferentSM->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 195   fOutputContainer->Add(fHmggDifferentSM);
 196
 197   fHOpeningAngle = new TH2F("hopang","2-cluster opening angle",100,0.,50.,100,0,10);
 198   fHOpeningAngle->SetXTitle("#alpha_{#gamma #gamma}");
 199   fHOpeningAngle->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 200   fOutputContainer->Add(fHOpeningAngle);
 201
 202   fHOpeningAngleDifferentSM = new TH2F("hopangDifferentSM","2-cluster opening angle, different SM",100,0,50.,100,0,10);
 203   fHOpeningAngleDifferentSM->SetXTitle("#alpha_{#gamma #gamma}");
 204   fHOpeningAngleDifferentSM->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 205   fOutputContainer->Add(fHOpeningAngleDifferentSM);
 206
 207   fHIncidentAngle = new TH2F("hinang","#gamma incident angle in SM",100,0.,20.,100,0,10);
 208   fHIncidentAngle->SetXTitle("#alpha_{#gamma SM center}");
 209   fHIncidentAngle->SetYTitle("p_{T #gamma} (GeV/c)");
 210   fOutputContainer->Add(fHIncidentAngle);
 211
 212   fHIncidentAngleDifferentSM = new TH2F("hinangDifferentSM","#gamma incident angle in SM, different SM pair",100,0,20.,100,0,10);
 213   fHIncidentAngleDifferentSM->SetXTitle("#alpha_{#gamma - SM center}");
 214   fHIncidentAngleDifferentSM->SetYTitle("p_{T #gamma} (GeV/c)");
 215   fOutputContainer->Add(fHIncidentAngleDifferentSM);
 216
 217   fHAsymmetry = new TH2F("hasym","2-cluster opening angle",100,0.,1.,100,0,10);
 218   fHAsymmetry->SetXTitle("a");
 219   fHAsymmetry->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 220   fOutputContainer->Add(fHAsymmetry);
 221
 222   fHAsymmetryDifferentSM = new TH2F("hasymDifferentSM","2-cluster opening angle, different SM",100,0,1.,100,0,10);
 223   fHAsymmetryDifferentSM->SetXTitle("a");
 224   fHAsymmetryDifferentSM->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 225   fOutputContainer->Add(fHAsymmetryDifferentSM);
 226
 227
 228   //TString pairname[] = {"A side (0-2)", "C side (1-3)","Row 0 (0-1)", "Row 1 (2-3)"};
 229
 230   fHmggMaskFrame = new TH2F("hmggMaskFrame","2-cluster invariant mass, frame masked",fNbins,fMinBin,fMaxBin,100,0,10);
 231   fHmggMaskFrame->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
 232   fHmggMaskFrame->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 233   fOutputContainer->Add(fHmggMaskFrame);
 234
 235   fHmggDifferentSMMaskFrame = new TH2F("hmggDifferentSMMaskFrame","2-cluster invariant mass, different SM, frame masked",
 236                                        fNbins,fMinBin,fMaxBin,100,0,10);
 237   fHmggDifferentSMMaskFrame->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
 238   fHmggDifferentSMMaskFrame->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 239   fOutputContainer->Add(fHmggDifferentSMMaskFrame);
 240
 241
 242   for(Int_t iSM = 0; iSM < nSM; iSM++) {
 243
 244     snprintf(hname, buffersize, "hmgg_SM%d",iSM);
 245     snprintf(htitl, buffersize, "Two-gamma inv. mass for super mod %d",iSM);
 246     fHmggSM[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
 247     fHmggSM[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
 248     fHmggSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 249     fOutputContainer->Add(fHmggSM[iSM]);
 250
 251     snprintf(hname, buffersize, "hmgg_SM%d_MaskFrame",iSM);
 252     snprintf(htitl, buffersize, "Two-gamma inv. mass for super mod %d",iSM);
 253     fHmggSMMaskFrame[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
 254     fHmggSMMaskFrame[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
 255     fHmggSMMaskFrame[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 256     fOutputContainer->Add(fHmggSMMaskFrame[iSM]);
 257
 258
 259     if(iSM < nSM/2){
 260       snprintf(hname,buffersize, "hmgg_PairSameSectorSM%d",iSM);
 261       snprintf(htitl,buffersize, "Two-gamma inv. mass for SM pair Sector: %d",iSM);
 262       fHmggPairSameSectorSM[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
 263       fHmggPairSameSectorSM[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
 264       fHmggPairSameSectorSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 265       fOutputContainer->Add(fHmggPairSameSectorSM[iSM]);
 266
 267       snprintf(hname,buffersize, "hmgg_PairSameSectorSM%d_MaskFrame",iSM);
 268       snprintf(htitl,buffersize, "Two-gamma inv. mass for SM pair Sector: %d",iSM);
 269       fHmggPairSameSectorSMMaskFrame[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
 270       fHmggPairSameSectorSMMaskFrame[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
 271       fHmggPairSameSectorSMMaskFrame[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 272       fOutputContainer->Add(fHmggPairSameSectorSMMaskFrame[iSM]);
 273
 274       fhClusterPairDiffTimeSameSector[iSM] = new TH2F(Form("hClusterPairDiffTimeSameSector%d",iSM),
 275                                                       Form("cluster pair time difference vs E, Sector %d",iSM),
 276                                                       100,0,10, 200,-100,100);
 277       fhClusterPairDiffTimeSameSector[iSM]->SetXTitle("E_{pair} (GeV)");
 278       fhClusterPairDiffTimeSameSector[iSM]->SetYTitle("#Delta t (ns)");
 279       fOutputContainer->Add(fhClusterPairDiffTimeSameSector[iSM]);
 280
 281
 282     }
 283
 284     if(iSM < nSM-2){
 285       snprintf(hname,buffersize, "hmgg_PairSameSideSM%d",iSM);
 286       snprintf(htitl,buffersize, "Two-gamma inv. mass for SM pair Sector: %d",iSM);
 287       fHmggPairSameSideSM[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
 288       fHmggPairSameSideSM[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
 289       fHmggPairSameSideSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 290       fOutputContainer->Add(fHmggPairSameSideSM[iSM]);
 291
 292       snprintf(hname,buffersize, "hmgg_PairSameSideSM%d_MaskFrame",iSM);
 293       snprintf(htitl,buffersize, "Two-gamma inv. mass for SM pair Sector: %d",iSM);
 294       fHmggPairSameSideSMMaskFrame[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
 295       fHmggPairSameSideSMMaskFrame[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
 296       fHmggPairSameSideSMMaskFrame[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 297       fOutputContainer->Add(fHmggPairSameSideSMMaskFrame[iSM]);
 298
 299       fhClusterPairDiffTimeSameSide[iSM] = new TH2F(Form("hClusterPairDiffTimeSameSide%d",iSM),
 300                                                     Form("cluster pair time difference vs E,  Side %d",iSM),
 301                                                     100,0,10, 200,-100,100);
 302       fhClusterPairDiffTimeSameSide[iSM]->SetXTitle("E_{pair} (GeV)");
 303       fhClusterPairDiffTimeSameSide[iSM]->SetYTitle("#Delta t (ns)");
 304       fOutputContainer->Add(fhClusterPairDiffTimeSameSide[iSM]);
 305
 306     }
 307
 308     snprintf(hname, buffersize, "hopang_SM%d",iSM);
 309     snprintf(htitl, buffersize, "Opening angle for super mod %d",iSM);
 310     fHOpeningAngleSM[iSM] = new TH2F(hname,htitl,100,0.,50.,100,0,10);
 311     fHOpeningAngleSM[iSM]->SetXTitle("#alpha_{#gamma #gamma} (deg)");
 312     fHOpeningAngleSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 313     fOutputContainer->Add(fHOpeningAngleSM[iSM]);
 314
 315     snprintf(hname,buffersize, "hopang_PairSM%d",iSM);
 316     snprintf(htitl,buffersize, "Opening angle for SM pair: %d",iSM);
 317     fHOpeningAnglePairSM[iSM] = new TH2F(hname,htitl,100,0.,50.,100,0,10);
 318     fHOpeningAnglePairSM[iSM]->SetXTitle("#alpha_{#gamma #gamma} (deg)");
 319     fHOpeningAnglePairSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 320     fOutputContainer->Add(fHOpeningAnglePairSM[iSM]);
 321
 322     snprintf(hname, buffersize, "hinang_SM%d",iSM);
 323     snprintf(htitl, buffersize, "Incident angle for super mod %d",iSM);
 324     fHIncidentAngleSM[iSM] = new TH2F(hname,htitl,100,0.,20.,100,0,10);
 325     fHIncidentAngleSM[iSM]->SetXTitle("#alpha_{#gamma - SM center} (deg)");
 326     fHIncidentAngleSM[iSM]->SetYTitle("p_{T #gamma} (GeV/c)");
 327     fOutputContainer->Add(fHIncidentAngleSM[iSM]);
 328
 329     snprintf(hname,buffersize, "hinang_PairSM%d",iSM);
 330     snprintf(htitl,buffersize, "Incident angle for SM pair: %d",iSM);
 331     fHIncidentAnglePairSM[iSM] = new TH2F(hname,htitl,100,0.,20.,100,0,10);
 332     fHIncidentAnglePairSM[iSM]->SetXTitle("#alpha_{#gamma - SM center} (deg)");
 333     fHIncidentAnglePairSM[iSM]->SetYTitle("p_{T #gamma} (GeV/c)");
 334     fOutputContainer->Add(fHIncidentAnglePairSM[iSM]);
 335
 336     snprintf(hname, buffersize, "hasym_SM%d",iSM);
 337     snprintf(htitl, buffersize, "Asymmetry for super mod %d",iSM);
 338     fHAsymmetrySM[iSM] = new TH2F(hname,htitl,100,0.,1.,100,0,10);
 339     fHAsymmetrySM[iSM]->SetXTitle("a");
 340     fHAsymmetrySM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 341     fOutputContainer->Add(fHAsymmetrySM[iSM]);
 342
 343     snprintf(hname,buffersize, "hasym_PairSM%d",iSM);
 344     snprintf(htitl,buffersize, "Asymmetry for SM pair: %d",iSM);
 345     fHAsymmetryPairSM[iSM] = new TH2F(hname,htitl,100,0.,1.,100,0,10);
 346     fHAsymmetryPairSM[iSM]->SetXTitle("a");
 347     fHAsymmetryPairSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
 348     fOutputContainer->Add(fHAsymmetryPairSM[iSM]);
 349
 350     Int_t colmax = 48;
 351     Int_t rowmax = 24;
 352
 353     fhTowerDecayPhotonHit[iSM]  = new TH2F (Form("hTowerDecPhotonHit_Mod%d",iSM),
 354                                             Form("Entries in grid of cells in Module %d",iSM),
 355                                             colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5);
 356     fhTowerDecayPhotonHit[iSM]->SetYTitle("row (phi direction)");
 357     fhTowerDecayPhotonHit[iSM]->SetXTitle("column (eta direction)");
 358     fOutputContainer->Add(fhTowerDecayPhotonHit[iSM]);
 359
 360     fhTowerDecayPhotonEnergy[iSM]  = new TH2F (Form("hTowerDecPhotonEnergy_Mod%d",iSM),
 361                                                Form("Accumulated energy in grid of cells in Module %d",iSM),
 362                                                colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5);
 363     fhTowerDecayPhotonEnergy[iSM]->SetYTitle("row (phi direction)");
 364     fhTowerDecayPhotonEnergy[iSM]->SetXTitle("column (eta direction)");
 365     fOutputContainer->Add(fhTowerDecayPhotonEnergy[iSM]);
 366
 367     fhTowerDecayPhotonAsymmetry[iSM]  = new TH2F (Form("hTowerDecPhotonAsymmetry_Mod%d",iSM),
 368                                                   Form("Accumulated asymmetry in grid of cells in Module %d",iSM),
 369                                                   colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5);
 370     fhTowerDecayPhotonAsymmetry[iSM]->SetYTitle("row (phi direction)");
 371     fhTowerDecayPhotonAsymmetry[iSM]->SetXTitle("column (eta direction)");
 372     fOutputContainer->Add(fhTowerDecayPhotonAsymmetry[iSM]);
 373
 374     fhTowerDecayPhotonHitMaskFrame[iSM]  = new TH2F (Form("hTowerDecPhotonHit_Mod%d_MaskFrame",iSM),Form("Entries in grid of cells in Module %d",iSM),
 375                                             colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5);
 376     fhTowerDecayPhotonHitMaskFrame[iSM]->SetYTitle("row (phi direction)");
 377     fhTowerDecayPhotonHitMaskFrame[iSM]->SetXTitle("column (eta direction)");
 378     fOutputContainer->Add(fhTowerDecayPhotonHitMaskFrame[iSM]);
 379
 380     fhClusterTimeSM[iSM] = new TH2F(Form("hClusterTime_SM%d",iSM),"cluster time vs E",100,0,10, 100,0,1000);
 381     fhClusterTimeSM[iSM]->SetXTitle("E (GeV)");
 382     fhClusterTimeSM[iSM]->SetYTitle("t (ns)");
 383     fOutputContainer->Add(fhClusterTimeSM[iSM]);
 384
 385     fhClusterPairDiffTimeSameSM[iSM] = new TH2F(Form("hClusterPairDiffTimeSameSM%d",iSM),
 386                                                 Form("cluster pair time difference vs E, SM %d",iSM),
 387                                                 100,0,10, 200,-100,100);
 388     fhClusterPairDiffTimeSameSM[iSM]->SetXTitle("E (GeV)");
 389     fhClusterPairDiffTimeSameSM[iSM]->SetYTitle("#Delta t (ns)");
 390     fOutputContainer->Add(fhClusterPairDiffTimeSameSM[iSM]);
 391
 392   }
 393
 394   fhClusterTime = new TH2F("hClusterTime","cluster time vs E",100,0,10, 100,0,1000);
 395   fhClusterTime->SetXTitle("E (GeV)");
 396   fhClusterTime->SetYTitle("t (ns)");
 397   fOutputContainer->Add(fhClusterTime);
 398
 399   fhClusterPairDiffTime = new TH2F("hClusterPairDiffTime","cluster pair time difference vs E",100,0,10, 200,-100,100);
 400   fhClusterPairDiffTime->SetXTitle("E_{pair} (GeV)");
 401   fhClusterPairDiffTime->SetYTitle("#Delta t (ns)");
 402   fOutputContainer->Add(fhClusterPairDiffTime);
 403
 404
 405   fhNEvents        = new TH1I("hNEvents", "Number of analyzed events"   , 1 , 0 , 1  ) ;
 406   fOutputContainer->Add(fhNEvents);
 407
 408   fOutputContainer->SetOwner(kTRUE);
 409
 410   //  fCalibData = new AliEMCALCalibData();
 411
 412   PostData(1,fOutputContainer);
 413
 414 }
 415
 416 //__________________________________________________
 417 Bool_t AliAnalysisTaskEMCALPi0CalibSelection::MaskFrameCluster(const Int_t iSM,  const Int_t ieta) const {
 418   //Check if cell is in one of the regions where we have significant amount of material in front of EMCAL
 419
 420   Int_t icol = ieta;
 421   if(iSM%2) icol+=48; // Impair SM, shift index [0-47] to [48-96]
 422
 423   if (fNMaskCellColumns && fMaskCellColumns) {
 424     for (Int_t imask = 0; imask < fNMaskCellColumns; imask++) {
 425       if(icol==fMaskCellColumns[imask]) return kTRUE;
 426     }
 427   }
 428
 429   return kFALSE;
 430
 431 }
 432
 433 //__________________________________________________
 434 void AliAnalysisTaskEMCALPi0CalibSelection::UserExec(Option_t* /* option */)
 435 {
 436   //Analysis per event.
 437
 438   if(fRecoUtils->GetParticleType()!=AliEMCALRecoUtils::kPhoton){
 439     printf("Wrong particle type for cluster position recalculation! = %d\n", fRecoUtils->GetParticleType());
 440     abort();
 441   }
 442
 443   if(!(((AliESDEvent*)InputEvent())->GetFiredTriggerClasses()).Contains(fTriggerName)) {
 444     //printf("Reject Event %d, FiredClass %s\n",(Int_t)Entry(),(((AliESDEvent*)InputEvent())->GetFiredTriggerClasses()).Data());
 445     return;
 446   }
 447
 448   fhNEvents->Fill(0); //Event analyzed
 449
 450   //Get the input event
 451   AliVEvent* event = 0;
 452   if(fFilteredInput) event = AODEvent();
 453   else               event = InputEvent();
 454
 455   if(!event) {
 456     printf("Input event not available!\n");
 457     return;
 458   }
 459
 460   if(DebugLevel() > 1)
 461     printf("AliAnalysisTaskEMCALPi0CalibSelection <<< %s: Event %d >>>\n",event->GetName(), (Int_t)Entry());
 462
 463
 464   //Get the primary vertex
 465   Double_t v[3];
 466   event->GetPrimaryVertex()->GetXYZ(v) ;
 467
 468   if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Vertex: (%.3f,%.3f,%.3f)\n",v[0],v[1],v[2]);
 469
 470   //Int_t runNum = aod->GetRunNumber();
 471   //if(DebugLevel() > 1) printf("Run number: %d\n",runNum);
 472
 473   Int_t nSM = (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules();
 474   //Get the matrix with geometry information
 475   if(fhNEvents->GetEntries()==1){
 476     if(fLoadMatrices){
 477       printf("AliAnalysisTaskEMCALPi0CalibSelection::UserExec() - Load user defined geometry matrices\n");
 478       for(Int_t mod=0; mod < nSM ; mod++){
 479         if(fMatrix[mod]){
 480           if(DebugLevel() > 1)
 481             fMatrix[mod]->Print();
 482           fEMCALGeo->SetMisalMatrix(fMatrix[mod],mod) ;
 483         }
 484       }//SM loop
 485     }//Load matrices
 486     else if(!gGeoManager){
 487       printf("AliAnalysisTaskEMCALPi0CalibSelection::UserExec() - Get geo matrices from data\n");
 488       //Still not implemented in AOD, just a workaround to be able to work at least with ESDs
 489       if(!strcmp(event->GetName(),"AliAODEvent")) {
 490         if(DebugLevel() > 1)
 491           printf("AliAnalysisTaskEMCALPi0CalibSelection Use ideal geometry, values geometry matrix not kept in AODs.\n");
 492       }//AOD
 493       else {
 494         if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Load Misaligned matrices. \n");
 495         AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event) ;
 496         if(!esd) {
 497           printf("AliAnalysisTaskEMCALPi0CalibSelection::UserExec() - This event does not contain ESDs?");
 498           return;
 499         }
 500         for(Int_t mod=0; mod < nSM; mod++){
 501           if(DebugLevel() > 1)
 502             esd->GetEMCALMatrix(mod)->Print();
 503           if(esd->GetEMCALMatrix(mod)) fEMCALGeo->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
 504         }
 505       }//ESD
 506     }//Load matrices from Data
 507   }//first event
 508
 509   if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Will use fLogWeight %.3f .\n",fLogWeight);
 510   Int_t absId1   = -1;
 511   Int_t iSupMod1 = -1;
 512   Int_t iphi1    = -1;
 513   Int_t ieta1    = -1;
 514   Int_t absId2   = -1;
 515   Int_t iSupMod2 = -1;
 516   Int_t iphi2    = -1;
 517   Int_t ieta2    = -1;
 518   Bool_t shared  = kFALSE;
 519
 520   TLorentzVector p1;
 521   TLorentzVector p2;
 522   TLorentzVector p12;
 523
 524   //Get the list of clusters
 525   TRefArray * caloClustersArr  = new TRefArray();
 526   event->GetEMCALClusters(caloClustersArr);
 527   const Int_t kNumberOfEMCALClusters   = caloClustersArr->GetEntries() ;
 528   if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection - N CaloClusters: %d \n", kNumberOfEMCALClusters);
 529
 530   // Get EMCAL cells
 531   AliVCaloCells *emCells = event->GetEMCALCells();
 532
 533   // loop over EMCAL clusters
 534   //----------------------------------------------------------
 535   // First recalibrate and recalculate energy and position
 536   Float_t pos[]={0,0,0};
 537   if(fCorrectClusters){
 538     for(Int_t iClu=0; iClu<kNumberOfEMCALClusters; iClu++) {
 539       AliVCluster *c1 = (AliVCluster *) caloClustersArr->At(iClu);
 540
 541       Float_t e1i = c1->E();   // cluster energy before correction
 542       if      (e1i < fEmin) continue;
 543       else if (e1i > fEmax) continue;
 544       else if (c1->GetNCells() < fMinNCells) continue;
 545
 546       if(fRecoUtils->ClusterContainsBadChannel(fEMCALGeo, c1->GetCellsAbsId(), c1->GetNCells())) continue;
 547
 548       if(DebugLevel() > 2)
 549       {
 550         printf("Std  : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",c1->GetID(),c1->E(),c1->GetDispersion(),c1->GetM02(),c1->GetM20());
 551         c1->GetPosition(pos);
 552         printf("Std  : i %d, x %f, y %f, z %f\n",c1->GetID(), pos[0], pos[1], pos[2]);
 553       }
 554
 555       //Correct cluster energy and position if requested, and not corrected previously, by default Off
 556       if(fRecoUtils->IsRecalibrationOn())       {
 557         fRecoUtils->RecalibrateClusterEnergy(fEMCALGeo, c1, emCells);
 558         fRecoUtils->RecalculateClusterShowerShapeParameters(fEMCALGeo, emCells,c1);
 559         fRecoUtils->RecalculateClusterPID(c1);
 560       }
 561       if(DebugLevel() > 2)
 562         printf("Energy: after recalibration %f; ",c1->E());
 563
 564       // Recalculate cluster position
 565       fRecoUtils->RecalculateClusterPosition(fEMCALGeo, emCells,c1);
 566
 567       // Correct Non-Linearity
 568
 569       c1->SetE(fRecoUtils->CorrectClusterEnergyLinearity(c1));
 570       if(DebugLevel() > 2)
 571         printf("after linearity correction %f\n",c1->E());
 572
 573       if(DebugLevel() > 2)
 574       {
 575         printf("Cor  : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",c1->GetID(),c1->E(),c1->GetDispersion(),c1->GetM02(),c1->GetM20());
 576         c1->GetPosition(pos);
 577         printf("Cor  : i %d, x %f, y %f, z %f\n",c1->GetID(), pos[0], pos[1], pos[2]);
 578       }
 579     }
 580   }
 581
 582   //----------------------------------------------------------
 583   //Now the invariant mass analysis with the corrected clusters
 584   for(Int_t iClu=0; iClu<kNumberOfEMCALClusters-1; iClu++) {
 585
 586     AliVCluster *c1 = (AliVCluster *) caloClustersArr->At(iClu);
 587     if(fRecoUtils->ClusterContainsBadChannel(fEMCALGeo, c1->GetCellsAbsId(), c1->GetNCells())) continue;
 588
 589     Float_t e1i = c1->E();   // cluster energy before correction
 590     if      (e1i < fEmin) continue;
 591     else if (e1i > fEmax) continue;
 592     else if (c1->GetNCells() < fMinNCells) continue;
 593
 594     if(DebugLevel() > 2)
 595     {
 596       printf("IMA  : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",c1->GetID(),e1i,c1->GetDispersion(),c1->GetM02(),c1->GetM20());
 597       c1->GetPosition(pos);
 598       printf("IMA  : i %d, x %f, y %f, z %f\n",c1->GetID(), pos[0], pos[1], pos[2]);
 599     }
 600
 601     fRecoUtils->GetMaxEnergyCell(fEMCALGeo, emCells,c1,absId1,iSupMod1,ieta1,iphi1,shared);
 602     c1->GetMomentum(p1,v);
 603
 604     //Check if cluster is in fidutial region, not too close to borders
 605     Bool_t in1 = fRecoUtils->CheckCellFiducialRegion(fEMCALGeo, c1, emCells);
 606     // Clusters not facing frame structures
 607     Bool_t mask1 = MaskFrameCluster(iSupMod1, ieta1);
 608     //if(mask1) printf("Reject eta %d SM %d\n",ieta1, iSupMod1);
 609
 610     Double_t time1 = c1->GetTOF()*1.e9;
 611     fhClusterTime            ->Fill(c1->E(),time1);
 612     fhClusterTimeSM[iSupMod1]->Fill(c1->E(),time1);
 613
 614     // Combine cluster with other clusters and get the invariant mass
 615     for (Int_t jClu=iClu+1; jClu<kNumberOfEMCALClusters; jClu++) {
 616       AliAODCaloCluster *c2 = (AliAODCaloCluster *) caloClustersArr->At(jClu);
 617
 618       if(fRecoUtils->ClusterContainsBadChannel(fEMCALGeo, c2->GetCellsAbsId(), c2->GetNCells())) continue;
 619
 620       Float_t e2i = c2->E();
 621       if      (e2i < fEmin) continue;
 622       else if (e2i > fEmax) continue;
 623       else if (c2->GetNCells() < fMinNCells) continue;
 624
 625
 626       fRecoUtils->GetMaxEnergyCell(fEMCALGeo, emCells,c2,absId2,iSupMod2,ieta2,iphi2,shared);
 627       c2->GetMomentum(p2,v);
 628
 629       p12 = p1+p2;
 630       Float_t invmass = p12.M()*1000;
 631       //printf("*** mass %f\n",invmass);
 632
 633       //Asimetry cut
 634       Float_t asym = TMath::Abs(p1.E()-p2.E())/(p1.E()+p2.E());
 635       //printf("asymmetry %f\n",asym);
 636       if(asym > fAsyCut) continue;
 637
 638       //Time cut
 639       Double_t time2 = c2->GetTOF()*1.e9;
 640       fhClusterPairDiffTime->Fill(p12.E(),time1-time2);
 641       if(TMath::Abs(time1-time2) > fDTimeCut) continue;
 642
 643       if(invmass < fMaxBin && invmass > fMinBin ){
 644
 645         //Check if cluster is in fidutial region, not too close to borders
 646         Bool_t in2 = fRecoUtils->CheckCellFiducialRegion(fEMCALGeo, c2, emCells);
 647
 648         // Clusters not facing frame structures
 649         Bool_t mask2 = MaskFrameCluster(iSupMod2, ieta2);
 650         //if(mask2) printf("Reject eta %d SM %d\n",ieta2, iSupMod2);
 651
 652         if(in1 && in2){
 653
 654           fHmgg->Fill(invmass,p12.Pt());
 655
 656           if(iSupMod1==iSupMod2) {
 657             fHmggSM[iSupMod1]->Fill(invmass,p12.Pt());
 658             fhClusterPairDiffTimeSameSM[iSupMod1]->Fill(p12.E(),time1-time2);
 659           }
 660           else
 661             fHmggDifferentSM ->Fill(invmass,p12.Pt());
 662
 663           // Same sector
 664           Int_t j=0;
 665           for(Int_t i = 0; i < nSM/2; i++){
 666             j=2*i;
 667             if((iSupMod1==j && iSupMod2==j+1) || (iSupMod1==j+1 && iSupMod2==j)) {
 668               fHmggPairSameSectorSM[i]->Fill(invmass,p12.Pt());
 669               fhClusterPairDiffTimeSameSector[i]->Fill(p12.E(),time1-time2);
 670             }
 671           }
 672
 673           // Same side
 674           for(Int_t i = 0; i < nSM-2; i++){
 675             if((iSupMod1==i && iSupMod2==i+2) || (iSupMod1==i+2 && iSupMod2==i)) {
 676               fHmggPairSameSideSM[i]->Fill(invmass,p12.Pt());
 677               fhClusterPairDiffTimeSameSide[i]->Fill(p12.E(),time1-time2);
 678             }
 679           }
 680
 681
 682           if(!mask1 && !mask2){
 683
 684             fHmggMaskFrame->Fill(invmass,p12.Pt());
 685
 686             if(iSupMod1==iSupMod2) fHmggSMMaskFrame[iSupMod1]->Fill(invmass,p12.Pt());
 687             else                   fHmggDifferentSMMaskFrame ->Fill(invmass,p12.Pt());
 688
 689             // Same sector
 690             j=0;
 691             for(Int_t i = 0; i < nSM/2; i++){
 692               j=2*i;
 693               if((iSupMod1==j && iSupMod2==j+1) || (iSupMod1==j+1 && iSupMod2==j)) fHmggPairSameSectorSMMaskFrame[i]->Fill(invmass,p12.Pt());
 694             }
 695
 696             // Same side
 697             for(Int_t i = 0; i < nSM-2; i++){
 698               if((iSupMod1==i && iSupMod2==i+2) || (iSupMod1==i+2 && iSupMod2==i)) fHmggPairSameSideSMMaskFrame[i]->Fill(invmass,p12.Pt());
 699             }
 700
 701           }// Pair not facing frame
 702
 703
 704           if(invmass > 100. && invmass < 160.){//restrict to clusters really close to pi0 peak
 705
 706             //Opening angle of 2 photons
 707             Float_t opangle = p1.Angle(p2.Vect())*TMath::RadToDeg();
 708             //printf("*******>>>>>>>> In PEAK pt %f, angle %f \n",p12.Pt(),opangle);
 709
 710             //Incident angle of each photon
 711             Float_t inangle1 =0., inangle2=0.;
 712             if(gGeoManager){
 713               Float_t posSM1cen[3]={0.,0.,0.};
 714               Float_t depth = fRecoUtils->GetDepth(p1.Energy(),fRecoUtils->GetParticleType(),iSupMod1);
 715               fEMCALGeo->RecalculateTowerPosition(11.5, 23.5, iSupMod1, depth, fRecoUtils->GetMisalTransShiftArray(),fRecoUtils->GetMisalRotShiftArray(),posSM1cen);
 716               Float_t posSM2cen[3]={0.,0.,0.};
 717               depth = fRecoUtils->GetDepth(p2.Energy(),fRecoUtils->GetParticleType(),iSupMod2);
 718               fEMCALGeo->RecalculateTowerPosition(11.5, 23.5, iSupMod2, depth, fRecoUtils->GetMisalTransShiftArray(),fRecoUtils->GetMisalRotShiftArray(),posSM2cen);
 719               //printf("SM1 %d pos (%2.3f,%2.3f,%2.3f) \n",iSupMod1,posSM1cen[0],posSM1cen[1],posSM1cen[2]);
 720               //printf("SM2 %d pos (%2.3f,%2.3f,%2.3f) \n",iSupMod2,posSM2cen[0],posSM2cen[1],posSM2cen[2]);
 721
 722               TVector3 vecSM1cen(posSM1cen[0]-v[0],posSM1cen[1]-v[1],posSM1cen[2]-v[2]);
 723               TVector3 vecSM2cen(posSM2cen[0]-v[0],posSM2cen[1]-v[1],posSM2cen[2]-v[2]);
 724               inangle1 = p1.Angle(vecSM1cen)*TMath::RadToDeg();
 725               inangle2 = p2.Angle(vecSM2cen)*TMath::RadToDeg();
 726               //printf("Incident angle: cluster 1 %2.3f; cluster 2 %2.3f\n",inangle1,inangle2);
 727             }
 728
 729             fHOpeningAngle ->Fill(opangle,p12.Pt());
 730             fHIncidentAngle->Fill(inangle1,p1.Pt());
 731             fHIncidentAngle->Fill(inangle2,p2.Pt());
 732             fHAsymmetry    ->Fill(asym,p12.Pt());
 733
 734             if(iSupMod1==iSupMod2) {
 735               fHOpeningAngleSM[iSupMod1] ->Fill(opangle,p12.Pt());
 736               fHIncidentAngleSM[iSupMod1]->Fill(inangle1,p1.Pt());
 737               fHIncidentAngleSM[iSupMod1]->Fill(inangle2,p2.Pt());
 738               fHAsymmetrySM[iSupMod1]    ->Fill(asym,p12.Pt());
 739             }
 740             else{
 741               fHOpeningAngleDifferentSM  ->Fill(opangle,p12.Pt());
 742               fHIncidentAngleDifferentSM ->Fill(inangle1,p1.Pt());
 743               fHIncidentAngleDifferentSM ->Fill(inangle2,p2.Pt());
 744               fHAsymmetryDifferentSM     ->Fill(asym,p12.Pt());
 745             }
 746
 747             if((iSupMod1==0 && iSupMod2==2) || (iSupMod1==2 && iSupMod2==0)) {
 748               fHOpeningAnglePairSM[0] ->Fill(opangle,p12.Pt());
 749               fHIncidentAnglePairSM[0]->Fill(inangle1,p1.Pt());
 750               fHIncidentAnglePairSM[0]->Fill(inangle2,p2.Pt());
 751               fHAsymmetryPairSM[0]    ->Fill(asym,p12.Pt());
 752
 753             }
 754             if((iSupMod1==1 && iSupMod2==3) || (iSupMod1==3 && iSupMod2==1)) {
 755               fHOpeningAnglePairSM[1] ->Fill(opangle,p12.Pt());
 756               fHIncidentAnglePairSM[1]->Fill(inangle1,p1.Pt());
 757               fHIncidentAnglePairSM[1]->Fill(inangle2,p2.Pt());
 758               fHAsymmetryPairSM[1]    ->Fill(asym,p12.Pt());
 759
 760             }
 761
 762             if((iSupMod1==0 && iSupMod2==1) || (iSupMod1==1 && iSupMod2==0)) {
 763               fHOpeningAnglePairSM[2] ->Fill(opangle,p12.Pt());
 764               fHIncidentAnglePairSM[2]->Fill(inangle1,p1.Pt());
 765               fHIncidentAnglePairSM[2]->Fill(inangle2,p2.Pt());
 766               fHAsymmetryPairSM[2]    ->Fill(asym,p12.Pt());
 767
 768
 769             }
 770             if((iSupMod1==2 && iSupMod2==3) || (iSupMod1==3 && iSupMod2==2)) {
 771               fHOpeningAnglePairSM[3] ->Fill(opangle,p12.Pt());
 772               fHIncidentAnglePairSM[3]->Fill(inangle1,p1.Pt());
 773               fHIncidentAnglePairSM[3]->Fill(inangle2,p2.Pt());
 774               fHAsymmetryPairSM[3]    ->Fill(asym,p12.Pt());
 775             }
 776
 777           }// pair in 100 < mass < 160
 778
 779         }//in acceptance cuts
 780
 781         //In case of filling only channels with second cluster in same SM
 782         if(fSameSM && iSupMod1!=iSupMod2) continue;
 783
 784         if (fGroupNCells == 0){
 785           fHmpi0[iSupMod1][ieta1][iphi1]->Fill(invmass);
 786           fHmpi0[iSupMod2][ieta2][iphi2]->Fill(invmass);
 787
 788           if(invmass > 100. && invmass < 160.){//restrict to clusters really close to pi0 peak
 789             fhTowerDecayPhotonHit      [iSupMod1]->Fill(ieta1,iphi1);
 790             fhTowerDecayPhotonEnergy   [iSupMod1]->Fill(ieta1,iphi1,p1.E());
 791             fhTowerDecayPhotonAsymmetry[iSupMod1]->Fill(ieta1,iphi1,asym);
 792
 793             fhTowerDecayPhotonHit      [iSupMod2]->Fill(ieta2,iphi2);
 794             fhTowerDecayPhotonEnergy   [iSupMod2]->Fill(ieta2,iphi2,p2.E());
 795             fhTowerDecayPhotonAsymmetry[iSupMod2]->Fill(ieta2,iphi2,asym);
 796
 797             if(!mask1)fhTowerDecayPhotonHitMaskFrame[iSupMod1]->Fill(ieta1,iphi1);
 798             if(!mask2)fhTowerDecayPhotonHitMaskFrame[iSupMod2]->Fill(ieta2,iphi2);
 799
 800           }// pair in mass of pi0
 801         }
 802         else  {
 803           //printf("Regroup N %d, eta1 %d, phi1 %d, eta2 %d, phi2 %d \n",fGroupNCells, ieta1, iphi1, ieta2, iphi2);
 804           for (Int_t i = -fGroupNCells; i < fGroupNCells+1; i++) {
 805             for (Int_t j = -fGroupNCells; j < fGroupNCells+1; j++) {
 806               //printf("\t i %d, j %d\n",i,j);
 807               if((ieta1+i >= 0) && (iphi1+j >= 0) && (ieta1+i < 48) && (iphi1+j < 24)){
 808                 //printf("\t \t eta1+i %d, phi1+j %d\n", ieta1+i, iphi1+j);
 809                 fHmpi0[iSupMod1][ieta1+i][iphi1+j]->Fill(invmass);
 810               }
 811               if((ieta2+i >= 0) && (iphi2+j >= 0) && (ieta2+i < 48) && (iphi2+j < 24)){
 812                 //printf("\t \t eta2+i %d, phi2+j %d\n", ieta2+i, iphi2+j);
 813                 fHmpi0[iSupMod2][ieta2+i][iphi2+j]->Fill(invmass);
 814               }
 815             }// j loop
 816           }//i loop
 817         }//group cells
 818
 819         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",
 820                                     iSupMod1,iphi1,ieta1,iSupMod2,iphi2,ieta2,p12.M(),e1i,c1->E(),e2i,c2->E());
 821       }
 822
 823     }
 824
 825   } // end of loop over EMCAL clusters
 826
 827   delete caloClustersArr;
 828
 829   PostData(1,fOutputContainer);
 830
 831 }
 832
 833 //_____________________________________________________
 834 void AliAnalysisTaskEMCALPi0CalibSelection::PrintInfo(){
 835
 836   //Print settings
 837   printf("Cluster cuts: %2.2f < E < %2.2f GeV; number of cells > %d; Assymetry < %1.2f, pair time diff < %2.2f\n",
 838          fEmin,fEmax, fMinNCells, fAsyCut, fDTimeCut) ;
 839   printf("Group %d cells\n", fGroupNCells) ;
 840   printf("Cluster maximal cell away from border at least %d cells\n", fRecoUtils->GetNumberOfCellsFromEMCALBorder()) ;
 841   printf("Histograms: bins %d; energy range: %2.2f < E < %2.2f GeV\n",fNbins,fMinBin,fMaxBin) ;
 842   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",
 843          fRecoUtils->IsBadChannelsRemovalSwitchedOn(),fFilteredInput,fCorrectClusters, fSameSM) ;
 844   printf("EMCAL Geometry name: < %s >, Load Matrices %d\n",fEMCALGeoName.Data(), fLoadMatrices) ;
 845   if(fLoadMatrices) {for(Int_t ism = 0; ism < AliEMCALGeoParams::fgkEMCALModules; ism++) fMatrix[ism]->Print() ; }
 846
 847
 848 }
 849