--- /dev/null
+/**************************************************************************
+ * 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
+
+}
+
+
+
--- /dev/null
+/* $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());
+}
+