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