Many improvements in EMCAL calibration class:
authorgconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 23 Jul 2010 10:12:19 +0000 (10:12 +0000)
committergconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 23 Jul 2010 10:12:19 +0000 (10:12 +0000)
  - Added possibility to remove bad channels/clusters by hand
  - Added possibility to recalibrate channels on the flight
  - Fill always max energy channel invariant mass for 2 clusters of the pair
  - Added possiblity to group cells
  - Set histograms range and binning via setters
  - Publish parameters in output file

Changes in PROOF*/SETUP.C to ensure that the new par file compiled library is the one used.

PWG4/CaloCalib/AliAnalysisTaskEMCALPi0CalibSelection.cxx
PWG4/CaloCalib/AliAnalysisTaskEMCALPi0CalibSelection.h
PWG4/CaloCalib/macros/anaEMCALCalib.C [new file with mode: 0644]
PWG4/PROOF-INF.PWG4CaloCalib/SETUP.C
PWG4/PROOF-INF.PWG4PartCorrBase/SETUP.C
PWG4/PROOF-INF.PWG4PartCorrDep/SETUP.C

index 719d2a2..291b056 100644 (file)
 #include "AliEMCALGeometry.h"
 #include "AliAODCaloCluster.h"
 #include "AliAODCaloCells.h"
-#include "AliEMCALAodCluster.h"
-#include "AliEMCALCalibData.h"
+//#include "AliEMCALAodCluster.h"
+//#include "AliEMCALCalibData.h"
 
 ClassImp(AliAnalysisTaskEMCALPi0CalibSelection)
 
 
 //__________________________________________________
 AliAnalysisTaskEMCALPi0CalibSelection::AliAnalysisTaskEMCALPi0CalibSelection(const char* name) :
-  AliAnalysisTaskSE(name),fEMCALGeo(0x0),fCalibData(0x0), fEmin(1.), fEmax(15.), fMinNCells(2), fLogWeight(4.5), fCopyAOD(kFALSE), 
-  fEMCALGeoName("EMCAL_FIRSTYEAR"), fOutputContainer(0x0),fHmgg(0x0)
+  AliAnalysisTaskSE(name),fEMCALGeo(0x0),//fCalibData(0x0), 
+  fEmin(0.5), fEmax(15.), fMinNCells(2), fGroupNCells(0),
+  fLogWeight(4.5), fCopyAOD(kFALSE), fEMCALGeoName("EMCAL_FIRSTYEAR"),     
+  fRemoveBadChannels(kFALSE),fEMCALBadChannelMap(0x0),
+  fRecalibration(kFALSE),fEMCALRecalibrationFactors(), 
+  fNbins(300), fMinBin(0.), fMaxBin(300.),
+  fOutputContainer(0x0),fHmgg(0x0), fhNEvents(0x0),fCuts(0x0)
 {
   //Named constructor which should be used.
   
-  DefineOutput(1,TList::Class());
-  
   for(Int_t iMod=0; iMod < 12; iMod++) {
     for(Int_t iX=0; iX<24; iX++) {
       for(Int_t iZ=0; iZ<48; iZ++) {
@@ -62,8 +65,9 @@ AliAnalysisTaskEMCALPi0CalibSelection::AliAnalysisTaskEMCALPi0CalibSelection(con
       }
     } 
   }
-       
+  
   DefineOutput(1, TList::Class());
+  DefineOutput(2, TList::Class());  // will contain cuts or local params
 
 }
 
@@ -77,11 +81,46 @@ AliAnalysisTaskEMCALPi0CalibSelection::~AliAnalysisTaskEMCALPi0CalibSelection()
     delete fOutputContainer ;
   }
        
-  if(fCalibData)  delete fCalibData;
+  //if(fCalibData)  delete fCalibData;
   if(fEMCALGeo)   delete fEMCALGeo;
        
+       if(fEMCALBadChannelMap) { 
+               fEMCALBadChannelMap->Clear();
+               delete  fEMCALBadChannelMap;
+       }
+       
+       
+       if(fEMCALRecalibrationFactors) { 
+               fEMCALRecalibrationFactors->Clear();
+               delete  fEMCALRecalibrationFactors;
+       }       
 }
 
+//_____________________________________________________
+void AliAnalysisTaskEMCALPi0CalibSelection::LocalInit()
+{
+       // Local Initialization
+       
+       // Create cuts/param objects and publish to slot
+       
+       char onePar[255] ;
+       fCuts = new TList();
+
+       sprintf(onePar,"Custer cuts: %2.2f < E < %2.2f GeV; min number of cells %d;", fEmin,fEmax, fMinNCells) ;
+       fCuts->Add(new TObjString(onePar));
+       sprintf(onePar,"Group %d cells;", fGroupNCells) ;
+       fCuts->Add(new TObjString(onePar));
+       sprintf(onePar,"Histograms: bins %d; energy range: %2.2f < E < %2.2f GeV;",fNbins,fMinBin,fMaxBin) ;
+       fCuts->Add(new TObjString(onePar));
+       sprintf(onePar,"Switchs: Remove Bad Channels? %d; Copy AODs? %d;  Recalibrate %d?",fRemoveBadChannels,fCopyAOD,fRecalibration) ;
+       fCuts->Add(new TObjString(onePar));
+       sprintf(onePar,"EMCAL Geometry name: < %s >",fEMCALGeoName.Data()) ;
+       fCuts->Add(new TObjString(onePar));
+
+       // Post Data
+       PostData(2, fCuts);
+       
+}
 
 //__________________________________________________
 void AliAnalysisTaskEMCALPi0CalibSelection::CreateAODFromAOD()
@@ -154,11 +193,13 @@ void AliAnalysisTaskEMCALPi0CalibSelection::CreateAODFromAOD()
     //Check if it is a EMCAL cluster
     if(!cluster->IsEMCALCluster())  continue ;
     
