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