Update in AliEMCALRecoUtils, some of its new functionalities moved from AliCalorimete...
authorgconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 22 Oct 2010 16:56:23 +0000 (16:56 +0000)
committergconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 22 Oct 2010 16:56:23 +0000 (16:56 +0000)
PWG4/CaloCalib/AliAnalysisTaskEMCALPi0CalibSelection.cxx
PWG4/CaloCalib/AliAnalysisTaskEMCALPi0CalibSelection.h
PWG4/CaloCalib/macros/anaEMCALCalib.C
PWG4/PartCorrBase/AliCaloTrackReader.cxx
PWG4/PartCorrBase/AliCalorimeterUtils.cxx
PWG4/PartCorrBase/AliCalorimeterUtils.h

index c459327..15cabfb 100644 (file)
@@ -53,7 +53,6 @@ AliAnalysisTaskEMCALPi0CalibSelection::AliAnalysisTaskEMCALPi0CalibSelection(con
   fLogWeight(4.5), fCopyAOD(kFALSE), fSameSM(kFALSE), fOldAOD(kFALSE),
   fEMCALGeoName("EMCAL_FIRSTYEAR"), fNCellsFromEMCALBorder(0),
   fRemoveBadChannels(kFALSE),fEMCALBadChannelMap(0x0),
-  fRecalibration(kFALSE),fEMCALRecalibrationFactors(),
   fRecoUtils(new AliEMCALRecoUtils),
   fNbins(300), fMinBin(0.), fMaxBin(300.),fOutputContainer(0x0),
   fHmgg(0x0),           fHmggDifferentSM(0x0), 
@@ -109,12 +108,6 @@ AliAnalysisTaskEMCALPi0CalibSelection::~AliAnalysisTaskEMCALPi0CalibSelection()
                delete  fEMCALBadChannelMap;
        }
        
-       
-       if(fEMCALRecalibrationFactors) { 
-               fEMCALRecalibrationFactors->Clear();
-               delete  fEMCALRecalibrationFactors;
-       }       
-    
   if(fRecoUtils) delete fRecoUtils ;
 
 }
@@ -138,7 +131,7 @@ void AliAnalysisTaskEMCALPi0CalibSelection::LocalInit()
        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,fRecalibration, fOldAOD, fSameSM) ;
+            fRemoveBadChannels,fCopyAOD,fRecoUtils->IsRecalibrationOn(), fOldAOD, fSameSM) ;
        fCuts->Add(new TObjString(onePar));
        snprintf(onePar,buffersize, "EMCAL Geometry name: < %s >",fEMCALGeoName.Data()) ;
        fCuts->Add(new TObjString(onePar));
@@ -613,6 +606,11 @@ 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;
@@ -679,10 +677,11 @@ void AliAnalysisTaskEMCALPi0CalibSelection::UserExec(Option_t* /* option */)
   //  }
   
   if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Will use fLogWeight %.3f .\n",fLogWeight);
-  
+  Int_t absId1   = -1;
   Int_t iSupMod1 = -1;
   Int_t iphi1    = -1;
   Int_t ieta1    = -1;
+  Int_t absId2   = -1;
   Int_t iSupMod2 = -1;
   Int_t iphi2    = -1;
   Int_t ieta2    = -1;
@@ -727,24 +726,14 @@ void AliAnalysisTaskEMCALPi0CalibSelection::UserExec(Option_t* /* option */)
     //clu1.Recalibrate(fCalibData, emCells, fEMCALGeoName);
     //clu1.EvalEnergy();
     //clu1.EvalAll(fLogWeight, fEMCALGeoName);
-    if(IsRecalibrationOn())    {
-      Float_t energy = RecalibrateClusterEnergy(c1, emCells);
-      //clu1.SetE(energy);
-      c1->SetE(energy);
-    }
+    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, clu1.GetDispersion(),clu1.GetM02(),clu1.GetM20()); 
       printf("Recal: i %d, E %f, dispersion %f, m02 %f, m20 %f\n",iClu,e1ii, c1->GetDispersion(),c1->GetM02(),c1->GetM20());    
-      Float_t pos2[]={0,0,0};
-      //clu1.GetPosition(pos2);
-      c1->GetPosition(pos2);
-      printf("Recal: i %d, x %f, y %f, z %f\n",iClu, pos2[0], pos2[1], pos2[2]);
-    }
+    
     
     //clu1.GetMomentum(p1,v);
     
@@ -759,7 +748,9 @@ void AliAnalysisTaskEMCALPi0CalibSelection::UserExec(Option_t* /* option */)
     //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]);