+       if(ClusterContainsBadChannel(cluster->GetCellsAbsId(), cluster->GetNCells())) continue; 
+
+         
     Int_t id       = cluster->GetID();
     Float_t energy = cluster->E();
     cluster->GetPosition(posF);
     Char_t ttype   = cluster->GetType(); 
-    printf("set pos1\n");
     AliAODCaloCluster *caloCluster = new(caloClusters[jClusters++]) 
       AliAODCaloCluster(id,
                        0,
@@ -167,7 +208,7 @@ void AliAnalysisTaskEMCALPi0CalibSelection::CreateAODFromAOD()
                        posF,
                        NULL,
                        ttype);
-    printf("set pos2\n");
+
     caloCluster->SetCaloCluster(cluster->GetDistToBadChannel(),
                                cluster->GetDispersion(),
                                cluster->GetM20(), cluster->GetM02(),
@@ -178,17 +219,11 @@ void AliAnalysisTaskEMCALPi0CalibSelection::CreateAODFromAOD()
     caloCluster->SetNCells(cluster->GetNCells());
     caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
        
-       caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
-
-//     else{
-//             Double_t *fraction = cluster->GetCellsAmplitudeFraction();
-//             for(Int_t i = 0; i < cluster->GetNCells() ; i++) fraction[i] = 1;
-//             caloCluster->SetCellsAmplitudeFraction(fraction);
-//     }
+    caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
     
   } 
 
-  caloClusters.Expand(jClusters); // resize TObjArray to 'remove' slots for pseudo clusters     
+  caloClusters.Expand(jClusters); // resize TObjArray   
   // end of loop on calo clusters
   
   // fill EMCAL cell info
@@ -201,8 +236,6 @@ void AliAnalysisTaskEMCALPi0CalibSelection::CreateAODFromAOD()
     aodEMcells.SetType(AliAODCaloCells::kEMCAL);
        
        Double_t calibFactor = 1;
-       //if(fOldData) calibFactor = 0.0153;
-
     for (Int_t iCell = 0; iCell < nEMcell; iCell++) {      
       aodEMcells.SetCell(iCell,aodinEMcells.GetCellNumber(iCell),aodinEMcells.GetAmplitude(iCell)*calibFactor);
     }
@@ -284,7 +317,8 @@ void AliAnalysisTaskEMCALPi0CalibSelection::CreateAODFromESD()
     
     //Check which calorimeter information we want to keep.
     if(!cluster->IsEMCAL())  continue ;
-    
+       if(ClusterContainsBadChannel(cluster->GetCellsAbsId(), cluster->GetNCells())) continue; 
+         
     Int_t id       = cluster->GetID();
     Float_t energy = cluster->E();
     cluster->GetPosition(posF);
@@ -307,16 +341,11 @@ void AliAnalysisTaskEMCALPi0CalibSelection::CreateAODFromESD()
     caloCluster->SetPIDFromESD(cluster->GetPid());
     caloCluster->SetNCells(cluster->GetNCells());
     caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
-//     if(!fOldData){
-                 caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
-//     }
-//     else{
-//               Double_t *fraction = cluster->GetCellsAmplitudeFraction();
-//               for(Int_t i = 0; i < cluster->GetNCells() ; i++) fraction[i] = 1;
-//               caloCluster->SetCellsAmplitudeFraction(fraction);
-//       }    
+    caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
+
   } 
-  caloClusters.Expand(jClusters); // resize TObjArray to 'remove' slots for pseudo clusters     
+
+  caloClusters.Expand(jClusters); // resize TObjArray
   // end of loop on calo clusters
   
   // fill EMCAL cell info
@@ -329,9 +358,7 @@ void AliAnalysisTaskEMCALPi0CalibSelection::CreateAODFromESD()
     aodEMcells.CreateContainer(nEMcell);
     aodEMcells.SetType(AliAODCaloCells::kEMCAL);  
          
-       Double_t calibFactor = 1;
-       //if(fOldData) calibFactor = 0.0153;
-   
+       Double_t calibFactor = 1;   
        for (Int_t iCell = 0; iCell < nEMcell; iCell++) {      
       aodEMcells.SetCell(iCell,esdEMcells.GetCellNumber(iCell),esdEMcells.GetAmplitude(iCell)*calibFactor);
     }
@@ -346,32 +373,31 @@ void AliAnalysisTaskEMCALPi0CalibSelection::UserCreateOutputObjects()
 {
   //Create output container, init geometry and calibration
        
-  printf("--- Create histograms for calibration with pi0 ---\n");
   fEMCALGeo =  AliEMCALGeometry::GetInstance(fEMCALGeoName) ;  
        
   fOutputContainer = new TList();
   
   char hname[128], htitl[128];
   
-       for(Int_t iMod=0; iMod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); iMod++) {
+  for(Int_t iMod=0; iMod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); iMod++) {
     for(Int_t iRow=0; iRow<AliEMCALGeoParams::fgkEMCALRows; iRow++) {
       for(Int_t iCol=0; iCol<AliEMCALGeoParams::fgkEMCALCols; iCol++) {
        sprintf(hname,"%d_%d_%d",iMod,iCol,iRow);
        sprintf(htitl,"Two-gamma inv. mass for super mod %d, cell(col,row)=(%d,%d)",iMod,iCol,iRow);
-       fHmpi0[iMod][iCol][iRow] = new TH1F(hname,htitl,100,0.,300.);
+       fHmpi0[iMod][iCol][iRow] = new TH1F(hname,htitl,fNbins,fMinBin,fMaxBin);
        fOutputContainer->Add(fHmpi0[iMod][iCol][iRow]);
       }
     }
   }
 
