Move exotic cell rejection methods from AliAnalysisTaskEMCALClusterize to AliEMCALRec...
authorgconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 24 Oct 2011 10:04:24 +0000 (10:04 +0000)
committergconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 24 Oct 2011 10:04:24 +0000 (10:04 +0000)
IsExoticCluster now checks if the cell with highest energy is exotic
EMCALPi0CalibSelection: Add time calibration histograms, and bad cluster rejection
AliCaloTrackReader: Add getter for list of clusters under use

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

index ad2866e..04dcd89 100644 (file)
@@ -73,12 +73,14 @@ AliEMCALRecoUtils::AliEMCALRecoUtils():
   fUseRunCorrectionFactors(kFALSE),       fRunCorrectionFactorsSet(kFALSE),
   fRemoveBadChannels(kFALSE),             fRecalDistToBadChannels(kFALSE),        fEMCALBadChannelMap(),
   fNCellsFromEMCALBorder(0),              fNoEMCALBorderAtEta0(kTRUE),
-  fRejectExoticCluster(kFALSE),           fPIDUtils(),                            fAODFilterMask(32),
+  fRejectExoticCluster(kFALSE),           fRejectExoticCells(kFALSE), 
+  fExoticCellFraction(0.97),              fExoticCellDiffTime(1e6),               fExoticCellMinAmplitude(0.5),
+  fPIDUtils(),                            fAODFilterMask(32),
   fMatchedTrackIndex(0x0),                fMatchedClusterIndex(0x0), 
   fResidualEta(0x0), fResidualPhi(0x0),   fCutEtaPhiSum(kTRUE),                   fCutEtaPhiSeparate(kFALSE), 
   fCutR(0.1),                             fCutEta(0.025),                         fCutPhi(0.05),
-  fClusterWindow(100),
-  fMass(0.139),                           fStepSurface(10.),                      fStepCluster(5.),
+  fClusterWindow(100),                    fMass(0.139),                           
+  fStepSurface(10.),                      fStepCluster(5.),
   fTrackCutsType(kLooseCut),              fCutMinTrackPt(0),                      fCutMinNClusterTPC(-1), 
   fCutMinNClusterITS(-1),                 fCutMaxChi2PerClusterTPC(1e10),         fCutMaxChi2PerClusterITS(1e10),
   fCutRequireTPCRefit(kFALSE),            fCutRequireITSRefit(kFALSE),            fCutAcceptKinkDaughters(kFALSE),
@@ -142,8 +144,10 @@ AliEMCALRecoUtils::AliEMCALRecoUtils(const AliEMCALRecoUtils & reco)
   fRemoveBadChannels(reco.fRemoveBadChannels),               fRecalDistToBadChannels(reco.fRecalDistToBadChannels),
   fEMCALBadChannelMap(reco.fEMCALBadChannelMap),
   fNCellsFromEMCALBorder(reco.fNCellsFromEMCALBorder),       fNoEMCALBorderAtEta0(reco.fNoEMCALBorderAtEta0),
-  fRejectExoticCluster(reco.fRejectExoticCluster),           fPIDUtils(reco.fPIDUtils), 
-  fAODFilterMask(reco.fAODFilterMask),
+  fRejectExoticCluster(reco.fRejectExoticCluster),           fRejectExoticCells(reco.fRejectExoticCells), 
+  fExoticCellFraction(reco.fExoticCellFraction),             fExoticCellDiffTime(reco.fExoticCellDiffTime),               
+  fExoticCellMinAmplitude(reco.fExoticCellMinAmplitude),
+  fPIDUtils(reco.fPIDUtils),                                 fAODFilterMask(reco.fAODFilterMask),
   fMatchedTrackIndex(  reco.fMatchedTrackIndex?  new TArrayI(*reco.fMatchedTrackIndex):0x0),
   fMatchedClusterIndex(reco.fMatchedClusterIndex?new TArrayI(*reco.fMatchedClusterIndex):0x0),
   fResidualEta(        reco.fResidualEta?        new TArrayF(*reco.fResidualEta):0x0),