-    GetMaxEnergyCellPosAndClusterPos(emCells,c1,iSupMod1,ieta1,iphi1);
+    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]);
@@ -780,20 +771,17 @@ void AliAnalysisTaskEMCALPi0CalibSelection::UserExec(Option_t* /* option */)
       //clu2.Recalibrate(fCalibData, emCells,fEMCALGeoName);
       //clu2.EvalEnergy();
       //clu2.EvalAll(fLogWeight,fEMCALGeoName);
-      if(IsRecalibrationOn())  {
-        Float_t energy = RecalibrateClusterEnergy(c2, emCells);
-        //clu2.SetE(energy);
-        c2->SetE(energy);
-      }
-      
+      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    
-      GetMaxEnergyCellPosAndClusterPos(emCells,c2,iSupMod2,ieta2,iphi2);
-  
+      fRecoUtils->RecalculateClusterPosition(fEMCALGeo, emCells,c2);
+      fRecoUtils->GetMaxEnergyCell          (fEMCALGeo, emCells,c2,absId2,iSupMod2,ieta2,iphi2);
+      
       c2->GetMomentum(p2,v);
 
       p12 = p1+p2;
@@ -832,11 +820,11 @@ void AliAnalysisTaskEMCALPi0CalibSelection::UserExec(Option_t* /* option */)
             //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.};
-            fEMCALGeo->RecalculateTowerPosition(11.5, 23.5, iSupMod1, p1.E(), 0,
-                                                fRecoUtils->GetMisalTransShiftArray(),fRecoUtils->GetMisalRotShiftArray(),posSM1cen); 
+            Float_t depth = fRecoUtils->GetDepth(p1.Energy(),fRecoUtils->GetParticleType(),iSupMod1);
+            fEMCALGeo->RecalculateTowerPosition(11.5, 23.5, iSupMod1, depth, fRecoUtils->GetMisalTransShiftArray(),fRecoUtils->GetMisalRotShiftArray(),posSM1cen); 
             Float_t posSM2cen[3]={0.,0.,0.}; 
-            fEMCALGeo->RecalculateTowerPosition(11.5, 23.5, iSupMod2, p2.E(), 0,
-                                                fRecoUtils->GetMisalTransShiftArray(),fRecoUtils->GetMisalRotShiftArray(),posSM2cen); 
+            depth = fRecoUtils->GetDepth(p2.Energy(),fRecoUtils->GetParticleType(),iSupMod2);
+            fEMCALGeo->RecalculateTowerPosition(11.5, 23.5, iSupMod2, depth, fRecoUtils->GetMisalTransShiftArray(),fRecoUtils->GetMisalRotShiftArray(),posSM2cen); 
             //printf("SM1 %d pos (%2.3f,%2.3f,%2.3f) \n",iSupMod1,posSM1cen[0],posSM1cen[1],posSM1cen[2]);
             //printf("SM2 %d pos (%2.3f,%2.3f,%2.3f) \n",iSupMod2,posSM2cen[0],posSM2cen[1],posSM2cen[2]);
             
@@ -1016,88 +1004,6 @@ Bool_t AliAnalysisTaskEMCALPi0CalibSelection::CheckCellFiducialRegion(AliVCluste
        
 }      
 
