MFT classes modified to include MUON+MFT tracks in ESD&AOD. Clusterization algorithm...
[u/mrichter/AliRoot.git] / MFT / AODtrain.C
1 // ### Settings that make sense when using the Alien plugin
2 //==============================================================================
3 Int_t       runOnData          = 0;       // Set to 1 if processing real data
4 Int_t       iCollision         = 0;       // 0=pp, 1=Pb-Pb
5 //==============================================================================
6 Bool_t      usePhysicsSelection = kTRUE; // use physics selection
7 Bool_t      useTender           = kFALSE; // use tender wagon
8 Bool_t      useCentrality       = kTRUE; // centrality
9 Bool_t      useV0tender         = kFALSE;  // use V0 correction in tender
10 Bool_t      useDBG              = kTRUE;  // activate debugging
11 Bool_t      useMC               = kTRUE;  // use MC info
12 Bool_t      useKFILTER          = kTRUE;  // use Kinematics filter
13 Bool_t      useTR               = kTRUE;  // use track references
14 Bool_t      useCORRFW           = kFALSE; // do not change
15 Bool_t      useAODTAGS          = kFALSE; // use AOD tags
16 Bool_t      useSysInfo          = kFALSE; // use sys info
17
18 // ### Analysis modules to be included. Some may not be yet fully implemented.
19 //==============================================================================
20 Int_t       iAODhandler        = 1;      // Analysis produces an AOD or dAOD's
21 Int_t       iESDMCLabelAddition= 0;
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             = 0;      // Jet analysis (PWG4)
25 Int_t       iJETANdelta        = 0;      // Jet delta AODs
26 Int_t       iPWGHFvertexing     = 0;      // Vertexing HF task (PWG3)
27 Int_t       iPWGDQJPSIfilter    = 0;      // JPSI filtering (PWG3)
28 Int_t       iPWGHFd2h           = 0;      // D0->2 hadrons (PWG3)
29 Bool_t doPIDResponse  = 1;
30 Bool_t doPIDqa        = 1; //new
31
32 // ### Configuration macros used for each module
33 //==============================================================================
34  TString configPWGHFd2h = (iCollision==0)?"$ALICE_ROOT/PWGHF/vertexingHF/ConfigVertexingHF.C"
35                           :"$ALICE_ROOT/PWGHF/vertexingHF/ConfigVertexingHF_Pb_AllCent.C";
36                                                  
37 // Temporaries.
38 class AliOADBPhysicsSelection;                                                                                                                  
39 AliOADBPhysicsSelection *CreateOADBphysicsSelection();
40 void AODmerge();
41 void AddAnalysisTasks(Int_t);
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    AliLog::SetGlobalDebugLevel(3);
52   if (merge) {
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(merge);
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(Int_t merge){                                                                                                                                          
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    if (usePhysicsSelection) {
152    // Physics selection task
153       gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C");
154       mgr->RegisterExtraFile("event_stat.root");
155       AliPhysicsSelectionTask *physSelTask = AddTaskPhysicsSelection(kFALSE);
156 //      AliOADBPhysicsSelection * oadbDefaultPbPb = CreateOADBphysicsSelection();      
157 //      physSelTask->GetPhysicsSelection()->SetCustomOADBObjects(oadbDefaultPbPb,0,0);
158 //      if (!merge) mgr->AddStatisticsTask(AliVEvent::kAny);
159    }
160    // Centrality (only Pb-Pb)
161    if (iCollision && useCentrality) {
162       gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskCentrality.C");
163       AliCentralitySelectionTask *taskCentrality = AddTaskCentrality();
164       taskCentrality->SelectCollisionCandidates(AliVEvent::kAny);
165    }
166
167         if(iESDMCLabelAddition) {
168                 gROOT->LoadMacro("$ALICE_ROOT/PWG/muondep/AddTaskESDMCLabelAddition.C");
169                 AliAnalysisTaskESDMCLabelAddition *esdmclabel = AddTaskESDMCLabelAddition(kFALSE);
170         }
171         
172    if (iESDfilter) {
173       //  ESD filter task configuration.
174       gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskESDFilter.C");       
175       if (iMUONcopyAOD) {
176          printf("Registering delta AOD file\n");
177          mgr->RegisterExtraFile("AliAOD.Muons.root");
178          mgr->RegisterExtraFile("AliAOD.Dimuons.root");
179          AliAnalysisTaskESDfilter *taskesdfilter = AddTaskESDFilter(useKFILTER, kTRUE, kFALSE, kFALSE /*usePhysicsSelection*/,kFALSE,kTRUE,kTRUE,kTRUE,1100,1);
180       } else {
181     AliAnalysisTaskESDfilter *taskesdfilter = AddTaskESDFilter(useKFILTER, kFALSE, kFALSE, kFALSE /*usePhysicsSelection*/,kFALSE,kTRUE,kTRUE,kTRUE,1100,1);
182       }   
183    }   
184
185 // ********** PWG3 wagons ******************************************************           
186    // PWGHF vertexing
187    if (iPWGHFvertexing) {
188       gROOT->LoadMacro("$ALICE_ROOT/PWGHF/vertexingHF/macros/AddTaskVertexingHF.C");
189       if (!iPWGHFd2h) TFile::Cp(gSystem->ExpandPathName(configPWGHFd2h.Data()), "file:ConfigVertexingHF.C");
190       AliAnalysisTaskSEVertexingHF *taskvertexingHF = AddTaskVertexingHF();
191       if (!taskvertexingHF) ::Warning("AnalysisTrainNew", "AliAnalysisTaskSEVertexingHF cannot run for this train conditions - EXCLUDED");
192       else mgr->RegisterExtraFile("AliAOD.VertexingHF.root");
193       taskvertexingHF->SelectCollisionCandidates(0);
194    }   
195       
196    // PWGDQ JPSI filtering (only pp)
197    if (iPWGDQJPSIfilter && (iCollision==0)) {
198       gROOT->LoadMacro("$ALICE_ROOT/PWGDQ/dielectron/macros/AddTaskJPSIFilter.C");
199       AliAnalysisTaskSE *taskJPSIfilter = AddTaskJPSIFilter();
200       if (!taskJPSIfilter) ::Warning("AnalysisTrainNew", "AliAnalysisTaskDielectronFilter cannot run for this train conditions - EXCLUDED");
201       else mgr->RegisterExtraFile("AliAOD.Dielectron.root");
202       taskJPSIfilter->SelectCollisionCandidates(0);
203    }   
204
205    // PWGHF D2h
206    if (iPWGHFd2h) {   
207      gROOT->LoadMacro("$ALICE_ROOT/PWGHF/vertexingHF/AddD2HTrain.C");
208      TFile::Cp(gSystem->ExpandPathName(configPWGHFd2h.Data()), "file:ConfigVertexingHF.C");
209      AddD2HTrain(kFALSE, 1,0,0,0,0,0,0,0,0,0,0);                                 
210    }
211    
212    // ********** PWG4 wagons ******************************************************
213    // Jet analysis
214
215    // Configurations flags, move up?
216    TString kDeltaAODJetName = "AliAOD.Jets.root"; //
217    Bool_t  kIsPbPb = (iCollision==0)?false:true; // can be more intlligent checking the name of the data set
218    TString kDefaultJetBackgroundBranch = "";
219    TString kJetSubtractBranches = "";
220    UInt_t kHighPtFilterMask = 128;// from esd filter
221    UInt_t iPhysicsSelectionFlag = AliVEvent::kMB;
222    if (iJETAN) {
223      gROOT->LoadMacro("$ALICE_ROOT/PWG4/macros/AddTaskJets.C");
224      // Default jet reconstructor running on ESD's
225      AliAnalysisTaskJets *taskjets = AddTaskJets("AOD","UA1",0.4,kHighPtFilterMask,1.,0); // no background subtraction     
226      if (!taskjets) ::Fatal("AnalysisTrainNew", "AliAnalysisTaskJets cannot run for this train conditions - EXCLUDED");
227      if(kDeltaAODJetName.Length()>0) taskjets->SetNonStdOutputFile(kDeltaAODJetName.Data());
228      if (iJETANdelta) {
229         //            AddTaskJetsDelta("AliAOD.Jets.root"); // need to modify this accordingly in the add task jets
230         mgr->RegisterExtraFile(kDeltaAODJetName.Data());
231         TString cTmp("");
232         if(kIsPbPb){
233           // UA1 intrinsic background subtraction
234           taskjets = AddTaskJets("AOD","UA1",0.4,kHighPtFilterMask,1.,2); // background subtraction
235           if(kDeltaAODJetName.Length()>0)taskjets->SetNonStdOutputFile(kDeltaAODJetName.Data());
236        }
237        // SICONE 
238        taskjets = AddTaskJets("AOD","SISCONE",0.4,kHighPtFilterMask,0.15,0); //no background subtration to be done later....                                                                                  
239        if(kDeltaAODJetName.Length()>0)taskjets->SetNonStdOutputFile(kDeltaAODJetName.Data());
240        cTmp = taskjets->GetNonStdBranch();
241        if(cTmp.Length()>0)kJetSubtractBranches += Form("%s ",cTmp.Data());
242          
243        // Add the clusters..
244        gROOT->LoadMacro("$ALICE_ROOT/PWG4/macros/AddTaskJetCluster.C");
245        AliAnalysisTaskJetCluster *taskCl = 0;
246        Float_t fCenUp = 0;
247        Float_t fCenLo = 0;
248        Float_t fTrackEtaWindow = 0.9;
249        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                                                                                 
250        taskCl->SetBackgroundCalc(kTRUE);
251        taskCl->SetNRandomCones(10);
252        taskCl->SetCentralityCut(fCenLo,fCenUp);
253        taskCl->SetGhostEtamax(fTrackEtaWindow);
254        kDefaultJetBackgroundBranch = Form("%s_%s",AliAODJetEventBackground::StdBranchName(),taskCl->GetJetOutputBranch());
255
256        taskCl = AddTaskJetCluster("AOD","",kHighPtFilterMask,iPhysicsSelectionFlag,"ANTIKT",0.4,2,1,kDeltaAODJetName.Data(),0.15);
257        taskCl->SetCentralityCut(fCenLo,fCenUp);
258        if(kIsPbPb)taskCl->SetBackgroundBranch(kDefaultJetBackgroundBranch.Data());
259        taskCl->SetNRandomCones(10);
260        kJetSubtractBranches += Form("%s ",taskCl->GetJetOutputBranch());
261
262        taskCl = AddTaskJetCluster("AOD","",kHighPtFilterMask,iPhysicsSelectionFlag,"ANTIKT",0.2,0,1,kDeltaAODJetName.Data(),0.15);
263        taskCl->SetCentralityCut(fCenLo,fCenUp);
264        if(kIsPbPb)taskCl->SetBackgroundBranch(kDefaultJetBackgroundBranch.Data());
265        kJetSubtractBranches += Form("%s ",taskCl->GetJetOutputBranch());
266          
267        // DO THE BACKGROUND SUBTRACTION
268        if(kIsPbPb&&kJetSubtractBranches.Length()){
269          gROOT->LoadMacro("$ALICE_ROOT/PWG4/macros/AddTaskJetBackgroundSubtract.C");
270          AliAnalysisTaskJetBackgroundSubtract *taskSubtract = 0;
271          taskSubtract = AddTaskJetBackgroundSubtract(kJetSubtractBranches,1,"B0","B%d");
272          taskSubtract->SetBackgroundBranch(kDefaultJetBackgroundBranch.Data());
273          if(kDeltaAODJetName.Length()>0)taskSubtract->SetNonStdOutputFile(kDeltaAODJetName.Data());
274        }
275      } 
276    }
277 }
278
279 //______________________________________________________________________________
280 Bool_t LoadCommonLibraries()
281 {
282 // Load common analysis libraries.
283    if (!gSystem->Getenv("ALICE_ROOT")) {
284       ::Error("AnalysisTrainNew.C::LoadCommonLibraries", "Analysis train requires that analysis libraries are compiled with a local AliRoot"); 
285       return kFALSE;
286    }   
287    Bool_t success = kTRUE;
288    // Load framework classes. Par option ignored here.
289    success &= LoadLibrary("libSTEERBase.so");
290    success &= LoadLibrary("libESD.so");
291    success &= LoadLibrary("libAOD.so");
292    success &= LoadLibrary("libANALYSIS.so");
293    success &= LoadLibrary("libOADB.so");
294    success &= LoadLibrary("libANALYSISalice.so");
295    success &= LoadLibrary("libCORRFW.so");
296    gROOT->ProcessLine(".include $ALICE_ROOT/include");
297    if (success) {
298       ::Info("AnalysisTrainNew.C::LoadCommodLibraries", "Load common libraries:    SUCCESS");
299       ::Info("AnalysisTrainNew.C::LoadCommodLibraries", "Include path for Aclic compilation:\n%s",
300               gSystem->GetIncludePath());
301    } else {           
302       ::Info("AnalysisTrainNew.C::LoadCommodLibraries", "Load common libraries:    FAILED");
303    }   
304    return success;
305 }
306
307 //______________________________________________________________________________
308 Bool_t LoadAnalysisLibraries()
309 {
310 // Load common analysis libraries.
311    if (useTender) {
312       if (!LoadLibrary("TENDER") ||
313           !LoadLibrary("TENDERSupplies")) return kFALSE;
314    }       
315    if (iESDfilter || iPWGMuonTrain) {
316       if (!LoadLibrary("PWGmuon")) return kFALSE;
317    }   
318    if (iESDMCLabelAddition) {
319      if (!LoadLibrary("PWGmuondep")) return kFALSE;
320    }
321    // JETAN
322    if (iJETAN) {
323       if (!LoadLibrary("JETAN")) return kFALSE;
324    }
325    if (iJETANdelta) {
326       if (!LoadLibrary("JETAN") ||
327           !LoadLibrary("CGAL") ||
328           !LoadLibrary("fastjet") ||
329           !LoadLibrary("siscone") ||
330           !LoadLibrary("SISConePlugin") ||
331           !LoadLibrary("FASTJETAN")) return kFALSE;
332    }     
333    // PWG3 Vertexing HF
334    if (iPWGHFvertexing || iPWGHFd2h) {
335       if (!LoadLibrary("PWGflowBase") ||
336           !LoadLibrary("PWGflowTasks") ||
337           !LoadLibrary("PWGHFvertexingHF")) return kFALSE;
338    }    
339  //    if (iPWGHFvertexing || iPWG3d2h) {
340  //     if (!LoadLibrary("PWG3base") ||
341  //         !LoadLibrary("PWGHFvertexingHF")) return kFALSE;
342  //  }   
343    // PWG3 dielectron
344    if (iPWGDQJPSIfilter) {
345       if (!LoadLibrary("PWGDQdielectron")) return kFALSE;
346    }   
347    
348    ::Info("AnalysisTrainNew.C::LoadAnalysisLibraries", "Load other libraries:   SUCCESS");
349    return kTRUE;
350 }
351
352 //______________________________________________________________________________
353 Bool_t LoadLibrary(const char *module)
354 {
355 // Load a module library in a given mode. Reports success.
356    Int_t result;
357    TString mod(module);
358    if (!mod.Length()) {
359       ::Error("AnalysisTrainNew.C::LoadLibrary", "Empty module name");
360       return kFALSE;
361    }   
362    // If a library is specified, just load it
363    if (mod.EndsWith(".so")) {
364       mod.Remove(mod.Index(".so"));
365       result = gSystem->Load(mod);
366       if (result < 0) {
367          ::Error("AnalysisTrainNew.C::LoadLibrary", "Could not load library %s", module);
368          return kFALSE;
369       }
370       return kTRUE;
371    } 
372    // Check if the library is already loaded
373    if (strlen(gSystem->GetLibraries(Form("%s.so", module), "", kFALSE)) > 0) return kTRUE;    
374    result = gSystem->Load(Form("lib%s.so", module));
375    if (result < 0) {
376       ::Error("AnalysisTrainNew.C::LoadLibrary", "Could not load module %s", module);
377       return kFALSE;
378    }
379    return kTRUE;
380 }           
381
382
383 //______________________________________________________________________________
384 TChain *CreateChain()
385 {
386 // Create the input chain
387    chain = new TChain("esdTree");
388    if (gSystem->AccessPathName("AliESDs.root")) 
389       ::Error("AnalysisTrainNew.C::CreateChain", "File: AliESDs.root not in ./data dir");
390    else 
391        chain->Add("AliESDs.root");
392    if (chain->GetNtrees()) return chain;
393    return NULL;
394 }   
395
396 //______________________________________________________________________________
397 void AODmerge()
398 {
399 // Merging method. No staging and no terminate phase.
400   TStopwatch timer;
401   timer.Start();
402   TString outputDir = "wn.xml";
403   TString outputFiles = "EventStat_temp.root,AliAOD.root,AliAOD.Muons.root";
404   TString mergeExcludes = "";
405   TObjArray *list = outputFiles.Tokenize(",");
406   TIter *iter = new TIter(list);
407   TObjString *str;
408   TString outputFile;
409   Bool_t merged = kTRUE;
410   while((str=(TObjString*)iter->Next())) {
411     outputFile = str->GetString();
412     // Skip already merged outputs
413     if (!gSystem->AccessPathName(outputFile)) {
414        printf("Output file <%s> found. Not merging again.",outputFile.Data());
415        continue;
416     }
417     if (mergeExcludes.Contains(outputFile.Data())) continue;
418     merged = AliAnalysisAlien::MergeOutput(outputFile, outputDir, 10, 0);
419     if (!merged) {
420        printf("ERROR: Cannot merge %s\n", outputFile.Data());
421        return;
422     }
423   }
424   // all outputs merged, validate
425   ofstream out;
426   out.open("outputs_valid", ios::out);
427   out.close();
428   timer.Print();
429 }