@@ -206,8 +210,13 @@ AliEMCALRecoUtils & AliEMCALRecoUtils::operator = (const AliEMCALRecoUtils & rec
   
   fNCellsFromEMCALBorder     = reco.fNCellsFromEMCALBorder;
   fNoEMCALBorderAtEta0       = reco.fNoEMCALBorderAtEta0;
+  
   fRejectExoticCluster       = reco.fRejectExoticCluster;           
-
+  fRejectExoticCells         = reco.fRejectExoticCells; 
+  fExoticCellFraction        = reco.fExoticCellFraction;
+  fExoticCellDiffTime        = reco.fExoticCellDiffTime;              
+  fExoticCellMinAmplitude    = reco.fExoticCellMinAmplitude;
+  
   fPIDUtils                  = reco.fPIDUtils;
 
   fAODFilterMask             = reco.fAODFilterMask;
@@ -221,7 +230,6 @@ AliEMCALRecoUtils & AliEMCALRecoUtils::operator = (const AliEMCALRecoUtils & rec
   fMass                      = reco.fMass;
   fStepSurface               = reco.fStepSurface;
   fStepCluster               = reco.fStepCluster;
-  fRejectExoticCluster       = reco.fRejectExoticCluster;
 
   fTrackCutsType             = reco.fTrackCutsType;
   fCutMinTrackPt             = reco.fCutMinTrackPt;
@@ -288,7 +296,7 @@ AliEMCALRecoUtils & AliEMCALRecoUtils::operator = (const AliEMCALRecoUtils & rec
 }
 
 
-//__________________________________________________
+//_____________________________________
 AliEMCALRecoUtils::~AliEMCALRecoUtils()
 {
   //Destructor.
@@ -317,6 +325,40 @@ AliEMCALRecoUtils::~AliEMCALRecoUtils()
   InitTrackCuts();
 }
 
+//_______________________________________________________________________________
+Bool_t AliEMCALRecoUtils::AcceptCalibrateCell(const Int_t absID, const Int_t bc,
+                                              Float_t  & amp,    Double_t & time, 
+                                              AliVCaloCells* cells) 
+{
+  // Reject cell if criteria not passed and calibrate it
+  
+  AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
+  
+  if(absID < 0 || absID >= 24*48*geom->GetNumberOfSuperModules()) return kFALSE;
+  
+  Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1; 
+  geom->GetCellIndex(absID,imod,iTower,iIphi,iIeta); 
+  geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);       
+  
+  // Do not include bad channels found in analysis,
+  if( IsBadChannelsRemovalSwitchedOn() && GetEMCALChannelStatus(imod, ieta, iphi)) {
+    return kFALSE;
+  }
+  
+  //Recalibrate energy
+  amp  = cells->GetCellAmplitude(absID);
+  if(IsRecalibrationOn())
+    amp *= GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
+  
+  
+  // Recalibrate time
+  time = cells->GetCellTime(absID);
+  
+  RecalibrateCellTime(absID,bc,time);
+  
+  return kTRUE;
+}
+
 //_______________________________________________________________
 Bool_t AliEMCALRecoUtils::CheckCellFiducialRegion(AliEMCALGeometry* geom, AliVCluster* cluster, AliVCaloCells* cells) 
 {
@@ -415,22 +457,89 @@ Bool_t AliEMCALRecoUtils::ClusterContainsBadChannel(AliEMCALGeometry* geom, USho
        
 }
 
+//_______________________________________________________________________
+Bool_t AliEMCALRecoUtils::IsExoticCell(const Int_t absID, AliVCaloCells* cells, const Int_t bc)
+{
+  // Look to cell neighbourhood and reject if it seems exotic
+  // Do before recalibrating the cells
+  if(!fRejectExoticCells) return kFALSE;
+    
+  AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
+  
+  Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1; 
+  geom->GetCellIndex(absID,imod,iTower,iIphi,iIeta); 
+  geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);       
+  
+  //Get close cells index, energy and time, not in corners
+  Int_t absID1 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi+1, ieta);
+  Int_t absID2 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi-1, ieta);
+  Int_t absID3 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
+  Int_t absID4 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1);
+  
+  Float_t  ecell  = 0, ecell1  = 0, ecell2  = 0, ecell3  = 0, ecell4  = 0;
+  Double_t tcell  = 0, tcell1  = 0, tcell2  = 0, tcell3  = 0, tcell4  = 0;
+  Bool_t   accept = 0, accept1 = 0, accept2 = 0, accept3 = 0, accept4 = 0;
+  
+  accept  = AcceptCalibrateCell(absID, bc, ecell ,tcell ,cells); 
+    
+  if(!accept) return kTRUE; // reject this cell
+  
+  if(ecell < fExoticCellMinAmplitude) return kFALSE; // do not reject low energy cells
+  
+  accept1 = AcceptCalibrateCell(absID1,bc, ecell1,tcell1,cells); 
+  accept2 = AcceptCalibrateCell(absID2,bc, ecell2,tcell2,cells); 
+  accept3 = AcceptCalibrateCell(absID3,bc, ecell3,tcell3,cells); 
+  accept4 = AcceptCalibrateCell(absID4,bc, ecell4,tcell4,cells); 
+  
+  /*
+    printf("Cell absID %d \n",absID);
+    printf("\t  accept1 %d, accept2 %d, accept3 %d, accept4 %d\n",
+           accept1,accept2,accept3,accept4);
+    printf("\t id %d: id1 %d, id2 %d, id3 %d, id4 %d\n",
+           absID,absID1,absID2,absID3,absID4);
+    printf("\t e %f: e1 %f, e2 %f, e3 %f, e4 %f\n",
+           ecell,ecell1,ecell2,ecell3,ecell4);
+    printf("\t t %f: t1 %f, t2 %f, t3 %f, t4 %f;\n dt1 %f, dt2 %f, dt3 %f, dt4 %f\n",
+           tcell*1.e9,tcell1*1.e9,tcell2*1.e9,tcell3*1.e9,tcell4*1.e9,
+           TMath::Abs(tcell-tcell1)*1.e9, TMath::Abs(tcell-tcell2)*1.e9, TMath::Abs(tcell-tcell3)*1.e9, TMath::Abs(tcell-tcell4)*1.e9);
+  */
+  
+  if(TMath::Abs(tcell-tcell1)*1.e9 > fExoticCellDiffTime) ecell1 = 0 ;
+  if(TMath::Abs(tcell-tcell2)*1.e9 > fExoticCellDiffTime) ecell2 = 0 ;
+  if(TMath::Abs(tcell-tcell3)*1.e9 > fExoticCellDiffTime) ecell3 = 0 ;
+  if(TMath::Abs(tcell-tcell4)*1.e9 > fExoticCellDiffTime) ecell4 = 0 ;
+
+  Float_t eCross = ecell1+ecell2+ecell3+ecell4;
+
+  //printf("\t eCell %f, eCross %f, 1-eCross/eCell %f\n",ecell,eCross,1-eCross/ecell);
+  
+  if(1-eCross/ecell > fExoticCellFraction) {
+    AliDebug(2,Form("AliEMCALRecoUtils::IsExoticCell() - EXOTIC CELL id %d, eCell %f, eCross %f, 1-eCross/eCell %f\n",absID,ecell,eCross,1-eCross/ecell));
+    return kTRUE;
+  }
+  
+  return kFALSE;
+  
+}
+
 //_________________________________________________
