Move some methods from EMCALcalibration task and CalorimterUtils class to AliEMCALRec...
authorgconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 24 Oct 2010 19:14:33 +0000 (19:14 +0000)
committergconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 24 Oct 2010 19:14:33 +0000 (19:14 +0000)
Calibration task: remove AOD filtering, reorganize class first correct clusters, if requested, and then do invariant mass.
Filter task, updated with EMCAL reconstruction utilities access and other stuff.

PWG4/CaloCalib/AliAnalysisTaskCaloFilter.cxx
PWG4/CaloCalib/AliAnalysisTaskCaloFilter.h
PWG4/CaloCalib/AliAnalysisTaskEMCALPi0CalibSelection.cxx
PWG4/CaloCalib/AliAnalysisTaskEMCALPi0CalibSelection.h
PWG4/CaloCalib/macros/anaEMCALCalib.C
PWG4/PartCorrBase/AliCalorimeterUtils.cxx
PWG4/PartCorrBase/AliCalorimeterUtils.h

index 82f8c05..b228f3a 100644 (file)
 #include "AliAnalysisTaskCaloFilter.h"
 #include "AliESDEvent.h"
 #include "AliAODEvent.h"
-#include "AliESDCaloCluster.h"
-#include "AliESDCaloCells.h"
 #include "AliLog.h"
+#include "AliVCluster.h"
+#include "AliVCaloCells.h"
+#include "AliEMCALRecoUtils.h"
+#include "AliEMCALGeometry.h"
+#include "AliVEventHandler.h"
+#include "AliAnalysisManager.h"
+#include "AliInputEventHandler.h"
 
 ClassImp(AliAnalysisTaskCaloFilter)
   
 ////////////////////////////////////////////////////////////////////////
 
 AliAnalysisTaskCaloFilter::AliAnalysisTaskCaloFilter():
-  AliAnalysisTaskSE(),
-  fCalorimeter("EMCAL PHOS")
+  AliAnalysisTaskSE(), //fCuts(0x0),
+  fCaloFilter(0), fCorrect(kFALSE), 
+  fEMCALGeo(0x0),fEMCALGeoName("EMCAL_FIRSTYEAR"), 
+  fEMCALRecoUtils(new AliEMCALRecoUtils)
 {
   // Default constructor
 }
 
 //__________________________________________________
 AliAnalysisTaskCaloFilter::AliAnalysisTaskCaloFilter(const char* name):
-  AliAnalysisTaskSE(name),
-  fCalorimeter("EMCAL PHOS")
+  AliAnalysisTaskSE(name), //fCuts(0x0),
+  fCaloFilter(0), fCorrect(kFALSE),
+  fEMCALGeo(0x0),fEMCALGeoName("EMCAL_FIRSTYEAR"), 
+  fEMCALRecoUtils(new AliEMCALRecoUtils)
 {
   // Constructor
 }
 
 //__________________________________________________
-void AliAnalysisTaskCaloFilter::CreateAODFromAOD()
+AliAnalysisTaskCaloFilter::~AliAnalysisTaskCaloFilter()
 {
-
-  // Copy AOD header, vertex, CaloClusters and CaloCells to output AOD
-  
-  AliAODEvent* aod = dynamic_cast<AliAODEvent*>(InputEvent());
-  if(!aod) {
-    printf("AliAnalysisTaskCaloFilter::CreateAODFromAOD() - This event does not contain AODs?");
-    return;
-  }
-  
-  // set arrays and pointers
-  Float_t posF[3];
-  Double_t pos[3];
-  
-  Double_t covVtx[6];
-  
-  for (Int_t i = 0; i < 6; i++)  covVtx[i] = 0.;
-  
-  // Update the header
-  AliAODHeader* headerin = aod->GetHeader();
-  AliAODHeader* header = AODEvent()->GetHeader();
-  header->SetRunNumber(headerin->GetRunNumber());
-  header->SetBunchCrossNumber(headerin->GetBunchCrossNumber());
-  header->SetOrbitNumber(headerin->GetOrbitNumber());
-  header->SetPeriodNumber(headerin->GetPeriodNumber());
-  header->SetEventType(headerin->GetEventType());
-  header->SetMuonMagFieldScale(headerin->GetMuonMagFieldScale()); // FIXME
-  header->SetCentrality(headerin->GetCentrality());        // FIXME
+  //Destructor.
+       
+  if(fEMCALGeo)       delete fEMCALGeo;        
+  if(fEMCALRecoUtils) delete fEMCALRecoUtils ;
   
+}
+
+//__________________________________________________
+void AliAnalysisTaskCaloFilter::UserCreateOutputObjects()
+{
+  // Init geometry 
+       
+  fEMCALGeo =  AliEMCALGeometry::GetInstance(fEMCALGeoName) ;  
   
-  header->SetTriggerMask(headerin->GetTriggerMask()); 
-  header->SetTriggerCluster(headerin->GetTriggerCluster());
-  header->SetMagneticField(headerin->GetMagneticField());
-  header->SetZDCN1Energy(headerin->GetZDCN1Energy());
-  header->SetZDCP1Energy(headerin->GetZDCP1Energy());
-  header->SetZDCN2Energy(headerin->GetZDCN2Energy());
-  header->SetZDCP2Energy(headerin->GetZDCP2Energy());
-  header->SetZDCEMEnergy(headerin->GetZDCEMEnergy(0),headerin->GetZDCEMEnergy(1));
-  Float_t diamxy[2]={aod->GetDiamondX(),aod->GetDiamondY()};
-  Float_t diamcov[3]; aod->GetDiamondCovXY(diamcov);
-  header->SetDiamond(diamxy,diamcov);
-  //
+}  
+
+//__________________________________________________
+void AliAnalysisTaskCaloFilter::UserExec(Option_t */*option*/)
+{
+  // Execute analysis for current event
   //
-  Int_t nVertices = 1 ;/* = prim. vtx*/;
-  Int_t nCaloClus = aod->GetNumberOfCaloClusters();
-  
-  AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
-  
-  // Access to the AOD container of vertices
-  TClonesArray &vertices = *(AODEvent()->GetVertices());
-  Int_t jVertices=0;
   
-  // Add primary vertex. The primary tracks will be defined
-  // after the loops on the composite objects (V0, cascades, kinks)
-  const AliAODVertex *vtx = aod->GetPrimaryVertex();
+  //Magic line to write events to file
+  AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kTRUE);
   
-  vtx->GetXYZ(pos); // position
-  vtx->GetCovMatrix(covVtx); //covariance matrix
+  if (fDebug > 0)  
+    printf("CaloFilter: Analysing event # %d\n", (Int_t)Entry());
   
-  AliAODVertex * primary = new(vertices[jVertices++])
-    AliAODVertex(pos, covVtx, vtx->GetChi2perNDF(), NULL, -1, AliAODVertex::kPrimary);
-  primary->SetName(vtx->GetName());
-  primary->SetTitle(vtx->GetTitle());
+  // Copy input ESD or AOD header, vertex, CaloClusters and CaloCells to output AOD
   
-  // Access to the AOD container of clusters
-  TClonesArray &caloClusters = *(AODEvent()->GetCaloClusters());
-  Int_t jClusters=0;
+  AliVEvent* event = InputEvent();
+  Bool_t bAOD = kFALSE;
+  if(!strcmp(event->GetName(),"AliAODEvent")) bAOD=kTRUE;
+  Bool_t bESD = kFALSE;
+  if(!strcmp(event->GetName(),"AliESDEvent")) bESD=kTRUE;
   
