First implementation of the time calibration for analysis in AliEMCALRecoUtils
authorgconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 27 Sep 2011 17:04:08 +0000 (17:04 +0000)
committergconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 27 Sep 2011 17:04:08 +0000 (17:04 +0000)
All the other classes addapt to this change
AliAnaCalorimeterQA: Some minor fixes, put back the track matching related histograms
AliCalTrackReader: Add a swich to recalculate the clusters or not, even if the AliEMCALRecoUtils switch are on, avoid multiple aplication of calibration parameters in case of using clusterization task.
AliAnalysisTaskEMCALClusterizer: Possibility to add matched tracks in AODCaloCluster added, this needs a patch in AliAODCaloCluster not committed yet.

EMCAL/AliEMCALRecoUtils.cxx
EMCAL/AliEMCALRecoUtils.h
PWG4/CaloCalib/AliAnalysisTaskCaloFilter.cxx
PWG4/CaloCalib/AliAnalysisTaskEMCALClusterize.cxx
PWG4/PartCorrBase/AliCaloTrackReader.cxx
PWG4/PartCorrBase/AliCaloTrackReader.h
PWG4/PartCorrDep/AliAnaCalorimeterQA.cxx
PWG4/PartCorrDep/AliAnaCalorimeterQA.h

index 8903777..7a22a37 100644 (file)
@@ -57,7 +57,7 @@
 #include "AliEMCALRecoUtils.h"
 #include "AliEMCALGeometry.h"
 #include "AliEMCALTrack.h"
-#include "AliEMCALCalibTimeDepCorrection.h"
+#include "AliEMCALCalibTimeDepCorrection.h" // Run dependent
 #include "AliEMCALPIDUtils.h"
 
 ClassImp(AliEMCALRecoUtils)
@@ -67,8 +67,9 @@ AliEMCALRecoUtils::AliEMCALRecoUtils():
   fParticleType(kPhoton),                 fPosAlgo(kUnchanged),                   fW0(4.), 
   fNonLinearityFunction(kNoCorrection),   fNonLinearThreshold(30),
   fSmearClusterEnergy(kFALSE),            fRandom(),
-  fRecalibration(kFALSE),                 fEMCALRecalibrationFactors(),
-  fUseTimeCorrectionFactors(kFALSE),      fTimeCorrectionFactorsSet(kFALSE),
+  fCellsRecalibrated(kFALSE),             fRecalibration(kFALSE),                 fEMCALRecalibrationFactors(),
+  fTimeRecalibration(kFALSE),             fEMCALTimeRecalibrationFactors(),
+  fUseRunCorrectionFactors(kFALSE),       fRunCorrectionFactorsSet(kFALSE),
   fRemoveBadChannels(kFALSE),             fRecalDistToBadChannels(kFALSE),        fEMCALBadChannelMap(),
   fNCellsFromEMCALBorder(0),              fNoEMCALBorderAtEta0(kTRUE),
   fRejectExoticCluster(kFALSE),           fPIDUtils(),                            fAODFilterMask(32),
@@ -132,8 +133,10 @@ AliEMCALRecoUtils::AliEMCALRecoUtils(const AliEMCALRecoUtils & reco)
   fParticleType(reco.fParticleType),                         fPosAlgo(reco.fPosAlgo),     fW0(reco.fW0),
   fNonLinearityFunction(reco.fNonLinearityFunction),         fNonLinearThreshold(reco.fNonLinearThreshold),
   fSmearClusterEnergy(reco.fSmearClusterEnergy),             fRandom(),
+  fCellsRecalibrated(reco.fCellsRecalibrated),
   fRecalibration(reco.fRecalibration),                       fEMCALRecalibrationFactors(reco.fEMCALRecalibrationFactors),
-  fUseTimeCorrectionFactors(reco.fUseTimeCorrectionFactors), fTimeCorrectionFactorsSet(reco.fTimeCorrectionFactorsSet),
+  fTimeRecalibration(reco.fTimeRecalibration),               fEMCALTimeRecalibrationFactors(reco.fEMCALTimeRecalibrationFactors),
+  fUseRunCorrectionFactors(reco.fUseRunCorrectionFactors),   fRunCorrectionFactorsSet(reco.fRunCorrectionFactorsSet),
   fRemoveBadChannels(reco.fRemoveBadChannels),               fRecalDistToBadChannels(reco.fRecalDistToBadChannels),
   fEMCALBadChannelMap(reco.fEMCALBadChannelMap),
   fNCellsFromEMCALBorder(reco.fNCellsFromEMCALBorder),       fNoEMCALBorderAtEta0(reco.fNoEMCALBorderAtEta0),
@@ -184,10 +187,15 @@ AliEMCALRecoUtils & AliEMCALRecoUtils::operator = (const AliEMCALRecoUtils & rec
   fNonLinearThreshold        = reco.fNonLinearThreshold;
   fSmearClusterEnergy        = reco.fSmearClusterEnergy;
 
+  fCellsRecalibrated         = reco.fCellsRecalibrated;
   fRecalibration             = reco.fRecalibration;
   fEMCALRecalibrationFactors = reco.fEMCALRecalibrationFactors;
-  fUseTimeCorrectionFactors  = reco.fUseTimeCorrectionFactors;
-  fTimeCorrectionFactorsSet  = reco.fTimeCorrectionFactorsSet;
+
+  fTimeRecalibration             = reco.fTimeRecalibration;
+  fEMCALTimeRecalibrationFactors = reco.fEMCALTimeRecalibrationFactors;
+
+  fUseRunCorrectionFactors   = reco.fUseRunCorrectionFactors;
+  fRunCorrectionFactorsSet   = reco.fRunCorrectionFactorsSet;
   
   fRemoveBadChannels         = reco.fRemoveBadChannels;
   fRecalDistToBadChannels    = reco.fRecalDistToBadChannels;
@@ -285,6 +293,11 @@ AliEMCALRecoUtils::~AliEMCALRecoUtils()
                delete  fEMCALRecalibrationFactors;
        }       
   
+  if(fEMCALTimeRecalibrationFactors) { 
+               fEMCALTimeRecalibrationFactors->Clear();
+               delete  fEMCALTimeRecalibrationFactors;
+       }       
+  
   if(fEMCALBadChannelMap) { 
                fEMCALBadChannelMap->Clear();
                delete  fEMCALBadChannelMap;
@@ -716,6 +729,31 @@ void AliEMCALRecoUtils::InitEMCALRecalibrationFactors(){
        TH1::AddDirectory(oldStatus);           
 }
 
+//________________________________________________________________
+void AliEMCALRecoUtils::InitEMCALTimeRecalibrationFactors(){
+       //Init EMCAL recalibration factors
+       AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
+       //In order to avoid rewriting the same histograms
+       Bool_t oldStatus = TH1::AddDirectoryStatus();
+       TH1::AddDirectory(kFALSE);
+  
+       fEMCALTimeRecalibrationFactors = new TObjArray(4);
+       for (int i = 0; i < 4; i++) 
+    fEMCALTimeRecalibrationFactors->Add(new TH1F(Form("hAllTimeAvBC%d",i),
+                                                 Form("hAllTimeAvBC%d",i),  
+                                                 11521,0.,11521)          );
+       //Init the histograms with 1
+       for (Int_t bc = 0; bc < 4; bc++) {
+               
+      SetEMCALChannelTimeRecalibrationFactor(bc,0.);
+  }
+
+       fEMCALTimeRecalibrationFactors->SetOwner(kTRUE);
+       fEMCALTimeRecalibrationFactors->Compress();
+       
+       //In order to avoid rewriting the same histograms
+       TH1::AddDirectory(oldStatus);           
+}
 
 //________________________________________________________________
 void AliEMCALRecoUtils::InitEMCALBadChannelStatusMap(){
@@ -730,9 +768,7 @@ void AliEMCALRecoUtils::InitEMCALBadChannelStatusMap(){
        for (int i = 0; i < 10; i++) {
                fEMCALBadChannelMap->Add(new TH2I(Form("EMCALBadChannelMap_Mod%d",i),Form("EMCALBadChannelMap_Mod%d",i), 48, 0, 48, 24, 0, 24));
        }
-       
-       //delete hTemp;
-       
+               
        fEMCALBadChannelMap->SetOwner(kTRUE);
        fEMCALBadChannelMap->Compress();
        
@@ -741,8 +777,10 @@ void AliEMCALRecoUtils::InitEMCALBadChannelStatusMap(){
 }
 
 //________________________________________________________________
-void AliEMCALRecoUtils::RecalibrateClusterEnergy(AliEMCALGeometry* geom, AliVCluster * cluster, AliVCaloCells * cells){
-       // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
+void AliEMCALRecoUtils::RecalibrateClusterEnergy(AliEMCALGeometry* geom, AliVCluster * cluster, AliVCaloCells * cells, const Int_t bc){
+       // Recalibrate the cluster energy and Time, considering the recalibration map 
+  // and the energy of the cells and time that compose the cluster.
+  // bc= bunch crossing number returned by esdevent->GetBunchCrossNumber();
        
   if(!cluster){
     AliInfo("Cluster pointer null!");
@@ -756,34 +794,160 @@ void AliEMCALRecoUtils::RecalibrateClusterEnergy(AliEMCALGeometry* geom, AliVClu
        
        //Initialize some used variables
        Float_t energy = 0;
-       Int_t absId    = -1;
-  Int_t icol = -1, irow = -1, imod=1;
+       Int_t   absId  =-1;
+  Int_t   icol   =-1, irow =-1, imod=1;
        Float_t factor = 1, frac = 0;
-       
+  Int_t   absIdMax = -1;
+  Float_t emax     = 0;
+  
        //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
        for(Int_t icell = 0; icell < ncells; icell++){
                absId = index[icell];
                frac =  fraction[icell];
                if(frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off
-               Int_t iTower = -1, iIphi = -1, iIeta = -1; 
-               geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta); 
-               if(fEMCALRecalibrationFactors->GetEntries() <= imod) continue;
-               geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);                  
-               factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
-    AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n",
-             imod,icol,irow,frac,factor,cells->GetCellAmplitude(absId)));
-               
+    
+    if(!fCellsRecalibrated && IsRecalibrationOn()){
+      
+      // Energy  
+      Int_t iTower = -1, iIphi = -1, iIeta = -1; 
+      geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta); 
+      if(fEMCALRecalibrationFactors->GetEntries() <= imod) continue;
+      geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);                   
+      factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
+      
+      AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n",
+                      imod,icol,irow,frac,factor,cells->GetCellAmplitude(absId)));
+      
+    } 
+    
                energy += cells->GetCellAmplitude(absId)*factor*frac;
+    
+    if(emax < cells->GetCellAmplitude(absId)*factor*frac){
+      emax     = cells->GetCellAmplitude(absId)*factor*frac;
+      absIdMax = absId;
+    }
+
        }
        
-       
-               AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Energy before %f, after %f\n",cluster->E(),energy));
-       
-       cluster->SetE(energy);
-       
+  cluster->SetE(energy);
+
+  AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Energy before %f, after %f\n",cluster->E(),energy));
+
+       // Recalculate time of cluster only for ESDs
+  if(!strcmp("AliESDCaloCluster",Form("%s",cluster->ClassName()))){
+    
+    // Time
+    Double_t weightedTime = 0;
+    Double_t weight       = 0;
+    Double_t weightTot    = 0;
+    Double_t maxcellTime  = 0;
+    for(Int_t icell = 0; icell < ncells; icell++){
+      absId = index[icell];
+      frac =  fraction[icell];
+      if(frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off
+
+      Double_t celltime = cells->GetCellTime(absId);
+      RecalibrateCellTime(absId, bc, celltime);
+      if(absId == absIdMax) maxcellTime = celltime;
+
+      if(!fCellsRecalibrated){
+      
+        Int_t iTower = -1, iIphi = -1, iIeta = -1; 
+        geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta); 
+        if(fEMCALRecalibrationFactors->GetEntries() <= imod) continue;
+        geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);                 
+        factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
+        
+        AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n",
+                        imod,icol,irow,frac,factor,cells->GetCellTime(absId)));
+        
+      } 
+      
+      weight        = GetCellWeight(cells->GetCellAmplitude(absId)*factor*frac , energy );
+      weightTot    += weight;
+      weightedTime += celltime * weight;
+      
+    }
+    
+    if(weightTot > 0)
+      cluster->SetTOF(weightedTime/weightTot);
+    else 
+      cluster->SetTOF(maxcellTime);
+    
+  }
+}
+
+//________________________________________________________________
+void AliEMCALRecoUtils::RecalibrateCells(AliEMCALGeometry* geom, AliVCaloCells * cells, Int_t bc){
+       // Recalibrate the cells time and energy, considering the recalibration map and the energy 
+  // of the cells that compose the cluster.
+  // bc= bunch crossing number returned by esdevent->GetBunchCrossNumber();
+
+  if(!IsRecalibrationOn()) return;
+  
+  if(!cells){
+    AliInfo("Cells pointer null!");
+    return;
+  }  
+  
+  fCellsRecalibrated = kTRUE;
+  
+  Int_t absId  =-1;
+  Int_t icol   =-1, irow  =-1, imod  = 1;
+  Int_t iTower =-1, iIeta =-1, iIphi =-1;
+
+  Int_t nEMcell = cells->GetNumberOfCells() ;
+  
+  for (Int_t iCell = 0; iCell < nEMcell; iCell++) { 
+    
+    absId = cells->GetCellNumber(iCell);
+    
+    // Energy
+    Float_t factor = 1;
+    if(IsRecalibrationOn()){
+      geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta); 
+      if(fEMCALRecalibrationFactors->GetEntries() <= imod) continue;
+      geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);   
+      factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
+               }
+    
+    Float_t cellE      = cells->GetAmplitude(iCell) * factor ;
+    
+    //Time
+    Double_t celltime = cells->GetCellTime(absId);
+    RecalibrateCellTime(absId, bc, celltime);
+    
+    //Set new values
+    cells->SetCell(iCell,cells->GetCellNumber(iCell),cellE, celltime);
+    
+  }
+  
 }
 
+//________________________________________________________________
+void AliEMCALRecoUtils::RecalibrateCellTime(const Int_t absId, const Int_t bc, Double_t & celltime){
+       // Recalibrate time of cell with absID  considering the recalibration map 
+  // bc= bunch crossing number returned by esdevent->GetBunchCrossNumber();
+    
+  if(!fCellsRecalibrated && IsTimeRecalibrationOn()){
+//    printf("cell time org %g, ",celltime);
 
+    Double_t timeBCoffset = 0.;
+    if( bc%4 ==0 || bc%4==1) timeBCoffset = 100.*1.e-9; //in ns        
+    
+    Double_t celloffset = GetEMCALChannelTimeRecalibrationFactor(bc%4,absId)*1.e-9; 
+    
+//    printf("absId %d, time %f bc %d-%d: bc0 %f, bc1 %f, bc2 %f, bc3 %f \n", absId, celltime*1.e9,bc, bc%4, 
+//           GetEMCALChannelTimeRecalibrationFactor(0,absId),GetEMCALChannelTimeRecalibrationFactor(1,absId),
+//           GetEMCALChannelTimeRecalibrationFactor(2,absId),GetEMCALChannelTimeRecalibrationFactor(3,absId));
+    
+    celltime -= timeBCoffset ;
+    celltime -= celloffset   ;  
+//    printf("new %g\n",celltime);
+  }
+  
+}
+  
 //__________________________________________________
 void AliEMCALRecoUtils::RecalculateClusterPosition(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
 {
@@ -825,20 +989,26 @@ void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerGlobal(AliEMCALGeomet
   //printf("** Cluster energy %f, ncells %d, depth %f\n",clEnergy,clu->GetNCells(),depth);
   
   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
-    absId = clu->GetCellAbsId(iDig);
-    fraction  = clu->GetCellAmplitudeFraction(iDig);
-    if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
-    geom->GetCellIndex(absId,iSM,iTower,iIphi,iIeta); 
-    geom->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);                      
     
-    if(IsRecalibrationOn()) {
-      recalFactor = GetEMCALChannelRecalibrationFactor(iSM,ieta,iphi);
+      absId = clu->GetCellAbsId(iDig);
+      fraction  = clu->GetCellAmplitudeFraction(iDig);
+      if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
+    
+    if(!fCellsRecalibrated){
+
+      geom->GetCellIndex(absId,iSM,iTower,iIphi,iIeta); 
+      geom->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);                    
+      
+      if(IsRecalibrationOn()) {
+        recalFactor = GetEMCALChannelRecalibrationFactor(iSM,ieta,iphi);
+      }
     }
+    
     eCell  = cells->GetCellAmplitude(absId)*fraction*recalFactor;
     
     weight = GetCellWeight(eCell,clEnergy);
-    //printf("cell energy %f, weight %f\n",eCell,weight);
     totalWeight += weight;
+    
     geom->RelPosCellInSModule(absId,depth,pLocal[0],pLocal[1],pLocal[2]);
     //printf("pLocal (%f,%f,%f), SM %d, absId %d\n",pLocal[0],pLocal[1],pLocal[2],iSupModMax,absId);
     geom->GetGlobal(pLocal,pGlobal,iSupModMax);
@@ -906,17 +1076,24 @@ void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerIndex(AliEMCALGeometr
     absId = clu->GetCellAbsId(iDig);
     fraction  = clu->GetCellAmplitudeFraction(iDig);
     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
-    geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta); 
-    geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);                  
-    
+
     if     (iDig==0)  startingSM = iSupMod;
     else if(iSupMod != startingSM) areInSameSM = kFALSE;
 
     eCell  = cells->GetCellAmplitude(absId);
     
-    if(IsRecalibrationOn()) {
-      recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
+    geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta); 
+    geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);          
+    
+    if(!fCellsRecalibrated){
+      
+      if(IsRecalibrationOn()) {
+                
+        recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
+        
+      }
     }