-Bool_t AliEMCALRecoUtils::IsExoticCluster(AliVCluster *cluster) const {
-  // Check if the cluster has high energy  but small number of cells
-  // The criteria comes from Gustavo's study
-  //
+Bool_t AliEMCALRecoUtils::IsExoticCluster(AliVCluster *cluster, AliVCaloCells *cells, const Int_t bc) {
+  // Check if the cluster highest energy tower is exotic
   
   if(!cluster){
     AliInfo("Cluster pointer null!");
     return kFALSE;
   }
   
-  Int_t nc = cluster->GetNCells() ;
+  if(!fRejectExoticCluster) return kFALSE;
   
-  if      ( nc > 8 )                   return kFALSE ; // Good cluster, needed for 3x3 clusterizer  
-  else if ( nc < 1 + cluster->E()/3. ) return kTRUE  ; // Bad cluster
-  else                                 return kFALSE ; // Good cluster
+  // Get highest energy tower
+  AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
+  Int_t iSupMod = -1, absId = -1, ieta = -1, iphi = -1;
+  Bool_t shared = kFALSE;
+  GetMaxEnergyCell(geom, cells, cluster, absId, iSupMod, ieta, iphi, shared);
+
+  return IsExoticCell(absId,cells,bc);
   
 }
 
@@ -885,7 +994,7 @@ void AliEMCALRecoUtils::RecalibrateClusterEnergy(AliEMCALGeometry* geom, AliVClu
 }
 
 //________________________________________________________________
-void AliEMCALRecoUtils::RecalibrateCells(AliEMCALGeometry* geom, AliVCaloCells * cells, Int_t bc){
+void AliEMCALRecoUtils::RecalibrateCells(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();
@@ -899,33 +1008,24 @@ void AliEMCALRecoUtils::RecalibrateCells(AliEMCALGeometry* geom, AliVCaloCells *
   
   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() ;
+  Int_t    absId  =-1;
+  Bool_t   accept = kFALSE;
+  Float_t  ecell  = 0;
+  Double_t tcell  = 0;
   
+  Int_t    nEMcell  = cells->GetNumberOfCells() ;  
   for (Int_t iCell = 0; iCell < nEMcell; iCell++) { 
     
-    absId = cells->GetCellNumber(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);
+    accept = AcceptCalibrateCell(absId, bc, ecell ,tcell ,cells); 
+    if(!accept) {
+      ecell = 0;
+      tcell = 0;
+    }
     
     //Set new values
-    cells->SetCell(iCell,cells->GetCellNumber(iCell),cellE, celltime);
+    cells->SetCell(iCell,absId,ecell, tcell);
     
   }
   
@@ -937,7 +1037,7 @@ void AliEMCALRecoUtils::RecalibrateCellTime(const Int_t absId, const Int_t bc, D
        // Recalibrate time of cell with absID  considering the recalibration map 
   // bc= bunch crossing number returned by esdevent->GetBunchCrossNumber();
   
-  if(!fCellsRecalibrated && IsTimeRecalibrationOn()){
+  if(!fCellsRecalibrated && IsTimeRecalibrationOn() && bc >= 0){
     
     celltime -= GetEMCALChannelTimeRecalibrationFactor(bc%4,absId)*1.e-9;    ;  
     
@@ -1797,17 +1897,22 @@ UInt_t AliEMCALRecoUtils::FindMatchedPosForTrack(Int_t trkIndex) const
   return pos;
 }
 
-//__________________________________________________________
-Bool_t AliEMCALRecoUtils::IsGoodCluster(AliVCluster *cluster, AliEMCALGeometry *geom, AliVCaloCells* cells)
+//___________________________________________________________________________________
+Bool_t AliEMCALRecoUtils::IsGoodCluster(AliVCluster *cluster, AliEMCALGeometry *geom, 
+                                        AliVCaloCells* cells,const Int_t bc)
 {
   // check if the cluster survives some quality cut
   //
   //
   Bool_t isGood=kTRUE;
-  if(!cluster || !cluster->IsEMCAL()) return kFALSE;
+  
+  if(!cluster || !cluster->IsEMCAL())              return kFALSE;
+  
   if(ClusterContainsBadChannel(geom,cluster->GetCellsAbsId(),cluster->GetNCells())) return kFALSE;
+  
   if(!CheckCellFiducialRegion(geom,cluster,cells)) return kFALSE;
-  if(fRejectExoticCluster && IsExoticCluster(cluster)) return kFALSE;
+  
+  if(IsExoticCluster(cluster, cells,bc))           return kFALSE;
 
   return isGood;
 }
index 6552e07..9644ee4 100644 (file)
@@ -98,7 +98,7 @@ public:
   void     SetW0(Float_t w0)                             { fW0  = w0                ; }
 
   //-----------------------------------------------------
-  //Non Linearity
+  // Non Linearity
   //-----------------------------------------------------
 
   Float_t  CorrectClusterEnergyLinearity(AliVCluster* clu) ;
@@ -127,14 +127,15 @@ public:
   Bool_t   IsClusterEnergySmeared()                const { return fSmearClusterEnergy          ; }   
   void     SetSmearingParameters(Int_t i, Float_t param) { if(i < 3){ fSmearClusterParam[i] = param ; }
                                                            else     { AliInfo(Form("Index %d larger than 2, do nothing\n",i)) ; } }
-  
   //-----------------------------------------------------
-  // Energy Recalibration
+  // Recalibration
   //-----------------------------------------------------
-  
-  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   AcceptCalibrateCell(const Int_t absId, const Int_t bc,
+                               Float_t & amp, Double_t & time, AliVCaloCells* cells) ; // Energy and Time
+  void     RecalibrateCells(AliVCaloCells * cells, Int_t bc) ; // Energy and Time
+  void     RecalibrateClusterEnergy(AliEMCALGeometry* geom, AliVCluster* cluster, AliVCaloCells * cells, const Int_t bc=-1) ; // Energy and time
 
+  // Energy recalibration
   Bool_t   IsRecalibrationOn()                     const { return fRecalibration ; }
   void     SwitchOffRecalibration()                      { fRecalibration = kFALSE ; }
   void     SwitchOnRecalibration()                       { fRecalibration = kTRUE  ; 
@@ -160,10 +161,7 @@ public:
                                                            SwitchOnRecalibration()           ; }
   void     SetRunDependentCorrections(Int_t runnumber);
       
-  //-----------------------------------------------------
-  // Time Recalibration
-  //-----------------------------------------------------
-  
+  // Time Recalibration  
   void     RecalibrateCellTime(const Int_t absId, const Int_t bc, Double_t & time);
   
   Bool_t   IsTimeRecalibrationOn()                 const { return fTimeRecalibration   ; }
@@ -172,25 +170,25 @@ public:
     if(!fEMCALTimeRecalibrationFactors)InitEMCALTimeRecalibrationFactors() ; }
   void     InitEMCALTimeRecalibrationFactors() ;
   
-  Float_t  GetEMCALChannelTimeRecalibrationFactor(Int_t bc, Int_t absID) const { 
+  Float_t  GetEMCALChannelTimeRecalibrationFactor(const Int_t bc, const Int_t absID) const { 
     if(fEMCALTimeRecalibrationFactors) 
       return (Float_t) ((TH1F*)fEMCALTimeRecalibrationFactors->At(bc))->GetBinContent(absID); 
-    else return 1 ; } 
+    else return 0 ; } 
        
-  void     SetEMCALChannelTimeRecalibrationFactor(Int_t bc,Int_t absID, Double_t c = 1) { 
+  void     SetEMCALChannelTimeRecalibrationFactor(const Int_t bc, const Int_t absID, Double_t c = 0) { 
     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)          ; }
+  TH1F *   GetEMCALChannelTimeRecalibrationFactors(const Int_t bc)const       { return (TH1F*)fEMCALTimeRecalibrationFactors->At(bc) ; }       
+  void     SetEMCALChannelTimeRecalibrationFactors(TObjArray *map)            { fEMCALTimeRecalibrationFactors = map                 ; }
+  void     SetEMCALChannelTimeRecalibrationFactors(const Int_t bc , TH1F* h)  { fEMCALTimeRecalibrationFactors->AddAt(h,bc)          ; }
   
   //-----------------------------------------------------
   // Modules fiducial region, remove clusters in borders
   //-----------------------------------------------------
 
   Bool_t   CheckCellFiducialRegion(AliEMCALGeometry* geom, AliVCluster* cluster, AliVCaloCells* cells) ;
-  void     SetNumberOfCellsFromEMCALBorder(Int_t n)      { fNCellsFromEMCALBorder = n      ; }
+  void     SetNumberOfCellsFromEMCALBorder(const Int_t n){ fNCellsFromEMCALBorder = n      ; }
   Int_t    GetNumberOfCellsFromEMCALBorder()      const  { return fNCellsFromEMCALBorder   ; }
     
   void     SwitchOnNoFiducialBorderInEMCALEta0()         { fNoEMCALBorderAtEta0 = kTRUE    ; }
@@ -285,14 +283,24 @@ public:
   void     SetStep(Double_t step)                     { fStepCluster = step           ; }
   void     SetStepSurface(Double_t step)              { fStepSurface = step           ; }
  
-  //Cluster cut
-  Bool_t   IsGoodCluster(AliVCluster *cluster, AliEMCALGeometry *geom, AliVCaloCells* cells);
-  Bool_t   IsExoticCluster(AliVCluster *cluster) const ;
-
-  void     SwitchOnRejectExoticCluster()              { fRejectExoticCluster=kTRUE     ; }
-  void     SwitchOffRejectExoticCluster()             { fRejectExoticCluster=kFALSE    ; }
+  // Exotic cells / clusters
+  
+  Bool_t   IsExoticCell(const Int_t absId, AliVCaloCells* cells, const Int_t bc =-1) ;
+  void     SwitchOnRejectExoticCell()                 { fRejectExoticCells = kTRUE     ; }
+  void     SwitchOffRejectExoticCell()                { fRejectExoticCells = kFALSE    ; } 
+   
+  void     SetExoticCellFractionCut(Float_t f)        { fExoticCellFraction     = f    ; }
+  void     SetExoticCellDiffTimeCut(Float_t dt)       { fExoticCellDiffTime     = dt   ; }
+  void     SetExoticCellMinAmplitudeCut(Float_t ma)   { fExoticCellMinAmplitude = ma   ; }
+  
+  Bool_t   IsExoticCluster(AliVCluster *cluster, AliVCaloCells* cells, const Int_t bc=0) ;
+  void     SwitchOnRejectExoticCluster()              { fRejectExoticCluster = kTRUE   ;
+                                                        fRejectExoticCells   = kTRUE   ; }
+  void     SwitchOffRejectExoticCluster()             { fRejectExoticCluster = kFALSE  ; }
   Bool_t   IsRejectExoticCluster()              const { return fRejectExoticCluster    ; }
-
+  
+  //Cluster cut
+  Bool_t   IsGoodCluster(AliVCluster *cluster, AliEMCALGeometry *geom, AliVCaloCells* cells, const Int_t bc =-1);
 
   //Track Cuts 
   Bool_t   IsAccepted(AliESDtrack *track);
@@ -368,8 +376,12 @@ private:
   Int_t      fNCellsFromEMCALBorder;     // Number of cells from EMCAL border the cell with maximum amplitude has to be.
   Bool_t     fNoEMCALBorderAtEta0;       // Do fiducial cut in EMCAL region eta = 0?
   
-  // Cluster cuts
+  // Exotic cell / cluster
   Bool_t     fRejectExoticCluster;       // Switch on or off exotic cluster rejection
+  Bool_t     fRejectExoticCells;         // Remove exotic cells
+  Float_t    fExoticCellFraction;        // Good cell if fraction < 1-ecross/ecell
+  Float_t    fExoticCellDiffTime;        // If time of candidate to exotic and close cell is too different (in ns), it must be noisy, set amp to 0
+  Float_t    fExoticCellMinAmplitude;    // Check for exotic only if amplitud is larger than this value
   
   // PID
   AliEMCALPIDUtils * fPIDUtils;          // Recalculate PID parameters
@@ -404,7 +416,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.
   
-  ClassDef(AliEMCALRecoUtils, 16)
+  ClassDef(AliEMCALRecoUtils, 17)
   
 };
 
index 1b2f2b6..7b3fa67 100644 (file)
@@ -429,7 +429,7 @@ void AliAnalysisTaskCaloFilter::UserExec(Option_t */*option*/)
       //        continue;
       //      }
 
-      if(fEMCALRecoUtils->IsRejectExoticCluster() && fEMCALRecoUtils->IsExoticCluster(cluster)) continue;      
+      if(fEMCALRecoUtils->IsExoticCluster(cluster, InputEvent()->GetEMCALCells(),InputEvent()->GetBunchCrossNumber())) continue;       
 
       fEMCALRecoUtils->GetMatchedResiduals(cluster->GetID(),dR,dZ);
       cluster->SetTrackDistance(dR,dZ);
@@ -483,7 +483,7 @@ void AliAnalysisTaskCaloFilter::UserExec(Option_t */*option*/)
       if(TMath::Abs(dR) < 990 && TMath::Abs(dZ) < 990) { //Default value in PHOS 999, in EMCAL 1024, why?
         if(DebugLevel() > 2) 
           printf("*** Cluster Track-Matched *** dR %f, dZ %f\n",caloCluster->GetEmcCpvDistance(),caloCluster->Chi2());
-        caloCluster->AddTrackMatched(0x0); 
+        caloCluster->AddTrackMatched(new AliAODTrack); 
       }// fill the array with one entry to signal a possible match
       //TArrayI* matchedT =    ((AliESDCaloCluster*)cluster)->GetTracksMatched();
       //if (InputEvent()->GetNumberOfTracks() > 0 && matchedT && cluster->GetTrackMatchedIndex() >= 0) {       
index 55fec4f..665f4b5 100644 (file)
@@ -32,7 +32,6 @@
 #include "TROOT.h"
 #include "TInterpreter.h"
 #include "TFile.h"
-//#include "string.h"
 
 // --- AliRoot Analysis Steering
 #include "AliAnalysisTask.h"
@@ -83,8 +82,7 @@ AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize(const char *name)
 , fCellLabels(),          fCellSecondLabels(),        fCellTime()
 , fMaxEvent(1000000000),  fDoTrackMatching(kFALSE)
 , fSelectCell(kFALSE),    fSelectCellMinE(0.005),     fSelectCellMinFrac(0.001)
-, fRemoveLEDEvents(kFALSE), fRemoveExoticEvents(kFALSE), fRemoveExoticCells(kFALSE)
-, fExoticCellFraction(0.9), fExoticCellDiffTime(1e6),    fExoticCellMinAmplitude(0.5)
+, fRemoveLEDEvents(kFALSE), fRemoveExoticEvents(kFALSE)
 
 {
   //ctor
@@ -122,8 +120,7 @@ AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize()
 , fCellLabels(),            fCellSecondLabels(),        fCellTime()
 , fMaxEvent(1000000000),    fDoTrackMatching(kFALSE)
 , fSelectCell(kFALSE),      fSelectCellMinE(0.005),     fSelectCellMinFrac(0.001)
-, fRemoveLEDEvents(kFALSE), fRemoveExoticEvents(kFALSE), fRemoveExoticCells(kFALSE)
-, fExoticCellFraction(0.9), fExoticCellDiffTime(1e6),    fExoticCellMinAmplitude(0.5)
+, fRemoveLEDEvents(kFALSE), fRemoveExoticEvents(kFALSE)
 {
   // Constructor
   for(Int_t i = 0; i < 10;    i++)  fGeomMatrix[i] =  0;
@@ -167,52 +164,6 @@ AliAnalysisTaskEMCALClusterize::~AliAnalysisTaskEMCALClusterize()
   
 }
 
-//_________________________________________________________________________________
-Bool_t AliAnalysisTaskEMCALClusterize::AcceptCalibrateCell(const Int_t absID,
-                                                           Float_t  & amp, 
-                                                           Double_t & time, 
-                                                           AliVCaloCells* cells) 
-{
-  // Reject cell if criteria not passed and calibrate it
-  
-  if(absID < 0 || absID >= 24*48*fGeom->GetNumberOfSuperModules()) return kFALSE;
-  
-  Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1; 
-  fGeom->GetCellIndex(absID,imod,iTower,iIphi,iIeta); 
-  fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);      
-  
-  // Do not include bad channels found in analysis,
-  if( fRecoUtils->IsBadChannelsRemovalSwitchedOn() && 
-     fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi)) {
-    return kFALSE;
-  }
-  //Recalibrate energy
-  amp  = cells->GetCellAmplitude(absID);
-  if(fRecoUtils->IsRecalibrationOn()){ 
-    amp *=fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
-  }
-  
-  // Do not include cells with too low energy, nor exotic cell
-  if(amp < fRecParam->GetMinECut() ) {
-    amp = 0.;
-    return kFALSE;
-  }
-  
-  // Recalibrate time
-  time = cells->GetCellTime(absID);
-  Int_t bc = InputEvent()->GetBunchCrossNumber();
-  
-  // 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[absID];
-  
-  fRecoUtils->RecalibrateCellTime(absID,bc,time);
-    
-  //printf("AliAnalysisTaskEMCALClusterize::AcceptCalibrateCell() - Accepted Cell id %d, amp %f, t %f\n",absID,amp,time*1.e9);
-  
-  return kTRUE;
-}
-
-
 //_________________________________________________
 Bool_t AliAnalysisTaskEMCALClusterize::AccessOCDB()
 {
@@ -396,14 +347,27 @@ void AliAnalysisTaskEMCALClusterize::ClusterizeCells()
     
   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++)
   {
     
     // Get cell values, recalibrate and not include bad channels found in analysis, nor cells with too low energy, nor exotic cell
     id = cells->GetCellNumber(icell);
-    Bool_t accept =  AcceptCalibrateCell(id,amp,time,cells);
-    if(  accept && IsExoticCell(id,amp,time,cells)) accept = kFALSE;
+    Bool_t accept = fRecoUtils->AcceptCalibrateCell(id,bc,amp,time,cells);
+    
+    // Do not include cells with too low energy, nor exotic cell
+    if(amp < fRecParam->GetMinECut() ) accept = kFALSE;
+    
+    // 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];
+      fRecoUtils->RecalibrateCellTime(id,bc,time);
+    }
+    
+    if(  accept && fRecoUtils->IsExoticCell(id,cells,bc)){
+      accept = kFALSE;
+    }
+
     if( !accept ){
       fCellLabels[id]      =-1; //reset the entry in the array for next event
       fCellSecondLabels[id]=-1; //reset the entry in the array for next event
@@ -809,71 +773,6 @@ void AliAnalysisTaskEMCALClusterize::InitGeometry()
   
 }
 
-//___________________________________________________________________
-Bool_t AliAnalysisTaskEMCALClusterize::IsExoticCell(const Int_t absID, 
-                                                    const Float_t ecell, 
-                                                    const Float_t tcell, 
-                                                    AliVCaloCells* cells)
-{
-  //Look to cell neighbourhood and reject if it seems exotic
-  
-  if(!fRemoveExoticCells) return kFALSE;
-  
-  if(ecell < fExoticCellMinAmplitude) return kFALSE;
-  
-  Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1; 
-  fGeom->GetCellIndex(absID,imod,iTower,iIphi,iIeta); 
-  fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);      
-  
-  //Get close cells index, energy and time, not in corners
-  Int_t absID1 = fGeom-> GetAbsCellIdFromCellIndexes(imod, iphi+1, ieta);
-  Int_t absID2 = fGeom-> GetAbsCellIdFromCellIndexes(imod, iphi-1, ieta);
-  Int_t absID3 = fGeom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
-  Int_t absID4 = fGeom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1);
-  
-  Float_t  ecell1  = 0, ecell2  = 0, ecell3  = 0, ecell4  = 0;
-  Double_t tcell1  = 0, tcell2  = 0, tcell3  = 0, tcell4  = 0;
-  Bool_t   accept1 = 0, accept2 = 0, accept3 = 0, accept4 = 0;
-  
-  accept1 = AcceptCalibrateCell(absID1,ecell1,tcell1,cells); 
-  accept2 = AcceptCalibrateCell(absID2,ecell2,tcell2,cells); 
-  accept3 = AcceptCalibrateCell(absID3,ecell3,tcell3,cells); 
-  accept4 = AcceptCalibrateCell(absID4,ecell4,tcell4,cells); 
-
-  if(DebugLevel() > 2 ){
-    printf("AliAnalysisTaksEMCALClusterize::IsExoticCell() -  Cell absID %d \n",absID);
-    printf("\t  accept1 %d, accept2 %d, accept3 %d, accept4 %d\n",
-         accept1,accept2,accept3,accept4);
-    printf("\t id %d: id1 %d, id2 %d, id3 %d, id4 %d\n",
-           absID,absID1,absID2,absID3,absID4);
-    printf("\t e %f: e1 %f, e2 %f, e3 %f, e4 %f\n",
-           ecell,ecell1,ecell2,ecell3,ecell4);
-    printf("\t t %f: t1 %f, t2 %f, t3 %f, t4 %f;\n dt1 %f, dt2 %f, dt3 %f, dt4 %f\n",
-           tcell*1.e9,tcell1*1.e9,tcell2*1.e9,tcell3*1.e9,tcell4*1.e9,
-           TMath::Abs(tcell-tcell1)*1.e9, TMath::Abs(tcell-tcell2)*1.e9, TMath::Abs(tcell-tcell3)*1.e9, TMath::Abs(tcell-tcell4)*1.e9);
-  }
-  
-  if(TMath::Abs(tcell-tcell1)*1.e9 > fExoticCellDiffTime) ecell1 = 0 ;
-  if(TMath::Abs(tcell-tcell2)*1.e9 > fExoticCellDiffTime) ecell2 = 0 ;
-  if(TMath::Abs(tcell-tcell3)*1.e9 > fExoticCellDiffTime) ecell3 = 0 ;
-  if(TMath::Abs(tcell-tcell4)*1.e9 > fExoticCellDiffTime) ecell4 = 0 ;
-  
-  
-  Float_t eCross = ecell1+ecell2+ecell3+ecell4;
-  if(DebugLevel() > 2 ){  
-    printf("\t eCell %f, eCross %f, 1-eCross/eCell %f\n",ecell,eCross,1-eCross/ecell);
-  }
-  
-  if(1-eCross/ecell > fExoticCellFraction) {
-    //cells->SetCell(icell,absID,0,1e6);
-    if(DebugLevel() > 2 ) printf("AliAnalysisTaskEMCALClusterize::IsExoticCell() - EXOTIC CELL REMOVED id %d, eCell %f, eCross %f, 1-eCross/eCell %f\n",absID,ecell,eCross,1-eCross/ecell);
-    return kTRUE;
-  }
-  
-  return kFALSE;
-  
-}
-
 //____________________________________________________
 Bool_t AliAnalysisTaskEMCALClusterize::IsExoticEvent()
 {
@@ -886,14 +785,15 @@ Bool_t AliAnalysisTaskEMCALClusterize::IsExoticEvent()
   // Loop on cells
   AliVCaloCells * cells = fEvent->GetEMCALCells();
   Float_t totCellE = 0;
+  Int_t bc = InputEvent()->GetBunchCrossNumber();
   for(Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++){
     
     Float_t  ecell = 0 ;
     Double_t tcell = 0 ;
     
     Int_t absID   = cells->GetCellNumber(icell);
-    Bool_t accept = AcceptCalibrateCell(absID,ecell,tcell,cells);
-    if(accept && !IsExoticCell(absID, ecell,tcell,cells)) totCellE += ecell;
+    Bool_t accept = fRecoUtils->AcceptCalibrateCell(absID,bc,ecell,tcell,cells);
+    if(accept && !fRecoUtils->IsExoticCell(absID,cells,bc)) totCellE += ecell;
   }
   
   //  TString triggerclasses = "";
index d2b5f6e..7c65f97 100644 (file)
@@ -44,20 +44,9 @@ class AliAnalysisTaskEMCALClusterize : public AliAnalysisTaskSE {
   void           SwitchOnExoticEventsRemoval()                  { fRemoveExoticEvents= kTRUE   ; }
   void           SwitchOffExoticEventsRemoval()                 { fRemoveExoticEvents= kFALSE  ; } 
   
-  Bool_t         IsExoticCell(const Int_t absId, const Float_t ecell, 
-                              const Float_t tcell, AliVCaloCells* cells);
-  void           SwitchOnExoticCellRemoval()                    { fRemoveExoticCells = kTRUE   ; }
-  void           SwitchOffExoticCellRemoval()                   { fRemoveExoticCells = kFALSE  ; } 
-  
-  void           SetExoticCellFractionCut(Float_t f)            { fExoticCellFraction = f      ; }
-  void           SetExoticCellDiffTimeCut(Float_t dt)           { fExoticCellDiffTime = dt     ; }
-  void           SetExoticCellMinAmplitudeCut(Float_t ma)       { fExoticCellMinAmplitude = ma ; }
-  
   Bool_t         IsLEDEvent();
   void           SwitchOnLEDEventsRemoval()                     { fRemoveLEDEvents   = kTRUE   ; }
   void           SwitchOffLEDEventsRemoval()                    { fRemoveLEDEvents   = kFALSE  ; } 
-
-  Bool_t         AcceptCalibrateCell(const Int_t absId, Float_t & amp, Double_t & time, AliVCaloCells* cells) ;
   
   //OCDB
   Bool_t         AccessOCDB();
@@ -160,12 +149,8 @@ class AliAnalysisTaskEMCALClusterize : public AliAnalysisTaskSE {
   Float_t                fSelectCellMinFrac;       // Min fraction of cell energy after unfolding cut
   Bool_t                 fRemoveLEDEvents;         // Remove LED events, use only for LHC11a 
   Bool_t                 fRemoveExoticEvents;      // Remove exotic events
-  Bool_t                 fRemoveExoticCells;       // Remove exotic cells
-  Float_t                fExoticCellFraction;      // Good cell if fraction < 1-ecross/ecell
-  Float_t                fExoticCellDiffTime;      // If time of candidate to exotic and close cell is too different, it must be noisy, set amp to 0
-  Float_t                fExoticCellMinAmplitude;  // Check for exotic only if amplitud is larger than this value
   
-  ClassDef(AliAnalysisTaskEMCALClusterize, 14);
+  ClassDef(AliAnalysisTaskEMCALClusterize, 15);
 
 };
 
index 444d118..f7a1622 100644 (file)
@@ -47,15 +47,18 @@ ClassImp(AliAnalysisTaskEMCALPi0CalibSelection)
 AliAnalysisTaskEMCALPi0CalibSelection::AliAnalysisTaskEMCALPi0CalibSelection(const char* name) :
 AliAnalysisTaskSE(name),fEMCALGeo(0x0), 
 fEmin(0.5),               fEmax(15.),      
-fL0min(0.01),             fL0max(0.5),     fDTimeCut(100.), 
-fAsyCut(1.),              fMinNCells(2),   fGroupNCells(0),
-fLogWeight(4.5),          fSameSM(kFALSE), fFilteredInput(kFALSE),
+fL0min(0.01),             fL0max(0.5),              fDTimeCut(100.), 
+fAsyCut(1.),              fMinNCells(2),            fGroupNCells(0),
+fLogWeight(4.5),          fSameSM(kFALSE),          fFilteredInput(kFALSE),
 fCorrectClusters(kFALSE), fEMCALGeoName("EMCAL_COMPLETEV1"), 
 fTriggerName("EMC"),      fRecoUtils(new AliEMCALRecoUtils), 
 fCuts(0x0),               fLoadMatrices(0),
 fNMaskCellColumns(11),    fMaskCellColumns(0x0),
+fInvMassCutMin(110.),     fInvMassCutMax(160.),
 //Histograms
-fNbins(300), fMinBin(0.), fMaxBin(300.),   fOutputContainer(0x0),
+fOutputContainer(0x0),    fNbins(300),              
+fMinBin(0.),              fMaxBin(300.),   
+fNTimeBins(1000),         fMinTimeBin(0.),          fMaxTimeBin(1000.),   
 fHmgg(0x0),               fHmggDifferentSM(0x0), 
 fHmggMaskFrame(0x0),      fHmggDifferentSMMaskFrame(0x0), 
 fHOpeningAngle(0x0),      fHOpeningAngleDifferentSM(0x0),  
@@ -69,11 +72,17 @@ fhClusterTime(0x0),       fhClusterPairDiffTime(0x0)
   for(Int_t iMod=0; iMod < AliEMCALGeoParams::fgkEMCALModules; iMod++) {
     for(Int_t iX=0; iX<24; iX++) {
       for(Int_t iZ=0; iZ<48; iZ++) {
-        fHmpi0[iMod][iZ][iX]=0;
+        fHmpi0[iMod][iZ][iX]   = 0 ;
       }
     } 
   }
   
+  fHTpi0[0]= 0 ;
+  fHTpi0[1]= 0 ;
+  fHTpi0[2]= 0 ;
+  fHTpi0[3]= 0 ;
+  
+  
   fMaskCellColumns = new Int_t[fNMaskCellColumns];
   fMaskCellColumns[0] = 6 ;  fMaskCellColumns[1] = 7 ;  fMaskCellColumns[2] = 8 ; 
   fMaskCellColumns[3] = 35;  fMaskCellColumns[4] = 36;  fMaskCellColumns[5] = 37; 
@@ -141,14 +150,16 @@ void AliAnalysisTaskEMCALPi0CalibSelection::LocalInit()
   char onePar[buffersize] ;
   fCuts = new TList();
   
-  snprintf(onePar,buffersize, "Custer cuts: %2.2f < E < %2.2f GeV; %2.2f < Lambda0_2 < %2.2f GeV; min number of cells %d; Assymetry cut %1.2f, time1-time2 < %2.2f", 
-           fEmin,fEmax, fL0min, fL0max, fMinNCells, fAsyCut, fDTimeCut) ;
+  snprintf(onePar,buffersize, "Custer cuts: %2.2f < E < %2.2f GeV; %2.2f < Lambda0_2 < %2.2f GeV; min number of cells %d; Assymetry cut %1.2f, time1-time2 < %2.2f; %3.1f < Mass < %3.1f", 
+           fEmin,fEmax, fL0min, fL0max, fMinNCells, fAsyCut, fDTimeCut, fInvMassCutMin, fInvMassCutMax) ;
   fCuts->Add(new TObjString(onePar));
   snprintf(onePar,buffersize, "Group %d cells;", fGroupNCells) ;
   fCuts->Add(new TObjString(onePar));
   snprintf(onePar,buffersize, "Cluster maximal cell away from border at least %d cells;", fRecoUtils->GetNumberOfCellsFromEMCALBorder()) ;
   fCuts->Add(new TObjString(onePar));
-  snprintf(onePar,buffersize, "Histograms: bins %d; energy range: %2.2f < E < %2.2f GeV;",fNbins,fMinBin,fMaxBin) ;
+  snprintf(onePar,buffersize, "Histograms, Mass bins %d; energy range: %2.2f < E < %2.2f GeV;",fNbins,fMinBin,fMaxBin) ;
+  fCuts->Add(new TObjString(onePar));
+  snprintf(onePar,buffersize, "Histograms, Time bins %d; energy range: %2.2f < E < %2.2f GeV;",fNTimeBins,fMinTimeBin,fMaxTimeBin) ;
   fCuts->Add(new TObjString(onePar));
   snprintf(onePar,buffersize, "Switchs: Remove Bad Channels? %d; Use filtered input? %d;  Correct Clusters? %d, Mass per channel same SM clusters? %d ",
            fRecoUtils->IsBadChannelsRemovalSwitchedOn(),fFilteredInput,fCorrectClusters, fSameSM) ;
@@ -181,11 +192,21 @@ void AliAnalysisTaskEMCALPi0CalibSelection::UserCreateOutputObjects()
         snprintf(hname,buffersize, "%d_%d_%d",iMod,iCol,iRow);
         snprintf(htitl,buffersize, "Two-gamma inv. mass for super mod %d, cell(col,row)=(%d,%d)",iMod,iCol,iRow);
         fHmpi0[iMod][iCol][iRow] = new TH1F(hname,htitl,fNbins,fMinBin,fMaxBin);
+        fHmpi0[iMod][iCol][iRow]->SetXTitle("mass (MeV/c^{2})");
         fOutputContainer->Add(fHmpi0[iMod][iCol][iRow]);
       }
     }
   }
   