-  fHmgg = new TH1F("hmgg","2-cluster invariant mass",100,0.,300.);
+  fHmgg = new TH1F("hmgg","2-cluster invariant mass",fNbins,fMinBin,fMaxBin);
   fOutputContainer->Add(fHmgg);
-       
-  fCalibData = new AliEMCALCalibData();
-       
-  if(!fCalibData)
-    AliFatal("Calibration parameters not found in CDB!");
-       
+  
+  fhNEvents        = new TH1I("hNEvents", "Number of analyzed events"   , 1 , 0 , 1  ) ;
+  fOutputContainer->Add(fhNEvents);
+
+//  fCalibData = new AliEMCALCalibData();
+               
   PostData(1,fOutputContainer);
 
 }
@@ -382,14 +408,21 @@ void AliAnalysisTaskEMCALPi0CalibSelection::UserExec(Option_t* /* option */)
   //Analysis per event.
   if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection <<< Event %d >>>\n",(Int_t)Entry());
   
+  fhNEvents->Fill(0); //Event analyzed
+  
   AliAODEvent* aod = 0x0;
-  if(!strcmp(InputEvent()->GetName(),"AliAODEvent")) {
+  Bool_t kAOD = kFALSE;
+  if(!strcmp(InputEvent()->GetName(),"AliAODEvent")) kAOD=kTRUE;
+  Bool_t kESD = kFALSE;
+  if(!strcmp(InputEvent()->GetName(),"AliESDEvent")) kESD=kTRUE;
+
+  if(kAOD){
     //Input are ESDs
     aod = dynamic_cast<AliAODEvent*>(InputEvent());
     // Create new AOD with only CaloClusters and CaloCells
     if(fCopyAOD) CreateAODFromAOD();
   }
-  else  if(!strcmp(InputEvent()->GetName(),"AliESDEvent")) {
+  else  if(kESD) {
     //Input are ESDs
     aod = AODEvent();
     // Create AOD with CaloClusters and use it as input.
@@ -401,7 +434,7 @@ void AliAnalysisTaskEMCALPi0CalibSelection::UserExec(Option_t* /* option */)
     abort();
   }    
   
-       Double_t v[3];// = {aod->GetVertex(0)->GetX(),aod->GetVertex(0)->GetY(),aod->GetVertex(0)->GetZ()}; //to check!!
+  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); 
   
@@ -411,30 +444,30 @@ void AliAnalysisTaskEMCALPi0CalibSelection::UserExec(Option_t* /* option */)
   Int_t runNum = aod->GetRunNumber();
   if(DebugLevel() > 1) printf("Run number: %d\n",runNum);
   
+  //FIXME Not neede the matrices for the moment MEFIX
   //Get the matrix with geometry information
   //Still not implemented in AOD, just a workaround to be able to work at least with ESDs      
-  if(!strcmp(InputEvent()->GetName(),"AliAODEvent")) {
-   // if(DebugLevel() > 1) 
-      printf("AliAnalysisTaskEMCALPi0CalibSelection Use ideal geometry, values geometry matrix not kept in AODs.\n");
-  }
-  else{        
-    if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Load Misaligned matrices. \n");
-    AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent()) ;
-    for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){ 
-      if(esd->GetEMCALMatrix(mod)) fEMCALGeo->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
-    }
-  }
+//  if(!strcmp(InputEvent()->GetName(),"AliAODEvent")) {
+//    if(DebugLevel() > 1) 
+//      printf("AliAnalysisTaskEMCALPi0CalibSelection Use ideal geometry, values geometry matrix not kept in AODs.\n");
+//  }
+//  else{      
+//    if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Load Misaligned matrices. \n");
+//    AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent()) ;
+//    for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){ 
+//      if(esd->GetEMCALMatrix(mod)) fEMCALGeo->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
+//    }
+//  }
   
   if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Will use fLogWeight %.3f .\n",fLogWeight);
     
-  Int_t maxId;
-  Int_t iSupMod = -1;
-  Int_t iTower  = -1;
-  Int_t iIphi   = -1;
-  Int_t iIeta   = -1;
-  Int_t iphi   = -1;
-  Int_t ieta   = -1;
-  
+  Int_t iSupMod1 = -1;
+  Int_t iphi1    = -1;
+  Int_t ieta1    = -1;
+  Int_t iSupMod2 = -1;
+  Int_t iphi2    = -1;
+  Int_t ieta2    = -1;
+       
   TLorentzVector p1;
   TLorentzVector p2;
   TLorentzVector p12;
@@ -447,117 +480,303 @@ void AliAnalysisTaskEMCALPi0CalibSelection::UserExec(Option_t* /* option */)
   
   // EMCAL cells
   AliAODCaloCells *emCells = aod->GetEMCALCells();
-       
-       
+  
   // loop over EMCAL clusters
   for(Int_t iClu=0; iClu<kNumberOfEMCALClusters; iClu++) {
     
     AliAODCaloCluster *c1 = (AliAODCaloCluster *) caloClustersArr->At(iClu);
-    if(!c1->IsEMCALCluster()) continue; // EMCAL cluster!
-                 
-       Float_t e1i = c1->E();   // cluster energy before correction   
-       if(e1i < fEmin) continue;
-       else if(e1i > fEmax) continue;
-       else if (c1->GetNCells() < fMinNCells) continue; 
-         
-       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};
-               c1->GetPosition(pos);
-               printf("Std  : i %d, x %f, y %f, z %f\n",iClu, pos[0], pos[1], pos[2]);
-       }
-         
-    AliEMCALAodCluster clu1(*c1);
-    clu1.Recalibrate(fCalibData, emCells, fEMCALGeoName);
-       clu1.EvalEnergy();
+    if(!fCopyAOD && kESD && ClusterContainsBadChannel(c1->GetCellsAbsId(), c1->GetNCells())) continue; 
+    
+    Float_t e1i = c1->E();   // cluster energy before correction   
+    if(e1i < fEmin) continue;
+    else if(e1i > fEmax) continue;
+    else if (c1->GetNCells() < fMinNCells) continue; 
+    
+    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};
+       c1->GetPosition(pos);
+       printf("Std  : i %d, x %f, y %f, z %f\n",iClu, pos[0], pos[1], pos[2]);
+      }
+    
+    //AliEMCALAodCluster clu1(*c1);
+    //clu1.Recalibrate(fCalibData, emCells, fEMCALGeoName);
+    //clu1.EvalEnergy();
     //clu1.EvalAll(fLogWeight, fEMCALGeoName);
