PWGGA/CaloTrackCorrelations/AliAnaCalorimeterQA.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   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 // Class to check results from simulations or reconstructed real data.
  18 // Fill few histograms and do some checking plots
  19 //
  20 //-- Author: Gustavo Conesa (INFN-LNF)
  21 //_________________________________________________________________________
  22
  23
  24 // --- ROOT system ---
  25 #include <TObjArray.h>
  26 #include <TParticle.h>
  27 #include <TDatabasePDG.h>
  28 #include <TH3F.h>
  29 #include <TObjString.h>
  30
  31 //---- AliRoot system ----
  32 #include "AliAnaCalorimeterQA.h"
  33 #include "AliCaloTrackReader.h"
  34 #include "AliStack.h"
  35 #include "AliVCaloCells.h"
  36 #include "AliFiducialCut.h"
  37 #include "AliVCluster.h"
  38 #include "AliVTrack.h"
  39 #include "AliVEvent.h"
  40 #include "AliVEventHandler.h"
  41 #include "AliAODMCParticle.h"
  42 #include "AliMCAnalysisUtils.h"
  43
  44 // --- Detectors ---
  45 #include "AliPHOSGeoUtils.h"
  46 #include "AliEMCALGeometry.h"
  47
  48 ClassImp(AliAnaCalorimeterQA)
  49
  50 //________________________________________
  51 AliAnaCalorimeterQA::AliAnaCalorimeterQA() :
  52 AliAnaCaloTrackCorrBaseClass(),             fCalorimeter(""),
  53
  54 //Switches
  55 fFillAllCellTimeHisto(kTRUE),
  56 fFillAllPosHisto(kFALSE),              fFillAllPosHisto2(kTRUE),
  57 fFillAllTH12(kFALSE),                  fFillAllTH3(kTRUE),
  58 fFillAllTMHisto(kTRUE),                fFillAllPi0Histo(kTRUE),
  59 fCorrelate(kTRUE),                     fStudyBadClusters(kFALSE),
  60 fStudyClustersAsymmetry(kFALSE),       fStudyExotic(kFALSE),
  61 fStudyWeight(kFALSE),
  62
  63 //Parameters and cuts
  64 fNModules(12),                         fNRCU(2),
  65 fNMaxCols(48),                         fNMaxRows(24),
  66 fTimeCutMin(-10000),                   fTimeCutMax(10000),
  67 fEMCALCellAmpMin(0),                   fPHOSCellAmpMin(0),
  68
  69 // Exotic
  70 fExoNECrossCuts(0),                    fExoECrossCuts(),
  71 fExoNDTimeCuts(0),                     fExoDTimeCuts(),
  72
  73 //Histograms
  74 fhE(0),                                fhPt(0),
  75 fhPhi(0),                              fhEta(0),                               fhEtaPhiE(0),
  76 fhECharged(0),                         fhPtCharged(0),
  77 fhPhiCharged(0),                       fhEtaCharged(0),                        fhEtaPhiECharged(0),
  78
  79 //Invariant mass
  80 fhIM(0 ),                              fhAsym(0),
  81
  82 fhNCellsPerCluster(0),                 fhNCellsPerClusterNoCut(0),             fhNClusters(0),
  83
  84 //Timing
  85 fhClusterTimeEnergy(0),                fhCellTimeSpreadRespectToCellMax(0),
  86 fhCellIdCellLargeTimeSpread(0),        fhClusterPairDiffTimeE(0),
  87 fhClusterMaxCellCloseCellRatio(0),     fhClusterMaxCellCloseCellDiff(0),
  88 fhClusterMaxCellDiff(0),               fhClusterMaxCellDiffNoCut(0),
  89 fhClusterMaxCellDiffAverageTime(0),    fhClusterMaxCellDiffWeightedTime(0),
  90 fhClusterMaxCellECross(0),
  91 fhLambda0(0),                          fhLambda1(0),                           fhDispersion(0),
  92
  93 //bad clusters
  94 fhBadClusterEnergy(0),                 fhBadClusterTimeEnergy(0),
  95 fhBadClusterPairDiffTimeE(0),          fhBadCellTimeSpreadRespectToCellMax(0),
  96 fhBadClusterMaxCellCloseCellRatio(0),  fhBadClusterMaxCellCloseCellDiff(0),    fhBadClusterMaxCellDiff(0),
  97 fhBadClusterMaxCellDiffAverageTime(0), fhBadClusterMaxCellDiffWeightedTime(0),
  98 fhBadClusterMaxCellECross(0),
  99 fhBadClusterDeltaIEtaDeltaIPhiE0(0),   fhBadClusterDeltaIEtaDeltaIPhiE2(0),
 100 fhBadClusterDeltaIEtaDeltaIPhiE6(0),   fhBadClusterDeltaIA(0),
 101
 102 //Position
 103 fhRNCells(0),                          fhXNCells(0),
 104 fhYNCells(0),                          fhZNCells(0),
 105 fhRE(0),                               fhXE(0),
 106 fhYE(0),                               fhZE(0),
 107 fhXYZ(0),
 108 fhRCellE(0),                           fhXCellE(0),
 109 fhYCellE(0),                           fhZCellE(0),
 110 fhXYZCell(0),
 111 fhDeltaCellClusterRNCells(0),          fhDeltaCellClusterXNCells(0),
 112 fhDeltaCellClusterYNCells(0),          fhDeltaCellClusterZNCells(0),
 113 fhDeltaCellClusterRE(0),               fhDeltaCellClusterXE(0),
 114 fhDeltaCellClusterYE(0),               fhDeltaCellClusterZE(0),
 115
 116 // Cells
 117 fhNCells(0),                           fhAmplitude(0),
 118 fhAmpId(0),                            fhEtaPhiAmp(0),
 119 fhTime(0),                             fhTimeVz(0),
 120 fhTimeId(0),                           fhTimeAmp(0),
 121 fhCellECross(0),
 122 fhCaloCorrNClusters(0),                fhCaloCorrEClusters(0),
 123 fhCaloCorrNCells(0),                   fhCaloCorrECells(0),
 124 fhCaloV0SCorrNClusters(0),             fhCaloV0SCorrEClusters(0),
 125 fhCaloV0SCorrNCells(0),                fhCaloV0SCorrECells(0),
 126 fhCaloV0MCorrNClusters(0),             fhCaloV0MCorrEClusters(0),
 127 fhCaloV0MCorrNCells(0),                fhCaloV0MCorrECells(0),
 128 fhCaloTrackMCorrNClusters(0),          fhCaloTrackMCorrEClusters(0),
 129 fhCaloTrackMCorrNCells(0),             fhCaloTrackMCorrECells(0),
 130
 131 //Super-Module dependent histgrams
 132 fhEMod(0),                             fhAmpMod(0),                            fhTimeMod(0),
 133 fhNClustersMod(0),                     fhNCellsMod(0),
 134 fhNCellsPerClusterMod(0),              fhNCellsPerClusterModNoCut(0),
 135
 136 fhGridCells(0),                        fhGridCellsE(0),                        fhGridCellsTime(0),
 137 fhTimeAmpPerRCU(0),                    fhIMMod(0),
 138
 139 // Weight studies
 140 fhECellClusterRatio(0),                fhECellClusterLogRatio(0),
 141 fhEMaxCellClusterRatio(0),             fhEMaxCellClusterLogRatio(0),
 142
 143 // MC and reco
 144 fhRecoMCE(),                           fhRecoMCPhi(),                          fhRecoMCEta(),
 145 fhRecoMCDeltaE(),                      fhRecoMCRatioE(),
 146 fhRecoMCDeltaPhi(),                    fhRecoMCDeltaEta(),
 147
 148 // MC only
 149 fhGenMCE(),                            fhGenMCEtaPhi(),
 150 fhGenMCAccE(),                         fhGenMCAccEtaPhi(),
 151
 152 //matched MC
 153 fhEMVxyz(0),                           fhEMR(0),
 154 fhHaVxyz(0),                           fhHaR(0),
 155 fh1EOverP(0),                          fh2dR(0),
 156 fh2EledEdx(0),                         fh2MatchdEdx(0),
 157 fhMCEle1EOverP(0),                     fhMCEle1dR(0),                          fhMCEle2MatchdEdx(0),
 158 fhMCChHad1EOverP(0),                   fhMCChHad1dR(0),                        fhMCChHad2MatchdEdx(0),
 159 fhMCNeutral1EOverP(0),                 fhMCNeutral1dR(0),                      fhMCNeutral2MatchdEdx(0), fh1EOverPR02(0),
 160 fhMCEle1EOverPR02(0),                  fhMCChHad1EOverPR02(0),                 fhMCNeutral1EOverPR02(0)
 161 {
 162   //Default Ctor
 163
 164   //Weight studies
 165   for(Int_t i =0; i < 14; i++){
 166     fhLambda0ForW0[i] = 0;
 167     //fhLambda1ForW0[i] = 0;
 168
 169     for(Int_t j = 0; j < 5; j++){
 170       fhLambda0ForW0MC[i][j] = 0;
 171       //fhLambda1ForW0MC[i][j] = 0;
 172     }
 173
 174   }
 175
 176   //Cluster size
 177   fhDeltaIEtaDeltaIPhiE0[0] = 0 ;         fhDeltaIEtaDeltaIPhiE2[0] = 0;          fhDeltaIEtaDeltaIPhiE6[0] = 0;
 178   fhDeltaIEtaDeltaIPhiE0[1] = 0 ;         fhDeltaIEtaDeltaIPhiE2[1] = 0;          fhDeltaIEtaDeltaIPhiE6[1] = 0;
 179   fhDeltaIA[0]              = 0 ;         fhDeltaIAL0[0]            = 0;          fhDeltaIAL1[0]            = 0;
 180   fhDeltaIA[1]              = 0 ;         fhDeltaIAL0[1]            = 0;          fhDeltaIAL1[1]            = 0;
 181   fhDeltaIANCells[0]        = 0 ;         fhDeltaIANCells[1]        = 0;
 182   fhDeltaIAMC[0]            = 0 ;         fhDeltaIAMC[1]            = 0;
 183   fhDeltaIAMC[2]            = 0 ;         fhDeltaIAMC[3]            = 0;
 184
 185   // Exotic
 186   for (Int_t ie = 0; ie < 10 ; ie++)
 187   {
 188     fhExoDTime[ie] = 0;
 189     for (Int_t idt = 0; idt < 5 ; idt++)
 190     {
 191       fhExoNCell    [ie][idt] = 0;
 192       fhExoL0       [ie][idt] = 0;
 193       fhExoECross   [ie][idt] = 0;
 194       fhExoTime     [ie][idt] = 0;
 195     }
 196   }
 197
 198   // MC
 199
 200   for(Int_t i = 0; i < 6; i++){
 201
 202     fhRecoMCE[i][0]         = 0; fhRecoMCE[i][1]        = 0;
 203     fhRecoMCPhi[i][0]       = 0; fhRecoMCPhi[i][1]      = 0;
 204     fhRecoMCEta[i][0]       = 0; fhRecoMCEta[i][1]      = 0;
 205     fhRecoMCDeltaE[i][0]    = 0; fhRecoMCDeltaE[i][1]   = 0;
 206     fhRecoMCRatioE[i][0]    = 0; fhRecoMCRatioE[i][1]   = 0;
 207     fhRecoMCDeltaPhi[i][0]  = 0; fhRecoMCDeltaPhi[i][1] = 0;
 208     fhRecoMCDeltaEta[i][0]  = 0; fhRecoMCDeltaEta[i][1] = 0;
 209
 210   }
 211
 212   //Initialize parameters
 213   InitParameters();
 214 }
 215
 216 //____________________________________________________________________________________________________________________________
 217 void AliAnaCalorimeterQA::BadClusterHistograms(AliVCluster* clus, const TObjArray *caloClusters, AliVCaloCells * cells,
 218                                                const Int_t absIdMax, const Double_t maxCellFraction, const Float_t eCrossFrac,
 219                                                const Double_t tmax,  Double_t timeAverages[2]
 220                                                )
 221 {
 222   //Bad cluster histograms
 223
 224   //  printf("AliAnaCalorimeterQA::BadClusterHistograms() - Event %d - Calorimeter %s \n \t  E %f, n cells %d, max cell absId %d, maxCellFrac %f\n",
 225   //         GetReader()->GetEventNumber(), fCalorimeter.Data(),
 226   //         clus->E(),clus->GetNCells(),absIdMax,maxCellFraction);
 227
 228   fhBadClusterEnergy     ->Fill(clus->E());
 229   Double_t tof = clus->GetTOF()*1.e9;
 230   fhBadClusterTimeEnergy   ->Fill(clus->E(),tof);
 231   fhBadClusterMaxCellDiff  ->Fill(clus->E(),maxCellFraction);
 232   fhBadClusterMaxCellECross->Fill(clus->E(),eCrossFrac);
 233
 234   if(fStudyClustersAsymmetry) ClusterAsymmetryHistograms(clus,absIdMax,kFALSE);
 235
 236   //Clusters in event time differencem bad minus good
 237
 238   for(Int_t iclus2 = 0; iclus2 < caloClusters->GetEntriesFast(); iclus2++ ){
 239
 240     AliVCluster* clus2 =  (AliVCluster*)caloClusters->At(iclus2);
 241
 242     if(clus->GetID()==clus2->GetID()) continue;
 243
 244     Float_t maxCellFraction2 = 0.;
 245     Int_t absIdMax2 = GetCaloUtils()->GetMaxEnergyCell(cells, clus2,maxCellFraction2);
 246     if(IsGoodCluster(absIdMax2,cells)){
 247       Double_t tof2   = clus2->GetTOF()*1.e9;
 248       fhBadClusterPairDiffTimeE  ->Fill(clus->E(), (tof-tof2));
 249     }
 250
 251   } // loop
 252
 253   // Max cell compared to other cells in cluster
 254   if(fFillAllCellTimeHisto)
 255   {
 256     fhBadClusterMaxCellDiffAverageTime      ->Fill(clus->E(),tmax-timeAverages[0]);
 257     fhBadClusterMaxCellDiffWeightedTime     ->Fill(clus->E(),tmax-timeAverages[1]);
 258   }
 259
 260   for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
 261   {
 262     Int_t absId  = clus->GetCellsAbsId()[ipos];
 263     if(absId!=absIdMax && cells->GetCellAmplitude(absIdMax) > 0.01)
 264     {
 265       Float_t frac = cells->GetCellAmplitude(absId)/cells->GetCellAmplitude(absIdMax);
 266
 267       fhBadClusterMaxCellCloseCellRatio->Fill(clus->E(),frac);
 268       fhBadClusterMaxCellCloseCellDiff ->Fill(clus->E(),cells->GetCellAmplitude(absIdMax)-cells->GetCellAmplitude(absId));
 269
 270       if(fFillAllCellTimeHisto)
 271       {
 272         Double_t time  = cells->GetCellTime(absId);
 273         GetCaloUtils()->RecalibrateCellTime(time, fCalorimeter, absId,GetReader()->GetInputEvent()->GetBunchCrossNumber());
 274
 275         Float_t diff = (tmax-time*1e9);
 276         fhBadCellTimeSpreadRespectToCellMax->Fill(clus->E(), diff);
 277
 278       }
 279     }// Not max
 280   }//loop
 281
 282 }
 283
 284 //______________________________________________________________________
 285 void AliAnaCalorimeterQA::CalculateAverageTime(AliVCluster *clus,
 286                                                AliVCaloCells* cells,
 287                                                Double_t timeAverages[2])
 288 {
 289   // Calculate time averages and weights
 290
 291   // First recalculate energy in case non linearity was applied
 292   Float_t  energy = 0;
 293   Float_t  ampMax = 0, amp = 0;
 294   Int_t    absIdMax =-1;
 295   for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
 296   {
 297     Int_t id       = clus->GetCellsAbsId()[ipos];
 298
 299     //Recalibrate cell energy if needed
 300     amp = cells->GetCellAmplitude(id);
 301     GetCaloUtils()->RecalibrateCellAmplitude(amp,fCalorimeter, id);
 302
 303     energy    += amp;
 304
 305     if(amp> ampMax)
 306     {
 307       ampMax   = amp;
 308       absIdMax = id;
 309     }
 310
 311   } // energy loop
 312
 313   // Calculate average time of cells in cluster and weighted average
 314   Double_t aTime  = 0;
 315   Double_t wTime  = 0;
 316   Float_t  wTot   = 0;
 317   Double_t time   = 0;
 318   Int_t    id     =-1;
 319   Double_t w      = 0;
 320   Int_t    ncells = clus->GetNCells();
 321   for (Int_t ipos = 0; ipos < ncells; ipos++)
 322   {
 323     id   = clus ->GetCellsAbsId()[ipos];
 324     amp  = cells->GetCellAmplitude(id);
 325     time = cells->GetCellTime(id);
 326
 327     //Recalibrate energy and time
 328     GetCaloUtils()->RecalibrateCellAmplitude(amp , fCalorimeter, id);
 329     GetCaloUtils()->RecalibrateCellTime     (time, fCalorimeter, id, GetReader()->GetInputEvent()->GetBunchCrossNumber());
 330
 331     w      = GetCaloUtils()->GetEMCALRecoUtils()->GetCellWeight(cells->GetCellAmplitude(id),energy);
 332     aTime += time*1e9;
 333     wTime += time*1e9 * w;
 334     wTot  += w;
 335
 336   }
 337
 338   if(ncells > 0) aTime /= ncells;
 339   else           aTime  = 0;
 340
 341   if(wTot   > 0) wTime /= wTot;
 342   else           wTime  = 0;
 343
 344   timeAverages[0] = aTime;
 345   timeAverages[1] = wTime;
 346
 347 }
 348
 349 //____________________________________________________________
 350 void AliAnaCalorimeterQA::CellHistograms(AliVCaloCells *cells)
 351 {
 352   // Plot histograms related to cells only
 353
 354   Int_t ncells = cells->GetNumberOfCells();
 355
 356   if(GetDebug() > 0)
 357     printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - %s cell entries %d\n", fCalorimeter.Data(), ncells );
 358
 359   //Init arrays and used variables
 360   Int_t *nCellsInModule = new Int_t[fNModules];
 361   for(Int_t imod = 0; imod < fNModules; imod++ ) nCellsInModule[imod] = 0;
 362
 363   Int_t    icol   = -1;
 364   Int_t    irow   = -1;
 365   Int_t    iRCU   = -1;
 366   Float_t  amp    = 0.;
 367   Double_t time   = 0.;
 368   Int_t    id     = -1;
 369   Float_t  recalF = 1.;
 370   Int_t    bc     = GetReader()->GetInputEvent()->GetBunchCrossNumber();
 371
 372   for (Int_t iCell = 0; iCell < cells->GetNumberOfCells(); iCell++) {
 373     if(GetDebug() > 2)
 374       printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - Cell : amp %f, absId %d \n", cells->GetAmplitude(iCell), cells->GetCellNumber(iCell));
 375     Int_t nModule = GetModuleNumberCellIndexes(cells->GetCellNumber(iCell),fCalorimeter, icol, irow, iRCU);
 376     if(GetDebug() > 2)
 377       printf("\t module %d, column %d, row %d \n", nModule,icol,irow);
 378
 379     if(nModule < fNModules)
 380     {
 381       //Check if the cell is a bad channel
 382       if(GetCaloUtils()->IsBadChannelsRemovalSwitchedOn()){
 383         if(fCalorimeter=="EMCAL")
 384         {
 385           if(GetCaloUtils()->GetEMCALChannelStatus(nModule,icol,irow)) continue;
 386         }
 387         else
 388         {
 389           if(GetCaloUtils()->GetPHOSChannelStatus(nModule,icol,irow) ) continue;
 390         }
 391       } // use bad channel map
 392
 393       amp     = cells->GetAmplitude(iCell)*recalF;
 394       time    = cells->GetTime(iCell);
 395       id      = cells->GetCellNumber(iCell);
 396
 397       // Amplitude recalibration if set
 398       GetCaloUtils()->RecalibrateCellAmplitude(amp,  fCalorimeter, id);
 399
 400       // Time recalibration if set
 401       GetCaloUtils()->RecalibrateCellTime     (time, fCalorimeter, id, GetReader()->GetInputEvent()->GetBunchCrossNumber());
 402
 403       //Transform time to ns
 404       time *= 1.0e9;
 405
 406       if(time < fTimeCutMin || time > fTimeCutMax)
 407       {
 408           if(GetDebug() > 0 )
 409             printf("AliAnaCalorimeterQA - Remove cell with Time %f\n",time);
 410           continue;
 411       }
 412
 413       // Remove exotic cells, defined only for EMCAL
 414       if(fCalorimeter=="EMCAL" &&
 415          GetCaloUtils()->GetEMCALRecoUtils()->IsExoticCell(id, cells, bc)) continue;
 416
 417
 418       fhAmplitude->Fill(amp);
 419       fhAmpId    ->Fill(amp,id);
 420       fhAmpMod   ->Fill(amp,nModule);
 421
 422       if ((fCalorimeter=="EMCAL" && amp > fEMCALCellAmpMin) ||
 423           (fCalorimeter=="PHOS"  && amp > fPHOSCellAmpMin )   )
 424       {
 425
 426         //E cross for exotic cells
 427         if(amp > 0.01) fhCellECross->Fill(amp,1-GetECross(id,cells)/amp);
 428
 429         nCellsInModule[nModule]++ ;
 430
 431         Int_t icols = icol;
 432         Int_t irows = irow;
 433
 434         if(fCalorimeter=="EMCAL")
 435         {
 436           icols = (nModule % 2) ? icol + fNMaxCols : icol;
 437           if(nModule < 10 )
 438             irows = irow + fNMaxRows       * Int_t(nModule / 2);
 439           else // 1/3 SM
 440             irows = irow + (fNMaxRows / 3) * Int_t(nModule / 2);
 441         }
 442         else
 443         {
 444           irows = irow + fNMaxRows * nModule;
 445         }
 446
 447         fhGridCells ->Fill(icols,irows);
 448         fhGridCellsE->Fill(icols,irows,amp);
 449
 450         if(fFillAllCellTimeHisto)
 451         {
 452           //printf("%s: time %g\n",fCalorimeter.Data(), time);
 453
 454           Double_t v[3] = {0,0,0}; //vertex ;
 455           GetReader()->GetVertex(v);
 456           if(amp > 0.5) fhTimeVz   ->Fill(TMath::Abs(v[2]),time);
 457
 458           fhTime     ->Fill(time);
 459           fhTimeId   ->Fill(time,id);
 460           fhTimeAmp  ->Fill(amp,time);
 461           fhGridCellsTime->Fill(icols,irows,time);
 462           fhTimeMod  ->Fill(time,nModule);
 463           fhTimeAmpPerRCU  [nModule*fNRCU+iRCU]->Fill(amp, time);
 464
 465         }
 466       }
 467
 468       //Get Eta-Phi position of Cell
 469       if(fFillAllPosHisto)
 470       {
 471         if(fCalorimeter=="EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet()){
 472           Float_t celleta = 0.;
 473           Float_t cellphi = 0.;
 474           GetEMCALGeometry()->EtaPhiFromIndex(id, celleta, cellphi);
 475
 476           fhEtaPhiAmp->Fill(celleta,cellphi,amp);
 477           Double_t cellpos[] = {0, 0, 0};
 478           GetEMCALGeometry()->GetGlobal(id, cellpos);
 479           fhXCellE->Fill(cellpos[0],amp)  ;
 480           fhYCellE->Fill(cellpos[1],amp)  ;
 481           fhZCellE->Fill(cellpos[2],amp)  ;
 482           Float_t rcell = TMath::Sqrt(cellpos[0]*cellpos[0]+cellpos[1]*cellpos[1]);//+cellpos[2]*cellpos[2]);
 483           fhRCellE->Fill(rcell,amp)  ;
 484           fhXYZCell->Fill(cellpos[0],cellpos[1],cellpos[2])  ;
 485         }//EMCAL Cells
 486         else if(fCalorimeter=="PHOS" && GetCaloUtils()->IsPHOSGeoMatrixSet()){
 487           TVector3 xyz;
 488           Int_t relId[4], module;
 489           Float_t xCell, zCell;
 490
 491           GetPHOSGeometry()->AbsToRelNumbering(id,relId);
 492           module = relId[0];
 493           GetPHOSGeometry()->RelPosInModule(relId,xCell,zCell);
 494           GetPHOSGeometry()->Local2Global(module,xCell,zCell,xyz);
 495           Float_t rcell = TMath::Sqrt(xyz.X()*xyz.X()+xyz.Y()*xyz.Y());
 496           fhXCellE ->Fill(xyz.X(),amp)  ;
 497           fhYCellE ->Fill(xyz.Y(),amp)  ;
 498           fhZCellE ->Fill(xyz.Z(),amp)  ;
 499           fhRCellE ->Fill(rcell  ,amp)  ;
 500           fhXYZCell->Fill(xyz.X(),xyz.Y(),xyz.Z())  ;
 501         }//PHOS cells
 502       }//fill cell position histograms
 503
 504       if     (fCalorimeter=="EMCAL" && amp > fEMCALCellAmpMin) ncells ++ ;
 505       else if(fCalorimeter=="PHOS"  && amp > fPHOSCellAmpMin)  ncells ++ ;
 506       //else
 507       //  printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - no %s CELLS passed the analysis cut\n",fCalorimeter.Data());
 508     }//nmodules
 509   }//cell loop
 510
 511   if(ncells > 0 )fhNCells->Fill(ncells) ; //fill the cells after the cut
 512
 513   //Number of cells per module
 514   for(Int_t imod = 0; imod < fNModules; imod++ ) {
 515
 516     if(GetDebug() > 1)
 517       printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - module %d calo %s cells %d\n", imod, fCalorimeter.Data(), nCellsInModule[imod]);
 518
 519     fhNCellsMod->Fill(nCellsInModule[imod],imod) ;
 520
 521   }
 522
 523   delete [] nCellsInModule;
 524
 525 }
 526
 527 //__________________________________________________________________________
 528 void AliAnaCalorimeterQA::CellInClusterPositionHistograms(AliVCluster* clus)
 529 {
 530   // Fill histograms releated to cell position
 531
 532
 533   Int_t nCaloCellsPerCluster = clus->GetNCells();
 534   UShort_t * indexList = clus->GetCellsAbsId();
 535   Float_t pos[3];
 536   clus->GetPosition(pos);
 537   Float_t clEnergy = clus->E();
 538
 539   //Loop on cluster cells
 540   for (Int_t ipos = 0; ipos < nCaloCellsPerCluster; ipos++) {
 541
 542     //  printf("Index %d\n",ipos);
 543     Int_t absId  = indexList[ipos];
 544
 545     //Get position of cell compare to cluster
 546
 547     if(fCalorimeter=="EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet()){
 548
 549       Double_t cellpos[] = {0, 0, 0};
 550       GetEMCALGeometry()->GetGlobal(absId, cellpos);
 551
 552       fhDeltaCellClusterXNCells->Fill(pos[0]-cellpos[0],nCaloCellsPerCluster) ;
 553       fhDeltaCellClusterYNCells->Fill(pos[1]-cellpos[1],nCaloCellsPerCluster) ;
 554       fhDeltaCellClusterZNCells->Fill(pos[2]-cellpos[2],nCaloCellsPerCluster) ;
 555
 556       fhDeltaCellClusterXE->Fill(pos[0]-cellpos[0],clEnergy)  ;
 557       fhDeltaCellClusterYE->Fill(pos[1]-cellpos[1],clEnergy)  ;
 558       fhDeltaCellClusterZE->Fill(pos[2]-cellpos[2],clEnergy)  ;
 559
 560       Float_t r     = TMath::Sqrt(pos[0]    *pos[0]     + pos[1]    * pos[1]    );
 561       Float_t rcell = TMath::Sqrt(cellpos[0]*cellpos[0] + cellpos[1]* cellpos[1]);
 562
 563       fhDeltaCellClusterRNCells->Fill(r-rcell, nCaloCellsPerCluster) ;
 564       fhDeltaCellClusterRE     ->Fill(r-rcell, clEnergy)  ;
 565
 566     }//EMCAL and its matrices are available
 567     else if(fCalorimeter=="PHOS" && GetCaloUtils()->IsPHOSGeoMatrixSet()){
 568       TVector3 xyz;
 569       Int_t relId[4], module;
 570       Float_t xCell, zCell;
 571
 572       GetPHOSGeometry()->AbsToRelNumbering(absId,relId);
 573       module = relId[0];
 574       GetPHOSGeometry()->RelPosInModule(relId,xCell,zCell);
 575       GetPHOSGeometry()->Local2Global(module,xCell,zCell,xyz);
 576
 577       fhDeltaCellClusterXNCells->Fill(pos[0]-xyz.X(),nCaloCellsPerCluster) ;
 578       fhDeltaCellClusterYNCells->Fill(pos[1]-xyz.Y(),nCaloCellsPerCluster) ;
 579       fhDeltaCellClusterZNCells->Fill(pos[2]-xyz.Z(),nCaloCellsPerCluster) ;
 580
 581       fhDeltaCellClusterXE->Fill(pos[0]-xyz.X(),clEnergy)  ;
 582       fhDeltaCellClusterYE->Fill(pos[1]-xyz.Y(),clEnergy)  ;
 583       fhDeltaCellClusterZE->Fill(pos[2]-xyz.Z(),clEnergy)  ;
 584
 585       Float_t r     = TMath::Sqrt(pos[0]  * pos[0]  + pos[1]  * pos[1] );
 586       Float_t rcell = TMath::Sqrt(xyz.X() * xyz.X() + xyz.Y() * xyz.Y());
 587
 588       fhDeltaCellClusterRNCells->Fill(r-rcell, nCaloCellsPerCluster) ;
 589       fhDeltaCellClusterRE     ->Fill(r-rcell, clEnergy)  ;
 590
 591     }//PHOS and its matrices are available
 592   }// cluster cell loop
 593 }
 594
 595 //___________________________________________________________________________________________
 596 void AliAnaCalorimeterQA::ClusterAsymmetryHistograms(AliVCluster* clus, const Int_t absIdMax,
 597                                                      const Bool_t goodCluster)
 598 {
 599   // Study the shape of the cluster in cell units terms
 600
 601   //No use to study clusters with less than 4 cells
 602   if(clus->GetNCells() <=3 ) return;
 603
 604   Int_t dIeta = 0;
 605   Int_t dIphi = 0;
 606
 607   Int_t ietaMax=-1; Int_t iphiMax = 0; Int_t rcuMax = 0;
 608   Int_t smMax = GetModuleNumberCellIndexes(absIdMax,fCalorimeter, ietaMax, iphiMax, rcuMax);
 609
 610   for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++) {
 611
 612     Int_t absId = clus->GetCellsAbsId()[ipos];
 613
 614     Int_t ieta=-1; Int_t iphi = 0; Int_t rcu = 0;
 615     Int_t sm = GetModuleNumberCellIndexes(absId,fCalorimeter, ieta, iphi, rcu);
 616
 617     if(dIphi < TMath::Abs(iphi-iphiMax)) dIphi = TMath::Abs(iphi-iphiMax);
 618
 619     if(smMax==sm){
 620       if(dIeta < TMath::Abs(ieta-ietaMax)) dIeta = TMath::Abs(ieta-ietaMax);
 621     }
 622     else {
 623       Int_t ietaShift    = ieta;
 624       Int_t ietaMaxShift = ietaMax;
 625       if (ieta > ietaMax)  ietaMaxShift+=48;
 626       else                 ietaShift   +=48;
 627       if(dIeta < TMath::Abs(ietaShift-ietaMaxShift)) dIeta = TMath::Abs(ietaShift-ietaMaxShift);
 628     }
 629
 630   }// fill cell-cluster histogram loop
 631
 632
 633   Float_t dIA = 1.*(dIphi-dIeta)/(dIeta+dIphi);
 634
 635   if(goodCluster)
 636   {
 637
 638     // Was cluster matched?
 639     Bool_t matched = GetCaloPID()->IsTrackMatched(clus,GetCaloUtils(),GetReader()->GetInputEvent());
 640
 641     if     (clus->E() < 2 ) fhDeltaIEtaDeltaIPhiE0[matched]->Fill(dIeta,dIphi);
 642     else if(clus->E() < 6 ) fhDeltaIEtaDeltaIPhiE2[matched]->Fill(dIeta,dIphi);
 643     else                    fhDeltaIEtaDeltaIPhiE6[matched]->Fill(dIeta,dIphi);
 644
 645     fhDeltaIA[matched]->Fill(clus->E(),dIA);
 646
 647     if(clus->E() > 0.5){
 648
 649       fhDeltaIAL0[matched]    ->Fill(clus->GetM02(),dIA);
 650       fhDeltaIAL1[matched]    ->Fill(clus->GetM20(),dIA);
 651       fhDeltaIANCells[matched]->Fill(clus->GetNCells(),dIA);
 652
 653     }
 654
 655     // Origin of  clusters
 656     Int_t  nLabel = clus->GetNLabels();
 657     Int_t* labels = clus->GetLabels();
 658     if(IsDataMC()){
 659       Int_t tag = GetMCAnalysisUtils()->CheckOrigin(labels,nLabel, GetReader(),0);
 660       if(   GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton) &&
 661          !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0)    &&
 662          !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta)    &&
 663          !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)        ){
 664         fhDeltaIAMC[0]->Fill(clus->E(),dIA);//Pure Photon
 665       }
 666       else if ( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCElectron) &&
 667                !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)  ){
 668         fhDeltaIAMC[1]->Fill(clus->E(),dIA);//Pure electron
 669       }
 670       else if ( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton) &&
 671                GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)  ){
 672         fhDeltaIAMC[2]->Fill(clus->E(),dIA);//Converted cluster
 673       }
 674       else if(!GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton)){
 675         fhDeltaIAMC[3]->Fill(clus->E(),dIA);//Hadrons
 676       }
 677
 678     }  // MC
 679   } // good cluster
 680   else
 681   {
 682     if     (clus->E() < 2 ) fhBadClusterDeltaIEtaDeltaIPhiE0->Fill(dIeta,dIphi);
 683     else if(clus->E() < 6 ) fhBadClusterDeltaIEtaDeltaIPhiE2->Fill(dIeta,dIphi);
 684     else                    fhBadClusterDeltaIEtaDeltaIPhiE6->Fill(dIeta,dIphi);
 685
 686     fhBadClusterDeltaIA->Fill(clus->E(),dIA);
 687
 688   }
 689 }
 690
 691 //_________________________________________________________________________________________________________________________
 692 void AliAnaCalorimeterQA::ClusterHistograms(AliVCluster* clus,const TObjArray *caloClusters, AliVCaloCells * cells,
 693                                             const Int_t absIdMax, const Double_t maxCellFraction, const Float_t eCrossFrac,
 694                                             const Double_t tmax,  Double_t timeAverages[2])
 695 {
 696   //Fill CaloCluster related histograms
 697
 698   Double_t tof = clus->GetTOF()*1.e9;
 699
 700   fhLambda0             ->Fill(clus->E(),clus->GetM02());
 701   fhLambda1             ->Fill(clus->E(),clus->GetM20());
 702   fhDispersion          ->Fill(clus->E(),clus->GetDispersion());
 703
 704   fhClusterMaxCellDiff  ->Fill(clus->E(),maxCellFraction);
 705   fhClusterMaxCellECross->Fill(clus->E(),eCrossFrac);
 706   fhClusterTimeEnergy   ->Fill(clus->E(),tof);
 707
 708   if(fStudyClustersAsymmetry) ClusterAsymmetryHistograms(clus,absIdMax,kTRUE);
 709
 710   //Clusters in event time difference
 711   for(Int_t iclus2 = 0; iclus2 < caloClusters->GetEntriesFast(); iclus2++ ){
 712
 713     AliVCluster* clus2 =  (AliVCluster*) caloClusters->At(iclus2);
 714
 715     if(clus->GetID()==clus2->GetID()) continue;
 716
 717     if(clus->GetM02() > 0.01 && clus2->GetM02() > 0.01)
 718     {
 719       Double_t tof2   = clus2->GetTOF()*1.e9;
 720       fhClusterPairDiffTimeE  ->Fill(clus->E(), tof-tof2);
 721     }
 722   }
 723
 724   Int_t    nModule = GetModuleNumber(clus);
 725   Int_t    nCaloCellsPerCluster = clus->GetNCells();
 726
 727   if(nCaloCellsPerCluster > 1){
 728
 729     // check time of cells respect to max energy cell
 730
 731     if(fFillAllCellTimeHisto)
 732     {
 733       fhClusterMaxCellDiffAverageTime      ->Fill(clus->E(),tmax-timeAverages[0]);
 734       fhClusterMaxCellDiffWeightedTime     ->Fill(clus->E(),tmax-timeAverages[1]);
 735     }
 736
 737     for (Int_t ipos = 0; ipos < nCaloCellsPerCluster; ipos++)
 738     {
 739       Int_t absId  = clus->GetCellsAbsId()[ipos];
 740       if(absId == absIdMax || cells->GetCellAmplitude(absIdMax) < 0.01) continue;
 741
 742       Float_t frac = cells->GetCellAmplitude(absId)/cells->GetCellAmplitude(absIdMax);
 743       fhClusterMaxCellCloseCellRatio->Fill(clus->E(),frac);
 744       fhClusterMaxCellCloseCellDiff ->Fill(clus->E(),cells->GetCellAmplitude(absIdMax)-cells->GetCellAmplitude(absId));
 745
 746       if(fFillAllCellTimeHisto)
 747       {
 748         Double_t time  = cells->GetCellTime(absId);
 749         GetCaloUtils()->RecalibrateCellTime(time, fCalorimeter, absId,GetReader()->GetInputEvent()->GetBunchCrossNumber());
 750
 751         Float_t diff = (tmax-time*1.0e9);
 752         fhCellTimeSpreadRespectToCellMax->Fill(clus->E(), diff);
 753         if(TMath::Abs(TMath::Abs(diff) > 100) && clus->E() > 1 ) fhCellIdCellLargeTimeSpread->Fill(absId);
 754       }
 755
 756     }// fill cell-cluster histogram loop
 757
 758   }//check time and energy of cells respect to max energy cell if cluster of more than 1 cell
 759
 760
 761   // Get vertex for photon momentum calculation and event selection
 762   Double_t v[3] = {0,0,0}; //vertex ;
 763   //GetReader()->GetVertex(v); //
 764
 765   TLorentzVector mom  ;
 766   clus->GetMomentum(mom,v);
 767
 768   Float_t e   = mom.E();
 769   Float_t pt  = mom.Pt();
 770   Float_t eta = mom.Eta();
 771   Float_t phi = mom.Phi();
 772   if(phi < 0) phi +=TMath::TwoPi();
 773
 774   if(GetDebug() > 0) {
 775     printf("AliAnaCalorimeterQA::ClusterHistograms() - cluster: E %2.3f, pT %2.3f, eta %2.3f, phi %2.3f \n",e,pt,eta,phi*TMath::RadToDeg());
 776   }
 777
 778   fhE     ->Fill(e);
 779   if(nModule >=0 && nModule < fNModules) fhEMod->Fill(e,nModule);
 780   if(fFillAllTH12){
 781     fhPt     ->Fill(pt);
 782     fhPhi    ->Fill(phi);
 783     fhEta    ->Fill(eta);
 784   }
 785
 786   if(fFillAllTH3)
 787     fhEtaPhiE->Fill(eta,phi,e);
 788
 789   //Cells per cluster
 790   fhNCellsPerCluster   ->Fill(e, nCaloCellsPerCluster);
 791
 792   //Position
 793   if(fFillAllPosHisto2){
 794
 795     Float_t pos[3] ;
 796     clus->GetPosition(pos);
 797
 798     fhXE     ->Fill(pos[0],e);
 799     fhYE     ->Fill(pos[1],e);
 800     fhZE     ->Fill(pos[2],e);
 801     if(fFillAllTH3)
 802       fhXYZ    ->Fill(pos[0], pos[1],pos[2]);
 803
 804     fhXNCells->Fill(pos[0],nCaloCellsPerCluster);
 805     fhYNCells->Fill(pos[1],nCaloCellsPerCluster);
 806     fhZNCells->Fill(pos[2],nCaloCellsPerCluster);
 807     Float_t rxyz = TMath::Sqrt(pos[0]*pos[0]+pos[1]*pos[1]);//+pos[2]*pos[2]);
 808     fhRE     ->Fill(rxyz,e);
 809     fhRNCells->Fill(rxyz  ,nCaloCellsPerCluster);
 810   }
 811
 812   if(nModule >=0 && nModule < fNModules) fhNCellsPerClusterMod[nModule]->Fill(e, nCaloCellsPerCluster);
 813
 814 }
 815
 816 //____________________________________________________________________________
 817 void AliAnaCalorimeterQA::ClusterLoopHistograms(const TObjArray *caloClusters,
 818                                                 AliVCaloCells* cells)
 819 {
 820   // Fill clusters related histograms
 821
 822   TLorentzVector mom  ;
 823   Int_t  nLabel                = 0  ;
 824   Int_t *labels                = 0x0;
 825   Int_t  nCaloClusters         = caloClusters->GetEntriesFast() ;
 826   Int_t  nCaloClustersAccepted = 0  ;
 827   Int_t  nCaloCellsPerCluster  = 0  ;
 828   Bool_t matched               = kFALSE;
 829   Int_t  nModule               =-1  ;
 830
 831   // Get vertex for photon momentum calculation and event selection
 832   Double_t v[3] = {0,0,0}; //vertex ;
 833   //GetReader()->GetVertex(v);
 834
 835   Int_t *nClustersInModule     = new Int_t[fNModules];
 836   for(Int_t imod = 0; imod < fNModules; imod++ ) nClustersInModule[imod] = 0;
 837
 838   if(GetDebug() > 0)
 839     printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - In %s there are %d clusters \n", fCalorimeter.Data(), nCaloClusters);
 840
 841   // Loop over CaloClusters
 842   for(Int_t iclus = 0; iclus < nCaloClusters; iclus++){
 843
 844     if(GetDebug() > 0)
 845       printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - cluster: %d/%d, data %d \n",
 846              iclus+1,nCaloClusters,GetReader()->GetDataType());
 847
 848     AliVCluster* clus =  (AliVCluster*)caloClusters->At(iclus);
 849
 850     // Get the fraction of the cluster energy that carries the cell with highest energy and its absId
 851     Float_t maxCellFraction = 0.;
 852     Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(cells, clus,maxCellFraction);
 853
 854     //Cut on time of clusters
 855     Double_t tof = clus->GetTOF()*1.e9;
 856     if(tof < fTimeCutMin || tof > fTimeCutMax)
 857     {
 858       if(GetDebug() > 0 )printf("AliAnaCalorimeterQA - Remove cluster with TOF %f\n",tof);
 859       continue;
 860     }
 861
 862     // Get cluster kinematics
 863     clus->GetMomentum(mom,v);
 864
 865     // Check only certain regions
 866     Bool_t in = kTRUE;
 867     if(IsFiducialCutOn()) in =  GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
 868     if(!in) continue;
 869
 870     // MC labels
 871     nLabel = clus->GetNLabels();
 872     labels = clus->GetLabels();
 873
 874     // SuperModule number of cluster
 875     nModule = GetModuleNumber(clus);
 876
 877     // Cells per cluster
 878     nCaloCellsPerCluster = clus->GetNCells();
 879
 880     // Cluster mathed with track?
 881     matched = GetCaloPID()->IsTrackMatched(clus,GetCaloUtils(), GetReader()->GetInputEvent());
 882
 883     // Get some time averages
 884     Double_t averTime[4] = {0.,0.,0.,0.};
 885     CalculateAverageTime(clus, cells, averTime);
 886
 887     //Get time of max cell
 888     Double_t tmax  = cells->GetCellTime(absIdMax);
 889     GetCaloUtils()->RecalibrateCellTime(tmax, fCalorimeter, absIdMax,GetReader()->GetInputEvent()->GetBunchCrossNumber());
 890     tmax*=1.e9;
 891
 892     // Fill histograms related to single cluster
 893
 894
 895     // Fill some histograms before applying the exotic cell / bad map cut
 896     fhNCellsPerClusterNoCut  ->Fill(clus->E(), nCaloCellsPerCluster);
 897     if(nModule >=0 && nModule < fNModules) fhNCellsPerClusterModNoCut[nModule]->Fill(clus->E(), nCaloCellsPerCluster);
 898
 899     fhClusterMaxCellDiffNoCut->Fill(clus->E(),maxCellFraction);
 900
 901     Float_t ampMax = cells->GetCellAmplitude(absIdMax);
 902     GetCaloUtils()->RecalibrateCellAmplitude(ampMax,fCalorimeter, absIdMax);
 903
 904     if(fStudyExotic) ExoticHistograms(absIdMax, ampMax, clus, cells);
 905
 906     //Check bad clusters if requested and rejection was not on
 907     Bool_t goodCluster = IsGoodCluster(absIdMax, cells);
 908
 909     Float_t eCrossFrac = 0;
 910     if(ampMax > 0.01) eCrossFrac = 1-GetECross(absIdMax,cells)/ampMax;
 911
 912     if(!goodCluster)
 913     {
 914       BadClusterHistograms(clus, caloClusters, cells, absIdMax,
 915                            maxCellFraction, eCrossFrac, tmax, averTime);
 916       continue;
 917     }
 918
 919     ClusterHistograms(clus, caloClusters, cells, absIdMax,
 920                       maxCellFraction, eCrossFrac, tmax, averTime);
 921
 922     nCaloClustersAccepted++;
 923     nModule = GetModuleNumber(clus);
 924     if(nModule >=0 && nModule < fNModules)
 925     {
 926       if     (fCalorimeter=="EMCAL" && mom.E() > 2*fEMCALCellAmpMin)  nClustersInModule[nModule]++;
 927       else if(fCalorimeter=="PHOS"  && mom.E() > 2*fPHOSCellAmpMin )  nClustersInModule[nModule]++;
 928     }
 929
 930     // Cluster weights
 931     if(fStudyWeight) WeightHistograms(clus, cells);
 932
 933     // Cells in cluster position
 934     if(fFillAllPosHisto) CellInClusterPositionHistograms(clus);
 935
 936     // Fill histograms related to single cluster, mc vs data
 937     Int_t  mcOK = kFALSE;
 938     Int_t  pdg  = -1;
 939     if(IsDataMC() && nLabel > 0 && labels)
 940       mcOK = ClusterMCHistograms(mom, matched, labels, nLabel, pdg);
 941
 942     // Matched clusters with tracks, also do some MC comparison, needs input from ClusterMCHistograms
 943     if( matched &&  fFillAllTMHisto)
 944       ClusterMatchedWithTrackHistograms(clus,mom,mcOK,pdg);
 945
 946     // Invariant mass
 947     // Try to reduce background with a mild shower shape cut and no more than 1 maxima
 948     // in cluster and remove low energy clusters
 949     if(fFillAllPi0Histo && nCaloClusters > 1 && nCaloCellsPerCluster > 1 &&
 950        GetCaloUtils()->GetNumberOfLocalMaxima(clus,cells) == 1 &&
 951        clus->GetM02() < 0.5 && clus->E() > 0.3)
 952       InvariantMassHistograms(iclus, mom, nModule, caloClusters,cells);
 953
 954   }//cluster loop
 955
 956   // Number of clusters histograms
 957   if(nCaloClustersAccepted > 0) fhNClusters->Fill(nCaloClustersAccepted);
 958
 959   // Number of clusters per module
 960   for(Int_t imod = 0; imod < fNModules; imod++ )
 961   {
 962     if(GetDebug() > 1)
 963       printf("AliAnaCalorimeterQA::ClusterLoopHistograms() - module %d calo %s clusters %d\n", imod, fCalorimeter.Data(), nClustersInModule[imod]);
 964     fhNClustersMod->Fill(nClustersInModule[imod],imod);
 965   }
 966
 967   delete [] nClustersInModule;
 968
 969 }
 970
 971 //______________________________________________________________________________________________________
 972 Bool_t AliAnaCalorimeterQA::ClusterMCHistograms(const TLorentzVector mom, const Bool_t matched,
 973                                                 const Int_t * labels, const Int_t nLabels, Int_t & pdg )
 974 {
 975
 976   //Fill histograms only possible when simulation
 977
 978   if(!labels || nLabels<=0)
 979   {
 980     if(GetDebug() > 1) printf("AliAnaCalorimeterQA::ClusterMCHistograms() - Strange, labels array %p, n labels %d \n", labels,nLabels);
 981     return kFALSE;
 982   }
 983
 984   if(GetDebug() > 1)
 985   {
 986     printf("AliAnaCalorimeterQA::ClusterMCHistograms() - Primaries: nlabels %d\n",nLabels);
 987   }
 988
 989   Float_t e   = mom.E();
 990   Float_t eta = mom.Eta();
 991   Float_t phi = mom.Phi();
 992   if(phi < 0) phi +=TMath::TwoPi();
 993
 994   AliAODMCParticle * aodprimary  = 0x0;
 995   TParticle * primary = 0x0;
 996
 997   //Play with the MC stack if available
 998   Int_t label = labels[0];
 999
1000   if(label < 0)
1001   {
1002     if(GetDebug() >= 0) printf("AliAnaCalorimeterQA::ClusterMCHistograms() *** bad label ***:  label %d \n", label);
1003     return kFALSE;
1004   }
1005
1006   Int_t pdg0  =-1;Int_t status = -1; Int_t iMother = -1; Int_t iParent = -1;
1007   Float_t vxMC= 0; Float_t vyMC = 0;
1008   Float_t eMC = 0; Float_t ptMC= 0; Float_t phiMC =0; Float_t etaMC = 0;
1009   Int_t charge = 0;
1010
1011   //Check the origin.
1012   Int_t tag = GetMCAnalysisUtils()->CheckOrigin(labels,nLabels, GetReader(),0);
1013
1014   if     ( GetReader()->ReadStack() &&
1015           !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCUnknown))
1016   { //it MC stack and known tag
1017
1018     if( label >= GetMCStack()->GetNtrack())
1019     {
1020       if(GetDebug() >= 0) printf("AliAnaCalorimeterQA::ClusterMCHistograms() *** large label ***:  label %d, n tracks %d \n", label, GetMCStack()->GetNtrack());
1021       return kFALSE;
1022     }
1023
1024     primary = GetMCStack()->Particle(label);
1025     iMother = label;
1026     pdg0    = TMath::Abs(primary->GetPdgCode());
1027     pdg     = pdg0;
1028     status  = primary->GetStatusCode();
1029     vxMC    = primary->Vx();
1030     vyMC    = primary->Vy();
1031     iParent = primary->GetFirstMother();
1032
1033     if(GetDebug() > 1 )
1034     {
1035       printf("AliAnaCalorimeterQA::ClusterMCHistograms() - Cluster most contributing mother: \n");
1036       printf("\t Mother label %d, pdg %d, %s, status %d, parent %d \n",iMother, pdg0, primary->GetName(),status, iParent);
1037     }
1038
1039     //Get final particle, no conversion products
1040     if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion))
1041     {
1042       //Get the parent
1043       primary = GetMCStack()->Particle(iParent);
1044       pdg = TMath::Abs(primary->GetPdgCode());
1045
1046       if(GetDebug() > 1 ) printf("AliAnaCalorimeterQA::ClusterMCHistograms() - Converted cluster!. Find before conversion: \n");
1047
1048       while((pdg == 22 || pdg == 11) && status != 1)
1049       {
1050         Int_t iMotherOrg = iMother;
1051         iMother = iParent;
1052         primary = GetMCStack()->Particle(iMother);
1053         status  = primary->GetStatusCode();
1054         pdg     = TMath::Abs(primary->GetPdgCode());
1055         iParent = primary->GetFirstMother();
1056
1057         // If gone too back and non stable, assign the decay photon/electron
1058         // there are other possible decays, ignore them for the moment
1059         if(pdg==111 || pdg==221)
1060         {
1061           primary = GetMCStack()->Particle(iMotherOrg);
1062           break;
1063         }
1064
1065         if( iParent < 0 )
1066         {
1067           iParent = iMother;
1068           printf("break\n");
1069           break;
1070         }
1071
1072         if(GetDebug() > 1 )printf("\t pdg %d, index %d, %s, status %d \n",pdg, iMother,  primary->GetName(),status);
1073       }
1074
1075       if(GetDebug() > 1 )
1076       {
1077         printf("AliAnaCalorimeterQA::ClusterHistograms() - Converted Cluster mother before conversion: \n");
1078         printf("\t Mother label %d, pdg %d, %s, status %d, parent %d \n",iMother, pdg, primary->GetName(), status, iParent);
1079       }
1080
1081     }
1082
1083     //Overlapped pi0 (or eta, there will be very few), get the meson
1084     if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) ||
1085        GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta))
1086     {
1087       if(GetDebug() > 1 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped Meson decay!, Find it: \n");
1088
1089       while(pdg != 111 && pdg != 221)
1090       {
1091         //printf("iMother %d, pdg %d, iParent %d, pdg %d\n",iMother,pdg,iParent,GetMCStack()->Particle(iParent)->GetPdgCode());
1092         iMother = iParent;
1093         primary = GetMCStack()->Particle(iMother);
1094         status  = primary->GetStatusCode();
1095         pdg     = TMath::Abs(primary->GetPdgCode());
1096         iParent = primary->GetFirstMother();
1097
1098         if( iParent < 0 )break;
1099
1100         if(GetDebug() > 1 ) printf("\t pdg %d, %s, index %d\n",pdg,  primary->GetName(),iMother);
1101
1102         if(iMother==-1)
1103         {
1104           printf("AliAnaCalorimeterQA::ClusterHistograms() - Tagged as Overlapped photon but meson not found, why?\n");
1105           //break;
1106         }
1107       }
1108
1109       if(GetDebug() > 2 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped %s decay, label %d \n",
1110                                  primary->GetName(),iMother);
1111     }
1112
1113     eMC    = primary->Energy();
1114     ptMC   = primary->Pt();
1115     phiMC  = primary->Phi();
1116     etaMC  = primary->Eta();
1117     pdg    = TMath::Abs(primary->GetPdgCode());
1118     charge = (Int_t) TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
1119
1120   }
1121   else if( GetReader()->ReadAODMCParticles() &&
1122           !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCUnknown))
1123   {//it MC AOD and known tag
1124     //Get the list of MC particles
1125     if(!GetReader()->GetAODMCParticles(0))
1126       AliFatal("MCParticles not available!");
1127
1128     aodprimary = (AliAODMCParticle*) (GetReader()->GetAODMCParticles(0))->At(label);
1129     iMother = label;
1130     pdg0    = TMath::Abs(aodprimary->GetPdgCode());
1131     pdg     = pdg0;
1132     status  = aodprimary->IsPrimary();
1133     vxMC    = aodprimary->Xv();
1134     vyMC    = aodprimary->Yv();
1135     iParent = aodprimary->GetMother();
1136
1137     if(GetDebug() > 1 )
1138     {
1139       printf("AliAnaCalorimeterQA::ClusterHistograms() - Cluster most contributing mother: \n");
1140       printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d, parent %d \n",
1141              iMother, pdg0, aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary(), iParent);
1142     }
1143
1144     //Get final particle, no conversion products
1145     if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion))
1146     {
1147       if(GetDebug() > 1 )
1148         printf("AliAnaCalorimeterQA::ClusterHistograms() - Converted cluster!. Find before conversion: \n");
1149       //Get the parent
1150       aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles(0))->At(iParent);
1151       pdg = TMath::Abs(aodprimary->GetPdgCode());
1152       while ((pdg == 22 || pdg == 11) && !aodprimary->IsPhysicalPrimary())
1153       {
1154         Int_t iMotherOrg = iMother;
1155         iMother    = iParent;
1156         aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles(0))->At(iMother);
1157         status     = aodprimary->IsPrimary();
1158         iParent    = aodprimary->GetMother();
1159         pdg        = TMath::Abs(aodprimary->GetPdgCode());
1160
1161         // If gone too back and non stable, assign the decay photon/electron
1162         // there are other possible decays, ignore them for the moment
1163         if(pdg==111 || pdg==221)
1164         {
1165           aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles(0))->At(iMotherOrg);
1166           break;
1167         }
1168
1169         if(iParent < 0 )
1170         {
1171           iParent = iMother;
1172           break;
1173         }
1174
1175         if(GetDebug() > 1 )
1176           printf("\t pdg %d, index %d, Primary? %d, Physical Primary? %d \n",
1177                  pdg, iMother, aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary());
1178       }
1179
1180       if(GetDebug() > 1 )
1181       {
1182         printf("AliAnaCalorimeterQA::ClusterHistograms() - Converted Cluster mother before conversion: \n");
1183         printf("\t Mother label %d, pdg %d, parent %d, Primary? %d, Physical Primary? %d \n",
1184                iMother, pdg, iParent, aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary());
1185       }
1186
1187     }
1188
1189     //Overlapped pi0 (or eta, there will be very few), get the meson
1190     if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) ||
1191        GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta))
1192     {
1193       if(GetDebug() > 1 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped Meson decay!, Find it: PDG %d, mom %d \n",pdg, iMother);
1194       while(pdg != 111 && pdg != 221)
1195       {
1196         iMother    = iParent;
1197         aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles(0))->At(iMother);
1198         status     = aodprimary->IsPrimary();
1199         iParent    = aodprimary->GetMother();
1200         pdg        = TMath::Abs(aodprimary->GetPdgCode());
1201
1202         if(iParent < 0 ) break;
1203
1204         if(GetDebug() > 1 ) printf("\t pdg %d, index %d\n",pdg, iMother);
1205
1206         if(iMother==-1)
1207         {
1208           printf("AliAnaCalorimeterQA::ClusterHistograms() - Tagged as Overlapped photon but meson not found, why?\n");
1209           //break;
1210         }
1211       }
1212
1213       if(GetDebug() > 2 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped %s decay, label %d \n",
1214                                  aodprimary->GetName(),iMother);
1215     }
1216
1217     status = aodprimary->IsPrimary();
1218     eMC    = aodprimary->E();
1219     ptMC   = aodprimary->Pt();
1220     phiMC  = aodprimary->Phi();
1221     etaMC  = aodprimary->Eta();
1222     pdg    = TMath::Abs(aodprimary->GetPdgCode());
1223     charge = aodprimary->Charge();
1224
1225   }
1226
1227   //Float_t vz = primary->Vz();
1228   Float_t rVMC = TMath::Sqrt(vxMC*vxMC + vyMC*vyMC);
1229   if((pdg == 22 || TMath::Abs(pdg)==11) && status!=1)
1230   {
1231     fhEMVxyz   ->Fill(vxMC,vyMC);//,vz);
1232     fhEMR      ->Fill(e,rVMC);
1233   }
1234
1235   //printf("reco e %f, pt %f, phi %f, eta %f \n", e, pt, phi, eta);
1236   //printf("prim e %f, pt %f, phi %f, eta %f \n", eMC,ptMC,phiMC ,etaMC );
1237   //printf("vertex: vx %f, vy %f, vz %f, r %f \n", vxMC, vyMC, vz, r);
1238
1239   //Overlapped pi0 (or eta, there will be very few)
1240   if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0))
1241   {
1242     fhRecoMCE  [kmcPi0][matched]     ->Fill(e,eMC);
1243     if(e > 0.5 && eMC > 0.5) fhRecoMCEta[kmcPi0][(matched)]->Fill(eta,etaMC);
1244     if(e > 0.5 && eMC > 0.5) fhRecoMCPhi[kmcPi0][(matched)]->Fill(phi,phiMC);
1245     if(eMC > 0) fhRecoMCRatioE  [kmcPi0][(matched)]->Fill(e,e/eMC);
1246     fhRecoMCDeltaE  [kmcPi0][(matched)]->Fill(e,eMC-e);
1247     fhRecoMCDeltaPhi[kmcPi0][(matched)]->Fill(e,phiMC-phi);
1248     fhRecoMCDeltaEta[kmcPi0][(matched)]->Fill(e,etaMC-eta);
1249   }//Overlapped pizero decay
1250   else     if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta))
1251   {
1252     fhRecoMCE  [kmcEta][(matched)]     ->Fill(e,eMC);
1253     if(e > 0.5 && eMC > 0.5) fhRecoMCEta[kmcEta][(matched)]->Fill(eta,etaMC);
1254     if(e > 0.5 && eMC > 0.5) fhRecoMCPhi[kmcEta][(matched)]->Fill(phi,phiMC);
1255     if(eMC > 0) fhRecoMCRatioE  [kmcEta][(matched)]->Fill(e,e/eMC);
1256     fhRecoMCDeltaE  [kmcEta][(matched)]->Fill(e,eMC-e);
1257     fhRecoMCDeltaPhi[kmcEta][(matched)]->Fill(e,phiMC-phi);
1258     fhRecoMCDeltaEta[kmcEta][(matched)]->Fill(e,etaMC-eta);
1259   }//Overlapped eta decay
1260   else if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton) &&
1261           !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion))
1262   {
1263     fhRecoMCE  [kmcPhoton][(matched)]     ->Fill(e,eMC);
1264     if(e > 0.5 && eMC > 0.5) fhRecoMCEta[kmcPhoton][(matched)]->Fill(eta,etaMC);
1265     if(e > 0.5 && eMC > 0.5) fhRecoMCPhi[kmcPhoton][(matched)]->Fill(phi,phiMC);
1266     if(eMC > 0) fhRecoMCRatioE  [kmcPhoton][(matched)]->Fill(e,e/eMC);
1267
1268     fhRecoMCDeltaE  [kmcPhoton][(matched)]->Fill(e,eMC-e);
1269     fhRecoMCDeltaPhi[kmcPhoton][(matched)]->Fill(e,phiMC-phi);
1270     fhRecoMCDeltaEta[kmcPhoton][(matched)]->Fill(e,etaMC-eta);
1271   }//photon
1272   else if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCElectron) &&
1273           !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion))
1274   {
1275     fhRecoMCE  [kmcElectron][(matched)]     ->Fill(e,eMC);
1276     if(e > 0.5 && eMC > 0.5) fhRecoMCEta[kmcElectron][(matched)]->Fill(eta,etaMC);
1277     if(e > 0.5 && eMC > 0.5) fhRecoMCPhi[kmcElectron][(matched)]->Fill(phi,phiMC);
1278     if(eMC > 0) fhRecoMCRatioE  [kmcElectron][(matched)]->Fill(e,e/eMC);
1279     fhRecoMCDeltaE  [kmcElectron][(matched)]->Fill(e,eMC-e);
1280     fhRecoMCDeltaPhi[kmcElectron][(matched)]->Fill(e,phiMC-phi);
1281     fhRecoMCDeltaEta[kmcElectron][(matched)]->Fill(e,etaMC-eta);
1282     fhEMVxyz   ->Fill(vxMC,vyMC);//,vz);
1283     fhEMR      ->Fill(e,rVMC);
1284   }
1285   else if(charge == 0)
1286   {
1287     fhRecoMCE  [kmcNeHadron][(matched)]     ->Fill(e,eMC);
1288     if(e > 0.5 && eMC > 0.5) fhRecoMCEta[kmcNeHadron][(matched)]->Fill(eta,etaMC);
1289     if(e > 0.5 && eMC > 0.5) fhRecoMCPhi[kmcNeHadron][(matched)]->Fill(phi,phiMC);
1290     if(eMC > 0) fhRecoMCRatioE  [kmcNeHadron][(matched)]->Fill(e,e/eMC);
1291     fhRecoMCDeltaE  [kmcNeHadron][(matched)]->Fill(e,eMC-e);
1292     fhRecoMCDeltaPhi[kmcNeHadron][(matched)]->Fill(e,phiMC-phi);
1293     fhRecoMCDeltaEta[kmcNeHadron][(matched)]->Fill(e,etaMC-eta);
1294     fhHaVxyz     ->Fill(vxMC,vyMC);//,vz);
1295     fhHaR        ->Fill(e,rVMC);
1296   }
1297   else if(charge!=0)
1298   {
1299     fhRecoMCE  [kmcChHadron][(matched)]     ->Fill(e,eMC);
1300     if(e > 0.5 && eMC > 0.5) fhRecoMCEta[kmcChHadron][(matched)]->Fill(eta,etaMC);
1301     if(e > 0.5 && eMC > 0.5) fhRecoMCPhi[kmcChHadron][(matched)]->Fill(phi,phiMC);
1302     if(eMC > 0) fhRecoMCRatioE  [kmcChHadron][(matched)]->Fill(e,e/eMC);
1303     fhRecoMCDeltaE  [kmcChHadron][(matched)]->Fill(e,eMC-e);
1304     fhRecoMCDeltaPhi[kmcChHadron][(matched)]->Fill(e,phiMC-phi);
1305     fhRecoMCDeltaEta[kmcChHadron][(matched)]->Fill(e,etaMC-eta);
1306     fhHaVxyz     ->Fill(vxMC,vyMC);//,vz);
1307     fhHaR        ->Fill(e,rVMC);
1308   }
1309
1310   if(primary || aodprimary) return kTRUE ;
1311   else                      return kFALSE;
1312
1313 }
1314
1315 //________________________________________________________________________________________________
1316 void AliAnaCalorimeterQA::ClusterMatchedWithTrackHistograms(AliVCluster *clus, TLorentzVector mom,
1317                                                             const Bool_t okPrimary, const Int_t pdg)
1318 {
1319   //Histograms for clusters matched with tracks
1320
1321   Float_t e   = mom.E();
1322   Float_t pt  = mom.Pt();
1323   Float_t eta = mom.Eta();
1324   Float_t phi = mom.Phi();
1325   if(phi < 0) phi +=TMath::TwoPi();
1326
1327   if(fFillAllTH12)
1328   {
1329     fhECharged      ->Fill(e);
1330     fhPtCharged     ->Fill(pt);
1331     fhPhiCharged    ->Fill(phi);
1332     fhEtaCharged    ->Fill(eta);
1333   }
1334
1335   //Study the track and matched cluster if track exists.
1336
1337   AliVTrack *track = GetCaloUtils()->GetMatchedTrack(clus, GetReader()->GetInputEvent());
1338
1339   if(!track) return ;
1340
1341   Double_t tpt   = track->Pt();
1342   Double_t tmom  = track->P();
1343   Double_t dedx  = track->GetTPCsignal();
1344   Int_t    nITS  = track->GetNcls(0);
1345   Int_t    nTPC  = track->GetNcls(1);
1346
1347   // Residuals
1348   Float_t deta  = clus->GetTrackDz();
1349   Float_t dphi  = clus->GetTrackDx();
1350   Double_t  dR  = TMath::Sqrt(dphi*dphi + deta*deta);
1351
1352   Double_t eOverP = e/tmom;
1353
1354   fh1EOverP->Fill(tpt, eOverP);
1355   if(dR < 0.02) fh1EOverPR02->Fill(tpt,eOverP);
1356
1357   fh2dR->Fill(e,dR);
1358   fh2MatchdEdx->Fill(tmom,dedx);
1359
1360   if(IsDataMC() && okPrimary)
1361   {
1362     Double_t  charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
1363
1364     if(TMath::Abs(pdg) == 11)
1365     {
1366       fhMCEle1EOverP->Fill(tpt,eOverP);
1367       fhMCEle1dR->Fill(dR);
1368       fhMCEle2MatchdEdx->Fill(tmom,dedx);
1369       if(dR < 0.02) fhMCEle1EOverPR02->Fill(tpt,eOverP);
1370     }
1371     else if(charge!=0)
1372     {
1373       fhMCChHad1EOverP->Fill(tpt,eOverP);
1374       fhMCChHad1dR->Fill(dR);
1375       fhMCChHad2MatchdEdx->Fill(tmom,dedx);
1376       if(dR < 0.02) fhMCChHad1EOverPR02->Fill(tpt,eOverP);
1377     }
1378     else if(charge == 0)
1379     {
1380       fhMCNeutral1EOverP->Fill(tpt,eOverP);
1381       fhMCNeutral1dR->Fill(dR);
1382       fhMCNeutral2MatchdEdx->Fill(tmom,dedx);
1383       if(dR < 0.02) fhMCNeutral1EOverPR02->Fill(tpt,eOverP);
1384     }
1385   }//DataMC
1386
1387   if(dR < 0.02 && eOverP > 0.6 && eOverP< 1.2
1388      && clus->GetNCells() > 1 && nITS > 3 && nTPC > 20)
1389   {
1390     fh2EledEdx->Fill(tmom,dedx);
1391   }
1392
1393 }
1394
1395 //___________________________________
1396 void AliAnaCalorimeterQA::Correlate()
1397 {
1398   // Correlate information from PHOS and EMCAL and with V0 and track multiplicity
1399
1400   //Clusters
1401   TObjArray * caloClustersEMCAL = GetEMCALClusters();
1402   TObjArray * caloClustersPHOS  = GetPHOSClusters();
1403
1404   Int_t nclEMCAL = caloClustersEMCAL->GetEntriesFast();
1405   Int_t nclPHOS  = caloClustersPHOS ->GetEntriesFast();
1406
1407   Float_t sumClusterEnergyEMCAL = 0;
1408   Float_t sumClusterEnergyPHOS  = 0;
1409   Int_t iclus = 0;
1410   for(iclus = 0 ; iclus <  caloClustersEMCAL->GetEntriesFast() ; iclus++)
1411     sumClusterEnergyEMCAL += ((AliVCluster*)caloClustersEMCAL->At(iclus))->E();
1412   for(iclus = 0 ; iclus <  caloClustersPHOS->GetEntriesFast(); iclus++)
1413     sumClusterEnergyPHOS += ((AliVCluster*)caloClustersPHOS->At(iclus))->E();
1414
1415
1416   //Cells
1417
1418   AliVCaloCells * cellsEMCAL = GetEMCALCells();
1419   AliVCaloCells * cellsPHOS  = GetPHOSCells();
1420
1421   Int_t ncellsEMCAL = cellsEMCAL->GetNumberOfCells();
1422   Int_t ncellsPHOS  = cellsPHOS ->GetNumberOfCells();
1423
1424   Float_t sumCellEnergyEMCAL = 0;
1425   Float_t sumCellEnergyPHOS  = 0;
1426   Int_t icell = 0;
1427   for(icell = 0 ; icell < cellsEMCAL->GetNumberOfCells()  ; icell++)
1428     sumCellEnergyEMCAL += cellsEMCAL->GetAmplitude(icell);
1429   for(icell = 0 ; icell <  cellsPHOS->GetNumberOfCells(); icell++)
1430     sumCellEnergyPHOS += cellsPHOS->GetAmplitude(icell);
1431
1432
1433   //Fill Histograms
1434   fhCaloCorrNClusters->Fill(nclEMCAL,nclPHOS);
1435   fhCaloCorrEClusters->Fill(sumClusterEnergyEMCAL,sumClusterEnergyPHOS);
1436   fhCaloCorrNCells   ->Fill(ncellsEMCAL,ncellsPHOS);
1437   fhCaloCorrECells   ->Fill(sumCellEnergyEMCAL,sumCellEnergyPHOS);
1438
1439   Int_t v0S = GetV0Signal(0)+GetV0Signal(1);
1440   Int_t v0M = GetV0Multiplicity(0)+GetV0Multiplicity(1);
1441   Int_t trM = GetTrackMultiplicity();
1442   if(fCalorimeter=="PHOS"){
1443     fhCaloV0MCorrNClusters   ->Fill(v0M,nclPHOS);
1444     fhCaloV0MCorrEClusters   ->Fill(v0M,sumClusterEnergyPHOS);
1445     fhCaloV0MCorrNCells      ->Fill(v0M,ncellsPHOS);
1446     fhCaloV0MCorrECells      ->Fill(v0M,sumCellEnergyPHOS);
1447
1448     fhCaloV0SCorrNClusters   ->Fill(v0S,nclPHOS);
1449     fhCaloV0SCorrEClusters   ->Fill(v0S,sumClusterEnergyPHOS);
1450     fhCaloV0SCorrNCells      ->Fill(v0S,ncellsPHOS);
1451     fhCaloV0SCorrECells      ->Fill(v0S,sumCellEnergyPHOS);
1452
1453     fhCaloTrackMCorrNClusters->Fill(trM,nclPHOS);
1454     fhCaloTrackMCorrEClusters->Fill(trM,sumClusterEnergyPHOS);
1455     fhCaloTrackMCorrNCells   ->Fill(trM,ncellsPHOS);
1456     fhCaloTrackMCorrECells   ->Fill(trM,sumCellEnergyPHOS);
1457   }
1458   else{
1459     fhCaloV0MCorrNClusters   ->Fill(v0M,nclEMCAL);
1460     fhCaloV0MCorrEClusters   ->Fill(v0M,sumClusterEnergyEMCAL);
1461     fhCaloV0MCorrNCells      ->Fill(v0M,ncellsEMCAL);
1462     fhCaloV0MCorrECells      ->Fill(v0M,sumCellEnergyEMCAL);
1463
1464     fhCaloV0SCorrNClusters   ->Fill(v0S,nclEMCAL);
1465     fhCaloV0SCorrEClusters   ->Fill(v0S,sumClusterEnergyEMCAL);
1466     fhCaloV0SCorrNCells      ->Fill(v0S,ncellsEMCAL);
1467     fhCaloV0SCorrECells      ->Fill(v0S,sumCellEnergyEMCAL);
1468
1469     fhCaloTrackMCorrNClusters->Fill(trM,nclEMCAL);
1470     fhCaloTrackMCorrEClusters->Fill(trM,sumClusterEnergyEMCAL);
1471     fhCaloTrackMCorrNCells   ->Fill(trM,ncellsEMCAL);
1472     fhCaloTrackMCorrECells   ->Fill(trM,sumCellEnergyEMCAL);
1473   }
1474
1475   if(GetDebug() > 0 )
1476   {
1477     printf("AliAnaCalorimeterQA::Correlate(): \n");
1478     printf("\t EMCAL: N cells %d, N clusters  %d, summed E cells %f, summed E clusters %f \n",
1479            ncellsEMCAL,nclEMCAL, sumCellEnergyEMCAL,sumClusterEnergyEMCAL);
1480     printf("\t PHOS : N cells %d, N clusters  %d, summed E cells %f, summed E clusters %f \n",
1481            ncellsPHOS,nclPHOS,sumCellEnergyPHOS,sumClusterEnergyPHOS);
1482     printf("\t V0 : Signal %d, Multiplicity  %d, Track Multiplicity %d \n", v0S,v0M,trM);
1483   }
1484
1485 }
1486
1487 //__________________________________________________
1488 TObjString * AliAnaCalorimeterQA::GetAnalysisCuts()
1489 {
1490   //Save parameters used for analysis
1491   TString parList ; //this will be list of parameters used for this analysis.
1492   const Int_t buffersize = 255;
1493   char onePar[buffersize] ;
1494
1495   snprintf(onePar,buffersize,"--- AliAnaCalorimeterQA ---\n") ;
1496   parList+=onePar ;
1497   snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
1498   parList+=onePar ;
1499   snprintf(onePar,buffersize,"Time Cut : %2.2f < T < %2.2f ns  \n",fTimeCutMin, fTimeCutMax) ;
1500   parList+=onePar ;
1501   snprintf(onePar,buffersize,"PHOS Cell Amplitude > %2.2f GeV, EMCAL Cell Amplitude > %2.2f GeV  \n",fPHOSCellAmpMin, fEMCALCellAmpMin) ;
1502   parList+=onePar ;
1503   //Get parameters set in base class.
1504   //parList += GetBaseParametersList() ;
1505
1506   //Get parameters set in FiducialCut class (not available yet)
1507   //parlist += GetFidCut()->GetFidCutParametersList()
1508
1509   return new TObjString(parList) ;
1510 }
1511
1512 //___________________________________________________________________________________
1513 void AliAnaCalorimeterQA::ExoticHistograms(const Int_t absIdMax, const Float_t ampMax,
1514                                            AliVCluster *clus, AliVCaloCells* cells)
1515 {
1516   // Calculate weights
1517
1518   if(ampMax < 0.01)
1519   {
1520     printf("AliAnaCalorimeterQA::ExoticHistograms()- Low amplitude energy %f\n",ampMax);
1521     return;
1522   }
1523
1524   Float_t l0 = clus->GetM02();
1525   Float_t en = clus->E();
1526   Int_t   nc = clus->GetNCells();
1527
1528   Double_t tmax  = cells->GetCellTime(absIdMax);
1529   GetCaloUtils()->RecalibrateCellTime(tmax, fCalorimeter, absIdMax,GetReader()->GetInputEvent()->GetBunchCrossNumber());
1530
1531   for(Int_t ie = 0; ie < fExoNECrossCuts; ie++)
1532   {
1533     for(Int_t idt = 0; idt < fExoNDTimeCuts; idt++)
1534     {
1535       Float_t eCrossFrac = 1-GetECross(absIdMax,cells, fExoDTimeCuts[idt])/ampMax;
1536
1537       if(eCrossFrac > fExoECrossCuts[ie])
1538       {
1539         //Exotic
1540         fhExoL0    [ie][idt]->Fill(en,l0  );
1541         fhExoTime  [ie][idt]->Fill(en,tmax);
1542
1543         // Diff time, do for one cut in e cross
1544         if(ie == 0)
1545         {
1546           for (Int_t icell = 0; icell < clus->GetNCells(); icell++)
1547           {
1548             Int_t absId  = clus->GetCellsAbsId()[icell];
1549             Double_t time  = cells->GetCellTime(absId);
1550             GetCaloUtils()->RecalibrateCellTime(time, fCalorimeter, absId,GetReader()->GetInputEvent()->GetBunchCrossNumber());
1551
1552             Float_t diff = (tmax-time)*1e9;
1553             fhExoDTime[idt]->Fill(en, diff);
1554           }
1555         }
1556       }
1557       else
1558       {
1559         fhExoECross[ie][idt]->Fill(en,eCrossFrac);
1560         fhExoNCell[ie][idt]->Fill(en,nc);
1561       }
1562     } // D time cut loop
1563   } // e cross cut loop
1564 }
1565
1566 //____________________________________________________
1567 TList * AliAnaCalorimeterQA::GetCreateOutputObjects()
1568 {
1569   // Create histograms to be saved in output file and
1570   // store them in outputContainer
1571
1572   TList * outputContainer = new TList() ;
1573   outputContainer->SetName("QAHistos") ;
1574
1575   //Histograms
1576   Int_t nptbins     = GetHistogramRanges()->GetHistoPtBins();           Float_t ptmax     = GetHistogramRanges()->GetHistoPtMax();           Float_t ptmin     = GetHistogramRanges()->GetHistoPtMin();
1577   Int_t nfineptbins = GetHistogramRanges()->GetHistoFinePtBins();           Float_t ptfinemax = GetHistogramRanges()->GetHistoFinePtMax();       Float_t ptfinemin = GetHistogramRanges()->GetHistoFinePtMin();
1578   Int_t nphibins    = GetHistogramRanges()->GetHistoPhiBins();              Float_t phimax    = GetHistogramRanges()->GetHistoPhiMax();          Float_t phimin    = GetHistogramRanges()->GetHistoPhiMin();
1579   Int_t netabins    = GetHistogramRanges()->GetHistoEtaBins();          Float_t etamax    = GetHistogramRanges()->GetHistoEtaMax();          Float_t etamin    = GetHistogramRanges()->GetHistoEtaMin();
1580   Int_t nmassbins   = GetHistogramRanges()->GetHistoMassBins();         Float_t massmax   = GetHistogramRanges()->GetHistoMassMax();           Float_t massmin   = GetHistogramRanges()->GetHistoMassMin();
1581   Int_t nasymbins   = GetHistogramRanges()->GetHistoAsymmetryBins();    Float_t asymmax   = GetHistogramRanges()->GetHistoAsymmetryMax();    Float_t asymmin   = GetHistogramRanges()->GetHistoAsymmetryMin();
1582   Int_t nPoverEbins = GetHistogramRanges()->GetHistoPOverEBins();       Float_t EOverPmax = GetHistogramRanges()->GetHistoPOverEMax();       Float_t EOverPmin = GetHistogramRanges()->GetHistoPOverEMin();
1583   Int_t ndedxbins   = GetHistogramRanges()->GetHistodEdxBins();         Float_t dedxmax   = GetHistogramRanges()->GetHistodEdxMax();         Float_t dedxmin   = GetHistogramRanges()->GetHistodEdxMin();
1584   Int_t ndRbins     = GetHistogramRanges()->GetHistodRBins();           Float_t dRmax     = GetHistogramRanges()->GetHistodRMax();           Float_t dRmin     = GetHistogramRanges()->GetHistodRMin();
1585   Int_t ntimebins   = GetHistogramRanges()->GetHistoTimeBins();         Float_t timemax   = GetHistogramRanges()->GetHistoTimeMax();         Float_t timemin   = GetHistogramRanges()->GetHistoTimeMin();
1586   Int_t nclbins     = GetHistogramRanges()->GetHistoNClustersBins();    Int_t   nclmax    = GetHistogramRanges()->GetHistoNClustersMax();    Int_t   nclmin    = GetHistogramRanges()->GetHistoNClustersMin();
1587   Int_t ncebins     = GetHistogramRanges()->GetHistoNCellsBins();       Int_t   ncemax    = GetHistogramRanges()->GetHistoNCellsMax();       Int_t   ncemin    = GetHistogramRanges()->GetHistoNCellsMin();
1588   Int_t nceclbins   = GetHistogramRanges()->GetHistoNClusterCellBins(); Int_t   nceclmax  = GetHistogramRanges()->GetHistoNClusterCellMax(); Int_t   nceclmin  = GetHistogramRanges()->GetHistoNClusterCellMin();
1589   Int_t nvdistbins  = GetHistogramRanges()->GetHistoVertexDistBins();   Float_t vdistmax  = GetHistogramRanges()->GetHistoVertexDistMax();   Float_t vdistmin  = GetHistogramRanges()->GetHistoVertexDistMin();
1590   Int_t rbins       = GetHistogramRanges()->GetHistoRBins();            Float_t rmax      = GetHistogramRanges()->GetHistoRMax();            Float_t rmin      = GetHistogramRanges()->GetHistoRMin();
1591   Int_t xbins       = GetHistogramRanges()->GetHistoXBins();            Float_t xmax      = GetHistogramRanges()->GetHistoXMax();            Float_t xmin      = GetHistogramRanges()->GetHistoXMin();
1592   Int_t ybins       = GetHistogramRanges()->GetHistoYBins();            Float_t ymax      = GetHistogramRanges()->GetHistoYMax();            Float_t ymin      = GetHistogramRanges()->GetHistoYMin();
1593   Int_t zbins       = GetHistogramRanges()->GetHistoZBins();            Float_t zmax      = GetHistogramRanges()->GetHistoZMax();            Float_t zmin      = GetHistogramRanges()->GetHistoZMin();
1594   Int_t ssbins      = GetHistogramRanges()->GetHistoShowerShapeBins();  Float_t ssmax     = GetHistogramRanges()->GetHistoShowerShapeMax();  Float_t ssmin     = GetHistogramRanges()->GetHistoShowerShapeMin();
1595   Int_t tdbins      = GetHistogramRanges()->GetHistoDiffTimeBins() ;    Float_t tdmax     = GetHistogramRanges()->GetHistoDiffTimeMax();     Float_t tdmin     = GetHistogramRanges()->GetHistoDiffTimeMin();
1596
1597   Int_t nv0sbins    = GetHistogramRanges()->GetHistoV0SignalBins();          Int_t nv0smax = GetHistogramRanges()->GetHistoV0SignalMax();          Int_t nv0smin = GetHistogramRanges()->GetHistoV0SignalMin();
1598   Int_t nv0mbins    = GetHistogramRanges()->GetHistoV0MultiplicityBins();    Int_t nv0mmax = GetHistogramRanges()->GetHistoV0MultiplicityMax();    Int_t nv0mmin = GetHistogramRanges()->GetHistoV0MultiplicityMin();
1599   Int_t ntrmbins    = GetHistogramRanges()->GetHistoTrackMultiplicityBins(); Int_t ntrmmax = GetHistogramRanges()->GetHistoTrackMultiplicityMax(); Int_t ntrmmin = GetHistogramRanges()->GetHistoTrackMultiplicityMin();
1600
1601   //EMCAL
1602   fNMaxCols = 48;
1603   fNMaxRows = 24;
1604   fNRCU     = 2 ;
1605   //PHOS
1606   if(fCalorimeter=="PHOS")
1607   {
1608     fNMaxCols = 56;
1609     fNMaxRows = 64;
1610     fNRCU     = 4 ;
1611   }
1612
1613   fhE  = new TH1F ("hE","E reconstructed clusters ", nptbins*5,ptmin,ptmax*5);
1614   fhE->SetXTitle("E (GeV)");
1615   outputContainer->Add(fhE);
1616
1617   if(fFillAllTH12)
1618   {
1619     fhPt  = new TH1F ("hPt","p_{T} reconstructed clusters", nptbins,ptmin,ptmax);
1620     fhPt->SetXTitle("p_{T} (GeV/c)");
1621     outputContainer->Add(fhPt);
1622
1623     fhPhi  = new TH1F ("hPhi","#phi reconstructed clusters ",nphibins,phimin,phimax);
1624     fhPhi->SetXTitle("#phi (rad)");
1625     outputContainer->Add(fhPhi);
1626
1627     fhEta  = new TH1F ("hEta","#eta reconstructed clusters ",netabins,etamin,etamax);
1628     fhEta->SetXTitle("#eta ");
1629     outputContainer->Add(fhEta);
1630   }
1631
1632   if(fFillAllTH3)
1633   {
1634     fhEtaPhiE  = new TH3F ("hEtaPhiE","#eta vs #phi vs energy, reconstructed clusters",
1635                            netabins,etamin,etamax,nphibins,phimin,phimax,nptbins,ptmin,ptmax);
1636     fhEtaPhiE->SetXTitle("#eta ");
1637     fhEtaPhiE->SetYTitle("#phi (rad)");
1638     fhEtaPhiE->SetZTitle("E (GeV) ");
1639     outputContainer->Add(fhEtaPhiE);
1640   }
1641
1642   fhClusterTimeEnergy  = new TH2F ("hClusterTimeEnergy","energy vs TOF, reconstructed clusters",
1643                                    nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
1644   fhClusterTimeEnergy->SetXTitle("E (GeV) ");
1645   fhClusterTimeEnergy->SetYTitle("TOF (ns)");
1646   outputContainer->Add(fhClusterTimeEnergy);
1647
1648   fhClusterPairDiffTimeE = new TH2F("hClusterPairDiffTimeE","cluster pair time difference vs E, only good clusters",
1649                                     nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1650   fhClusterPairDiffTimeE->SetXTitle("E_{cluster} (GeV)");
1651   fhClusterPairDiffTimeE->SetYTitle("#Delta t (ns)");
1652   outputContainer->Add(fhClusterPairDiffTimeE);
1653
1654   fhLambda0  = new TH2F ("hLambda0","shower shape, #lambda^{2}_{0} vs E",
1655                          nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1656   fhLambda0->SetXTitle("E_{cluster}");
1657   fhLambda0->SetYTitle("#lambda^{2}_{0}");
1658   outputContainer->Add(fhLambda0);
1659
1660   fhLambda1  = new TH2F ("hLambda1","shower shape, #lambda^{2}_{1} vs E for bad cluster ",
1661                          nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1662   fhLambda1->SetXTitle("E_{cluster}");
1663   fhLambda1->SetYTitle("#lambda^{2}_{1}");
1664   outputContainer->Add(fhLambda1);
1665
1666   fhDispersion  = new TH2F ("hDispersion","shower shape, Dispersion^{2} vs E for bad cluster ",
1667                             nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1668   fhDispersion->SetXTitle("E_{cluster}");
1669   fhDispersion->SetYTitle("Dispersion");
1670   outputContainer->Add(fhDispersion);
1671
1672   fhClusterMaxCellCloseCellRatio  = new TH2F ("hClusterMaxCellCloseCellRatio","energy vs ratio of max cell / neighbour cell, reconstructed clusters",
1673                                               nptbins,ptmin,ptmax, 100,0,1.);
1674   fhClusterMaxCellCloseCellRatio->SetXTitle("E_{cluster} (GeV) ");
1675   fhClusterMaxCellCloseCellRatio->SetYTitle("E_{cell i}/E_{cell max}");
1676   outputContainer->Add(fhClusterMaxCellCloseCellRatio);
1677
1678   fhClusterMaxCellCloseCellDiff  = new TH2F ("hClusterMaxCellCloseCellDiff","energy vs ratio of max cell / neighbour cell, reconstructed clusters",
1679                                              nptbins,ptmin,ptmax, 500,0,100.);
1680   fhClusterMaxCellCloseCellDiff->SetXTitle("E_{cluster} (GeV) ");
1681   fhClusterMaxCellCloseCellDiff->SetYTitle("E_{cell max}-E_{cell i} (GeV)");
1682   outputContainer->Add(fhClusterMaxCellCloseCellDiff);
1683
1684   fhClusterMaxCellDiff  = new TH2F ("hClusterMaxCellDiff","energy vs difference of cluster energy - max cell energy / cluster energy, good clusters",
1685                                     nptbins,ptmin,ptmax, 500,0,1.);
1686   fhClusterMaxCellDiff->SetXTitle("E_{cluster} (GeV) ");
1687   fhClusterMaxCellDiff->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
1688   outputContainer->Add(fhClusterMaxCellDiff);
1689
1690   fhClusterMaxCellDiffNoCut  = new TH2F ("hClusterMaxCellDiffNoCut","energy vs difference of cluster energy - max cell energy / cluster energy",
1691                                          nptbins,ptmin,ptmax, 500,0,1.);
1692   fhClusterMaxCellDiffNoCut->SetXTitle("E_{cluster} (GeV) ");
1693   fhClusterMaxCellDiffNoCut->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
1694   outputContainer->Add(fhClusterMaxCellDiffNoCut);
1695
1696   fhClusterMaxCellECross  = new TH2F ("hClusterMaxCellECross","1 - Energy in cross around max energy cell / max energy cell vs cluster energy, good clusters",
1697                                        nptbins,ptmin,ptmax, 400,-1,1.);
1698   fhClusterMaxCellECross->SetXTitle("E_{cluster} (GeV) ");
1699   fhClusterMaxCellECross->SetYTitle("1- E_{cross}/E_{cell max}");
1700   outputContainer->Add(fhClusterMaxCellECross);
1701
1702   fhNCellsPerClusterNoCut  = new TH2F ("hNCellsPerClusterNoCut","# cells per cluster vs energy, no bad clusters cut",
1703                                        nptbins,ptmin,ptmax, nceclbins,nceclmin,nceclmax);
1704   fhNCellsPerClusterNoCut->SetXTitle("E (GeV)");
1705   fhNCellsPerClusterNoCut->SetYTitle("n cells");
1706   outputContainer->Add(fhNCellsPerClusterNoCut);
1707
1708   fhNCellsPerCluster  = new TH2F ("hNCellsPerCluster","# cells per cluster vs energy",nptbins,ptmin,ptmax, nceclbins,nceclmin,nceclmax);
1709   fhNCellsPerCluster->SetXTitle("E (GeV)");
1710   fhNCellsPerCluster->SetYTitle("n cells");
1711   outputContainer->Add(fhNCellsPerCluster);
1712
1713   fhNClusters  = new TH1F ("hNClusters","# clusters", nclbins,nclmin,nclmax);
1714   fhNClusters->SetXTitle("number of clusters");
1715   outputContainer->Add(fhNClusters);
1716
1717   if(fStudyBadClusters)
1718   {
1719     fhBadClusterEnergy  = new TH1F ("hBadClusterEnergy","Bad cluster energy", nptbins,ptmin,ptmax);
1720     fhBadClusterEnergy->SetXTitle("E_{cluster} (GeV) ");
1721     outputContainer->Add(fhBadClusterEnergy);
1722
1723     fhBadClusterMaxCellCloseCellRatio  = new TH2F ("hBadClusterMaxCellCloseCellRatio","energy vs ratio of max cell / neighbour cell constributing cell, reconstructed bad clusters",
1724                                                    nptbins,ptmin,ptmax, 100,0,1.);
1725     fhBadClusterMaxCellCloseCellRatio->SetXTitle("E_{cluster} (GeV) ");
1726     fhBadClusterMaxCellCloseCellRatio->SetYTitle("ratio");
1727     outputContainer->Add(fhBadClusterMaxCellCloseCellRatio);
1728
1729     fhBadClusterMaxCellCloseCellDiff  = new TH2F ("hBadClusterMaxCellCloseCellDiff","energy vs ratio of max cell - neighbour cell constributing cell, reconstructed bad clusters",
1730                                                   nptbins,ptmin,ptmax, 500,0,100);
1731     fhBadClusterMaxCellCloseCellDiff->SetXTitle("E_{cluster} (GeV) ");
1732     fhBadClusterMaxCellCloseCellDiff->SetYTitle("E_{cell max} - E_{cell i} (GeV)");
1733     outputContainer->Add(fhBadClusterMaxCellCloseCellDiff);
1734
1735     fhBadClusterMaxCellDiff  = new TH2F ("hBadClusterMaxCellDiff","energy vs difference of cluster energy - max cell energy / cluster energy for bad clusters",
1736                                          nptbins,ptmin,ptmax, 500,0,1.);
1737     fhBadClusterMaxCellDiff->SetXTitle("E_{cluster} (GeV) ");
1738     fhBadClusterMaxCellDiff->SetYTitle("(E_{cluster} - E_{cell max}) / E_{cluster}");
1739     outputContainer->Add(fhBadClusterMaxCellDiff);
1740
1741     fhBadClusterTimeEnergy  = new TH2F ("hBadClusterTimeEnergy","energy vs TOF of reconstructed bad clusters",
1742                                         nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
1743     fhBadClusterTimeEnergy->SetXTitle("E_{cluster} (GeV) ");
1744     fhBadClusterTimeEnergy->SetYTitle("TOF (ns)");
1745     outputContainer->Add(fhBadClusterTimeEnergy);
1746
1747     fhBadClusterPairDiffTimeE = new TH2F("hBadClusterPairDiffTimeE","cluster pair time difference (bad - good) vs E from bad cluster",nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1748     fhBadClusterPairDiffTimeE->SetXTitle("E_{bad cluster} (GeV)");
1749     fhBadClusterPairDiffTimeE->SetYTitle("#Delta t (ns)");
1750     outputContainer->Add(fhBadClusterPairDiffTimeE);
1751
1752     fhBadClusterMaxCellECross  = new TH2F ("hBadClusterMaxCellECross","1 - Energy in cross around max energy cell / max energy cell vs cluster energy, bad clusters",
1753                                         nptbins,ptmin,ptmax, 400,-1,1.);
1754     fhBadClusterMaxCellECross->SetXTitle("E_{cluster} (GeV) ");
1755     fhBadClusterMaxCellECross->SetYTitle("1- E_{cross}/E_{cell max}");
1756     outputContainer->Add(fhBadClusterMaxCellECross);
1757
1758     if(fFillAllCellTimeHisto)
1759     {
1760       fhBadCellTimeSpreadRespectToCellMax = new TH2F ("hBadCellTimeSpreadRespectToCellMax","t_{cell max}-t_{cell i} from bad cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1761       fhBadCellTimeSpreadRespectToCellMax->SetXTitle("E (GeV)");
1762       fhBadCellTimeSpreadRespectToCellMax->SetYTitle("#Delta t_{cell max - i} (ns)");
1763       outputContainer->Add(fhBadCellTimeSpreadRespectToCellMax);
1764
1765       fhBadClusterMaxCellDiffAverageTime = new TH2F ("hBadClusterMaxCellDiffAverageTime","t_{cell max}-t_{average} from bad cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1766       fhBadClusterMaxCellDiffAverageTime->SetXTitle("E (GeV)");
1767       fhBadClusterMaxCellDiffAverageTime->SetYTitle("#Delta t_{cell max - average} (ns)");
1768       outputContainer->Add(fhBadClusterMaxCellDiffAverageTime);
1769
1770       fhBadClusterMaxCellDiffWeightedTime = new TH2F ("hBadClusterMaxCellDiffWeightedTime","t_{cell max}-t_{weighted} from bad cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1771       fhBadClusterMaxCellDiffWeightedTime->SetXTitle("E (GeV)");
1772       fhBadClusterMaxCellDiffWeightedTime->SetYTitle("#Delta t_{cell max - weighted} (ns)");
1773       outputContainer->Add(fhBadClusterMaxCellDiffWeightedTime);
1774
1775     }
1776
1777   }
1778
1779   if(fStudyExotic)
1780   {
1781     for(Int_t ie = 0; ie <fExoNECrossCuts; ie++)
1782     {
1783
1784       fhExoDTime[ie]  = new TH2F(Form("hExoDTime_ECross%d",ie),
1785                                  Form("#Delta time = t_{max}-t_{cells} vs E_{cluster} for exotic, 1-E_{+}/E_{max} < %2.2f",fExoECrossCuts[ie]),
1786                                  nptbins,ptmin,ptmax,tdbins,tdmin,tdmax);
1787       fhExoDTime[ie] ->SetYTitle("#Delta t (ns)");
1788       fhExoDTime[ie] ->SetXTitle("E (GeV)");
1789       outputContainer->Add(fhExoDTime[ie]) ;
1790
1791       for(Int_t idt = 0; idt < fExoNDTimeCuts; idt++)
1792       {
1793         fhExoNCell[ie][idt]  = new TH2F(Form("hExoNCell_ECross%d_DT%d",ie,idt),
1794                                      Form("N cells per cluster vs E cluster, 1-E_{+}/E_{max} < %2.2f, #Delta t < %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
1795                                      nptbins,ptmin,ptmax,nceclbins,nceclmin,nceclmax);
1796         fhExoNCell[ie][idt] ->SetYTitle("N cells");
1797         fhExoNCell[ie][idt] ->SetXTitle("E (GeV)");
1798         outputContainer->Add(fhExoNCell[ie][idt]) ;
1799
1800         fhExoL0   [ie][idt]  = new TH2F(Form("hExoL0_ECross%d_DT%d",ie,idt),
1801                                      Form("M02 vs E cluster for exotic, 1-E_{+}/E_{max} < %2.2f, #Delta t = %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
1802                                      nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1803         fhExoL0   [ie][idt] ->SetYTitle("#lambda^{2}_{0}");
1804         fhExoL0   [ie][idt] ->SetXTitle("E (GeV)");
1805         outputContainer->Add(fhExoL0[ie][idt]) ;
1806
1807         fhExoECross[ie][idt]  = new TH2F(Form("hExoECross_ECross%d_DT%d",ie,idt),
1808                                       Form("E cross for cells vs E cell, 1-E_{+}/E_{max} < %2.2f, #Delta t < %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
1809                                       nptbins,ptmin,ptmax,400,0,1);
1810         fhExoECross[ie][idt] ->SetYTitle("1-E_{+}/E_{cell max}");
1811         fhExoECross[ie][idt] ->SetXTitle("E_{cell} (GeV)");
1812         outputContainer->Add(fhExoECross[ie][idt]) ;
1813
1814         fhExoTime  [ie][idt]  = new TH2F(Form("hExoTime_ECross%d_DT%d",ie,idt),
1815                                         Form("Time of cluster (max cell) vs E cluster for exotic, 1-E_{+}/E_{max} < %2.2f, #Delta t = %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
1816                                         nptbins,ptmin,ptmax,ntimebins,timemin,timemax);
1817         fhExoTime  [ie][idt] ->SetYTitle("time_{max} (ns)");
1818         fhExoTime  [ie][idt] ->SetXTitle("E (GeV)");
1819         outputContainer->Add(fhExoTime[ie][idt]) ;
1820
1821       }
1822     }
1823   }
1824
1825   // Cluster size in terms of cells
1826   if(fStudyClustersAsymmetry)
1827   {
1828     fhDeltaIEtaDeltaIPhiE0[0]  = new TH2F ("hDeltaIEtaDeltaIPhiE0"," Cluster size in columns vs rows for E < 2 GeV, n cells > 3",
1829                                            50,0,50,50,0,50);
1830     fhDeltaIEtaDeltaIPhiE0[0]->SetXTitle("#Delta Column");
1831     fhDeltaIEtaDeltaIPhiE0[0]->SetYTitle("#Delta Row");
1832     outputContainer->Add(fhDeltaIEtaDeltaIPhiE0[0]);
1833
1834     fhDeltaIEtaDeltaIPhiE2[0]  = new TH2F ("hDeltaIEtaDeltaIPhiE2"," Cluster size in columns vs rows for 2 <E < 6 GeV, n cells > 3",
1835                                            50,0,50,50,0,50);
1836     fhDeltaIEtaDeltaIPhiE2[0]->SetXTitle("#Delta Column");
1837     fhDeltaIEtaDeltaIPhiE2[0]->SetYTitle("#Delta Row");
1838     outputContainer->Add(fhDeltaIEtaDeltaIPhiE2[0]);
1839
1840     fhDeltaIEtaDeltaIPhiE6[0]  = new TH2F ("hDeltaIEtaDeltaIPhiE6"," Cluster size in columns vs rows for E > 6 GeV, n cells > 3",
1841                                            50,0,50,50,0,50);
1842     fhDeltaIEtaDeltaIPhiE6[0]->SetXTitle("#Delta Column");
1843     fhDeltaIEtaDeltaIPhiE6[0]->SetYTitle("#Delta Row");
1844     outputContainer->Add(fhDeltaIEtaDeltaIPhiE6[0]);
1845
1846     fhDeltaIA[0]  = new TH2F ("hDeltaIA"," Cluster *asymmetry* in cell units vs E",
1847                               nptbins,ptmin,ptmax,21,-1.05,1.05);
1848     fhDeltaIA[0]->SetXTitle("E_{cluster}");
1849     fhDeltaIA[0]->SetYTitle("A_{cell in cluster}");
1850     outputContainer->Add(fhDeltaIA[0]);
1851
1852     fhDeltaIAL0[0]  = new TH2F ("hDeltaIAL0"," Cluster *asymmetry* in cell units vs #lambda^{2}_{0}",
1853                                 ssbins,ssmin,ssmax,21,-1.05,1.05);
1854     fhDeltaIAL0[0]->SetXTitle("#lambda^{2}_{0}");
1855     fhDeltaIAL0[0]->SetYTitle("A_{cell in cluster}");
1856     outputContainer->Add(fhDeltaIAL0[0]);
1857
1858     fhDeltaIAL1[0]  = new TH2F ("hDeltaIAL1"," Cluster *asymmetry* in cell units vs #lambda^{2}_{1}",
1859                                 ssbins,ssmin,ssmax,21,-1.05,1.05);
1860     fhDeltaIAL1[0]->SetXTitle("#lambda^{2}_{1}");
1861     fhDeltaIAL1[0]->SetYTitle("A_{cell in cluster}");
1862     outputContainer->Add(fhDeltaIAL1[0]);
1863
1864     fhDeltaIANCells[0]  = new TH2F ("hDeltaIANCells"," Cluster *asymmetry* in cell units vs N cells in cluster",
1865                                     nceclbins,nceclmin,nceclmax,21,-1.05,1.05);
1866     fhDeltaIANCells[0]->SetXTitle("N_{cell in cluster}");
1867     fhDeltaIANCells[0]->SetYTitle("A_{cell in cluster}");
1868     outputContainer->Add(fhDeltaIANCells[0]);
1869
1870
1871     fhDeltaIEtaDeltaIPhiE0[1]  = new TH2F ("hDeltaIEtaDeltaIPhiE0Charged"," Cluster size in columns vs rows for E < 2 GeV, n cells > 3, matched with track",
1872                                            50,0,50,50,0,50);
1873     fhDeltaIEtaDeltaIPhiE0[1]->SetXTitle("#Delta Column");
1874     fhDeltaIEtaDeltaIPhiE0[1]->SetYTitle("#Delta Row");
1875     outputContainer->Add(fhDeltaIEtaDeltaIPhiE0[1]);
1876
1877     fhDeltaIEtaDeltaIPhiE2[1]  = new TH2F ("hDeltaIEtaDeltaIPhiE2Charged"," Cluster size in columns vs rows for 2 <E < 6 GeV, n cells > 3, matched with track",
1878                                            50,0,50,50,0,50);
1879     fhDeltaIEtaDeltaIPhiE2[1]->SetXTitle("#Delta Column");
1880     fhDeltaIEtaDeltaIPhiE2[1]->SetYTitle("#Delta Row");
1881     outputContainer->Add(fhDeltaIEtaDeltaIPhiE2[1]);
1882
1883     fhDeltaIEtaDeltaIPhiE6[1]  = new TH2F ("hDeltaIEtaDeltaIPhiE6Charged"," Cluster size in columns vs rows for E > 6 GeV, n cells > 3, matched with track",
1884                                            50,0,50,50,0,50);
1885     fhDeltaIEtaDeltaIPhiE6[1]->SetXTitle("#Delta Column");
1886     fhDeltaIEtaDeltaIPhiE6[1]->SetYTitle("#Delta Row");
1887     outputContainer->Add(fhDeltaIEtaDeltaIPhiE6[1]);
1888
1889     fhDeltaIA[1]  = new TH2F ("hDeltaIACharged"," Cluster *asymmetry* in cell units vs E, matched with track",
1890                               nptbins,ptmin,ptmax,21,-1.05,1.05);
1891     fhDeltaIA[1]->SetXTitle("E_{cluster}");
1892     fhDeltaIA[1]->SetYTitle("A_{cell in cluster}");
1893     outputContainer->Add(fhDeltaIA[1]);
1894
1895     fhDeltaIAL0[1]  = new TH2F ("hDeltaIAL0Charged"," Cluster *asymmetry* in cell units vs #lambda^{2}_{0}, matched with track",
1896                                 ssbins,ssmin,ssmax,21,-1.05,1.05);
1897     fhDeltaIAL0[1]->SetXTitle("#lambda^{2}_{0}");
1898     fhDeltaIAL0[1]->SetYTitle("A_{cell in cluster}");
1899     outputContainer->Add(fhDeltaIAL0[1]);
1900
1901     fhDeltaIAL1[1]  = new TH2F ("hDeltaIAL1Charged"," Cluster *asymmetry* in cell units vs #lambda^{2}_{1}, matched with track",
1902                                 ssbins,ssmin,ssmax,21,-1.05,1.05);
1903     fhDeltaIAL1[1]->SetXTitle("#lambda^{2}_{1}");
1904     fhDeltaIAL1[1]->SetYTitle("A_{cell in cluster}");
1905     outputContainer->Add(fhDeltaIAL1[1]);
1906
1907     fhDeltaIANCells[1]  = new TH2F ("hDeltaIANCellsCharged"," Cluster *asymmetry* in cell units vs N cells in cluster, matched with track",
1908                                     nceclbins,nceclmin,nceclmax,21,-1.05,1.05);
1909     fhDeltaIANCells[1]->SetXTitle("N_{cell in cluster}");
1910     fhDeltaIANCells[1]->SetYTitle("A_{cell in cluster}");
1911     outputContainer->Add(fhDeltaIANCells[1]);
1912
1913     if(IsDataMC()){
1914       TString particle[]={"Photon","Electron","Conversion","Hadron"};
1915       for (Int_t iPart = 0; iPart < 4; iPart++) {
1916
1917         fhDeltaIAMC[iPart]  = new TH2F (Form("hDeltaIA_MC%s",particle[iPart].Data()),Form(" Cluster *asymmetry* in cell units vs E, from %s",particle[iPart].Data()),
1918                                         nptbins,ptmin,ptmax,21,-1.05,1.05);
1919         fhDeltaIAMC[iPart]->SetXTitle("E_{cluster}");
1920         fhDeltaIAMC[iPart]->SetYTitle("A_{cell in cluster}");
1921         outputContainer->Add(fhDeltaIAMC[iPart]);
1922       }
1923     }
1924
1925     if(fStudyBadClusters)
1926     {
1927       fhBadClusterDeltaIEtaDeltaIPhiE0  = new TH2F ("hBadClusterDeltaIEtaDeltaIPhiE0"," Cluster size in columns vs rows for E < 2 GeV, n cells > 3",
1928                                                     50,0,50,50,0,50);
1929       fhBadClusterDeltaIEtaDeltaIPhiE0->SetXTitle("#Delta Column");
1930       fhBadClusterDeltaIEtaDeltaIPhiE0->SetYTitle("#Delta Row");
1931       outputContainer->Add(fhBadClusterDeltaIEtaDeltaIPhiE0);
1932
1933       fhBadClusterDeltaIEtaDeltaIPhiE2  = new TH2F ("hBadClusterDeltaIEtaDeltaIPhiE2"," Cluster size in columns vs rows for 2 <E < 6 GeV, n cells > 3",
1934                                                     50,0,50,50,0,50);
1935       fhBadClusterDeltaIEtaDeltaIPhiE2->SetXTitle("#Delta Column");
1936       fhBadClusterDeltaIEtaDeltaIPhiE2->SetYTitle("#Delta Row");
1937       outputContainer->Add(fhBadClusterDeltaIEtaDeltaIPhiE2);
1938
1939       fhBadClusterDeltaIEtaDeltaIPhiE6  = new TH2F ("hBadClusterDeltaIEtaDeltaIPhiE6"," Cluster size in columns vs rows for E > 6 GeV, n cells > 3",
1940                                                     50,0,50,50,0,50);
1941       fhBadClusterDeltaIEtaDeltaIPhiE6->SetXTitle("#Delta Column");
1942       fhBadClusterDeltaIEtaDeltaIPhiE6->SetYTitle("#Delta Row");
1943       outputContainer->Add(fhBadClusterDeltaIEtaDeltaIPhiE6);
1944
1945       fhBadClusterDeltaIA  = new TH2F ("hBadClusterDeltaIA"," Cluster *asymmetry* in cell units vs E",
1946                                        nptbins,ptmin,ptmax,21,-1.05,1.05);
1947       fhBadClusterDeltaIA->SetXTitle("E_{cluster}");
1948       fhBadClusterDeltaIA->SetYTitle("A_{cell in cluster}");
1949       outputContainer->Add(fhBadClusterDeltaIA);
1950     }
1951   }
1952
1953   if(fStudyWeight)
1954   {
1955     fhECellClusterRatio  = new TH2F ("hECellClusterRatio"," cell energy / cluster energy vs cluster energy",
1956                                      nptbins,ptmin,ptmax, 100,0,1.);
1957     fhECellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
1958     fhECellClusterRatio->SetYTitle("E_{cell i}/E_{cluster}");
1959     outputContainer->Add(fhECellClusterRatio);
1960
1961     fhECellClusterLogRatio  = new TH2F ("hECellClusterLogRatio"," Log(cell energy / cluster energy) vs cluster energy",
1962                                         nptbins,ptmin,ptmax, 100,-10,0);
1963     fhECellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
1964     fhECellClusterLogRatio->SetYTitle("Log(E_{cell i}/E_{cluster})");
1965     outputContainer->Add(fhECellClusterLogRatio);
1966
1967     fhEMaxCellClusterRatio  = new TH2F ("hEMaxCellClusterRatio"," max cell energy / cluster energy vs cluster energy",
1968                                         nptbins,ptmin,ptmax, 100,0,1.);
1969     fhEMaxCellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
1970     fhEMaxCellClusterRatio->SetYTitle("E_{max cell}/E_{cluster}");
1971     outputContainer->Add(fhEMaxCellClusterRatio);
1972
1973     fhEMaxCellClusterLogRatio  = new TH2F ("hEMaxCellClusterLogRatio"," Log(max cell energy / cluster energy) vs cluster energy",
1974                                            nptbins,ptmin,ptmax, 100,-10,0);
1975     fhEMaxCellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
1976     fhEMaxCellClusterLogRatio->SetYTitle("Log (E_{max cell}/E_{cluster})");
1977     outputContainer->Add(fhEMaxCellClusterLogRatio);
1978
1979     for(Int_t iw = 0; iw < 14; iw++){
1980       fhLambda0ForW0[iw]  = new TH2F (Form("hLambda0ForW0%d",iw),Form("shower shape, #lambda^{2}_{0} vs E, w0 = %1.1f",1+0.5*iw),
1981                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1982       fhLambda0ForW0[iw]->SetXTitle("E_{cluster}");
1983       fhLambda0ForW0[iw]->SetYTitle("#lambda^{2}_{0}");
1984       outputContainer->Add(fhLambda0ForW0[iw]);
1985
1986 //      fhLambda1ForW0[iw]  = new TH2F (Form("hLambda1ForW0%d",iw),Form("shower shape, #lambda^{2}_{1} vs E, w0 = %1.1f",1+0.5*iw),
1987 //                                      nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1988 //      fhLambda1ForW0[iw]->SetXTitle("E_{cluster}");
1989 //      fhLambda1ForW0[iw]->SetYTitle("#lambda^{2}_{1}");
1990 //      outputContainer->Add(fhLambda1ForW0[iw]);
1991
1992       if(IsDataMC()){
1993         TString mcnames[] = {"Photon", "Electron","Conversion","Pi0","Hadron"};
1994         for(Int_t imc = 0; imc < 5; imc++){
1995           fhLambda0ForW0MC[iw][imc]  = new TH2F (Form("hLambda0ForW0%d_MC%s",iw,mcnames[imc].Data()),
1996                                                  Form("shower shape, #lambda^{2}_{0} vs E, w0 = %1.1f, for MC %s",1+0.5*iw,mcnames[imc].Data()),
1997                                                  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1998           fhLambda0ForW0MC[iw][imc]->SetXTitle("E_{cluster}");
1999           fhLambda0ForW0MC[iw][imc]->SetYTitle("#lambda^{2}_{0}");
2000           outputContainer->Add(fhLambda0ForW0MC[iw][imc]);
2001
2002 //          fhLambda1ForW0MC[iw][imc]  = new TH2F (Form("hLambda1ForW0%d_MC%s",iw,mcnames[imc].Data()),
2003 //                                                 Form("shower shape, #lambda^{2}_{1} vs E, w0 = %1.1f, for MC %s",1+0.5*iw,mcnames[imc].Data()),
2004 //                                                 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2005 //          fhLambda1ForW0MC[iw][imc]->SetXTitle("E_{cluster}");
2006 //          fhLambda1ForW0MC[iw][imc]->SetYTitle("#lambda^{2}_{1}");
2007 //          outputContainer->Add(fhLambda1ForW0MC[iw][imc]);
2008         }
2009       }
2010     }
2011   }
2012
2013   //Track Matching
2014   if(fFillAllTMHisto)
2015   {
2016     if(fFillAllTH12)
2017     {
2018       fhECharged  = new TH1F ("hECharged","E reconstructed clusters, matched with track", nptbins,ptmin,ptmax);
2019       fhECharged->SetXTitle("E (GeV)");
2020       outputContainer->Add(fhECharged);
2021
2022       fhPtCharged  = new TH1F ("hPtCharged","p_{T} reconstructed clusters, matched with track", nptbins,ptmin,ptmax);
2023       fhPtCharged->SetXTitle("p_{T} (GeV/c)");
2024       outputContainer->Add(fhPtCharged);
2025
2026       fhPhiCharged  = new TH1F ("hPhiCharged","#phi reconstructed clusters, matched with track",nphibins,phimin,phimax);
2027       fhPhiCharged->SetXTitle("#phi (rad)");
2028       outputContainer->Add(fhPhiCharged);
2029
2030       fhEtaCharged  = new TH1F ("hEtaCharged","#eta reconstructed clusters, matched with track",netabins,etamin,etamax);
2031       fhEtaCharged->SetXTitle("#eta ");
2032       outputContainer->Add(fhEtaCharged);
2033     }
2034     if(fFillAllTH3){
2035       fhEtaPhiECharged  = new TH3F ("hEtaPhiECharged","#eta vs #phi, reconstructed clusters, matched with track",
2036                                     netabins,etamin,etamax,nphibins,phimin,phimax,nptbins,ptmin,ptmax);
2037       fhEtaPhiECharged->SetXTitle("#eta ");
2038       fhEtaPhiECharged->SetYTitle("#phi ");
2039       fhEtaPhiECharged->SetZTitle("E (GeV) ");
2040       outputContainer->Add(fhEtaPhiECharged);
2041     }
2042
2043     fh1EOverP = new TH2F("h1EOverP","TRACK matches E/p",nptbins,ptmin,ptmax, nPoverEbins,EOverPmin,EOverPmax);
2044     fh1EOverP->SetYTitle("E/p");
2045     fh1EOverP->SetXTitle("p_{T} (GeV/c)");
2046     outputContainer->Add(fh1EOverP);
2047
2048     fh2dR = new TH2F("h2dR","TRACK matches dR",nptbins,ptmin,ptmax,ndRbins,dRmin,dRmax);
2049     fh2dR->SetXTitle("#Delta R (rad)");
2050     fh2dR->SetXTitle("E cluster (GeV)");
2051     outputContainer->Add(fh2dR) ;
2052
2053     fh2MatchdEdx = new TH2F("h2MatchdEdx","dE/dx vs. p for all matches",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
2054     fh2MatchdEdx->SetXTitle("p (GeV/c)");
2055     fh2MatchdEdx->SetYTitle("<dE/dx>");
2056     outputContainer->Add(fh2MatchdEdx);
2057
2058     fh2EledEdx = new TH2F("h2EledEdx","dE/dx vs. p for electrons",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
2059     fh2EledEdx->SetXTitle("p (GeV/c)");
2060     fh2EledEdx->SetYTitle("<dE/dx>");
2061     outputContainer->Add(fh2EledEdx) ;
2062
2063     fh1EOverPR02 = new TH2F("h1EOverPR02","TRACK matches E/p, all",nptbins,ptmin,ptmax, nPoverEbins,EOverPmin,EOverPmax);
2064     fh1EOverPR02->SetYTitle("E/p");
2065     fh1EOverPR02->SetXTitle("p_{T} (GeV/c)");
2066     outputContainer->Add(fh1EOverPR02);
2067   }
2068
2069   if(fFillAllPi0Histo)
2070   {
2071     fhIM  = new TH2F ("hIM","Cluster pairs Invariant mass vs reconstructed pair energy, ncell > 1",nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
2072     fhIM->SetXTitle("p_{T, cluster pairs} (GeV) ");
2073     fhIM->SetYTitle("M_{cluster pairs} (GeV/c^{2})");
2074     outputContainer->Add(fhIM);
2075
2076     fhAsym  = new TH2F ("hAssym","Cluster pairs Asymmetry vs reconstructed pair energy",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax);
2077     fhAsym->SetXTitle("p_{T, cluster pairs} (GeV) ");
2078     fhAsym->SetYTitle("Asymmetry");
2079     outputContainer->Add(fhAsym);
2080   }
2081
2082
2083   if(fFillAllPosHisto2)
2084   {
2085     if(fFillAllTH3)
2086     {
2087       fhXYZ  = new TH3F ("hXYZ","Cluster: x vs y vs z",xbins,xmin,xmax,ybins,ymin,ymax,zbins,zmin,zmax);
2088       fhXYZ->SetXTitle("x (cm)");
2089       fhXYZ->SetYTitle("y (cm)");
2090       fhXYZ->SetZTitle("z (cm) ");
2091       outputContainer->Add(fhXYZ);
2092     }
2093
2094     fhXNCells  = new TH2F ("hXNCells","Cluster X position vs N Cells per Cluster",xbins,xmin,xmax,nceclbins,nceclmin,nceclmax);
2095     fhXNCells->SetXTitle("x (cm)");
2096     fhXNCells->SetYTitle("N cells per cluster");
2097     outputContainer->Add(fhXNCells);
2098
2099     fhZNCells  = new TH2F ("hZNCells","Cluster Z position vs N Cells per Cluster",zbins,zmin,zmax,nceclbins,nceclmin,nceclmax);
2100     fhZNCells->SetXTitle("z (cm)");
2101     fhZNCells->SetYTitle("N cells per cluster");
2102     outputContainer->Add(fhZNCells);
2103
2104     fhXE  = new TH2F ("hXE","Cluster X position vs cluster energy",xbins,xmin,xmax,nptbins,ptmin,ptmax);
2105     fhXE->SetXTitle("x (cm)");
2106     fhXE->SetYTitle("E (GeV)");
2107     outputContainer->Add(fhXE);
2108
2109     fhZE  = new TH2F ("hZE","Cluster Z position vs cluster energy",zbins,zmin,zmax,nptbins,ptmin,ptmax);
2110     fhZE->SetXTitle("z (cm)");
2111     fhZE->SetYTitle("E (GeV)");
2112     outputContainer->Add(fhZE);
2113
2114     fhRNCells  = new TH2F ("hRNCells","Cluster R position vs N Cells per Cluster",rbins,rmin,rmax,nceclbins,nceclmin,nceclmax);
2115     fhRNCells->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
2116     fhRNCells->SetYTitle("N cells per cluster");
2117     outputContainer->Add(fhRNCells);
2118
2119
2120     fhYNCells  = new TH2F ("hYNCells","Cluster Y position vs N Cells per Cluster",ybins,ymin,ymax,nceclbins,nceclmin,nceclmax);
2121     fhYNCells->SetXTitle("y (cm)");
2122     fhYNCells->SetYTitle("N cells per cluster");
2123     outputContainer->Add(fhYNCells);
2124
2125     fhRE  = new TH2F ("hRE","Cluster R position vs cluster energy",rbins,rmin,rmax,nptbins,ptmin,ptmax);
2126     fhRE->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
2127     fhRE->SetYTitle("E (GeV)");
2128     outputContainer->Add(fhRE);
2129
2130     fhYE  = new TH2F ("hYE","Cluster Y position vs cluster energy",ybins,ymin,ymax,nptbins,ptmin,ptmax);
2131     fhYE->SetXTitle("y (cm)");
2132     fhYE->SetYTitle("E (GeV)");
2133     outputContainer->Add(fhYE);
2134   }
2135
2136   if(fFillAllPosHisto)
2137   {
2138     fhRCellE  = new TH2F ("hRCellE","Cell R position vs cell energy",rbins,rmin,rmax,nptbins,ptmin,ptmax);
2139     fhRCellE->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
2140     fhRCellE->SetYTitle("E (GeV)");
2141     outputContainer->Add(fhRCellE);
2142
2143     fhXCellE  = new TH2F ("hXCellE","Cell X position vs cell energy",xbins,xmin,xmax,nptbins,ptmin,ptmax);
2144     fhXCellE->SetXTitle("x (cm)");
2145     fhXCellE->SetYTitle("E (GeV)");
2146     outputContainer->Add(fhXCellE);
2147
2148     fhYCellE  = new TH2F ("hYCellE","Cell Y position vs cell energy",ybins,ymin,ymax,nptbins,ptmin,ptmax);
2149     fhYCellE->SetXTitle("y (cm)");
2150     fhYCellE->SetYTitle("E (GeV)");
2151     outputContainer->Add(fhYCellE);
2152
2153     fhZCellE  = new TH2F ("hZCellE","Cell Z position vs cell energy",zbins,zmin,zmax,nptbins,ptmin,ptmax);
2154     fhZCellE->SetXTitle("z (cm)");
2155     fhZCellE->SetYTitle("E (GeV)");
2156     outputContainer->Add(fhZCellE);
2157
2158     fhXYZCell  = new TH3F ("hXYZCell","Cell : x vs y vs z",xbins,xmin,xmax,ybins,ymin,ymax,zbins,zmin,zmax);
2159     fhXYZCell->SetXTitle("x (cm)");
2160     fhXYZCell->SetYTitle("y (cm)");
2161     fhXYZCell->SetZTitle("z (cm)");
2162     outputContainer->Add(fhXYZCell);
2163
2164
2165     Float_t dx = TMath::Abs(xmin)+TMath::Abs(xmax);
2166     Float_t dy = TMath::Abs(ymin)+TMath::Abs(ymax);
2167     Float_t dz = TMath::Abs(zmin)+TMath::Abs(zmax);
2168     Float_t dr = TMath::Abs(rmin)+TMath::Abs(rmax);
2169
2170     fhDeltaCellClusterRNCells  = new TH2F ("hDeltaCellClusterRNCells","Cluster-Cell R position vs N Cells per Cluster",rbins*2,-dr,dr,nceclbins,nceclmin,nceclmax);
2171     fhDeltaCellClusterRNCells->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
2172     fhDeltaCellClusterRNCells->SetYTitle("N cells per cluster");
2173     outputContainer->Add(fhDeltaCellClusterRNCells);
2174
2175     fhDeltaCellClusterXNCells  = new TH2F ("hDeltaCellClusterXNCells","Cluster-Cell X position vs N Cells per Cluster",xbins*2,-dx,dx,nceclbins,nceclmin,nceclmax);
2176     fhDeltaCellClusterXNCells->SetXTitle("x (cm)");
2177     fhDeltaCellClusterXNCells->SetYTitle("N cells per cluster");
2178     outputContainer->Add(fhDeltaCellClusterXNCells);
2179
2180     fhDeltaCellClusterYNCells  = new TH2F ("hDeltaCellClusterYNCells","Cluster-Cell Y position vs N Cells per Cluster",ybins*2,-dy,dy,nceclbins,nceclmin,nceclmax);
2181     fhDeltaCellClusterYNCells->SetXTitle("y (cm)");
2182     fhDeltaCellClusterYNCells->SetYTitle("N cells per cluster");
2183     outputContainer->Add(fhDeltaCellClusterYNCells);
2184
2185     fhDeltaCellClusterZNCells  = new TH2F ("hDeltaCellClusterZNCells","Cluster-Cell Z position vs N Cells per Cluster",zbins*2,-dz,dz,nceclbins,nceclmin,nceclmax);
2186     fhDeltaCellClusterZNCells->SetXTitle("z (cm)");
2187     fhDeltaCellClusterZNCells->SetYTitle("N cells per cluster");
2188     outputContainer->Add(fhDeltaCellClusterZNCells);
2189
2190     fhDeltaCellClusterRE  = new TH2F ("hDeltaCellClusterRE","Cluster-Cell R position vs cluster energy",rbins*2,-dr,dr,nptbins,ptmin,ptmax);
2191     fhDeltaCellClusterRE->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
2192     fhDeltaCellClusterRE->SetYTitle("E (GeV)");
2193     outputContainer->Add(fhDeltaCellClusterRE);
2194
2195     fhDeltaCellClusterXE  = new TH2F ("hDeltaCellClusterXE","Cluster-Cell X position vs cluster energy",xbins*2,-dx,dx,nptbins,ptmin,ptmax);
2196     fhDeltaCellClusterXE->SetXTitle("x (cm)");
2197     fhDeltaCellClusterXE->SetYTitle("E (GeV)");
2198     outputContainer->Add(fhDeltaCellClusterXE);
2199
2200     fhDeltaCellClusterYE  = new TH2F ("hDeltaCellClusterYE","Cluster-Cell Y position vs cluster energy",ybins*2,-dy,dy,nptbins,ptmin,ptmax);
2201     fhDeltaCellClusterYE->SetXTitle("y (cm)");
2202     fhDeltaCellClusterYE->SetYTitle("E (GeV)");
2203     outputContainer->Add(fhDeltaCellClusterYE);
2204
2205     fhDeltaCellClusterZE  = new TH2F ("hDeltaCellClusterZE","Cluster-Cell Z position vs cluster energy",zbins*2,-dz,dz,nptbins,ptmin,ptmax);
2206     fhDeltaCellClusterZE->SetXTitle("z (cm)");
2207     fhDeltaCellClusterZE->SetYTitle("E (GeV)");
2208     outputContainer->Add(fhDeltaCellClusterZE);
2209
2210     fhEtaPhiAmp  = new TH3F ("hEtaPhiAmp","Cell #eta vs cell #phi vs cell energy",netabins,etamin,etamax,nphibins,phimin,phimax,nptbins,ptmin,ptmax);
2211     fhEtaPhiAmp->SetXTitle("#eta ");
2212     fhEtaPhiAmp->SetYTitle("#phi (rad)");
2213     fhEtaPhiAmp->SetZTitle("E (GeV) ");
2214     outputContainer->Add(fhEtaPhiAmp);
2215
2216   }
2217
2218   //Calo cells
2219   fhNCells  = new TH1F ("hNCells","# cells", ncebins,ncemin+0.5,ncemax);
2220   fhNCells->SetXTitle("n cells");
2221   outputContainer->Add(fhNCells);
2222
2223   fhAmplitude  = new TH1F ("hAmplitude","Cell Energy", nptbins*2,ptmin,ptmax);
2224   fhAmplitude->SetXTitle("Cell Energy (GeV)");
2225   outputContainer->Add(fhAmplitude);
2226
2227   fhAmpId  = new TH2F ("hAmpId","Cell Energy", nfineptbins,ptfinemin,ptfinemax,fNMaxRows*fNMaxCols*fNModules,0,fNMaxRows*fNMaxCols*fNModules);
2228   fhAmpId->SetXTitle("Cell Energy (GeV)");
2229   outputContainer->Add(fhAmpId);
2230
2231   if(fFillAllCellTimeHisto)
2232   {
2233     fhCellTimeSpreadRespectToCellMax = new TH2F ("hCellTimeSpreadRespectToCellMax","t_{cell max}-t_{cell i} per cluster", nptbins,ptmin,ptmax,tdbins,tdmin,tdmax);
2234     fhCellTimeSpreadRespectToCellMax->SetXTitle("E (GeV)");
2235     fhCellTimeSpreadRespectToCellMax->SetYTitle("#Delta t_{cell max-i} (ns)");
2236     outputContainer->Add(fhCellTimeSpreadRespectToCellMax);
2237
2238     fhClusterMaxCellDiffAverageTime = new TH2F ("hClusterMaxCellDiffAverageTime","t_{cell max}-t_{average} per cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
2239     fhClusterMaxCellDiffAverageTime->SetXTitle("E (GeV)");
2240     fhClusterMaxCellDiffAverageTime->SetYTitle("#Delta t_{cell max - average} (ns)");
2241     outputContainer->Add(fhClusterMaxCellDiffAverageTime);
2242
2243     fhClusterMaxCellDiffWeightedTime = new TH2F ("hClusterMaxCellDiffWeightedTime","t_{cell max}-t_{weighted} per cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
2244     fhClusterMaxCellDiffWeightedTime->SetXTitle("E (GeV)");
2245     fhClusterMaxCellDiffWeightedTime->SetYTitle("#Delta t_{cell max - weighted} (ns)");
2246     outputContainer->Add(fhClusterMaxCellDiffWeightedTime);
2247
2248     fhCellIdCellLargeTimeSpread= new TH1F ("hCellIdCellLargeTimeSpread","Cells with time 100 ns larger than cell max in cluster ",
2249                                            fNMaxCols*fNMaxRows*fNModules,0,fNMaxCols*fNMaxRows*fNModules);
2250     fhCellIdCellLargeTimeSpread->SetXTitle("Absolute Cell Id");
2251     outputContainer->Add(fhCellIdCellLargeTimeSpread);
2252
2253     fhTime  = new TH1F ("hTime","Cell Time",ntimebins,timemin,timemax);
2254     fhTime->SetXTitle("Cell Time (ns)");
2255     outputContainer->Add(fhTime);
2256
2257     fhTimeVz  = new TH2F ("hTimeVz","Cell Time vs vertex, amplitude > 0.5 GeV",100, 0, 50,ntimebins,timemin,timemax);
2258     fhTimeVz->SetXTitle("|v_{z}| (cm)");
2259     fhTimeVz->SetYTitle("Cell Time (ns)");
2260     outputContainer->Add(fhTimeVz);
2261
2262     fhTimeId  = new TH2F ("hTimeId","Cell Time vs Absolute Id",
2263                           ntimebins,timemin,timemax,fNMaxRows*fNMaxCols*fNModules,0,fNMaxRows*fNMaxCols*fNModules);
2264     fhTimeId->SetXTitle("Cell Time (ns)");
2265     fhTimeId->SetYTitle("Cell Absolute Id");
2266     outputContainer->Add(fhTimeId);
2267
2268     fhTimeAmp  = new TH2F ("hTimeAmp","Cell Time vs Cell Energy",nptbins*2,ptmin,ptmax,ntimebins,timemin,timemax);
2269     fhTimeAmp->SetYTitle("Cell Time (ns)");
2270     fhTimeAmp->SetXTitle("Cell Energy (GeV)");
2271     outputContainer->Add(fhTimeAmp);
2272
2273   }
2274
2275   fhCellECross  = new TH2F ("hCellECross","1 - Energy in cross around cell /  cell energy",
2276                             nptbins,ptmin,ptmax, 400,-1,1.);
2277   fhCellECross->SetXTitle("E_{cell} (GeV) ");
2278   fhCellECross->SetYTitle("1- E_{cross}/E_{cell}");
2279   outputContainer->Add(fhCellECross);
2280
2281
2282   if(fCorrelate)
2283   {
2284     //PHOS vs EMCAL
2285     fhCaloCorrNClusters  = new TH2F ("hCaloCorrNClusters","# clusters in EMCAL vs PHOS", nclbins,nclmin,nclmax,nclbins,nclmin,nclmax);
2286     fhCaloCorrNClusters->SetXTitle("number of clusters in EMCAL");
2287     fhCaloCorrNClusters->SetYTitle("number of clusters in PHOS");
2288     outputContainer->Add(fhCaloCorrNClusters);
2289
2290     fhCaloCorrEClusters  = new TH2F ("hCaloCorrEClusters","summed energy of clusters in EMCAL vs PHOS", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2291     fhCaloCorrEClusters->SetXTitle("#Sigma E of clusters in EMCAL (GeV)");
2292     fhCaloCorrEClusters->SetYTitle("#Sigma E of clusters in PHOS (GeV)");
2293     outputContainer->Add(fhCaloCorrEClusters);
2294
2295     fhCaloCorrNCells  = new TH2F ("hCaloCorrNCells","# Cells in EMCAL vs PHOS", ncebins,ncemin,ncemax, ncebins,ncemin,ncemax);
2296     fhCaloCorrNCells->SetXTitle("number of Cells in EMCAL");
2297     fhCaloCorrNCells->SetYTitle("number of Cells in PHOS");
2298     outputContainer->Add(fhCaloCorrNCells);
2299
2300     fhCaloCorrECells  = new TH2F ("hCaloCorrECells","summed energy of Cells in EMCAL vs PHOS", nptbins*2,ptmin,ptmax*2,nptbins*2,ptmin,ptmax*2);
2301     fhCaloCorrECells->SetXTitle("#Sigma E of Cells in EMCAL (GeV)");
2302     fhCaloCorrECells->SetYTitle("#Sigma E of Cells in PHOS (GeV)");
2303     outputContainer->Add(fhCaloCorrECells);
2304
2305     //Calorimeter VS V0 signal
2306     fhCaloV0SCorrNClusters  = new TH2F ("hCaloV0SNClusters",Form("# clusters in %s vs V0 signal",fCalorimeter.Data()), nv0sbins,nv0smin,nv0smax,nclbins,nclmin,nclmax);
2307     fhCaloV0SCorrNClusters->SetXTitle("V0 signal");
2308     fhCaloV0SCorrNClusters->SetYTitle(Form("number of clusters in %s",fCalorimeter.Data()));
2309     outputContainer->Add(fhCaloV0SCorrNClusters);
2310
2311     fhCaloV0SCorrEClusters  = new TH2F ("hCaloV0SEClusters",Form("summed energy of clusters in %s vs V0 signal",fCalorimeter.Data()), nv0sbins,nv0smin,nv0smax,nptbins,ptmin,ptmax);
2312     fhCaloV0SCorrEClusters->SetXTitle("V0 signal");
2313     fhCaloV0SCorrEClusters->SetYTitle(Form("#Sigma E of clusters in %s (GeV)",fCalorimeter.Data()));
2314     outputContainer->Add(fhCaloV0SCorrEClusters);
2315
2316     fhCaloV0SCorrNCells  = new TH2F ("hCaloV0SNCells",Form("# Cells in %s vs V0 signal",fCalorimeter.Data()), nv0sbins,nv0smin,nv0smax, ncebins,ncemin,ncemax);
2317     fhCaloV0SCorrNCells->SetXTitle("V0 signal");
2318     fhCaloV0SCorrNCells->SetYTitle(Form("number of Cells in %s",fCalorimeter.Data()));
2319     outputContainer->Add(fhCaloV0SCorrNCells);
2320
2321     fhCaloV0SCorrECells  = new TH2F ("hCaloV0SECells",Form("summed energy of Cells in %s vs V0 signal",fCalorimeter.Data()), nv0sbins,nv0smin,nv0smax,nptbins,ptmin,ptmax);
2322     fhCaloV0SCorrECells->SetXTitle("V0 signal");
2323     fhCaloV0SCorrECells->SetYTitle(Form("#Sigma E of Cells in %s (GeV)",fCalorimeter.Data()));
2324     outputContainer->Add(fhCaloV0SCorrECells);
2325
2326     //Calorimeter VS V0 multiplicity
2327     fhCaloV0MCorrNClusters  = new TH2F ("hCaloV0MNClusters",Form("# clusters in %s vs V0 signal",fCalorimeter.Data()), nv0mbins,nv0mmin,nv0mmax,nclbins,nclmin,nclmax);
2328     fhCaloV0MCorrNClusters->SetXTitle("V0 signal");
2329     fhCaloV0MCorrNClusters->SetYTitle(Form("number of clusters in %s",fCalorimeter.Data()));
2330     outputContainer->Add(fhCaloV0MCorrNClusters);
2331
2332     fhCaloV0MCorrEClusters  = new TH2F ("hCaloV0MEClusters",Form("summed energy of clusters in %s vs V0 signal",fCalorimeter.Data()), nv0mbins,nv0mmin,nv0mmax,nptbins,ptmin,ptmax);
2333     fhCaloV0MCorrEClusters->SetXTitle("V0 signal");
2334     fhCaloV0MCorrEClusters->SetYTitle(Form("#Sigma E of clusters in %s (GeV)",fCalorimeter.Data()));
2335     outputContainer->Add(fhCaloV0MCorrEClusters);
2336
2337     fhCaloV0MCorrNCells  = new TH2F ("hCaloV0MNCells",Form("# Cells in %s vs V0 signal",fCalorimeter.Data()), nv0mbins,nv0mmin,nv0mmax, ncebins,ncemin,ncemax);
2338     fhCaloV0MCorrNCells->SetXTitle("V0 signal");
2339     fhCaloV0MCorrNCells->SetYTitle(Form("number of Cells in %s",fCalorimeter.Data()));
2340     outputContainer->Add(fhCaloV0MCorrNCells);
2341
2342     fhCaloV0MCorrECells  = new TH2F ("hCaloV0MECells",Form("summed energy of Cells in %s vs V0 signal",fCalorimeter.Data()), nv0mbins,nv0mmin,nv0mmax,nptbins,ptmin,ptmax);
2343     fhCaloV0MCorrECells->SetXTitle("V0 signal");
2344     fhCaloV0MCorrECells->SetYTitle(Form("#Sigma E of Cells in %s (GeV)",fCalorimeter.Data()));
2345     outputContainer->Add(fhCaloV0MCorrECells);
2346
2347     //Calorimeter VS Track multiplicity
2348     fhCaloTrackMCorrNClusters  = new TH2F ("hCaloTrackMNClusters",Form("# clusters in %s vs # tracks",fCalorimeter.Data()), ntrmbins,ntrmmin,ntrmmax,nclbins,nclmin,nclmax);
2349     fhCaloTrackMCorrNClusters->SetXTitle("# tracks");
2350     fhCaloTrackMCorrNClusters->SetYTitle(Form("number of clusters in %s",fCalorimeter.Data()));
2351     outputContainer->Add(fhCaloTrackMCorrNClusters);
2352
2353     fhCaloTrackMCorrEClusters  = new TH2F ("hCaloTrackMEClusters",Form("summed energy of clusters in %s vs # tracks",fCalorimeter.Data()), ntrmbins,ntrmmin,ntrmmax,nptbins,ptmin,ptmax);
2354     fhCaloTrackMCorrEClusters->SetXTitle("# tracks");
2355     fhCaloTrackMCorrEClusters->SetYTitle(Form("#Sigma E of clusters in %s (GeV)",fCalorimeter.Data()));
2356     outputContainer->Add(fhCaloTrackMCorrEClusters);
2357
2358     fhCaloTrackMCorrNCells  = new TH2F ("hCaloTrackMNCells",Form("# Cells in %s vs # tracks",fCalorimeter.Data()), ntrmbins,ntrmmin,ntrmmax, ncebins,ncemin,ncemax);
2359     fhCaloTrackMCorrNCells->SetXTitle("# tracks");
2360     fhCaloTrackMCorrNCells->SetYTitle(Form("number of Cells in %s",fCalorimeter.Data()));
2361     outputContainer->Add(fhCaloTrackMCorrNCells);
2362
2363     fhCaloTrackMCorrECells  = new TH2F ("hCaloTrackMECells",Form("summed energy of Cells in %s vs # tracks",fCalorimeter.Data()), ntrmbins,ntrmmin,ntrmmax,nptbins,ptmin,ptmax);
2364     fhCaloTrackMCorrECells->SetXTitle("# tracks");
2365     fhCaloTrackMCorrECells->SetYTitle(Form("#Sigma E of Cells in %s (GeV)",fCalorimeter.Data()));
2366     outputContainer->Add(fhCaloTrackMCorrECells);
2367
2368
2369   }//correlate calorimeters
2370
2371   //Module histograms
2372
2373   fhEMod  = new TH2F ("hE_Mod","Cluster reconstructed Energy in each present Module",nptbins,ptmin,ptmax,fNModules,0,fNModules);
2374   fhEMod->SetXTitle("E (GeV)");
2375   fhEMod->SetYTitle("Module");
2376   outputContainer->Add(fhEMod);
2377
2378   fhAmpMod  = new TH2F ("hAmp_Mod","Cell energy in each present Module",nptbins,ptmin,ptmax,fNModules,0,fNModules);
2379   fhAmpMod->SetXTitle("E (GeV)");
2380   fhAmpMod->SetYTitle("Module");
2381   outputContainer->Add(fhAmpMod);
2382
2383   if(fFillAllCellTimeHisto)
2384   {
2385     fhTimeMod  = new TH2F ("hTime_Mod","Cell time in each present Module",ntimebins,timemin,timemax,fNModules,0,fNModules);
2386     fhTimeMod->SetXTitle("t (ns)");
2387     fhTimeMod->SetYTitle("Module");
2388     outputContainer->Add(fhTimeMod);
2389   }
2390
2391   fhNClustersMod  = new TH2F ("hNClusters_Mod","# clusters vs Module", nclbins,nclmin+0.5,nclmax,fNModules,0,fNModules);
2392   fhNClustersMod->SetXTitle("number of clusters");
2393   fhNClustersMod->SetYTitle("Module");
2394   outputContainer->Add(fhNClustersMod);
2395
2396   fhNCellsMod  = new TH2F ("hNCells_Mod","# cells vs Module", ncebins,ncemin+0.5,ncemax,fNModules,0,fNModules);
2397   fhNCellsMod->SetXTitle("n cells");
2398   fhNCellsMod->SetYTitle("Module");
2399   outputContainer->Add(fhNCellsMod);
2400
2401   Int_t colmaxs = fNMaxCols;
2402   Int_t rowmaxs = fNMaxRows;
2403   if(fCalorimeter=="EMCAL")
2404   {
2405     colmaxs=2*fNMaxCols;
2406     rowmaxs=Int_t(fNModules/2)*fNMaxRows;
2407   }
2408   else
2409   {
2410     rowmaxs=fNModules*fNMaxRows;
2411   }
2412
2413   fhGridCells  = new TH2F ("hGridCells",Form("Entries in grid of cells"),
2414                            colmaxs+2,-1.5,colmaxs+0.5, rowmaxs+2,-1.5,rowmaxs+0.5);
2415   fhGridCells->SetYTitle("row (phi direction)");
2416   fhGridCells->SetXTitle("column (eta direction)");
2417   outputContainer->Add(fhGridCells);
2418
2419   fhGridCellsE  = new TH2F ("hGridCellsE","Accumulated energy in grid of cells",
2420                             colmaxs+2,-1.5,colmaxs+0.5, rowmaxs+2,-1.5,rowmaxs+0.5);
2421   fhGridCellsE->SetYTitle("row (phi direction)");
2422   fhGridCellsE->SetXTitle("column (eta direction)");
2423   outputContainer->Add(fhGridCellsE);
2424
2425   if(fFillAllCellTimeHisto)
2426   {
2427     fhGridCellsTime  = new TH2F ("hGridCellsTime","Accumulated time in grid of cells",
2428                                  colmaxs+2,-1.5,colmaxs+0.5, rowmaxs+2,-1.5,rowmaxs+0.5);
2429     fhGridCellsTime->SetYTitle("row (phi direction)");
2430     fhGridCellsTime->SetXTitle("column (eta direction)");
2431     outputContainer->Add(fhGridCellsTime);
2432   }
2433
2434   fhNCellsPerClusterMod      = new TH2F*[fNModules];
2435   fhNCellsPerClusterModNoCut = new TH2F*[fNModules];
2436   fhIMMod                    = new TH2F*[fNModules];
2437   if(fFillAllCellTimeHisto) fhTimeAmpPerRCU = new TH2F*[fNModules*fNRCU];
2438
2439   for(Int_t imod = 0; imod < fNModules; imod++)
2440   {
2441     fhNCellsPerClusterMod[imod]  = new TH2F (Form("hNCellsPerCluster_Mod%d",imod),
2442                                              Form("# cells per cluster vs cluster energy in Module %d",imod),
2443                                              nptbins,ptmin,ptmax, nceclbins,nceclmin,nceclmax);
2444     fhNCellsPerClusterMod[imod]->SetXTitle("E (GeV)");
2445     fhNCellsPerClusterMod[imod]->SetYTitle("n cells");
2446     outputContainer->Add(fhNCellsPerClusterMod[imod]);
2447
2448     fhNCellsPerClusterModNoCut[imod]  = new TH2F (Form("hNCellsPerClusterNoCut_Mod%d",imod),
2449                                                   Form("# cells per cluster vs cluster energy in Module %d, no cut",imod),
2450                                                   nptbins,ptmin,ptmax, nceclbins,nceclmin,nceclmax);
2451     fhNCellsPerClusterModNoCut[imod]->SetXTitle("E (GeV)");
2452     fhNCellsPerClusterModNoCut[imod]->SetYTitle("n cells");
2453     outputContainer->Add(fhNCellsPerClusterModNoCut[imod]);
2454
2455     if(fFillAllCellTimeHisto)
2456     {
2457       for(Int_t ircu = 0; ircu < fNRCU; ircu++)
2458       {
2459         fhTimeAmpPerRCU[imod*fNRCU+ircu]  = new TH2F (Form("hTimeAmp_Mod%d_RCU%d",imod,ircu),
2460                                                       Form("Cell Energy vs Cell Time in Module %d, RCU %d ",imod,ircu),
2461                                                       nptbins,ptmin,ptmax,ntimebins,timemin,timemax);
2462         fhTimeAmpPerRCU[imod*fNRCU+ircu]->SetXTitle("E (GeV)");
2463         fhTimeAmpPerRCU[imod*fNRCU+ircu]->SetYTitle("time (ns)");
2464         outputContainer->Add(fhTimeAmpPerRCU[imod*fNRCU+ircu]);
2465
2466       }
2467     }
2468
2469     if(fFillAllPi0Histo)
2470     {
2471       fhIMMod[imod]  = new TH2F (Form("hIM_Mod%d",imod),
2472                                  Form("Cluster pairs Invariant mass vs reconstructed pair energy in Module %d, n cell > 1",imod),
2473                                  nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
2474       fhIMMod[imod]->SetXTitle("p_{T, cluster pairs} (GeV) ");
2475       fhIMMod[imod]->SetYTitle("M_{cluster pairs} (GeV/c^{2})");
2476       outputContainer->Add(fhIMMod[imod]);
2477
2478     }
2479   }
2480
2481   // Monte Carlo Histograms
2482
2483   TString particleName[] = { "Photon", "Pi0", "Eta", "Electron", "NeutralHadron", "ChargedHadron" };
2484
2485   if(IsDataMC())
2486   {
2487     for(Int_t iPart = 0; iPart < 6; iPart++)
2488     {
2489       for(Int_t iCh = 0; iCh < 2; iCh++)
2490       {
2491         fhRecoMCRatioE[iPart][iCh]  = new TH2F (Form("hRecoMCRatioE_%s_Match%d",particleName[iPart].Data(),iCh),
2492                                                 Form("Reco/Gen E, %s, Matched %d",particleName[iPart].Data(),iCh),
2493                                                 nptbins, ptmin, ptmax, 200,0,2);
2494         fhRecoMCRatioE[iPart][iCh]->SetXTitle("E_{reconstructed}/E_{generated}");
2495         outputContainer->Add(fhRecoMCRatioE[iPart][iCh]);
2496
2497
2498         fhRecoMCDeltaE[iPart][iCh]  = new TH2F (Form("hRecoMCDeltaE_%s_Match%d",particleName[iPart].Data(),iCh),
2499                                                 Form("MC - Reco E, %s, Matched %d",particleName[iPart].Data(),iCh),
2500                                                 nptbins, ptmin, ptmax, nptbins*2,-ptmax,ptmax);
2501         fhRecoMCDeltaE[iPart][iCh]->SetXTitle("#Delta E (GeV)");
2502         outputContainer->Add(fhRecoMCDeltaE[iPart][iCh]);
2503
2504         fhRecoMCDeltaPhi[iPart][iCh]  = new TH2F (Form("hRecoMCDeltaPhi_%s_Match%d",particleName[iPart].Data(),iCh),
2505                                                   Form("MC - Reco #phi, %s, Matched %d",particleName[iPart].Data(),iCh),
2506                                                   nptbins, ptmin, ptmax, nphibins*2,-phimax,phimax);
2507         fhRecoMCDeltaPhi[iPart][iCh]->SetXTitle("#Delta #phi (rad)");
2508         outputContainer->Add(fhRecoMCDeltaPhi[iPart][iCh]);
2509
2510         fhRecoMCDeltaEta[iPart][iCh]  = new TH2F (Form("hRecoMCDeltaEta_%s_Match%d",particleName[iPart].Data(),iCh),
2511                                                   Form("MC- Reco #eta, %s, Matched %d",particleName[iPart].Data(),iCh),
2512                                                   nptbins, ptmin, ptmax,netabins*2,-etamax,etamax);
2513         fhRecoMCDeltaEta[iPart][iCh]->SetXTitle("#Delta #eta ");
2514         outputContainer->Add(fhRecoMCDeltaEta[iPart][iCh]);
2515
2516         fhRecoMCE[iPart][iCh]  = new TH2F (Form("hRecoMCE_%s_Match%d",particleName[iPart].Data(),iCh),
2517                                            Form("E distribution, reconstructed vs generated, %s, Matched %d",particleName[iPart].Data(),iCh),
2518                                            nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2519         fhRecoMCE[iPart][iCh]->SetXTitle("E_{rec} (GeV)");
2520         fhRecoMCE[iPart][iCh]->SetYTitle("E_{gen} (GeV)");
2521         outputContainer->Add(fhRecoMCE[iPart][iCh]);
2522
2523         fhRecoMCPhi[iPart][iCh]  = new TH2F (Form("hRecoMCPhi_%s_Match%d",particleName[iPart].Data(),iCh),
2524                                              Form("#phi distribution, reconstructed vs generated, %s, Matched %d",particleName[iPart].Data(),iCh),
2525                                              nphibins,phimin,phimax, nphibins,phimin,phimax);
2526         fhRecoMCPhi[iPart][iCh]->SetXTitle("#phi_{rec} (rad)");
2527         fhRecoMCPhi[iPart][iCh]->SetYTitle("#phi_{gen} (rad)");
2528         outputContainer->Add(fhRecoMCPhi[iPart][iCh]);
2529
2530         fhRecoMCEta[iPart][iCh]  = new TH2F (Form("hRecoMCEta_%s_Match%d",particleName[iPart].Data(),iCh),
2531                                              Form("#eta distribution, reconstructed vs generated, %s, Matched %d",particleName[iPart].Data(),iCh),
2532                                              netabins,etamin,etamax,netabins,etamin,etamax);
2533         fhRecoMCEta[iPart][iCh]->SetXTitle("#eta_{rec} ");
2534         fhRecoMCEta[iPart][iCh]->SetYTitle("#eta_{gen} ");
2535         outputContainer->Add(fhRecoMCEta[iPart][iCh]);
2536       }
2537     }
2538
2539     //Pure MC
2540     for(Int_t iPart = 0; iPart < 4; iPart++)
2541     {
2542       fhGenMCE[iPart]     = new TH1F(Form("hGenMCE_%s",particleName[iPart].Data()) ,
2543                                      Form("p_{T} of generated %s",particleName[iPart].Data()),
2544                                      nptbins,ptmin,ptmax);
2545       fhGenMCEtaPhi[iPart] = new TH2F(Form("hGenMCEtaPhi_%s",particleName[iPart].Data()),
2546                                       Form("Y vs #phi of generated %s",particleName[iPart].Data()),
2547                                       netabins,etamin,etamax,nphibins,phimin,phimax);
2548
2549       fhGenMCE[iPart]    ->SetXTitle("p_{T} (GeV/c)");
2550       fhGenMCEtaPhi[iPart]->SetXTitle("#eta");
2551       fhGenMCEtaPhi[iPart]->SetYTitle("#phi (rad)");
2552
2553       outputContainer->Add(fhGenMCE[iPart]);
2554       outputContainer->Add(fhGenMCEtaPhi[iPart]);
2555
2556
2557       fhGenMCAccE[iPart]     = new TH1F(Form("hGenMCAccE_%s",particleName[iPart].Data()) ,
2558                                         Form("p_{T} of generated %s",particleName[iPart].Data()),
2559                                         nptbins,ptmin,ptmax);
2560       fhGenMCAccEtaPhi[iPart] = new TH2F(Form("hGenMCAccEtaPhi_%s",particleName[iPart].Data()),
2561                                          Form("Y vs #phi of generated %s",particleName[iPart].Data()),
2562                                          netabins,etamin,etamax,nphibins,phimin,phimax);
2563
2564       fhGenMCAccE[iPart]    ->SetXTitle("p_{T} (GeV/c)");
2565       fhGenMCAccEtaPhi[iPart]->SetXTitle("#eta");
2566       fhGenMCAccEtaPhi[iPart]->SetYTitle("#phi (rad)");
2567
2568       outputContainer->Add(fhGenMCAccE[iPart]);
2569       outputContainer->Add(fhGenMCAccEtaPhi[iPart]);
2570
2571     }
2572
2573     //Vertex of generated particles
2574
2575     fhEMVxyz  = new TH2F ("hEMVxyz","Production vertex of reconstructed ElectroMagnetic particles",nvdistbins,vdistmin,vdistmax,nvdistbins,vdistmin,vdistmax);//,100,0,500);
2576     fhEMVxyz->SetXTitle("v_{x}");
2577     fhEMVxyz->SetYTitle("v_{y}");
2578     //fhEMVxyz->SetZTitle("v_{z}");
2579     outputContainer->Add(fhEMVxyz);
2580
2581     fhHaVxyz  = new TH2F ("hHaVxyz","Production vertex of reconstructed hadrons",nvdistbins,vdistmin,vdistmax,nvdistbins,vdistmin,vdistmax);//,100,0,500);
2582     fhHaVxyz->SetXTitle("v_{x}");
2583     fhHaVxyz->SetYTitle("v_{y}");
2584     //fhHaVxyz->SetZTitle("v_{z}");
2585     outputContainer->Add(fhHaVxyz);
2586
2587     fhEMR  = new TH2F ("hEMR","Distance to production vertex of reconstructed ElectroMagnetic particles vs E rec",nptbins,ptmin,ptmax,nvdistbins,vdistmin,vdistmax);
2588     fhEMR->SetXTitle("E (GeV)");
2589     fhEMR->SetYTitle("TMath::Sqrt(v_{x}^{2}+v_{y}^{2})");
2590     outputContainer->Add(fhEMR);
2591
2592     fhHaR  = new TH2F ("hHaR","Distance to production vertex of reconstructed Hadrons vs E rec",nptbins,ptmin,ptmax,nvdistbins,vdistmin,vdistmax);
2593     fhHaR->SetXTitle("E (GeV)");
2594     fhHaR->SetYTitle("TMath::Sqrt(v_{x}^{2}+v_{y}^{2})");
2595     outputContainer->Add(fhHaR);
2596
2597
2598     //Track Matching
2599
2600     fhMCEle1EOverP = new TH2F("hMCEle1EOverP","TRACK matches E/p, MC electrons",nptbins,ptmin,ptmax, nPoverEbins,EOverPmin,EOverPmax);
2601     fhMCEle1EOverP->SetYTitle("E/p");
2602     fhMCEle1EOverP->SetXTitle("p_{T} (GeV/c)");
2603     outputContainer->Add(fhMCEle1EOverP);
2604
2605     fhMCEle1dR = new TH1F("hMCEle1dR","TRACK matches dR, MC electrons",ndRbins,dRmin,dRmax);
2606     fhMCEle1dR->SetXTitle("#Delta R (rad)");
2607     outputContainer->Add(fhMCEle1dR) ;
2608
2609     fhMCEle2MatchdEdx = new TH2F("hMCEle2MatchdEdx","dE/dx vs. p for all matches, MC electrons",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
2610     fhMCEle2MatchdEdx->SetXTitle("p (GeV/c)");
2611     fhMCEle2MatchdEdx->SetYTitle("<dE/dx>");
2612     outputContainer->Add(fhMCEle2MatchdEdx);
2613
2614     fhMCChHad1EOverP = new TH2F("hMCChHad1EOverP","TRACK matches E/p, MC charged hadrons",nptbins,ptmin,ptmax, nPoverEbins,EOverPmin,EOverPmax);
2615     fhMCChHad1EOverP->SetYTitle("E/p");
2616     fhMCChHad1EOverP->SetXTitle("p_{T} (GeV/c)");
2617     outputContainer->Add(fhMCChHad1EOverP);
2618
2619     fhMCChHad1dR = new TH1F("hMCChHad1dR","TRACK matches dR, MC charged hadrons",ndRbins,dRmin,dRmax);
2620     fhMCChHad1dR->SetXTitle("#Delta R (rad)");
2621     outputContainer->Add(fhMCChHad1dR) ;
2622
2623     fhMCChHad2MatchdEdx = new TH2F("hMCChHad2MatchdEdx","dE/dx vs. p for all matches, MC charged hadrons",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
2624     fhMCChHad2MatchdEdx->SetXTitle("p (GeV/c)");
2625     fhMCChHad2MatchdEdx->SetYTitle("<dE/dx>");
2626     outputContainer->Add(fhMCChHad2MatchdEdx);
2627
2628     fhMCNeutral1EOverP = new TH2F("hMCNeutral1EOverP","TRACK matches E/p, MC neutrals",nptbins,ptmin,ptmax, nPoverEbins,EOverPmin,EOverPmax);
2629     fhMCNeutral1EOverP->SetYTitle("E/p");
2630     fhMCNeutral1EOverP->SetXTitle("p_{T} (GeV/c)");
2631     outputContainer->Add(fhMCNeutral1EOverP);
2632
2633     fhMCNeutral1dR = new TH1F("hMCNeutral1dR","TRACK matches dR, MC neutrals",ndRbins,dRmin,dRmax);
2634     fhMCNeutral1dR->SetXTitle("#Delta R (rad)");
2635     outputContainer->Add(fhMCNeutral1dR) ;
2636
2637     fhMCNeutral2MatchdEdx = new TH2F("hMCNeutral2MatchdEdx","dE/dx vs. p for all matches, MC neutrals",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
2638     fhMCNeutral2MatchdEdx->SetXTitle("p (GeV/c)");
2639     fhMCNeutral2MatchdEdx->SetYTitle("<dE/dx>");
2640     outputContainer->Add(fhMCNeutral2MatchdEdx);
2641
2642     fhMCEle1EOverPR02 = new TH2F("hMCEle1EOverPR02","TRACK matches E/p, MC electrons",nptbins,ptmin,ptmax, nPoverEbins,EOverPmin,EOverPmax);
2643     fhMCEle1EOverPR02->SetYTitle("E/p");
2644     fhMCEle1EOverPR02->SetXTitle("p_{T} (GeV/c)");
2645     outputContainer->Add(fhMCEle1EOverPR02);
2646
2647     fhMCChHad1EOverPR02 = new TH2F("hMCChHad1EOverPR02","TRACK matches E/p, MC charged hadrons",nptbins,ptmin,ptmax, nPoverEbins,EOverPmin,EOverPmax);
2648     fhMCChHad1EOverPR02->SetYTitle("E/p");
2649     fhMCChHad1EOverPR02->SetXTitle("p_{T} (GeV/c)");
2650     outputContainer->Add(fhMCChHad1EOverPR02);
2651
2652     fhMCNeutral1EOverPR02 = new TH2F("hMCNeutral1EOverPR02","TRACK matches E/p, MC neutrals",nptbins,ptmin,ptmax, nPoverEbins,EOverPmin,EOverPmax);
2653     fhMCNeutral1EOverPR02->SetYTitle("E/p");
2654     fhMCNeutral1EOverPR02->SetXTitle("p_{T} (GeV/c)");
2655     outputContainer->Add(fhMCNeutral1EOverPR02);
2656   }
2657
2658   //  for(Int_t i = 0; i < outputContainer->GetEntries() ; i++)
2659   //    printf("i=%d, name= %s\n",i,outputContainer->At(i)->GetName());
2660
2661   return outputContainer;
2662 }
2663
2664 //__________________________________________________________________________________________________
2665 Float_t AliAnaCalorimeterQA::GetECross(const Int_t absID, AliVCaloCells* cells, const Float_t dtcut)
2666 {
2667   // Get energy in cross axis around maximum cell, for EMCAL only
2668
2669   Int_t icol =-1, irow=-1,iRCU = -1;
2670   Int_t imod = GetModuleNumberCellIndexes(absID, fCalorimeter, icol, irow, iRCU);
2671
2672   if(fCalorimeter=="EMCAL")
2673   {
2674     //Get close cells index, energy and time, not in corners
2675
2676     Int_t absID1 = -1;
2677     Int_t absID2 = -1;
2678
2679     if( irow < AliEMCALGeoParams::fgkEMCALRows-1) absID1 = GetCaloUtils()->GetEMCALGeometry()->GetAbsCellIdFromCellIndexes(imod, irow+1, icol);
2680     if( irow > 0 )                                absID2 = GetCaloUtils()->GetEMCALGeometry()->GetAbsCellIdFromCellIndexes(imod, irow-1, icol);
2681
2682     // In case of cell in eta = 0 border, depending on SM shift the cross cell index
2683     Int_t absID3 = -1;
2684     Int_t absID4 = -1;
2685
2686     if     ( icol == AliEMCALGeoParams::fgkEMCALCols - 1 && !(imod%2) )
2687     {
2688       absID3 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod+1, irow, 0);
2689       absID4 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod  , irow, icol-1);
2690     }
2691     else if( icol == 0 && imod%2 )
2692     {
2693       absID3 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod  , irow, icol+1);
2694       absID4 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod-1, irow, AliEMCALGeoParams::fgkEMCALCols-1);
2695     }
2696     else
2697     {
2698       if( icol < AliEMCALGeoParams::fgkEMCALCols-1 )
2699         absID3 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod, irow, icol+1);
2700       if( icol > 0 )
2701         absID4 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod, irow, icol-1);
2702     }
2703
2704     //Recalibrate cell energy if needed
2705     //Float_t  ecell = cells->GetCellAmplitude(absID);
2706     //GetCaloUtils()->RecalibrateCellAmplitude(ecell,fCalorimeter, absID);
2707     Double_t tcell = cells->GetCellTime(absID);
2708     GetCaloUtils()->RecalibrateCellTime(tcell, fCalorimeter, absID,GetReader()->GetInputEvent()->GetBunchCrossNumber());
2709
2710     Float_t  ecell1  = 0, ecell2  = 0, ecell3  = 0, ecell4  = 0;
2711     Double_t tcell1  = 0, tcell2  = 0, tcell3  = 0, tcell4  = 0;
2712
2713     if(absID1 >0 )
2714     {
2715       ecell1 = cells->GetCellAmplitude(absID1);
2716       GetCaloUtils()->RecalibrateCellAmplitude(ecell1, fCalorimeter, absID1);
2717       tcell1 = cells->GetCellTime(absID1);
2718       GetCaloUtils()->RecalibrateCellTime     (tcell1, fCalorimeter, absID1,GetReader()->GetInputEvent()->GetBunchCrossNumber());
2719     }
2720     if(absID2 >0 )
2721     {
2722       ecell2 = cells->GetCellAmplitude(absID2);
2723       GetCaloUtils()->RecalibrateCellAmplitude(ecell2, fCalorimeter, absID2);
2724       tcell2 = cells->GetCellTime(absID2);
2725       GetCaloUtils()->RecalibrateCellTime     (tcell2, fCalorimeter, absID2, GetReader()->GetInputEvent()->GetBunchCrossNumber());
2726     }
2727     if(absID3 >0 )
2728     {
2729       ecell3 = cells->GetCellAmplitude(absID3);
2730       GetCaloUtils()->RecalibrateCellAmplitude(ecell3, fCalorimeter, absID3);
2731       tcell3 = cells->GetCellTime(absID3);
2732       GetCaloUtils()->RecalibrateCellTime     (tcell3, fCalorimeter, absID3, GetReader()->GetInputEvent()->GetBunchCrossNumber());
2733     }
2734     if(absID4 >0 )
2735     {
2736       ecell4 = cells->GetCellAmplitude(absID4);
2737       GetCaloUtils()->RecalibrateCellAmplitude(ecell4, fCalorimeter, absID4);
2738       tcell4 = cells->GetCellTime(absID4);
2739       GetCaloUtils()->RecalibrateCellTime     (tcell4, fCalorimeter, absID4, GetReader()->GetInputEvent()->GetBunchCrossNumber());
2740     }
2741
2742     if(TMath::Abs(tcell-tcell1)*1.e9 > dtcut) ecell1 = 0 ;
2743     if(TMath::Abs(tcell-tcell2)*1.e9 > dtcut) ecell2 = 0 ;
2744     if(TMath::Abs(tcell-tcell3)*1.e9 > dtcut) ecell3 = 0 ;
2745     if(TMath::Abs(tcell-tcell4)*1.e9 > dtcut) ecell4 = 0 ;
2746
2747     return ecell1+ecell2+ecell3+ecell4;
2748   }
2749   else //PHOS
2750   {
2751
2752     Int_t absId1 = -1, absId2 = -1, absId3 = -1, absId4 = -1;
2753
2754     Int_t relId1[] = { imod+1, 0, irow+1, icol   };
2755     Int_t relId2[] = { imod+1, 0, irow-1, icol   };
2756     Int_t relId3[] = { imod+1, 0, irow  , icol+1 };
2757     Int_t relId4[] = { imod+1, 0, irow  , icol-1 };
2758
2759     GetCaloUtils()->GetPHOSGeometry()->RelToAbsNumbering(relId1, absId1);
2760     GetCaloUtils()->GetPHOSGeometry()->RelToAbsNumbering(relId2, absId2);
2761     GetCaloUtils()->GetPHOSGeometry()->RelToAbsNumbering(relId3, absId3);
2762     GetCaloUtils()->GetPHOSGeometry()->RelToAbsNumbering(relId4, absId4);
2763
2764     Float_t  ecell1  = 0, ecell2  = 0, ecell3  = 0, ecell4  = 0;
2765
2766     if(absId1 > 0 ) ecell1 = cells->GetCellAmplitude(absId1);
2767     if(absId2 > 0 ) ecell2 = cells->GetCellAmplitude(absId2);
2768     if(absId3 > 0 ) ecell3 = cells->GetCellAmplitude(absId3);
2769     if(absId4 > 0 ) ecell4 = cells->GetCellAmplitude(absId4);
2770
2771     return ecell1+ecell2+ecell3+ecell4;
2772
2773   }
2774
2775 }
2776
2777 //__________________________________________________________________________________________________
2778 void AliAnaCalorimeterQA::InvariantMassHistograms(const Int_t iclus,   const TLorentzVector mom,
2779                                                   const Int_t nModule, const TObjArray* caloClusters,
2780                                                   AliVCaloCells * cells)
2781 {
2782   // Fill Invariant mass histograms
2783
2784   if(GetDebug()>1) printf("AliAnaCalorimeterQA::InvariantMassHistograms() - Start \n");
2785
2786   //Get vertex for photon momentum calculation and event selection
2787   Double_t v[3] = {0,0,0}; //vertex ;
2788   //GetReader()->GetVertex(v);
2789
2790   Int_t nModule2      = -1;
2791   TLorentzVector mom2 ;
2792   Int_t nCaloClusters = caloClusters->GetEntriesFast();
2793
2794   for(Int_t jclus = iclus + 1 ; jclus < nCaloClusters ; jclus++)
2795   {
2796     AliVCluster* clus2 =  (AliVCluster*)caloClusters->At(jclus);
2797
2798     Float_t maxCellFraction = 0.;
2799     Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(cells, clus2,maxCellFraction);
2800
2801     // Try to rediuce background with a mild shower shape cut and no more than 1 maxima
2802     // in cluster and remove low energy clusters
2803     if( clus2->GetNCells() <= 1 || !IsGoodCluster(absIdMax,cells) ||
2804        GetCaloUtils()->GetNumberOfLocalMaxima(clus2,cells) > 1 ||
2805        clus2->GetM02() > 0.5 || clus2->E() < 0.3) continue;
2806
2807     //Get cluster kinematics
2808     clus2->GetMomentum(mom2,v);
2809
2810     //Check only certain regions
2811     Bool_t in2 = kTRUE;
2812     if(IsFiducialCutOn()) in2 =  GetFiducialCut()->IsInFiducialCut(mom2,fCalorimeter) ;
2813     if(!in2) continue;
2814
2815     //Get module of cluster
2816     nModule2 = GetModuleNumber(clus2);
2817
2818     //Fill histograms
2819
2820     //All modules
2821     fhIM  ->Fill((mom+mom2).Pt(),(mom+mom2).M());
2822
2823     //Single module
2824     if(nModule == nModule2 && nModule >=0 && nModule < fNModules)
2825       fhIMMod[nModule]->Fill((mom+mom2).Pt(),(mom+mom2).M());
2826
2827
2828     //Asymetry histograms
2829     fhAsym->Fill((mom+mom2).Pt(),TMath::Abs((mom.E()-mom2.E())/(mom.E()+mom2.E())));
2830
2831   }// 2nd cluster loop
2832
2833 }
2834
2835 //______________________________
2836 void AliAnaCalorimeterQA::Init()
2837 {
2838   //Check if the data or settings are ok
2839
2840   if(fCalorimeter != "PHOS" && fCalorimeter !="EMCAL")
2841     AliFatal(Form("Wrong calorimeter name <%s>", fCalorimeter.Data()));
2842
2843   if(GetReader()->GetDataType()== AliCaloTrackReader::kMC)
2844     AliFatal("Analysis of reconstructed data, MC reader not aplicable");
2845
2846 }
2847
2848 //________________________________________
2849 void AliAnaCalorimeterQA::InitParameters()
2850 {
2851   //Initialize the parameters of the analysis.
2852   AddToHistogramsName("AnaCaloQA_");
2853
2854   fCalorimeter     = "EMCAL"; //or PHOS
2855   fNModules        = 12; // set maximum to maximum number of EMCAL modules
2856   fNRCU            = 2;  // set maximum number of RCU in EMCAL per SM
2857   fTimeCutMin      = -9999999;
2858   fTimeCutMax      =  9999999;
2859   fEMCALCellAmpMin = 0.2;
2860   fPHOSCellAmpMin  = 0.2;
2861
2862   // Exotic studies
2863   fExoNECrossCuts  = 10 ;
2864   fExoNDTimeCuts   = 4  ;
2865
2866   fExoDTimeCuts [0] = 1.e4 ; fExoDTimeCuts [1] = 50.0 ; fExoDTimeCuts [2] = 25.0 ; fExoDTimeCuts [3] = 10.0 ;
2867   fExoECrossCuts[0] = 0.80 ; fExoECrossCuts[1] = 0.85 ; fExoECrossCuts[2] = 0.90 ; fExoECrossCuts[3] = 0.92 ; fExoECrossCuts[4] = 0.94 ;
2868   fExoECrossCuts[5] = 0.95 ; fExoECrossCuts[6] = 0.96 ; fExoECrossCuts[7] = 0.97 ; fExoECrossCuts[8] = 0.98 ; fExoECrossCuts[9] = 0.99 ;
2869
2870 }
2871
2872 //___________________________________________________________________________________
2873 Bool_t AliAnaCalorimeterQA::IsGoodCluster(const Int_t absIdMax, AliVCaloCells* cells)
2874 {
2875   //Identify cluster as exotic or not
2876
2877   if(!fStudyBadClusters) return kTRUE;
2878
2879   if(fCalorimeter=="EMCAL")
2880   {
2881     if(!GetCaloUtils()->GetEMCALRecoUtils()->IsRejectExoticCluster())
2882     {
2883       return !( GetCaloUtils()->GetEMCALRecoUtils()->IsExoticCell(absIdMax,cells,(GetReader()->GetInputEvent())->GetBunchCrossNumber()) );
2884     }
2885     else
2886     {
2887       return kTRUE;
2888     }
2889   }
2890   else // PHOS
2891   {
2892     Float_t ampMax = cells->GetCellAmplitude(absIdMax);
2893     GetCaloUtils()->RecalibrateCellAmplitude(ampMax, fCalorimeter, absIdMax);
2894
2895     if(ampMax < 0.01) return kFALSE;
2896
2897     if(1-GetECross(absIdMax,cells)/ampMax > 0.95) return kFALSE;
2898     else                                          return kTRUE;
2899   }
2900
2901 }
2902
2903 //_________________________________________________________
2904 void AliAnaCalorimeterQA::Print(const Option_t * opt) const
2905 {
2906   //Print some relevant parameters set for the analysis
2907   if(! opt)
2908     return;
2909
2910   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
2911   AliAnaCaloTrackCorrBaseClass::Print(" ");
2912
2913   printf("Select Calorimeter %s \n",fCalorimeter.Data());
2914   printf("Time Cut: %3.1f < TOF  < %3.1f\n", fTimeCutMin, fTimeCutMax);
2915   printf("EMCAL Min Amplitude   : %2.1f GeV/c\n", fEMCALCellAmpMin) ;
2916   printf("PHOS Min Amplitude    : %2.1f GeV/c\n", fPHOSCellAmpMin) ;
2917
2918 }
2919
2920 //_____________________________________________________
2921 void  AliAnaCalorimeterQA::MakeAnalysisFillHistograms()
2922 {
2923   //Fill Calorimeter QA histograms
2924
2925   //Play with the MC stack if available
2926   if(IsDataMC()) MCHistograms();
2927
2928   //Get List with CaloClusters
2929   TObjArray * caloClusters = NULL;
2930   if      (fCalorimeter == "PHOS")  caloClusters = GetPHOSClusters();
2931   else if (fCalorimeter == "EMCAL") caloClusters = GetEMCALClusters();
2932   else
2933     AliFatal(Form("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - Wrong calorimeter name <%s>, END\n", fCalorimeter.Data()));
2934
2935   // Do not do anything if there are no clusters
2936   if(caloClusters->GetEntriesFast() == 0) return;
2937
2938   //Get List with CaloCells
2939   AliVCaloCells * cells = 0x0;
2940   if(fCalorimeter == "PHOS") cells =  GetPHOSCells();
2941   else                                   cells =  GetEMCALCells();
2942
2943   if(!caloClusters || !cells) {
2944     AliFatal(Form("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - No CaloClusters or CaloCells available\n"));
2945     return; // trick coverity
2946   }
2947
2948   //printf("QA: N cells %d, N clusters %d \n",cells->GetNumberOfCells(),caloClusters->GetEntriesFast());
2949
2950   // Correlate Calorimeters and V0 and track Multiplicity
2951   if(fCorrelate)        Correlate();
2952
2953   // Clusters
2954   ClusterLoopHistograms(caloClusters,cells);
2955
2956   // Cells
2957   CellHistograms(cells);
2958
2959   if(GetDebug() > 0)
2960     printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - End \n");
2961
2962 }
2963
2964 //______________________________________
2965 void AliAnaCalorimeterQA::MCHistograms()
2966 {
2967   //Get the MC arrays and do some checks before filling MC histograms
2968
2969   TLorentzVector mom  ;
2970
2971   if(GetReader()->ReadStack()){
2972
2973     if(!GetMCStack())
2974       AliFatal("Stack not available, is the MC handler called?\n");
2975
2976     //Fill some pure MC histograms, only primaries.
2977     for(Int_t i=0 ; i<GetMCStack()->GetNprimary(); i++)
2978     {//Only primary particles, for all MC transport put GetNtrack()
2979       TParticle *primary = GetMCStack()->Particle(i) ;
2980
2981       if (primary->GetStatusCode() > 11) continue; //Working for PYTHIA and simple generators, check for HERWIG
2982       primary->Momentum(mom);
2983       MCHistograms(mom,TMath::Abs(primary->GetPdgCode()));
2984     } //primary loop
2985   }
2986   else if(GetReader()->ReadAODMCParticles()){
2987
2988     if(!GetReader()->GetAODMCParticles(0))
2989       AliFatal("AODMCParticles not available!");
2990
2991     //Fill some pure MC histograms, only primaries.
2992     for(Int_t i=0 ; i < (GetReader()->GetAODMCParticles(0))->GetEntriesFast(); i++)
2993     {
2994       AliAODMCParticle *aodprimary = (AliAODMCParticle*) (GetReader()->GetAODMCParticles(0))->At(i) ;
2995
2996       if (!aodprimary->IsPrimary()) continue; //accept all which is not MC transport generated. Don't know how to avoid partons
2997
2998       mom.SetPxPyPzE(aodprimary->Px(),aodprimary->Py(),aodprimary->Pz(),aodprimary->E());
2999       MCHistograms(mom,TMath::Abs(aodprimary->GetPdgCode()));
3000     } //primary loop
3001
3002   }
3003 }
3004
3005 //_______________________________________________________________________________
3006 void AliAnaCalorimeterQA::MCHistograms(const TLorentzVector mom, const Int_t pdg)
3007 {
3008   //Fill pure monte carlo related histograms
3009
3010   Float_t eMC    = mom.E();
3011   Float_t phiMC  = mom.Phi();
3012   if(phiMC < 0)
3013     phiMC  += TMath::TwoPi();
3014   Float_t etaMC  = mom.Eta();
3015
3016   if (TMath::Abs(etaMC) > 1) return;
3017
3018   Bool_t in = kFALSE;
3019
3020   //Rough stimate of acceptance for pi0, Eta and electrons
3021   if(fCalorimeter == "PHOS")
3022   {
3023     if(GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter))
3024       in = kTRUE ;
3025     if(GetDebug() > 2) printf("AliAnaCalorimeterQA::MCHistograms() - In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),in);
3026
3027   }
3028   else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet())
3029   {
3030     if(GetEMCALGeometry())
3031     {
3032       Int_t absID=0;
3033       GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(mom.Eta(),mom.Phi(),absID);
3034
3035       if( absID >= 0)
3036         in = kTRUE;
3037
3038       if(GetDebug() > 2) printf("AliAnaCalorimeterQA::MCHistograms() - In %s Real acceptance? %d\n",fCalorimeter.Data(),in);
3039     }
3040     else
3041     {
3042       if(GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter))
3043         in = kTRUE ;
3044       if(GetDebug() > 2) printf("AliAnaCalorimeterQA::MCHistograms() - In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),in);
3045     }
3046   }
3047
3048   if (pdg==22)
3049   {
3050     fhGenMCE[kmcPhoton]    ->Fill(eMC);
3051     if(eMC > 0.5) fhGenMCEtaPhi[kmcPhoton]->Fill(etaMC,phiMC);
3052     if(in)
3053     {
3054       fhGenMCAccE[kmcPhoton]    ->Fill(eMC);
3055       if(eMC > 0.5) fhGenMCAccEtaPhi[kmcPhoton]->Fill(etaMC,phiMC);
3056     }
3057   }
3058   else if (pdg==111)
3059   {
3060     fhGenMCE[kmcPi0]    ->Fill(eMC);
3061     if(eMC > 0.5) fhGenMCEtaPhi[kmcPi0]->Fill(etaMC,phiMC);
3062     if(in)
3063     {
3064       fhGenMCAccE[kmcPi0]    ->Fill(eMC);
3065       if(eMC > 0.5) fhGenMCAccEtaPhi[kmcPi0]->Fill(etaMC,phiMC);
3066     }
3067   }
3068   else if (pdg==221)
3069   {
3070     fhGenMCE[kmcEta]    ->Fill(eMC);
3071     if(eMC > 0.5) fhGenMCEtaPhi[kmcEta]->Fill(etaMC,phiMC);
3072     if(in)
3073     {
3074       fhGenMCAccE[kmcEta]    ->Fill(eMC);
3075       if(eMC > 0.5) fhGenMCAccEtaPhi[kmcEta]->Fill(etaMC,phiMC);
3076     }
3077   }
3078   else if (TMath::Abs(pdg)==11)
3079   {
3080     fhGenMCE[kmcElectron]    ->Fill(eMC);
3081     if(eMC > 0.5) fhGenMCEtaPhi[kmcElectron]->Fill(etaMC,phiMC);
3082     if(in)
3083     {
3084       fhGenMCAccE[kmcElectron]    ->Fill(eMC);
3085       if(eMC > 0.5) fhGenMCAccEtaPhi[kmcElectron]->Fill(etaMC,phiMC);
3086     }
3087   }
3088 }
3089
3090 //_________________________________________________________________________________
3091 void AliAnaCalorimeterQA::WeightHistograms(AliVCluster *clus, AliVCaloCells* cells)
3092 {
3093   // Calculate weights
3094
3095   // First recalculate energy in case non linearity was applied
3096   Float_t  energy = 0;
3097   Float_t  ampMax = 0;
3098   for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
3099   {
3100     Int_t id       = clus->GetCellsAbsId()[ipos];
3101
3102     //Recalibrate cell energy if needed
3103     Float_t amp = cells->GetCellAmplitude(id);
3104     GetCaloUtils()->RecalibrateCellAmplitude(amp, fCalorimeter, id);
3105
3106     energy    += amp;
3107
3108     if(amp> ampMax)
3109       ampMax = amp;
3110
3111   } // energy loop
3112
3113   if(energy <=0 )
3114   {
3115     printf("AliAnaCalorimeterQA::WeightHistograms()- Wrong calculated energy %f\n",energy);
3116     return;
3117   }
3118
3119   fhEMaxCellClusterRatio   ->Fill(energy,ampMax/energy);
3120   fhEMaxCellClusterLogRatio->Fill(energy,TMath::Log(ampMax/energy));
3121
3122   //Get the ratio and log ratio to all cells in cluster
3123   for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
3124   {
3125     Int_t id       = clus->GetCellsAbsId()[ipos];
3126
3127     //Recalibrate cell energy if needed
3128     Float_t amp = cells->GetCellAmplitude(id);
3129     GetCaloUtils()->RecalibrateCellAmplitude(amp, fCalorimeter, id);
3130
3131     fhECellClusterRatio   ->Fill(energy,amp/energy);
3132     fhECellClusterLogRatio->Fill(energy,TMath::Log(amp/energy));
3133   }
3134
3135   //Recalculate shower shape for different W0
3136   if(fCalorimeter=="EMCAL")
3137   {
3138     Float_t l0org = clus->GetM02();
3139     Float_t l1org = clus->GetM20();
3140     Float_t dorg  = clus->GetDispersion();
3141
3142     for(Int_t iw = 0; iw < 14; iw++){
3143       GetCaloUtils()->GetEMCALRecoUtils()->SetW0(1+iw*0.5);
3144       GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), cells, clus);
3145
3146       fhLambda0ForW0[iw]->Fill(energy,clus->GetM02());
3147       //fhLambda1ForW0[iw]->Fill(energy,clus->GetM20());
3148
3149       if(IsDataMC()){
3150
3151         Int_t tag = GetMCAnalysisUtils()->CheckOrigin(clus->GetLabels(),clus->GetNLabels(), GetReader(),0);
3152
3153         if(   GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton) &&
3154            !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0)      &&
3155            !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta)      &&
3156            !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)        ){
3157           fhLambda0ForW0MC[iw][0]->Fill(energy,clus->GetM02());
3158           //fhLambda1ForW0MC[iw][0]->Fill(energy,clus->GetM20());
3159         } // Pure Photon
3160         else if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCElectron) &&
3161                 !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)        ){
3162           fhLambda0ForW0MC[iw][1]->Fill(energy,clus->GetM02());
3163           //fhLambda1ForW0MC[iw][1]->Fill(energy,clus->GetM20());
3164         } // Electron
3165         else if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)        ){
3166           fhLambda0ForW0MC[iw][2]->Fill(energy,clus->GetM02());
3167           //fhLambda1ForW0MC[iw][2]->Fill(energy,clus->GetM20());
3168         } // Conversion
3169         else if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) ){
3170           fhLambda0ForW0MC[iw][3]->Fill(energy,clus->GetM02());
3171           //fhLambda1ForW0MC[iw][3]->Fill(energy,clus->GetM20());
3172         }// Pi0
3173         else if(!GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta) &&
3174                 !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton) ){
3175           fhLambda0ForW0MC[iw][4]->Fill(energy,clus->GetM02());
3176           //fhLambda1ForW0MC[iw][4]->Fill(energy,clus->GetM20());
3177         }// Hadron
3178
3179       }// Is MC
3180     } // w0 loop
3181
3182     // Set the original values back
3183     clus->SetM02(l0org);
3184     clus->SetM20(l1org);
3185     clus->SetDispersion(dorg);
3186
3187   }// EMCAL
3188
3189 }
3190
3191
3192