+  Int_t nchannels = nSM*AliEMCALGeoParams::fgkEMCALRows*AliEMCALGeoParams::fgkEMCALCols;
+  for(Int_t ibc = 0; ibc < 4; ibc++){
+    fHTpi0[ibc] = new TH2F(Form("hTime_BC%d",ibc),Form("Time of cell clusters under pi0 peak, bunch crossing %d",ibc),
+                           nchannels,0,nchannels, fNTimeBins,fMinTimeBin,fMaxTimeBin);
+    fOutputContainer->Add(fHTpi0[ibc]);       
+    fHTpi0[ibc]->SetYTitle("time (ns)");
+    fHTpi0[ibc]->SetXTitle("abs. Id. ");
+  }
+  
   fHmgg = new TH2F("hmgg","2-cluster invariant mass",fNbins,fMinBin,fMaxBin,100,0,10);
   fHmgg->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
   fHmgg->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
@@ -398,7 +419,7 @@ void AliAnalysisTaskEMCALPi0CalibSelection::UserCreateOutputObjects()
   fhClusterTime->SetYTitle("t (ns)");
   fOutputContainer->Add(fhClusterTime);
   
-  fhClusterPairDiffTime = new TH2F("hClusterPairDiffTime","cluster pair time difference vs E",100,0,10, 200,-100,100);
+  fhClusterPairDiffTime = new TH2F("hClusterPairDiffTime","cluster pair time difference vs E",100,0,10, 800,-400,400);
   fhClusterPairDiffTime->SetXTitle("E_{pair} (GeV)");
   fhClusterPairDiffTime->SetYTitle("#Delta t (ns)");
   fOutputContainer->Add(fhClusterPairDiffTime);
