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