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