+    
     eCell  = cells->GetCellAmplitude(absId)*fraction*recalFactor;
     
     weight = GetCellWeight(eCell,clEnergy);
@@ -1091,9 +1268,15 @@ void AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(AliEMCALGeometry
     //Get the cell energy, if recalibration is on, apply factors
     fraction  = cluster->GetCellAmplitudeFraction(iDigit);
     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
-    if(IsRecalibrationOn()) {
-      recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
+    
+    if(!fCellsRecalibrated){
+      
+      if(IsRecalibrationOn()) {
+        recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
+      }
+      
     }
+    
     eCell  = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
     
     if(cluster->E() > 0 && eCell > 0){
@@ -1751,14 +1934,14 @@ void AliEMCALRecoUtils::Print(const Option_t *) const
 }
 
 //_____________________________________________________________________
-void AliEMCALRecoUtils::SetTimeDependentCorrections(Int_t runnumber){
+void AliEMCALRecoUtils::SetRunDependentCorrections(Int_t runnumber){
   //Get EMCAL time dependent corrections from file and put them in the recalibration histograms
   //Do it only once and only if it is requested
   
-  if(!fUseTimeCorrectionFactors) return;
-  if(fTimeCorrectionFactorsSet)  return;
+  if(!fUseRunCorrectionFactors) return;
+  if(fRunCorrectionFactorsSet)  return;
   
-  printf("AliEMCALRecoUtils::GetTimeDependentCorrections() - Get Correction Factors for Run number %d\n",runnumber);
+  AliInfo(Form("AliEMCALRecoUtils::GetRunDependentCorrections() - Get Correction Factors for Run number %d\n",runnumber));
  
   AliEMCALCalibTimeDepCorrection  *corr =  new AliEMCALCalibTimeDepCorrection();
   corr->ReadRootInfo(Form("CorrectionFiles/Run%d_Correction.root",runnumber));
@@ -1776,6 +1959,6 @@ void AliEMCALRecoUtils::SetTimeDependentCorrections(Int_t runnumber){
       }
     }
   }
-   fTimeCorrectionFactorsSet = kTRUE;
+   fRunCorrectionFactorsSet = kTRUE;
 }
 
index f2292ef..21b2cb2 100644 (file)
@@ -129,10 +129,11 @@ public:
                                                            else     { AliInfo(Form("Index %d larger than 2, do nothing\n",i)) ; } }
   
   //-----------------------------------------------------
-  //Recalibration
+  // Energy Recalibration
   //-----------------------------------------------------
-
-  void RecalibrateClusterEnergy(AliEMCALGeometry* geom, AliVCluster* cluster, AliVCaloCells * cells) ;
+  
+  void     RecalibrateCells(AliEMCALGeometry* geom, AliVCaloCells * cells, Int_t bc) ; // Energy and Time
+  void     RecalibrateClusterEnergy(AliEMCALGeometry* geom, AliVCluster* cluster, AliVCaloCells * cells, const Int_t bc=0) ; // Energy and time
 
   Bool_t   IsRecalibrationOn()                     const { return fRecalibration ; }
   void     SwitchOffRecalibration()                      { fRecalibration = kFALSE ; }
@@ -140,27 +141,52 @@ public:
                                                            if(!fEMCALRecalibrationFactors)InitEMCALRecalibrationFactors() ; }
   void     InitEMCALRecalibrationFactors() ;
 
-  //Recalibrate channels with time dependent corrections
-  void     SwitchOffTimeDepCorrection()                  { fUseTimeCorrectionFactors = kFALSE ; }
-  void     SwitchOnTimeDepCorrection()                   { fUseTimeCorrectionFactors = kTRUE  ; 
-                                                           SwitchOnRecalibration()            ; }
-  void     SetTimeDependentCorrections(Int_t runnumber);
-    
+  TH2F *   GetEMCALChannelRecalibrationFactors(Int_t iSM)     const { return (TH2F*)fEMCALRecalibrationFactors->At(iSM) ; }    
+  void     SetEMCALChannelRecalibrationFactors(TObjArray *map)      { fEMCALRecalibrationFactors = map                  ; }
+  void     SetEMCALChannelRecalibrationFactors(Int_t iSM , TH2F* h) { fEMCALRecalibrationFactors->AddAt(h,iSM)          ; }
+  
   Float_t  GetEMCALChannelRecalibrationFactor(Int_t iSM , Int_t iCol, Int_t iRow) const { 
-                                                           if(fEMCALRecalibrationFactors) 
-                                                             return (Float_t) ((TH2F*)fEMCALRecalibrationFactors->At(iSM))->GetBinContent(iCol,iRow); 
-                                                           else return 1 ; } 
+    if(fEMCALRecalibrationFactors) 
+      return (Float_t) ((TH2F*)fEMCALRecalibrationFactors->At(iSM))->GetBinContent(iCol,iRow); 
+    else return 1 ; } 
        
   void     SetEMCALChannelRecalibrationFactor(Int_t iSM , Int_t iCol, Int_t iRow, Double_t c = 1) { 
-                                                           if(!fEMCALRecalibrationFactors) InitEMCALRecalibrationFactors() ;
-                                                           ((TH2F*)fEMCALRecalibrationFactors->At(iSM))->SetBinContent(iCol,iRow,c) ; }  
+    if(!fEMCALRecalibrationFactors) InitEMCALRecalibrationFactors() ;
+    ((TH2F*)fEMCALRecalibrationFactors->At(iSM))->SetBinContent(iCol,iRow,c) ; }
+  
+  //Recalibrate channels energy with run dependent corrections
+  void     SwitchOffRunDepCorrection()                   { fUseRunCorrectionFactors = kFALSE ; }
+  void     SwitchOnRunDepCorrection()                    { fUseRunCorrectionFactors = kTRUE  ; 
+                                                           SwitchOnRecalibration()           ; }
+  void     SetRunDependentCorrections(Int_t runnumber);
+      
+  //-----------------------------------------------------
+  // Time Recalibration
+  //-----------------------------------------------------
+  
+  void     RecalibrateCellTime(const Int_t absId, const Int_t bc, Double_t & time);
+  
+  Bool_t   IsTimeRecalibrationOn()                 const { return fTimeRecalibration   ; }
+  void     SwitchOffTimeRecalibration()                  { fTimeRecalibration = kFALSE ; }
+  void     SwitchOnTimeRecalibration()                   { fTimeRecalibration = kTRUE  ; 
+    if(!fEMCALTimeRecalibrationFactors)InitEMCALTimeRecalibrationFactors() ; }
+  void     InitEMCALTimeRecalibrationFactors() ;
+  
+  Float_t  GetEMCALChannelTimeRecalibrationFactor(Int_t bc, Int_t absID) const { 
+    if(fEMCALTimeRecalibrationFactors) 
+      return (Float_t) ((TH1F*)fEMCALTimeRecalibrationFactors->At(bc))->GetBinContent(absID); 
+    else return 1 ; } 
+       
+  void     SetEMCALChannelTimeRecalibrationFactor(Int_t bc,Int_t absID, Double_t c = 1) { 
+    if(!fEMCALTimeRecalibrationFactors) InitEMCALTimeRecalibrationFactors() ;
+    ((TH1F*)fEMCALTimeRecalibrationFactors->At(bc))->SetBinContent(absID,c) ; }  
+  
+  TH1F *   GetEMCALChannelTimeRecalibrationFactors(Int_t bc)      const { return (TH1F*)fEMCALTimeRecalibrationFactors->At(bc) ; }     
+  void     SetEMCALChannelTimeRecalibrationFactors(TObjArray *map)      { fEMCALTimeRecalibrationFactors = map                 ; }
+  void     SetEMCALChannelTimeRecalibrationFactors(Int_t bc , TH1F* h)  { fEMCALTimeRecalibrationFactors->AddAt(h,bc)          ; }
   
-  TH2F *   GetEMCALChannelRecalibrationFactors(Int_t iSM)     const { return (TH2F*)fEMCALRecalibrationFactors->At(iSM) ; }    
-  void     SetEMCALChannelRecalibrationFactors(TObjArray *map)      { fEMCALRecalibrationFactors = map                  ; }
-  void     SetEMCALChannelRecalibrationFactors(Int_t iSM , TH2F* h) { fEMCALRecalibrationFactors->AddAt(h,iSM)          ; }
-
   //-----------------------------------------------------
-  //Modules fiducial region, remove clusters in borders
+  // Modules fiducial region, remove clusters in borders
   //-----------------------------------------------------
 
   Bool_t   CheckCellFiducialRegion(AliEMCALGeometry* geom, AliVCluster* cluster, AliVCaloCells* cells) ;
@@ -312,13 +338,18 @@ private:
   Float_t    fSmearClusterParam[3];      // Smearing parameters
   TRandom3   fRandom;                    // Random generator
     
-  // Recalibration 
+  // Energy Recalibration 
+  Bool_t     fCellsRecalibrated;         // Internal bool to check if cells (time/energy) where recalibrated and not recalibrate them when recalculating different things
   Bool_t     fRecalibration;             // Switch on or off the recalibration
   TObjArray* fEMCALRecalibrationFactors; // Array of histograms with map of recalibration factors, EMCAL
     
-  // Recalibrate with run dependent corrections
-  Bool_t     fUseTimeCorrectionFactors;  // Use Time Dependent Correction
-  Bool_t     fTimeCorrectionFactorsSet;  // Time Correction set at leat once
+  // Time Recalibration 
+  Bool_t     fTimeRecalibration;             // Switch on or off the time recalibration
+  TObjArray* fEMCALTimeRecalibrationFactors; // Array of histograms with map of time recalibration factors, EMCAL
+  
+  // Recalibrate with run dependent corrections, energy
+  Bool_t     fUseRunCorrectionFactors;   // Use Run Dependent Correction
+  Bool_t     fRunCorrectionFactorsSet;   // Run Correction set at leat once
     
   // Bad Channels
   Bool_t     fRemoveBadChannels;         // Check the channel status provided and remove clusters with bad channels
@@ -363,7 +394,7 @@ private:
   Float_t    fCutMaxDCAToVertexZ;        // Track-to-vertex cut in max absolute distance in z-plane
   Bool_t     fCutDCAToVertex2D;          // If true a 2D DCA cut is made. Tracks are accepted if sqrt((DCAXY / fCutMaxDCAToVertexXY)^2 + (DCAZ / fCutMaxDCAToVertexZ)^2) < 1 AND sqrt((DCAXY / fCutMinDCAToVertexXY)^2 + (DCAZ / fCutMinDCAToVertexZ)^2) > 1
   
-  ClassDef(AliEMCALRecoUtils, 13)
+  ClassDef(AliEMCALRecoUtils, 14)
   
 };
 
index fbccca2..1b2f2b6 100644 (file)
@@ -344,7 +344,7 @@ void AliAnalysisTaskCaloFilter::UserExec(Option_t */*option*/)
       }//Load matrices from Data 
       
       //Recover time dependent corrections, put then in recalibration histograms. Do it once
-      fEMCALRecoUtils->SetTimeDependentCorrections(InputEvent()->GetRunNumber());
+      fEMCALRecoUtils->SetRunDependentCorrections(InputEvent()->GetRunNumber());
 
     }//first event
     
index ed234ae..f1d5c1e 100644 (file)
@@ -459,7 +459,7 @@ void AliAnalysisTaskEMCALClusterize::UserExec(Option_t *)
     }//Load matrices from Data 
     
     //Recover time dependent corrections, put then in recalibration histograms. Do it once
-    fRecoUtils->SetTimeDependentCorrections(InputEvent()->GetRunNumber());
+    fRecoUtils->SetRunDependentCorrections(InputEvent()->GetRunNumber());
     
   }//first event
   