-//__________________________________________________
-void AliAnalysisTaskEMCALPi0CalibSelection::GetMaxEnergyCellPosAndClusterPos(AliVCaloCells* cells, AliVCluster* clu, Int_t& iSupMod, Int_t& ieta, Int_t& iphi)
-{
-  //For a given CaloCluster calculates the absId of the cell 
-  //with maximum energy deposit.
-  
-  Double_t eMax       = -1.;
-  Double_t eCell      = -1.;
-  Float_t  fraction   = 1.;
-  Int_t    cellAbsId  = -1;
-  Float_t recalFactor = 1.;
-       
-  Int_t maxId   = -1;
-  Int_t iTower  = -1;
-  Int_t iIphi   = -1;
-  Int_t iIeta   = -1;
-       
-  Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
-  Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
-  Bool_t  areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now
-  Int_t   startingSM = -1;
-  
-  for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
-    cellAbsId = clu->GetCellAbsId(iDig);
-    fraction  = clu->GetCellAmplitudeFraction(iDig);
-    if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
-    Int_t imodrc   = -1, iphirc  = -1, ietarc  =-1;
-    Int_t iTowerrc = -1, iIphirc = -1, iIetarc =-1;
-    fEMCALGeo->GetCellIndex(cellAbsId,imodrc,iTowerrc,iIphirc,iIetarc); 
-    fEMCALGeo->GetCellPhiEtaIndexInSModule(imodrc,iTowerrc,iIphirc, iIetarc,iphirc,ietarc);                    
-    if     (iDig==0)  startingSM = imodrc;
-    else if(imodrc != startingSM) areInSameSM = kFALSE;
-
-    if(IsRecalibrationOn()) {
-      recalFactor = GetEMCALChannelRecalibrationFactor(imodrc,ietarc,iphirc);
-    }
-    eCell  = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
-    
-    weight = TMath::Log(eCell/clEnergy) + 4;
-    if(weight < 0) weight = 0;
-    totalWeight += weight;
-    weightedCol += ietarc*weight;
-    weightedRow += iphirc*weight;
-    
-    //printf("Max cell? cell %d, amplitude org %f, fraction %f, recalibration %f, amplitude new %f \n",cellAbsId, cells->GetCellAmplitude(cellAbsId), fraction, recalFactor, eCell) ;
-    
-    if(eCell > eMax)  { 
-      eMax  = eCell; 
-      maxId = cellAbsId;
-      //printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell);
-    }
-  }// cell loop
-  
-  //Get from the absid the supermodule, tower and eta/phi numbers
-  fEMCALGeo->GetCellIndex(maxId,iSupMod,iTower,iIphi,iIeta); 
-  //Gives SuperModule and Tower numbers
-  fEMCALGeo->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
-                                        iIphi, iIeta,iphi,ieta); 
-  
-  Float_t xyzNew[3];
-  if(areInSameSM == kTRUE) {
-    //printf("In Same SM\n");
-    weightedCol = weightedCol/totalWeight;
-    weightedRow = weightedRow/totalWeight;
-    
-    //Float_t *xyzNew = RecalculatePosition(weightedRow, weightedCol, clEnergy, 0, iSupMod); //1 = electrons, 0 photons
-    fEMCALGeo->RecalculateTowerPosition(weightedRow, weightedCol, iSupMod, clEnergy, 0, //1 = electrons, 0 photons
-                                        fRecoUtils->GetMisalTransShiftArray(), fRecoUtils->GetMisalRotShiftArray(), xyzNew);
-  }
-  else {
-    //printf("In Different SM\n");
-    //Float_t *xyzNew = RecalculatePosition(iphi,        ieta,        clEnergy, 0, iSupMod); //1 = electrons, 0 photons
-    fEMCALGeo->RecalculateTowerPosition(iphi, ieta, iSupMod, clEnergy, 0, //1 = electrons, 0 photons
-                                        fRecoUtils->GetMisalTransShiftArray(), fRecoUtils->GetMisalRotShiftArray(), xyzNew);
-    
-  }
-
-  clu->SetPosition(xyzNew);
-
-  //printf("\t Max : cell %d, iSupMod %d, ieta %d, iphi %d \n",maxId,iSupMod, ieta,iphi);
-  
-}
 
 //__________________________________________________
 //void AliAnalysisTaskEMCALPi0CalibSelection::SetCalibCorrections(AliEMCALCalibData* const cdata)
@@ -1160,74 +1066,3 @@ Bool_t AliAnalysisTaskEMCALPi0CalibSelection::ClusterContainsBadChannel(UShort_t
        
 }
 
-
-//________________________________________________________________
-void AliAnalysisTaskEMCALPi0CalibSelection::InitEMCALRecalibrationFactors(){
-       //Init EMCAL recalibration factors
-       if(DebugLevel() > 0 )printf("AliAnalysisTaskEMCALPi0CalibSelection::InitEMCALRecalibrationFactors()\n");
-       //In order to avoid rewriting the same histograms
-       Bool_t oldStatus = TH1::AddDirectoryStatus();
-       TH1::AddDirectory(kFALSE);
-       
-       fEMCALRecalibrationFactors = new TObjArray(12);
-       for (int i = 0; i < 12; i++) fEMCALRecalibrationFactors->Add(new TH2F(Form("EMCALRecalFactors_SM%d",i),Form("EMCALRecalFactors_SM%d",i),  48, 0, 48, 24, 0, 24));
-       //Init the histograms with 1
-       for (Int_t sm = 0; sm < 12; sm++) {
-               for (Int_t i = 0; i < 48; i++) {
-                       for (Int_t j = 0; j < 24; j++) {
-                               SetEMCALChannelRecalibrationFactor(sm,i,j,1.);
-                       }
-               }
-       }
-       fEMCALRecalibrationFactors->SetOwner(kTRUE);
-       fEMCALRecalibrationFactors->Compress();
-       
-       //In order to avoid rewriting the same histograms
-       TH1::AddDirectory(oldStatus);           
-}
-
-//________________________________________________________________
-Float_t AliAnalysisTaskEMCALPi0CalibSelection::RecalibrateClusterEnergy(AliAODCaloCluster * cluster, AliAODCaloCells * cells){
-       // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
-       // AOD case
-       
-       if(!cells) {
-               printf("AliAnalysisTaskEMCALPi0CalibSelection::RecalibrateClusterEnergy(AOD) - Cells pointer does not exist, stop!");
-               abort();
-       }
-       
-       //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
-       UShort_t * index    = cluster->GetCellsAbsId() ;
-       Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
-       Int_t ncells = cluster->GetNCells();
-       
-       //Initialize some used variables
-       Float_t energy = 0;
-       Int_t absId    = -1;
-    Int_t icol = -1, irow = -1, imod=1;
-       Float_t factor = 1, frac = 0;
-       
-       //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
-       for(Int_t icell = 0; icell < ncells; icell++){
-               absId = index[icell];
-               frac =  fraction[icell];
-               if(frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off
-               Int_t iTower = -1, iIphi = -1, iIeta = -1; 
-               fEMCALGeo->GetCellIndex(absId,imod,iTower,iIphi,iIeta); 
-               if(fEMCALRecalibrationFactors->GetEntries() <= imod) continue;
-               fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);                     
-               factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
-               if(DebugLevel()>2)
-                       printf("AliAnalysisTaskEMCALPi0CalibSelection::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n",
-                                  imod,icol,irow,frac,factor,cells->GetCellAmplitude(absId));
-               
-               energy += cells->GetCellAmplitude(absId)*factor*frac;
-       }
-       
-       if(DebugLevel()>1)
-               printf("AliAnalysisTaskEMCALPi0CalibSelection::RecalibrateClusterEnergy - Energy before %f, after %f\n",cluster->E(),energy);
-       
-       return energy;
-       
-}
-
index 94377a3..87bbf4d 100644 (file)
@@ -91,29 +91,7 @@ public:
   void   SetEMCALChannelStatusMap(TObjArray *map)  {fEMCALBadChannelMap = map;}
        
   Bool_t ClusterContainsBadChannel(UShort_t* cellList, Int_t nCells);