-
-       if(DebugLevel() > 2)
-       { 
-               printf("Recal: i %d, E %f, dispersion %f, m02 %f, m20 %f\n",iClu,clu1.E(), clu1.GetDispersion(),clu1.GetM02(),clu1.GetM20());    
-               Float_t pos2[]={0,0,0};
-               clu1.GetPosition(pos2);
-               printf("Recal: i %d, x %f, y %f, z %f\n",iClu, pos2[0], pos2[1], pos2[2]);
-       }
-         
-       clu1.GetMomentum(p1,v);
-
-    //Get tower with maximum energy and fill in the end the pi0 histogram for this cell.
-       MaxEnergyCellPos(emCells,&clu1,maxId);
+    if(IsRecalibrationOn())    {
+      Float_t energy = RecalibrateClusterEnergy(c1, emCells);
+      //clu1.SetE(energy);
+      c1->SetE(energy);
+    }
     
-    //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 e1ii = clu1.E(); // cluster energy after correction
+    Float_t e1ii = c1->E(); // cluster energy after correction
     
-    Float_t e1ii = clu1.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);
+    c1->GetMomentum(p1,v);
+    
+    //Get tower with maximum energy and fill in the end the pi0 histogram for this cell.
+    //MaxEnergyCellPos(emCells,&clu1,iSupMod1,ieta1,iphi1);
+    MaxEnergyCellPos(emCells,c1,iSupMod1,ieta1,iphi1);
     
     for (Int_t jClu=iClu; jClu<kNumberOfEMCALClusters; jClu++) {
       AliAODCaloCluster *c2 = (AliAODCaloCluster *) caloClustersArr->At(jClu);
-      if(!c2->IsEMCALCluster())   continue; // EMCAL cluster!
       if(c2->IsEqual(c1)) continue;
+      if(!fCopyAOD && kESD && ClusterContainsBadChannel(c2->GetCellsAbsId(), c2->GetNCells())) continue;       
       
       Float_t e2i = c2->E();
-//      if(e2i<fEmin) continue;
-               if(e2i > fEmax) continue;
-               else if (c2->GetNCells() < fMinNCells) continue; 
-
-      AliEMCALAodCluster clu2(*c2);
-      clu2.Recalibrate(fCalibData, emCells,fEMCALGeoName);
-         clu2.EvalEnergy();
+      if(e2i < fEmin) continue;
+      else if (e2i > fEmax) continue;
+      else if (c2->GetNCells() < fMinNCells) continue; 
+      
+      //AliEMCALAodCluster clu2(*c2);
+      //clu2.Recalibrate(fCalibData, emCells,fEMCALGeoName);
+      //clu2.EvalEnergy();
       //clu2.EvalAll(fLogWeight,fEMCALGeoName);
-
-      clu2.GetMomentum(p2,v);
-      Float_t e2ii = clu2.E();
+      if(IsRecalibrationOn())  {
+       Float_t energy = RecalibrateClusterEnergy(c2, emCells);
+       //clu2.SetE(energy);
+       c2->SetE(energy);
+      }
+      
+      Float_t e2ii = c2->E();
+      
+      //clu2.GetMomentum(p2,v);
+      c2->GetMomentum(p2,v);
+      
+      //Get tower with maximum energy and fill in the end the pi0 histogram for this cell.
+      //MaxEnergyCellPos(emCells,&clu2,iSupMod2,ieta2,iphi2);
+      MaxEnergyCellPos(emCells,c2,iSupMod2,ieta2,iphi2);
       
       p12 = p1+p2;
       Float_t invmass = p12.M()*1000; 
-
-      fHmgg->Fill(invmass); 
-
-      if(invmass < 300)
-                 fHmpi0[iSupMod][ieta][iphi]->Fill(invmass);
       
-      if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Mass in (SM%d,%d,%d): %.3f GeV  E1_i=%f E1_ii=%f  E2_i=%f E2_ii=%f\n",
-                                 iSupMod,iphi,ieta,p12.M(),e1i,e1ii,e2i,e2ii);
+      if(invmass < fMaxBin && invmass > fMinBin){
+       fHmgg->Fill(invmass); 
+       
+       if (fGroupNCells == 0){
+         fHmpi0[iSupMod1][ieta1][iphi1]->Fill(invmass);
+         fHmpi0[iSupMod2][ieta2][iphi2]->Fill(invmass);
+       }       
+       else  {
+         //printf("Regroup N %d, eta1 %d, phi1 %d, eta2 %d, phi2 %d \n",fGroupNCells, ieta1, iphi1, ieta2, iphi2);
+         for (Int_t i = -fGroupNCells; i < fGroupNCells+1; i++) {
+           for (Int_t j = -fGroupNCells; j < fGroupNCells+1; j++) {
+             //printf("\t i %d, j %d\n",i,j);
+             if((ieta1+i >= 0) && (iphi1+j >= 0) && (ieta1+i < 48) && (iphi1+j < 24)){
+               //printf("\t \t eta1+i %d, phi1+j %d\n", ieta1+i, iphi1+j);
+               fHmpi0[iSupMod1][ieta1+i][iphi1+j]->Fill(invmass);
+             }
+             if((ieta2+i >= 0) && (iphi2+j >= 0) && (ieta2+i < 48) && (iphi2+j < 24)){
+               //printf("\t \t eta2+i %d, phi2+j %d\n", ieta2+i, iphi2+j);
+               fHmpi0[iSupMod2][ieta2+i][iphi2+j]->Fill(invmass);
+             }
+           }// j loop
+         }//i loop
+       }//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);
+      }
+      
     }
     
   } // end of loop over EMCAL clusters
   
   delete caloClustersArr;