@@ -583,6 +583,8 @@ void AliAnalysisTaskEMCALClusterize::UserExec(Option_t *)
         for(Int_t icell=0; icell < clus->GetNCells(); icell++ ){
           fCellLabels[index[icell]]       = label;
           fCellSecondLabels[index[icell]] = label2;
+          //printf("Clusterizer in : TOF %g\n",clus->GetTOF()*1.e9);
+
           fCellTime[icell]                = clus->GetTOF();        
 
           //printf("1) absID %d, label[0] %d label[1] %d\n",index[icell], fCellLabels[index[icell]],fCellSecondLabels[index[icell]]);
@@ -593,14 +595,16 @@ void AliAnalysisTaskEMCALClusterize::UserExec(Option_t *)
     
     // Create digits 
     //.................
-    Int_t idigit =  0;
-    Int_t id     = -1;
-    Float_t amp  = -1; 
-    Float_t time = -1; 
+    Int_t    idigit =  0;
+    Int_t    id     = -1;
+    Float_t  amp    = -1; 
+    Double_t time   = -1; 
     
     TTree *digitsTree = new TTree("digitstree","digitstree");
     digitsTree->Branch("EMCAL","TClonesArray", &fDigitsArr, 32000);
     
+    Int_t bc = InputEvent()->GetBunchCrossNumber();
+    
     for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
     {
       if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
@@ -629,10 +633,17 @@ void AliAnalysisTaskEMCALClusterize::UserExec(Option_t *)
         //printf("CalibFactor %f times %f for id %d\n",fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi),amp,id);
         amp *=fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
       }
-      
+            
       // In case of AOD analysis cell time is 0, approximate replacing by time of the cluster the digit belongs.
       if (time*1e9 < 1.) time = fCellTime[id];
       
+      // Recalibrate time
+      fRecoUtils->RecalibrateCellTime(id,bc,time);
+      
+//      printf("Clusterizer: Id %d, Time org %e, Time new %e; Amp org %f, Amp new %f\n",
+//             id, cells->GetTime(icell),time, cells->GetAmplitude(icell),amp);
+
+      
       //Create the digit, put a fake primary deposited energy to trick the clusterizer when checking the most likely primary
       new((*fDigitsArr)[idigit]) AliEMCALDigit( fCellLabels[id], fCellLabels[id],id, amp, time,AliEMCALDigit::kHG,idigit, 0, 0, 1); 
       //if(fCellLabels[id]>=0)printf("2) Digit cell %d, label %d\n",id,fCellLabels[id]) ;
@@ -705,8 +716,8 @@ void AliAnalysisTaskEMCALClusterize::UserExec(Option_t *)
     AliAODCaloCluster *newCluster = (AliAODCaloCluster *) fCaloClusterArr->At(i);
     //if(Entry()==0) Info("UserExec","newCluster E %f\n", newCluster->E());
     
-    //Add matched track, if any, only with ESDs
-    if(esdevent && fDoTrackMatching){
+    //Add matched track
+    if(fDoTrackMatching){
       Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
       if(trackIndex >= 0){
         newCluster->AddTrackMatched(event->GetTrack(trackIndex));
index 4a3beeb..045a21a 100755 (executable)
@@ -51,25 +51,31 @@ ClassImp(AliCaloTrackReader)
   
 //____________________________________________________________________________
   AliCaloTrackReader::AliCaloTrackReader() : 
-    TObject(), fEventNumber(-1), //fCurrentFileName(""),
-    fDataType(0), fDebug(0), 
-    fFiducialCut(0x0), fCheckFidCut(kFALSE), fComparePtHardAndJetPt(kFALSE), fPtHardAndJetPtFactor(7),
-    fCTSPtMin(0), fEMCALPtMin(0),fPHOSPtMin(0), 
-    fCTSPtMax(1000), fEMCALPtMax(1000),fPHOSPtMax(1000), fAODBranchList(new TList ),
+    TObject(),                   fEventNumber(-1), //fCurrentFileName(""),
+    fDataType(0),                fDebug(0), 
+    fFiducialCut(0x0),           fCheckFidCut(kFALSE), 
+    fComparePtHardAndJetPt(0),   fPtHardAndJetPtFactor(7),
+    fCTSPtMin(0),                fEMCALPtMin(0),                  fPHOSPtMin(0), 
+    fCTSPtMax(1000),             fEMCALPtMax(1000),               fPHOSPtMax(1000), 
+    fAODBranchList(new TList ),
     fCTSTracks(new TObjArray()), fEMCALClusters(new TObjArray()), fPHOSClusters(new TObjArray()),
-    fEMCALCells(0x0), fPHOSCells(0x0),
-    fInputEvent(0x0), fOutputEvent(0x0),fMC(0x0),
-    fFillCTS(0),fFillEMCAL(0),fFillPHOS(0),
-    fFillEMCALCells(0),fFillPHOSCells(0),  fSelectEmbeddedClusters(kFALSE),
-    fTrackStatus(0), fTrackFilterMask(0), fESDtrackCuts(0), fTrackMult(0), fTrackMultEtaCut(0.8),
-    fReadStack(kFALSE), fReadAODMCParticles(kFALSE), 
-    fDeltaAODFileName("deltaAODPartCorr.root"),fFiredTriggerClassName(""),
-    fAnaLED(kFALSE),fTaskName(""),fCaloUtils(0x0), 
-    fMixedEvent(NULL), fNMixedEvent(1), fVertex(NULL), 
-    fWriteOutputDeltaAOD(kFALSE),fOldAOD(kFALSE),fCaloFilterPatch(kFALSE),
-    fEMCALClustersListName(""),fZvtxCut(0.), 
-    fDoEventSelection(kFALSE),   fDoV0ANDEventSelection(kFALSE), fUseEventsWithPrimaryVertex(kFALSE),
-    fTriggerAnalysis (new AliTriggerAnalysis), fCentralityClass("V0M"),fCentralityOpt(10),
+    fEMCALCells(0x0),            fPHOSCells(0x0),
+    fInputEvent(0x0),            fOutputEvent(0x0),fMC(0x0),
+    fFillCTS(0),                 fFillEMCAL(0),                   fFillPHOS(0),
+    fFillEMCALCells(0),          fFillPHOSCells(0), 
+    fRecalculateClusters(kFALSE),fSelectEmbeddedClusters(kFALSE),
+    fTrackStatus(0),             fTrackFilterMask(0),             fESDtrackCuts(0), 
+    fTrackMult(0),               fTrackMultEtaCut(0.8),
+    fReadStack(kFALSE),          fReadAODMCParticles(kFALSE), 
+    fDeltaAODFileName("deltaAODPartCorr.root"),
+    fFiredTriggerClassName(""),  fAnaLED(kFALSE),
+    fTaskName(""),               fCaloUtils(0x0), 
+    fMixedEvent(NULL),           fNMixedEvent(1),                 fVertex(NULL), 
+    fWriteOutputDeltaAOD(kFALSE),fOldAOD(kFALSE),                 fCaloFilterPatch(kFALSE),
+    fEMCALClustersListName(""),  fZvtxCut(0.), 
+    fDoEventSelection(kFALSE),   fDoV0ANDEventSelection(kFALSE),  fUseEventsWithPrimaryVertex(kFALSE),
+    fTriggerAnalysis (new AliTriggerAnalysis), 
+    fCentralityClass("V0M"),     fCentralityOpt(10),
     fEventPlaneMethod("Q")
    
 {
@@ -322,7 +328,8 @@ void AliCaloTrackReader::Print(const Option_t * opt) const
   printf("Track filter mask (AODs) =  %d\n", (Int_t) fTrackFilterMask) ;
   printf("Track Mult Eta Cut =  %d\n", (Int_t) fTrackMultEtaCut) ;
   printf("Write delta AOD =     %d\n", fWriteOutputDeltaAOD) ;
-
+  printf("Recalculate Clusters = %d\n", fRecalculateClusters) ;
+  
   if(fComparePtHardAndJetPt)
          printf("Compare jet pt and pt hard to accept event, factor = %2.2f",fPtHardAndJetPtFactor);
                
@@ -681,12 +688,12 @@ void AliCaloTrackReader::FillInputEMCALAlgorithm(AliVCluster * clus, const Int_t
   
   //Reject clusters with bad channels, close to borders and exotic;
   if(!GetCaloUtils()->GetEMCALRecoUtils()->IsGoodCluster(clus,GetCaloUtils()->GetEMCALGeometry(),GetEMCALCells())) return;
-//  //Check if the cluster contains any bad channel and if close to calorimeter borders
-//  if(GetCaloUtils()->ClusterContainsBadChannel("EMCAL",clus->GetCellsAbsId(), clus->GetNCells())) 
-//    return;
-//  if(!GetCaloUtils()->CheckCellFiducialRegion(clus, (AliVCaloCells*)fInputEvent->GetEMCALCells(), fInputEvent, vindex)) 
-//    return;
-//  
+  //  //Check if the cluster contains any bad channel and if close to calorimeter borders
+  //  if(GetCaloUtils()->ClusterContainsBadChannel("EMCAL",clus->GetCellsAbsId(), clus->GetNCells())) 
+  //    return;
+  //  if(!GetCaloUtils()->CheckCellFiducialRegion(clus, (AliVCaloCells*)fInputEvent->GetEMCALCells(), fInputEvent, vindex)) 
+  //    return;
+  //  
   
   //Mask all cells in collumns facing ALICE thick material if requested
   if(GetCaloUtils()->GetNMaskCellColumns()){
@@ -721,35 +728,37 @@ void AliCaloTrackReader::FillInputEMCALAlgorithm(AliVCluster * clus, const Int_t
     //clus->GetPosition(pos);
     //printf("Before Corrections: e %f, x %f, y %f, z %f\n",clus->E(),pos[0],pos[1],pos[2]);
     
-    //Recalibrate the cluster energy 
-    if(GetCaloUtils()->IsRecalibrationOn()) {
-      Float_t energy = GetCaloUtils()->RecalibrateClusterEnergy(clus, GetEMCALCells());
-      clus->SetE(energy);
-      //printf("Recalibrated Energy %f\n",clus->E());  
-      GetCaloUtils()->RecalculateClusterShowerShapeParameters(GetEMCALCells(),clus);
-      GetCaloUtils()->RecalculateClusterPID(clus);
-    }
-    
-    //Recalculate distance to bad channels, if new list of bad channels provided
-    GetCaloUtils()->RecalculateClusterDistanceToBadChannel(GetEMCALCells(),clus);
-    
-    //Recalculate cluster position
-    if(GetCaloUtils()->IsRecalculationOfClusterPositionOn()){
-      GetCaloUtils()->RecalculateClusterPosition(GetEMCALCells(),clus); 
-      //clus->GetPosition(pos);
-      //printf("After  Corrections: e %f, x %f, y %f, z %f\n",clus->E(),pos[0],pos[1],pos[2]);
+    if(fRecalculateClusters){
+      //Recalibrate the cluster energy 
+      if(GetCaloUtils()->IsRecalibrationOn()) {
+        Float_t energy = GetCaloUtils()->RecalibrateClusterEnergy(clus, GetEMCALCells());
+        clus->SetE(energy);
+        //printf("Recalibrated Energy %f\n",clus->E());  
+        GetCaloUtils()->RecalculateClusterShowerShapeParameters(GetEMCALCells(),clus);
+        GetCaloUtils()->RecalculateClusterPID(clus);
+      }
+      
+      //Recalculate distance to bad channels, if new list of bad channels provided
+      GetCaloUtils()->RecalculateClusterDistanceToBadChannel(GetEMCALCells(),clus);
+      
+      //Recalculate cluster position
+      if(GetCaloUtils()->IsRecalculationOfClusterPositionOn()){
+        GetCaloUtils()->RecalculateClusterPosition(GetEMCALCells(),clus); 
+        //clus->GetPosition(pos);
+        //printf("After  Corrections: e %f, x %f, y %f, z %f\n",clus->E(),pos[0],pos[1],pos[2]);
+      }
     }
     
     //Correct non linearity
     if(GetCaloUtils()->IsCorrectionOfClusterEnergyOn()){
       GetCaloUtils()->CorrectClusterEnergy(clus) ;
       //printf("Linearity Corrected Energy %f\n",clus->E());  
-    }          
-    
-    //In case of MC analysis, to match resolution/calibration in real data
-    Float_t rdmEnergy = GetCaloUtils()->GetEMCALRecoUtils()->SmearClusterEnergy(clus);
-    // printf("\t Energy %f, smeared %f\n", clus->E(),rdmEnergy);
-    clus->SetE(rdmEnergy);
+      
+      //In case of MC analysis, to match resolution/calibration in real data
+      Float_t rdmEnergy = GetCaloUtils()->GetEMCALRecoUtils()->SmearClusterEnergy(clus);
+      // printf("\t Energy %f, smeared %f\n", clus->E(),rdmEnergy);
+      clus->SetE(rdmEnergy);
+    }
     
     if (fMixedEvent) 
       clus->SetID(iclus) ; 
@@ -764,6 +773,12 @@ void AliCaloTrackReader::FillInputEMCAL() {
   
   if(fDebug > 2 ) printf("AliCaloTrackReader::FillInputEMCAL()\n");
   
+  // First recalibrate cells, time or energy
+  //  if(GetCaloUtils()->IsRecalibrationOn())
+  //    GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCells(GetCaloUtils()->GetEMCALGeometry(), 
+  //                                                          GetEMCALCells(), 
+  //                                                          fInputEvent->GetBunchCrossNumber());
+  
   //Loop to select clusters in fiducial cut and fill container with aodClusters
   if(fEMCALClustersListName==""){
     Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
@@ -845,10 +860,14 @@ void AliCaloTrackReader::FillInputPHOS() {
             printf("AliCaloTrackReader::FillInputPHOS() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
                    momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
           
+          if(fRecalculateClusters){
+            
             //Recalibrate the cluster energy 
-          if(GetCaloUtils()->IsRecalibrationOn()) {
-            Float_t energy = GetCaloUtils()->RecalibrateClusterEnergy(clus, (AliAODCaloCells*)GetPHOSCells());
-            clus->SetE(energy);
+            if(GetCaloUtils()->IsRecalibrationOn()) {
+              Float_t energy = GetCaloUtils()->RecalibrateClusterEnergy(clus, (AliAODCaloCells*)GetPHOSCells());
+              clus->SetE(energy);
+            }
+            
           }
           
           if (fMixedEvent) {
index 2a2d843..e62e9cf 100755 (executable)
@@ -156,6 +156,10 @@ public:
   void             SwitchOnPHOSCells()                     { fFillPHOSCells = kTRUE        ; }
   void             SwitchOffPHOSCells()                    { fFillPHOSCells = kFALSE       ; }
 
+  Bool_t           AreClustersRecalculated()         const { return fRecalculateClusters   ; }
+  void             SwitchOnClusterRecalculation()          { fRecalculateClusters = kTRUE  ; }
+  void             SwitchOffClusterRecalculation()         { fRecalculateClusters = kFALSE ; }  
+  
   Bool_t           IsEmbeddedClusterSelectionOn()    const { return fSelectEmbeddedClusters   ; }
   void             SwitchOnEmbeddedClustersSelection()     { fSelectEmbeddedClusters = kTRUE  ; }
   void             SwitchOffEmbeddedClustersSelection()    { fSelectEmbeddedClusters = kFALSE ; }
@@ -353,7 +357,8 @@ public:
   Bool_t           fFillPHOS;       // use data from PHOS
   Bool_t           fFillEMCALCells; // use data from EMCAL
   Bool_t           fFillPHOSCells;  // use data from PHOS
-  Bool_t           fSelectEmbeddedClusters;   // Use only simulated clusters that come from embedding.
+  Bool_t           fRecalculateClusters;    // Correct clusters, recalculate them if recalibration parameters is given
+  Bool_t           fSelectEmbeddedClusters; // Use only simulated clusters that come from embedding.
   
   ULong_t          fTrackStatus        ; // Track selection bit, select tracks refitted in TPC, ITS ...
   ULong_t          fTrackFilterMask    ; // Track selection bit, for AODs (any difference with track status?)
@@ -396,7 +401,7 @@ public:
   Int_t            fCentralityBin[2];    // Minimum and maximum value of the centrality for the analysis
   TString          fEventPlaneMethod;    // Name of event plane method, by default "Q"
   
-  ClassDef(AliCaloTrackReader,32)
+  ClassDef(AliCaloTrackReader,33)
 } ;
 
 
index b653d96..9d1be34 100755 (executable)
@@ -144,12 +144,12 @@ fhGenMCE(),                fhGenMCEtaPhi(),
 fhGenMCAccE(),             fhGenMCAccEtaPhi(),   
 
 //matched MC
-fhEMVxyz(0),               fhEMR(0),                   fhHaVxyz(0),             fhHaR(0)
-//fh1pOverE(0),              fh1dR(0),                   fh2EledEdx(0),           fh2MatchdEdx(0),
-//fhMCEle1pOverE(0),         fhMCEle1dR(0),              fhMCEle2MatchdEdx(0),
-//fhMCChHad1pOverE(0),       fhMCChHad1dR(0),            fhMCChHad2MatchdEdx(0),
-//fhMCNeutral1pOverE(0),     fhMCNeutral1dR(0),          fhMCNeutral2MatchdEdx(0),fh1pOverER02(0),           
-//fhMCEle1pOverER02(0),      fhMCChHad1pOverER02(0),     fhMCNeutral1pOverER02(0)
+fhEMVxyz(0),               fhEMR(0),                   fhHaVxyz(0),             fhHaR(0),
+fh1pOverE(0),              fh1dR(0),                   fh2EledEdx(0),           fh2MatchdEdx(0),
+fhMCEle1pOverE(0),         fhMCEle1dR(0),              fhMCEle2MatchdEdx(0),
+fhMCChHad1pOverE(0),       fhMCChHad1dR(0),            fhMCChHad2MatchdEdx(0),
+fhMCNeutral1pOverE(0),     fhMCNeutral1dR(0),          fhMCNeutral2MatchdEdx(0),fh1pOverER02(0),           
+fhMCEle1pOverER02(0),      fhMCChHad1pOverER02(0),     fhMCNeutral1pOverER02(0)
 {
   //Default Ctor
   
@@ -199,9 +199,9 @@ TList *  AliAnaCalorimeterQA::GetCreateOutputObjects()
   Int_t netabins    = GetHistoEtaBins();          Float_t etamax    = GetHistoEtaMax();          Float_t etamin    = GetHistoEtaMin(); 
   Int_t nmassbins   = GetHistoMassBins();         Float_t massmax   = GetHistoMassMax();              Float_t massmin   = GetHistoMassMin();
   Int_t nasymbins   = GetHistoAsymmetryBins();    Float_t asymmax   = GetHistoAsymmetryMax();    Float_t asymmin   = GetHistoAsymmetryMin();
-  //Int_t nPoverEbins = GetHistoPOverEBins();       Float_t pOverEmax = GetHistoPOverEMax();       Float_t pOverEmin = GetHistoPOverEMin();
-  //Int_t ndedxbins   = GetHistodEdxBins();         Float_t dedxmax   = GetHistodEdxMax();         Float_t dedxmin   = GetHistodEdxMin();
-  //Int_t ndRbins     = GetHistodRBins();           Float_t dRmax     = GetHistodRMax();           Float_t dRmin     = GetHistodRMin();
+  Int_t nPoverEbins = GetHistoPOverEBins();       Float_t pOverEmax = GetHistoPOverEMax();       Float_t pOverEmin = GetHistoPOverEMin();
+  Int_t ndedxbins   = GetHistodEdxBins();         Float_t dedxmax   = GetHistodEdxMax();         Float_t dedxmin   = GetHistodEdxMin();
+  Int_t ndRbins     = GetHistodRBins();           Float_t dRmax     = GetHistodRMax();           Float_t dRmin     = GetHistodRMin();
   Int_t ntimebins   = GetHistoTimeBins();         Float_t timemax   = GetHistoTimeMax();         Float_t timemin   = GetHistoTimeMin();       
   Int_t nclbins     = GetHistoNClustersBins();    Int_t   nclmax    = GetHistoNClustersMax();    Int_t   nclmin    = GetHistoNClustersMin(); 
   Int_t ncebins     = GetHistoNCellsBins();       Int_t   ncemax    = GetHistoNCellsMax();       Int_t   ncemin    = GetHistoNCellsMin(); 
@@ -548,29 +548,29 @@ TList *  AliAnaCalorimeterQA::GetCreateOutputObjects()
       outputContainer->Add(fhEtaPhiECharged);  
     }
     
-//    fh1pOverE = new TH2F("h1pOverE","TRACK matches p/E",nptbins,ptmin,ptmax, nPoverEbins,pOverEmin,pOverEmax);
-//    fh1pOverE->SetYTitle("p/E");
-//    fh1pOverE->SetXTitle("p_{T} (GeV/c)");
-//    outputContainer->Add(fh1pOverE);
-//    
-//    fh1dR = new TH1F("h1dR","TRACK matches dR",ndRbins,dRmin,dRmax);
-//    fh1dR->SetXTitle("#Delta R (rad)");
-//    outputContainer->Add(fh1dR) ;
-//    
-//    fh2MatchdEdx = new TH2F("h2MatchdEdx","dE/dx vs. p for all matches",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
-//    fh2MatchdEdx->SetXTitle("p (GeV/c)");
-//    fh2MatchdEdx->SetYTitle("<dE/dx>");
-//    outputContainer->Add(fh2MatchdEdx);
-//    
-//    fh2EledEdx = new TH2F("h2EledEdx","dE/dx vs. p for electrons",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
-//    fh2EledEdx->SetXTitle("p (GeV/c)");
-//    fh2EledEdx->SetYTitle("<dE/dx>");
-//    outputContainer->Add(fh2EledEdx) ;
-//    
-//    fh1pOverER02 = new TH2F("h1pOverER02","TRACK matches p/E, all",nptbins,ptmin,ptmax, nPoverEbins,pOverEmin,pOverEmax);
-//    fh1pOverER02->SetYTitle("p/E");
-//    fh1pOverER02->SetXTitle("p_{T} (GeV/c)");
-//    outputContainer->Add(fh1pOverER02);      
+    fh1pOverE = new TH2F("h1pOverE","TRACK matches p/E",nptbins,ptmin,ptmax, nPoverEbins,pOverEmin,pOverEmax);
+    fh1pOverE->SetYTitle("p/E");
+    fh1pOverE->SetXTitle("p_{T} (GeV/c)");
+    outputContainer->Add(fh1pOverE);
+    
+    fh1dR = new TH1F("h1dR","TRACK matches dR",ndRbins,dRmin,dRmax);
+    fh1dR->SetXTitle("#Delta R (rad)");
+    outputContainer->Add(fh1dR) ;
+    
+    fh2MatchdEdx = new TH2F("h2MatchdEdx","dE/dx vs. p for all matches",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
+    fh2MatchdEdx->SetXTitle("p (GeV/c)");
+    fh2MatchdEdx->SetYTitle("<dE/dx>");
+    outputContainer->Add(fh2MatchdEdx);
+    
+    fh2EledEdx = new TH2F("h2EledEdx","dE/dx vs. p for electrons",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
+    fh2EledEdx->SetXTitle("p (GeV/c)");
+    fh2EledEdx->SetYTitle("<dE/dx>");
+    outputContainer->Add(fh2EledEdx) ;
+    
+    fh1pOverER02 = new TH2F("h1pOverER02","TRACK matches p/E, all",nptbins,ptmin,ptmax, nPoverEbins,pOverEmin,pOverEmax);
+    fh1pOverER02->SetYTitle("p/E");
+    fh1pOverER02->SetXTitle("p_{T} (GeV/c)");
+    outputContainer->Add(fh1pOverER02);        
   }
   
   if(fFillAllPi0Histo){
@@ -1113,62 +1113,62 @@ TList *  AliAnaCalorimeterQA::GetCreateOutputObjects()
     
     //Track Matching 
     
-//    fhMCEle1pOverE = new TH2F("hMCEle1pOverE","TRACK matches p/E, MC electrons",nptbins,ptmin,ptmax, nPoverEbins,pOverEmin,pOverEmax);
-//    fhMCEle1pOverE->SetYTitle("p/E");
-//    fhMCEle1pOverE->SetXTitle("p_{T} (GeV/c)");
-//    outputContainer->Add(fhMCEle1pOverE);
-//    
-//    fhMCEle1dR = new TH1F("hMCEle1dR","TRACK matches dR, MC electrons",ndRbins,dRmin,dRmax);
-//    fhMCEle1dR->SetXTitle("#Delta R (rad)");
-//    outputContainer->Add(fhMCEle1dR) ;
-//    
-//    fhMCEle2MatchdEdx = new TH2F("hMCEle2MatchdEdx","dE/dx vs. p for all matches, MC electrons",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
-//    fhMCEle2MatchdEdx->SetXTitle("p (GeV/c)");
-//    fhMCEle2MatchdEdx->SetYTitle("<dE/dx>");
-//    outputContainer->Add(fhMCEle2MatchdEdx);
-//    
-//    fhMCChHad1pOverE = new TH2F("hMCChHad1pOverE","TRACK matches p/E, MC charged hadrons",nptbins,ptmin,ptmax, nPoverEbins,pOverEmin,pOverEmax);
-//    fhMCChHad1pOverE->SetYTitle("p/E");
-//    fhMCChHad1pOverE->SetXTitle("p_{T} (GeV/c)");
-//    outputContainer->Add(fhMCChHad1pOverE);
-//    
-//    fhMCChHad1dR = new TH1F("hMCChHad1dR","TRACK matches dR, MC charged hadrons",ndRbins,dRmin,dRmax);
-//    fhMCChHad1dR->SetXTitle("#Delta R (rad)");
-//    outputContainer->Add(fhMCChHad1dR) ;
-//    
-//    fhMCChHad2MatchdEdx = new TH2F("hMCChHad2MatchdEdx","dE/dx vs. p for all matches, MC charged hadrons",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
-//    fhMCChHad2MatchdEdx->SetXTitle("p (GeV/c)");
-//    fhMCChHad2MatchdEdx->SetYTitle("<dE/dx>");
-//    outputContainer->Add(fhMCChHad2MatchdEdx);
-//    
-//    fhMCNeutral1pOverE = new TH2F("hMCNeutral1pOverE","TRACK matches p/E, MC neutrals",nptbins,ptmin,ptmax, nPoverEbins,pOverEmin,pOverEmax);
-//    fhMCNeutral1pOverE->SetYTitle("p/E");
-//    fhMCNeutral1pOverE->SetXTitle("p_{T} (GeV/c)");
-//    outputContainer->Add(fhMCNeutral1pOverE);
-//    
-//    fhMCNeutral1dR = new TH1F("hMCNeutral1dR","TRACK matches dR, MC neutrals",ndRbins,dRmin,dRmax);
-//    fhMCNeutral1dR->SetXTitle("#Delta R (rad)");
-//    outputContainer->Add(fhMCNeutral1dR) ;
-//    
-//    fhMCNeutral2MatchdEdx = new TH2F("hMCNeutral2MatchdEdx","dE/dx vs. p for all matches, MC neutrals",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
-//    fhMCNeutral2MatchdEdx->SetXTitle("p (GeV/c)");
-//    fhMCNeutral2MatchdEdx->SetYTitle("<dE/dx>");
-//    outputContainer->Add(fhMCNeutral2MatchdEdx);
-//    
-//    fhMCEle1pOverER02 = new TH2F("hMCEle1pOverER02","TRACK matches p/E, MC electrons",nptbins,ptmin,ptmax, nPoverEbins,pOverEmin,pOverEmax);
-//    fhMCEle1pOverER02->SetYTitle("p/E");
-//    fhMCEle1pOverER02->SetXTitle("p_{T} (GeV/c)");
-//    outputContainer->Add(fhMCEle1pOverER02);
-//    
-//    fhMCChHad1pOverER02 = new TH2F("hMCChHad1pOverER02","TRACK matches p/E, MC charged hadrons",nptbins,ptmin,ptmax, nPoverEbins,pOverEmin,pOverEmax);
-//    fhMCChHad1pOverER02->SetYTitle("p/E");
-//    fhMCChHad1pOverER02->SetXTitle("p_{T} (GeV/c)");
-//    outputContainer->Add(fhMCChHad1pOverER02);
-//    
-//    fhMCNeutral1pOverER02 = new TH2F("hMCNeutral1pOverER02","TRACK matches p/E, MC neutrals",nptbins,ptmin,ptmax, nPoverEbins,pOverEmin,pOverEmax);
-//    fhMCNeutral1pOverER02->SetYTitle("p/E");
-//    fhMCNeutral1pOverER02->SetXTitle("p_{T} (GeV/c)");
-//    outputContainer->Add(fhMCNeutral1pOverER02);
+    fhMCEle1pOverE = new TH2F("hMCEle1pOverE","TRACK matches p/E, MC electrons",nptbins,ptmin,ptmax, nPoverEbins,pOverEmin,pOverEmax);
+    fhMCEle1pOverE->SetYTitle("p/E");
+    fhMCEle1pOverE->SetXTitle("p_{T} (GeV/c)");
+    outputContainer->Add(fhMCEle1pOverE);
+    
+    fhMCEle1dR = new TH1F("hMCEle1dR","TRACK matches dR, MC electrons",ndRbins,dRmin,dRmax);
+    fhMCEle1dR->SetXTitle("#Delta R (rad)");
+    outputContainer->Add(fhMCEle1dR) ;
+    
+    fhMCEle2MatchdEdx = new TH2F("hMCEle2MatchdEdx","dE/dx vs. p for all matches, MC electrons",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
+    fhMCEle2MatchdEdx->SetXTitle("p (GeV/c)");
+    fhMCEle2MatchdEdx->SetYTitle("<dE/dx>");
+    outputContainer->Add(fhMCEle2MatchdEdx);
+    
+    fhMCChHad1pOverE = new TH2F("hMCChHad1pOverE","TRACK matches p/E, MC charged hadrons",nptbins,ptmin,ptmax, nPoverEbins,pOverEmin,pOverEmax);
+    fhMCChHad1pOverE->SetYTitle("p/E");
+    fhMCChHad1pOverE->SetXTitle("p_{T} (GeV/c)");
+    outputContainer->Add(fhMCChHad1pOverE);
+    
+    fhMCChHad1dR = new TH1F("hMCChHad1dR","TRACK matches dR, MC charged hadrons",ndRbins,dRmin,dRmax);
+    fhMCChHad1dR->SetXTitle("#Delta R (rad)");
+    outputContainer->Add(fhMCChHad1dR) ;
+    
+    fhMCChHad2MatchdEdx = new TH2F("hMCChHad2MatchdEdx","dE/dx vs. p for all matches, MC charged hadrons",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
+    fhMCChHad2MatchdEdx->SetXTitle("p (GeV/c)");
+    fhMCChHad2MatchdEdx->SetYTitle("<dE/dx>");
+    outputContainer->Add(fhMCChHad2MatchdEdx);
+    
+    fhMCNeutral1pOverE = new TH2F("hMCNeutral1pOverE","TRACK matches p/E, MC neutrals",nptbins,ptmin,ptmax, nPoverEbins,pOverEmin,pOverEmax);
+    fhMCNeutral1pOverE->SetYTitle("p/E");
+    fhMCNeutral1pOverE->SetXTitle("p_{T} (GeV/c)");
+    outputContainer->Add(fhMCNeutral1pOverE);
+    
+    fhMCNeutral1dR = new TH1F("hMCNeutral1dR","TRACK matches dR, MC neutrals",ndRbins,dRmin,dRmax);
+    fhMCNeutral1dR->SetXTitle("#Delta R (rad)");
+    outputContainer->Add(fhMCNeutral1dR) ;
+    
+    fhMCNeutral2MatchdEdx = new TH2F("hMCNeutral2MatchdEdx","dE/dx vs. p for all matches, MC neutrals",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
+    fhMCNeutral2MatchdEdx->SetXTitle("p (GeV/c)");
+    fhMCNeutral2MatchdEdx->SetYTitle("<dE/dx>");
+    outputContainer->Add(fhMCNeutral2MatchdEdx);
+    
+    fhMCEle1pOverER02 = new TH2F("hMCEle1pOverER02","TRACK matches p/E, MC electrons",nptbins,ptmin,ptmax, nPoverEbins,pOverEmin,pOverEmax);
+    fhMCEle1pOverER02->SetYTitle("p/E");
+    fhMCEle1pOverER02->SetXTitle("p_{T} (GeV/c)");
+    outputContainer->Add(fhMCEle1pOverER02);
+    
+    fhMCChHad1pOverER02 = new TH2F("hMCChHad1pOverER02","TRACK matches p/E, MC charged hadrons",nptbins,ptmin,ptmax, nPoverEbins,pOverEmin,pOverEmax);
+    fhMCChHad1pOverER02->SetYTitle("p/E");
+    fhMCChHad1pOverER02->SetXTitle("p_{T} (GeV/c)");
+    outputContainer->Add(fhMCChHad1pOverER02);
+    
+    fhMCNeutral1pOverER02 = new TH2F("hMCNeutral1pOverER02","TRACK matches p/E, MC neutrals",nptbins,ptmin,ptmax, nPoverEbins,pOverEmin,pOverEmax);
+    fhMCNeutral1pOverER02->SetYTitle("p/E");
+    fhMCNeutral1pOverER02->SetXTitle("p_{T} (GeV/c)");
+    outputContainer->Add(fhMCNeutral1pOverER02);
   }
   
   //  for(Int_t i = 0; i < outputContainer->GetEntries() ; i++)
@@ -1228,22 +1228,25 @@ void AliAnaCalorimeterQA::Print(const Option_t * opt) const
 void  AliAnaCalorimeterQA::MakeAnalysisFillHistograms() 
 {
   //Fill Calorimeter QA histograms
+  
+  
   TLorentzVector mom  ;
   TLorentzVector mom2 ;
   TObjArray * caloClusters = NULL;
-  Int_t  nLabel = 0  ;
-  Int_t *labels = 0x0;
-  Int_t  nCaloClusters         = 0;
-  Int_t  nCaloClustersAccepted = 0;
-  Int_t  nCaloCellsPerCluster  = 0;
-  Bool_t matched    = kFALSE;
-  Int_t  nModule    = -1;
-  
+  Int_t  nLabel                = 0  ;
+  Int_t *labels                = 0x0;
+  Int_t  nCaloClusters         = 0  ;
+  Int_t  nCaloClustersAccepted = 0  ;
+  Int_t  nCaloCellsPerCluster  = 0  ;
+  Bool_t matched               = kFALSE;
+  Int_t  nModule               =-1  ;
+  Int_t  bc                    = GetReader()->GetInputEvent()->GetBunchCrossNumber();
+
   //Get vertex for photon momentum calculation and event selection
   Double_t v[3] = {0,0,0}; //vertex ;
   GetReader()->GetVertex(v);
   if (TMath::Abs(v[2]) > GetZvertexCut()) return ;  
-  
+
   //Play with the MC stack if available        
   //Get the MC arrays and do some checks
   if(IsDataMC()){
@@ -1255,7 +1258,7 @@ void  AliAnaCalorimeterQA::MakeAnalysisFillHistograms()
       //Fill some pure MC histograms, only primaries.
       for(Int_t i=0 ; i<GetMCStack()->GetNprimary(); i++){//Only primary particles, for all MC transport put GetNtrack()
         TParticle *primary = GetMCStack()->Particle(i) ;
-        //printf("i %d, %s: status = %d, primary? %d\n",i, primary->GetName(), primary->GetStatusCode(), primary->IsPrimary());
+        
         if (primary->GetStatusCode() > 11) continue; //Working for PYTHIA and simple generators, check for HERWIG 
         primary->Momentum(mom);
         MCHistograms(mom,TMath::Abs(primary->GetPdgCode()));
@@ -1269,11 +1272,9 @@ void  AliAnaCalorimeterQA::MakeAnalysisFillHistograms()
       //Fill some pure MC histograms, only primaries.
       for(Int_t i=0 ; i < (GetReader()->GetAODMCParticles(0))->GetEntriesFast(); i++){
         AliAODMCParticle *aodprimary = (AliAODMCParticle*) (GetReader()->GetAODMCParticles(0))->At(i) ;
-        //printf("i %d, %s: primary? %d physical primary? %d, flag %d\n",
-        //        i,(TDatabasePDG::Instance()->GetParticle(aodprimary->GetPdgCode()))->GetName(), 
-        //        aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary(), aodprimary->GetFlag());
+        
         if (!aodprimary->IsPrimary()) continue; //accept all which is not MC transport generated. Don't know how to avoid partons
-        //aodprimary->Momentum(mom);
+        
         mom.SetPxPyPzE(aodprimary->Px(),aodprimary->Py(),aodprimary->Pz(),aodprimary->E());
         MCHistograms(mom,TMath::Abs(aodprimary->GetPdgCode()));
       } //primary loop
@@ -1281,436 +1282,503 @@ void  AliAnaCalorimeterQA::MakeAnalysisFillHistograms()
     }
   }// is data and MC   
   
-  
   //Get List with CaloClusters  
   if      (fCalorimeter == "PHOS")  caloClusters = GetPHOSClusters();
   else if (fCalorimeter == "EMCAL") caloClusters = GetEMCALClusters();
   else 
     AliFatal(Form("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - Wrong calorimeter name <%s>, END\n", fCalorimeter.Data()));
   
-  
-  if(!caloClusters) {
-    AliFatal(Form("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - No CaloClusters available\n"));
+  //Get List with CaloCells
+  AliVCaloCells * cell = 0x0; 
+  if(fCalorimeter == "PHOS") cell =  GetPHOSCells();
+  else                                  cell =  GetEMCALCells();
+    
+  if(!caloClusters || !cell) {
+    AliFatal(Form("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - No CaloClusters or CaloCells available\n"));
+    return; // trick coverity
   }
-  else{
+  
+  //----------------------------------------------------------
+  //Correlate Calorimeters and V0 and track Multiplicity
+  //----------------------------------------------------------
+  
+  if(fCorrelate)       Correlate();
+  
+  //----------------------------------------------------------
+  // CALOCLUSTERS
+  //----------------------------------------------------------
+  
+  nCaloClusters = caloClusters->GetEntriesFast() ; 
+  Int_t *nClustersInModule = new Int_t[fNModules];
+  for(Int_t imod = 0; imod < fNModules; imod++ ) nClustersInModule[imod] = 0;
+  
+  if(GetDebug() > 0)
+    printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - In %s there are %d clusters \n", fCalorimeter.Data(), nCaloClusters);
+  
+  AliVTrack * track = 0x0;
+  Float_t pos[3] ;
+  Double_t tof = 0;
+  
+  //Loop over CaloClusters
+  //if(nCaloClusters > 0)printf("QA  : Vertex Cut passed %f, cut %f, entries %d, %s\n",v[2], 40., nCaloClusters, fCalorimeter.Data());
+  for(Int_t iclus = 0; iclus < nCaloClusters; iclus++){
     
-    //----------------------------------------------------------
-    //Correlate Calorimeters and V0 and track Multiplicity
-    //----------------------------------------------------------
-
-    if(fCorrelate)     Correlate();
+    if(GetDebug() > 0) printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - cluster: %d/%d, data %d \n",
+                              iclus+1,nCaloClusters,GetReader()->GetDataType());
     
-    //----------------------------------------------------------
-    // CALOCLUSTERS
-    //----------------------------------------------------------
+    AliVCluster* clus =  (AliVCluster*)caloClusters->At(iclus);
+        
+    //Get cluster kinematics
+    clus->GetPosition(pos);
+    clus->GetMomentum(mom,v);
     
-    nCaloClusters = caloClusters->GetEntriesFast() ; 
-    Int_t *nClustersInModule = new Int_t[fNModules];
-    for(Int_t imod = 0; imod < fNModules; imod++ ) nClustersInModule[imod] = 0;
+    //Check only certain regions
+    Bool_t in = kTRUE;
+    if(IsFiducialCutOn()) in =  GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
+    if(!in) continue;
     
+    //MC labels
+    nLabel = clus->GetNLabels();
+    labels = clus->GetLabels();
     
-    if(GetDebug() > 0)
-      printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - In %s there are %d clusters \n", fCalorimeter.Data(), nCaloClusters);
+    //Cells per cluster
+    nCaloCellsPerCluster = clus->GetNCells();
+    //if(mom.E() > 10 && nCaloCellsPerCluster == 1 ) printf("%s:************** E = %f ********** ncells = %d\n",fCalorimeter.Data(), mom.E(),nCaloCellsPerCluster);
     
-    //AliVTrack * track = 0x0;
-    Float_t pos[3] ;
-    Double_t tof = 0;
-    //Loop over CaloClusters
-    //if(nCaloClusters > 0)printf("QA  : Vertex Cut passed %f, cut %f, entries %d, %s\n",v[2], 40., nCaloClusters, fCalorimeter.Data());
-    for(Int_t iclus = 0; iclus < nCaloClusters; iclus++){
-      
-      if(GetDebug() > 0) printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - cluster: %d/%d, data %d \n",
-                                iclus+1,nCaloClusters,GetReader()->GetDataType());
-      
-      AliVCluster* clus =  (AliVCluster*)caloClusters->At(iclus);
-      AliVCaloCells * cell = 0x0; 
-      if(fCalorimeter == "PHOS") cell =  GetPHOSCells();
-      else                                      cell =  GetEMCALCells();
-      
-      //Get cluster kinematics
-      clus->GetPosition(pos);
-      clus->GetMomentum(mom,v);
-      tof = clus->GetTOF()*1e9;
-      if(tof < fTimeCutMin || tof > fTimeCutMax) continue;
-      
-      //Check only certain regions
-      Bool_t in = kTRUE;
-      if(IsFiducialCutOn()) in =  GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
-      if(!in) continue;
-      
-      //MC labels
-      nLabel = clus->GetNLabels();
-      labels = clus->GetLabels();
-      
-      //Cells per cluster
-      nCaloCellsPerCluster = clus->GetNCells();
-      //if(mom.E() > 10 && nCaloCellsPerCluster == 1 ) printf("%s:************** E = %f ********** ncells = %d\n",fCalorimeter.Data(), mom.E(),nCaloCellsPerCluster);
-      // Cluster mathed with track?
-      matched = GetCaloPID()->IsTrackMatched(clus,GetCaloUtils());
-
-      
-      //======================
-      //Cells in cluster
-      //======================
-      
-      //Get list of contributors
-      UShort_t * indexList = clus->GetCellsAbsId() ;
-      
-      if(fFillAllPosHisto) FillCellPositionHistograms(nCaloCellsPerCluster,indexList,pos,mom.E());
-      
-      // Get the fraction of the cluster energy that carries the cell with highest energy
-      Float_t maxCellFraction = 0.;
-      Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(cell, clus,maxCellFraction);
-      Int_t smMax =0; Int_t ietaaMax=-1; Int_t iphiiMax = 0; Int_t rcuMax = 0;
-      smMax = GetModuleNumberCellIndexes(absIdMax,fCalorimeter, ietaaMax, iphiiMax, rcuMax);
-      Int_t dIeta = 0;
-      Int_t dIphi = 0;
-      Double_t tmax  = cell->GetCellTime(absIdMax)*1e9;
-      //Float_t  emax  = cell->GetCellAmplitude(absIdMax);
-      
-      if     (clus->E() < 2.){
-        fhLambda0vsClusterMaxCellDiffE0->Fill(clus->GetM02(),      maxCellFraction);
-        fhNCellsvsClusterMaxCellDiffE0 ->Fill(nCaloCellsPerCluster,maxCellFraction);
-      }
-      else if(clus->E() < 6.){
-        fhLambda0vsClusterMaxCellDiffE2->Fill(clus->GetM02(),      maxCellFraction);
-        fhNCellsvsClusterMaxCellDiffE2 ->Fill(nCaloCellsPerCluster,maxCellFraction);
-      }
-      else{
-        fhLambda0vsClusterMaxCellDiffE6->Fill(clus->GetM02(),      maxCellFraction);  
-        fhNCellsvsClusterMaxCellDiffE6 ->Fill(nCaloCellsPerCluster,maxCellFraction);
+    // Cluster mathed with track?
+    matched = GetCaloPID()->IsTrackMatched(clus,GetCaloUtils());
+    
+    
+    //======================
+    //Cells in cluster
+    //======================
+    
+    //Get list of contributors
+    UShort_t * indexList = clus->GetCellsAbsId() ;
+    
+    if(fFillAllPosHisto) FillCellPositionHistograms(nCaloCellsPerCluster,indexList,pos,mom.E());
+    
+    // Get the fraction of the cluster energy that carries the cell with highest energy
+    Float_t maxCellFraction = 0.;
+    Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(cell, clus,maxCellFraction);
+    Int_t smMax =0; Int_t ietaaMax=-1; Int_t iphiiMax = 0; Int_t rcuMax = 0;
+    smMax = GetModuleNumberCellIndexes(absIdMax,fCalorimeter, ietaaMax, iphiiMax, rcuMax);
+    Int_t dIeta = 0;
+    Int_t dIphi = 0;
+    
+    if     (clus->E() < 2.){
+      fhLambda0vsClusterMaxCellDiffE0->Fill(clus->GetM02(),      maxCellFraction);
+      fhNCellsvsClusterMaxCellDiffE0 ->Fill(nCaloCellsPerCluster,maxCellFraction);
+    }
+    else if(clus->E() < 6.){
+      fhLambda0vsClusterMaxCellDiffE2->Fill(clus->GetM02(),      maxCellFraction);
+      fhNCellsvsClusterMaxCellDiffE2 ->Fill(nCaloCellsPerCluster,maxCellFraction);
+    }
+    else{
+      fhLambda0vsClusterMaxCellDiffE6->Fill(clus->GetM02(),      maxCellFraction);  
+      fhNCellsvsClusterMaxCellDiffE6 ->Fill(nCaloCellsPerCluster,maxCellFraction);
+    }
+    
+    fhNCellsPerClusterNoCut  ->Fill(clus->E(), nCaloCellsPerCluster);
+    nModule = GetModuleNumber(clus);
+    if(nModule >=0 && nModule < fNModules) fhNCellsPerClusterModNoCut[nModule]->Fill(clus->E(), nCaloCellsPerCluster);
+    
+    fhClusterMaxCellDiffNoCut->Fill(clus->E(),maxCellFraction);
+    
+    // Look inside the cluster    
+    // Calculate average time of cells in cluster and weighted average
+    Float_t  rawEnergy         = 0;
+    Float_t  factor            = 1;
+    Double_t time              = 0;
+    Double_t tmax              = cell->GetCellTime(absIdMax);
+    Double_t averTime          = 0;
+    Double_t weightedTime      = 0;
+    Double_t weight            = 0;
+    Double_t averTimeNoMax     = 0;
+    Double_t weightedTimeNoMax = 0;
+    Double_t weightNoMax       = 0;
+    
+    for (Int_t ipos = 0; ipos < nCaloCellsPerCluster; ipos++) {
+      Int_t icol     = -1;
+      Int_t irow     = -1;
+      Int_t iRCU     = -1;
+      Int_t id       = indexList[ipos];
+      Int_t nModuleCell = GetModuleNumberCellIndexes(id,fCalorimeter, icol, irow, iRCU);
+      
+      //Recalibrate energy
+      if(GetCaloUtils()->IsRecalibrationOn()){
+        if(fCalorimeter == "PHOS") {
+          factor = GetCaloUtils()->GetPHOSChannelRecalibrationFactor(nModuleCell,icol,irow);
+        }
+        else                              {
+          factor = GetCaloUtils()->GetEMCALChannelRecalibrationFactor(nModuleCell,icol,irow);
+          
+        }
       }
       
-      fhNCellsPerClusterNoCut  ->Fill(clus->E(), nCaloCellsPerCluster);
-      nModule = GetModuleNumber(clus);
-      if(nModule >=0 && nModule < fNModules) fhNCellsPerClusterModNoCut[nModule]->Fill(clus->E(), nCaloCellsPerCluster);
+      //Recalibrate time
+      time      = cell->GetCellTime(id);
+      if(fCalorimeter=="EMCAL" && GetCaloUtils()->GetEMCALRecoUtils()->IsTimeRecalibrationOn()){
+        GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTime(id,bc,time);
+        if(id==absIdMax) tmax = time;
+      }
       
-      fhClusterMaxCellDiffNoCut->Fill(clus->E(),maxCellFraction);
+//      printf("QA  cluster       : Id %d, Time org %e, Time new %e; Amp org %f, Amp new %f\n",
+//             id, cell->GetCellTime(id),time, cell->GetCellAmplitude(id),cell->GetCellAmplitude(id)*factor);
       
-      // Calculate average time of cells in cluster and weighted average
-      Double_t averTime     = 0;
-      Double_t weightedTime = 0;
-      Double_t weight       = 0;
-      Double_t averTimeNoMax     = 0;
-      Double_t weightedTimeNoMax = 0;
-      Double_t weightNoMax       = 0;
+      // Get the energy of the cluster without the non linearity correction
+      rawEnergy    += cell->GetCellAmplitude(id)*factor;
       
-      //if(GetReader()->GetDataType()==AliCaloTrackReader::kESD) {
-        Float_t rawEnergy = 0;
-        for (Int_t ipos = 0; ipos < nCaloCellsPerCluster; ipos++) 
-          rawEnergy += cell->GetCellAmplitude(indexList[ipos]);
+      Double_t w    = GetCaloUtils()->GetEMCALRecoUtils()->GetCellWeight(cell->GetCellAmplitude(id),rawEnergy);
+      averTime     += time*1e9;
+      weightedTime += time*1e9 * w;
+      weight       += w;
+      if(id != absIdMax){
+        averTimeNoMax     += time*1e9;
+        weightedTimeNoMax += time*1e9 * w;
+        weightNoMax       += w;
+      }
+    }        
+    
+    tmax*=1.e9;
+    averTime      /= nCaloCellsPerCluster;
+    if(nCaloCellsPerCluster > 1 ) averTimeNoMax /= (nCaloCellsPerCluster-1);
+    
+    if(weight > 0 )     tof = weightedTime / weight;        
+    else { 
+
+      tof = clus->GetTOF()*1e9;
+
+      printf("AliAnaCalorimeterQA:: Null weight! Investigate: E %f GeV, ncells %d, time max cell %f ns, average time %f ns, absIdMax %d, SM %d\n",
+             rawEnergy,nCaloCellsPerCluster, tmax, averTime,absIdMax,GetModuleNumber(clus));
+    }
+    
+    //printf("QA: E %f, tof org %g, tof new %g, ncells %d\n",clus->E(),clus->GetTOF()*1.e9,tof, clus->GetNCells());
+    
+    if(weightNoMax > 0 ) weightedTimeNoMax /= weightNoMax;
+    
+    //Cut on time of clusters
+    if(tof < fTimeCutMin || tof > fTimeCutMax){ 
+      printf("Remove cluster with TOF %f\n",tof);
+      continue;
+    }
+    
+    //======================
+    //Bad clusters selection
+    //======================
+    //Check bad clusters if rejection was not on
+    
+    Bool_t badCluster = kFALSE;
+    if(fCalorimeter=="EMCAL" && !GetCaloUtils()->GetEMCALRecoUtils()->IsRejectExoticCluster()){
+      //Bad clusters histograms
+      Float_t minNCells = 1;
+      if(clus->E() > 7) minNCells = TMath::Max(1,TMath::FloorNint(1 + TMath::Log(clus->E() - 7 )*1.7 ));
+      if(nCaloCellsPerCluster <= minNCells) {
+        //if(clus->GetM02() < 0.05) {
         
-        for (Int_t ipos = 0; ipos < nCaloCellsPerCluster; ipos++) {
-          Double_t w    = TMath::Max( 0., 4.5 + TMath::Log(cell->GetCellAmplitude(indexList[ipos]) / rawEnergy ));
-          averTime     += cell->GetCellTime(indexList[ipos])*1e9;
-          weightedTime += cell->GetCellTime(indexList[ipos])*1e9 * w;
-          weight       += w;
-          if(indexList[ipos]!=absIdMax){
-            averTimeNoMax     += cell->GetCellTime(indexList[ipos])*1e9;
-            weightedTimeNoMax += cell->GetCellTime(indexList[ipos])*1e9 * w;
-            weightNoMax       += w;
-          }
-        }        
+        badCluster = kTRUE;
         
-        averTime      /= nCaloCellsPerCluster;
-        if(nCaloCellsPerCluster > 1 ) averTimeNoMax /= (nCaloCellsPerCluster-1);
-
-       if(weight > 0 )      weightedTime /= weight;        
-       else {         
-         printf("AliAnaCalorimeterQA:: Null weight! Investigate: E %f GeV, ncells %d, time max cell %f ns, average time %f ns, absIdMax %d, SM %d\n",
-                   rawEnergy,nCaloCellsPerCluster, tmax, averTime,absIdMax,GetModuleNumber(clus));
-       }
+        fhBadClusterEnergy     ->Fill(clus->E());
+        fhBadClusterTimeEnergy ->Fill(clus->E(),tof);
+        fhBadClusterMaxCellDiff->Fill(clus->E(),maxCellFraction);
         
-       if(weightNoMax > 0 ) weightedTimeNoMax /= weightNoMax;
-
-      //} // only possible in ESDs
-      
-      //======================
-      //Bad clusters selection
-      //======================
-      //Check bad clusters if rejection was not on
-      
-      Bool_t badCluster = kFALSE;
-      if(fCalorimeter=="EMCAL" && !GetCaloUtils()->GetEMCALRecoUtils()->IsRejectExoticCluster()){
-        //Bad clusters histograms
-        Float_t minNCells = 1;
-        if(clus->E() > 7) minNCells = TMath::Max(1,TMath::FloorNint(1 + TMath::Log(clus->E() - 7 )*1.7 ));
-        if(nCaloCellsPerCluster <= minNCells) {
-          //if(clus->GetM02() < 0.05) {
+        if(clus->GetM02() > 0 || TMath::Abs(clus->GetM20()) > 0 || clus->GetDispersion() > 0){
           
-          badCluster = kTRUE;
+          fhBadClusterL0         ->Fill(clus->E(),clus->GetM02());
+          fhBadClusterL1         ->Fill(clus->E(),clus->GetM20());
+          fhBadClusterD          ->Fill(clus->E(),clus->GetDispersion());
           
-          fhBadClusterEnergy     ->Fill(clus->E());
-          fhBadClusterTimeEnergy ->Fill(clus->E(),tof);
-          fhBadClusterMaxCellDiff->Fill(clus->E(),maxCellFraction);
-
-          if(clus->GetM02() > 0 || TMath::Abs(clus->GetM20()) > 0 || clus->GetDispersion() > 0){
-            
-            fhBadClusterL0         ->Fill(clus->E(),clus->GetM02());
-            fhBadClusterL1         ->Fill(clus->E(),clus->GetM20());
-            fhBadClusterD          ->Fill(clus->E(),clus->GetDispersion());
-            
-          }
+        }
+        
+        //Clusters in event time difference
+        
+        for(Int_t iclus2 = 0; iclus2 < nCaloClusters; iclus2++ ){
           
-          //Clusters in event time difference
+          AliVCluster* clus2 =  (AliVCluster*)caloClusters->At(iclus2);
           
-          for(Int_t iclus2 = 0; iclus2 < nCaloClusters; iclus2++ ){
-            
-            AliVCluster* clus2 =  (AliVCluster*)caloClusters->At(iclus2);
-            
-            if(clus->GetID()==clus2->GetID()) continue;
-            
-            if(clus->GetM02() > 0.01) 
-              fhBadClusterPairDiffTimeE  ->Fill(clus->E(), tof-clus2->GetTOF()*1.e9);
-            
-          } // loop
+          if(clus->GetID()==clus2->GetID()) continue;
           
-          // Max cell compared to other cells in cluster
-          if(GetReader()->GetDataType()==AliCaloTrackReader::kESD) {
-            fhBadClusterMaxCellDiffAverageTime->Fill(clus->E(),tmax-averTime);
-            fhBadClusterMaxCellDiffWeightTime ->Fill(clus->E(),tmax-weightedTime);
-            fhBadClusterDiffWeightAverTime    ->Fill(clus->E(),weightedTime-averTime);
-            fhBadClusterMaxCellDiffAverageNoMaxTime->Fill(clus->E(),tmax-averTimeNoMax);
-            fhBadClusterMaxCellDiffWeightNoMaxTime ->Fill(clus->E(),tmax-weightedTimeNoMax);
-            //printf("E %f, tmax %f, aver %f, weigh %f, averNoMax %f, weightNoMax %f\n",
-            //       clus->E(),tmax,averTime,weightedTime,averTimeNoMax,weightedTimeNoMax);
-          }           
-
-          if( weight > 0 ) fhBadClusterNoMaxCellWeight ->Fill(clus->E(),weightNoMax/weight);
-          //if( weight > 0 && weightNoMax > 0.0)printf("weight %f, weightNoMax %f\n",weight,weightNoMax);
-
-          for (Int_t ipos = 0; ipos < nCaloCellsPerCluster; ipos++) {
-            Int_t absId  = indexList[ipos]; 
-            if(absId!=absIdMax){
-              Float_t frac = cell->GetCellAmplitude(absId)/cell->GetCellAmplitude(absIdMax);
-              fhBadClusterMaxCellCloseCellRatio->Fill(mom.E(),frac);
-              fhBadClusterMaxCellCloseCellDiff ->Fill(mom.E(),cell->GetCellAmplitude(absIdMax)-cell->GetCellAmplitude(absId));
-              if(GetReader()->GetDataType()==AliCaloTrackReader::kESD) {
-                Float_t diff = (tmax-cell->GetCellTime(absId)*1e9);
-                fhBadCellTimeSpreadRespectToCellMax->Fill(mom.E(), diff);
+          if(clus->GetM02() > 0.01 && clus2->GetM02() > 0.01) {
+            Double_t tof2   = clus2->GetTOF(); 
+            // approximation in case of time recalibration
+            if(fCalorimeter=="EMCAL" && GetCaloUtils()->GetEMCALRecoUtils()->IsTimeRecalibrationOn()){
+              Float_t maxCellFraction2 = -1;
+              Int_t absIdMax2 = GetCaloUtils()->GetMaxEnergyCell(cell, clus2,maxCellFraction2); 
+              tof2   = cell->GetCellTime(absIdMax2);
+              GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTime(absIdMax2,bc,tof2);
+            }            
+            fhBadClusterPairDiffTimeE  ->Fill(clus->E(), tof-tof2*1.e9);
+          }
+        } // loop
+        
+        // Max cell compared to other cells in cluster
+        if(GetReader()->GetDataType()==AliCaloTrackReader::kESD) {
+          fhBadClusterMaxCellDiffAverageTime     ->Fill(clus->E(),tmax-averTime);
+          fhBadClusterMaxCellDiffWeightTime      ->Fill(clus->E(),tmax-weightedTime);
+          fhBadClusterDiffWeightAverTime         ->Fill(clus->E(),weightedTime-averTime);
+          fhBadClusterMaxCellDiffAverageNoMaxTime->Fill(clus->E(),tmax-averTimeNoMax);
+          fhBadClusterMaxCellDiffWeightNoMaxTime ->Fill(clus->E(),tmax-weightedTimeNoMax);
+          //printf("E %f, tmax %f, aver %f, weigh %f, averNoMax %f, weightNoMax %f\n",
+          //       clus->E(),tmax,averTime,weightedTime,averTimeNoMax,weightedTimeNoMax);
+        }           
+        
+        if( weight > 0 ) fhBadClusterNoMaxCellWeight ->Fill(clus->E(),weightNoMax/weight);
+        //if( weight > 0 && weightNoMax > 0.0)printf("weight %f, weightNoMax %f\n",weight,weightNoMax);
+        
+        for (Int_t ipos = 0; ipos < nCaloCellsPerCluster; ipos++) {
+          Int_t absId  = indexList[ipos]; 
+          if(absId!=absIdMax){
+            Float_t frac = cell->GetCellAmplitude(absId)/cell->GetCellAmplitude(absIdMax);
+            fhBadClusterMaxCellCloseCellRatio->Fill(mom.E(),frac);
+            fhBadClusterMaxCellCloseCellDiff ->Fill(mom.E(),cell->GetCellAmplitude(absIdMax)-cell->GetCellAmplitude(absId));
+            time      = cell->GetCellTime(absId);
+            if(GetReader()->GetDataType()==AliCaloTrackReader::kESD) {
+              
+              if(fCalorimeter=="EMCAL" && GetCaloUtils()->GetEMCALRecoUtils()->IsTimeRecalibrationOn()){
+                GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTime(absId,bc,time);
               }
-            }// Not max
-          }//loop
-        }//Bad cluster
-      }
+              Float_t diff = (tmax-time*1e9);
+              fhBadCellTimeSpreadRespectToCellMax->Fill(mom.E(), diff);
+              
+            } // ESD
+          }// Not max
+        }//loop
+      }//Bad cluster
+    }
+    
+    if(!badCluster){
+      
+      fhClusterMaxCellDiff->Fill(clus->E(),maxCellFraction);        
+      fhClusterTimeEnergy ->Fill(mom.E(),tof);
       
-      if(!badCluster){
-                
-        fhClusterMaxCellDiff->Fill(clus->E(),maxCellFraction);        
-        fhClusterTimeEnergy ->Fill(mom.E(),tof);
+      //Clusters in event time difference
+      for(Int_t iclus2 = 0; iclus2 < nCaloClusters; iclus2++ ){
         
-        //Clusters in event time difference
-        for(Int_t iclus2 = 0; iclus2 < nCaloClusters; iclus2++ ){
+        AliVCluster* clus2 =  (AliVCluster*) caloClusters->At(iclus2);
+        
+        if(clus->GetID()==clus2->GetID()) continue;
+        
+        if(clus->GetM02() > 0.01 && clus2->GetM02() > 0.01) {
           
-          AliVCluster* clus2 =  (AliVCluster*) caloClusters->At(iclus2);
+          Double_t tof2   = clus2->GetTOF(); 
+          // approximation in case of time recalibration
+          if(fCalorimeter=="EMCAL" && GetCaloUtils()->GetEMCALRecoUtils()->IsTimeRecalibrationOn()){
+            Float_t maxCellFraction2 = -1;
+            Int_t absIdMax2 = GetCaloUtils()->GetMaxEnergyCell(cell, clus2,maxCellFraction2);
+            tof2   = cell->GetCellTime(absIdMax2);
+            GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTime(absIdMax2,bc,tof2);
+          }            
           
-          if(clus->GetID()==clus2->GetID()) continue;
+          fhClusterPairDiffTimeE  ->Fill(clus->E(), tof-tof2*1.e9);
+        }
+      }        
+      
+      if(nCaloCellsPerCluster > 1){
+        
+        // check time of cells respect to max energy cell
+        
+        if(GetReader()->GetDataType()==AliCaloTrackReader::kESD) {
+          fhClusterMaxCellDiffAverageTime     ->Fill(clus->E(),tmax-averTime);
+          fhClusterMaxCellDiffWeightTime      ->Fill(clus->E(),tmax-weightedTime);
+          fhClusterDiffWeightAverTime         ->Fill(clus->E(), weightedTime-averTime);
+          fhClusterMaxCellDiffAverageNoMaxTime->Fill(clus->E(),tmax-averTimeNoMax);
+          fhClusterMaxCellDiffWeightNoMaxTime ->Fill(clus->E(),tmax-weightedTimeNoMax);
           
-          if(clus->GetM02() > 0.01) {
-            fhClusterPairDiffTimeE  ->Fill(clus->E(), tof-clus2->GetTOF()*1.e9);
-          }
-        }        
+        }
         
-        if(nCaloCellsPerCluster > 1){
+        if( weight > 0 )fhClusterNoMaxCellWeight ->Fill(clus->E(),weightNoMax/weight);
+        
+        for (Int_t ipos = 0; ipos < nCaloCellsPerCluster; ipos++) {
           
-          // check time of cells respect to max energy cell
+          Int_t absId  = indexList[ipos];             
+          if(absId == absIdMax) continue;
+          
+          Float_t frac = cell->GetCellAmplitude(absId)/cell->GetCellAmplitude(absIdMax);            
+          fhClusterMaxCellCloseCellRatio->Fill(mom.E(),frac);
+          fhClusterMaxCellCloseCellDiff ->Fill(mom.E(),cell->GetCellAmplitude(absIdMax)-cell->GetCellAmplitude(absId));
           
           if(GetReader()->GetDataType()==AliCaloTrackReader::kESD) {
-            fhClusterMaxCellDiffAverageTime->Fill(clus->E(),tmax-averTime);
-            fhClusterMaxCellDiffWeightTime ->Fill(clus->E(),tmax-weightedTime);
-            fhClusterDiffWeightAverTime ->Fill(clus->E(), weightedTime-averTime);
-            fhClusterMaxCellDiffAverageNoMaxTime->Fill(clus->E(),tmax-averTimeNoMax);
-            fhClusterMaxCellDiffWeightNoMaxTime ->Fill(clus->E(),tmax-weightedTimeNoMax);
+            
+            time      = cell->GetCellTime(absId);
+            if(fCalorimeter=="EMCAL" && GetCaloUtils()->GetEMCALRecoUtils()->IsTimeRecalibrationOn()){
+              GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTime(absId,bc,time);
+            }
+            
+            Float_t diff = (tmax-time*1.0e9);
+            fhCellTimeSpreadRespectToCellMax->Fill(mom.E(), diff);
+            if(TMath::Abs(TMath::Abs(diff) > 100) && mom.E() > 1 ) fhCellIdCellLargeTimeSpread->Fill(absId);
           }
           
-          if( weight > 0 )fhClusterNoMaxCellWeight ->Fill(clus->E(),weightNoMax/weight);
-
-          for (Int_t ipos = 0; ipos < nCaloCellsPerCluster; ipos++) {
+          Int_t sm =0; Int_t ietaa=-1; Int_t iphii = 0; Int_t rcu = 0;
+          sm = GetModuleNumberCellIndexes(absId,fCalorimeter, ietaa, iphii, rcu);
+          if(dIphi < TMath::Abs(iphii-iphiiMax)) dIphi = TMath::Abs(iphii-iphiiMax);
+          if(smMax==sm){
+            if(dIeta < TMath::Abs(ietaa-ietaaMax)) dIeta = TMath::Abs(ietaa-ietaaMax);
+          }
+          else {
+            Int_t ietaaShift    = ietaa;
+            Int_t ietaaMaxShift = ietaaMax;
+            if (ietaa > ietaaMax)  ietaaMaxShift+=48;
+            else                   ietaaShift   +=48;
+            if(dIeta < TMath::Abs(ietaaShift-ietaaMaxShift)) dIeta = TMath::Abs(ietaaShift-ietaaMaxShift);
+          }
+          
+          //if(TMath::Abs(clus->GetM20()) < 0.0001 && clus->GetNCells() > 3){
+          //  printf("Good : E %f, mcells %d, l0 %f, l1 %f, d %f, cell max t %f, cluster TOF %f, sm %d, icol %d, irow %d; Max icol %d, irow %d \n", 
+          //            clus->E(), clus->GetNCells(),clus->GetM02(), clus->GetM20(), clus->GetDispersion(),tmax, tof,sm,ietaa,iphii, ietaaMax, iphiiMax);
+          //}
+          
+          
+        }// fill cell-cluster histogram loop
+        
+        if(nCaloCellsPerCluster > 3){
+          
+          if     (mom.E() < 2 ) fhDeltaIEtaDeltaIPhiE0[matched]->Fill(dIeta,dIphi);
+          else if(mom.E() < 6 ) fhDeltaIEtaDeltaIPhiE2[matched]->Fill(dIeta,dIphi);
+          else                  fhDeltaIEtaDeltaIPhiE6[matched]->Fill(dIeta,dIphi);
+          
+          Float_t dIA    = 1.*(dIphi-dIeta)/(dIeta+dIphi);
+          fhDeltaIA[matched]->Fill(mom.E(),dIA);
+          
+          if(mom.E() > 0.5){
             
-            Int_t absId  = indexList[ipos];             
-            if(absId == absIdMax) continue;
+            fhDeltaIAL0[matched]->Fill(clus->GetM02(),dIA);
+            fhDeltaIAL1[matched]->Fill(clus->GetM20(),dIA);
+            fhDeltaIANCells[matched]->Fill(nCaloCellsPerCluster,dIA);
             
-            Float_t frac = cell->GetCellAmplitude(absId)/cell->GetCellAmplitude(absIdMax);            
-            fhClusterMaxCellCloseCellRatio->Fill(mom.E(),frac);
-            fhClusterMaxCellCloseCellDiff ->Fill(mom.E(),cell->GetCellAmplitude(absIdMax)-cell->GetCellAmplitude(absId));
-
-            if(GetReader()->GetDataType()==AliCaloTrackReader::kESD) {
-              Float_t diff = (tmax-cell->GetCellTime(absId)*1e9);
-              fhCellTimeSpreadRespectToCellMax->Fill(mom.E(), diff);
-              if(TMath::Abs(TMath::Abs(diff) > 100) && mom.E() > 1 ) fhCellIdCellLargeTimeSpread->Fill(absId);
+          }
+          
+          if(IsDataMC()){
+            Int_t tag = GetMCAnalysisUtils()->CheckOrigin(labels,nLabel, GetReader(),0);
+            if( (GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton) || 
+                 GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0)    || 
+                 GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta)         ) &&
+               !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)        ){
+              fhDeltaIAMC[0]->Fill(mom.E(),dIA);//Pure Photon
             }
-            
-            Int_t sm =0; Int_t ietaa=-1; Int_t iphii = 0; Int_t rcu = 0;
-            sm = GetModuleNumberCellIndexes(absId,fCalorimeter, ietaa, iphii, rcu);
-            if(dIphi < TMath::Abs(iphii-iphiiMax)) dIphi = TMath::Abs(iphii-iphiiMax);
-            if(smMax==sm){
-              if(dIeta < TMath::Abs(ietaa-ietaaMax)) dIeta = TMath::Abs(ietaa-ietaaMax);
+            else if ( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCElectron) &&
+                     !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)  ){
+              fhDeltaIAMC[1]->Fill(mom.E(),dIA);//Pure electron
             }
-            else {
-              Int_t ietaaShift    = ietaa;
-              Int_t ietaaMaxShift = ietaaMax;
-              if (ietaa > ietaaMax)  ietaaMaxShift+=48;
-              else                   ietaaShift   +=48;
-              if(dIeta < TMath::Abs(ietaaShift-ietaaMaxShift)) dIeta = TMath::Abs(ietaaShift-ietaaMaxShift);
+            else if ( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)  ){
+              fhDeltaIAMC[2]->Fill(mom.E(),dIA);//Converted cluster
             }
-            
-            //if(TMath::Abs(clus->GetM20()) < 0.0001 && clus->GetNCells() > 3){
-            //  printf("Good : E %f, mcells %d, l0 %f, l1 %f, d %f, cell max t %f, cluster TOF %f, sm %d, icol %d, irow %d; Max icol %d, irow %d \n", 
-            //            clus->E(), clus->GetNCells(),clus->GetM02(), clus->GetM20(), clus->GetDispersion(),tmax, tof,sm,ietaa,iphii, ietaaMax, iphiiMax);
-            //}
-            
-            
-          }// fill cell-cluster histogram loop
-
-          if(nCaloCellsPerCluster > 3){
-            
-            if     (mom.E() < 2 ) fhDeltaIEtaDeltaIPhiE0[matched]->Fill(dIeta,dIphi);
-            else if(mom.E() < 6 ) fhDeltaIEtaDeltaIPhiE2[matched]->Fill(dIeta,dIphi);
-            else                  fhDeltaIEtaDeltaIPhiE6[matched]->Fill(dIeta,dIphi);
-            
-            Float_t dIA    = 1.*(dIphi-dIeta)/(dIeta+dIphi);
-            fhDeltaIA[matched]->Fill(mom.E(),dIA);
-            
-            if(mom.E() > 0.5){
-                     
-              fhDeltaIAL0[matched]->Fill(clus->GetM02(),dIA);
-              fhDeltaIAL1[matched]->Fill(clus->GetM20(),dIA);
-              fhDeltaIANCells[matched]->Fill(nCaloCellsPerCluster,dIA);
-              
+            else{ 
+              fhDeltaIAMC[3]->Fill(mom.E(),dIA);//Hadrons
             }
             
-            if(IsDataMC()){
-              Int_t tag = GetMCAnalysisUtils()->CheckOrigin(labels,nLabel, GetReader(),0);
-              if( (GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton) || 
-                   GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0)    || 
-                   GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta)         ) &&
-                 !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)        ){
-                fhDeltaIAMC[0]->Fill(mom.E(),dIA);//Pure Photon
-              }
-              else if ( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCElectron) &&
-                       !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)  ){
-                fhDeltaIAMC[1]->Fill(mom.E(),dIA);//Pure electron
-              }
-              else if ( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)  ){
-                fhDeltaIAMC[2]->Fill(mom.E(),dIA);//Converted cluster
-              }
-              else{ 
-                fhDeltaIAMC[3]->Fill(mom.E(),dIA);//Hadrons
-              }
-              
-            }
-          }// 4 cells at least in cluster for size study
-          
-        }//check time and energy of cells respect to max energy cell if cluster of more than 1 cell
-        
+          }
+        }// 4 cells at least in cluster for size study
         
-        //Get module of cluster
-
-        nCaloClustersAccepted++;
-        if(nModule >=0 && nModule < fNModules) {
-          if     (fCalorimeter=="EMCAL" && mom.E() > 2*fEMCALCellAmpMin)  nClustersInModule[nModule]++;
-          else if(fCalorimeter=="PHOS"  && mom.E() > 2*fPHOSCellAmpMin )  nClustersInModule[nModule]++;
+      }//check time and energy of cells respect to max energy cell if cluster of more than 1 cell
+      
+      
+      //Get module of cluster
+      
+      nCaloClustersAccepted++;
+      if(nModule >=0 && nModule < fNModules) {
+        if     (fCalorimeter=="EMCAL" && mom.E() > 2*fEMCALCellAmpMin)  nClustersInModule[nModule]++;
+        else if(fCalorimeter=="PHOS"  && mom.E() > 2*fPHOSCellAmpMin )  nClustersInModule[nModule]++;
+      }
+      
+      //-----------------------------------------------------------
+      //Fill histograms related to single cluster or track matching
+      //-----------------------------------------------------------
+      
+      // Get Track matched
+      if(matched){
+      
+        if(!strcmp("AliESDCaloCluster",Form("%s",clus->ClassName())))
+        { 
+          track = dynamic_cast<AliVTrack*> (GetReader()->GetInputEvent()->GetTrack(clus->GetTrackMatchedIndex())); 
         }
+        else //AOD
+        {
+          track = dynamic_cast<AliVTrack*> (clus->GetTrackMatched(0)); 
+        }
+      }
+      
+      ClusterHistograms(mom, pos, nCaloCellsPerCluster, nModule, matched,track, labels, nLabel);       
+      
+      
+      //-----------------------------------------------------------
+      //Invariant mass
+      //-----------------------------------------------------------
+      if(fFillAllPi0Histo){
+        if(GetDebug()>1) printf("Invariant mass \n");
         
-        //-----------------------------------------------------------
-        //Fill histograms related to single cluster or track matching
-        //-----------------------------------------------------------
-        
-        ClusterHistograms(mom, pos, nCaloCellsPerCluster, nModule, matched,0x0/*track*/, labels, nLabel);      
-        
-        
-        //-----------------------------------------------------------
-        //Invariant mass
-        //-----------------------------------------------------------
-        if(fFillAllPi0Histo){
-          if(GetDebug()>1) printf("Invariant mass \n");
-          
-          //do not do for bad vertex
-          // Float_t fZvtxCut = 40. ;  
-          if(v[2]<-GetZvertexCut() || v[2]> GetZvertexCut()) continue ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
-          
-          Int_t nModule2 = -1;
-          if (nCaloClusters > 1 && nCaloCellsPerCluster > 1) {
-            for(Int_t jclus = iclus + 1 ; jclus < nCaloClusters ; jclus++) {
-              AliVCluster* clus2 =  (AliVCluster*)caloClusters->At(jclus);
-              if( clus2->GetNCells() > 1 ){
-                
-                //Get cluster kinematics
-                clus2->GetMomentum(mom2,v);
-                
-                //Check only certain regions
-                Bool_t in2 = kTRUE;
-                if(IsFiducialCutOn()) in2 =  GetFiducialCut()->IsInFiducialCut(mom2,fCalorimeter) ;
-                if(!in2) continue;     
-                
-                //Get module of cluster
-                nModule2 = GetModuleNumber(clus2);
-                
-                //Fill invariant mass histograms
-                
-                //All modules
-                fhIM  ->Fill((mom+mom2).Pt(),(mom+mom2).M());
-                
-                //Single module
-                if(nModule == nModule2 && nModule >=0 && nModule < fNModules)
-                  fhIMMod[nModule]->Fill((mom+mom2).Pt(),(mom+mom2).M());
-                
-                //Asymetry histograms
-                fhAsym->Fill((mom+mom2).Pt(),TMath::Abs((mom.E()-mom2.E())/(mom.E()+mom2.E())));
-              } 
-            }// 2nd cluster loop
-          } // At least one cell in cluster and one cluster in the array
-        }//Fill Pi0
+        //do not do for bad vertex
+        // Float_t fZvtxCut = 40. ;    
+        if(v[2]<-GetZvertexCut() || v[2]> GetZvertexCut()) continue ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
         
-      }//good cluster
+        Int_t nModule2 = -1;
+        if (nCaloClusters > 1 && nCaloCellsPerCluster > 1) {
+          for(Int_t jclus = iclus + 1 ; jclus < nCaloClusters ; jclus++) {
+            AliVCluster* clus2 =  (AliVCluster*)caloClusters->At(jclus);
+            if( clus2->GetNCells() > 1 ){
+              
+              //Get cluster kinematics
+              clus2->GetMomentum(mom2,v);
+              
+              //Check only certain regions
+              Bool_t in2 = kTRUE;
+              if(IsFiducialCutOn()) in2 =  GetFiducialCut()->IsInFiducialCut(mom2,fCalorimeter) ;
+              if(!in2) continue;       
+              
+              //Get module of cluster
+              nModule2 = GetModuleNumber(clus2);
+              
+              //Fill invariant mass histograms
+              
+              //All modules
+              fhIM  ->Fill((mom+mom2).Pt(),(mom+mom2).M());
+              
+              //Single module
+              if(nModule == nModule2 && nModule >=0 && nModule < fNModules)
+                fhIMMod[nModule]->Fill((mom+mom2).Pt(),(mom+mom2).M());
+              
+              //Asymetry histograms
+              fhAsym->Fill((mom+mom2).Pt(),TMath::Abs((mom.E()-mom2.E())/(mom.E()+mom2.E())));
+            } 
+          }// 2nd cluster loop
+        } // At least one cell in cluster and one cluster in the array
+      }//Fill Pi0
       
-    }//cluster loop
+    }//good cluster
     
-    //Number of clusters histograms
-    if(nCaloClustersAccepted > 0) fhNClusters->Fill(nCaloClustersAccepted);
-    //  Number of clusters per module
-    for(Int_t imod = 0; imod < fNModules; imod++ ){ 
-      if(GetDebug() > 1) 
-        printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - module %d calo %s clusters %d\n", imod, fCalorimeter.Data(), nClustersInModule[imod]); 
-      fhNClustersMod->Fill(nClustersInModule[imod],imod);
-    }
-    delete [] nClustersInModule;
-    //delete caloClusters;
-  }// calo clusters array exists
+  }//cluster loop
+  
+  //Number of clusters histograms
+  if(nCaloClustersAccepted > 0) fhNClusters->Fill(nCaloClustersAccepted);
+  //  Number of clusters per module
+  for(Int_t imod = 0; imod < fNModules; imod++ ){ 
+    if(GetDebug() > 1) 
+      printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - module %d calo %s clusters %d\n", imod, fCalorimeter.Data(), nClustersInModule[imod]); 
+    fhNClustersMod->Fill(nClustersInModule[imod],imod);
+  }
+  delete [] nClustersInModule;
+  //delete caloClusters;
   
   //----------------------------------------------------------
   // CALOCELLS
   //----------------------------------------------------------
-  
-  AliVCaloCells * cell = 0x0; 
-  Int_t ncells = 0;
-  if(fCalorimeter == "PHOS") 
-    cell = GetPHOSCells();
-  else                       
-    cell = GetEMCALCells();
-  
-  if(!cell){ 
-    AliFatal(Form("No %s CELLS available for analysis",fCalorimeter.Data()));
-    return; // just to trick coverity
-  }
+    
+  Int_t ncells = cell->GetNumberOfCells();
   
   if(GetDebug() > 0) 
-    printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - %s cell entries %d\n", fCalorimeter.Data(), cell->GetNumberOfCells());    
+    printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - %s cell entries %d\n", fCalorimeter.Data(), ncells );    
   
   //Init arrays and used variables
   Int_t *nCellsInModule = new Int_t[fNModules];
   for(Int_t imod = 0; imod < fNModules; imod++ ) nCellsInModule[imod] = 0;
-  Int_t icol     = -1;
-  Int_t irow     = -1;
-  Int_t iRCU     = -1;
-  Float_t amp    = 0.;
-  Float_t time   = 0.;
-  Int_t id       = -1;
-  Float_t recalF = 1.;  
-  
+  Int_t    icol   = -1;
+  Int_t    irow   = -1;
+  Int_t    iRCU   = -1;
+  Float_t  amp    = 0.;
+  Double_t time   = 0.;
+  Int_t    id     = -1;
+  Float_t  recalF = 1.;  
+
   for (Int_t iCell = 0; iCell < cell->GetNumberOfCells(); iCell++) {      
     if(GetDebug() > 2)  
       printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - Cell : amp %f, absId %d \n", cell->GetAmplitude(iCell), cell->GetCellNumber(iCell));
@@ -1733,27 +1801,40 @@ void  AliAnaCalorimeterQA::MakeAnalysisFillHistograms()
         }
       } // use bad channel map
       
+      amp     = cell->GetAmplitude(iCell)*recalF;
+      time    = cell->GetTime(iCell);
+      id      = cell->GetCellNumber(iCell);
+      
       //Get Recalibration factor if set
       if (GetCaloUtils()->IsRecalibrationOn()) {
-        if(fCalorimeter == "PHOS") recalF = GetCaloUtils()->GetPHOSChannelRecalibrationFactor(nModule,icol,irow);
-        else                              recalF = GetCaloUtils()->GetEMCALChannelRecalibrationFactor(nModule,icol,irow);
+        if(fCalorimeter == "PHOS") {
+          amp *= GetCaloUtils()->GetPHOSChannelRecalibrationFactor(nModule,icol,irow);
+        }
+        else                              {
+          amp *= GetCaloUtils()->GetEMCALChannelRecalibrationFactor(nModule,icol,irow);
+          GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTime(id,bc,time);
+//          printf("QA         : Id %d, Time org %e, Time new %e; Amp org %f, Amp new %f\n",
+//                          id, cell->GetTime(iCell),time, cell->GetAmplitude(iCell),amp);
+        }
         //if(fCalorimeter == "PHOS")printf("Recalibration factor (sm,row,col)=(%d,%d,%d) -  %f\n",nModule,icol,irow,recalF);
       }
       
-      amp     = cell->GetAmplitude(iCell)*recalF;
-      time    = cell->GetTime(iCell)*1e9;//transform time to ns
+      //Transform time to ns
+      time *= 1.0e9;
       
       //Remove noisy channels, only possible in ESDs
       if(GetReader()->GetDataType() == AliCaloTrackReader::kESD){
-        if(time < fTimeCutMin || time > fTimeCutMax) continue;
+        if(time < fTimeCutMin || time > fTimeCutMax){
+          printf("Remove cluster with TOF %f\n",time);
+          continue;
+        }
       }
       //if(amp > 3 && fCalorimeter=="EMCAL") printf("Amp = %f, time = %f, (mod, col, row)= (%d,%d,%d)\n",
       //                                                                                  amp,time,nModule,icol,irow);
       
-      id      = cell->GetCellNumber(iCell);
       fhAmplitude->Fill(amp);
       fhAmpId    ->Fill(amp,id);
-            
+      
       if ((fCalorimeter=="EMCAL" && amp > fEMCALCellAmpMin) ||
           (fCalorimeter=="PHOS"  && amp > fPHOSCellAmpMin))   {
         
@@ -1768,7 +1849,7 @@ void  AliAnaCalorimeterQA::MakeAnalysisFillHistograms()
         else {
           irows = irow + fNMaxRows * fNModules;
         }
-
+        
         fhGridCellsMod ->Fill(icols,irows);
         fhGridCellsEMod->Fill(icols,irows,amp);
         
@@ -1778,9 +1859,9 @@ void  AliAnaCalorimeterQA::MakeAnalysisFillHistograms()
           fhTimeId   ->Fill(time,id);
           fhTimeAmp  ->Fill(amp,time);
           fhGridCellsTimeMod->Fill(icols,irows,time);
-
+          
           fhTimeAmpPerRCU  [nModule*fNRCU+iRCU]->Fill(amp, time);
-    
+          
         }
       }
       
@@ -1845,7 +1926,7 @@ void  AliAnaCalorimeterQA::MakeAnalysisFillHistograms()
 //_____________________________________________________________________________________________
 void AliAnaCalorimeterQA::ClusterHistograms(const TLorentzVector mom, 
                                             Float_t *pos, const Int_t nCaloCellsPerCluster,const Int_t nModule,
-                                            const Bool_t matched,  const AliVTrack * /*track*/,  
+                                            const Bool_t matched,  const AliVTrack * track,  
                                             const Int_t * labels, const Int_t nLabels){
   //Fill CaloCluster related histograms
        
@@ -2158,118 +2239,118 @@ void AliAnaCalorimeterQA::ClusterHistograms(const TLorentzVector mom,
     }
     //printf("track index %d ntracks %d\n", esd->GetNumberOfTracks());
     
-//    //Study the track and matched cluster if track exists.
-//    if(!track) return;
-//    Double_t emcpos[3] = {0.,0.,0.};
-//    Double_t emcmom[3] = {0.,0.,0.};
-//    Double_t radius    = 441.0; //[cm] EMCAL radius +13cm
-//    Double_t bfield    = 0.;
-//    Double_t tphi      = 0;
-//    Double_t teta      = 0;
-//    Double_t tmom      = 0;
-//    Double_t tpt       = 0;
-//    Double_t tmom2     = 0;
-//    Double_t tpcSignal = 0;
-//    Bool_t okpos = kFALSE;
-//    Bool_t okmom = kFALSE;
-//    Bool_t okout = kFALSE;
-//    Int_t nITS   = 0;
-//    Int_t nTPC   = 0;
-//    
-//    //In case of ESDs get the parameters in this way
-//    if(GetReader()->GetDataType()==AliCaloTrackReader::kESD) {
-//      if (track->GetOuterParam() ) {
-//        okout = kTRUE;
-//        
-//        bfield = GetReader()->GetInputEvent()->GetMagneticField();
-//        okpos = track->GetOuterParam()->GetXYZAt(radius,bfield,emcpos);
-//        okmom = track->GetOuterParam()->GetPxPyPzAt(radius,bfield,emcmom);
-//        if(!(okpos && okmom)) return;
-//        
-//        TVector3 position(emcpos[0],emcpos[1],emcpos[2]);
-//        TVector3 momentum(emcmom[0],emcmom[1],emcmom[2]);
-//        tphi = position.Phi();
-//        teta = position.Eta();
-//        tmom = momentum.Mag();
-//        
-//        tpt       = track->Pt();
-//        tmom2     = track->P();
-//        tpcSignal = track->GetTPCsignal();
-//        
-//        nITS = track->GetNcls(0);
-//        nTPC = track->GetNcls(1);
-//      }//Outer param available 
-//    }// ESDs
-//    else if(GetReader()->GetDataType()==AliCaloTrackReader::kAOD) {
-//      AliAODPid* pid = (AliAODPid*) ((AliAODTrack *) track)->GetDetPid();
-//      if (pid) {
-//        okout = kTRUE;
-//        pid->GetEMCALPosition(emcpos);
-//        pid->GetEMCALMomentum(emcmom);       
-//        
-//        TVector3 position(emcpos[0],emcpos[1],emcpos[2]);
-//        TVector3 momentum(emcmom[0],emcmom[1],emcmom[2]);
-//        tphi = position.Phi();
-//        teta = position.Eta();
-//        tmom = momentum.Mag();
-//        
-//        tpt       = track->Pt();
-//        tmom2     = track->P();
-//        tpcSignal = pid->GetTPCsignal();
-//        
-//       }//pid 
-//    }//AODs
-//             
-//    if(okout){
-//      printf("okout\n");
-//      Double_t deta = teta - eta;
-//      Double_t dphi = tphi - phi;
-//      if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();
-//      if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();
-//      Double_t dR = sqrt(dphi*dphi + deta*deta);
-//                     
-//      Double_t pOverE = tmom/e;
-//                     
-//      fh1pOverE->Fill(tpt, pOverE);
-//      if(dR < 0.02) fh1pOverER02->Fill(tpt,pOverE);
-//                     
-//      fh1dR->Fill(dR);
-//      fh2MatchdEdx->Fill(tmom2,tpcSignal);
-//                     
-//      if(IsDataMC() && primary){ 
-//        Int_t pdg = primary->GetPdgCode();
-//        Double_t  charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
-//                             
-//        if(TMath::Abs(pdg) == 11){
-//          fhMCEle1pOverE->Fill(tpt,pOverE);
-//          fhMCEle1dR->Fill(dR);
-//          fhMCEle2MatchdEdx->Fill(tmom2,tpcSignal);          
-//          if(dR < 0.02) fhMCEle1pOverER02->Fill(tpt,pOverE);
-//        }
-//        else if(charge!=0){
-//          fhMCChHad1pOverE->Fill(tpt,pOverE);
-//          fhMCChHad1dR->Fill(dR);
-//          fhMCChHad2MatchdEdx->Fill(tmom2,tpcSignal);        
-//          if(dR < 0.02) fhMCChHad1pOverER02->Fill(tpt,pOverE);
-//        }
-//        else if(charge == 0){
-//          fhMCNeutral1pOverE->Fill(tpt,pOverE);
-//          fhMCNeutral1dR->Fill(dR);
-//          fhMCNeutral2MatchdEdx->Fill(tmom2,tpcSignal);      
-//          if(dR < 0.02) fhMCNeutral1pOverER02->Fill(tpt,pOverE);
-//        }
-//      }//DataMC
-//      
-//      if(dR < 0.02 && pOverE > 0.5 && pOverE < 1.5
-//         && nCaloCellsPerCluster > 1 && nITS > 3 && nTPC > 20) {
-//        fh2EledEdx->Fill(tmom2,tpcSignal);
-//      }
-//    }
-//    else{//no ESD external param or AODPid
-//
-//      if(GetDebug() >= 0) printf("No ESD external param or AliAODPid \n");
-//      
-//    }//No out params
+    //Study the track and matched cluster if track exists.
+    if(!track) return;
+    Double_t emcpos[3] = {0.,0.,0.};
+    Double_t emcmom[3] = {0.,0.,0.};
+    Double_t radius    = 441.0; //[cm] EMCAL radius +13cm
+    Double_t bfield    = 0.;
+    Double_t tphi      = 0;
+    Double_t teta      = 0;
+    Double_t tmom      = 0;
+    Double_t tpt       = 0;
+    Double_t tmom2     = 0;
+    Double_t tpcSignal = 0;
+    Bool_t okpos = kFALSE;
+    Bool_t okmom = kFALSE;
+    Bool_t okout = kFALSE;
+    Int_t nITS   = 0;
+    Int_t nTPC   = 0;
+    
+    //In case of ESDs get the parameters in this way
+    if(GetReader()->GetDataType()==AliCaloTrackReader::kESD) {
+      if (track->GetOuterParam() ) {
+        okout = kTRUE;
+        
+        bfield = GetReader()->GetInputEvent()->GetMagneticField();
+        okpos = track->GetOuterParam()->GetXYZAt(radius,bfield,emcpos);
+        okmom = track->GetOuterParam()->GetPxPyPzAt(radius,bfield,emcmom);
+        if(!(okpos && okmom)) return;
+        
+        TVector3 position(emcpos[0],emcpos[1],emcpos[2]);
+        TVector3 momentum(emcmom[0],emcmom[1],emcmom[2]);
+        tphi = position.Phi();
+        teta = position.Eta();
+        tmom = momentum.Mag();
+        
+        tpt       = track->Pt();
+        tmom2     = track->P();
+        tpcSignal = track->GetTPCsignal();
+        
+        nITS = track->GetNcls(0);
+        nTPC = track->GetNcls(1);
+      }//Outer param available 
+    }// ESDs
+    else if(GetReader()->GetDataType()==AliCaloTrackReader::kAOD) {
+      AliAODPid* pid = (AliAODPid*) ((AliAODTrack *) track)->GetDetPid();
+      if (pid) {
+        okout = kTRUE;
+        pid->GetEMCALPosition(emcpos);
+        pid->GetEMCALMomentum(emcmom); 
+        
+        TVector3 position(emcpos[0],emcpos[1],emcpos[2]);
+        TVector3 momentum(emcmom[0],emcmom[1],emcmom[2]);
+        tphi = position.Phi();
+        teta = position.Eta();
+        tmom = momentum.Mag();
+        
+        tpt       = track->Pt();
+        tmom2     = track->P();
+        tpcSignal = pid->GetTPCsignal();
+        
+       }//pid 
+    }//AODs
+               
+    if(okout){
+      //printf("okout\n");
+      Double_t deta = teta - eta;
+      Double_t dphi = tphi - phi;
+      if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();
+      if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();
+      Double_t dR = sqrt(dphi*dphi + deta*deta);
+                       
+      Double_t pOverE = tmom/e;
+                       
+      fh1pOverE->Fill(tpt, pOverE);
+      if(dR < 0.02) fh1pOverER02->Fill(tpt,pOverE);
+                       
+      fh1dR->Fill(dR);
+      fh2MatchdEdx->Fill(tmom2,tpcSignal);
+                       
+      if(IsDataMC() && primary){ 
+        Int_t pdg = primary->GetPdgCode();
+        Double_t  charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
+                               
+        if(TMath::Abs(pdg) == 11){
+          fhMCEle1pOverE->Fill(tpt,pOverE);
+          fhMCEle1dR->Fill(dR);
+          fhMCEle2MatchdEdx->Fill(tmom2,tpcSignal);            
+          if(dR < 0.02) fhMCEle1pOverER02->Fill(tpt,pOverE);
+        }
+        else if(charge!=0){
+          fhMCChHad1pOverE->Fill(tpt,pOverE);
+          fhMCChHad1dR->Fill(dR);
+          fhMCChHad2MatchdEdx->Fill(tmom2,tpcSignal);  
+          if(dR < 0.02) fhMCChHad1pOverER02->Fill(tpt,pOverE);
+        }
+        else if(charge == 0){
+          fhMCNeutral1pOverE->Fill(tpt,pOverE);
+          fhMCNeutral1dR->Fill(dR);
+          fhMCNeutral2MatchdEdx->Fill(tmom2,tpcSignal);        
+          if(dR < 0.02) fhMCNeutral1pOverER02->Fill(tpt,pOverE);
+        }
+      }//DataMC
+      
+      if(dR < 0.02 && pOverE > 0.5 && pOverE < 1.5
+         && nCaloCellsPerCluster > 1 && nITS > 3 && nTPC > 20) {
+        fh2EledEdx->Fill(tmom2,tpcSignal);
+      }
+    }
+    else{//no ESD external param or AODPid
+
+      if(GetDebug() >= 0) printf("No ESD external param or AliAODPid \n");
+      
+    }//No out params
   }//matched clusters with tracks
   
 }// Clusters
index 612a790..74f0132 100755 (executable)
@@ -291,27 +291,27 @@ public:
   TH2F *   fhHaR       ;                      //! Hadron distance to vertex vs rec energy  
        
   //Histograms for MC track-matching
-//  TH2F *   fh1pOverE;                         //! p/E for track-cluster matches
-//  TH1F *   fh1dR;                             //! distance between projected track and cluster
-//  TH2F *   fh2EledEdx;                        //! dE/dx vs. momentum for electron candidates
-//  TH2F *   fh2MatchdEdx;                      //! dE/dx vs. momentum for all matches
-//     
-//  TH2F *   fhMCEle1pOverE;                    //! p/E for track-cluster matches, MC electrons
-//  TH1F *   fhMCEle1dR;                        //! distance between projected track and cluster, MC electrons
-//  TH2F *   fhMCEle2MatchdEdx;                 //! dE/dx vs. momentum for all matches, MC electrons   
-//     
-//  TH2F *   fhMCChHad1pOverE;                  //! p/E for track-cluster matches, MC charged hadrons
-//  TH1F *   fhMCChHad1dR;                      //! distance between projected track and cluster, MC charged hadrons
-//  TH2F *   fhMCChHad2MatchdEdx;               //! dE/dx vs. momentum for all matches, MC charged
-//     
-//  TH2F *   fhMCNeutral1pOverE;                //! p/E for track-cluster matches, MC neutral
-//  TH1F *   fhMCNeutral1dR;                    //! distance between projected track and cluster, MC neutral
-//  TH2F *   fhMCNeutral2MatchdEdx;             //! dE/dx vs. momentum for all matches, MC neutral     
-//     
-//  TH2F *   fh1pOverER02;                      //! p/E for track-cluster matches, dR > 0.2    
-//  TH2F *   fhMCEle1pOverER02;                 //! p/E for track-cluster matches, dR > 0.2, MC electrons
-//  TH2F *   fhMCChHad1pOverER02;               //! p/E for track-cluster matches, dR > 0.2, MC charged hadrons
-//  TH2F *   fhMCNeutral1pOverER02;             //! p/E for track-cluster matches, dR > 0.2, MC neutral
+  TH2F *   fh1pOverE;                         //! p/E for track-cluster matches
+  TH1F *   fh1dR;                             //! distance between projected track and cluster
+  TH2F *   fh2EledEdx;                        //! dE/dx vs. momentum for electron candidates
+  TH2F *   fh2MatchdEdx;                      //! dE/dx vs. momentum for all matches
+       
+  TH2F *   fhMCEle1pOverE;                    //! p/E for track-cluster matches, MC electrons
+  TH1F *   fhMCEle1dR;                        //! distance between projected track and cluster, MC electrons
+  TH2F *   fhMCEle2MatchdEdx;                 //! dE/dx vs. momentum for all matches, MC electrons     
+       
+  TH2F *   fhMCChHad1pOverE;                  //! p/E for track-cluster matches, MC charged hadrons
+  TH1F *   fhMCChHad1dR;                      //! distance between projected track and cluster, MC charged hadrons
+  TH2F *   fhMCChHad2MatchdEdx;               //! dE/dx vs. momentum for all matches, MC charged
+       
+  TH2F *   fhMCNeutral1pOverE;                //! p/E for track-cluster matches, MC neutral
+  TH1F *   fhMCNeutral1dR;                    //! distance between projected track and cluster, MC neutral
+  TH2F *   fhMCNeutral2MatchdEdx;             //! dE/dx vs. momentum for all matches, MC neutral       
+       
+  TH2F *   fh1pOverER02;                      //! p/E for track-cluster matches, dR > 0.2      
+  TH2F *   fhMCEle1pOverER02;                 //! p/E for track-cluster matches, dR > 0.2, MC electrons
+  TH2F *   fhMCChHad1pOverER02;               //! p/E for track-cluster matches, dR > 0.2, MC charged hadrons
+  TH2F *   fhMCNeutral1pOverER02;             //! p/E for track-cluster matches, dR > 0.2, MC neutral
        
   ClassDef(AliAnaCalorimeterQA,19)
 } ;