-       
-  // Recalibration
-  Bool_t IsRecalibrationOn() const { return fRecalibration  ; }
-  void SwitchOnRecalibration()     {fRecalibration = kTRUE  ; InitEMCALRecalibrationFactors();}
-  void SwitchOffRecalibration()    {fRecalibration = kFALSE ; }
-       
-  void InitEMCALRecalibrationFactors() ;
-  
-  Float_t GetEMCALChannelRecalibrationFactor(Int_t iSM , Int_t iCol, Int_t iRow) const { 
-       if(fEMCALRecalibrationFactors) return (Float_t) ((TH2F*)fEMCALRecalibrationFactors->At(iSM))->GetBinContent(iCol,iRow); 
-       else return 1;}
-       
-  void SetEMCALChannelRecalibrationFactor(Int_t iSM , Int_t iCol, Int_t iRow, Double_t c = 1) { 
-       if(!fEMCALRecalibrationFactors) InitEMCALRecalibrationFactors();
-       ((TH2F*)fEMCALRecalibrationFactors->At(iSM))->SetBinContent(iCol,iRow,c);}
-       
-  void SetEMCALChannelRecalibrationFactors(Int_t iSM , TH2F* h) {fEMCALRecalibrationFactors->AddAt(h,iSM);}
-       
-  TH2F * GetEMCALChannelRecalibrationFactors(Int_t iSM) const {return (TH2F*)fEMCALRecalibrationFactors->At(iSM);}
-       
-  void SetEMCALChannelRecalibrationFactors(TObjArray *map) {fEMCALRecalibrationFactors = map;}
-  Float_t RecalibrateClusterEnergy(AliAODCaloCluster* cluster, AliAODCaloCells * cells);
-       
+                 
   void SetEMCALRecoUtils(AliEMCALRecoUtils * ru) {fRecoUtils = ru;}
   AliEMCALRecoUtils* GetEMCALRecoUtils() const   {return fRecoUtils;}
   
@@ -141,9 +119,6 @@ private:
 
   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
-  Bool_t     fRecalibration;             // Switch on or off the recalibration
-  TObjArray *fEMCALRecalibrationFactors; // Array of histograms with map of recalibration factors, EMCAL                 
   AliEMCALRecoUtils * fRecoUtils;  // Access to reconstruction utilities
   
   //Output histograms  
@@ -181,7 +156,7 @@ private:
   TH1I*   fhNEvents;        //! Number of events counter histogram
   TList * fCuts ;           //! List with analysis cuts
 
-  ClassDef(AliAnalysisTaskEMCALPi0CalibSelection,8);
+  ClassDef(AliAnalysisTaskEMCALPi0CalibSelection,9);
 
 };
 
index c86a039..f44f4c6 100644 (file)
@@ -46,8 +46,8 @@ 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.
   //-------------------------------------------------------------------------------------------------
@@ -57,10 +57,10 @@ void anaEMCALCalib(Int_t mode=mLocal)
     cout<<"Wrong  data type "<<kInputData<<endl;
     break;
   }
-  
+
   TChain *chain       = new TChain(kTreeName) ;
   CreateChain(mode, chain);  