-  for (Int_t iClust=0; iClust<nCaloClus; ++iClust) {
-    
-    AliAODCaloCluster * cluster = aod->GetCaloCluster(iClust);
-    
-    //Check which calorimeter information we want to keep.
-    if     (fCalorimeter.Contains("PHOS")  && !fCalorimeter.Contains("EMCAL") && cluster->IsEMCAL()) continue ;
-    else if(fCalorimeter.Contains("EMCAL") && !fCalorimeter.Contains("PHOS")  && cluster->IsPHOS())  continue ;
-    
-    Int_t id       = cluster->GetID();
-    Float_t energy = cluster->E();
-    cluster->GetPosition(posF);
-    Char_t ttype   = cluster->GetType(); 
-    
-    
-    AliAODCaloCluster *caloCluster = new(caloClusters[jClusters++]) 
-      AliAODCaloCluster(id,
-                       0,
-                       0x0,
-                       energy,
-                       posF,
-                       NULL,
-                       ttype);
-    
-    caloCluster->SetCaloCluster(cluster->GetDistanceToBadChannel(),
-                               cluster->GetDispersion(),
-                               cluster->GetM20(), cluster->GetM02(),
-                               cluster->GetEmcCpvDistance(),  
-                               cluster->GetNExMax(),cluster->GetTOF()) ;
-    
-    caloCluster->SetPIDFromESD(cluster->GetPID());
-    caloCluster->SetNCells(cluster->GetNCells());
-    caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
-    caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
-    
-  } 
-  caloClusters.Expand(jClusters); // resize TObjArray to 'remove' slots for pseudo clusters     
-  // end of loop on calo clusters
-  
-  // fill EMCAL cell info
-  if (fCalorimeter.Contains("EMCAL") && aod->GetEMCALCells()) { // protection against missing ESD information
-    AliAODCaloCells &aodinEMcells = *(aod->GetEMCALCells());
-    Int_t nEMcell = aodinEMcells.GetNumberOfCells() ;
-    
-    AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
-    aodEMcells.CreateContainer(nEMcell);
-    aodEMcells.SetType(AliVCaloCells::kEMCALCell);
-    for (Int_t iCell = 0; iCell < nEMcell; iCell++) {      
-      aodEMcells.SetCell(iCell,aodinEMcells.GetCellNumber(iCell),aodinEMcells.GetAmplitude(iCell));
-    }
-    aodEMcells.Sort();
-  }
-  
-  // fill PHOS cell info
-  if (fCalorimeter.Contains("PHOS") && aod->GetPHOSCells()) { // protection against missing ESD information
-    AliAODCaloCells &aodinPHcells = *(aod->GetPHOSCells());
-    Int_t nPHcell = aodinPHcells.GetNumberOfCells() ;
-    
-    AliAODCaloCells &aodPHcells = *(AODEvent()->GetPHOSCells());
-    aodPHcells.CreateContainer(nPHcell);
-    aodPHcells.SetType(AliVCaloCells::kPHOSCell);
-    for (Int_t iCell = 0; iCell < nPHcell; iCell++) {      
-      aodPHcells.SetCell(iCell,aodinPHcells.GetCellNumber(iCell),aodinPHcells.GetAmplitude(iCell));
-    }
-    aodPHcells.Sort();
-  }
-  
-  
-  return;
-}
-
-//__________________________________________________
-void AliAnalysisTaskCaloFilter::CreateAODFromESD()
-{
-
-  // Copy ESD header, vertex, CaloClusters and CaloCells to output AOD
-  
-  AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
-  if(!esd) {
-    printf("AliAnalysisTaskCaloFilter::CreateAODFromESD() - This event does not contain AODs?");
+  if(!event) {
+    printf("AliAnalysisTaskCaloFilter::CreateAODFromESD() - This event does not contain Input?");
     return;
   }
   
@@ -210,34 +111,42 @@ void AliAnalysisTaskCaloFilter::CreateAODFromESD()
   Double_t covVtx[6];
   
   for (Int_t i = 0; i < 6; i++)  covVtx[i] = 0.;
-  
-  // Update the header
-  
+      
   AliAODHeader* header = AODEvent()->GetHeader();
-  header->SetRunNumber(esd->GetRunNumber());
-  header->SetBunchCrossNumber(esd->GetBunchCrossNumber());
-  header->SetOrbitNumber(esd->GetOrbitNumber());
-  header->SetPeriodNumber(esd->GetPeriodNumber());
-  header->SetEventType(esd->GetEventType());
+  
+  header->SetRunNumber(event->GetRunNumber());
+  header->SetOfflineTrigger(fInputHandler->IsEventSelected()); // propagate the decision of the physics selection
+  if(bESD)
+    header->SetFiredTriggerClasses(((AliESDEvent*)event)->GetFiredTriggerClasses());
+  header->SetTriggerMask(event->GetTriggerMask()); 
+  header->SetTriggerCluster(event->GetTriggerCluster());
+  
+  header->SetBunchCrossNumber(event->GetBunchCrossNumber());
+  header->SetOrbitNumber(event->GetOrbitNumber());
+  header->SetPeriodNumber(event->GetPeriodNumber());
+  header->SetEventType(event->GetEventType());
   header->SetMuonMagFieldScale(-999.); // FIXME
   header->SetCentrality(-999.);        // FIXME
   
   
-  header->SetTriggerMask(esd->GetTriggerMask()); 
-  header->SetTriggerCluster(esd->GetTriggerCluster());
-  header->SetMagneticField(esd->GetMagneticField());
-  header->SetZDCN1Energy(esd->GetZDCN1Energy());
-  header->SetZDCP1Energy(esd->GetZDCP1Energy());
-  header->SetZDCN2Energy(esd->GetZDCN2Energy());
-  header->SetZDCP2Energy(esd->GetZDCP2Energy());
-  header->SetZDCEMEnergy(esd->GetZDCEMEnergy(0),esd->GetZDCEMEnergy(1));
-  Float_t diamxy[2]={esd->GetDiamondX(),esd->GetDiamondY()};
-  Float_t diamcov[3]; esd->GetDiamondCovXY(diamcov);
+  header->SetTriggerMask(event->GetTriggerMask()); 
+  header->SetTriggerCluster(event->GetTriggerCluster());
+  header->SetMagneticField(event->GetMagneticField());
+  header->SetZDCN1Energy(event->GetZDCN1Energy());
+  header->SetZDCP1Energy(event->GetZDCP1Energy());
+  header->SetZDCN2Energy(event->GetZDCN2Energy());
+  header->SetZDCP2Energy(event->GetZDCP2Energy());
+  header->SetZDCEMEnergy(event->GetZDCEMEnergy(0),event->GetZDCEMEnergy(1));
+  Float_t diamxy[2]={event->GetDiamondX(),event->GetDiamondY()};
+  Float_t diamcov[3]; event->GetDiamondCovXY(diamcov);
   header->SetDiamond(diamxy,diamcov);
+  if(bESD){
+    header->SetDiamondZ(((AliESDEvent*)event)->GetDiamondZ(),((AliESDEvent*)event)->GetSigma2DiamondZ());
+  }
   //
   //
   Int_t nVertices = 1 ;/* = prim. vtx*/;
-  Int_t nCaloClus = esd->GetNumberOfCaloClusters();
+  Int_t nCaloClus = event->GetNumberOfCaloClusters();
   
   AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
   
@@ -247,15 +156,21 @@ void AliAnalysisTaskCaloFilter::CreateAODFromESD()
   
   // Add primary vertex. The primary tracks will be defined
   // after the loops on the composite objects (V0, cascades, kinks)
-  const AliESDVertex *vtx = esd->GetPrimaryVertex();
-  
-  vtx->GetXYZ(pos); // position
-  vtx->GetCovMatrix(covVtx); //covariance matrix
+  event->GetPrimaryVertex()->GetXYZ(pos);
+  Float_t chi = 0;
+  if      (bESD){
+    ((AliESDEvent*)event)->GetPrimaryVertex()->GetCovMatrix(covVtx);
+    chi = ((AliESDEvent*)event)->GetPrimaryVertex()->GetChi2toNDF();
+  }
+  else if (bAOD){
+    ((AliAODEvent*)event)->GetPrimaryVertex()->GetCovMatrix(covVtx);
+    chi = ((AliAODEvent*)event)->GetPrimaryVertex()->GetChi2perNDF();//Different from ESD?
+  }
   
   AliAODVertex * primary = new(vertices[jVertices++])
-    AliAODVertex(pos, covVtx, vtx->GetChi2toNDF(), NULL, -1, AliAODVertex::kPrimary);
-  primary->SetName(vtx->GetName());
-  primary->SetTitle(vtx->GetTitle());
+    AliAODVertex(pos, covVtx, chi, NULL, -1, AliAODVertex::kPrimary);
+  primary->SetName(event->GetPrimaryVertex()->GetName());
+  primary->SetTitle(event->GetPrimaryVertex()->GetTitle());
   
   // Access to the AOD container of clusters
   TClonesArray &caloClusters = *(AODEvent()->GetCaloClusters());
@@ -263,12 +178,48 @@ void AliAnalysisTaskCaloFilter::CreateAODFromESD()
   
   for (Int_t iClust=0; iClust<nCaloClus; ++iClust) {
     
-    AliESDCaloCluster * cluster = esd->GetCaloCluster(iClust);
+    AliVCluster * cluster = event->GetCaloCluster(iClust);
     
     //Check which calorimeter information we want to keep.
-    if     (fCalorimeter.Contains("PHOS")  && !fCalorimeter.Contains("EMCAL") && cluster->IsEMCAL()) continue ;
-    else if(fCalorimeter.Contains("EMCAL") && !fCalorimeter.Contains("PHOS")  && cluster->IsPHOS())  continue ;
     
+    if(fCaloFilter!=kBoth){
+      if     (fCaloFilter==kPHOS  && cluster->IsEMCAL()) continue ;
+      else if(fCaloFilter==kEMCAL && cluster->IsPHOS())  continue ;
+    }  
+    
+    //--------------------------------------------------------------
+    //If EMCAL, and requested, correct energy, position ...
+    //printf("Hello IsEMCAL? %d, Correct? %d\n",cluster->IsEMCAL(), fCorrect);
+    if(cluster->IsEMCAL() && fCorrect){
+      Float_t position[]={0,0,0};
+      if(DebugLevel() > 2)
+        printf("Check cluster %d for bad channels and close to border\n",cluster->GetID());
+      if(fEMCALRecoUtils->ClusterContainsBadChannel(fEMCALGeo,cluster->GetCellsAbsId(), cluster->GetNCells())) continue;       
+      if(fEMCALRecoUtils->CheckCellFiducialRegion(fEMCALGeo, cluster, event->GetEMCALCells())) continue;
+      if(DebugLevel() > 2)
+      { 
+        printf("Filter, before  : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",iClust,cluster->E(),
+               cluster->GetDispersion(),cluster->GetM02(),cluster->GetM20());
+        cluster->GetPosition(position);
+        printf("Filter, before  : i %d, x %f, y %f, z %f\n",cluster->GetID(), position[0], position[1], position[2]);
+      }
+            
+      if(fEMCALRecoUtils->IsRecalibrationOn()) fEMCALRecoUtils->RecalibrateClusterEnergy(fEMCALGeo, cluster, event->GetEMCALCells());
+      cluster->SetE(fEMCALRecoUtils->CorrectClusterEnergyLinearity(cluster));
+      fEMCALRecoUtils->RecalculateClusterPosition(fEMCALGeo, event->GetEMCALCells(),cluster);
+      
+      if(DebugLevel() > 2)
+      { 
+        printf("Filter, after   : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",cluster->GetID(),cluster->E(),
+               cluster->GetDispersion(),cluster->GetM02(),cluster->GetM20());
+        cluster->GetPosition(position);
+        printf("Filter, after   : i %d, x %f, y %f, z %f\n",cluster->GetID(), position[0], position[1], position[2]);
+      }    
+      
+    }
+    //--------------------------------------------------------------
+
+    //Now fill AODs
     Int_t id       = cluster->GetID();
     Float_t energy = cluster->E();
     cluster->GetPosition(posF);
@@ -293,34 +244,49 @@ void AliAnalysisTaskCaloFilter::CreateAODFromESD()
     caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
     caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
     
+    if(DebugLevel() > 2)
+    { 
+      printf("Filter, aod       : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",caloCluster->GetID(),caloCluster->E(),
+             caloCluster->GetDispersion(),caloCluster->GetM02(),caloCluster->GetM20());
+      caloCluster->GetPosition(posF);
+      printf("Filter, aod       : i %d, x %f, y %f, z %f\n",caloCluster->GetID(), posF[0], posF[1], posF[2]);
+    }    
+    
   } 
   caloClusters.Expand(jClusters); // resize TObjArray to 'remove' slots for pseudo clusters     
   // end of loop on calo clusters
   
   // fill EMCAL cell info
-  if (fCalorimeter.Contains("EMCAL") && esd->GetEMCALCells()) { // protection against missing ESD information
-    AliESDCaloCells &esdEMcells = *(esd->GetEMCALCells());
-    Int_t nEMcell = esdEMcells.GetNumberOfCells() ;
+  if ((fCaloFilter==kBoth ||  fCaloFilter==kEMCAL) && event->GetEMCALCells()) { // protection against missing ESD information
+    AliVCaloCells &eventEMcells = *(event->GetEMCALCells());
+    Int_t nEMcell = eventEMcells.GetNumberOfCells() ;
     
     AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
     aodEMcells.CreateContainer(nEMcell);
     aodEMcells.SetType(AliVCaloCells::kEMCALCell);
-    for (Int_t iCell = 0; iCell < nEMcell; iCell++) {      
-      aodEMcells.SetCell(iCell,esdEMcells.GetCellNumber(iCell),esdEMcells.GetAmplitude(iCell));
+    Double_t calibFactor = 1.;   
+    for (Int_t iCell = 0; iCell < nEMcell; iCell++) { 
+      if(fCorrect && fEMCALRecoUtils->IsRecalibrationOn()){ 
+        Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1; 
+        fEMCALGeo->GetCellIndex(eventEMcells.GetCellNumber(iCell),imod,iTower,iIphi,iIeta); 
+        fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);    
+        calibFactor = fEMCALRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
+      }
+      aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),eventEMcells.GetAmplitude(iCell)*calibFactor);
     }
     aodEMcells.Sort();
   }
   
   // fill PHOS cell info
-  if (fCalorimeter.Contains("PHOS") && esd->GetPHOSCells()) { // protection against missing ESD information
-    AliESDCaloCells &esdPHcells = *(esd->GetPHOSCells());
-    Int_t nPHcell = esdPHcells.GetNumberOfCells() ;
+  if ((fCaloFilter==kBoth ||  fCaloFilter==kPHOS) && event->GetPHOSCells()) { // protection against missing ESD information
+    AliVCaloCells &eventPHcells = *(event->GetPHOSCells());
+    Int_t nPHcell = eventPHcells.GetNumberOfCells() ;
     
     AliAODCaloCells &aodPHcells = *(AODEvent()->GetPHOSCells());
     aodPHcells.CreateContainer(nPHcell);
     aodPHcells.SetType(AliVCaloCells::kPHOSCell);
     for (Int_t iCell = 0; iCell < nPHcell; iCell++) {      
-      aodPHcells.SetCell(iCell,esdPHcells.GetCellNumber(iCell),esdPHcells.GetAmplitude(iCell));
+      aodPHcells.SetCell(iCell,eventPHcells.GetCellNumber(iCell),eventPHcells.GetAmplitude(iCell));
     }
     aodPHcells.Sort();
   }
@@ -329,39 +295,26 @@ void AliAnalysisTaskCaloFilter::CreateAODFromESD()
   return;
 }
 
-//__________________________________________________
-void AliAnalysisTaskCaloFilter::Init()
-{
-  // Initialization
-  if (fDebug > 1) AliInfo("Init() \n");
-  
-}
-
-//__________________________________________________
-void AliAnalysisTaskCaloFilter::UserCreateOutputObjects()
-{
-  // Create the output container
-}
-
-//__________________________________________________
-void AliAnalysisTaskCaloFilter::UserExec(Option_t */*option*/)
-{
-  // Execute analysis for current event
-  //
-  
-  Long64_t ientry = Entry();
-  if (fDebug > 0)  printf("CaloFilter: Analysing event # %5d\n", (Int_t) ientry);
-
-  const char * dataevent =InputEvent()->GetName();
-  if(!strcmp(dataevent,"AliAODEvent"))   CreateAODFromAOD();  
-  
-  else if(!strcmp(dataevent,"AliESDEvent"))  CreateAODFromESD();
-  else {
-    AliFatal(Form("Unknown event type %s, ABORT",dataevent));
-    
-  }
+//_____________________________________________________
+//void AliAnalysisTaskCaloFilter::LocalInit()
+//{
+//     // Local Initialization
+//     
+//     // Create cuts/param objects and publish to slot
+//     const Int_t buffersize = 255;
+//     char onePar[buffersize] ;
+//     fCuts = new TList();
+//  
+//     snprintf(onePar,buffersize, "Calorimeter Filtering Option %d", fCaloFilter) ;
+//     fCuts->Add(new TObjString(onePar));
+//     snprintf(onePar,buffersize, "Not only filter but correct? %d cells;", fCorrect) ;
+//     fCuts->Add(new TObjString(onePar));
+//  
+//     // Post Data
+//     PostData(2, fCuts);
+//     
+//}
 
-}
 
 //__________________________________________________
 void AliAnalysisTaskCaloFilter::Terminate(Option_t */*option*/)
index 201364e..6237d2a 100644 (file)
 // Author: Gustavo Conesa Balbastre (INFN - Frascati)
 //////////////////////////////////////////////////////////
 
+class TList;
+
 #include "AliAnalysisTaskSE.h"
+class AliEMCALRecoUtils;
+class AliEMCALGeometry;
 
 class AliAnalysisTaskCaloFilter : public AliAnalysisTaskSE
 {
  public:
   AliAnalysisTaskCaloFilter();
   AliAnalysisTaskCaloFilter(const char* name);
-  virtual ~AliAnalysisTaskCaloFilter() {;}
+  virtual ~AliAnalysisTaskCaloFilter() ;
   
-  // Implementation of interface methods
+private:
+  AliAnalysisTaskCaloFilter(const AliAnalysisTaskCaloFilter&);
+  AliAnalysisTaskCaloFilter& operator=(const AliAnalysisTaskCaloFilter&);
+  
+public:
   virtual void   UserCreateOutputObjects();
-  virtual void   Init();
-  virtual void   LocalInit() {Init();}
+  //virtual void   Init();
+  //virtual void   LocalInit() ;
   virtual void   UserExec(Option_t *option);
   virtual void   Terminate(Option_t *option);
   
-  void SetCalorimeter(TString calo) {fCalorimeter = calo;}
-  TString GetCalorimeter() const {return fCalorimeter;}
-  
-  void CreateAODFromESD();
-  void CreateAODFromAOD();
+  enum caloFilter {kBoth = 0, kEMCAL = 1, kPHOS=2};
+  void SetCaloFilter(Int_t calo) {fCaloFilter = calo;}
+  TString GetCaloFilter() const  {return fCaloFilter;}
 
- private:
-  AliAnalysisTaskCaloFilter(const AliAnalysisTaskCaloFilter&);
-  AliAnalysisTaskCaloFilter& operator=(const AliAnalysisTaskCaloFilter&);
+  void SetGeometryName(TString name) { fEMCALGeoName = name ; }
+  TString GeometryName() const       { return fEMCALGeoName ; }
   
-  TString fCalorimeter; //Calorimeter to filter
+  void SetEMCALRecoUtils(AliEMCALRecoUtils * ru) {fEMCALRecoUtils = ru;}
+  AliEMCALRecoUtils* GetEMCALRecoUtils() const   {return fEMCALRecoUtils;}
+
+  void SwitchOnClusterCorrection()  {fCorrect = kTRUE ;}
+  void SwitchOffClusterCorrection() {fCorrect = kFALSE;}
+
+private:
   
-  ClassDef(AliAnalysisTaskCaloFilter, 1); // Analysis task for standard ESD filtering
+  //TList* fCuts ;      //! List with analysis cuts
+  Int_t  fCaloFilter; // Calorimeter to filter
+  Int_t  fCorrect;    // Recalibrate or recalculate different cluster parameters
+  //EMCAL specific
+  AliEMCALGeometry  * fEMCALGeo;       //! EMCAL geometry
+  TString             fEMCALGeoName;   // Name of geometry to use.
+  AliEMCALRecoUtils * fEMCALRecoUtils; // Pointer to EMCAL utilities for clusterization
+
+  ClassDef(AliAnalysisTaskCaloFilter, 2); // Analysis task for standard ESD filtering
 };
 
 #endif
index 15cabfb..78df0f4 100644 (file)
 //#include <Riostream.h>
 // Root 
 #include "TLorentzVector.h"
-//#include "TVector3.h"
 #include "TRefArray.h"
 #include "TList.h"
 #include "TH1F.h"
+#include <TGeoManager.h>
 
 // AliRoot
 #include "AliAnalysisTaskEMCALPi0CalibSelection.h"
@@ -39,7 +39,6 @@
 #include "AliVCluster.h"
 #include "AliVCaloCells.h"
 #include "AliEMCALRecoUtils.h"
-
 //#include "AliEMCALAodCluster.h"
 //#include "AliEMCALCalibData.h"
 
@@ -50,9 +49,8 @@ ClassImp(AliAnalysisTaskEMCALPi0CalibSelection)
 AliAnalysisTaskEMCALPi0CalibSelection::AliAnalysisTaskEMCALPi0CalibSelection(const char* name) :
   AliAnalysisTaskSE(name),fEMCALGeo(0x0),//fCalibData(0x0), 
   fEmin(0.5), fEmax(15.), fAsyCut(1.),fMinNCells(2), fGroupNCells(0),
-  fLogWeight(4.5), fCopyAOD(kFALSE), fSameSM(kFALSE), fOldAOD(kFALSE),
-  fEMCALGeoName("EMCAL_FIRSTYEAR"), fNCellsFromEMCALBorder(0),
-  fRemoveBadChannels(kFALSE),fEMCALBadChannelMap(0x0),
+  fLogWeight(4.5), fSameSM(kFALSE), fOldAOD(kFALSE), fFilteredInput(kFALSE),
+  fCorrectClusters(kFALSE), fEMCALGeoName("EMCAL_FIRSTYEAR"), 
   fRecoUtils(new AliEMCALRecoUtils),
   fNbins(300), fMinBin(0.), fMaxBin(300.),fOutputContainer(0x0),
   fHmgg(0x0),           fHmggDifferentSM(0x0), 
@@ -103,11 +101,6 @@ AliAnalysisTaskEMCALPi0CalibSelection::~AliAnalysisTaskEMCALPi0CalibSelection()
   //if(fCalibData)  delete fCalibData;
   if(fEMCALGeo)   delete fEMCALGeo;
        
-       if(fEMCALBadChannelMap) { 
-               fEMCALBadChannelMap->Clear();
-               delete  fEMCALBadChannelMap;
-       }
-       
   if(fRecoUtils) delete fRecoUtils ;
 
 }
