new task for reclusterization in EMCAL
authorgconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 5 Jan 2011 14:05:48 +0000 (14:05 +0000)
committergconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 5 Jan 2011 14:05:48 +0000 (14:05 +0000)
PWG4/CaloCalib/AliAnalysisTaskEMCALClusterize.cxx [new file with mode: 0644]
PWG4/CaloCalib/macros/emcalReclusterize.C [new file with mode: 0644]
PWG4/PWG4CaloCalibLinkDef.h
PWG4/libPWG4CaloCalib.pkg

diff --git a/PWG4/CaloCalib/AliAnalysisTaskEMCALClusterize.cxx b/PWG4/CaloCalib/AliAnalysisTaskEMCALClusterize.cxx
new file mode 100644 (file)
index 0000000..2c35fc0
--- /dev/null
@@ -0,0 +1,490 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+/* $Id: $ */
+
+//_________________________________________________________________________
+// This analysis provides a new list of clusters to be used in other analysis
+//
+// Author: Gustavo Conesa Balbastre,
+//         Adapted from analysis class from Deepa Thomas
+//
+//
+//_________________________________________________________________________
+
+// --- Root ---
+#include "TString.h"
+#include "TRefArray.h"
+#include "TClonesArray.h"
+#include "TTree.h"
+#include "TGeoManager.h"
+
+// --- AliRoot Analysis Steering
+#include "AliAnalysisTask.h"
+#include "AliAnalysisManager.h"
+#include "AliESDEvent.h"
+#include "AliGeomManager.h"
+#include "AliVCaloCells.h"
+#include "AliAODCaloCluster.h"
+#include "AliCDBManager.h"
+#include "AliCDBStorage.h"
+#include "AliCDBEntry.h"
+#include "AliLog.h"
+#include "AliVEventHandler.h"
+
+// --- EMCAL
+#include "AliEMCALRecParam.h"
+#include "AliEMCALAfterBurnerUF.h"
+#include "AliEMCALGeometry.h"
+#include "AliEMCALClusterizerNxN.h"
+#include "AliEMCALClusterizerv1.h"
+#include "AliEMCALRecPoint.h"
+#include "AliEMCALDigit.h"
+#include "AliCaloCalibPedestal.h"
+#include "AliEMCALCalibData.h"
+
+#include "AliAnalysisTaskEMCALClusterize.h"
+
+ClassImp(AliAnalysisTaskEMCALClusterize)
+
+//________________________________________________________________________
+AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize(const char *name) 
+  : AliAnalysisTaskSE(name)
+  , fGeom(0), fGeomName("EMCAL_FIRSTYEARV1"),  fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
+  , fCalibData(0),       fPedestalData(0),     fOCDBpath("raw://")
+  , fDigitsArr(0),       fClusterArr(0),       fCaloClusterArr(0)
+  , fRecParam(0),        fClusterizer(0),      fUnfolder(0),           fJustUnfold(kFALSE) 
+  , fOutputAODBranch(0), fOutputAODBranchName("newEMCALClusters"),     fFillAODFile(kTRUE)
+  
+  {
+  //ctor
+  for(Int_t i = 0; i < 10; i++) fGeomMatrix[i] = 0 ;
+  fDigitsArr       = new TClonesArray("AliEMCALDigit",200);
+  fClusterArr      = new TObjArray(100);
+  fCaloClusterArr  = new TObjArray(100);
+  fRecParam        = new AliEMCALRecParam;
+}
+
+//________________________________________________________________________
+AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize() 
+  : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskEMCALClusterize")
+  , fGeom(0), fGeomName("EMCAL_FIRSTYEARV1"),  fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
+  , fCalibData(0),       fPedestalData(0),     fOCDBpath("raw://")
+  , fDigitsArr(0),       fClusterArr(0),       fCaloClusterArr(0)
+  , fRecParam(0),        fClusterizer(0),      fUnfolder(0),           fJustUnfold(kFALSE) 
+  , fOutputAODBranch(0), fOutputAODBranchName("newEMCALClusters"),     fFillAODFile(kFALSE)
+{
+  // Constructor
+  for(Int_t i = 0; i < 10; i++) fGeomMatrix[i] = 0 ;
+  fDigitsArr       = new TClonesArray("AliEMCALDigit",200);
+  fClusterArr      = new TObjArray(100);
+  fCaloClusterArr  = new TObjArray(100);
+  fRecParam        = new AliEMCALRecParam;
+}
+
+//________________________________________________________________________
+AliAnalysisTaskEMCALClusterize::~AliAnalysisTaskEMCALClusterize()
+{
+  //dtor 
+  
+  if (fDigitsArr){
+    fDigitsArr->Clear("C");
+    delete fDigitsArr; fDigitsArr=0;
+  }
+  
+  if (fClusterArr){
+    fClusterArr->Delete();
+    delete fClusterArr; fClusterArr=0;
+  }
+  
+  if (fCaloClusterArr){
+    fCaloClusterArr->Delete();
+    delete fCaloClusterArr; fCaloClusterArr=0;
+  }
+
+  if(fClusterizer) {delete fClusterizer; fClusterizer = 0;}
+  if(fGeom)        {delete fGeom;        fGeom        = 0;}
+  if(fUnfolder)    {delete fUnfolder;    fUnfolder    = 0;}
+
+}
+
+//-------------------------------------------------------------------
+void AliAnalysisTaskEMCALClusterize::UserCreateOutputObjects()
+{
+  // Init geometry, create list of output clusters
+  
+  fGeom =  AliEMCALGeometry::GetInstance(fGeomName) ;  
+
+  fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
+  fOutputAODBranch->SetName(fOutputAODBranchName);
+  AddAODBranch("TClonesArray", &fOutputAODBranch);
+
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskEMCALClusterize::UserExec(Option_t *) 
+{
+  // Main loop
+  // Called for each event
+  
+  //printf("___ Event __ %d __\n",(Int_t)Entry());
+  
+  //Remove the contents of output list set in the previous event 
+  fOutputAODBranch->Clear("C");
+  
+  AliVEvent * event = InputEvent();
+  if (!event) {
+    printf("ERROR: event not available\n");
+    return;
+  }
+  
+  //Magic line to write events to AOD file
+  AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(fFillAODFile);
+  
+  //-------------------------------------------------------------------------------------
+  //Set the geometry matrix, for the first event, skip the rest
+  //-------------------------------------------------------------------------------------
+  if(!fGeomMatrixSet){
+    if(fLoadGeomMatrices){
+      for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
+        if(fGeomMatrix[mod]){
+          if(DebugLevel() > 1) 
+            fGeomMatrix[mod]->Print();
+          fGeom->SetMisalMatrix(fGeomMatrix[mod],mod) ;  
+        }
+        fGeomMatrixSet=kTRUE;
+      }//SM loop
+    }//Load matrices
+    else if(!gGeoManager){
+      printf("AliAnalysisTaskEMCALClusterize::UserExec() - Get geo matrices from data\n");
+      //Still not implemented in AOD, just a workaround to be able to work at least with ESDs  
+      if(!strcmp(event->GetName(),"AliAODEvent")) {
+        if(DebugLevel() > 1) 
+          printf("AliAnalysisTaskCaloFilter Use ideal geometry, values geometry matrix not kept in AODs.\n");
+      }//AOD
+      else {   
+        if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALClusterize Load Misaligned matrices. \n");
+        AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event) ;
+        if(!esd) {
+          printf("AliAnalysisTaskCaloFilter::UserExec() - This event does not contain ESDs?");
+          return;
+        }
+        for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
+          if(DebugLevel() > 1) 
+            esd->GetEMCALMatrix(mod)->Print();
+          if(esd->GetEMCALMatrix(mod)) fGeom->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
+        } 
+        fGeomMatrixSet=kTRUE;
+      }//ESD
+    }//Load matrices from Data 
+  }//first event
+
+  //Get the list of cells needed for unfolding and reclustering
+  AliVCaloCells *cells= event->GetEMCALCells();
+  Int_t nClustersOrg = 0;
+  //-------------------------------------------
+  //---------Unfolding clusters----------------
+  //-------------------------------------------
+  if (fJustUnfold) {
+    
+    //Fill the array with the EMCAL clusters, copy them
+    for (Int_t i = 0; i < event->GetNumberOfCaloClusters(); i++)
+    {
+      AliVCluster *clus = event->GetCaloCluster(i);
+      if(clus->IsEMCAL()){        
+        //printf("Org Cluster %d, E %f \n",i, clus->E());
+        if(dynamic_cast<AliESDCaloCluster*> (clus)){
+          fCaloClusterArr->Add( new AliESDCaloCluster(*(dynamic_cast<AliESDCaloCluster*> (clus))) );   
+        }//ESD
+        else{
+          fCaloClusterArr->Add( new AliAODCaloCluster(*(dynamic_cast<AliAODCaloCluster*> (clus))) );   
+        }//AOD
+        nClustersOrg++;
+      }
+    }
+    
+    //Do the unfolding
+    fUnfolder->UnfoldClusters(fCaloClusterArr, cells);
+    
+    //CLEAN-UP
+    fUnfolder->Clear();
+    
+    //printf("nClustersOrg %d\n",nClustersOrg);
+  }
+  //-------------------------------------------
+  //---------- Recluster cells ----------------
+  //-------------------------------------------
+  
+  else{
+    //-------------------------------------------------------------------------------------
+    //Transform CaloCells into Digits
+    //-------------------------------------------------------------------------------------
+    Int_t idigit =  0;
+    Int_t id     = -1;
+    Float_t amp  = -1; 
+    Float_t time = -1; 
+    
+    Double_t cellAmplitude = 0;
+    Double_t cellTime      = 0;
+    Short_t  cellNumber    = 0;
+        
+    TTree *digitsTree = new TTree("digitstree","digitstree");
+    digitsTree->Branch("EMCAL","TClonesArray", &fDigitsArr, 32000);
+    
+    for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
+    {
+      if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
+        break;
+      
+      time = cellTime;
+      amp  = cellAmplitude;
+      id   = cellNumber;
+      
+      AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->New(idigit);
+      digit->SetId(id);
+      digit->SetAmplitude(amp);
+      digit->SetTime(time);
+      digit->SetTimeR(time);
+      digit->SetIndexInList(idigit);
+      digit->SetType(AliEMCALDigit::kHG);
+      //if(Entry()==0) printf("Digit: Id %d, amp %f, time %e, index %d\n",id, amp,time,idigit);
+      idigit++;
+    }
+    
+    //Fill the tree with digits
+    digitsTree->Fill();
+    
+    //-------------------------------------------------------------------------------------
+    //Do the clusterization
+    //-------------------------------------------------------------------------------------        
+    TTree *clustersTree = new TTree("clustertree","clustertree");
+    
+    fClusterizer->SetInput(digitsTree);
+    fClusterizer->SetOutput(clustersTree);
+    fClusterizer->Digits2Clusters("");
+    
+    //-------------------------------------------------------------------------------------
+    //Transform the recpoints into AliVClusters
+    //-------------------------------------------------------------------------------------
+    
+    clustersTree->SetBranchStatus("*",0); //disable all branches
+    clustersTree->SetBranchStatus("EMCALECARP",1); //Enable only the branch we need
+    
+    TBranch *branch = clustersTree->GetBranch("EMCALECARP");
+    branch->SetAddress(&fClusterArr);
+    branch->GetEntry(0);
+    
+    RecPoints2Clusters(fDigitsArr, fClusterArr, fCaloClusterArr);
+    
+    //---CLEAN UP-----
+    fClusterizer->Clear();
+    fDigitsArr  ->Clear("C");
+    fClusterArr ->Delete(); // Do not Clear(), it leaks, why?
+
+    clustersTree->Delete("all");
+    digitsTree  ->Delete("all");
+    
+  }
+  
+  //-------------------------------------------------------------------------------------
+  //Put the new clusters in the AOD list
+  //-------------------------------------------------------------------------------------
+  
+  Int_t kNumberOfCaloClusters   = fCaloClusterArr->GetEntries();
+  //printf("New clusters %d\n",kNumberOfCaloClusters);
+  //if(nClustersOrg!=kNumberOfCaloClusters) printf("Different number: Org %d, New %d\n",nClustersOrg,kNumberOfCaloClusters);
+  for(Int_t i = 0; i < kNumberOfCaloClusters; i++){
+    AliAODCaloCluster *newCluster = (AliAODCaloCluster *) fCaloClusterArr->At(i);
+    //if(Entry()==0)printf("newCluster E %f\n", newCluster->E());
+    new((*fOutputAODBranch)[i])  AliAODCaloCluster(*newCluster);
+  }
+  
+  //---CLEAN UP-----
+  fCaloClusterArr->Delete(); // Do not Clear(), it leaks, why?
+  
+ }      
+
+
+//_____________________________________________________________________
+Bool_t AliAnalysisTaskEMCALClusterize::UserNotify()
+{
+  //Access to OCDB stuff
+  if(DebugLevel() > 1 )printf(" AliAnalysisTaskEMCALClusterize::UserNotify() - Begin \n");
+  AliVEvent * event = InputEvent();
+  if (event)
+  {
+    fGeom = AliEMCALGeometry::GetInstance(fGeomName);
+    AliCDBManager *cdb = AliCDBManager::Instance();
+    cdb->SetDefaultStorage(fOCDBpath.Data());
+    cdb->SetRun(event->GetRunNumber());
+    //
+    // EMCAL from RAW OCDB
+    if (fOCDBpath.Contains("alien:"))
+    {
+      cdb->SetSpecificStorage("EMCAL/Calib/Data","alien://Folder=/alice/data/2010/OCDB");
+      cdb->SetSpecificStorage("EMCAL/Calib/Pedestals","alien://Folder=/alice/data/2010/OCDB");
+    }
+    TString path = cdb->GetDefaultStorage()->GetBaseFolder();
+    
+    // init parameters:
+    //Get calibration parameters       
+    if(!fCalibData)
+    {
+      AliCDBEntry *entry = (AliCDBEntry*) 
+           AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
+      if (entry) fCalibData =  (AliEMCALCalibData*) entry->GetObject();
+    }
+    
+    if(!fCalibData)
+      AliFatal("Calibration parameters not found in CDB!");
+    
+    //Get calibration parameters       
+    if(!fPedestalData)
+    {
+      AliCDBEntry *entry = (AliCDBEntry*) 
+           AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals");
+      if (entry) fPedestalData =  (AliCaloCalibPedestal*) entry->GetObject();
+    }
+    
+    if(!fPedestalData)
+      AliFatal("Dead map not found in CDB!");
+    
+    //      cout << "[i] Change of run number: " << fAOD->GetRunNumber() << endl;
+    InitClusterization();
+    
+  }
+  else
+  {
+    printf(" AliAnalysisTaskEMCALClusterize::UserNotify() - Event not available!!! \n");
+  }
+
+  return kTRUE;
+  
+}
+
+//________________________________________________________________________________________
+void AliAnalysisTaskEMCALClusterize::InitClusterization()
+{
+  //Select clusterization/unfolding algorithm and set all the needed parameters
+  
+  if(!fJustUnfold){
+    
+    //First init the clusterizer
+    if     (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
+      fClusterizer = new AliEMCALClusterizerv1 (fGeom, fCalibData, fPedestalData);
+    else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN)
+      fClusterizer = new AliEMCALClusterizerNxN(fGeom, fCalibData, fPedestalData);
+    else{
+      printf("AliAnalysisTaskEMCALClusterize::InitClusterization() - Clusterizer < %d > not available",
+             fRecParam->GetClusterizerFlag());
+      abort();
+    }
+    
+    //Now set the parameters
+    fClusterizer->SetECAClusteringThreshold( fRecParam->GetClusteringThreshold() );
+    fClusterizer->SetECALogWeight          ( fRecParam->GetW0()                  );
+    fClusterizer->SetMinECut               ( fRecParam->GetMinECut()             );    
+    fClusterizer->SetUnfolding             ( fRecParam->GetUnfold()              );
+    fClusterizer->SetECALocalMaxCut        ( fRecParam->GetLocMaxCut()           );
+    fClusterizer->SetTimeCut               ( fRecParam->GetTimeCut()             );
+    fClusterizer->SetTimeMin               ( fRecParam->GetTimeMin()             );
+    fClusterizer->SetTimeMax               ( fRecParam->GetTimeMax()             );
+    fClusterizer->SetInputCalibrated       ( kTRUE                               );
+    
+    //In case of unfolding after clusterization is requested, set the corresponding parameters
+    if(fRecParam->GetUnfold()){
+      
+      Int_t i=0;
+      for (i = 0; i < 8; i++) {
+        fClusterizer->SetSSPars(i, fRecParam->GetSSPars(i));
+      }//end of loop over parameters
+      for (i = 0; i < 3; i++) {
+        fClusterizer->SetPar5  (i, fRecParam->GetPar5(i));
+        fClusterizer->SetPar6  (i, fRecParam->GetPar6(i));
+      }//end of loop over parameters
+      
+      fClusterizer->InitClusterUnfolding();
+            
+    }// to unfold
+    
+    
+  }else{
+    //Now init the unfolding afterburner 
+    fUnfolder =  new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut());
+ }
+}
+//________________________________________________________________________________________
+void AliAnalysisTaskEMCALClusterize::RecPoints2Clusters(TClonesArray *digitsArr, TObjArray *recPoints, TObjArray *clusArray)
+{
+  // Restore clusters from recPoints
+  // Cluster energy, global position, cells and their amplitude fractions are restored
+  
+  for(Int_t i = 0; i < recPoints->GetEntriesFast(); i++)
+  {
+    AliEMCALRecPoint *recPoint = (AliEMCALRecPoint*) recPoints->At(i);
+    
+    const Int_t ncells = recPoint->GetMultiplicity();
+    Int_t ncells_true = 0;
+    
+    // cells and their amplitude fractions
+    UShort_t   absIds[ncells];  
+    Double32_t ratios[ncells];
+    
+    for (Int_t c = 0; c < ncells; c++) {
+      AliEMCALDigit *digit = (AliEMCALDigit*) digitsArr->At(recPoint->GetDigitsList()[c]);
+      
+      absIds[ncells_true] = digit->GetId();
+      ratios[ncells_true] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude();
+      
+      if (ratios[ncells_true] > 0.001) ncells_true++;
+    }
+    
+    if (ncells_true < 1) {
+      Warning("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters", "skipping cluster with no cells");
+      continue;
+    }
+    
+    TVector3 gpos;
+    Float_t g[3];
+    
+    // calculate new cluster position
+    recPoint->EvalGlobalPosition(fRecParam->GetW0(), digitsArr);
+    recPoint->GetGlobalPosition(gpos);
+    gpos.GetXYZ(g);
+    
+    // create a new cluster
+    AliAODCaloCluster *clus = new AliAODCaloCluster();
+    clus->SetType(AliVCluster::kEMCALClusterv1);
+    clus->SetE(recPoint->GetEnergy());
+    clus->SetPosition(g);
+    clus->SetNCells(ncells_true);
+    clus->SetCellsAbsId(absIds);
+    clus->SetCellsAmplitudeFraction(ratios);
+    clus->SetDispersion(recPoint->GetDispersion());
+    clus->SetChi2(-1); //not yet implemented
+    clus->SetTOF(recPoint->GetTime()) ; //time-of-fligh
+    clus->SetNExMax(recPoint->GetNExMax()); //number of local maxima
+    Float_t elipAxis[2];
+    recPoint->GetElipsAxis(elipAxis);
+    clus->SetM02(elipAxis[0]*elipAxis[0]) ;
+    clus->SetM20(elipAxis[1]*elipAxis[1]) ;
+
+    clusArray->Add(clus);
+
+  } // recPoints loop
+  
+}
+
+
+
diff --git a/PWG4/CaloCalib/macros/emcalReclusterize.C b/PWG4/CaloCalib/macros/emcalReclusterize.C
new file mode 100644 (file)
index 0000000..781b23b
--- /dev/null
@@ -0,0 +1,467 @@
+/* $Id:  $ */
+//--------------------------------------------------
+// Example macro to do Calorimeters filtering
+// copy ESDs into AODs
+//
+// 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 = "ESD"; //ESD, AOD, MC
+TString kTreeName = "esdTree";
+Bool_t kUsePhysSel = kFALSE;
+
+void emcalReclusterize(Int_t mode=mLocal)
+{
+  // Main
+  char cmd[200] ; 
+  sprintf(cmd, ".! rm -rf aod.root") ; 
+  gROOT->ProcessLine(cmd) ; 
+  //--------------------------------------------------------------------
+  // 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
+    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();
+    
+    // ESD physics selection task
+    if(kInputData == "ESD" && kUsePhysSel){
+      gROOT->LoadMacro("AddTaskPhysicsSelection.C");
+      AliPhysicsSelectionTask* physSelTask = AddTaskPhysicsSelection();
+    }
+    
+    // Create containers for input/output
+    AliAnalysisDataContainer *cinput1  = mgr->GetCommonInputContainer();
+    AliAnalysisDataContainer *coutput1 = mgr->GetCommonOutputContainer();
+    AliAnalysisTaskEMCALClusterize * clusterize = new AliAnalysisTaskEMCALClusterize("EMCALClusterize");
+    if(kUsePhysSel)clusterize->SelectCollisionCandidates();
+    //clusterize->SetOCDBPath("local://$ALICE_ROOT/OCDB");
+    clusterize->FillAODFile(kFALSE); // fill aod.root with clusters?, not really needed for analysis.
+    clusterize->JustUnfold(kTRUE); // if TRUE, do just unfolding, do not recluster cells
+    AliEMCALRecParam * params = clusterize->GetRecParam();
+    params->SetClusterizerFlag(AliEMCALRecParam::kClusterizerNxN);
+    params->SetClusteringThreshold(0.1); // 100 MeV                                             
+    params->SetMinECut(0.01);  //10 MeV    
+    params->SetUnfold(kFALSE);
+    params->SetW0(4.5);
+    params->SetTimeCut(-1);///Open this cut for AODs
+    params->SetTimeMin(-1);//Open this cut for AODs
+    params->SetTimeMax(1e6);//Open this cut for AODs
+
+//    TGeoHMatrix *matrix[4];
+    
+//     double rotationMatrix[4][9] = {-0.014585, -0.999892, -0.002031, 0.999892, -0.014589,  0.001950, -0.001979, -0.002003,  0.999996,
+//                                -0.014585,  0.999892,  0.002031, 0.999892,  0.014589, -0.001950, -0.001979,  0.002003, -0.999996,
+//                                -0.345861, -0.938280, -0.003412, 0.938281, -0.345869,  0.001950, -0.003010, -0.002527,  0.999992,
+//                                -0.345861,  0.938280,  0.003412, 0.938281,  0.345869, -0.001950, -0.003010,  0.002527, -0.999992};
+    
+//     double translationMatrix[4][3] = {0.367264,    446.508738,  175.97185+0.3,
+//                                   1.078181,    445.826258, -174.026758+0.3,
+//                                   -153.843916, 418.304256,  175.956905+0.8,
+//                                   -152.649580, 417.621779, -174.040392+0.8};
+
+
+//     double rotationMatrix[4][9] = {-0.014587, -0.999892, -0.002031, 0.999892, -0.014591,  0.001979, -0.002009, -0.002002,  0.999996,
+//                              -0.014587,  0.999892,  0.002031, 0.999892,  0.014591, -0.001979, -0.002009,  0.002002, -0.999996,
+//                              -0.345864, -0.938278, -0.003412, 0.938276, -0.345874,  0.003010, -0.004004, -0.002161,  0.999990,
+//                              -0.345864,  0.938278,  0.003412, 0.938276,  0.345874, -0.003010, -0.004004,  0.002161, -0.999990};
+    
+//     double translationMatrix[4][3] = {0.351659,    447.576446,  176.269742,
+//                                   1.062577,    446.893974, -173.728870,
+//                                   -154.213287, 419.306156,  176.753692,
+//                                   -153.018950, 418.623681, -173.243605};
+
+//     for(int j=0; j<4; j++)
+//       {
+//     matrix[j] = new TGeoHMatrix();
+//     matrix[j]->SetRotation(rotationMatrix[j]);
+//     matrix[j]->SetTranslation(translationMatrix[j]);
+//     matrix[j]->Print();
+//     clusterize->SetGeometryMatrixInSM(matrix[j],j);
+//       }
+    
+    
+//     clusterize->SwitchOnLoadOwnEMCALGeometryMatrices();
+    
+        
+    mgr->AddTask(clusterize);
+    mgr->ConnectInput  (clusterize,  0, cinput1);
+    mgr->ConnectOutput (clusterize, 0, coutput1 );
+   
+
+    //-----------------------
+    // 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");
+  gSystem->Load("libMatrix.so");
+  gSystem->Load("libPhysics.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("libEMCALUtils.so");  
+
+    //SetupPar("EMCALUtils"); 
+    //SetupPar("PWG4CaloCalib"); 
+    
+    TGeoManager::Import("geometry.root") ; //need file "geometry.root" in local dir!!!!
+    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 = "AliAODs.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 d56bb3d2225f435541e83c389f63e0f0d59b9a52..4034b4e3d828f0a52f668c3fdabe0a6797d0fa33 100644 (file)
@@ -7,5 +7,6 @@
 #pragma link C++ class AliAnalysisTaskCaloFilter+;
 #pragma link C++ class AliAnalysisTaskPHOSPi0CalibSelection+;
 #pragma link C++ class AliAnalysisTaskEMCALPi0CalibSelection+;
+#pragma link C++ class AliAnalysisTaskEMCALClusterize+;
 
 #endif
index 33211d80e93f9179ab096aaab88734c0e419c9c0..12d2bdb749d411b5e4c706ce3e701a0a1331ccc3 100644 (file)
@@ -1,7 +1,6 @@
 #-*- Mode: Makefile -*-
 
-SRCS =  CaloCalib/AliAnalysisTaskCaloFilter.cxx CaloCalib/AliAnalysisTaskPHOSPi0CalibSelection.cxx CaloCalib/AliAnalysisTaskEMCALPi0CalibSelection.cxx 
-      
+SRCS =  CaloCalib/AliAnalysisTaskCaloFilter.cxx CaloCalib/AliAnalysisTaskPHOSPi0CalibSelection.cxx CaloCalib/AliAnalysisTaskEMCALPi0CalibSelection.cxx CaloCalib/AliAnalysisTaskEMCALClusterize.cxx 
 HDRS:= $(SRCS:.cxx=.h) 
 
 DHDR:= PWG4CaloCalibLinkDef.h
@@ -12,6 +11,6 @@ EINCLUDE= PHOS EMCAL PWG4/CaloCalib
 
 ifeq (win32gcc,$(ALICE_TARGET))
 PACKSOFLAGS:= $(SOFLAGS) -L$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET) -lSTEERBase \
-                         -lESD -lAOD -lANALYSIS -lANALYSISalice -lPHOSbase -lEMCALbase \
+                         -lESD -lAOD -lANALYSIS -lANALYSISalice -lPHOSbase -lEMCALUtils -lEMCALbase -lEMCALrec\
                          -L$(ROOTLIBDIR) -lEG
 endif