-       
+  
   PostData(1,fOutputContainer);
-       
+  
 }
 
 //__________________________________________________
-void AliAnalysisTaskEMCALPi0CalibSelection::MaxEnergyCellPos(AliAODCaloCells* const cells, AliAODCaloCluster* const clu, Int_t& maxId)
+void AliAnalysisTaskEMCALPi0CalibSelection::MaxEnergyCellPos(AliAODCaloCells* const cells, AliAODCaloCluster* const 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 = -111;
-  
+  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;
+       
   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
-    Int_t cellAbsId = clu->GetCellAbsId(iDig);
-    Double_t eCell  = cells->GetCellAmplitude(cellAbsId)*clu->GetCellAmplitudeFraction(iDig);
-    if(eCell>eMax)  { 
-      eMax = eCell; 
+    cellAbsId = clu->GetCellAbsId(iDig);
+    fraction  = clu->GetCellAmplitudeFraction(iDig);
+    if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
+    if(IsRecalibrationOn()) {
+      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);                  
+      recalFactor = GetEMCALChannelRecalibrationFactor(imodrc,ietarc,iphirc);
+    }
+    eCell  = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
+    //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);
     }
   }
-
+  
+  //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);    
+  
+  //printf("\t Max : cell %d, iSupMod %d, ieta %d, iphi %d \n",maxId,iSupMod, ieta,iphi);
+  
 }
 
 //__________________________________________________
-void AliAnalysisTaskEMCALPi0CalibSelection::SetCalibCorrections(AliEMCALCalibData* const cdata)
-{
-  //Set new correction factors (~1) to calibration coefficients, delete previous.
+//void AliAnalysisTaskEMCALPi0CalibSelection::SetCalibCorrections(AliEMCALCalibData* const cdata)
+//{
+//  //Set new correction factors (~1) to calibration coefficients, delete previous.
+//
+//   if(fCalibData) delete fCalibData;
+//   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;
+       
+}
 
-   if(fCalibData) delete fCalibData;
-   fCalibData = cdata;
+
+//________________________________________________________________
+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 ad9ab1f..5477c28 100644 (file)
 
 // Root includes
 class TH1F;
+#include "TH2I.h"
+#include "TObjArray.h"
 
 // AliRoot includes
 #include "AliAnalysisTaskSE.h"
 class AliEMCALGeometry;
 class AliAODCaloCluster;
 class AliAODCaloCells;
-class AliEMCALCalibData ;
+//class AliEMCALCalibData ;
 #include "AliEMCALGeoParams.h"
 
 class AliAnalysisTaskEMCALPi0CalibSelection : public AliAnalysisTaskSE
@@ -32,13 +34,15 @@ public:
   // Implementation of interface methods
   virtual void UserCreateOutputObjects();
   virtual void UserExec(Option_t * opt);
-  
+  virtual void LocalInit() ;
+
   void SetClusterMinEnergy(Float_t emin) {fEmin=emin;}
   void SetClusterMaxEnergy(Float_t emax) {fEmax=emax;}
-  void SetClusterMinNCells(Float_t n)    {fMinNCells=n;}
+  void SetClusterMinNCells(Int_t n)      {fMinNCells=n;}
+  void SetNCellsGroup(Int_t n)           {fGroupNCells=n;}
 
   void SetLogWeight(Float_t weight) {fLogWeight=weight;}
-  void SetCalibCorrections(AliEMCALCalibData* const cdata);
+  //void SetCalibCorrections(AliEMCALCalibData* const cdata);
   void CreateAODFromESD();
   void CreateAODFromAOD();     
 
@@ -47,29 +51,87 @@ public:
        
   void SetGeometryName(TString name)   { fEMCALGeoName = name ; }
   TString GeometryName() const { return fEMCALGeoName ; }
-       
+  // 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);
+       
+  // 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 SetInvariantMassHistoBinRange(Int_t nBins, Float_t minbin, Float_t maxbin){
+       fNbins = nBins; fMinBin = minbin; fMaxBin = maxbin; }
+       
 private:
 
-  void MaxEnergyCellPos(AliAODCaloCells* const cells, AliAODCaloCluster* const clu, Int_t& maxId);
+  void MaxEnergyCellPos(AliAODCaloCells* const cells, AliAODCaloCluster* const clu, Int_t& iSM, Int_t& ieta, Int_t& iphi);
 
 private:
 
   AliEMCALGeometry * fEMCALGeo;  //! EMCAL geometry
-  AliEMCALCalibData* fCalibData; // corrections to CC from the previous iteration
+  //AliEMCALCalibData* fCalibData; // corrections to CC from the previous iteration
        
   Float_t fEmin;          // min. cluster energy
   Float_t fEmax;          // max. cluster energy
   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?
   TString fEMCALGeoName;  // Name of geometry to use.
-                         
+       
+  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                 
   //Output histograms  
+  Int_t   fNbins;  // N       mass bins of invariant mass histograms
+  Float_t fMinBin; // Minimum mass bins of invariant mass histograms
+  Float_t fMaxBin; // Maximum mass bins of invariant mass histograms
+
   TList*  fOutputContainer; //!histogram container
   TH1F*   fHmpi0[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows];//! two-cluster inv. mass assigned to each cell.
   TH1F*   fHmgg;            //! two-cluster inv.mass