@@ -587,7 +608,8 @@ void AliAnalysisTaskEMCALPi0CalibSelection::UserExec(Option_t* /* option */)
   }
   
   //----------------------------------------------------------
-  //Now the invariant mass analysis with the corrected clusters  
+  //Now the invariant mass analysis with the corrected clusters
+  Int_t bc = event->GetBunchCrossNumber();
   for(Int_t iClu=0; iClu<kNumberOfEMCALClusters-1; iClu++) {
     
     AliVCluster *c1 = (AliVCluster *) caloClustersArr->At(iClu);
@@ -596,6 +618,7 @@ void AliAnalysisTaskEMCALPi0CalibSelection::UserExec(Option_t* /* option */)
     Float_t e1i = c1->E();   // cluster energy before correction   
     if      (e1i < fEmin) continue;
     else if (e1i > fEmax) continue;
+    else if (!fRecoUtils->IsGoodCluster(c1,fEMCALGeo,emCells,bc));
     else if (c1->GetNCells() < fMinNCells) continue; 
     else if (c1->GetM02() < fL0min || c1->GetM02() > fL0max) continue;
 
@@ -623,11 +646,11 @@ void AliAnalysisTaskEMCALPi0CalibSelection::UserExec(Option_t* /* option */)
     for (Int_t jClu=iClu+1; jClu<kNumberOfEMCALClusters; jClu++) {
       AliAODCaloCluster *c2 = (AliAODCaloCluster *) caloClustersArr->At(jClu);
       
-      if(fRecoUtils->ClusterContainsBadChannel(fEMCALGeo, c2->GetCellsAbsId(), c2->GetNCells())) continue;     
-      
+
       Float_t e2i = c2->E();
       if      (e2i < fEmin) continue;
       else if (e2i > fEmax) continue;
+      else if (!fRecoUtils->IsGoodCluster(c2,fEMCALGeo,emCells,bc))continue;
       else if (c2->GetNCells() < fMinNCells) continue; 
       else if (c2->GetM02() < fL0min || c2->GetM02() > fL0max) continue;
 
@@ -710,9 +733,21 @@ void AliAnalysisTaskEMCALPi0CalibSelection::UserExec(Option_t* /* option */)
           }// Pair not facing frame
           
           
-          if(invmass > 100. && invmass < 160.){//restrict to clusters really close to pi0 peak
+          if(invmass > fInvMassCutMin && invmass < fInvMassCutMax){//restrict to clusters really close to pi0 peak
+            
+            
+            // Check time of cells in both clusters, and fill time histogram
+            for(Int_t icell = 0; icell < c1->GetNCells(); icell++){
+              Int_t absID = c1->GetCellAbsId(icell);   
+              fHTpi0[bc%4]->Fill(absID, emCells->GetCellTime(absID)*1.e9);  
+            }
+            
+            for(Int_t icell = 0; icell < c2->GetNCells(); icell++){
+              Int_t absID = c2->GetCellAbsId(icell);   
+              fHTpi0[bc%4]->Fill(absID, emCells->GetCellTime(absID)*1.e9);  
+            }
             
-            //Opening angle of 2 photons
+             //Opening angle of 2 photons
             Float_t opangle = p1.Angle(p2.Vect())*TMath::RadToDeg();
             //printf("*******>>>>>>>> In PEAK pt %f, angle %f \n",p12.Pt(),opangle);
             
index 31fd9de..61de935 100644 (file)
@@ -57,6 +57,9 @@ public:
   void    SetNCellsGroup(Int_t n)                        { fGroupNCells = n          ; }
   void    SetLogWeight(Float_t w)                        { fLogWeight   = w          ; }
          
+  void    SetPairMinMassCut(Float_t min)                 { fInvMassCutMin = min      ; }
+  void    SetPairMaxMassCut(Float_t max)                 { fInvMassCutMax = max      ; }
+  
   void    SwitchOnSameSM()                               { fSameSM = kTRUE           ; }
   void    SwitchOffSameSM()                              { fSameSM = kFALSE          ; }
   
@@ -80,8 +83,12 @@ public:
   AliEMCALRecoUtils* GetEMCALRecoUtils()    const        { return fRecoUtils         ; }
   
   void    SetInvariantMassHistoBinRange(Int_t nBins, Float_t minbin, Float_t maxbin){
-                                        fNbins = nBins; fMinBin = minbin; fMaxBin = maxbin; }
+                                        fNbins     = nBins ; fMinBin     = minbin ; fMaxBin     = maxbin ; }
+
+  void    SetTimeHistoBinRange         (Int_t nBins, Float_t minbin, Float_t maxbin){
+                                        fNTimeBins = nBins ; fMinTimeBin = minbin ; fMaxTimeBin = maxbin ; }
 
+  
   // Mask clusters
   void    SetNMaskCellColumns(Int_t n) {
     if(n > fNMaskCellColumns){ delete [] fMaskCellColumns ; fMaskCellColumns = new Int_t[n] ; }
@@ -120,12 +127,22 @@ private:
   Int_t   fNMaskCellColumns;   // Number of masked columns
   Int_t*  fMaskCellColumns;    //[fNMaskCellColumns] list of masked cell collumn
   
+  // Pi0 clusters selection
+  
+  Float_t fInvMassCutMin;      // Min mass cut for clusters to fill time or other histograms
+  Float_t fInvMassCutMax;      // Mas mass cut for clusters to fill time or other histograms
+  
   //Output histograms  
+
+  TList*  fOutputContainer;    //!histogram container
+
   Int_t   fNbins;              // N       mass bins of invariant mass histograms
   Float_t fMinBin;             // Minimum mass bins of invariant mass histograms
   Float_t fMaxBin;             // Maximum mass bins of invariant mass histograms
-
-  TList*  fOutputContainer;    //!histogram container
+  
+  Int_t   fNTimeBins;          // N       time bins of invariant mass histograms
+  Float_t fMinTimeBin;         // Minimum time bins of invariant mass histograms
+  Float_t fMaxTimeBin;         // Maximum time bins of invariant mass histograms
   
   TH1F*   fHmpi0[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows];//! two-cluster inv. mass assigned to each cell.
 
@@ -164,6 +181,7 @@ private:
   TH1I*   fhNEvents;                                                             //! Number of events counter histogram
  
   //Time
+  TH2F*   fHTpi0[4];                                                             //! Time of cell under pi0 mass, for 4 bunch crossings
   TH2F*   fhClusterTime ;                                                        //! Timing of clusters vs energy
   TH2F*   fhClusterTimeSM[AliEMCALGeoParams::fgkEMCALModules] ;                  //! Timing of clusters vs energy per SM
   TH2F*   fhClusterPairDiffTime;                                                 //! Diference in time of clusters
@@ -171,7 +189,7 @@ private:
   TH2F*   fhClusterPairDiffTimeSameSector[AliEMCALGeoParams::fgkEMCALModules/2]; //! Diference in time of clusters same sector
   TH2F*   fhClusterPairDiffTimeSameSide[AliEMCALGeoParams::fgkEMCALModules-2];   //! Diference in time of clusters same side
 
-  ClassDef(AliAnalysisTaskEMCALPi0CalibSelection,16);
+  ClassDef(AliAnalysisTaskEMCALPi0CalibSelection,17);
 
 };
 
index b7038cd..441a923 100755 (executable)
@@ -735,7 +735,7 @@ 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;
+  if(!GetCaloUtils()->GetEMCALRecoUtils()->IsGoodCluster(clus,GetCaloUtils()->GetEMCALGeometry(),GetEMCALCells(),fInputEvent->GetBunchCrossNumber())) return;
   //  //Check if the cluster contains any bad channel and if close to calorimeter borders
   //  if(GetCaloUtils()->ClusterContainsBadChannel("EMCAL",clus->GetCellsAbsId(), clus->GetNCells())) 
   //    return;
index 8822f4e..f334588 100755 (executable)
@@ -177,7 +177,8 @@ public:
   Int_t            GetV0Signal(Int_t i)              const { return fV0ADC[i]              ; }
   Int_t            GetV0Multiplicity(Int_t i)        const { return fV0Mul[i]              ; }
   
-  void             SetEMCALClusterListName(TString &name)  {fEMCALClustersListName = name  ; }
+  void             SetEMCALClusterListName(TString &name)  { fEMCALClustersListName = name ; }
+  TString          GetEMCALClusterListName()               { return fEMCALClustersListName ; }
 
   // Arrayes with clusters/track/cells access method
   virtual TObjArray*     GetCTSTracks()              const { return fCTSTracks             ; }