replaced THXF to THX in many function prototypes
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / run.C
1 void run(Int_t runWhat, const Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_t aDebug = kFALSE, Int_t aProof = kFALSE, Bool_t mc = kTRUE, const char* option = "")
2 {
3   // runWhat options: 0 = AlidNdEtaTask
4   //                  1 = AlidNdEtaCorrectionTask
5   //
6   // aProof option: 0 no proof
7   //                1 proof with chain
8   //                2 proof with dataset
9
10   TString taskName;
11   if (runWhat == 0)
12   {
13     taskName = "AlidNdEtaTask";
14   }
15   else if (runWhat == 1)
16   {
17     taskName = "AlidNdEtaCorrectionTask";
18     if (!mc)
19     {
20       Printf("%s needs MC. Exiting...", taskName.Data());
21       return;
22     }
23   }
24   else
25   {
26     Printf("Do not know what to run. Exiting...");
27     return;
28   }
29
30   Printf("Processing task: %s", taskName.Data());
31
32   if (nRuns < 0)
33     nRuns = 1234567890;
34
35   if (aProof)
36   {
37     TProof::Open("lxb6046");
38     //gProof->SetParallel(1);
39
40     // Enable the needed package
41     if (1)
42     {
43       gProof->UploadPackage("$ALICE_ROOT/STEERBase");
44       gProof->EnablePackage("$ALICE_ROOT/STEERBase");
45       gProof->UploadPackage("$ALICE_ROOT/ESD");
46       gProof->EnablePackage("$ALICE_ROOT/ESD");
47       gProof->UploadPackage("$ALICE_ROOT/AOD");
48       gProof->EnablePackage("$ALICE_ROOT/AOD");
49       gProof->UploadPackage("$ALICE_ROOT/ANALYSIS");
50       gProof->EnablePackage("$ALICE_ROOT/ANALYSIS");
51       gProof->UploadPackage("$ALICE_ROOT/ANALYSISalice");
52       gProof->EnablePackage("$ALICE_ROOT/ANALYSISalice");
53     }
54     else
55     {
56       gProof->UploadPackage("/afs/cern.ch/alice/caf/sw/ALICE/PARs/v4-14-Release/AF-v4-14");
57       gProof->EnablePackage("/afs/cern.ch/alice/caf/sw/ALICE/PARs/v4-14-Release/AF-v4-14");
58     }
59
60     gProof->UploadPackage("$ALICE_ROOT/PWG0base");
61     gProof->EnablePackage("$ALICE_ROOT/PWG0base");
62   }
63   else
64   {
65     gSystem->Load("libVMC");
66     gSystem->Load("libTree");
67     gSystem->Load("libSTEERBase");
68     gSystem->Load("libESD");
69     gSystem->Load("libAOD");
70     gSystem->Load("libANALYSIS");
71     gSystem->Load("libANALYSISalice");
72     gSystem->Load("libPWG0base");
73   }
74
75   // Create the analysis manager
76   mgr = new AliAnalysisManager;
77
78   TString compileTaskName;
79   compileTaskName.Form("%s.cxx+", taskName.Data());
80   if (aDebug)
81     compileTaskName += "+g";
82
83   // Create, add task
84   if (aProof) {
85     gProof->Load(compileTaskName);
86   } else
87     gROOT->Macro(compileTaskName);
88
89   AliPWG0Helper::AnalysisMode analysisMode = AliPWG0Helper::kSPD;
90   AliPWG0Helper::Trigger      trigger      = AliPWG0Helper::kMB1;
91
92   AliPWG0Helper::PrintConf(analysisMode, trigger);
93
94   AliESDtrackCuts* esdTrackCuts = 0;
95   if (analysisMode != AliPWG0Helper::kSPD)
96   {
97     // selection of esd tracks
98     gROOT->ProcessLine(".L ../CreateStandardCuts.C");
99     esdTrackCuts = CreateTrackCuts(analysisMode);
100     if (!esdTrackCuts)
101     {
102       printf("ERROR: esdTrackCuts could not be created\n");
103       return;
104     }
105     esdTrackCuts->SetHistogramsOn(kTRUE);
106   }
107
108   if (runWhat == 0)
109   {
110     task = new AlidNdEtaTask(option);
111
112     if (mc)
113       task->SetReadMC();
114
115     //task->SetUseMCVertex();
116     //task->SetUseMCKine();
117   }
118   else if (runWhat == 1)
119   {
120     task = new AlidNdEtaCorrectionTask(option);
121
122     //task->SetOnlyPrimaries();
123   }
124
125   task->SetTrigger(trigger);
126   task->SetAnalysisMode(analysisMode);
127   task->SetTrackCuts(esdTrackCuts);
128
129   mgr->AddTask(task);
130
131   if (mc) {
132     // Enable MC event handler
133     AliMCEventHandler* handler = new AliMCEventHandler;
134     handler->SetReadTR(kFALSE);
135     mgr->SetMCtruthEventHandler(handler);
136   }
137
138   // Add ESD handler
139   AliESDInputHandler* esdH = new AliESDInputHandler;
140   esdH->SetInactiveBranches("*");
141   mgr->SetInputEventHandler(esdH);
142
143   // Attach input
144   cInput  = mgr->CreateContainer("cInput", TChain::Class(), AliAnalysisManager::kInputContainer);
145   mgr->ConnectInput(task, 0, cInput);
146
147   // Attach output
148   cOutput = mgr->CreateContainer("cOutput", TList::Class(), AliAnalysisManager::kOutputContainer);
149   mgr->ConnectOutput(task, 0, cOutput);
150
151   // Enable debug printouts
152   if (aDebug)
153   {
154     mgr->SetDebugLevel(2);
155     AliLog::SetClassDebugLevel(taskName, AliLog::kDebug+2);
156   }
157   else
158     AliLog::SetClassDebugLevel(taskName, AliLog::kWarning);
159
160   // Run analysis
161   mgr->InitAnalysis();
162   mgr->PrintStatus();
163
164   if (aProof == 2)
165   {
166     // process dataset
167
168     mgr->StartAnalysis("proof", data, nRuns, offset);
169   }
170   else
171   {
172     // Create chain of input files
173     gROOT->LoadMacro("../CreateESDChain.C");
174     chain = CreateESDChain(data, nRuns, offset);
175     //chain = CreateChain("TE", data, nRuns, offset);
176
177     mgr->StartAnalysis((aProof > 0) ? "proof" : "local", chain);
178   }
179 }
180
181 void loadlibs()
182 {
183   gSystem->Load("libTree");
184   gSystem->Load("libVMC");
185
186   gSystem->Load("libSTEERBase");
187   gSystem->Load("libANALYSIS");
188   gSystem->Load("libPWG0base");
189 }
190
191 void FinishAnalysisAll(const char* dataInput = "analysis_esd_raw.root", const char* dataOutput = "analysis_esd.root", const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
192 {
193   loadlibs();
194
195   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
196   TFile::Open(correctionMapFile);
197   dNdEtaCorrection->LoadHistograms();
198
199   TFile* file = TFile::Open(dataInput);
200
201   if (!file)
202   {
203     cout << "Error. File not found" << endl;
204     return;
205   }
206
207     // Note: the last parameter does not define which analysis is going to happen, the histograms will be overwritten when loading from the f
208   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaNSD", "dndetaNSD");
209   fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
210   fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3, AlidNdEtaCorrection::kNSD, "ESD -> NSD");
211   //fdNdEtaAnalysis->DrawHistograms(kTRUE);
212   TFile* file2 = TFile::Open(dataOutput, "RECREATE");
213   fdNdEtaAnalysis->SaveHistograms();
214
215   file->cd();
216   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
217   fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
218   fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3, AlidNdEtaCorrection::kINEL, "ESD -> full inelastic");
219   //fdNdEtaAnalysis->DrawHistograms(kTRUE);
220   file2->cd();
221   fdNdEtaAnalysis->SaveHistograms();
222
223   file->cd();
224   fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTr", "dndetaTr");
225   fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
226   fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3, AlidNdEtaCorrection::kVertexReco, "ESD -> minimum bias");
227   //fdNdEtaAnalysis->DrawHistograms(kTRUE);
228   file2->cd();
229   fdNdEtaAnalysis->SaveHistograms();
230
231   file->cd();
232   fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx");
233   fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
234   fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3, AlidNdEtaCorrection::kTrack2Particle, "ESD -> MB with vertex");
235   //fdNdEtaAnalysis->DrawHistograms(kTRUE);
236   file2->cd();
237   fdNdEtaAnalysis->SaveHistograms();
238
239   file->cd();
240   fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTracks", "dndetaTracks");
241   fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
242   fdNdEtaAnalysis->Finish(0, 0.3, AlidNdEtaCorrection::kNone, "ESD raw with pt cut");
243   //fdNdEtaAnalysis->DrawHistograms(kTRUE);
244   file2->cd();
245   fdNdEtaAnalysis->SaveHistograms();
246
247   file->cd();
248   fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTracksAll", "dndetaTracksAll");
249   fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
250   fdNdEtaAnalysis->Finish(0, -1, AlidNdEtaCorrection::kNone, "ESD raw");
251   //fdNdEtaAnalysis->DrawHistograms(kTRUE);
252   file2->cd();
253   fdNdEtaAnalysis->SaveHistograms();
254
255   file2->Write();
256   file2->Close();
257 }
258
259 void* FinishAnalysis(const char* analysisFile = "analysis_esd_raw.root", const char* analysisDir = "fdNdEtaAnalysisESD", const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction", Bool_t useUncorrected = kFALSE, Bool_t simple = kFALSE)
260 {
261   loadlibs();
262
263   TFile* file = TFile::Open(analysisFile);
264
265   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis(analysisDir, analysisDir);
266   fdNdEtaAnalysis->LoadHistograms();
267
268   if (correctionMapFile)
269   {
270     AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
271     TFile::Open(correctionMapFile);
272     dNdEtaCorrection->LoadHistograms();
273
274     fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3, AlidNdEtaCorrection::kINEL);
275     //fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0, AlidNdEtaCorrection::kINEL);
276     //fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0, AlidNdEtaCorrection::kTrack2Particle);
277   }
278   else
279     fdNdEtaAnalysis->Finish(0, 0.3, AlidNdEtaCorrection::kNone);
280
281   fdNdEtaAnalysis->DrawHistograms(simple);
282
283   TH1* hist = fdNdEtaAnalysis->GetdNdEtaHistogram(1);
284   Int_t binLeft = hist->GetXaxis()->FindBin(-0.5);
285   Int_t binRight = hist->GetXaxis()->FindBin(0.5);
286   Float_t value1 = hist->Integral(binLeft, binRight);
287
288   hist = fdNdEtaAnalysis->GetdNdEtaHistogram(2);
289   Float_t value2 = hist->Integral(binLeft, binRight);
290
291   if (value2 > 0)
292     printf("Ratio is %f, values are %f %f\n", value1 / value2, value1, value2);
293
294   return fdNdEtaAnalysis;
295 }