-  
+
   if(chain){
     AliLog::SetGlobalLogLevel(AliLog::kError);//Minimum prints on screen
     
@@ -68,61 +68,61 @@ void anaEMCALCalib(Int_t mode=mLocal)
     // Make the analysis manager
     //-------------------------------------
     AliAnalysisManager *mgr  = new AliAnalysisManager("Manager", "Manager");
-    
+
     // AOD output handler
     if(copy){
-      AliAODHandler* aodoutHandler   = new AliAODHandler();
-      aodoutHandler->SetOutputFileName("aod.root");
-      ////aodoutHandler->SetCreateNonStandardAOD();
-      mgr->SetOutputEventHandler(aodoutHandler);
+       AliAODHandler* aodoutHandler   = new AliAODHandler();
+       aodoutHandler->SetOutputFileName("aod.root");
+       ////aodoutHandler->SetCreateNonStandardAOD();
+       mgr->SetOutputEventHandler(aodoutHandler);
     }
-    
+
     //input
     if(kInputData == "ESD")
-    {
-      // ESD handler
-      AliESDInputHandler *esdHandler = new AliESDInputHandler();
-      mgr->SetInputEventHandler(esdHandler);
-      esdHandler->SetReadFriends(kFALSE);
-      cout<<"ESD handler "<<mgr->GetInputEventHandler()<<endl;
-    }
+      {
+       // ESD handler
+       AliESDInputHandler *esdHandler = new AliESDInputHandler();
+       mgr->SetInputEventHandler(esdHandler);
+       esdHandler->SetReadFriends(kFALSE);
+       cout<<"ESD handler "<<mgr->GetInputEventHandler()<<endl;
+      }
     if(kInputData == "AOD")
-    {
-      // AOD handler
-      AliAODInputHandler *aodHandler = new AliAODInputHandler();
-      mgr->SetInputEventHandler(aodHandler);
-      cout<<"AOD handler "<<mgr->GetInputEventHandler()<<endl;
-      
-    }
+      {
+       // AOD handler
+       AliAODInputHandler *aodHandler = new AliAODInputHandler();
+       mgr->SetInputEventHandler(aodHandler);
+       cout<<"AOD handler "<<mgr->GetInputEventHandler()<<endl;
+       
+      }
     
-    //mgr->SetDebugLevel(1);
+    // mgr->SetDebugLevel(1);
     
     //-------------------------------------------------------------------------
     //Define task, put here any other task that you want to use.
     //-------------------------------------------------------------------------
     AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer();
     AliAnalysisDataContainer *coutput1 = mgr->GetCommonOutputContainer();
-    
+
     
     // ESD filter task
     if(kInputData == "ESD"){
-      
+
       gROOT->LoadMacro("AddTaskPhysicsSelection.C");
       AliPhysicsSelectionTask* physSelTask = AddTaskPhysicsSelection();
       if(!copy){
-        gROOT->LoadMacro("AddTaskESDFilter.C");
-        AliAnalysisTaskESDfilter *esdfilter = AddTaskESDFilter(kFALSE);
+       gROOT->LoadMacro("AddTaskESDFilter.C");
+       AliAnalysisTaskESDfilter *esdfilter = AddTaskESDFilter(kFALSE);
       }
     }
     
-    
-    
+
+   
     AliAnalysisTaskEMCALPi0CalibSelection * pi0calib = new AliAnalysisTaskEMCALPi0CalibSelection ("EMCALPi0Calibration");
     //pi0calib->SetDebugLevel(10); 
     pi0calib->CopyAOD(copy);
     pi0calib->SetClusterMinEnergy(0.5);
     pi0calib->SetClusterMaxEnergy(15.);
-    pi0calib->SetAsymmetryCut(0.7);
+    //pi0calib->SetAsymmetryCut(0.7);
     pi0calib->SetClusterMinNCells(0);
     pi0calib->SetNCellsGroup(0);
     pi0calib->SwitchOnBadChannelsRemoval();
@@ -130,25 +130,46 @@ void anaEMCALCalib(Int_t mode=mLocal)
     pi0calib->SwitchOnOldAODs();
     pi0calib->SetNumberOfCellsFromEMCALBorder(1);
     AliEMCALRecoUtils * reco = pi0calib->GetEMCALRecoUtils();
-    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->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->Print("");
-    
+
+//     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->Print("");
+
     // 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);   
@@ -166,35 +187,23 @@ void anaEMCALCalib(Int_t mode=mLocal)
     pi0calib->SetEMCALChannelStatus(2,19,22);
     //SM3
     pi0calib->SetEMCALChannelStatus(3,4,7);
-    
+
     mgr->AddTask(pi0calib);
     
-    pi0calib->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();
-    
-    pi0calib->SetEMCALChannelRecalibrationFactors(0,h0);
-    pi0calib->SetEMCALChannelRecalibrationFactors(1,h1);
-    pi0calib->SetEMCALChannelRecalibrationFactors(2,h2);
-    pi0calib->SetEMCALChannelRecalibrationFactors(3,h3);
-    
     // 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);
-    
+
     //-----------------------
     // Run the analysis
     //-----------------------    
@@ -209,9 +218,9 @@ void anaEMCALCalib(Int_t mode=mLocal)
     mgr->InitAnalysis();
     mgr->PrintStatus();
     mgr->StartAnalysis(smode.Data(),chain);
-    
-    cout <<" Analysis ended sucessfully "<< endl ;
-    
+
+cout <<" Analysis ended sucessfully "<< endl ;
+
   }
   else cout << "Chain was not produced ! "<<endl;
   
