ANALYSIS/macros/AODtrain.C

   1 // ### Settings that make sense when using the Alien plugin
   2 //==============================================================================
   3 Int_t       runOnData          = 1;       // Set to 1 if processing real data
   4 Int_t       iCollision         = 1;       // 0=pp, 1=Pb-Pb
   5 Int_t       run_flag           = 1100;    // year (2011 = 1100)
   6 //==============================================================================
   7 Bool_t      usePhysicsSelection = kTRUE; // use physics selection
   8 Bool_t      useTender           = kFALSE; // use tender wagon
   9 Bool_t      useCentrality       = kTRUE; // centrality
  10 Bool_t      useV0tender         = kFALSE;  // use V0 correction in tender
  11 Bool_t      useDBG              = kTRUE;  // activate debugging
  12 Bool_t      useMC               = kFALSE;  // use MC info
  13 Bool_t      useKFILTER          = kFALSE;  // use Kinematics filter
  14 Bool_t      useTR               = kFALSE;  // use track references
  15 Bool_t      useCORRFW           = kFALSE; // do not change
  16 Bool_t      useAODTAGS          = kFALSE; // use AOD tags
  17 Bool_t      useSysInfo          = kFALSE; // use sys info
  18
  19 // ### Analysis modules to be included. Some may not be yet fully implemented.
  20 //==============================================================================
  21 Int_t       iAODhandler        = 1;      // Analysis produces an AOD or dAOD's
  22 Int_t       iESDfilter         = 1;      // ESD to AOD filter (barrel + muon tracks)
  23 Int_t       iMUONcopyAOD       = 1;      // Task that copies only muon events in a separate AOD (PWG3)
  24 Int_t       iJETAN             = 1;      // Jet analysis (PWG4)
  25 Int_t       iJETANdelta        = 1;      // Jet delta AODs
  26 Int_t       iPWGHFvertexing     = 1;      // Vertexing HF task (PWG3)
  27 Int_t       iPWGDQJPSIfilter    = 1;      // JPSI filtering (PWG3)
  28 Int_t       iPWGHFd2h           = 1;      // D0->2 hadrons (PWG3)
  29
  30 // ### Configuration macros used for each module
  31 //==============================================================================
  32  TString configPWGHFd2h = (iCollision==0)?"$ALICE_ROOT/PWGHF/vertexingHF/ConfigVertexingHF.C"
  33                           :"$ALICE_ROOT/PWGHF/vertexingHF/ConfigVertexingHF_highmult.C";
  34
  35 // Temporaries.
  36 void AODmerge();
  37 void AddAnalysisTasks();
  38 Bool_t LoadCommonLibraries();
  39 Bool_t LoadAnalysisLibraries();
  40 Bool_t LoadLibrary(const char *);
  41 TChain *CreateChain();
  42
  43 //______________________________________________________________________________
  44 void AODtrain(Int_t merge=0)
  45 {
  46 // Main analysis train macro.
  47
  48   if (merge) {
  49     TGrid::Connect("alien://");
  50     if (!gGrid || !gGrid->IsConnected()) {
  51       ::Error("QAtrain", "No grid connection");
  52       return;
  53     }
  54   }
  55   // Set temporary merging directory to current one
  56   gSystem->Setenv("TMPDIR", gSystem->pwd());
  57   // Set temporary compilation directory to current one
  58   gSystem->SetBuildDir(gSystem->pwd(), kTRUE);
  59    printf("==================================================================\n");
  60    printf("===========    RUNNING FILTERING TRAIN   ==========\n");
  61    printf("==================================================================\n");
  62    printf("=  Configuring analysis train for:                               =\n");
  63    if (usePhysicsSelection)   printf("=  Physics selection                                                =\n");
  64    if (useTender)    printf("=  TENDER                                                        =\n");
  65    if (iESDfilter)   printf("=  ESD filter                                                    =\n");
  66    if (iMUONcopyAOD) printf("=  MUON copy AOD                                                 =\n");
  67    if (iJETAN)       printf("=  Jet analysis                                                  =\n");
  68    if (iJETANdelta)  printf("=     Jet delta AODs                                             =\n");
  69    if (iPWGHFvertexing) printf("=  PWGHF vertexing                                                =\n");
  70    if (iPWGDQJPSIfilter) printf("=  PWGDQ j/psi filter                                             =\n");
  71    if (iPWGHFd2h) printf("=  PWGHF D0->2 hadrons QA                                     =\n");
  72
  73    // Load common libraries and set include path
  74    if (!LoadCommonLibraries()) {
  75       ::Error("AnalysisTrain", "Could not load common libraries");
  76       return;
  77    }
  78
  79    // Make the analysis manager and connect event handlers
  80    AliAnalysisManager *mgr  = new AliAnalysisManager("Analysis Train", "Production train");
  81    if (useSysInfo) mgr->SetNSysInfo(100);
  82    // Load analysis specific libraries
  83    if (!LoadAnalysisLibraries()) {
  84       ::Error("AnalysisTrain", "Could not load analysis libraries");
  85       return;
  86    }
  87
  88    // Create input handler (input container created automatically)
  89    // ESD input handler
  90    AliESDInputHandler *esdHandler = new AliESDInputHandler();
  91    mgr->SetInputEventHandler(esdHandler);
  92    // Monte Carlo handler
  93    if (useMC) {
  94       AliMCEventHandler* mcHandler = new AliMCEventHandler();
  95       mgr->SetMCtruthEventHandler(mcHandler);
  96       mcHandler->SetReadTR(useTR);
  97    }
  98    // AOD output container, created automatically when setting an AOD handler
  99    if (iAODhandler) {
 100       // AOD output handler
 101       AliAODHandler* aodHandler   = new AliAODHandler();
 102       aodHandler->SetOutputFileName("AliAOD.root");
 103       mgr->SetOutputEventHandler(aodHandler);
 104    }
 105    // Debugging if needed
 106    if (useDBG) mgr->SetDebugLevel(3);
 107
 108    AddAnalysisTasks();
 109    if (merge) {
 110       AODmerge();
 111       mgr->InitAnalysis();
 112       mgr->SetGridHandler(new AliAnalysisAlien);
 113       mgr->StartAnalysis("gridterminate",0);
 114       return;
 115    }
 116    // Run the analysis
 117    //
 118    TChain *chain = CreateChain();
 119    if (!chain) return;
 120
 121    TStopwatch timer;
 122    timer.Start();
 123    mgr->SetSkipTerminate(kTRUE);
 124    if (mgr->InitAnalysis()) {
 125       mgr->PrintStatus();
 126       mgr->StartAnalysis("local", chain);
 127    }
 128    timer.Print();
 129 }
 130
 131 //______________________________________________________________________________
 132 void AddAnalysisTasks(){
 133   // Add all analysis task wagons to the train
 134    AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
 135
 136   //
 137   // Tender and supplies. Needs to be called for every event.
 138   //
 139    AliAnalysisManager::SetCommonFileName("AODQA.root");
 140    if (useTender) {
 141       gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/TenderSupplies/AddTaskTender.C");
 142       // IF V0 tender needed, put kTRUE below
 143       AliAnalysisTaskSE *tender = AddTaskTender(useV0tender);
 144 //      tender->SetDebugLevel(2);
 145    }
 146
 147    if (usePhysicsSelection) {
 148    // Physics selection task
 149       gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C");
 150       mgr->RegisterExtraFile("event_stat.root");
 151       AliPhysicsSelectionTask *physSelTask = AddTaskPhysicsSelection(useMC);
 152 //      AliOADBPhysicsSelection * oadbDefaultPbPb = CreateOADBphysicsSelection();
 153 //      physSelTask->GetPhysicsSelection()->SetCustomOADBObjects(oadbDefaultPbPb,0,0);
 154       mgr->AddStatisticsTask(AliVEvent::kAny);
 155    }
 156    // Centrality (only Pb-Pb)
 157    if (iCollision && useCentrality) {
 158       gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskCentrality.C");
 159       AliCentralitySelectionTask *taskCentrality = AddTaskCentrality();
 160       taskCentrality->SelectCollisionCandidates(AliVEvent::kAny);
 161    }
 162
 163    if (iESDfilter) {
 164       //  ESD filter task configuration.
 165       gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskESDFilter.C");
 166       if (iMUONcopyAOD) {
 167          printf("Registering delta AOD file\n");
 168          mgr->RegisterExtraFile("AliAOD.Muons.root");
 169          mgr->RegisterExtraFile("AliAOD.Dimuons.root");
 170          AliAnalysisTaskESDfilter *taskesdfilter = AddTaskESDFilter(useKFILTER, kTRUE, kFALSE, kFALSE /*usePhysicsSelection*/,kTRUE,kFALSE,kFALSE,run_flag);
 171       } else {
 172            AliAnalysisTaskESDfilter *taskesdfilter = AddTaskESDFilter(useKFILTER, kFALSE, kFALSE, kFALSE /*usePhysicsSelection*/,kFALSE,kTRUE,kFALSE,kFALSE,run_flag); // others
 173       }
 174    }
 175
 176 // ********** PWG3 wagons ******************************************************
 177    // PWGHF vertexing
 178    if (iPWGHFvertexing) {
 179       gROOT->LoadMacro("$ALICE_ROOT/PWGHF/vertexingHF/macros/AddTaskVertexingHF.C");
 180       if (!iPWGHFd2h) TFile::Cp(gSystem->ExpandPathName(configPWGHFd2h.Data()), "file:ConfigVertexingHF.C");
 181       AliAnalysisTaskSEVertexingHF *taskvertexingHF = AddTaskVertexingHF();
 182       if (!taskvertexingHF) ::Warning("AnalysisTrainNew", "AliAnalysisTaskSEVertexingHF cannot run for this train conditions - EXCLUDED");
 183       else mgr->RegisterExtraFile("AliAOD.VertexingHF.root");
 184       taskvertexingHF->SelectCollisionCandidates(0);
 185    }
 186
 187    // PWGDQ JPSI filtering (only pp)
 188    if (iPWGDQJPSIfilter && (iCollision==0)) {
 189       gROOT->LoadMacro("$ALICE_ROOT/PWGDQ/dielectron/macros/AddTaskJPSIFilter.C");
 190       AliAnalysisTaskSE *taskJPSIfilter = AddTaskJPSIFilter();
 191       if (!taskJPSIfilter) ::Warning("AnalysisTrainNew", "AliAnalysisTaskDielectronFilter cannot run for this train conditions - EXCLUDED");
 192       else mgr->RegisterExtraFile("AliAOD.Dielectron.root");
 193       taskJPSIfilter->SelectCollisionCandidates(0);
 194    }
 195
 196    // PWGHF D2h
 197    if (iPWGHFd2h) {
 198      gROOT->LoadMacro("$ALICE_ROOT/PWGHF/vertexingHF/AddD2HTrain.C");
 199      TFile::Cp(gSystem->ExpandPathName(configPWGHFd2h.Data()), "file:ConfigVertexingHF.C");
 200      AddD2HTrain(kFALSE, 1,0,0,0,0,0,0,0,0,0,0);
 201    }
 202
 203    // ********** PWG4 wagons ******************************************************
 204    // Jet analysis
 205
 206    // Configurations flags, move up?
 207    TString kDeltaAODJetName = "AliAOD.Jets.root"; //
 208    Bool_t  kIsPbPb = (iCollision==0)?false:true; // can be more intlligent checking the name of the data set
 209    TString kDefaultJetBackgroundBranch = "";
 210    TString kJetSubtractBranches = "";
 211    UInt_t kHighPtFilterMask = 128;// from esd filter
 212    UInt_t iPhysicsSelectionFlag = AliVEvent::kMB;
 213    if (iJETAN) {
 214      gROOT->LoadMacro("$ALICE_ROOT/PWGJE/macros/AddTaskJets.C");
 215      // Default jet reconstructor running on ESD's
 216      AliAnalysisTaskJets *taskjets = AddTaskJets("AOD","UA1",0.4,kHighPtFilterMask,1.,0); // no background subtraction
 217      if (!taskjets) ::Fatal("AnalysisTrainNew", "AliAnalysisTaskJets cannot run for this train conditions - EXCLUDED");
 218      if(kDeltaAODJetName.Length()>0) taskjets->SetNonStdOutputFile(kDeltaAODJetName.Data());
 219      if (iJETANdelta) {
 220         //            AddTaskJetsDelta("AliAOD.Jets.root"); // need to modify this accordingly in the add task jets
 221         mgr->RegisterExtraFile(kDeltaAODJetName.Data());
 222         TString cTmp("");
 223         if(kIsPbPb){
 224           // UA1 intrinsic background subtraction
 225           taskjets = AddTaskJets("AOD","UA1",0.4,kHighPtFilterMask,1.,2); // background subtraction
 226           if(kDeltaAODJetName.Length()>0)taskjets->SetNonStdOutputFile(kDeltaAODJetName.Data());
 227        }
 228        // SICONE
 229        taskjets = AddTaskJets("AOD","SISCONE",0.4,kHighPtFilterMask,0.15,0); //no background subtration to be done later....
 230        if(kDeltaAODJetName.Length()>0)taskjets->SetNonStdOutputFile(kDeltaAODJetName.Data());
 231        cTmp = taskjets->GetNonStdBranch();
 232        if(cTmp.Length()>0)kJetSubtractBranches += Form("%s ",cTmp.Data());
 233
 234        // Add the clusters..
 235        gROOT->LoadMacro("$ALICE_ROOT/PWGJE/macros/AddTaskJetCluster.C");
 236        AliAnalysisTaskJetCluster *taskCl = 0;
 237        Float_t fCenUp = 0;
 238        Float_t fCenLo = 0;
 239        Float_t fTrackEtaWindow = 0.9;
 240        taskCl = AddTaskJetCluster("AOD","",kHighPtFilterMask,iPhysicsSelectionFlag,"KT",0.4,0,1, kDeltaAODJetName.Data(),0.15,fTrackEtaWindow,0); // this one is for the background and random jets, random cones with no skip
 241        taskCl->SetBackgroundCalc(kTRUE);
 242        taskCl->SetNRandomCones(10);
 243        taskCl->SetCentralityCut(fCenLo,fCenUp);
 244        taskCl->SetGhostEtamax(fTrackEtaWindow);
 245        kDefaultJetBackgroundBranch = Form("%s_%s",AliAODJetEventBackground::StdBranchName(),taskCl->GetJetOutputBranch());
 246
 247        taskCl = AddTaskJetCluster("AOD","",kHighPtFilterMask,iPhysicsSelectionFlag,"ANTIKT",0.4,2,1,kDeltaAODJetName.Data(),0.15);
 248        taskCl->SetCentralityCut(fCenLo,fCenUp);
 249        if(kIsPbPb)taskCl->SetBackgroundBranch(kDefaultJetBackgroundBranch.Data());
 250        taskCl->SetNRandomCones(10);
 251        kJetSubtractBranches += Form("%s ",taskCl->GetJetOutputBranch());
 252
 253        taskCl = AddTaskJetCluster("AOD","",kHighPtFilterMask,iPhysicsSelectionFlag,"ANTIKT",0.2,0,1,kDeltaAODJetName.Data(),0.15);
 254        taskCl->SetCentralityCut(fCenLo,fCenUp);
 255        if(kIsPbPb)taskCl->SetBackgroundBranch(kDefaultJetBackgroundBranch.Data());
 256        kJetSubtractBranches += Form("%s ",taskCl->GetJetOutputBranch());
 257
 258        // DO THE BACKGROUND SUBTRACTION
 259        if(kIsPbPb&&kJetSubtractBranches.Length()){
 260          gROOT->LoadMacro("$ALICE_ROOT/PWGJE/macros/AddTaskJetBackgroundSubtract.C");
 261          AliAnalysisTaskJetBackgroundSubtract *taskSubtract = 0;
 262          taskSubtract = AddTaskJetBackgroundSubtract(kJetSubtractBranches,1,"B0","B%d");
 263          taskSubtract->SetBackgroundBranch(kDefaultJetBackgroundBranch.Data());
 264          if(kDeltaAODJetName.Length()>0)taskSubtract->SetNonStdOutputFile(kDeltaAODJetName.Data());
 265        }
 266      }
 267    }
 268 }
 269
 270 //______________________________________________________________________________
 271 Bool_t LoadCommonLibraries()
 272 {
 273 // Load common analysis libraries.
 274    if (!gSystem->Getenv("ALICE_ROOT")) {
 275       ::Error("AnalysisTrainNew.C::LoadCommonLibraries", "Analysis train requires that analysis libraries are compiled with a local AliRoot");
 276       return kFALSE;
 277    }
 278    Bool_t success = kTRUE;
 279    // Load framework classes. Par option ignored here.
 280    success &= LoadLibrary("libSTEERBase.so");
 281    success &= LoadLibrary("libESD.so");
 282    success &= LoadLibrary("libAOD.so");
 283    success &= LoadLibrary("libANALYSIS.so");
 284    success &= LoadLibrary("libOADB.so");
 285    success &= LoadLibrary("libANALYSISalice.so");
 286    success &= LoadLibrary("libCORRFW.so");
 287    gROOT->ProcessLine(".include $ALICE_ROOT/include");
 288    if (success) {
 289       ::Info("AnalysisTrainNew.C::LoadCommodLibraries", "Load common libraries:    SUCCESS");
 290       ::Info("AnalysisTrainNew.C::LoadCommodLibraries", "Include path for Aclic compilation:\n%s",
 291               gSystem->GetIncludePath());
 292    } else {
 293       ::Info("AnalysisTrainNew.C::LoadCommodLibraries", "Load common libraries:    FAILED");
 294    }
 295    return success;
 296 }
 297
 298 //______________________________________________________________________________
 299 Bool_t LoadAnalysisLibraries()
 300 {
 301 // Load common analysis libraries.
 302    if (useTender) {
 303       if (!LoadLibrary("TENDER") ||
 304           !LoadLibrary("TENDERSupplies")) return kFALSE;
 305    }
 306    if (iESDfilter || iPWGMuonTrain) {
 307       if (!LoadLibrary("PWGmuon")) return kFALSE;
 308    }
 309    // JETAN
 310    if (iJETAN) {
 311       if (!LoadLibrary("JETAN")) return kFALSE;
 312    }
 313    if (iJETANdelta) {
 314       if (!LoadLibrary("JETAN") ||
 315           !LoadLibrary("CGAL") ||
 316           !LoadLibrary("fastjet") ||
 317           !LoadLibrary("siscone") ||
 318           !LoadLibrary("SISConePlugin") ||
 319           !LoadLibrary("FASTJETAN")) return kFALSE;
 320    }
 321    // PWG3 Vertexing HF
 322    if (iPWGHFvertexing || iPWG3d2h) {
 323       if (!LoadLibrary("PWG3base") ||
 324           !LoadLibrary("PWGHFvertexingHF")) return kFALSE;
 325    }
 326    // PWG3 dielectron
 327    if (iPWGDQJPSIfilter) {
 328       if (!LoadLibrary("PWGDQdielectron")) return kFALSE;
 329    }
 330
 331    ::Info("AnalysisTrainNew.C::LoadAnalysisLibraries", "Load other libraries:   SUCCESS");
 332    return kTRUE;
 333 }
 334
 335 //______________________________________________________________________________
 336 Bool_t LoadLibrary(const char *module)
 337 {
 338 // Load a module library in a given mode. Reports success.
 339    Int_t result;
 340    TString mod(module);
 341    if (!mod.Length()) {
 342       ::Error("AnalysisTrainNew.C::LoadLibrary", "Empty module name");
 343       return kFALSE;
 344    }
 345    // If a library is specified, just load it
 346    if (mod.EndsWith(".so")) {
 347       mod.Remove(mod.Index(".so"));
 348       result = gSystem->Load(mod);
 349       if (result < 0) {
 350          ::Error("AnalysisTrainNew.C::LoadLibrary", "Could not load library %s", module);
 351          return kFALSE;
 352       }
 353       return kTRUE;
 354    }
 355    // Check if the library is already loaded
 356    if (strlen(gSystem->GetLibraries(Form("%s.so", module), "", kFALSE)) > 0) return kTRUE;
 357    result = gSystem->Load(Form("lib%s.so", module));
 358    if (result < 0) {
 359       ::Error("AnalysisTrainNew.C::LoadLibrary", "Could not load module %s", module);
 360       return kFALSE;
 361    }
 362    return kTRUE;
 363 }
 364
 365
 366 //______________________________________________________________________________
 367 TChain *CreateChain()
 368 {
 369 // Create the input chain
 370    chain = new TChain("esdTree");
 371    if (gSystem->AccessPathName("AliESDs.root"))
 372       ::Error("AnalysisTrainNew.C::CreateChain", "File: AliESDs.root not in ./data dir");
 373    else
 374        chain->Add("AliESDs.root");
 375    if (chain->GetNtrees()) return chain;
 376    return NULL;
 377 }
 378
 379 //______________________________________________________________________________
 380 void AODmerge()
 381 {
 382 // Merging method. No staging and no terminate phase.
 383   TStopwatch timer;
 384   timer.Start();
 385   TString outputDir = "wn.xml";
 386   TString outputFiles = "EventStat_temp.root,AODQA.root,AliAOD.root,AliAOD.VertexingHF.root,";
 387   outputFiles +=        "AliAOD.Muons.root,AliAOD.Jets.root";
 388   TString mergeExcludes = "";
 389   TObjArray *list = outputFiles.Tokenize(",");
 390   TIter *iter = new TIter(list);
 391   TObjString *str;
 392   TString outputFile;
 393   Bool_t merged = kTRUE;
 394   while((str=(TObjString*)iter->Next())) {
 395     outputFile = str->GetString();
 396     // Skip already merged outputs
 397     if (!gSystem->AccessPathName(outputFile)) {
 398        printf("Output file <%s> found. Not merging again.",outputFile.Data());
 399        continue;
 400     }
 401     if (mergeExcludes.Contains(outputFile.Data())) continue;
 402     merged = AliAnalysisAlien::MergeOutput(outputFile, outputDir, 10, 0);
 403     if (!merged) {
 404        printf("ERROR: Cannot merge %s\n", outputFile.Data());
 405        return;
 406     }
 407   }
 408   // all outputs merged, validate
 409   ofstream out;
 410   out.open("outputs_valid", ios::out);
 411   out.close();
 412   timer.Print();
 413 }