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