22758e10ae06eb01ec76bb5c7d6412f55f1d7172
[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     /*gProof->UploadPackage("STEERBase");
42     gProof->EnablePackage("STEERBase");
43     gProof->UploadPackage("ESD");
44     gProof->EnablePackage("ESD");
45     gProof->UploadPackage("AOD");
46     gProof->EnablePackage("AOD");
47     gProof->UploadPackage("ANALYSIS");
48     gProof->EnablePackage("ANALYSIS");*/
49
50     gProof->UploadPackage("$ALICE_ROOT/AF-v4-12");
51     gProof->EnablePackage("$ALICE_ROOT/AF-v4-12");
52
53     gProof->UploadPackage("$ALICE_ROOT/PWG0base");
54     gProof->EnablePackage("$ALICE_ROOT/PWG0base");
55   }
56   else
57   {
58     gSystem->Load("libVMC");
59     gSystem->Load("libTree");
60     gSystem->Load("libSTEERBase");
61     gSystem->Load("libESD");
62     gSystem->Load("libANALYSIS");
63     gSystem->Load("libANALYSISalice");
64     gSystem->Load("libPWG0base");
65   }
66
67   // Create the analysis manager
68   mgr = new AliAnalysisManager;
69
70   TString compileTaskName;
71   compileTaskName.Form("%s.cxx+", taskName.Data());
72   if (aDebug)
73     compileTaskName += "+g";
74
75   // Create, add task
76   if (aProof) {
77     gProof->Load(compileTaskName);
78   } else
79     gROOT->Macro(compileTaskName);
80
81   AliPWG0Helper::AnalysisMode analysisMode = AliPWG0Helper::kSPD;
82   AliPWG0Helper::Trigger      trigger = AliPWG0Helper::kMB1;
83
84   AliPWG0Helper::PrintConf(analysisMode, trigger);
85
86   AliESDtrackCuts* esdTrackCuts = 0;
87   if (analysisMode != AliPWG0Helper::kSPD)
88   {
89     // selection of esd tracks
90     gROOT->ProcessLine(".L ../CreateStandardCuts.C");
91     esdTrackCuts = CreateTrackCuts(analysisMode);
92     if (!esdTrackCuts)
93     {
94       printf("ERROR: esdTrackCuts could not be created\n");
95       return;
96     }
97   }
98
99   if (runWhat == 0)
100   {
101     task = new AlidNdEtaTask(option);
102
103     if (mc)
104       task->SetReadMC();
105
106     //task->SetUseMCVertex();
107     //task->SetUseMCKine();
108   }
109   else if (runWhat == 1)
110   {
111     task = new AlidNdEtaCorrectionTask(option);
112
113     //task->SetOnlyPrimaries();
114   }
115
116   task->SetTrigger(trigger);
117   task->SetAnalysisMode(analysisMode);
118   task->SetTrackCuts(esdTrackCuts);
119
120   mgr->AddTask(task);
121
122   if (mc) {
123     // Enable MC event handler
124     AliMCEventHandler* handler = new AliMCEventHandler;
125     handler->SetReadTR(kFALSE);
126     mgr->SetMCtruthEventHandler(handler);
127   }
128
129   // Add ESD handler
130   AliESDInputHandler* esdH = new AliESDInputHandler;
131   mgr->SetInputEventHandler(esdH);
132
133   // Attach input
134   cInput  = mgr->CreateContainer("cInput", TChain::Class(), AliAnalysisManager::kInputContainer);
135   mgr->ConnectInput(task, 0, cInput);
136
137   // Attach output
138   cOutput = mgr->CreateContainer("cOutput", TList::Class(), AliAnalysisManager::kOutputContainer);
139   mgr->ConnectOutput(task, 0, cOutput);
140
141   // Enable debug printouts
142   if (aDebug)
143   {
144     mgr->SetDebugLevel(2);
145     AliLog::SetClassDebugLevel(taskName, AliLog::kDebug+2);
146   }
147   else
148     AliLog::SetClassDebugLevel(taskName, AliLog::kWarning);
149
150   // Run analysis
151   mgr->InitAnalysis();
152   mgr->PrintStatus();
153
154   if (aProof == 2)
155   {
156     // process dataset
157
158     mgr->StartAnalysis("proof", data, nRuns, offset);
159   }
160   else
161   {
162     // Create chain of input files
163     gROOT->LoadMacro("../CreateESDChain.C");
164     chain = CreateESDChain(data, nRuns, offset);
165
166     mgr->StartAnalysis((aProof > 0) ? "proof" : "local", chain);
167   }
168 }
169
170 void loadlibs()
171 {
172   gSystem->Load("libTree");
173   gSystem->Load("libVMC");
174
175   gSystem->Load("libSTEERBase");
176   gSystem->Load("libANALYSIS");
177   gSystem->Load("libPWG0base");
178 }
179
180 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")
181 {
182   loadlibs();
183
184   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
185   TFile::Open(correctionMapFile);
186   dNdEtaCorrection->LoadHistograms();
187
188   TFile* file = TFile::Open(dataInput);
189
190   if (!file)
191   {
192     cout << "Error. File not found" << endl;
193     return;
194   }
195
196   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaNSD", "dndetaNSD");
197   fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
198   fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3, AlidNdEtaCorrection::kNSD, "ESD -> NSD");
199   //fdNdEtaAnalysis->DrawHistograms(kTRUE);
200   TFile* file2 = TFile::Open(dataOutput, "RECREATE");
201   fdNdEtaAnalysis->SaveHistograms();
202
203   file->cd();
204   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
205   fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
206   fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3, AlidNdEtaCorrection::kINEL, "ESD -> full inelastic");
207   //fdNdEtaAnalysis->DrawHistograms(kTRUE);
208   file2->cd();
209   fdNdEtaAnalysis->SaveHistograms();
210
211   file->cd();
212   fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTr", "dndetaTr");
213   fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
214   fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3, AlidNdEtaCorrection::kVertexReco, "ESD -> minimum bias");
215   //fdNdEtaAnalysis->DrawHistograms(kTRUE);
216   file2->cd();
217   fdNdEtaAnalysis->SaveHistograms();
218
219   file->cd();
220   fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx");
221   fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
222   fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3, AlidNdEtaCorrection::kTrack2Particle, "ESD -> MB with vertex");
223   //fdNdEtaAnalysis->DrawHistograms(kTRUE);
224   file2->cd();
225   fdNdEtaAnalysis->SaveHistograms();
226
227   file->cd();
228   fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTracks", "dndetaTracks");
229   fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
230   fdNdEtaAnalysis->Finish(0, 0.3, AlidNdEtaCorrection::kNone, "ESD raw");
231   //fdNdEtaAnalysis->DrawHistograms(kTRUE);
232   file2->cd();
233   fdNdEtaAnalysis->SaveHistograms();
234
235   file2->Write();
236   file2->Close();
237 }
238
239 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)
240 {
241   loadlibs();
242
243   TFile* file = TFile::Open(analysisFile);
244
245   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis(analysisDir, analysisDir);
246   fdNdEtaAnalysis->LoadHistograms();
247
248   if (correctionMapFile)
249   {
250     AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
251     TFile::Open(correctionMapFile);
252     dNdEtaCorrection->LoadHistograms();
253
254     fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3, AlidNdEtaCorrection::kINEL);
255     //fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0, AlidNdEtaCorrection::kINEL);
256     //fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0, AlidNdEtaCorrection::kTrack2Particle);
257   }
258   else
259     fdNdEtaAnalysis->Finish(0, 0.3, AlidNdEtaCorrection::kNone);
260
261   fdNdEtaAnalysis->DrawHistograms(simple);
262
263   TH1* hist = fdNdEtaAnalysis->GetdNdEtaHistogram(1);
264   Int_t binLeft = hist->GetXaxis()->FindBin(-0.5);
265   Int_t binRight = hist->GetXaxis()->FindBin(0.5);
266   Float_t value1 = hist->Integral(binLeft, binRight);
267
268   hist = fdNdEtaAnalysis->GetdNdEtaHistogram(2);
269   Float_t value2 = hist->Integral(binLeft, binRight);
270
271   if (value2 > 0)
272     printf("Ratio is %f, values are %f %f\n", value1 / value2, value1, value2);
273
274   return fdNdEtaAnalysis;
275 }