@@ -126,12 +119,12 @@ void AliAnalysisTaskEMCALPi0CalibSelection::LocalInit()
        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;", fNCellsFromEMCALBorder) ;
+  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) ;
        fCuts->Add(new TObjString(onePar));
-       snprintf(onePar,buffersize, "Switchs: Remove Bad Channels? %d; Copy AODs? %d;  Recalibrate %d?, Analyze Old AODs? %d, Mass per channel same SM clusters? %d ",
-            fRemoveBadChannels,fCopyAOD,fRecoUtils->IsRecalibrationOn(), fOldAOD, fSameSM) ;
+       snprintf(onePar,buffersize, "Switchs: Remove Bad Channels? %d; Use filtered input? %d;  Correct Clusters? %d, Analyze Old AODs? %d, Mass per channel same SM clusters? %d ",
+            fRecoUtils->IsBadChannelsRemovalSwitchedOn(),fFilteredInput,fCorrectClusters, fOldAOD, fSameSM) ;
        fCuts->Add(new TObjString(onePar));
        snprintf(onePar,buffersize, "EMCAL Geometry name: < %s >",fEMCALGeoName.Data()) ;
        fCuts->Add(new TObjString(onePar));
@@ -141,264 +134,6 @@ void AliAnalysisTaskEMCALPi0CalibSelection::LocalInit()
        
 }
 