@@ -244,7 +253,7 @@ void  LoadLibraries(const anaModes mode) {
     gSystem->Load("libANALYSISalice.so");
     TGeoManager::Import("geometry.root") ; //need file "geometry.root" in local dir!!!!
     gSystem->Load("libPWG4CaloCalib.so");
-    
+    //SetupPar("PWG4CaloCalib");
     /*
       //     gSystem->Load("libPWG4omega3pi.so");
       //     gSystem->Load("libCORRFW.so");
index 85a5220..3a82904 100755 (executable)
@@ -648,8 +648,8 @@ void AliCaloTrackReader::FillInputEMCAL() {
             printf("AliCaloTrackReader::FillInputEMCAL() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
                    momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
           
-          Float_t pos[3];
-          clus->GetPosition(pos);
+          //Float_t pos[3];
+          //clus->GetPosition(pos);
           //printf("Before Corrections: e %f, x %f, y %f, z %f\n",clus->E(),pos[0],pos[1],pos[2]);
           
           //Recalibrate the cluster energy 
@@ -668,7 +668,7 @@ void AliCaloTrackReader::FillInputEMCAL() {
           //Recalculate cluster position
           if(GetCaloUtils()->IsRecalculationOfClusterPositionOn()){
             GetCaloUtils()->RecalculateClusterPosition(GetEMCALCells(),clus); 
-            clus->GetPosition(pos);
+            //clus->GetPosition(pos);
             //printf("After  Corrections: e %f, x %f, y %f, z %f\n",clus->E(),pos[0],pos[1],pos[2]);
           }
           
index d5e1636..7bd86c9 100755 (executable)
@@ -35,8 +35,6 @@
 #include "AliVCluster.h"
 #include "AliVCaloCells.h"
 #include "AliMixedEvent.h"
-#include "AliEMCALRecoUtils.h"
-
 
 ClassImp(AliCalorimeterUtils)
   
@@ -51,7 +49,7 @@ ClassImp(AliCalorimeterUtils)
     fEMCALBadChannelMap(0x0),fPHOSBadChannelMap(0x0), 
     fNCellsFromEMCALBorder(0), fNCellsFromPHOSBorder(0), 
     fNoEMCALBorderAtEta0(kFALSE),fRecalibration(kFALSE),
-    fEMCALRecalibrationFactors(), fPHOSRecalibrationFactors(),
+    fPHOSRecalibrationFactors(),
     fEMCALRecoUtils(new AliEMCALRecoUtils),fRecalculatePosition(kFALSE),fCorrectELinearity(kFALSE)
 
 {
@@ -77,10 +75,6 @@ AliCalorimeterUtils::~AliCalorimeterUtils() {
     delete  fPHOSBadChannelMap;
   }
        
-  if(fEMCALRecalibrationFactors) { 
-    fEMCALRecalibrationFactors->Clear();
-    delete  fEMCALRecalibrationFactors;
-  }
   if(fPHOSRecalibrationFactors) { 
     fPHOSRecalibrationFactors->Clear();
     delete  fPHOSRecalibrationFactors;
@@ -486,31 +480,6 @@ void AliCalorimeterUtils::InitPHOSBadChannelStatusMap(){
 }
 
 //________________________________________________________________
-void AliCalorimeterUtils::InitEMCALRecalibrationFactors(){
-       //Init EMCAL recalibration factors
-       if(fDebug > 0 )printf("AliCalorimeterUtils::InitEMCALRecalibrationFactors()\n");
-       //In order to avoid rewriting the same histograms
-       Bool_t oldStatus = TH1::AddDirectoryStatus();
-       TH1::AddDirectory(kFALSE);
-
-       fEMCALRecalibrationFactors = new TObjArray(12);
-       for (int i = 0; i < 12; i++) fEMCALRecalibrationFactors->Add(new TH2F(Form("EMCALRecalFactors_SM%d",i),Form("EMCALRecalFactors_SM%d",i),  48, 0, 48, 24, 0, 24));
-       //Init the histograms with 1
-       for (Int_t sm = 0; sm < 12; sm++) {
-               for (Int_t i = 0; i < 48; i++) {
-                       for (Int_t j = 0; j < 24; j++) {
-                               SetEMCALChannelRecalibrationFactor(sm,i,j,1.);
-                       }
-               }
-       }
-       fEMCALRecalibrationFactors->SetOwner(kTRUE);
-       fEMCALRecalibrationFactors->Compress();
-       
-       //In order to avoid rewriting the same histograms
-       TH1::AddDirectory(oldStatus);           
-}
-
-//________________________________________________________________
 void AliCalorimeterUtils::InitPHOSRecalibrationFactors(){
        //Init EMCAL recalibration factors
        if(fDebug > 0 )printf("AliCalorimeterUtils::InitPHOSRecalibrationFactors()\n");
@@ -679,76 +648,6 @@ void AliCalorimeterUtils::RecalculateClusterPosition(AliVCaloCells* cells, AliVC
   
   //Recalculate EMCAL cluster position
   
-  Double_t eMax       = -1.;
-  Double_t eCell      = -1.;
-  Float_t  fraction   = 1.;
-  Int_t    cellAbsId  = -1;
-  Float_t recalFactor = 1.;
-       
-  Int_t maxId   = -1;
-  Int_t imod   = -1, iphi  = -1, ieta  =-1;
-  Int_t iTower = -1, iIphi = -1, iIeta =-1;
-
-  Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
-  Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
-  Bool_t  areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now
-  Int_t   startingSM = -1;
-  
-  for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
-    cellAbsId = clu->GetCellAbsId(iDig);
-    fraction  = clu->GetCellAmplitudeFraction(iDig);
-    if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
-       fEMCALGeo->GetCellIndex(cellAbsId,imod,iTower,iIphi,iIeta); 
-    fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);                        
-    if     (iDig==0)  startingSM = imod;
-    else if(imod != startingSM) areInSameSM = kFALSE;
-    
-    if(IsRecalibrationOn()) {
-      recalFactor = GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
-    }
-    eCell  = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
-    
-    weight = TMath::Log(eCell/clEnergy) + 4;
-    if(weight < 0) weight = 0;
-    totalWeight += weight;
-    weightedCol += ieta*weight;
-    weightedRow += iphi*weight;
-    
-    //printf("Max cell? cell %d, amplitude org %f, fraction %f, recalibration %f, amplitude new %f \n",cellAbsId, cells->GetCellAmplitude(cellAbsId), fraction, recalFactor, eCell) ;
-    
-    if(eCell > eMax)  { 
-      eMax  = eCell; 
-      maxId = cellAbsId;
-      //printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell);
-    }
-  }// cell loop
-  
-  //Get from the absid the supermodule, tower and eta/phi numbers
-  fEMCALGeo->GetCellIndex(maxId,imod,iTower,iIphi,iIeta); 
-  //Gives SuperModule and Tower numbers
-  fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,
-                                         iIphi, iIeta,iphi,ieta); 
-  
-  Float_t xyzNew[3];
-  if(areInSameSM == kTRUE) {
-    //printf("In Same SM\n");
-    weightedCol = weightedCol/totalWeight;
-    weightedRow = weightedRow/totalWeight;
-    
-    //Float_t *xyzNew = RecalculatePosition(weightedRow, weightedCol, clEnergy, 0, iSupMod); //1 = electrons, 0 photons
-    fEMCALGeo->RecalculateTowerPosition(weightedRow, weightedCol, imod, clEnergy, 0, //1 = electrons, 0 photons
-                                        fEMCALRecoUtils->GetMisalTransShiftArray(), fEMCALRecoUtils->GetMisalRotShiftArray(), xyzNew);
-  }
-  else {
-    //printf("In Different SM\n");
-    //Float_t *xyzNew = RecalculatePosition(iphi,        ieta,        clEnergy, 0, iSupMod); //1 = electrons, 0 photons
-    fEMCALGeo->RecalculateTowerPosition(iphi, ieta, imod, clEnergy, 0, //1 = electrons, 0 photons
-                                        fEMCALRecoUtils->GetMisalTransShiftArray(), fEMCALRecoUtils->GetMisalRotShiftArray(), xyzNew);
-    
-  }
-  
-  clu->SetPosition(xyzNew);
-  
-  //printf("\t Max : cell %d, iSupMod %d, ieta %d, iphi %d \n",maxId,iSupMod, ieta,iphi);
-  
+  fEMCALRecoUtils->RecalculateClusterPosition((AliEMCALGeometry*)fEMCALGeo, cells,clu);
+
 }
index 2a08022..71582ed 100755 (executable)
@@ -26,7 +26,7 @@ class AliVCluster;
 class AliVCaloCells;
 #include "AliPHOSGeoUtils.h"
 #include "AliEMCALGeoUtils.h"
-class AliEMCALRecoUtils;
+#include "AliEMCALRecoUtils.h"
 
 class AliCalorimeterUtils : public TObject {
 
@@ -112,35 +112,31 @@ class AliCalorimeterUtils : public TObject {
        
   // Recalibration
   Bool_t IsRecalibrationOn()  const { return fRecalibration ; }
-  void SwitchOnRecalibration()    {fRecalibration = kTRUE ; InitEMCALRecalibrationFactors(); InitPHOSRecalibrationFactors();}
-  void SwitchOffRecalibration()   {fRecalibration = kFALSE ; }
+  void SwitchOnRecalibration()    {fRecalibration = kTRUE ; InitPHOSRecalibrationFactors(); fEMCALRecoUtils->SwitchOnRecalibration();}
+  void SwitchOffRecalibration()   {fRecalibration = kFALSE;fEMCALRecoUtils->SwitchOffRecalibration();}
        
-  void InitEMCALRecalibrationFactors() ;
   void InitPHOSRecalibrationFactors () ;
        
-  Float_t GetEMCALChannelRecalibrationFactor(Int_t iSM , Int_t iCol, Int_t iRow) const { 
-    if(fEMCALRecalibrationFactors) return (Float_t) ((TH2F*)fEMCALRecalibrationFactors->At(iSM))->GetBinContent(iCol,iRow); 
-    else return 1;}
+  Float_t GetEMCALChannelRecalibrationFactor(Int_t iSM , Int_t iCol, Int_t iRow) const { return fEMCALRecoUtils->GetEMCALChannelRecalibrationFactor(iSM , iCol, iRow);}
   
   Float_t GetPHOSChannelRecalibrationFactor (Int_t imod, Int_t iCol, Int_t iRow) const { 
     if(fPHOSRecalibrationFactors)return (Float_t) ((TH2F*)fPHOSRecalibrationFactors->At(imod))->GetBinContent(iCol,iRow); 
     else return 1;}
   
   void SetEMCALChannelRecalibrationFactor(Int_t iSM , Int_t iCol, Int_t iRow, Double_t c = 1) { 
-    if(!fEMCALRecalibrationFactors) InitEMCALRecalibrationFactors();
-    ((TH2F*)fEMCALRecalibrationFactors->At(iSM))->SetBinContent(iCol,iRow,c);}
+    fEMCALRecoUtils->SetEMCALChannelRecalibrationFactor(iSM,iCol,iRow,c);}
        
   void SetPHOSChannelRecalibrationFactor (Int_t imod, Int_t iCol, Int_t iRow, Double_t c = 1) {
     if(!fPHOSRecalibrationFactors)  InitPHOSRecalibrationFactors();
     ((TH2F*)fPHOSRecalibrationFactors->At(imod))->SetBinContent(iCol,iRow,c);}
     
-  void SetEMCALChannelRecalibrationFactors(Int_t iSM , TH2F* h) {fEMCALRecalibrationFactors->AddAt(h,iSM);}
+  void SetEMCALChannelRecalibrationFactors(Int_t iSM , TH2F* h) {fEMCALRecoUtils->SetEMCALChannelRecalibrationFactors(iSM,h);}
   void SetPHOSChannelRecalibrationFactors(Int_t imod , TH2F* h) {fPHOSRecalibrationFactors ->AddAt(h,imod);}
        
-  TH2F * GetEMCALChannelRecalibrationFactors(Int_t iSM) const {return (TH2F*)fEMCALRecalibrationFactors->At(iSM);}
+  TH2F * GetEMCALChannelRecalibrationFactors(Int_t iSM) const {return fEMCALRecoUtils->GetEMCALChannelRecalibrationFactors(iSM);}
   TH2F * GetPHOSChannelRecalibrationFactors(Int_t imod) const {return (TH2F*)fPHOSRecalibrationFactors->At(imod);}
        
-  void SetEMCALChannelRecalibrationFactors(TObjArray *map) {fEMCALRecalibrationFactors = map;}
+  void SetEMCALChannelRecalibrationFactors(TObjArray *map) {fEMCALRecoUtils->SetEMCALChannelRecalibrationFactors(map);}
   void SetPHOSChannelRecalibrationFactors (TObjArray *map) {fPHOSRecalibrationFactors  = map;}
 
   Float_t RecalibrateClusterEnergy(AliVCluster* cluster, AliVCaloCells * cells);
@@ -174,7 +170,6 @@ class AliCalorimeterUtils : public TObject {
   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        * fEMCALRecalibrationFactors; // Array of histograms with map of recalibration factors, EMCAL
   TObjArray        * fPHOSRecalibrationFactors;  // Array of histograms with map of recalibration factors, PHOS
   AliEMCALRecoUtils* fEMCALRecoUtils;        //  EMCAL utils for cluster rereconstruction
   Bool_t             fRecalculatePosition;   // Recalculate cluster position