fbdd39f88bd249010711c44587db2317d83fa6dd
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / run.C
1 void run(Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_t aDebug = kFALSE, Bool_t aProof = kFALSE, Bool_t mc = kTRUE, const char* option = "")
2 {
3   if (aProof)
4   {
5     TProof::Open("lxb6046");
6
7     // Enable the needed package
8     gProof->UploadPackage("STEERBase");
9     gProof->EnablePackage("STEERBase");
10     gProof->UploadPackage("ESD");
11     gProof->EnablePackage("ESD");
12     gProof->UploadPackage("AOD");
13     gProof->EnablePackage("AOD");
14     gProof->UploadPackage("ANALYSIS");
15     gProof->EnablePackage("ANALYSIS");
16     gProof->UploadPackage("PWG0base");
17     gProof->EnablePackage("PWG0base");
18   }
19   else
20   {
21     gSystem->Load("libVMC");
22     gSystem->Load("libTree");
23     gSystem->Load("libSTEERBase");
24     gSystem->Load("libESD");
25     gSystem->Load("libANALYSIS");
26     gSystem->Load("libPWG0base");
27   }
28
29   // Create chain of input files
30   gROOT->LoadMacro("../CreateESDChain.C");
31   chain = CreateESDChain(data, nRuns, offset);
32
33   // Create the analysis manager
34   mgr = new AliAnalysisManager;
35
36   // selection of esd tracks
37   gROOT->ProcessLine(".L CreateCuts.C");
38   AliESDtrackCuts* esdTrackCuts = CreateTrackCuts();
39   if (!esdTrackCuts)
40   {
41     printf("ERROR: esdTrackCuts could not be created\n");
42     return;
43   }
44
45   TString taskName("AlidNdEtaTask.cxx+");
46   if (aDebug)
47     taskName += "+g";
48
49   // Create, add task
50   if (aProof) {
51     gProof->Load(taskName);
52   } else
53     gROOT->Macro(taskName);
54
55   task = new AlidNdEtaTask(option);
56   task->SetTrackCuts(esdTrackCuts);
57   task->SetAnalysisMode(AlidNdEtaTask::kSPD);
58
59   if (mc)
60     task->SetReadMC();
61
62   mgr->AddTask(task);
63
64   if (mc) {
65     // Enable MC event handler
66     AliMCEventHandler* handler = new AliMCEventHandler;
67     handler->SetReadTR(kFALSE);
68     mgr->SetMCtruthEventHandler(handler);
69   }
70
71   // Add ESD handler
72   AliESDInputHandler* esdH = new AliESDInputHandler;
73   mgr->SetInputEventHandler(esdH);
74
75   // Attach input
76   cInput  = mgr->CreateContainer("cInput", TChain::Class(), AliAnalysisManager::kInputContainer);
77   mgr->ConnectInput(task, 0, cInput);
78
79   // Attach output
80   cOutput = mgr->CreateContainer("cOutput", TList::Class(), AliAnalysisManager::kOutputContainer);
81   mgr->ConnectOutput(task, 0, cOutput);
82
83   // Enable debug printouts
84   if (aDebug)
85     mgr->SetDebugLevel(2);
86
87   // Run analysis
88   mgr->InitAnalysis();
89   mgr->PrintStatus();
90
91   mgr->StartAnalysis((aProof) ? "proof" : "local", chain);
92 }
93
94 void loadlibs()
95 {
96   gSystem->Load("libTree");
97   gSystem->Load("libVMC");
98
99   gSystem->Load("libSTEERBase");
100   gSystem->Load("libANALYSIS");
101   gSystem->Load("libPWG0base");
102 }
103
104 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")
105 {
106   loadlibs();
107
108   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
109   TFile::Open(correctionMapFile);
110   dNdEtaCorrection->LoadHistograms();
111
112   TFile* file = TFile::Open(dataInput);
113
114   if (!file)
115   {
116     cout << "Error. File not found" << endl;
117     return;
118   }
119
120   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
121   fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
122   fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3, AlidNdEtaCorrection::kINEL);
123   fdNdEtaAnalysis->DrawHistograms(kTRUE);
124   TFile* file2 = TFile::Open(dataOutput, "RECREATE");
125   fdNdEtaAnalysis->SaveHistograms();
126
127   file->cd();
128   fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTr", "dndetaTr");
129   fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
130   fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3, AlidNdEtaCorrection::kVertexReco);
131   fdNdEtaAnalysis->DrawHistograms(kTRUE);
132   file2->cd();
133   fdNdEtaAnalysis->SaveHistograms();
134
135   file->cd();
136   fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx");
137   fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
138   fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3, AlidNdEtaCorrection::kTrack2Particle);
139   fdNdEtaAnalysis->DrawHistograms(kTRUE);
140   file2->cd();
141   fdNdEtaAnalysis->SaveHistograms();
142
143   file->cd();
144   fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTracks", "dndetaTracks");
145   fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
146   fdNdEtaAnalysis->Finish(0, 0.3, AlidNdEtaCorrection::kNone);
147   fdNdEtaAnalysis->DrawHistograms(kTRUE);
148   file2->cd();
149   fdNdEtaAnalysis->SaveHistograms();
150 }
151
152 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)
153 {
154   loadlibs();
155
156   TFile* file = TFile::Open(analysisFile);
157
158   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis(analysisDir, analysisDir);
159   fdNdEtaAnalysis->LoadHistograms();
160
161   if (correctionMapFile)
162   {
163     AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
164     TFile::Open(correctionMapFile);
165     dNdEtaCorrection->LoadHistograms();
166
167     fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3, AlidNdEtaCorrection::kINEL);
168     //fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0, AlidNdEtaCorrection::kINEL);
169     //fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0, AlidNdEtaCorrection::kTrack2Particle);
170   }
171   else
172     fdNdEtaAnalysis->Finish(0, 0, AlidNdEtaCorrection::kNone);
173
174   fdNdEtaAnalysis->DrawHistograms(simple);
175
176   TH1* hist = fdNdEtaAnalysis->GetdNdEtaHistogram(1);
177   Int_t binLeft = hist->GetXaxis()->FindBin(-0.5);
178   Int_t binRight = hist->GetXaxis()->FindBin(0.5);
179   Float_t value1 = hist->Integral(binLeft, binRight);
180
181   hist = fdNdEtaAnalysis->GetdNdEtaHistogram(2);
182   Float_t value2 = hist->Integral(binLeft, binRight);
183
184   if (value2 > 0)
185     printf("Ratio is %f, values are %f %f\n", value1 / value2, value1, value2);
186
187   return fdNdEtaAnalysis;
188 }