-//__________________________________________________
-void AliAnalysisTaskEMCALPi0CalibSelection::CreateAODFromAOD()
-{
-  // Copy AOD header, vertex, CaloClusters and CaloCells to output AOD
-  AliAODEvent* aod = dynamic_cast<AliAODEvent*>(InputEvent());
-  
-  if(!aod) {
-  printf("AliAnalysisTaskEMCALPi0CalibSelection::CreateAODFromAOD() - This event does not contain AODs?");
-    return;
-  }
-  
-  // set arrays and pointers
-  Float_t posF[3];
-  Double_t pos[3];
-  
-  Double_t covVtx[6];
-  
-  for (Int_t i = 0; i < 6; i++)  covVtx[i] = 0.;
-  
-  // Update the header
-  AliAODHeader*headerin = aod->GetHeader();
-  AliAODHeader* header = AODEvent()->GetHeader();
-  header->SetRunNumber(headerin->GetRunNumber());
-  header->SetBunchCrossNumber(headerin->GetBunchCrossNumber());
-  header->SetOrbitNumber(headerin->GetOrbitNumber());
-  header->SetPeriodNumber(headerin->GetPeriodNumber());
-  header->SetEventType(headerin->GetEventType());
-  header->SetMuonMagFieldScale(headerin->GetMuonMagFieldScale());
-  header->SetCentrality(headerin->GetCentrality()); 
-  
-  header->SetTriggerMask(headerin->GetTriggerMask()); 
-  header->SetTriggerCluster(headerin->GetTriggerCluster());
-  header->SetMagneticField(headerin->GetMagneticField());
-  header->SetZDCN1Energy(headerin->GetZDCN1Energy());
-  header->SetZDCP1Energy(headerin->GetZDCP1Energy());
-  header->SetZDCN2Energy(headerin->GetZDCN2Energy());
-  header->SetZDCP2Energy(headerin->GetZDCP2Energy());
-  header->SetZDCEMEnergy(headerin->GetZDCEMEnergy(0),headerin->GetZDCEMEnergy(1));
-  Float_t diamxy[2]={aod->GetDiamondX(),aod->GetDiamondY()};
-  Float_t diamcov[3]; aod->GetDiamondCovXY(diamcov);
-  header->SetDiamond(diamxy,diamcov);
-  //
-  //
-  Int_t nVertices = 1 ;/* = prim. vtx*/;
-  Int_t nCaloClus = aod->GetNumberOfCaloClusters();
-  
-  AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
-  
-  // Access to the AOD container of vertices
-  TClonesArray &vertices = *(AODEvent()->GetVertices());
-  Int_t jVertices=0;
-  
-  // Add primary vertex. The primary tracks will be defined
-  // after the loops on the composite objects (V0, cascades, kinks)
-  const AliAODVertex *vtx = aod->GetPrimaryVertex();
-  
-  vtx->GetXYZ(pos); // position
-  vtx->GetCovMatrix(covVtx); //covariance matrix
-  
-  AliAODVertex * primary = new(vertices[jVertices++])
-  AliAODVertex(pos, covVtx, vtx->GetChi2perNDF(), NULL, -1, AliAODVertex::kPrimary);
-  primary->SetName(vtx->GetName());
-  primary->SetTitle(vtx->GetTitle());
-  
-  // Access to the AOD container of clusters
-  TClonesArray &caloClusters = *(AODEvent()->GetCaloClusters());
-  Int_t jClusters=0;
-  //printf("nCaloClus %d\n",nCaloClus);
-  
-  for (Int_t iClust=0; iClust<nCaloClus; ++iClust) {
-    
-    AliAODCaloCluster * cluster = aod->GetCaloCluster(iClust);
-    
-    //Check if it is a EMCAL cluster
-    if(!IsEMCALCluster(cluster))  continue ;
-    //printf("EMCAL cluster %d, ncells %d\n",iClust, cluster->GetNCells());
-    //if(ClusterContainsBadChannel(cluster->GetCellsAbsId(), cluster->GetNCells())) continue;  
-    //printf("copy\n");
-    Int_t id       = cluster->GetID();
-    Float_t energy = cluster->E();
-    cluster->GetPosition(posF);
-    Char_t ttype   = cluster->GetType(); 
-    AliAODCaloCluster *caloCluster = new(caloClusters[jClusters++]) 
-    AliAODCaloCluster(id,
-                      0,
-                      0x0,
-                      energy,
-                      posF,
-                      NULL,
-                      ttype);
-    
-    caloCluster->SetCaloCluster(cluster->GetDistanceToBadChannel(),
-                                cluster->GetDispersion(),
-                                cluster->GetM20(), cluster->GetM02(),
-                                cluster->GetEmcCpvDistance(),  
-                                cluster->GetNExMax(),cluster->GetTOF()) ;
-    
-    caloCluster->SetPIDFromESD(cluster->GetPID());
-    caloCluster->SetNCells(cluster->GetNCells());
-    caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
-    
-    caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
-    
-  } 
-  
-  caloClusters.Expand(jClusters); // resize TObjArray   
-  // end of loop on calo clusters
-  
-  // fill EMCAL cell info
-  if (aod->GetEMCALCells()) { // protection against missing AOD information
-    AliAODCaloCells &aodinEMcells = *(aod->GetEMCALCells());
-    Int_t nEMcell = aodinEMcells.GetNumberOfCells() ;
-    
-    AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
-    aodEMcells.CreateContainer(nEMcell);
-    aodEMcells.SetType(AliAODCaloCells::kEMCALCell);
-    
-    Double_t calibFactor = 1;
-    for (Int_t iCell = 0; iCell < nEMcell; iCell++) {      
-      aodEMcells.SetCell(iCell,aodinEMcells.GetCellNumber(iCell),aodinEMcells.GetAmplitude(iCell)*calibFactor);
-    }
-    aodEMcells.Sort();
-         
-  }
-  
-}
-
-//__________________________________________________
-void AliAnalysisTaskEMCALPi0CalibSelection::CreateAODFromESD()
-{
-  
-  // Copy Header, Vertex, CaloClusters and CaloCells from ESDs to AODs
-  AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
-  
-  if(!esd) {
-    printf("AliAnalysisTaskEMCALPi0CalibSelection::CreateAODFromAOD() - This event does not contain AODs?");
-    return;
-  }
-  
-  // set arrays and pointers
-  Float_t posF[3];
-  Double_t pos[3];
-  
-  Double_t covVtx[6];
-  
-  for (Int_t i = 0; i < 6; i++)  covVtx[i] = 0.;
-  
-  // Update the header
-  
-  AliAODHeader* header = AODEvent()->GetHeader();
-  header->SetRunNumber(esd->GetRunNumber());
-  header->SetBunchCrossNumber(esd->GetBunchCrossNumber());
-  header->SetOrbitNumber(esd->GetOrbitNumber());
-  header->SetPeriodNumber(esd->GetPeriodNumber());
-  header->SetEventType(esd->GetEventType());
-  header->SetMuonMagFieldScale(-999.); // FIXME
-  header->SetCentrality(-999.);        // FIXME
-  
-  
-  header->SetTriggerMask(esd->GetTriggerMask()); 
-  header->SetTriggerCluster(esd->GetTriggerCluster());
-  header->SetMagneticField(esd->GetMagneticField());
-  header->SetZDCN1Energy(esd->GetZDCN1Energy());
-  header->SetZDCP1Energy(esd->GetZDCP1Energy());
-  header->SetZDCN2Energy(esd->GetZDCN2Energy());
-  header->SetZDCP2Energy(esd->GetZDCP2Energy());
-  header->SetZDCEMEnergy(esd->GetZDCEMEnergy(0),esd->GetZDCEMEnergy(1));
-  Float_t diamxy[2]={esd->GetDiamondX(),esd->GetDiamondY()};
-  Float_t diamcov[3]; esd->GetDiamondCovXY(diamcov);
-  header->SetDiamond(diamxy,diamcov);
-  //
-  //
-  Int_t nVertices = 1 ;/* = prim. vtx*/;
-  Int_t nCaloClus = esd->GetNumberOfCaloClusters();
-  
-  AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
-  
-  // Access to the AOD container of vertices
-  TClonesArray &vertices = *(AODEvent()->GetVertices());
-  Int_t jVertices=0;
-  
-  // Add primary vertex. The primary tracks will be defined
-  // after the loops on the composite objects (V0, cascades, kinks)
-  const AliESDVertex *vtx = esd->GetPrimaryVertex();
-  
-  vtx->GetXYZ(pos); // position
-  vtx->GetCovMatrix(covVtx); //covariance matrix
-  
-  AliAODVertex * primary = new(vertices[jVertices++])
-  AliAODVertex(pos, covVtx, vtx->GetChi2toNDF(), NULL, -1, AliAODVertex::kPrimary);
-  primary->SetName(vtx->GetName());
-  primary->SetTitle(vtx->GetTitle());
-  
-  // Access to the AOD container of clusters
-  TClonesArray &caloClusters = *(AODEvent()->GetCaloClusters());
-  Int_t jClusters=0;
-  //printf("nCaloClus %d\n",nCaloClus);
-  
-  for (Int_t iClust=0; iClust<nCaloClus; ++iClust) {
-    
-    AliESDCaloCluster * cluster = esd->GetCaloCluster(iClust);
-    
-    //Check which calorimeter information we want to keep.
-    if(!IsEMCALCluster(cluster))  continue ;
-    //printf("EMCAL cluster %d, ncells %d\n",iClust, cluster->GetNCells());
-    
-    if(ClusterContainsBadChannel(cluster->GetCellsAbsId(), cluster->GetNCells())) continue;    
-    //printf("copy\n");
-    
-    Int_t id       = cluster->GetID();
-    Float_t energy = cluster->E();
-    cluster->GetPosition(posF);
-    
-    AliAODCaloCluster *caloCluster = new(caloClusters[jClusters++]) 
-    AliAODCaloCluster(id,
-                      0,
-                      0x0,
-                      energy,
-                      posF,
-                      NULL,
-                      AliVCluster::kEMCALClusterv1);
-    
-    caloCluster->SetCaloCluster(cluster->GetDistanceToBadChannel(),
-                                cluster->GetDispersion(),
-                                cluster->GetM20(), cluster->GetM02(),
-                                cluster->GetEmcCpvDistance(),  
-                                cluster->GetNExMax(),cluster->GetTOF()) ;
-    
-    caloCluster->SetPIDFromESD(cluster->GetPID());
-    caloCluster->SetNCells(cluster->GetNCells());
-    caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
-    caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
-    
-  } 
-  
-  caloClusters.Expand(jClusters); // resize TObjArray
-  // end of loop on calo clusters
-  
-  // fill EMCAL cell info
-  
-  if( esd->GetEMCALCells()) { // protection against missing ESD information
-    AliESDCaloCells &esdEMcells = *(esd->GetEMCALCells());
-    Int_t nEMcell = esdEMcells.GetNumberOfCells() ;
-    
-    AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
-    aodEMcells.CreateContainer(nEMcell);
-    aodEMcells.SetType(AliAODCaloCells::kEMCALCell);  
-         
-    Double_t calibFactor = 1;   
-    for (Int_t iCell = 0; iCell < nEMcell; iCell++) {      
-      aodEMcells.SetCell(iCell,esdEMcells.GetCellNumber(iCell),esdEMcells.GetAmplitude(iCell)*calibFactor);
-    }
-    aodEMcells.Sort();
-         
-  }
-
-}
-
 //_________________________________________________________________
 Int_t AliAnalysisTaskEMCALPi0CalibSelection::GetEMCALClusters(AliVEvent * event, TRefArray *clusters) const
 {
@@ -444,7 +179,7 @@ Bool_t AliAnalysisTaskEMCALPi0CalibSelection::IsEMCALCluster(AliVCluster* cluste
 //__________________________________________________
 void AliAnalysisTaskEMCALPi0CalibSelection::UserCreateOutputObjects()
 {
-  //Create output container, init geometry and calibration
+  //Create output container, init geometry 
        
   fEMCALGeo =  AliEMCALGeometry::GetInstance(fEMCALGeoName) ;  
   
@@ -589,11 +324,11 @@ void AliAnalysisTaskEMCALPi0CalibSelection::UserCreateOutputObjects()
     
   }  
   
-  
-  
   fhNEvents        = new TH1I("hNEvents", "Number of analyzed events"   , 1 , 0 , 1  ) ;
   fOutputContainer->Add(fhNEvents);
-
+  
+  fOutputContainer->SetOwner(kTRUE);
+  
 //  fCalibData = new AliEMCALCalibData();
                
   PostData(1,fOutputContainer);
@@ -604,58 +339,38 @@ void AliAnalysisTaskEMCALPi0CalibSelection::UserCreateOutputObjects()
 void AliAnalysisTaskEMCALPi0CalibSelection::UserExec(Option_t* /* option */)
 {
   //Analysis per event.
-  if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection <<< Event %d >>>\n",(Int_t)Entry());
   
   if(fRecoUtils->GetParticleType()!=AliEMCALRecoUtils::kPhoton){
     printf("Wrong particle type for cluster position recalculation! = %d\n", fRecoUtils->GetParticleType());
     abort();
   }
-    
+  
   fhNEvents->Fill(0); //Event analyzed
   
-  AliAODEvent* aod = 0x0;
-  Bool_t kAOD = kFALSE;
-  if(!strcmp(InputEvent()->GetName(),"AliAODEvent")) kAOD=kTRUE;
-  Bool_t kESD = kFALSE;
-  if(!strcmp(InputEvent()->GetName(),"AliESDEvent")) kESD=kTRUE;
+  //Get the input event
+  AliVEvent* event = 0;
+  if(fFilteredInput) event = AODEvent();
+  else               event = InputEvent();
   
-  if(kAOD){
-    //Input are ESDs
-    aod = dynamic_cast<AliAODEvent*>(InputEvent());
-    if(!aod) {
-      printf("AliAnalysisTaskEMCALPi0CalibSelection::UserExec() - This event does not contain AODs?");
-      return;
-    }
-    
-    // Create new AOD with only CaloClusters and CaloCells
-    if(fCopyAOD) CreateAODFromAOD();
-  }
-  else  if(kESD) {
-    //Input are ESDs
-    aod = AODEvent();
-    if(!aod) {
-      printf("AliAnalysisTaskEMCALPi0CalibSelection::UserExec() - This event does not contain AODs?");
-      return;
-    }
-    
-    // Create AOD with CaloClusters and use it as input.
-    // If filtering task is already executed, this is not needed.
-    if(fCopyAOD) CreateAODFromESD();
+  
+  
+  if(!event) {
+    printf("Input event not available!\n");
+    return;
   }
-  else {
-    printf("AliAnalysisTaskEMCALPi0CalibSelection: Unknown event type, STOP!\n");
-    abort();
-  }    
   
-  Double_t v[3];// = {aod->GetVertex(0)->GetX(),aod->GetVertex(0)->GetY(),aod->GetVertex(0)->GetZ()}; //to check!!
-  aod->GetPrimaryVertex()->GetXYZ(v) ;
-  //TVector3 vtx(v); 
+  if(DebugLevel() > 1) 
+    printf("AliAnalysisTaskEMCALPi0CalibSelection <<< %s: Event %d >>>\n",event->GetName(), (Int_t)Entry());
+  
+  
+  //Get the primary vertex
+  Double_t v[3];
+  event->GetPrimaryVertex()->GetXYZ(v) ;
   
-  //if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Vertex: (%.3f,%.3f,%.3f)\n",vtx.X(),vtx.Y(),vtx.Z());
   if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Vertex: (%.3f,%.3f,%.3f)\n",v[0],v[1],v[2]);
   
-  Int_t runNum = aod->GetRunNumber();
-  if(DebugLevel() > 1) printf("Run number: %d\n",runNum);
+  //Int_t runNum = aod->GetRunNumber();
+  //if(DebugLevel() > 1) printf("Run number: %d\n",runNum);
   
   //FIXME Not need the matrices for the moment MEFIX
   //Get the matrix with geometry information
@@ -688,26 +403,61 @@ void AliAnalysisTaskEMCALPi0CalibSelection::UserExec(Option_t* /* option */)
        
   TLorentzVector p1;
   TLorentzVector p2;
-//  TLorentzVector p11;
-//  TLorentzVector p22;
-
   TLorentzVector p12;
   
+  //Get the list of clusters
   TRefArray * caloClustersArr  = new TRefArray();
-  if(!fOldAOD) aod->GetEMCALClusters(caloClustersArr);
-  else  GetEMCALClusters(aod,caloClustersArr);
-    
+  if(!fOldAOD) event->GetEMCALClusters(caloClustersArr);
+  else                GetEMCALClusters(event,caloClustersArr);
   const Int_t kNumberOfEMCALClusters   = caloClustersArr->GetEntries() ;
   if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection - N CaloClusters: %d \n", kNumberOfEMCALClusters);
   
-  // EMCAL cells
-  AliAODCaloCells *emCells = aod->GetEMCALCells();
-   
+  // Get EMCAL cells
+  AliVCaloCells *emCells = event->GetEMCALCells();
+  
   // loop over EMCAL clusters
-  for(Int_t iClu=0; iClu<kNumberOfEMCALClusters; iClu++) {
+  //----------------------------------------------------------
+  // First recalibrate and recalculate energy and position
+  Float_t pos[]={0,0,0};
+  if(fCorrectClusters){
+    for(Int_t iClu=0; iClu<kNumberOfEMCALClusters-1; iClu++) {
+      
+      AliVCluster *c1 = (AliVCluster *) caloClustersArr->At(iClu);
+      
+      if(fRecoUtils->ClusterContainsBadChannel(fEMCALGeo, c1->GetCellsAbsId(), c1->GetNCells())) continue;     
+      
+      if(DebugLevel() > 2)
+      { 
+        printf("Std  : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",c1->GetID(),c1->E(),c1->GetDispersion(),c1->GetM02(),c1->GetM20());
+        c1->GetPosition(pos);
+        printf("Std  : i %d, x %f, y %f, z %f\n",c1->GetID(), pos[0], pos[1], pos[2]);
+      }
+      
+      //Correct cluster energy and position if requested, and not corrected previously, by default Off
+      if(fRecoUtils->IsRecalibrationOn())      fRecoUtils->RecalibrateClusterEnergy(fEMCALGeo, c1, emCells);
+      if(DebugLevel() > 2) 
+        printf("Energy: after recalibration %f; ",c1->E());
+      
+      // Correct Non-Linearity
+      c1->SetE(fRecoUtils->CorrectClusterEnergyLinearity(c1));
+      if(DebugLevel() > 2) 
+        printf("after linearity correction %f\n",c1->E());
+      // Recalculate cluster position
+      fRecoUtils->RecalculateClusterPosition(fEMCALGeo, emCells,c1);
+      if(DebugLevel() > 2)
+      { 
+        printf("Cor  : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",c1->GetID(),c1->E(),c1->GetDispersion(),c1->GetM02(),c1->GetM20());
+        c1->GetPosition(pos);
+        printf("Cor  : i %d, x %f, y %f, z %f\n",c1->GetID(), pos[0], pos[1], pos[2]);
+      }    
+    }    
+  }
+  //----------------------------------------------------------
+  //Now the invariant mass analysis with the corrected clusters  
+  for(Int_t iClu=0; iClu<kNumberOfEMCALClusters-1; iClu++) {
     
-    AliAODCaloCluster *c1 = (AliAODCaloCluster *) caloClustersArr->At(iClu);
-    if(!fCopyAOD && ClusterContainsBadChannel(c1->GetCellsAbsId(), c1->GetNCells())) continue; 
+    AliVCluster *c1 = (AliVCluster *) caloClustersArr->At(iClu);
+    if(fRecoUtils->ClusterContainsBadChannel(fEMCALGeo, c1->GetCellsAbsId(), c1->GetNCells())) continue;       
     
     Float_t e1i = c1->E();   // cluster energy before correction   
     if(e1i < fEmin) continue;
@@ -716,74 +466,43 @@ void AliAnalysisTaskEMCALPi0CalibSelection::UserExec(Option_t* /* option */)
     
     if(DebugLevel() > 2)
     { 
-      printf("Std  : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",iClu,e1i, c1->GetDispersion(),c1->GetM02(),c1->GetM20());
-      Float_t pos[]={0,0,0};
+      printf("IMA  : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",c1->GetID(),e1i,c1->GetDispersion(),c1->GetM02(),c1->GetM20());
       c1->GetPosition(pos);
-      printf("Std  : i %d, x %f, y %f, z %f\n",iClu, pos[0], pos[1], pos[2]);
+      printf("IMA  : i %d, x %f, y %f, z %f\n",c1->GetID(), pos[0], pos[1], pos[2]);
     }
     
-    //AliEMCALAodCluster clu1(*c1);
+    //AliEMCALAodCluster newc1(*((AliAODCaloCluster*)c1));
+    //newc1.EvalAllFromRecoUtils(fEMCALGeo,fRecoUtils,emCells);
+    //printf("i %d, recal? %d\n",iClu,newc1.IsRecalibrated());
     //clu1.Recalibrate(fCalibData, emCells, fEMCALGeoName);
     //clu1.EvalEnergy();
     //clu1.EvalAll(fLogWeight, fEMCALGeoName);
-    if(fRecoUtils->IsRecalibrationOn())        fRecoUtils->RecalibrateClusterEnergy(fEMCALGeo, c1, emCells);
-    
-    //Float_t e1ii = clu1.E(); // cluster energy after correction
-    Float_t e1ii = c1->E(); // cluster energy after correction
-        
-    if(DebugLevel() > 2)
-      printf("Recal: i %d, E %f, dispersion %f, m02 %f, m20 %f\n",iClu,e1ii, c1->GetDispersion(),c1->GetM02(),c1->GetM20());    
-    
-    
-    //clu1.GetMomentum(p1,v);
     
-    // Correct Non-Linearity
-    c1->SetE(fRecoUtils->CorrectClusterEnergyLinearity(c1));
-    //printf("\t  e1 org %2.3f, e1 cor  %2.3f \n",e1ii,c1->E());
-
-    //c1->GetMomentum(p1,v);
-    //printf("\t cor: e %2.3f, pt %2.3f, px %2.3f, py %2.3f, pz %2.3f\n",p1.E(), p1.Pt(),p1.Px(),p1.Py(),p1.Pz());
-
-    //Get tower with maximum energy and fill in the end the pi0 histogram for this cell, recalculate cluster position and recalibrate    
-    //Float_t pos[3];
-    //c1->GetPosition(pos);
-    //printf("Before Alignment: e %2.4f, x %2.4f, y %2.4f , z %2.4f\n",c1->E(),pos[0], pos[1],pos[2]);
-    fRecoUtils->RecalculateClusterPosition(fEMCALGeo, emCells,c1);
-    fRecoUtils->GetMaxEnergyCell          (fEMCALGeo, emCells,c1,absId1,iSupMod1,ieta1,iphi1);
-
-    //printf("i1 %d, corr1 %2.3f, e1 %2.3f, , ecorr1 %2.3f, pt %2.3f, px %2.3f, py %2.3f, pz %2.3f,\n",iClu, 1./corrFac, e1, p1.E(), p1.Pt(),p1.Px(),p1.Py(),p1.Pz());
-    //c1->GetPosition(pos);
-    //printf("After Alignment: e %2.4f, x %2.4f, y %2.4f , z %2.4f\n",c1->E(),pos[0], pos[1],pos[2]);
-
+    fRecoUtils->GetMaxEnergyCell(fEMCALGeo, emCells,c1,absId1,iSupMod1,ieta1,iphi1);
     c1->GetMomentum(p1,v);
-
-    for (Int_t jClu=iClu; jClu<kNumberOfEMCALClusters; jClu++) {
+    //newc1.GetMomentum(p1,v);
+    
+    // Combine cluster with other clusters and get the invariant mass
+    for (Int_t jClu=iClu+1; jClu<kNumberOfEMCALClusters; jClu++) {
       AliAODCaloCluster *c2 = (AliAODCaloCluster *) caloClustersArr->At(jClu);
-      if(c2->IsEqual(c1)) continue;
-      if(!fCopyAOD && ClusterContainsBadChannel(c2->GetCellsAbsId(), c2->GetNCells())) continue;       
+      //if(c2->IsEqual(c1)) continue;
+      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 (c2->GetNCells() < fMinNCells) continue; 
       
-      //AliEMCALAodCluster clu2(*c2);
+      //AliEMCALAodCluster newc2(*((AliAODCaloCluster*)c2));
+      //newc2.EvalAllFromRecoUtils(fEMCALGeo,fRecoUtils,emCells);
+      //printf("\t j %d, recal? %d\n",jClu,newc2.IsRecalibrated());
       //clu2.Recalibrate(fCalibData, emCells,fEMCALGeoName);
       //clu2.EvalEnergy();
       //clu2.EvalAll(fLogWeight,fEMCALGeoName);
-      if(fRecoUtils->IsRecalibrationOn())      fRecoUtils->RecalibrateClusterEnergy(fEMCALGeo, c2, emCells);
-
-      Float_t e2ii = c2->E();
-      
-      //Correct Non-Linearity
-      c2->SetE(fRecoUtils->CorrectClusterEnergyLinearity(c2));
-
-      //Get tower with maximum energy and fill in the end the pi0 histogram for this cell, recalculate cluster position and recalibrate    
-      fRecoUtils->RecalculateClusterPosition(fEMCALGeo, emCells,c2);
-      fRecoUtils->GetMaxEnergyCell          (fEMCALGeo, emCells,c2,absId2,iSupMod2,ieta2,iphi2);
       
+      fRecoUtils->GetMaxEnergyCell(fEMCALGeo, emCells,c2,absId2,iSupMod2,ieta2,iphi2);
       c2->GetMomentum(p2,v);
-
+      //newc2.GetMomentum(p2,v);
       p12 = p1+p2;
       Float_t invmass = p12.M()*1000; 
       //printf("*** mass %f\n",invmass);
@@ -795,16 +514,16 @@ void AliAnalysisTaskEMCALPi0CalibSelection::UserExec(Option_t* /* option */)
       if(invmass < fMaxBin && invmass > fMinBin){
         
         //Check if cluster is in fidutial region, not too close to borders
-        Bool_t in1 = CheckCellFiducialRegion(c1, aod->GetEMCALCells());
-        Bool_t in2 = CheckCellFiducialRegion(c2, aod->GetEMCALCells());
-
+        Bool_t in1 = fRecoUtils->CheckCellFiducialRegion(fEMCALGeo, c1, emCells);
+        Bool_t in2 = fRecoUtils->CheckCellFiducialRegion(fEMCALGeo, c2, emCells);
+        
         if(in1 && in2){
           
           fHmgg->Fill(invmass,p12.Pt()); 
-        
+          
           if(iSupMod1==iSupMod2) fHmggSM[iSupMod1]->Fill(invmass,p12.Pt()); 
           else                   fHmggDifferentSM ->Fill(invmass,p12.Pt());
-        
+          
           if((iSupMod1==0 && iSupMod2==2) || (iSupMod1==2 && iSupMod2==0)) fHmggPairSM[0]->Fill(invmass,p12.Pt()); 
           if((iSupMod1==1 && iSupMod2==3) || (iSupMod1==3 && iSupMod2==1)) fHmggPairSM[1]->Fill(invmass,p12.Pt()); 
           if((iSupMod1==0 && iSupMod2==1) || (iSupMod1==1 && iSupMod2==0)) fHmggPairSM[2]->Fill(invmass,p12.Pt()); 
@@ -815,10 +534,8 @@ void AliAnalysisTaskEMCALPi0CalibSelection::UserExec(Option_t* /* option */)
             //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);
-
+            
             //Inciden angle of each photon
-            //Float_t * posSM1cen = RecalculatePosition(11.5, 23.5, p1.E(),0, iSupMod1); 
-            //Float_t * posSM2cen = RecalculatePosition(11.5, 23.5, p2.E(),0, iSupMod2); 
             Float_t posSM1cen[3]={0.,0.,0.};
             Float_t depth = fRecoUtils->GetDepth(p1.Energy(),fRecoUtils->GetParticleType(),iSupMod1);
             fEMCALGeo->RecalculateTowerPosition(11.5, 23.5, iSupMod1, depth, fRecoUtils->GetMisalTransShiftArray(),fRecoUtils->GetMisalRotShiftArray(),posSM1cen); 
@@ -838,7 +555,7 @@ void AliAnalysisTaskEMCALPi0CalibSelection::UserExec(Option_t* /* option */)
             fHIncidentAngle->Fill(inangle1,p1.Pt()); 
             fHIncidentAngle->Fill(inangle2,p2.Pt()); 
             fHAsymmetry    ->Fill(asym,p12.Pt()); 
-
+            
             if(iSupMod1==iSupMod2) {
               fHOpeningAngleSM[iSupMod1] ->Fill(opangle,p12.Pt());
               fHIncidentAngleSM[iSupMod1]->Fill(inangle1,p1.Pt());
@@ -857,14 +574,14 @@ void AliAnalysisTaskEMCALPi0CalibSelection::UserExec(Option_t* /* option */)
               fHIncidentAnglePairSM[0]->Fill(inangle1,p1.Pt());
               fHIncidentAnglePairSM[0]->Fill(inangle2,p2.Pt());
               fHAsymmetryPairSM[0]    ->Fill(asym,p12.Pt());
-
+              
             } 
             if((iSupMod1==1 && iSupMod2==3) || (iSupMod1==3 && iSupMod2==1)) {
               fHOpeningAnglePairSM[1] ->Fill(opangle,p12.Pt()); 
               fHIncidentAnglePairSM[1]->Fill(inangle1,p1.Pt());
               fHIncidentAnglePairSM[1]->Fill(inangle2,p2.Pt());
               fHAsymmetryPairSM[1]    ->Fill(asym,p12.Pt());
-
+              
             }
             
             if((iSupMod1==0 && iSupMod2==1) || (iSupMod1==1 && iSupMod2==0)) {
@@ -872,8 +589,8 @@ void AliAnalysisTaskEMCALPi0CalibSelection::UserExec(Option_t* /* option */)
               fHIncidentAnglePairSM[2]->Fill(inangle1,p1.Pt());
               fHIncidentAnglePairSM[2]->Fill(inangle2,p2.Pt());
               fHAsymmetryPairSM[2]    ->Fill(asym,p12.Pt());
-
-
+              
+              
             }
             if((iSupMod1==2 && iSupMod2==3) || (iSupMod1==3 && iSupMod2==2)) {
               fHOpeningAnglePairSM[3] ->Fill(opangle,p12.Pt()); 
@@ -881,28 +598,28 @@ void AliAnalysisTaskEMCALPi0CalibSelection::UserExec(Option_t* /* option */)
               fHIncidentAnglePairSM[3]->Fill(inangle2,p2.Pt());
               fHAsymmetryPairSM[3]    ->Fill(asym,p12.Pt());
             }
-              
+            
           }// pair in 100 < mass < 160
-        
+          
         }//in acceptance cuts
         
         //In case of filling only channels with second cluster in same SM
         if(fSameSM && iSupMod1!=iSupMod2) continue;
         
         if (fGroupNCells == 0){
-            fHmpi0[iSupMod1][ieta1][iphi1]->Fill(invmass);
-            fHmpi0[iSupMod2][ieta2][iphi2]->Fill(invmass);
+          fHmpi0[iSupMod1][ieta1][iphi1]->Fill(invmass);
+          fHmpi0[iSupMod2][ieta2][iphi2]->Fill(invmass);
           
-            if(invmass > 100. && invmass < 160.){//restrict to clusters really close to pi0 peak
-              fhTowerDecayPhotonHit      [iSupMod1]->Fill(ieta1,iphi1);
-              fhTowerDecayPhotonEnergy   [iSupMod1]->Fill(ieta1,iphi1,p1.E());
-              fhTowerDecayPhotonAsymmetry[iSupMod1]->Fill(ieta1,iphi1,asym);
-              
-              fhTowerDecayPhotonHit      [iSupMod2]->Fill(ieta2,iphi2);
-              fhTowerDecayPhotonEnergy   [iSupMod2]->Fill(ieta2,iphi2,p2.E());
-              fhTowerDecayPhotonAsymmetry[iSupMod2]->Fill(ieta2,iphi2,asym);
-
-            }// pair in mass of pi0
+          if(invmass > 100. && invmass < 160.){//restrict to clusters really close to pi0 peak
+            fhTowerDecayPhotonHit      [iSupMod1]->Fill(ieta1,iphi1);
+            fhTowerDecayPhotonEnergy   [iSupMod1]->Fill(ieta1,iphi1,p1.E());
+            fhTowerDecayPhotonAsymmetry[iSupMod1]->Fill(ieta1,iphi1,asym);
+            
+            fhTowerDecayPhotonHit      [iSupMod2]->Fill(ieta2,iphi2);
+            fhTowerDecayPhotonEnergy   [iSupMod2]->Fill(ieta2,iphi2,p2.E());
+            fhTowerDecayPhotonAsymmetry[iSupMod2]->Fill(ieta2,iphi2,asym);
+            
+          }// pair in mass of pi0
         }      
         else  {
           //printf("Regroup N %d, eta1 %d, phi1 %d, eta2 %d, phi2 %d \n",fGroupNCells, ieta1, iphi1, ieta2, iphi2);
@@ -922,7 +639,7 @@ void AliAnalysisTaskEMCALPi0CalibSelection::UserExec(Option_t* /* option */)
         }//group cells
         
         if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Mass in (SM%d,%d,%d) and  (SM%d,%d,%d): %.3f GeV  E1_i=%f E1_ii=%f  E2_i=%f E2_ii=%f\n",
-                                    iSupMod1,iphi1,ieta1,iSupMod2,iphi2,ieta2,p12.M(),e1i,e1ii,e2i,e2ii);
+                                    iSupMod1,iphi1,ieta1,iSupMod2,iphi2,ieta2,p12.M(),e1i,c1->E(),e2i,c2->E());
       }
       
     }
@@ -935,76 +652,6 @@ void AliAnalysisTaskEMCALPi0CalibSelection::UserExec(Option_t* /* option */)
   
 }
 
-
-//_______________________________________________________________
-Bool_t AliAnalysisTaskEMCALPi0CalibSelection::CheckCellFiducialRegion(AliVCluster* cluster, AliVCaloCells* cells) 
-{
-  
-       // Given the list of AbsId of the cluster, get the maximum cell and 
-       // check if there are fNCellsFromBorder from the calorimeter border
-       
-  //If the distance to the border is 0 or negative just exit accept all clusters
-       if(cells->GetType()==AliVCaloCells::kEMCALCell && fNCellsFromEMCALBorder <= 0 ) return kTRUE;
-  
-  Int_t absIdMax       = -1;
-       Float_t ampMax  = -1;
-  
-  for(Int_t i = 0; i < cluster->GetNCells() ; i++){
-    Int_t absId = cluster->GetCellAbsId(i) ;
-    Float_t amp        = cells->GetCellAmplitude(absId);
-    if(amp > ampMax){
-      ampMax   = amp;
-      absIdMax = absId;
-    }
-  }
-       
-       if(DebugLevel() > 1)
-               printf("AliAnalysisTaskEMCALPi0CalibSelection::CheckCellFiducialRegion() - Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f\n", 
-           absIdMax, ampMax, cluster->E());
-       
-       if(absIdMax==-1) return kFALSE;
-       
-       //Check if the cell is close to the borders:
-       Bool_t okrow = kFALSE;
-       Bool_t okcol = kFALSE;
-  
-  Int_t iTower = -1, iIphi = -1, iIeta = -1, iphi = -1, ieta = -1, iSM = -1; 
-  fEMCALGeo->GetCellIndex(absIdMax,iSM,iTower,iIphi,iIeta); 
-  fEMCALGeo->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
-  if(iSM < 0 || iphi < 0 || ieta < 0 ) {
-    Fatal("CheckCellFidutialRegion","Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",iSM,ieta,iphi);
-  }
-  
-  //Check rows/phi
-  if(iSM < 10){
-    if(iphi >= fNCellsFromEMCALBorder && iphi < 24-fNCellsFromEMCALBorder) okrow =kTRUE; 
-  }
-  else{
-    if(iphi >= fNCellsFromEMCALBorder && iphi < 12-fNCellsFromEMCALBorder) okrow =kTRUE; 
-  }
-  
-  //Check collumns/eta
-  if(iSM%2==0){
-    if(ieta >= fNCellsFromEMCALBorder)     okcol = kTRUE;      
-  }
-  else {
-    if(ieta <  48-fNCellsFromEMCALBorder)  okcol = kTRUE;      
-  }
-  
-  if(DebugLevel() > 1)
-  {
-    printf("AliAnalysisTaskEMCALPi0CalibSelection::CheckCellFiducialRegion() - EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d ?",
-           fNCellsFromEMCALBorder, ieta, iphi, iSM);
-    if (okcol && okrow ) printf(" YES \n");
-    else  printf(" NO: column ok? %d, row ok? %d \n",okcol,okrow);
-  }
-       
-       if (okcol && okrow) return kTRUE; 
-       else                return kFALSE;
-       
-}      
-
-
 //__________________________________________________
 //void AliAnalysisTaskEMCALPi0CalibSelection::SetCalibCorrections(AliEMCALCalibData* const cdata)
 //{
@@ -1014,55 +661,3 @@ Bool_t AliAnalysisTaskEMCALPi0CalibSelection::CheckCellFiducialRegion(AliVCluste
 //   fCalibData = cdata;
 //     
 //}
-
-//________________________________________________________________
-void AliAnalysisTaskEMCALPi0CalibSelection::InitEMCALBadChannelStatusMap(){
-       //Init EMCAL bad channels map
-       if(DebugLevel() > 0 )printf("AliAnalysisTaskEMCALPi0CalibSelection::InitEMCALBadChannelStatusMap()\n");
-       //In order to avoid rewriting the same histograms
-       Bool_t oldStatus = TH1::AddDirectoryStatus();
-       TH1::AddDirectory(kFALSE);
-       
-       fEMCALBadChannelMap = new TObjArray(12);
-       //TH2F * hTemp = new  TH2I("EMCALBadChannelMap","EMCAL SuperModule bad channel map", 48, 0, 48, 24, 0, 24);
-       for (int i = 0; i < 12; i++) {
-               fEMCALBadChannelMap->Add(new TH2I(Form("EMCALBadChannelMap_Mod%d",i),Form("EMCALBadChannelMap_Mod%d",i), 48, 0, 48, 24, 0, 24));
-               //fEMCALBadChannelMap->Add((TH2I*)hTemp->Clone(Form("EMCALBadChannelMap_Mod%d",i)));
-       }
-       
-       //delete hTemp;
-       
-       fEMCALBadChannelMap->SetOwner(kTRUE);
-       fEMCALBadChannelMap->Compress();
-       
-       //In order to avoid rewriting the same histograms
-       TH1::AddDirectory(oldStatus);           
-}
-
-
-//_________________________________________________________________________________________________________
-Bool_t AliAnalysisTaskEMCALPi0CalibSelection::ClusterContainsBadChannel(UShort_t* cellList, Int_t nCells){
-       // Check that in the cluster cells, there is no bad channel of those stored 
-       // in fEMCALBadChannelMap or fPHOSBadChannelMap
-       
-       if(!fRemoveBadChannels)  return kFALSE;
-       if(!fEMCALBadChannelMap) return kFALSE;
-       
-       Int_t icol = -1;
-       Int_t irow = -1;
-       Int_t imod = -1;
-       for(Int_t iCell = 0; iCell<nCells; iCell++){
-               
-               //Get the column and row
-                       Int_t iTower = -1, iIphi = -1, iIeta = -1; 
-                       fEMCALGeo->GetCellIndex(cellList[iCell],imod,iTower,iIphi,iIeta); 
-                       if(fEMCALBadChannelMap->GetEntries() <= imod) continue;
-                       fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);                     
-                       if(GetEMCALChannelStatus(imod, icol, irow))return kTRUE;
-               
-       }// cell cluster loop
-       
-       return kFALSE;
-       
-}
-
index 87bbf4d..9cdcb60 100644 (file)
@@ -17,8 +17,6 @@ class TH1F;
 // AliRoot includes
 #include "AliAnalysisTaskSE.h"
 class AliEMCALGeometry;
-class AliAODCaloCluster;
-class AliAODCaloCells;
 //class AliEMCALCalibData ;
 #include "AliEMCALGeoParams.h"
 class AliEMCALRecoUtils;
@@ -50,12 +48,10 @@ public:
   void SetLogWeight(Float_t w)           {fLogWeight   = w   ;}
   
   //void SetCalibCorrections(AliEMCALCalibData* const cdata);
-  void CreateAODFromESD();
-  void CreateAODFromAOD();     
-
-  void CopyAOD(Bool_t copy)  { fCopyAOD = copy ; }
-  Bool_t IsAODCopied() const { return fCopyAOD ; }
        
+  void SwitchOnClusterCorrection()    {fCorrectClusters = kTRUE  ; }
+  void SwitchOffClusterCorrection()   {fCorrectClusters = kFALSE ; }
+  
   void SwitchOnSameSM()    {fSameSM = kTRUE  ; }
   void SwitchOffSameSM()   {fSameSM = kFALSE ; }
   
@@ -66,32 +62,7 @@ public:
   
   void SetGeometryName(TString name) { fEMCALGeoName = name ; }
   TString GeometryName() const       { return fEMCALGeoName ; }
-  //Modules fiducial region
-  Bool_t CheckCellFiducialRegion(AliVCluster* cluster, AliVCaloCells* cells) ;
-  void   SetNumberOfCellsFromEMCALBorder(Int_t n) {fNCellsFromEMCALBorder = n; }
-  Int_t  GetNumberOfCellsFromEMCALBorder() const  {return fNCellsFromEMCALBorder; }
-  
-  // Bad channels, copy from PWG4/PartCorrBase/AliCalorimeterUtils
-  Bool_t IsBadChannelsRemovalSwitchedOn()  const { return fRemoveBadChannels ; }
-  void SwitchOnBadChannelsRemoval ()  {fRemoveBadChannels = kTRUE  ; InitEMCALBadChannelStatusMap();}
-  void SwitchOffBadChannelsRemoval()  {fRemoveBadChannels = kFALSE ; }
-       
-  void InitEMCALBadChannelStatusMap() ;
-       
-  Int_t GetEMCALChannelStatus(Int_t iSM , Int_t iCol, Int_t iRow) const { 
-       if(fEMCALBadChannelMap) return (Int_t) ((TH2I*)fEMCALBadChannelMap->At(iSM))->GetBinContent(iCol,iRow); 
-       else return 0;}//Channel is ok by default
-       
-  void SetEMCALChannelStatus(Int_t iSM , Int_t iCol, Int_t iRow, Double_t c = 1) { 
-       if(!fEMCALBadChannelMap)InitEMCALBadChannelStatusMap() ;
-       ((TH2I*)fEMCALBadChannelMap->At(iSM))->SetBinContent(iCol,iRow,c);}
-       
-  TH2I * GetEMCALChannelStatusMap(Int_t iSM) const {return (TH2I*)fEMCALBadChannelMap->At(iSM);}
-  void   SetEMCALChannelStatusMap(TObjArray *map)  {fEMCALBadChannelMap = map;}
-       
-  Bool_t ClusterContainsBadChannel(UShort_t* cellList, Int_t nCells);
-                 
+                 
   void SetEMCALRecoUtils(AliEMCALRecoUtils * ru) {fRecoUtils = ru;}
   AliEMCALRecoUtils* GetEMCALRecoUtils() const   {return fRecoUtils;}
   
@@ -100,25 +71,26 @@ public:
          
   void GetMaxEnergyCellPosAndClusterPos(AliVCaloCells* cells, AliVCluster* clu, Int_t& iSM, Int_t& ieta, Int_t& iphi);
 
+  void UseFilteredEventAsInput() {fFilteredInput = kTRUE;}
+  void UseNormalEventAsInput()   {fFilteredInput = kFALSE;}
+
 private:
 
   AliEMCALGeometry * fEMCALGeo;  //! EMCAL geometry
   //AliEMCALCalibData* fCalibData; // corrections to CC from the previous iteration
        
-  Float_t fEmin;          // min. cluster energy
-  Float_t fEmax;          // max. cluster energy
-  Float_t fAsyCut;        // Asymmetry cut
-  Int_t   fMinNCells;     // min. ncells in cluster
-  Int_t   fGroupNCells;   // group n cells
-  Float_t fLogWeight;     // log weight used in cluster recalibration
-  Bool_t  fCopyAOD;       // Copy calo information only to AOD?
-  Bool_t  fSameSM;        // Combine clusters in channels on same SM
-  Bool_t  fOldAOD;        // Reading Old AODs, created before release 4.20
-  TString fEMCALGeoName;  // Name of geometry to use.
-  Int_t   fNCellsFromEMCALBorder; //  Number of cells from EMCAL border the cell with maximum amplitude has to be.
-
-  Bool_t     fRemoveBadChannels;         // Check the channel status provided and remove clusters with bad channels
-  TObjArray *fEMCALBadChannelMap;        // Array of histograms with map of bad channels, EMCAL
+  Float_t fEmin;           // min. cluster energy
+  Float_t fEmax;           // max. cluster energy
+  Float_t fAsyCut;         // Asymmetry cut
+  Int_t   fMinNCells;      // min. ncells in cluster
+  Int_t   fGroupNCells;    // group n cells
+  Float_t fLogWeight;      // log weight used in cluster recalibration
+  Bool_t  fSameSM;         // Combine clusters in channels on same SM
+  Bool_t  fOldAOD;         // Reading Old AODs, created before release 4.20
+  Bool_t  fFilteredInput;  // Read input produced with filter.
+  Bool_t  fCorrectClusters;// Correct clusters energy, position etc.
+  TString fEMCALGeoName;   // Name of geometry to use.
+
   AliEMCALRecoUtils * fRecoUtils;  // Access to reconstruction utilities
   
   //Output histograms  
@@ -156,7 +128,7 @@ private:
   TH1I*   fhNEvents;        //! Number of events counter histogram
   TList * fCuts ;           //! List with analysis cuts
 
-  ClassDef(AliAnalysisTaskEMCALPi0CalibSelection,9);
+  ClassDef(AliAnalysisTaskEMCALPi0CalibSelection,10);
 
 };
 
index f44f4c6..0daa9e3 100644 (file)
@@ -27,9 +27,75 @@ Int_t kFile = 1; // Number of files
 char * kXML = "collection.xml";
 //---------------------------------------------------------------------------
 
-const TString kInputData = "AOD"; //ESD, AOD, MC
+const TString kInputData = "ESD"; //ESD, AOD, MC
 TString kTreeName = "esdTree";
-Bool_t copy = kFALSE;
+Bool_t copy = kTRUE;
+
+
+void SetRecoUtilsParams(AliEMCALRecoUtils* reco){
+
+    reco->SetParticleType(AliEMCALRecoUtils::kPhoton);
+    reco->SetW0(4.);
+
+    reco->SetPositionAlgorithm(AliEMCALRecoUtils::kPosTowerGlobal);
+    reco->SetMisalTransShift(0,1.134);   reco->SetMisalTransShift(1,8.2); reco->SetMisalTransShift(2,1.197);
+    reco->SetMisalTransShift(3,-3.093);  reco->SetMisalTransShift(5,6.82);reco->SetMisalTransShift(5,1.635);
+
+    //reco->SetPositionAlgorithm(AliEMCALRecoUtils::kPosTowerIndex);
+    //reco->SetMisalTransShift(0,1.08);   reco->SetMisalTransShift(1,8.35); reco->SetMisalTransShift(2,1.12);
+    //reco->SetMisalRotShift(3,-8.05);    reco->SetMisalRotShift(4,8.05);  
+    //reco->SetMisalTransShift(3,-0.42);  reco->SetMisalTransShift(5,1.55);
+
+    reco->SetNonLinearityFunction(AliEMCALRecoUtils::kPi0GammaGamma);
+    reco->SetNonLinearityParam(0,1.04);     reco->SetNonLinearityParam(1,-0.1445);
+    reco->SetNonLinearityParam(2,1.046);    
+
+    reco->SetNonLinearityFunction(AliEMCALRecoUtils::kPi0GammaConversion);
+    reco->SetNonLinearityParam(0,1.033);     reco->SetNonLinearityParam(1,0.0566186);
+    reco->SetNonLinearityParam(2,0.982133);    
+
+
+     reco->SetNonLinearityFunction(AliEMCALRecoUtils::kPi0MC);
+     reco->SetNonLinearityParam(0,1.001);     reco->SetNonLinearityParam(1,-0.01264);
+     reco->SetNonLinearityParam(2,-0.03632);    
+     reco->SetNonLinearityParam(3,0.1798);     reco->SetNonLinearityParam(4,-0.522);
+
+     reco->SwitchOnRecalibration();
+     TFile * f = new TFile("RecalibrationFactors.root","read");
+     TH2F * h0 = (TH2F*)f->Get("EMCALRecalFactors_SM0")->Clone();
+     TH2F * h1 = (TH2F*)f->Get("EMCALRecalFactors_SM1")->Clone();
+     TH2F * h2 = (TH2F*)f->Get("EMCALRecalFactors_SM2")->Clone();
+     TH2F * h3 = (TH2F*)f->Get("EMCALRecalFactors_SM3")->Clone();
+
+     reco->SetEMCALChannelRecalibrationFactors(0,h0);
+     reco->SetEMCALChannelRecalibrationFactors(1,h1);
+     reco->SetEMCALChannelRecalibrationFactors(2,h2);
+     reco->SetEMCALChannelRecalibrationFactors(3,h3);
+
+    reco->SwitchOnBadChannelsRemoval();
+    // SM0
+    reco->SetEMCALChannelStatus(0,3,13);  reco->SetEMCALChannelStatus(0,44,1); reco->SetEMCALChannelStatus(0,3,13); 
+    reco->SetEMCALChannelStatus(0,20,7);  reco->SetEMCALChannelStatus(0,38,2);   
+    // SM1
+    reco->SetEMCALChannelStatus(1,4,7);   reco->SetEMCALChannelStatus(1,4,13);  reco->SetEMCALChannelStatus(1,9,20); 
+    reco->SetEMCALChannelStatus(1,14,15); reco->SetEMCALChannelStatus(1,23,16); reco->SetEMCALChannelStatus(1,32,23); 
+    reco->SetEMCALChannelStatus(1,37,5);  reco->SetEMCALChannelStatus(1,40,1);  reco->SetEMCALChannelStatus(1,40,2);
+    reco->SetEMCALChannelStatus(1,40,5);  reco->SetEMCALChannelStatus(1,41,0);  reco->SetEMCALChannelStatus(1,41,1);
+    reco->SetEMCALChannelStatus(1,41,2);  reco->SetEMCALChannelStatus(1,41,4);
+    // SM2        
+    reco->SetEMCALChannelStatus(2,14,15); reco->SetEMCALChannelStatus(2,18,16); reco->SetEMCALChannelStatus(2,18,17); 
+    reco->SetEMCALChannelStatus(2,18,18); reco->SetEMCALChannelStatus(2,18,20); reco->SetEMCALChannelStatus(2,18,21); 
+    reco->SetEMCALChannelStatus(2,18,23); reco->SetEMCALChannelStatus(2,19,16); reco->SetEMCALChannelStatus(2,19,17); 
+    reco->SetEMCALChannelStatus(2,19,19); reco->SetEMCALChannelStatus(2,19,20); reco->SetEMCALChannelStatus(2,19,21); 
+    reco->SetEMCALChannelStatus(2,19,22);
+    //SM3
+    reco->SetEMCALChannelStatus(3,4,7);
+
+    reco->SetNumberOfCellsFromEMCALBorder(1);
+
+    //reco->Print("");
+
+}
 
 void anaEMCALCalib(Int_t mode=mLocal)
 {
@@ -46,7 +112,7 @@ void anaEMCALCalib(Int_t mode=mLocal)
   //gSystem->Unload("libPWG4CaloCalib.so");
   //Try to set the new library
   //gSystem->Load("./PWG4CaloCalib/libPWG4CaloCalib.so");
-  gSystem->ListLibraries();
+  //gSystem->ListLibraries();
 
   //-------------------------------------------------------------------------------------------------
   //Create chain from ESD and from cross sections files, look below for options.
@@ -109,98 +175,66 @@ void anaEMCALCalib(Int_t mode=mLocal)
 
       gROOT->LoadMacro("AddTaskPhysicsSelection.C");
       AliPhysicsSelectionTask* physSelTask = AddTaskPhysicsSelection();
-      if(!copy){
-       gROOT->LoadMacro("AddTaskESDFilter.C");
-       AliAnalysisTaskESDfilter *esdfilter = AddTaskESDFilter(kFALSE);
-      }
+//       if(!copy){
+//     gROOT->LoadMacro("AddTaskESDFilter.C");
+//     AliAnalysisTaskESDfilter *esdfilter = AddTaskESDFilter(kFALSE);
+//       }
     }
     
+    // Create containers for input/output
+    AliAnalysisDataContainer *cinput1  = mgr->GetCommonInputContainer();
+    if(copy) {
 
-   
-    AliAnalysisTaskEMCALPi0CalibSelection * pi0calib = new AliAnalysisTaskEMCALPi0CalibSelection ("EMCALPi0Calibration");
-    //pi0calib->SetDebugLevel(10); 
-    pi0calib->CopyAOD(copy);
-    pi0calib->SetClusterMinEnergy(0.5);
-    pi0calib->SetClusterMaxEnergy(15.);
-    //pi0calib->SetAsymmetryCut(0.7);
-    pi0calib->SetClusterMinNCells(0);
-    pi0calib->SetNCellsGroup(0);
-    pi0calib->SwitchOnBadChannelsRemoval();
-    pi0calib->SwitchOffSameSM();
-    pi0calib->SwitchOnOldAODs();
-    pi0calib->SetNumberOfCellsFromEMCALBorder(1);
-    AliEMCALRecoUtils * reco = pi0calib->GetEMCALRecoUtils();
-
-    reco->SetParticleType(AliEMCALRecoUtils::kPhoton);
-    reco->SetW0(4.);
-
-    reco->SetPositionAlgorithm(AliEMCALRecoUtils::kPosTowerGlobal);
-    reco->SetMisalTransShift(0,1.134);   reco->SetMisalTransShift(1,8.2); reco->SetMisalTransShift(2,1.197);
-    reco->SetMisalTransShift(3,-3.093);  reco->SetMisalTransShift(5,6.82);reco->SetMisalTransShift(5,1.635);
-
-    //reco->SetPositionAlgorithm(AliEMCALRecoUtils::kPosTowerIndex);
-    //reco->SetMisalTransShift(0,1.08);   reco->SetMisalTransShift(1,8.35); reco->SetMisalTransShift(2,1.12);
-    //reco->SetMisalRotShift(3,-8.05);    reco->SetMisalRotShift(4,8.05);  
-    //reco->SetMisalTransShift(3,-0.42);  reco->SetMisalTransShift(5,1.55);
-
-    reco->SetNonLinearityFunction(AliEMCALRecoUtils::kPi0GammaGamma);
-    reco->SetNonLinearityParam(0,1.04);     reco->SetNonLinearityParam(1,-0.1445);
-    reco->SetNonLinearityParam(2,1.046);    
+      AliAnalysisDataContainer *coutput1 = mgr->GetCommonOutputContainer();
 
-//     reco->SetNonLinearityFunction(AliEMCALRecoUtils::kPi0GammaConversion);
-//     reco->SetNonLinearityParam(0,1.033);     reco->SetNonLinearityParam(1,0.0566186);
-//     reco->SetNonLinearityParam(2,0.982133);    
+      AliAnalysisTaskCaloFilter * filter = new AliAnalysisTaskCaloFilter("EMCALFilter");
+      //filter->SelectCollisionCandidates(); 
+      filter->SetCaloFilter(AliAnalysisTaskCaloFilter::kEMCAL);
+      filter->SwitchOnClusterCorrection();
 
+      AliEMCALRecoUtils * reco = filter->GetEMCALRecoUtils();
+      SetRecoUtilsParams(reco);
+      reco->Print("");
 
-//      reco->SetNonLinearityFunction(AliEMCALRecoUtils::kPi0MC);
-//      reco->SetNonLinearityParam(0,1.001);     reco->SetNonLinearityParam(1,-0.01264);
-//      reco->SetNonLinearityParam(2,-0.03632);    
-//      reco->SetNonLinearityParam(3,0.1798);     reco->SetNonLinearityParam(4,-0.522);
+      mgr->AddTask(filter);
 
-     reco->SwitchOnRecalibration();
-     TFile * f = new TFile("RecalibrationFactors.root","read");
-     TH2F * h0 = (TH2F*)f->Get("EMCALRecalFactors_SM0")->Clone();
-     TH2F * h1 = (TH2F*)f->Get("EMCALRecalFactors_SM1")->Clone();
-     TH2F * h2 = (TH2F*)f->Get("EMCALRecalFactors_SM2")->Clone();
-     TH2F * h3 = (TH2F*)f->Get("EMCALRecalFactors_SM3")->Clone();
-
-     reco->SetEMCALChannelRecalibrationFactors(0,h0);
-     reco->SetEMCALChannelRecalibrationFactors(1,h1);
-     reco->SetEMCALChannelRecalibrationFactors(2,h2);
-     reco->SetEMCALChannelRecalibrationFactors(3,h3);
-     reco->Print("");
+   
+    //AliAnalysisDataContainer *cout_cuts2 = mgr->CreateContainer("Cuts", TList::Class(), 
+    //                                        AliAnalysisManager::kOutputContainer, "pi0calib.root");
+    
+    mgr->ConnectInput  (filter,  0, cinput1);
+    mgr->ConnectOutput (filter, 0, coutput1 );
+    //mgr->ConnectOutput (filter, 2, cout_cuts2);
 
-    // SM0
-    pi0calib->SetEMCALChannelStatus(0,3,13);  pi0calib->SetEMCALChannelStatus(0,44,1); pi0calib->SetEMCALChannelStatus(0,3,13); 
-    pi0calib->SetEMCALChannelStatus(0,20,7);  pi0calib->SetEMCALChannelStatus(0,38,2);   
-    // SM1
-    pi0calib->SetEMCALChannelStatus(1,4,7);   pi0calib->SetEMCALChannelStatus(1,4,13);  pi0calib->SetEMCALChannelStatus(1,9,20); 
-    pi0calib->SetEMCALChannelStatus(1,14,15); pi0calib->SetEMCALChannelStatus(1,23,16); pi0calib->SetEMCALChannelStatus(1,32,23); 
-    pi0calib->SetEMCALChannelStatus(1,37,5);  pi0calib->SetEMCALChannelStatus(1,40,1);  pi0calib->SetEMCALChannelStatus(1,40,2);
-    pi0calib->SetEMCALChannelStatus(1,40,5);  pi0calib->SetEMCALChannelStatus(1,41,0);  pi0calib->SetEMCALChannelStatus(1,41,1);
-    pi0calib->SetEMCALChannelStatus(1,41,2);  pi0calib->SetEMCALChannelStatus(1,41,4);
-    // SM2        
-    pi0calib->SetEMCALChannelStatus(2,14,15); pi0calib->SetEMCALChannelStatus(2,18,16); pi0calib->SetEMCALChannelStatus(2,18,17); 
-    pi0calib->SetEMCALChannelStatus(2,18,18); pi0calib->SetEMCALChannelStatus(2,18,20); pi0calib->SetEMCALChannelStatus(2,18,21); 
-    pi0calib->SetEMCALChannelStatus(2,18,23); pi0calib->SetEMCALChannelStatus(2,19,16); pi0calib->SetEMCALChannelStatus(2,19,17); 
-    pi0calib->SetEMCALChannelStatus(2,19,19); pi0calib->SetEMCALChannelStatus(2,19,20); pi0calib->SetEMCALChannelStatus(2,19,21); 
-    pi0calib->SetEMCALChannelStatus(2,19,22);
-    //SM3
-    pi0calib->SetEMCALChannelStatus(3,4,7);
 
-    mgr->AddTask(pi0calib);
+    }
+     AliAnalysisTaskEMCALPi0CalibSelection * pi0calib = new AliAnalysisTaskEMCALPi0CalibSelection ("EMCALPi0Calibration");
+     //pi0calib->SelectCollisionCandidates(); 
+//     //pi0calib->SetDebugLevel(10); 
+     if(copy)
+       pi0calib->UseFilteredEventAsInput();
+     pi0calib->SetClusterMinEnergy(0.5);
+     pi0calib->SetClusterMaxEnergy(15.);
+     //pi0calib->SetAsymmetryCut(0.7);
+     pi0calib->SetClusterMinNCells(0);
+     pi0calib->SetNCellsGroup(0);
+     pi0calib->SwitchOnSameSM();
+     pi0calib->SwitchOffClusterCorrection();
+     //pi0calib->SwitchOnOldAODs();
+//     AliEMCALRecoUtils * reco = pi0calib->GetEMCALRecoUtils();
+//     SetRecoUtilsParams(reco);
+//     reco->Print("");
+
+     mgr->AddTask(pi0calib);
     
-    // Create containers for input/output
-    AliAnalysisDataContainer *cinput1  = mgr->GetCommonInputContainer();
-    if(copy) AliAnalysisDataContainer *coutput1 = mgr->GetCommonOutputContainer();
+
     AliAnalysisDataContainer *coutput2 = 
-      mgr->CreateContainer("pi0calib", TList::Class(), AliAnalysisManager::kOutputContainer, "pi0calib.root");
+    mgr->CreateContainer("pi0calib", TList::Class(), AliAnalysisManager::kOutputContainer, "pi0calib.root");
     
     AliAnalysisDataContainer *cout_cuts = mgr->CreateContainer("Cuts", TList::Class(), 
-                                                              AliAnalysisManager::kOutputContainer, "pi0calib.root");
+                                                      AliAnalysisManager::kOutputContainer, "pi0calib.root");
     
     mgr->ConnectInput  (pi0calib,     0, cinput1);
-    if(copy) mgr->ConnectOutput (pi0calib, 0, coutput1 );
     mgr->ConnectOutput (pi0calib, 1, coutput2 );
     mgr->ConnectOutput (pi0calib, 2, cout_cuts);
 
index 7bd86c9..21bea80 100755 (executable)
@@ -45,13 +45,11 @@ ClassImp(AliCalorimeterUtils)
     fEMCALGeoName("EMCAL_FIRSTYEAR"),fPHOSGeoName("PHOSgeo"), 
     fEMCALGeo(0x0), fPHOSGeo(0x0), 
     fEMCALGeoMatrixSet(kFALSE), fPHOSGeoMatrixSet(kFALSE), 
-    fRemoveBadChannels(kFALSE),
-    fEMCALBadChannelMap(0x0),fPHOSBadChannelMap(0x0), 
-    fNCellsFromEMCALBorder(0), fNCellsFromPHOSBorder(0), 
-    fNoEMCALBorderAtEta0(kFALSE),fRecalibration(kFALSE),
+    fRemoveBadChannels(kFALSE),fPHOSBadChannelMap(0x0), 
+    fNCellsFromPHOSBorder(0), fRecalibration(kFALSE),
     fPHOSRecalibrationFactors(),
-    fEMCALRecoUtils(new AliEMCALRecoUtils),fRecalculatePosition(kFALSE),fCorrectELinearity(kFALSE)
-
+    fEMCALRecoUtils(new AliEMCALRecoUtils),
+    fRecalculatePosition(kFALSE),fCorrectELinearity(kFALSE)
 {
   //Ctor
   
@@ -65,11 +63,7 @@ AliCalorimeterUtils::~AliCalorimeterUtils() {
   
   //if(fPHOSGeo)  delete fPHOSGeo  ;
   if(fEMCALGeo) delete fEMCALGeo ;
-       
-  if(fEMCALBadChannelMap) { 
-    fEMCALBadChannelMap->Clear();
-    delete  fEMCALBadChannelMap;
-  }
+
   if(fPHOSBadChannelMap) { 
     fPHOSBadChannelMap->Clear();
     delete  fPHOSBadChannelMap;
@@ -91,7 +85,7 @@ Bool_t AliCalorimeterUtils::CheckCellFiducialRegion(AliVCluster* cluster, AliVCa
        // check if there are fNCellsFromBorder from the calorimeter border
        
     //If the distance to the border is 0 or negative just exit accept all clusters
-       if(cells->GetType()==AliVCaloCells::kEMCALCell && fNCellsFromEMCALBorder <= 0 ) return kTRUE;
+       if(cells->GetType()==AliVCaloCells::kEMCALCell && fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder() <= 0 ) return kTRUE;
        if(cells->GetType()==AliVCaloCells::kPHOSCell  && fNCellsFromPHOSBorder  <= 0 ) return kTRUE;
 
   Int_t absIdMax       = -1;
@@ -169,29 +163,30 @@ Bool_t AliCalorimeterUtils::CheckCellFiducialRegion(AliVCluster* cluster, AliVCa
     }
       
                //Check rows/phi
+    Int_t nborder = fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder();
                if(iSM < 10){
-                       if(iphi >= fNCellsFromEMCALBorder && iphi < 24-fNCellsFromEMCALBorder) okrow =kTRUE; 
+                       if(iphi >= nborder && iphi < 24-nborder) okrow =kTRUE; 
            }
                else{
-                       if(iphi >= fNCellsFromEMCALBorder && iphi < 12-fNCellsFromEMCALBorder) okrow =kTRUE; 
+                       if(iphi >= nborder && iphi < 12-nborder) okrow =kTRUE; 
                }
                
-               //Check collumns/eta
-               if(!fNoEMCALBorderAtEta0){
-                       if(ieta  > fNCellsFromEMCALBorder && ieta < 48-fNCellsFromEMCALBorder) okcol =kTRUE; 
+               //Check columns/eta
+               if(!fEMCALRecoUtils->IsEMCALNoBorderAtEta0()){
+                       if(ieta  > nborder && ieta < 48-nborder) okcol =kTRUE; 
                }
                else{
                        if(iSM%2==0){
-                               if(ieta >= fNCellsFromEMCALBorder)     okcol = kTRUE;   
+                               if(ieta >= nborder)     okcol = kTRUE;  
                        }
                        else {
-                               if(ieta <  48-fNCellsFromEMCALBorder)  okcol = kTRUE;   
+                               if(ieta <  48-nborder)  okcol = kTRUE;  
                        }
                }//eta 0 not checked
                if(fDebug > 1)
                {
                        printf("AliCalorimeterUtils::CheckCellFiducialRegion() - EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d ?",
-                                  fNCellsFromEMCALBorder, ieta, iphi, iSM);
+                                  nborder, ieta, iphi, iSM);
                        if (okcol && okrow ) printf(" YES \n");
                        else  printf(" NO: column ok? %d, row ok? %d \n",okcol,okrow);
                }
@@ -226,7 +221,7 @@ Bool_t AliCalorimeterUtils::ClusterContainsBadChannel(TString calorimeter,UShort
        
        if (!fRemoveBadChannels) return kFALSE;
        //printf("fEMCALBadChannelMap %p, fPHOSBadChannelMap %p \n",fEMCALBadChannelMap,fPHOSBadChannelMap);
-       if(calorimeter == "EMCAL" && !fEMCALBadChannelMap) return kFALSE;
+       if(calorimeter == "EMCAL" && !fEMCALRecoUtils->GetEMCALChannelStatusMap(0)) return kFALSE;
        if(calorimeter == "PHOS"  && !fPHOSBadChannelMap)  return kFALSE;
 
        Int_t icol = -1;
@@ -236,16 +231,7 @@ Bool_t AliCalorimeterUtils::ClusterContainsBadChannel(TString calorimeter,UShort
        
                //Get the column and row
                if(calorimeter == "EMCAL"){
-                       Int_t iTower = -1, iIphi = -1, iIeta = -1; 
-                       fEMCALGeo->GetCellIndex(cellList[iCell],imod,iTower,iIphi,iIeta); 
-                       if(fEMCALBadChannelMap->GetEntries() <= imod) continue;
-                       fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);             
-      
-      if(imod < 0 || irow < 0 || icol < 0 ) {
-        Fatal("ClusterContainsBadChannel","Negative value for super module: %d, or cell icol: %d, or cell irow: %d, check EMCAL geometry name\n",imod,icol,irow);
-      }
-      
-                       if(GetEMCALChannelStatus(imod, icol, irow))return kTRUE;
+      return fEMCALRecoUtils->ClusterContainsBadChannel((AliEMCALGeometry*)fEMCALGeo,cellList,nCells);
                }
                else if(calorimeter=="PHOS"){
                        Int_t    relId[4];
@@ -432,34 +418,9 @@ void AliCalorimeterUtils::InitParameters()
                
   fRemoveBadChannels   = kFALSE;
        
-  fNCellsFromEMCALBorder   = 0;
   fNCellsFromPHOSBorder    = 0;
-  fNoEMCALBorderAtEta0     = kFALSE;
 }
 
-//________________________________________________________________
-void AliCalorimeterUtils::InitEMCALBadChannelStatusMap(){
-  //Init EMCAL bad channels map
-   if(fDebug > 0 )printf("AliCalorimeterUtils::InitEMCALBadChannelStatusMap()\n");
-  //In order to avoid rewriting the same histograms
-  Bool_t oldStatus = TH1::AddDirectoryStatus();
-  TH1::AddDirectory(kFALSE);
-
-  fEMCALBadChannelMap = new TObjArray(12);
-  //TH2F * hTemp = new  TH2I("EMCALBadChannelMap","EMCAL SuperModule bad channel map", 48, 0, 48, 24, 0, 24);
-  for (int i = 0; i < 12; i++) {
-    fEMCALBadChannelMap->Add(new TH2I(Form("EMCALBadChannelMap_Mod%d",i),Form("EMCALBadChannelMap_Mod%d",i), 48, 0, 48, 24, 0, 24));
-    //fEMCALBadChannelMap->Add((TH2I*)hTemp->Clone(Form("EMCALBadChannelMap_Mod%d",i)));
-  }
-
-  //delete hTemp;
-
-  fEMCALBadChannelMap->SetOwner(kTRUE);
-  fEMCALBadChannelMap->Compress();
-
-  //In order to avoid rewriting the same histograms
-  TH1::AddDirectory(oldStatus);                
-}
 
 //________________________________________________________________
 void AliCalorimeterUtils::InitPHOSBadChannelStatusMap(){
@@ -544,8 +505,8 @@ void AliCalorimeterUtils::Print(const Option_t * opt) const
   printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
   printf("Remove Clusters with bad channels? %d\n",fRemoveBadChannels);
   printf("Remove Clusters with max cell at less than %d cells from EMCAL border and %d cells from PHOS border\n",
-        fNCellsFromEMCALBorder, fNCellsFromPHOSBorder);
-  if(fNoEMCALBorderAtEta0) printf("Do not remove EMCAL clusters at Eta = 0\n");
+        fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder(), fNCellsFromPHOSBorder);
+  if(fEMCALRecoUtils->IsEMCALNoBorderAtEta0()) printf("Do not remove EMCAL clusters at Eta = 0\n");
   printf("Recalibrate Clusters? %d\n",fRecalibration);
   printf("Recalculate Clusters Position? %d\n",fRecalculatePosition);
   printf("Recalculate Clusters Energy? %d\n",fCorrectELinearity);
index 71582ed..6a23242 100755 (executable)
@@ -64,32 +64,29 @@ class AliCalorimeterUtils : public TObject {
        
   // Bad channels
   Bool_t IsBadChannelsRemovalSwitchedOn()  const { return fRemoveBadChannels ; }
-  void SwitchOnBadChannelsRemoval ()  {fRemoveBadChannels = kTRUE  ; InitEMCALBadChannelStatusMap(); InitPHOSBadChannelStatusMap();}
-  void SwitchOffBadChannelsRemoval()  {fRemoveBadChannels = kFALSE ; }
+  void SwitchOnBadChannelsRemoval ()  {fRemoveBadChannels = kTRUE  ; fEMCALRecoUtils->SwitchOnBadChannelsRemoval(); InitPHOSBadChannelStatusMap();}
+  void SwitchOffBadChannelsRemoval()  {fRemoveBadChannels = kFALSE ; fEMCALRecoUtils->SwitchOffBadChannelsRemoval();}
        
-  void InitEMCALBadChannelStatusMap() ;
   void InitPHOSBadChannelStatusMap () ;
 
   Int_t GetEMCALChannelStatus(Int_t iSM , Int_t iCol, Int_t iRow) const { 
-    if(fEMCALBadChannelMap) return (Int_t) ((TH2I*)fEMCALBadChannelMap->At(iSM))->GetBinContent(iCol,iRow); 
-    else return 0;}//Channel is ok by default
+    return fEMCALRecoUtils->GetEMCALChannelStatus(iSM,iCol,iRow); }//Channel is ok by default
 
   Int_t GetPHOSChannelStatus (Int_t imod, Int_t iCol, Int_t iRow) const { 
     if(fPHOSBadChannelMap)return (Int_t) ((TH2I*)fPHOSBadChannelMap->At(imod))->GetBinContent(iCol,iRow); 
     else return 0;}//Channel is ok by default
   
   void SetEMCALChannelStatus(Int_t iSM , Int_t iCol, Int_t iRow, Double_t c = 1) { 
-    if(!fEMCALBadChannelMap)InitEMCALBadChannelStatusMap() ;
-    ((TH2I*)fEMCALBadChannelMap->At(iSM))->SetBinContent(iCol,iRow,c);}
+    fEMCALRecoUtils->SetEMCALChannelStatus(iSM,iCol,iRow,c);}
   
   void SetPHOSChannelStatus (Int_t imod, Int_t iCol, Int_t iRow, Double_t c = 1) {
        if(!fPHOSBadChannelMap) InitPHOSBadChannelStatusMap() ; 
        ((TH2I*)fPHOSBadChannelMap->At(imod))->SetBinContent(iCol,iRow,c);}
     
-  TH2I * GetEMCALChannelStatusMap(Int_t iSM) const {return (TH2I*)fEMCALBadChannelMap->At(iSM);}
+  TH2I * GetEMCALChannelStatusMap(Int_t iSM) const {return fEMCALRecoUtils->GetEMCALChannelStatusMap(iSM);}
   TH2I * GetPHOSChannelStatusMap(Int_t imod) const {return (TH2I*)fPHOSBadChannelMap->At(imod);}
 
-  void SetEMCALChannelStatusMap(TObjArray *map) {fEMCALBadChannelMap = map;}
+  void SetEMCALChannelStatusMap(TObjArray *map) {fEMCALRecoUtils->SetEMCALChannelStatusMap(map);}
   void SetPHOSChannelStatusMap (TObjArray *map) {fPHOSBadChannelMap  = map;}
        
   Bool_t ClusterContainsBadChannel(TString calorimeter,UShort_t* cellList, Int_t nCells);
@@ -101,15 +98,11 @@ class AliCalorimeterUtils : public TObject {
        
   //Modules fiducial region
   Bool_t CheckCellFiducialRegion(AliVCluster* cluster, AliVCaloCells* cells, AliVEvent * event, Int_t iev=0) const ;
-       
-  void   SetNumberOfCellsFromEMCALBorder(Int_t n) {fNCellsFromEMCALBorder = n; }
-  Int_t  GetNumberOfCellsFromEMCALBorder() const  {return fNCellsFromEMCALBorder; }
   void   SetNumberOfCellsFromPHOSBorder(Int_t n)  {fNCellsFromPHOSBorder = n; }
   Int_t  GetNumberOfCellsFromPHOSBorder() const   {return fNCellsFromPHOSBorder; }
-       
-  void   SwitchOnNoFiducialBorderInEMCALEta0()  {fNoEMCALBorderAtEta0 = kTRUE; }
-  void   SwitchOffNoFiducialBorderInEMCALEta0() {fNoEMCALBorderAtEta0 = kFALSE; }
-       
+  void   SetNumberOfCellsFromEMCALBorder(Int_t n) {fEMCALRecoUtils->SetNumberOfCellsFromEMCALBorder(n); }
+  Int_t  GetNumberOfCellsFromEMCALBorder() const  {return fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder(); }
+  
   // Recalibration
   Bool_t IsRecalibrationOn()  const { return fRecalibration ; }
   void SwitchOnRecalibration()    {fRecalibration = kTRUE ; InitPHOSRecalibrationFactors(); fEMCALRecoUtils->SwitchOnRecalibration();}
@@ -164,11 +157,8 @@ class AliCalorimeterUtils : public TObject {
   Bool_t             fEMCALGeoMatrixSet;     //  Check if the transformation matrix is set for EMCAL
   Bool_t             fPHOSGeoMatrixSet ;     //  Check if the transformation matrix is set for PHOS
   Bool_t             fRemoveBadChannels;     //  Check the channel status provided and remove clusters with bad channels
-  TObjArray        * fEMCALBadChannelMap;    //  Array of histograms with map of bad channels, EMCAL
   TObjArray        * fPHOSBadChannelMap;     //  Array of histograms with map of bad channels, PHOS
-  Int_t              fNCellsFromEMCALBorder; //  Number of cells from EMCAL border the cell with maximum amplitude has to be.
   Int_t              fNCellsFromPHOSBorder;  //  Number of cells from PHOS  border the cell with maximum amplitude has to be.
-  Bool_t             fNoEMCALBorderAtEta0;   //  Do fiducial cut in EMCAL region eta = 0?
   Bool_t             fRecalibration;         //  Switch on or off the recalibration
   TObjArray        * fPHOSRecalibrationFactors;  // Array of histograms with map of recalibration factors, PHOS
   AliEMCALRecoUtils* fEMCALRecoUtils;        //  EMCAL utils for cluster rereconstruction