Enabled vertexing QA
[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         = 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*/,kFALSE,AliESDpid::kTOF_T0,kTRUE,kFALSE,kFALSE,run_flag);
171       } else {
172            AliAnalysisTaskESDfilter *taskesdfilter = AddTaskESDFilter(useKFILTER, kFALSE, kFALSE, kFALSE /*usePhysicsSelection*/,kFALSE,AliESDpid::kTOF_T0,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 }