PWG4/CaloCalib/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 = "/user/data/files/";
  22 char * kPattern = ""; // Data are in files kInDir/kPattern+i
  23 Int_t kFile = 1; // Number of files
  24 //---------------------------------------------------------------------------
  25 //Collection file for grid analysis
  26 char * kXML = "collection.xml";
  27 //---------------------------------------------------------------------------
  28
  29 const TString kInputData = "ESD"; //ESD, AOD, MC
  30 TString kTreeName = "esdTree";
  31 Bool_t kUsePhysSel = kFALSE;
  32
  33 void emcalReclusterize(Int_t mode=mLocal)
  34 {
  35   // Main
  36   char cmd[200] ;
  37   sprintf(cmd, ".! rm -rf aod.root") ;
  38   gROOT->ProcessLine(cmd) ;
  39   //--------------------------------------------------------------------
  40   // Load analysis libraries
  41   // Look at the method below,
  42   // change whatever you need for your analysis case
  43   // ------------------------------------------------------------------
  44   LoadLibraries(mode) ;
  45   //gSystem->Unload("libPWG4CaloCalib.so");
  46   //Try to set the new library
  47   //gSystem->Load("./PWG4CaloCalib/libPWG4CaloCalib.so");
  48   //gSystem->ListLibraries();
  49
  50   //-------------------------------------------------------------------------------------------------
  51   //Create chain from ESD and from cross sections files, look below for options.
  52   //-------------------------------------------------------------------------------------------------
  53   if(kInputData == "ESD") kTreeName = "esdTree" ;
  54   else if(kInputData == "AOD") kTreeName = "aodTree" ;
  55   else {
  56     cout<<"Wrong  data type "<<kInputData<<endl;
  57     break;
  58   }
  59
  60   TChain *chain       = new TChain(kTreeName) ;
  61   CreateChain(mode, chain);
  62
  63   if(chain){
  64     AliLog::SetGlobalLogLevel(AliLog::kError);//Minimum prints on screen
  65
  66     //--------------------------------------
  67     // Make the analysis manager
  68     //-------------------------------------
  69     AliAnalysisManager *mgr  = new AliAnalysisManager("Manager", "Manager");
  70
  71     // AOD output handler
  72     AliAODHandler* aodoutHandler   = new AliAODHandler();
  73     aodoutHandler->SetOutputFileName("aod.root");
  74     ////aodoutHandler->SetCreateNonStandardAOD();
  75     mgr->SetOutputEventHandler(aodoutHandler);
  76
  77     //input
  78     if(kInputData == "ESD")
  79       {
  80         // ESD handler
  81         AliESDInputHandler *esdHandler = new AliESDInputHandler();
  82         mgr->SetInputEventHandler(esdHandler);
  83         esdHandler->SetReadFriends(kFALSE);
  84         cout<<"ESD handler "<<mgr->GetInputEventHandler()<<endl;
  85       }
  86     if(kInputData == "AOD")
  87       {
  88         // AOD handler
  89         AliAODInputHandler *aodHandler = new AliAODInputHandler();
  90         mgr->SetInputEventHandler(aodHandler);
  91         cout<<"AOD handler "<<mgr->GetInputEventHandler()<<endl;
  92
  93       }
  94
  95     mgr->SetDebugLevel(1);
  96
  97     //-------------------------------------------------------------------------
  98     //Define task, put here any other task that you want to use.
  99     //-------------------------------------------------------------------------
 100     AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer();
 101     AliAnalysisDataContainer *coutput1 = mgr->GetCommonOutputContainer();
 102
 103     // ESD physics selection task
 104     if(kInputData == "ESD" && kUsePhysSel){
 105       gROOT->LoadMacro("AddTaskPhysicsSelection.C");
 106       AliPhysicsSelectionTask* physSelTask = AddTaskPhysicsSelection();
 107     }
 108
 109     // Create containers for input/output
 110     AliAnalysisDataContainer *cinput1  = mgr->GetCommonInputContainer();
 111     AliAnalysisDataContainer *coutput1 = mgr->GetCommonOutputContainer();
 112
 113     AliAnalysisTaskEMCALClusterize * clusterize = new AliAnalysisTaskEMCALClusterize("EMCALClusterize");
 114     if(kUsePhysSel)clusterize->SelectCollisionCandidates();
 115     //clusterize->SetOCDBPath("local://$ALICE_ROOT/OCDB");
 116     clusterize->FillAODFile(kFALSE); // fill aod.root with clusters?, not really needed for analysis.
 117     clusterize->JustUnfold(kTRUE); // if TRUE, do just unfolding, do not recluster cells
 118     AliEMCALRecParam * params = clusterize->GetRecParam();
 119     params->SetClusterizerFlag(AliEMCALRecParam::kClusterizerNxN);
 120     params->SetClusteringThreshold(0.1); // 100 MeV
 121     params->SetMinECut(0.01);  //10 MeV
 122     params->SetUnfold(kFALSE);
 123     params->SetW0(4.5);
 124     params->SetTimeCut(-1);///Open this cut for AODs
 125     params->SetTimeMin(-1);//Open this cut for AODs
 126     params->SetTimeMax(1e6);//Open this cut for AODs
 127
 128 //    TGeoHMatrix *matrix[4];
 129
 130 //     double rotationMatrix[4][9] = {-0.014585, -0.999892, -0.002031, 0.999892, -0.014589,  0.001950, -0.001979, -0.002003,  0.999996,
 131 //                                 -0.014585,  0.999892,  0.002031, 0.999892,  0.014589, -0.001950, -0.001979,  0.002003, -0.999996,
 132 //                                 -0.345861, -0.938280, -0.003412, 0.938281, -0.345869,  0.001950, -0.003010, -0.002527,  0.999992,
 133 //                                 -0.345861,  0.938280,  0.003412, 0.938281,  0.345869, -0.001950, -0.003010,  0.002527, -0.999992};
 134
 135 //     double translationMatrix[4][3] = {0.367264,    446.508738,  175.97185+0.3,
 136 //                                    1.078181,    445.826258, -174.026758+0.3,
 137 //                                    -153.843916, 418.304256,  175.956905+0.8,
 138 //                                    -152.649580, 417.621779, -174.040392+0.8};
 139
 140
 141 //     double rotationMatrix[4][9] = {-0.014587, -0.999892, -0.002031, 0.999892, -0.014591,  0.001979, -0.002009, -0.002002,  0.999996,
 142 //                               -0.014587,  0.999892,  0.002031, 0.999892,  0.014591, -0.001979, -0.002009,  0.002002, -0.999996,
 143 //                               -0.345864, -0.938278, -0.003412, 0.938276, -0.345874,  0.003010, -0.004004, -0.002161,  0.999990,
 144 //                               -0.345864,  0.938278,  0.003412, 0.938276,  0.345874, -0.003010, -0.004004,  0.002161, -0.999990};
 145
 146 //     double translationMatrix[4][3] = {0.351659,    447.576446,  176.269742,
 147 //                                    1.062577,    446.893974, -173.728870,
 148 //                                    -154.213287, 419.306156,  176.753692,
 149 //                                    -153.018950, 418.623681, -173.243605};
 150
 151 //     for(int j=0; j<4; j++)
 152 //       {
 153 //      matrix[j] = new TGeoHMatrix();
 154 //      matrix[j]->SetRotation(rotationMatrix[j]);
 155 //      matrix[j]->SetTranslation(translationMatrix[j]);
 156 //      matrix[j]->Print();
 157 //      clusterize->SetGeometryMatrixInSM(matrix[j],j);
 158 //       }
 159
 160
 161 //     clusterize->SwitchOnLoadOwnEMCALGeometryMatrices();
 162
 163
 164     mgr->AddTask(clusterize);
 165     mgr->ConnectInput  (clusterize,  0, cinput1);
 166     mgr->ConnectOutput (clusterize, 0, coutput1 );
 167
 168
 169     //-----------------------
 170     // Run the analysis
 171     //-----------------------
 172     TString smode = "";
 173     if (mode==mLocal || mode == mLocalCAF)
 174       smode = "local";
 175     else if (mode==mPROOF)
 176       smode = "proof";
 177     else if (mode==mGRID)
 178       smode = "local";
 179
 180     mgr->InitAnalysis();
 181     mgr->PrintStatus();
 182     mgr->StartAnalysis(smode.Data(),chain);
 183
 184 cout <<" Analysis ended sucessfully "<< endl ;
 185
 186   }
 187   else cout << "Chain was not produced ! "<<endl;
 188
 189 }
 190
 191 void  LoadLibraries(const anaModes mode) {
 192
 193   //--------------------------------------
 194   // Load the needed libraries most of them already loaded by aliroot
 195   //--------------------------------------
 196   gSystem->Load("libTree.so");
 197   gSystem->Load("libGeom.so");
 198   gSystem->Load("libVMC.so");
 199   gSystem->Load("libXMLIO.so");
 200   gSystem->Load("libMatrix.so");
 201   gSystem->Load("libPhysics.so");
 202
 203   //----------------------------------------------------------
 204   // >>>>>>>>>>> Local mode <<<<<<<<<<<<<<
 205   //----------------------------------------------------------
 206   if (mode==mLocal || mode == mLocalCAF || mode == mGRID) {
 207     //--------------------------------------------------------
 208     // If you want to use already compiled libraries
 209     // in the aliroot distribution
 210     //--------------------------------------------------------
 211
 212     gSystem->Load("libSTEERBase.so");
 213     gSystem->Load("libESD.so");
 214     gSystem->Load("libAOD.so");
 215     gSystem->Load("libANALYSIS.so");
 216     gSystem->Load("libANALYSISalice.so");
 217     gSystem->Load("libANALYSISalice.so");
 218     gSystem->Load("libEMCALUtils.so");
 219
 220     //SetupPar("EMCALUtils");
 221     //SetupPar("PWG4CaloCalib");
 222
 223     TGeoManager::Import("geometry.root") ; //need file "geometry.root" in local dir!!!!
 224     gSystem->Load("libPWG4CaloCalib.so");
 225
 226     /*
 227       //     gSystem->Load("libPWG4omega3pi.so");
 228       //     gSystem->Load("libCORRFW.so");
 229       //     gSystem->Load("libPWG3base.so");
 230       //     gSystem->Load("libPWG3muon.so");
 231  */
 232     //--------------------------------------------------------
 233     //If you want to use root and par files from aliroot
 234     //--------------------------------------------------------
 235     /*
 236            SetupPar("STEERBase");
 237            SetupPar("ESD");
 238            SetupPar("AOD");
 239            SetupPar("ANALYSIS");
 240            SetupPar("ANALYSISalice");
 241            SetupPar("PHOSUtils");
 242            SetupPar("EMCALUtils");
 243            //Create Geometry
 244            TGeoManager::Import("geometry.root") ; //need file "geometry.root" in local dir!!!!
 245            SetupPar("PWG4CaloCalib");
 246 */
 247   }
 248
 249   //---------------------------------------------------------
 250   // <<<<<<<<<< PROOF mode >>>>>>>>>>>>
 251   //---------------------------------------------------------
 252   else if (mode==mPROOF) {
 253     //
 254     // Connect to proof
 255     // Put appropriate username here
 256     // TProof::Reset("proof://mgheata@lxb6046.cern.ch");
 257     TProof::Open("proof://mgheata@lxb6046.cern.ch");
 258
 259     //    gProof->ClearPackages();
 260     //    gProof->ClearPackage("ESD");
 261     //    gProof->ClearPackage("AOD");
 262     //    gProof->ClearPackage("ANALYSIS");
 263     //    gProof->ClearPackage("PWG4PartCorrBase");
 264     //    gProof->ClearPackage("PWG4PartCorrDep");
 265
 266     // Enable the STEERBase Package
 267     gProof->UploadPackage("STEERBase.par");
 268     gProof->EnablePackage("STEERBase");
 269     // Enable the ESD Package
 270     gProof->UploadPackage("ESD.par");
 271     gProof->EnablePackage("ESD");
 272     // Enable the AOD Package
 273     gProof->UploadPackage("AOD.par");
 274     gProof->EnablePackage("AOD");
 275     // Enable the Analysis Package
 276     gProof->UploadPackage("ANALYSIS.par");
 277     gProof->EnablePackage("ANALYSIS");
 278     // Enable the PHOS geometry Package
 279     //gProof->UploadPackage("PHOSUtils.par");
 280     //gProof->EnablePackage("PHOSUtils");
 281     // Enable PartCorr analysis
 282     gProof->UploadPackage("PWG4PartCorrBase.par");
 283     gProof->EnablePackage("PWG4PartCorrBase");
 284     gProof->UploadPackage("PWG4PartCorrDep.par");
 285     gProof->EnablePackage("PWG4PartCorrDep");
 286     gProof->ShowEnabledPackages();
 287   }
 288
 289 }
 290
 291 void SetupPar(char* pararchivename)
 292 {
 293   //Load par files, create analysis libraries
 294   //For testing, if par file already decompressed and modified
 295   //classes then do not decompress.
 296
 297   TString cdir(Form("%s", gSystem->WorkingDirectory() )) ;
 298   TString parpar(Form("%s.par", pararchivename)) ;
 299   if ( gSystem->AccessPathName(parpar.Data()) ) {
 300     gSystem->ChangeDirectory(gSystem->Getenv("ALICE_ROOT")) ;
 301     TString processline(Form(".! make %s", parpar.Data())) ;
 302     gROOT->ProcessLine(processline.Data()) ;
 303     gSystem->ChangeDirectory(cdir) ;
 304     processline = Form(".! mv $ALICE_ROOT/%s .", parpar.Data()) ;
 305     gROOT->ProcessLine(processline.Data()) ;
 306   }
 307   if ( gSystem->AccessPathName(pararchivename) ) {
 308     TString processline = Form(".! tar xvzf %s",parpar.Data()) ;
 309     gROOT->ProcessLine(processline.Data());
 310   }
 311
 312   TString ocwd = gSystem->WorkingDirectory();
 313   gSystem->ChangeDirectory(pararchivename);
 314
 315   // check for BUILD.sh and execute
 316   if (!gSystem->AccessPathName("PROOF-INF/BUILD.sh")) {
 317     printf("*******************************\n");
 318     printf("*** Building PAR archive    ***\n");
 319     cout<<pararchivename<<endl;
 320     printf("*******************************\n");
 321
 322     if (gSystem->Exec("PROOF-INF/BUILD.sh")) {
 323       Error("runProcess","Cannot Build the PAR Archive! - Abort!");
 324       return -1;
 325     }
 326   }
 327   // check for SETUP.C and execute
 328   if (!gSystem->AccessPathName("PROOF-INF/SETUP.C")) {
 329     printf("*******************************\n");
 330     printf("*** Setup PAR archive       ***\n");
 331     cout<<pararchivename<<endl;
 332     printf("*******************************\n");
 333     gROOT->Macro("PROOF-INF/SETUP.C");
 334   }
 335
 336   gSystem->ChangeDirectory(ocwd.Data());
 337   printf("Current dir: %s\n", ocwd.Data());
 338 }
 339
 340
 341
 342 void CreateChain(const anaModes mode, TChain * chain){
 343   //Fills chain with data
 344   TString ocwd = gSystem->WorkingDirectory();
 345
 346   //-----------------------------------------------------------
 347   //Analysis of CAF data locally and with PROOF
 348   //-----------------------------------------------------------
 349   if(mode ==mPROOF || mode ==mLocalCAF){
 350     // Chain from CAF
 351     gROOT->LoadMacro("$ALICE_ROOT/PWG0/CreateESDChain.C");
 352     // The second parameter is the number of input files in the chain
 353     chain = CreateESDChain("ESD12001.txt", 5);
 354   }
 355
 356   //---------------------------------------
 357   //Local files analysis
 358   //---------------------------------------
 359   else if(mode == mLocal){
 360     //If you want to add several ESD files sitting in a common directory INDIR
 361     //Specify as environmental variables the directory (INDIR), the number of files
 362     //to analyze (NFILES) and the pattern name of the directories with files (PATTERN)
 363
 364     if(gSystem->Getenv("INDIR"))
 365       kInDir = gSystem->Getenv("INDIR") ;
 366     else cout<<"INDIR not set, use default: "<<kInDir<<endl;
 367
 368     if(gSystem->Getenv("PATTERN"))
 369       kPattern = gSystem->Getenv("PATTERN") ;
 370     else  cout<<"PATTERN not set, use default: "<<kPattern<<endl;
 371
 372     if(gSystem->Getenv("NFILES"))
 373       kFile = atoi(gSystem->Getenv("NFILES")) ;
 374     else cout<<"NFILES not set, use default: "<<kFile<<endl;
 375
 376     //Check if env variables are set and are correct
 377     if ( kInDir  && kFile) {
 378       printf("Get %d files from directory %s\n",kFile,kInDir);
 379       if ( ! gSystem->cd(kInDir) ) {//check if ESDs directory exist
 380         printf("%s does not exist\n", kInDir) ;
 381         return ;
 382       }
 383
 384
 385       cout<<"INDIR   : "<<kInDir<<endl;
 386       cout<<"NFILES  : "<<kFile<<endl;
 387       cout<<"PATTERN : " <<kPattern<<endl;
 388
 389       TString datafile="";
 390       if(kInputData == "ESD") datafile = "AliESDs.root" ;
 391       else if(kInputData == "AOD") datafile = "AliAODs.root" ;
 392       else if(kInputData == "MC")  datafile = "galice.root" ;
 393
 394       //Loop on ESD files, add them to chain
 395       Int_t event =0;
 396       Int_t skipped=0 ;
 397       char file[120] ;
 398
 399       for (event = 0 ; event < kFile ; event++) {
 400         sprintf(file, "%s/%s%d/%s", kInDir,kPattern,event,datafile.Data()) ;
 401         TFile * fESD = 0 ;
 402         //Check if file exists and add it, if not skip it
 403         if ( fESD = TFile::Open(file)) {
 404           if ( fESD->Get(kTreeName) ) {
 405             printf("++++ Adding %s\n", file) ;
 406             chain->AddFile(file);
 407           }
 408         }
 409         else {
 410           printf("---- Skipping %s\n", file) ;
 411           skipped++ ;
 412         }
 413       }
 414       printf("number of entries # %lld, skipped %d\n", chain->GetEntries(), skipped*100) ;
 415     }
 416     else {
 417       TString input = "AliESDs.root" ;
 418       cout<<">>>>>> No list added, take a single file <<<<<<<<< "<<input<<endl;
 419       chain->AddFile(input);
 420     }
 421
 422   }// local files analysis
 423
 424   //------------------------------
 425   //GRID xml files
 426   //-----------------------------
 427   else if(mode == mGRID){
 428     //Get colection file. It is specified by the environmental
 429     //variable XML
 430
 431     if(gSystem->Getenv("XML") )
 432       kXML = gSystem->Getenv("XML");
 433     else
 434       sprintf(kXML, "collection.xml") ;
 435
 436     if (!TFile::Open(kXML)) {
 437       printf("No collection file with name -- %s -- was found\n",kXML);
 438       return ;
 439     }
 440     else cout<<"XML file "<<kXML<<endl;
 441
 442     //Load necessary libraries and connect to the GRID
 443     gSystem->Load("libNetx.so") ;
 444     gSystem->Load("libRAliEn.so");
 445     TGrid::Connect("alien://") ;
 446
 447     //Feed Grid with collection file
 448     //TGridCollection * collection =  (TGridCollection*)gROOT->ProcessLine(Form("TAlienCollection::Open(\"%s\", 0)", kXML));
 449     TGridCollection * collection = (TGridCollection*) TAlienCollection::Open(kXML);
 450     if (! collection) {
 451       AliError(Form("%s not found", kXML)) ;
 452       return kFALSE ;
 453     }
 454     TGridResult* result = collection->GetGridResult("",0 ,0);
 455
 456     // Makes the ESD chain
 457     printf("*** Getting the Chain       ***\n");
 458     for (Int_t index = 0; index < result->GetEntries(); index++) {
 459       TString alienURL = result->GetKey(index, "turl") ;
 460       cout << "================== " << alienURL << endl ;
 461       chain->Add(alienURL) ;
 462     }
 463   }// xml analysis
 464
 465   gSystem->ChangeDirectory(ocwd.Data());
 466 }
 467