+  TH1I*   fhNEvents;        //! Number of events counter histogram
+  TList * fCuts ;           //! List with analysis cuts
 
-  ClassDef(AliAnalysisTaskEMCALPi0CalibSelection,2);
+  ClassDef(AliAnalysisTaskEMCALPi0CalibSelection,3);
 
 };
 
diff --git a/PWG4/CaloCalib/macros/anaEMCALCalib.C b/PWG4/CaloCalib/macros/anaEMCALCalib.C
new file mode 100644 (file)
index 0000000..04c29a2
--- /dev/null
@@ -0,0 +1,446 @@
+/* $Id:  $ */
+//--------------------------------------------------
+// Example macro to do analysis with the 
+// analysis classes in PWG4PartCorr
+// Can be executed with Root and AliRoot
+//
+// Pay attention to the options and definitions
+// set in the lines below
+//
+//  Author : Gustavo Conesa Balbastre (INFN-LNF)
+//
+//-------------------------------------------------
+enum anaModes {mLocal, mLocalCAF,mPROOF,mGRID};
+//mLocal: Analyze locally files in your computer
+//mLocalCAF: Analyze locally CAF files
+//mPROOF: Analyze CAF files with PROOF
+
+//---------------------------------------------------------------------------
+//Settings to read locally several files, only for "mLocal" mode
+//The different values are default, they can be set with environmental 
+//variables: INDIR, PATTERN, NFILES, respectivelly
+char * kInDir = "/user/data/files/"; 
+char * kPattern = ""; // Data are in files kInDir/kPattern+i 
+Int_t kFile = 1; // Number of files
+//---------------------------------------------------------------------------
+//Collection file for grid analysis
+char * kXML = "collection.xml";
+//---------------------------------------------------------------------------
+
+const TString kInputData = "AOD"; //ESD, AOD, MC
+TString kTreeName = "aodTree";
+Bool_t copy = kFALSE;
+
+void anaEMCALCalib(Int_t mode=mLocal)
+{
+  // Main
+
+  //--------------------------------------------------------------------
+  // Load analysis libraries
+  // Look at the method below, 
+  // change whatever you need for your analysis case
+  // ------------------------------------------------------------------
+  LoadLibraries(mode) ;
+  //gSystem->Unload("libPWG4CaloCalib.so");
+  //Try to set the new library
+  //gSystem->Load("./PWG4CaloCalib/libPWG4CaloCalib.so");
+  //gSystem->ListLibraries();
+
+  //-------------------------------------------------------------------------------------------------
+  //Create chain from ESD and from cross sections files, look below for options.
+  //-------------------------------------------------------------------------------------------------
+  if(kInputData == "ESD") kTreeName = "esdTree" ;
+  else if(kInputData == "AOD") kTreeName = "aodTree" ;
+  else {
+    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
+    
+    //--------------------------------------
+    // 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);
+    }
+
+    //input
+    if(kInputData == "ESD")
+      {
+       // 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;
+       
+      }
+    
+    //mgr->SetDebugLevel(1);
+    
+    //-------------------------------------------------------------------------
+    //Define task, put here any other task that you want to use.
+    //-------------------------------------------------------------------------
+    AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer();
+    AliAnalysisDataContainer *coutput1 = mgr->GetCommonOutputContainer();
+    
+    gROOT->LoadMacro("AddTaskPhysicsSelection.C");
+    AliPhysicsSelectionTask* physSelTask = AddTaskPhysicsSelection();
+    
+    // ESD filter task
+    if(kInputData == "ESD" && !copy){
+      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->SetClusterMinNCells(0);
+    pi0calib->SetNCellsGroup(0);
+    pi0calib->SwitchOnBadChannelsRemoval();
+    // 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);
+
+    
+    // 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");
+    
+    AliAnalysisDataContainer *cout_cuts = mgr->CreateContainer("Cuts", TList::Class(), 
+                                                              AliAnalysisManager::kParamContainer, "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
+    //-----------------------    
+    TString smode = "";
+    if (mode==mLocal || mode == mLocalCAF) 
+      smode = "local";
+    else if (mode==mPROOF) 
+      smode = "proof";
+    else if (mode==mGRID) 
+      smode = "local";
+    
+    mgr->InitAnalysis();
+    mgr->PrintStatus();
+    mgr->StartAnalysis(smode.Data(),chain);
+
+cout <<" Analysis ended sucessfully "<< endl ;
+
+  }
+  else cout << "Chain was not produced ! "<<endl;
+  
+}
+
+void  LoadLibraries(const anaModes mode) {
+  
+  //--------------------------------------
+  // Load the needed libraries most of them already loaded by aliroot
+  //--------------------------------------
+  gSystem->Load("libTree.so");
+  gSystem->Load("libGeom.so");
+  gSystem->Load("libVMC.so");
+  gSystem->Load("libXMLIO.so");
+  
+  //----------------------------------------------------------
+  // >>>>>>>>>>> Local mode <<<<<<<<<<<<<< 
+  //----------------------------------------------------------
+  if (mode==mLocal || mode == mLocalCAF || mode == mGRID) {
+    //--------------------------------------------------------
+    // If you want to use already compiled libraries 
+    // in the aliroot distribution
+    //--------------------------------------------------------
+
+    gSystem->Load("libSTEERBase.so");
+    gSystem->Load("libESD.so");
+    gSystem->Load("libAOD.so");
+    gSystem->Load("libANALYSIS.so");
+    gSystem->Load("libANALYSISalice.so");
+    gSystem->Load("libANALYSISalice.so");
+    gSystem->Load("libPWG4CaloCalib.so");
+    
+    /*
+      //     gSystem->Load("libPWG4omega3pi.so");
+      //     gSystem->Load("libCORRFW.so");
+      //     gSystem->Load("libPWG3base.so");
+      //     gSystem->Load("libPWG3muon.so");
+ */
+    //--------------------------------------------------------
+    //If you want to use root and par files from aliroot
+    //--------------------------------------------------------  
+    /*     
+          SetupPar("STEERBase");
+          SetupPar("ESD");
+          SetupPar("AOD");
+          SetupPar("ANALYSIS");
+          SetupPar("ANALYSISalice");  
+          SetupPar("PHOSUtils");
+          SetupPar("EMCALUtils");
+          //Create Geometry
+          TGeoManager::Import("geometry.root") ; //need file "geometry.root" in local dir!!!!
+          SetupPar("PWG4CaloCalib");
+*/
+  }
+
+  //---------------------------------------------------------
+  // <<<<<<<<<< PROOF mode >>>>>>>>>>>>
+  //---------------------------------------------------------
+  else if (mode==mPROOF) {
+    //
+    // Connect to proof
+    // Put appropriate username here
+    // TProof::Reset("proof://mgheata@lxb6046.cern.ch"); 
+    TProof::Open("proof://mgheata@lxb6046.cern.ch");
+    
+    //    gProof->ClearPackages();
+    //    gProof->ClearPackage("ESD");
+    //    gProof->ClearPackage("AOD");
+    //    gProof->ClearPackage("ANALYSIS");   
+    //    gProof->ClearPackage("PWG4PartCorrBase");
+    //    gProof->ClearPackage("PWG4PartCorrDep");
+    
+    // Enable the STEERBase Package
+    gProof->UploadPackage("STEERBase.par");
+    gProof->EnablePackage("STEERBase");
+    // Enable the ESD Package
+    gProof->UploadPackage("ESD.par");
+    gProof->EnablePackage("ESD");
+    // Enable the AOD Package
+    gProof->UploadPackage("AOD.par");
+    gProof->EnablePackage("AOD");
+    // Enable the Analysis Package
+    gProof->UploadPackage("ANALYSIS.par");
+    gProof->EnablePackage("ANALYSIS");
+    // Enable the PHOS geometry Package
+    //gProof->UploadPackage("PHOSUtils.par");
+    //gProof->EnablePackage("PHOSUtils");
+    // Enable PartCorr analysis
+    gProof->UploadPackage("PWG4PartCorrBase.par");
+    gProof->EnablePackage("PWG4PartCorrBase");
+    gProof->UploadPackage("PWG4PartCorrDep.par");
+    gProof->EnablePackage("PWG4PartCorrDep");    
+    gProof->ShowEnabledPackages();
+  }  
+  
+}
+
+void SetupPar(char* pararchivename)
+{
+  //Load par files, create analysis libraries
+  //For testing, if par file already decompressed and modified
+  //classes then do not decompress.
+  TString cdir(Form("%s", gSystem->WorkingDirectory() )) ; 
+  TString parpar(Form("%s.par", pararchivename)) ; 
+  if ( gSystem->AccessPathName(parpar.Data()) ) {
+    gSystem->ChangeDirectory(gSystem->Getenv("ALICE_ROOT")) ;
+    TString processline(Form(".! make %s", parpar.Data())) ; 
+    gROOT->ProcessLine(processline.Data()) ;
+    gSystem->ChangeDirectory(cdir) ; 
+    processline = Form(".! mv $ALICE_ROOT/%s .", parpar.Data()) ;
+    gROOT->ProcessLine(processline.Data()) ;
+  } 
+  if ( gSystem->AccessPathName(pararchivename) ) {  
+    TString processline = Form(".! tar xvzf %s",parpar.Data()) ;
+    gROOT->ProcessLine(processline.Data());
+  }
+  
+  TString ocwd = gSystem->WorkingDirectory();
+  gSystem->ChangeDirectory(pararchivename);
+  
+  // check for BUILD.sh and execute
+  if (!gSystem->AccessPathName("PROOF-INF/BUILD.sh")) {
+    printf("*******************************\n");
+    printf("*** Building PAR archive    ***\n");
+    cout<<pararchivename<<endl;
+    printf("*******************************\n");
+    
+    if (gSystem->Exec("PROOF-INF/BUILD.sh")) {
+      Error("runProcess","Cannot Build the PAR Archive! - Abort!");
+      return -1;
+    }
+  }
+  // check for SETUP.C and execute
+  if (!gSystem->AccessPathName("PROOF-INF/SETUP.C")) {
+    printf("*******************************\n");
+    printf("*** Setup PAR archive       ***\n");
+    cout<<pararchivename<<endl;
+    printf("*******************************\n");
+    gROOT->Macro("PROOF-INF/SETUP.C");
+  }
+  
+  gSystem->ChangeDirectory(ocwd.Data());
+  printf("Current dir: %s\n", ocwd.Data());
+}
+
+
+
+void CreateChain(const anaModes mode, TChain * chain){
+  //Fills chain with data
+  TString ocwd = gSystem->WorkingDirectory();
+  
+  //-----------------------------------------------------------
+  //Analysis of CAF data locally and with PROOF
+  //-----------------------------------------------------------
+  if(mode ==mPROOF || mode ==mLocalCAF){
+    // Chain from CAF
+    gROOT->LoadMacro("$ALICE_ROOT/PWG0/CreateESDChain.C");
+    // The second parameter is the number of input files in the chain
+    chain = CreateESDChain("ESD12001.txt", 5);  
+  }
+  
+  //---------------------------------------
+  //Local files analysis
+  //---------------------------------------
+  else if(mode == mLocal){    
+    //If you want to add several ESD files sitting in a common directory INDIR
+    //Specify as environmental variables the directory (INDIR), the number of files 
+    //to analyze (NFILES) and the pattern name of the directories with files (PATTERN)
+
+    if(gSystem->Getenv("INDIR"))  
+      kInDir = gSystem->Getenv("INDIR") ; 
+    else cout<<"INDIR not set, use default: "<<kInDir<<endl;   
+    
+    if(gSystem->Getenv("PATTERN"))   
+      kPattern = gSystem->Getenv("PATTERN") ; 
+    else  cout<<"PATTERN not set, use default: "<<kPattern<<endl;
+    
+    if(gSystem->Getenv("NFILES"))
+      kFile = atoi(gSystem->Getenv("NFILES")) ;
+    else cout<<"NFILES not set, use default: "<<kFile<<endl;
+    
+    //Check if env variables are set and are correct
+    if ( kInDir  && kFile) {
+      printf("Get %d files from directory %s\n",kFile,kInDir);
+      if ( ! gSystem->cd(kInDir) ) {//check if ESDs directory exist
+       printf("%s does not exist\n", kInDir) ;
+       return ;
+      }
+
+  
+      cout<<"INDIR   : "<<kInDir<<endl;
+      cout<<"NFILES  : "<<kFile<<endl;
+      cout<<"PATTERN : " <<kPattern<<endl;
+
+      TString datafile="";
+      if(kInputData == "ESD") datafile = "AliESDs.root" ;
+      else if(kInputData == "AOD") datafile = "aod.root" ;
+      else if(kInputData == "MC")  datafile = "galice.root" ;
+      
+      //Loop on ESD files, add them to chain
+      Int_t event =0;
+      Int_t skipped=0 ; 
+      char file[120] ;
+      
+      for (event = 0 ; event < kFile ; event++) {
+       sprintf(file, "%s/%s%d/%s", kInDir,kPattern,event,datafile.Data()) ; 
+       TFile * fESD = 0 ; 
+       //Check if file exists and add it, if not skip it
+       if ( fESD = TFile::Open(file)) {
+         if ( fESD->Get(kTreeName) ) { 
+           printf("++++ Adding %s\n", file) ;
+           chain->AddFile(file);
+         }
+       }
+       else { 
+         printf("---- Skipping %s\n", file) ;
+         skipped++ ;
+       }
+      }
+      printf("number of entries # %lld, skipped %d\n", chain->GetEntries(), skipped*100) ;     
+    }
+    else {
+      TString input = "AliESDs.root" ;
+      cout<<">>>>>> No list added, take a single file <<<<<<<<< "<<input<<endl;
+      chain->AddFile(input);
+    }
+    
+  }// local files analysis
+  
+  //------------------------------
+  //GRID xml files
+  //-----------------------------
+  else if(mode == mGRID){
+    //Get colection file. It is specified by the environmental
+    //variable XML
+
+    if(gSystem->Getenv("XML") )
+      kXML = gSystem->Getenv("XML");
+    else
+      sprintf(kXML, "collection.xml") ; 
+    
+    if (!TFile::Open(kXML)) {
+      printf("No collection file with name -- %s -- was found\n",kXML);
+      return ;
+    }
+    else cout<<"XML file "<<kXML<<endl;
+
+    //Load necessary libraries and connect to the GRID
+    gSystem->Load("libNetx.so") ; 
+    gSystem->Load("libRAliEn.so"); 
+    TGrid::Connect("alien://") ;
+
+    //Feed Grid with collection file
+    //TGridCollection * collection =  (TGridCollection*)gROOT->ProcessLine(Form("TAlienCollection::Open(\"%s\", 0)", kXML));
+    TGridCollection * collection = (TGridCollection*) TAlienCollection::Open(kXML);
+    if (! collection) {
+      AliError(Form("%s not found", kXML)) ; 
+      return kFALSE ; 
+    }
+    TGridResult* result = collection->GetGridResult("",0 ,0);
+   
+    // Makes the ESD chain 
+    printf("*** Getting the Chain       ***\n");
+    for (Int_t index = 0; index < result->GetEntries(); index++) {
+      TString alienURL = result->GetKey(index, "turl") ; 
+      cout << "================== " << alienURL << endl ; 
+      chain->Add(alienURL) ; 
+    }
+  }// xml analysis
+  
+  gSystem->ChangeDirectory(ocwd.Data());
+}
+
index b08351e..d799937 100755 (executable)
@@ -2,7 +2,9 @@ void SETUP()
 {
 
    // Load the ESD library
-   gSystem->Load("libPWG4CaloCalib");
+   //gSystem->Load("libPWG4CaloCalib");
+   TString ocwd = gSystem->WorkingDirectory();
+   gSystem->Load(ocwd+"/libPWG4CaloCalib.so");
 
    // Set the Include paths
    gSystem->SetIncludePath("-I$ROOTSYS/include -IPWG4CaloCalib");
index c91b29d..9b36c37 100755 (executable)
@@ -2,7 +2,9 @@ void SETUP()
 {
 
    // Load the ESD library
-   gSystem->Load("libPWG4PartCorrBase");
+   //gSystem->Load("libPWG4PartCorrBase");
+   TString ocwd = gSystem->WorkingDirectory();
+   gSystem->Load(ocwd+"/libPWG4PartCorrBase.so");
 
    // Set the Include paths
    gSystem->SetIncludePath("-I$ROOTSYS/include -IPWG4PartCorrBase");
index fbc2d01..6b3886e 100755 (executable)
@@ -2,7 +2,9 @@ void SETUP()
 {
 
    // Load the ESD library
-   gSystem->Load("libPWG4PartCorrDep");
+   //gSystem->Load("libPWG4PartCorrDep");
+   TString ocwd = gSystem->WorkingDirectory();
+   gSystem->Load(ocwd+"/libPWG4PartCorrDep.so");
 
    // Set the Include paths
    gSystem->SetIncludePath("-I$ROOTSYS/include -IPWG4PartCorrDep");