PWGGA/EMCALTasks/macros/emcalReclusterize.C

   1 // $Id$
   2 //--------------------------------------------------
   3 // Example macro to do Calorimeters filtering
   4 // copy ESDs into AODs
   5 //
   6 // Pay attention to the options and definitions
   7 // set in the lines below
   8 //
   9 //  Author : Gustavo Conesa Balbastre (INFN-LNF)
  10 //
  11 //-------------------------------------------------
  12 enum anaModes {mLocal, mLocalCAF,mPROOF,mGRID};
  13 //mLocal: Analyze locally files in your computer
  14 //mLocalCAF: Analyze locally CAF files
  15 //mPROOF: Analyze CAF files with PROOF
  16
  17 //---------------------------------------------------------------------------
  18 //Settings to read locally several files, only for "mLocal" mode
  19 //The different values are default, they can be set with environmental
  20 //variables: INDIR, PATTERN, NFILES, respectivelly
  21 char * kInDir = "/Users/Gustavo/Work/data/134908/pass1";
  22 //char * kInDir = "/Users/Gustavo/Work/data/137366/";
  23 char * kPattern = ""; // Data are in files kInDir/kPattern+i
  24 Int_t kFile = 1; // Number of files
  25 //---------------------------------------------------------------------------
  26 //Collection file for grid analysis
  27 char * kXML = "collection.xml";
  28 //---------------------------------------------------------------------------
  29
  30 TString kInputData = "AOD"; //ESD, AOD, MC
  31 TString kTreeName = "esdTree";
  32 Bool_t kUsePhysSel = kFALSE;
  33 Bool_t kEmbed = kTRUE;
  34
  35 void emcalReclusterize(Int_t mode=mLocal)
  36 {
  37   // Main
  38   //char cmd[200] ;
  39   //sprintf(cmd, ".! rm -rf outputAOD.root") ;
  40   //gROOT->ProcessLine(cmd) ;
  41   //--------------------------------------------------------------------
  42   // Load analysis libraries
  43   // Look at the method below,
  44   // change whatever you need for your analysis case
  45   // ------------------------------------------------------------------
  46   LoadLibraries(mode) ;
  47
  48   //-------------------------------------------------------------------------------------------------
  49   //Create chain from ESD and from cross sections files, look below for options.
  50   //-------------------------------------------------------------------------------------------------
  51   if(kInputData == "ESD")      kTreeName = "esdTree" ;
  52   else if(kInputData == "AOD") kTreeName = "aodTree" ;
  53   else {
  54     cout<<"Wrong  data type "<<kInputData<<endl;
  55     break;
  56   }
  57
  58   TChain *chain       = new TChain(kTreeName) ;
  59   CreateChain(mode, chain);
  60
  61   if(chain){
  62     AliLog::SetGlobalLogLevel(AliLog::kError);//Minimum prints on screen
  63
  64     //--------------------------------------
  65     // Make the analysis manager
  66     //-------------------------------------
  67     AliAnalysisManager *mgr  = new AliAnalysisManager("Manager", "Manager");
  68
  69     // AOD output handler
  70     AliAODHandler* aodoutHandler   = new AliAODHandler();
  71     aodoutHandler->SetOutputFileName("outputAOD.root");
  72     if(kEmbed){
  73       aodoutHandler->SetCreateNonStandardAOD();
  74       kInputData = "AOD";
  75     }
  76     mgr->SetOutputEventHandler(aodoutHandler);
  77
  78
  79     //input
  80     if(kInputData == "ESD")
  81     {
  82       // ESD handler
  83       AliESDInputHandler *esdHandler = new AliESDInputHandler();
  84       mgr->SetInputEventHandler(esdHandler);
  85       esdHandler->SetReadFriends(kFALSE);
  86       cout<<"ESD handler "<<mgr->GetInputEventHandler()<<endl;
  87     }
  88     if(kInputData == "AOD")
  89     {
  90       // AOD handler
  91       AliAODInputHandler *aodHandler = new AliAODInputHandler();
  92       mgr->SetInputEventHandler(aodHandler);
  93       if(kEmbed){
  94         aodHandler->SetMergeEvents(kTRUE);
  95         aodHandler->AddFriend("AliAOD.root");
  96       }
  97
  98       cout<<"AOD handler "<<mgr->GetInputEventHandler()<<endl;
  99
 100     }
 101
 102     //mgr->SetDebugLevel(1);
 103
 104     //-------------------------------------------------------------------------
 105     //Define task, put here any other task that you want to use.
 106     //-------------------------------------------------------------------------
 107     AliAnalysisDataContainer *cinput1  = mgr->GetCommonInputContainer();
 108     AliAnalysisDataContainer *coutput1 = mgr->GetCommonOutputContainer();
 109
 110     // ESD physics selection task
 111     if(kInputData == "ESD" && kUsePhysSel){
 112       gROOT->LoadMacro("AddTaskPhysicsSelection.C");
 113       AliPhysicsSelectionTask* physSelTask = AddTaskPhysicsSelection();
 114     }
 115
 116     gROOT->LoadMacro("$ALICE_ROOT/PWGGA/EMCALTasks/macros/AddTaskEMCALClusterize.C");
 117
 118     Bool_t kMC = kFALSE; // extra checks in case of MC analysis
 119     Bool_t outAOD = kFALSE ; // Generate output AOD with new clusters
 120     Bool_t calibEE = kTRUE; // It is set automatically, but here we force to use it or not in any case
 121     Bool_t calibTT = kTRUE; // It is set automatically, but here we force to use it or not in any case
 122     if(kRun < 122195 || (kRun > 126437 && kRun < 136851) || kMC) calibTT=kFALSE ; // Recalibration parameters not available for LHC10a,b,c,e,f,g
 123     Bool_t badMap  = kTRUE; // It is set automatically, but here we force to use it or not in any case
 124     Bool_t  clTM      = kTRUE;
 125     Bool_t  exo       = kTRUE;  // Remove exotic cells
 126     Bool_t  clnonlin  = kTRUE;  // Apply non linearity (clusterization)
 127     Int_t   minEcell  = 50;     // 50  MeV (10 MeV used in reconstruction)
 128     Int_t   minEseed  = 100;    // 100 MeV
 129     Int_t   dTime     = 0;      // default, 250 ns
 130     Int_t   wTime     = 30;     // default 425 < T < 825 ns, careful if time calibration is on
 131     Int_t   unfMinE   = 15;     // Remove cells with less than 15 MeV from cluster after unfolding
 132     Int_t   unfFrac   = 1;      // Remove cells with less than 1% of cluster energy after unfolding
 133
 134     TString arrayName = "";
 135     TString clusterzerName = "V1"; // "V2", "NxN","V1Unfold"
 136
 137     AliAnalysisTaskEMCALClusterize * clv1 = AddTaskEMCALClusterize(arrayNameV1,outAOD,kMC,exo,clusterizerName,clTrigger, clTM,
 138                                                                    minEcell,minEseed,dTime,wTime,unfMinE,unfFrac,
 139                                                                    calibEnergy,badMap,calibTime,nonlinear);
 140
 141     //    AliAnalysisTaskEMCALClusterize * clusterize = new AliAnalysisTaskEMCALClusterize();
 142 //    clusterize->SetConfigFileName("ConfigEMCALClusterize.C");
 143 //    clusterize->SetOCDBPath("local://$ALICE_ROOT/OCDB");
 144 //    mgr->AddTask(clusterize);
 145 //
 146 //    // Create containers for input/output
 147 //    AliAnalysisDataContainer *cinput1  = mgr->GetCommonInputContainer()  ;
 148 //    AliAnalysisDataContainer *coutput1 = mgr->GetCommonOutputContainer() ;
 149 //
 150 //    mgr->ConnectInput  (clusterize, 0,  cinput1 );
 151 //    mgr->ConnectOutput (clusterize, 0, coutput1 );
 152
 153     //-----------------------
 154     // Run the analysis
 155     //-----------------------
 156     TString smode = "";
 157     if (mode==mLocal || mode == mLocalCAF)
 158       smode = "local";
 159     else if (mode==mPROOF)
 160       smode = "proof";
 161     else if (mode==mGRID)
 162       smode = "local";
 163
 164     mgr->InitAnalysis();
 165     mgr->PrintStatus();
 166     mgr->StartAnalysis(smode.Data(),chain);
 167
 168     cout <<" Analysis ended sucessfully "<< endl ;
 169
 170   }
 171   else cout << "Chain was not produced ! "<<endl;
 172
 173 }
 174
 175 void  LoadLibraries(const anaModes mode) {
 176
 177   //--------------------------------------
 178   // Load the needed libraries most of them already loaded by aliroot
 179   //--------------------------------------
 180   gSystem->Load("libTree");
 181   gSystem->Load("libGeom");
 182   gSystem->Load("libVMC");
 183   gSystem->Load("libXMLIO");
 184   gSystem->Load("libMatrix");
 185   gSystem->Load("libPhysics");
 186
 187   //----------------------------------------------------------
 188   // >>>>>>>>>>> Local mode <<<<<<<<<<<<<<
 189   //----------------------------------------------------------
 190   if (mode==mLocal || mode == mLocalCAF || mode == mGRID) {
 191     //--------------------------------------------------------
 192     // If you want to use already compiled libraries
 193     // in the aliroot distribution
 194     //--------------------------------------------------------
 195
 196     gSystem->Load("libSTEERBase");
 197     gSystem->Load("libESD");
 198     gSystem->Load("libAOD");
 199     gSystem->Load("libANALYSIS");
 200     gSystem->Load("libANALYSISalice");
 201     gSystem->Load("libANALYSISalice");
 202
 203     //SetupPar("EMCALUtils");
 204     //SetupPar("PWGEMCAL");
 205     //SetupPar("PWGGAEMCALTasks");
 206
 207     gSystem->Load("libEMCALUtils");
 208     gSystem->Load("libPWGEMCAL");
 209     gSystem->Load("libPWGGAEMCALTasks");
 210
 211   }
 212
 213   //---------------------------------------------------------
 214   // <<<<<<<<<< PROOF mode >>>>>>>>>>>>
 215   //---------------------------------------------------------
 216   else if (mode==mPROOF) {
 217     //
 218     // Connect to proof
 219     // Put appropriate username here
 220     // TProof::Reset("proof://mgheata@lxb6046.cern.ch");
 221     TProof::Open("proof://mgheata@lxb6046.cern.ch");
 222
 223     //    gProof->ClearPackages();
 224     //    gProof->ClearPackage("ESD");
 225     //    gProof->ClearPackage("AOD");
 226     //    gProof->ClearPackage("ANALYSIS");
 227
 228     // Enable the STEERBase Package
 229     gProof->UploadPackage("STEERBase.par");
 230     gProof->EnablePackage("STEERBase");
 231     // Enable the ESD Package
 232     gProof->UploadPackage("ESD.par");
 233     gProof->EnablePackage("ESD");
 234     // Enable the AOD Package
 235     gProof->UploadPackage("AOD.par");
 236     gProof->EnablePackage("AOD");
 237     // Enable the Analysis Package
 238     gProof->UploadPackage("ANALYSIS.par");
 239     gProof->EnablePackage("ANALYSIS");
 240     // Enable the PHOS geometry Package
 241     //gProof->UploadPackage("PHOSUtils.par");
 242     //gProof->EnablePackage("PHOSUtils");
 243     gProof->ShowEnabledPackages();
 244   }
 245
 246 }
 247
 248 void SetupPar(char* pararchivename)
 249 {
 250   //Load par files, create analysis libraries
 251   //For testing, if par file already decompressed and modified
 252   //classes then do not decompress.
 253
 254   TString cdir(Form("%s", gSystem->WorkingDirectory() )) ;
 255   TString parpar(Form("%s.par", pararchivename)) ;
 256   if ( gSystem->AccessPathName(parpar.Data()) ) {
 257     gSystem->ChangeDirectory(gSystem->Getenv("ALICE_ROOT")) ;
 258     TString processline(Form(".! make %s", parpar.Data())) ;
 259     gROOT->ProcessLine(processline.Data()) ;
 260     gSystem->ChangeDirectory(cdir) ;
 261     processline = Form(".! mv $ALICE_ROOT/%s .", parpar.Data()) ;
 262     gROOT->ProcessLine(processline.Data()) ;
 263   }
 264   if ( gSystem->AccessPathName(pararchivename) ) {
 265     TString processline = Form(".! tar xvzf %s",parpar.Data()) ;
 266     gROOT->ProcessLine(processline.Data());
 267   }
 268
 269   TString ocwd = gSystem->WorkingDirectory();
 270   gSystem->ChangeDirectory(pararchivename);
 271
 272   // check for BUILD.sh and execute
 273   if (!gSystem->AccessPathName("PROOF-INF/BUILD.sh")) {
 274     printf("*******************************\n");
 275     printf("*** Building PAR archive    ***\n");
 276     cout<<pararchivename<<endl;
 277     printf("*******************************\n");
 278
 279     if (gSystem->Exec("PROOF-INF/BUILD.sh")) {
 280       Error("runProcess","Cannot Build the PAR Archive! - Abort!");
 281       return -1;
 282     }
 283   }
 284   // check for SETUP.C and execute
 285   if (!gSystem->AccessPathName("PROOF-INF/SETUP.C")) {
 286     printf("*******************************\n");
 287     printf("*** Setup PAR archive       ***\n");
 288     cout<<pararchivename<<endl;
 289     printf("*******************************\n");
 290     gROOT->Macro("PROOF-INF/SETUP.C");
 291   }
 292
 293   gSystem->ChangeDirectory(ocwd.Data());
 294   printf("Current dir: %s\n", ocwd.Data());
 295 }
 296
 297
 298
 299 void CreateChain(const anaModes mode, TChain * chain){
 300   //Fills chain with data
 301   TString ocwd = gSystem->WorkingDirectory();
 302
 303   //-----------------------------------------------------------
 304   //Analysis of CAF data locally and with PROOF
 305   //-----------------------------------------------------------
 306   if(mode ==mPROOF || mode ==mLocalCAF){
 307     // Chain from CAF
 308     gROOT->LoadMacro("$ALICE_ROOT/PWG0/CreateESDChain.C");
 309     // The second parameter is the number of input files in the chain
 310     chain = CreateESDChain("ESD12001.txt", 5);
 311   }
 312
 313   //---------------------------------------
 314   //Local files analysis
 315   //---------------------------------------
 316   else if(mode == mLocal){
 317     //If you want to add several ESD files sitting in a common directory INDIR
 318     //Specify as environmental variables the directory (INDIR), the number of files
 319     //to analyze (NFILES) and the pattern name of the directories with files (PATTERN)
 320
 321     if(gSystem->Getenv("INDIR"))
 322       kInDir = gSystem->Getenv("INDIR") ;
 323     else cout<<"INDIR not set, use default: "<<kInDir<<endl;
 324
 325     if(gSystem->Getenv("PATTERN"))
 326       kPattern = gSystem->Getenv("PATTERN") ;
 327     else  cout<<"PATTERN not set, use default: "<<kPattern<<endl;
 328
 329     if(gSystem->Getenv("NFILES"))
 330       kFile = atoi(gSystem->Getenv("NFILES")) ;
 331     else cout<<"NFILES not set, use default: "<<kFile<<endl;
 332
 333     //Check if env variables are set and are correct
 334     if ( kInDir  && kFile) {
 335       printf("Get %d files from directory %s\n",kFile,kInDir);
 336       if ( ! gSystem->cd(kInDir) ) {//check if ESDs directory exist
 337         printf("%s does not exist\n", kInDir) ;
 338         return ;
 339       }
 340
 341
 342       cout<<"INDIR   : "<<kInDir<<endl;
 343       cout<<"NFILES  : "<<kFile<<endl;
 344       cout<<"PATTERN : " <<kPattern<<endl;
 345
 346       TString datafile="";
 347       if(kInputData == "ESD") datafile = "AliESDs.root" ;
 348       else if(kInputData == "AOD") datafile = "AliAODs.root" ;
 349       else if(kInputData == "MC")  datafile = "galice.root" ;
 350
 351       //Loop on ESD files, add them to chain
 352       Int_t event =0;
 353       Int_t skipped=0 ;
 354       char file[120] ;
 355
 356       for (event = 0 ; event < kFile ; event++) {
 357         sprintf(file, "%s/%s%d/%s", kInDir,kPattern,event,datafile.Data()) ;
 358         TFile * fESD = 0 ;
 359         //Check if file exists and add it, if not skip it
 360         if ( fESD = TFile::Open(file)) {
 361           if ( fESD->Get(kTreeName) ) {
 362             printf("++++ Adding %s\n", file) ;
 363             chain->AddFile(file);
 364           }
 365         }
 366         else {
 367           printf("---- Skipping %s\n", file) ;
 368           skipped++ ;
 369         }
 370       }
 371       printf("number of entries # %lld, skipped %d\n", chain->GetEntries(), skipped*100) ;
 372     }
 373     else {
 374       TString input = "AliESDs.root" ;
 375       cout<<">>>>>> No list added, take a single file <<<<<<<<< "<<input<<endl;
 376       chain->AddFile(input);
 377     }
 378
 379   }// local files analysis
 380
 381   //------------------------------
 382   //GRID xml files
 383   //-----------------------------
 384   else if(mode == mGRID){
 385     //Get colection file. It is specified by the environmental
 386     //variable XML
 387
 388     if(gSystem->Getenv("XML") )
 389       kXML = gSystem->Getenv("XML");
 390     else
 391       sprintf(kXML, "collection.xml") ;
 392
 393     if (!TFile::Open(kXML)) {
 394       printf("No collection file with name -- %s -- was found\n",kXML);
 395       return ;
 396     }
 397     else cout<<"XML file "<<kXML<<endl;
 398
 399     //Load necessary libraries and connect to the GRID
 400     gSystem->Load("libNetx") ;
 401     gSystem->Load("libRAliEn");
 402     TGrid::Connect("alien://") ;
 403
 404     //Feed Grid with collection file
 405     //TGridCollection * collection =  (TGridCollection*)gROOT->ProcessLine(Form("TAlienCollection::Open(\"%s\", 0)", kXML));
 406     TGridCollection * collection = (TGridCollection*) TAlienCollection::Open(kXML);
 407     if (! collection) {
 408       AliError(Form("%s not found", kXML)) ;
 409       return kFALSE ;
 410     }
 411     TGridResult* result = collection->GetGridResult("",0 ,0);
 412
 413     // Makes the ESD chain
 414     printf("*** Getting the Chain       ***\n");
 415     for (Int_t index = 0; index < result->GetEntries(); index++) {
 416       TString alienURL = result->GetKey(index, "turl") ;
 417       cout << "================== " << alienURL << endl ;
 418       chain->Add(alienURL) ;
 419     }
 420   }// xml analysis
 421
 422   gSystem->ChangeDirectory(ocwd.Data());
 423 }
 424