]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG2/FORWARD/analysis/AliFMDAnalysisTaskBFCorrelation.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Minor fixes to analysis code
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis / AliFMDAnalysisTaskBFCorrelation.cxx
CommitLineData
1c038f67 1#include <TROOT.h>
2#include <TSystem.h>
3#include <TInterpreter.h>
4#include <TChain.h>
5#include <TFile.h>
6#include <TList.h>
7#include <iostream>
4b8bdb60 8#include <string>
9#include <TCanvas.h>
1c038f67 10#include "TH1F.h"
11#include "TH2F.h"
4b8bdb60 12#include "TProfile.h"
1c038f67 13#include "AliFMDAnalysisTaskBFCorrelation.h"
14#include "AliAnalysisManager.h"
15#include "AliESDFMD.h"
16#include "AliESDEvent.h"
17#include "AliAODEvent.h"
18#include "AliAODHandler.h"
19#include "AliMCEventHandler.h"
20#include "AliStack.h"
21#include "AliLog.h"
22#include "AliESDVertex.h"
23#include "TMath.h"
24#include "AliFMDAnaParameters.h"
25//#include "AliFMDGeometry.h"
26#include "AliGenEventHeader.h"
27#include "AliGenPythiaEventHeader.h"
28#include "AliHeader.h"
29//#include "TDatabasePDG.h"
30//#include "TParticlePDG.h"
31#include "AliESDInputHandler.h"
1c038f67 32
4b8bdb60 33ClassImp(AliFMDAnalysisTaskBFCorrelation)
1c038f67 34
35AliFMDAnalysisTaskBFCorrelation::AliFMDAnalysisTaskBFCorrelation()
36: fDebug(0),
37 fOutputList(0),
38 fInputList(0),
4b8bdb60 39 fInternalList(0),
1c038f67 40 fVertexString(0x0),
4b8bdb60 41 fStandalone(kTRUE),
42 fEvent(0),
2869f0ae 43 fnBinsX(0),
44 fXmin(0),
45 fXmax(0),
46 fnBinsY(0),
4b8bdb60 47 fYmin(0),
2869f0ae 48 fYmax(0)
49
1c038f67 50{
51 // Default constructor
52 DefineInput (0, TList::Class());
53 DefineOutput(0, TList::Class());
54}
4b8bdb60 55
1c038f67 56//_____________________________________________________________________
57AliFMDAnalysisTaskBFCorrelation::AliFMDAnalysisTaskBFCorrelation(const char* name, Bool_t SE):
58 AliAnalysisTask(name,name),
4b8bdb60 59 fDebug(0),
60 fOutputList(0),
61 fInputList(0),
62 fInternalList(0),
63 fVertexString(0x0),
64 fStandalone(kTRUE),
65 fEvent(0),
2869f0ae 66 fnBinsX(0),
67 fXmin(0),
68 fXmax(0),
69 fnBinsY(0),
4b8bdb60 70 fYmin(0),
2869f0ae 71 fYmax(0)
1c038f67 72{
73 fStandalone = SE;
74 if(fStandalone) {
75 DefineInput (0, TList::Class());
76 DefineInput(1, TObjString::Class());
77 DefineOutput(0, TList::Class());
1c038f67 78 }
79}
4b8bdb60 80
1c038f67 81//_____________________________________________________________________
82void AliFMDAnalysisTaskBFCorrelation::CreateOutputObjects()
83{
2869f0ae 84
85 // Nomenclature for naming histograms
86 // SE - Used in histogram names if used for containing Single Event information
87 // Rebin - Used in histogram names if they have been rebinned for analysis purposes
88
4b8bdb60 89 // Setup the list for storing results, if it does not exist
90
1c038f67 91 if(!fOutputList) {
92 fOutputList = new TList();
93 fOutputList->SetName("BackgroundCorrected");
94 }
4b8bdb60 95
96 // Setup the list for temporary storage during the analysis
97
98 if(!fInternalList) {
99 fInternalList = new TList();
100 fInternalList->SetName("InternalBFList");
101 }
102
2869f0ae 103 // Get the bounds for the input histogram.
104 // This version is optimized for 200 bins in the eta range (-4, 6)
105
106 AliFMDAnaParameters *pars = AliFMDAnaParameters::Instance();
107 TH2F *hDefault = (TH2F*)pars->GetBackgroundCorrection(1, 'I', 1);
108
109 fnBinsX = hDefault->GetNbinsX();
110 fnBinsY = hDefault->GetNbinsY();
111 fXmin = hDefault->GetXaxis()->GetXmin();
112 fXmax = hDefault->GetXaxis()->GetXmax();
113 fYmin = hDefault->GetYaxis()->GetXmin();
114 fYmax = hDefault->GetYaxis()->GetXmax();
115
116 // Histogram to contain an event of MC data. Same dimension as ESD event
117
118 TH2F *hEtaPhiParticleMap = new TH2F("hEtaPhiParticleMap",
119 "Distributions of MC particles (Eta,Phi)",
120 fnBinsX, fXmin, fXmax,
121 fnBinsY, fYmin, fYmax);
122 hEtaPhiParticleMap->Sumw2();
123
124 fInternalList->Add(hEtaPhiParticleMap);
125
126 // Create histograms with same binning as input for Response analysis
127 // and control histograms. One temporary histogram is also created to
128 // avoid the new and delete call. Must be reset after use !
129
130 TH1D *hSEMultESD = new TH1D("hSEMultESD", "Multiplicity", fnBinsX, fXmin, fXmax);
131 TH1D *hSEMultMC = new TH1D("hSEMultMC", "Multiplicity", fnBinsX, fXmin, fXmax);
132
133 TH1D *hTemp = new TH1D("hTemp", "Temporary histogram", fnBinsX, fXmin, fXmax);
134
135 fInternalList->Add(hSEMultESD);
136 fInternalList->Add(hSEMultMC);
137 fInternalList->Add(hTemp);
138
139 TH1F *hMultvsEtaESD = new TH1F("hMultvsEtaESD", "Multiplicity vs Eta (ESD)", fnBinsX, fXmin, fXmax);
140 TH1F *hMultvsEtaMC = new TH1F("hMultvsEtaMC", "Multiplicity vs Eta (MC)", fnBinsX, fXmin, fXmax);
1c038f67 141
2869f0ae 142 TH2F *hResponseMatrix = new TH2F("hResponseMatrix", "Response matrix", 151, -0.5, 150.5, 151, -0.5, 150.5);
143
144 fOutputList->Add(hMultvsEtaESD);
145 fOutputList->Add(hMultvsEtaMC);
146 fOutputList->Add(hResponseMatrix);
147
148 // Set up histograms for analysis and storage. 3 different binnings
149
4b8bdb60 150 for (Int_t i = 1; i <= 4; i++) {
2869f0ae 151 if (i == 3) continue;
152
153 Int_t nBinsX = fnBinsX/(5*i);
154
155 // Histograms for the event-by-event analysis
156
157 TH1D *hSERebinMultESD = new TH1D(Form("hSERebinMultESD_binning%d", i),
158 Form("Multiplicity per event vs eta ESD (%d bins)", nBinsX),
159 nBinsX, fXmin, fXmax);
160
161 TH1D *hSERebinMultMC = new TH1D(Form("hSERebinMultMC_binning%d", i),
162 Form("Multiplicity per event vs eta MC-truth (%d bins)", nBinsX),
163 nBinsX, fXmin, fXmax);
164
165 TH1D *hSERebinMultMirrorESD = new TH1D(Form("hSERebinMultMirrorESD_binning%d", i),
166 Form("Multiplicity per event vs eta ESD Mirrored (%d bins)", nBinsX),
167 nBinsX, fXmin, fXmax);
168
169 TH1D *hSERebinMultMirrorMC = new TH1D(Form("hSERebinMultMirrorMC_binning%d", i),
170 Form("Multiplicity per event vs eta MC-truth Mirrored (%d bins)", nBinsX),
171 nBinsX, fXmin, fXmax);
4b8bdb60 172
2869f0ae 173
174 TH1D *hSERebinMultWESD = new TH1D(Form("hSERebinMultWESD_binning%d", i),
175 Form("Multiplicity per event vs eta ESD (%d bins)", nBinsX),
176 nBinsX, fXmin, fXmax);
177
178 TH1D *hSERebinMultWMC = new TH1D(Form("hSERebinMultWMC_binning%d", i),
179 Form("Multiplicity per event vs eta MC-truth (%d bins)", nBinsX),
180 nBinsX, fXmin, fXmax);
181
182 TH1D *hSERebinMultMirrorWESD = new TH1D(Form("hSERebinMultMirrorWESD_binning%d", i),
183 Form("Multiplicity per event vs eta ESD Mirrored (%d bins)", nBinsX),
184 nBinsX, fXmin, fXmax);
185
186 TH1D *hSERebinMultMirrorWMC = new TH1D(Form("hSERebinMultMirrorWMC_binning%d", i),
187 Form("Multiplicity per event vs eta MC-truth Mirrored (%d bins)", nBinsX),
188 nBinsX, fXmin, fXmax);
189
190 fInternalList->Add(hSERebinMultESD);
191 fInternalList->Add(hSERebinMultMC);
192 fInternalList->Add(hSERebinMultMirrorESD);
193 fInternalList->Add(hSERebinMultMirrorMC);
194
195 fInternalList->Add(hSERebinMultWESD);
196 fInternalList->Add(hSERebinMultWMC);
197 fInternalList->Add(hSERebinMultMirrorWESD);
198 fInternalList->Add(hSERebinMultMirrorWMC);
199
200 // Histograms for storing the acummulated parameters.
201
202 TH1F *hRebinnESD = new TH1F(Form("hRebinnESD_binning%d", i),
203 Form("Counts ESD (%d bins)", nBinsX),
204 nBinsX, fXmin, fXmax);
205
206 TH1F *hRebinnMirrorESD = new TH1F(Form("hRebinnMirrorESD_binning%d", i),
207 Form("Counts ESD Mirrored (%d bins)", nBinsX),
208 nBinsX, fXmin, fXmax);
209
210 TH1F *hRebinn2ESD = new TH1F(Form("hRebinn2ESD_binning%d", i),
211 Form("Counts^2 ESD (%d bins)", nBinsX),
212 nBinsX, fXmin, fXmax);
213
214 TH1F *hRebinn2MirrorESD = new TH1F(Form("hRebinn2MirrorESD_binning%d", i),
215 Form("Counts^2 ESD Mirrored (%d bins)", nBinsX),
216 nBinsX, fXmin, fXmax);
217
218 TH1F *hRebinnfbESD = new TH1F(Form("hRebinnfbESD_binning%d", i),
219 Form("Fwd*bwd ESD (%d bins)", nBinsX),
220 nBinsX, fXmin, fXmax);
221
222
223 TH1F *hRebinnMC = new TH1F(Form("hRebinnMC_binning%d", i),
224 Form("Counts MC (%d bins)", nBinsX),
225 nBinsX, fXmin, fXmax);
226
227 TH1F *hRebinnMirrorMC = new TH1F(Form("hRebinnMirrorMC_binning%d", i),
228 Form("Counts MC Mirrored (%d bins)", nBinsX),
229 nBinsX, fXmin, fXmax);
230
231 TH1F *hRebinn2MC = new TH1F(Form("hRebinn2MC_binning%d", i),
232 Form("Counts^2 MC (%d bins)", nBinsX),
233 nBinsX, fXmin, fXmax);
234
235 TH1F *hRebinn2MirrorMC = new TH1F(Form("hRebinn2MirrorMC_binning%d", i),
236 Form("Counts^2 MC Mirrored (%d bins)", nBinsX),
237 nBinsX, fXmin, fXmax);
238
239 TH1F *hRebinnfbMC = new TH1F(Form("hRebinnfbMC_binning%d", i),
240 Form("Fwd*bwd MC (%d bins)", nBinsX),
241 nBinsX, fXmin, fXmax);
242
243 fOutputList->Add(hRebinnESD);
244 fOutputList->Add(hRebinnMirrorESD);
245 fOutputList->Add(hRebinn2ESD);
246 fOutputList->Add(hRebinn2MirrorESD);
247 fOutputList->Add(hRebinnfbESD);
248
249 fOutputList->Add(hRebinnMC);
250 fOutputList->Add(hRebinnMirrorMC);
251 fOutputList->Add(hRebinn2MC);
252 fOutputList->Add(hRebinn2MirrorMC);
253 fOutputList->Add(hRebinnfbMC);
254
255 // Histograms for storing the weights for the acummulated parameters.
256
257 TH1F *hRebinnWESD = new TH1F(Form("hRebinnWESD_binning%d", i),
258 Form("Counts ESD Weights (%d bins)", nBinsX),
259 nBinsX, fXmin, fXmax);
260
261 TH1F *hRebinnMirrorWESD = new TH1F(Form("hRebinnMirrorWESD_binning%d", i),
262 Form("Counts ESD Weights (%d bins)", nBinsX),
263 nBinsX, fXmin, fXmax);
264
265 TH1F *hRebinn2WESD = new TH1F(Form("hRebinn2WESD_binning%d", i),
266 Form("Counts^2 ESD Weights (%d bins)", nBinsX),
267 nBinsX, fXmin, fXmax);
268
269 TH1F *hRebinn2MirrorWESD = new TH1F(Form("hRebinn2MirrorWESD_binning%d", i),
270 Form("Counts^2 ESD Weights (%d bins)", nBinsX),
271 nBinsX, fXmin, fXmax);
272
273 TH1F *hRebinnfbWESD = new TH1F(Form("hRebinnfbWESD_binning%d", i),
274 Form("Fwd*bwd ESD Weights (%d bins)", nBinsX),
275 nBinsX, fXmin, fXmax);
276
277
278 TH1F *hRebinnWMC = new TH1F(Form("hRebinnWMC_binning%d", i),
279 Form("Counts MC weights (%d bins)", nBinsX),
280 nBinsX, fXmin, fXmax);
281
282 TH1F *hRebinnMirrorWMC = new TH1F(Form("hRebinnMirrorWMC_binning%d", i),
283 Form("Counts MC weights (%d bins)", nBinsX),
284 nBinsX, fXmin, fXmax);
285
286 TH1F *hRebinn2WMC = new TH1F(Form("hRebinn2WMC_binning%d", i),
287 Form("Counts^2 MC weights (%d bins)", nBinsX),
288 nBinsX, fXmin, fXmax);
289
290 TH1F *hRebinn2MirrorWMC = new TH1F(Form("hRebinn2MirrorWMC_binning%d", i),
291 Form("Counts^2 MC weights (%d bins)", nBinsX),
292 nBinsX, fXmin, fXmax);
293
294 TH1F *hRebinnfbWMC = new TH1F(Form("hRebinnfbWMC_binning%d", i),
295 Form("Fwd*bwd MC weights (%d bins)", nBinsX),
296 nBinsX, fXmin, fXmax);
297
298 fOutputList->Add(hRebinnWESD);
299 fOutputList->Add(hRebinnMirrorWESD);
300 fOutputList->Add(hRebinn2WESD);
301 fOutputList->Add(hRebinn2MirrorWESD);
302 fOutputList->Add(hRebinnfbWESD);
303
304 fOutputList->Add(hRebinnWMC);
305 fOutputList->Add(hRebinnMirrorWMC);
306 fOutputList->Add(hRebinn2WMC);
307 fOutputList->Add(hRebinn2MirrorWMC);
308 fOutputList->Add(hRebinnfbWMC);
309
310 // Histograms for the final result
311 /*
312 TH1D *hBFFcor = new TH1D(Form("hBFFcor_binning%d", i),
313 Form("Forward - backward correlations F (%d bins)", nBinsXBF),
314 nBinsXBF, 0, fXmax);
315 hBFFcor->GetXaxis()->SetTitle("#eta");
316 hBFFcor->GetXaxis()->SetTitle("b");
317 TH1D *hBFBcor = new TH1D(Form("hBFBcor_binning%d", i),
318 Form("Forward - backward correlations B (%d bins)", nBinsXBF),
319 nBinsXBF, 0, fXmax);
320 hBFBcor->GetXaxis()->SetTitle("#eta");
321 hBFBcor->GetXaxis()->SetTitle("b");
322
323 TH1D *hBFFcor_MC = new TH1D(Form("hBFFcor_MC_binning%d", i),
324 Form("Forward - backward correlations F (%d bins) MC", nBinsXBF),
325 nBinsXBF, 0, fXmax);
326 hBFFcor_MC->GetXaxis()->SetTitle("#eta");
327 hBFFcor_MC->GetXaxis()->SetTitle("b");
328 TH1D *hBFBcor_MC = new TH1D(Form("hBFBcor_MC_binning%d", i),
329 Form("Forward - backward correlations B (%d bins) MC", nBinsXBF),
330 nBinsXBF, 0, fXmax);
331 hBFBcor_MC->GetXaxis()->SetTitle("#eta");
332 hBFBcor_MC->GetXaxis()->SetTitle("b");
333
334 fOutputList->Add(hBFFcor);
335 fOutputList->Add(hBFBcor);
336 fOutputList->Add(hBFFcor_MC);
337 fOutputList->Add(hBFBcor_MC);
338 */
339
340 // Temporary histogram to avoid new-delete
341
342 TH1D* hRebinTemp = new TH1D(Form("hRebinTemp_binning%d", i),
343 Form("Temporary histogram (%d bins)", nBinsX),
344 nBinsX, fXmin, fXmax);
345
346 fInternalList->Add(hRebinTemp);
4b8bdb60 347 }
1c038f67 348}
4b8bdb60 349
1c038f67 350//_____________________________________________________________________
351void AliFMDAnalysisTaskBFCorrelation::ConnectInputData(Option_t */*option*/)
352{
353 if(fStandalone) {
354 fInputList = (TList*)GetInputData(0);
355 fVertexString = (TObjString*)GetInputData(1);
356 }
357}
4b8bdb60 358
1c038f67 359//_____________________________________________________________________
2869f0ae 360void AliFMDAnalysisTaskBFCorrelation::Exec(Option_t */*option*/) {
361
4b8bdb60 362 fEvent++;
d5ce0c32 363 //if (fEvent % 1000 == 0)
364 // std::cout << "Event # " << fEvent << std::endl;
4b8bdb60 365
1c038f67 366 AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
367
368 fVertexString = (TObjString*)fInputList->At(0);
2869f0ae 369
370 //Int_t vtxbin = fVertexString->GetString().Atoi();
371 //if (vtxbin != 5) return;
372
373 ProjectAndMirror("ESD");
4b8bdb60 374 CalculateValues("ESD");
1c038f67 375
2869f0ae 376 if(pars->GetProcessPrimary()) {
377
1c038f67 378 ProcessPrimary();
2869f0ae 379
380 CreateResponseMatrix();
381 }
382
1c038f67 383
384 if(fStandalone) {
385 PostData(0, fOutputList);
386 }
4b8bdb60 387}
4b8bdb60 388
2869f0ae 389//_____________________________________________________________________
390void AliFMDAnalysisTaskBFCorrelation::ProjectAndMirror(TString sType) {
391
392 sType.ToUpper();
393
394 if (!sType.Contains("ESD") && !sType.Contains("MC")) {
395 std::cout << "Wrong type specification for 'ProjectAndMirror'" << std::endl;
396 return;
4b8bdb60 397 }
398
2869f0ae 399 // Get Single Event histograms for storing hits without rebinning
400
401 TH1D *hMult = dynamic_cast<TH1D*>(fInternalList->FindObject(Form("hSEMult%s", sType.Data())));
402 hMult->Reset();
403
404 // Create generic names for retrieving histograms
4b8bdb60 405
2869f0ae 406 TList *list = 0; // List for getting either ESD or MC "hit map"
4b8bdb60 407
2869f0ae 408 TString sEtaPhiMap;
4b8bdb60 409
2869f0ae 410 if (sType.Contains("ESD")) {
4b8bdb60 411
2869f0ae 412 list = fInputList;
4b8bdb60 413
2869f0ae 414 sEtaPhiMap = "dNdetadphiHistogramSPDTrVtx";
415 }
416
417 if (sType.Contains("MC")) {
418
419 list = fInternalList;
420
421 sEtaPhiMap = "hEtaPhiParticleMap";
4b8bdb60 422 }
423
2869f0ae 424 TString sMult("hSERebinMult");
425 TString sMultMirror("hSERebinMultMirror");
426 TString sMultW("hSERebinMultW");
427 TString sMultMirrorW("hSERebinMultMirrorW");
428
429 sType += "_binning";
430
431 sMult += sType;
432 sMultMirror += sType;
433 sMultW += sType;
434 sMultMirrorW += sType;
1c038f67 435
2869f0ae 436 // Get the 2D histogram containing the particles of the event being analyzed
437
438 TH2F *hEtaPhiMap = dynamic_cast<TH2F*>(list->FindObject(sEtaPhiMap));
439
440 TH1D *hProjection = hEtaPhiMap->ProjectionX("hTemporary");
441 hMult->Add(hProjection);
442
443 // Loop over the 3 binnings
444
445 for (Int_t i = 1; i<=4; i++) {
446 if (i == 3) continue;
447
448 // Project the 2D "hit map" onto the eta-axis and rebin
449
450 TH1D *hProjRebin = hEtaPhiMap->ProjectionX("hProjRebin");
451 hProjRebin->Rebin(i*5);
452
453 // Retrieve the histograms to store the Singe Event information and reset
454
1b418b63 455 TH1D *hSEMult = static_cast<TH1D*>(fInternalList->FindObject(Form("%s%d", sMult.Data(), i)));
2869f0ae 456 hSEMult->Reset();
457
1b418b63 458 TH1D *hSEMultMirror = static_cast<TH1D*>(fInternalList->FindObject(Form("%s%d", sMultMirror.Data(), i)));
2869f0ae 459 hSEMultMirror->Reset();
460
1b418b63 461 TH1D *hSEMultW = static_cast<TH1D*>(fInternalList->FindObject(Form("%s%d", sMultW.Data(), i)));
2869f0ae 462 hSEMultW->Reset();
463
1b418b63 464 TH1D *hSEMultMirrorW = static_cast<TH1D*>(fInternalList->FindObject(Form("%s%d", sMultMirrorW.Data(), i)));
2869f0ae 465 hSEMultMirrorW->Reset();
466
467 // Fill the histograms with the Single Event information
468
469 hSEMult->Add(hProjRebin);
470
471 for (Int_t bin = 1; bin <= hSEMult->GetNbinsX(); bin++) {
472
473 hSEMultMirror->SetBinContent(hSEMultMirror->FindBin(-hProjRebin->GetBinCenter(bin)), hProjRebin->GetBinContent(bin));
474 hSEMultMirror->SetBinError(hSEMultMirror->FindBin(-hProjRebin->GetBinCenter(bin)), hProjRebin->GetBinError(bin));
475
476 // hMultDist->Fill(bin, hMultTemp->GetBinContent(bin));
477 }
478 hProjRebin->Delete();
479 }
1c038f67 480}
4b8bdb60 481
2869f0ae 482//_____________________________________________________________________
483void AliFMDAnalysisTaskBFCorrelation::CalculateValues(TString sType) {
484
485 sType.ToUpper();
486
487 if (!sType.Contains("ESD") && !sType.Contains("MC")) {
488 std::cout << "Wrong type specification for 'CalculateValues'" << std::endl;
489 return;
490 }
491
492 TString sMult("hSERebinMult");
493 TString sMultMirror("hSERebinMultMirror");
494 // TString_t *sMultW("hSEMultW");
495 // TString_t *sMultMirrorW("hSEMultMirrorW");
496
497 TString sn("hRebinn");
498 TString snMirror("hRebinnMirror");
499 TString sn2("hRebinn2");
500 TString sn2Mirror("hRebinn2Mirror");
501 TString snfb("hRebinnfb");
4b8bdb60 502
2869f0ae 503 sType += "_binning";
4b8bdb60 504
2869f0ae 505 sMult += sType;
506 sMultMirror += sType;
4b8bdb60 507
2869f0ae 508 sn += sType;
509 snMirror += sType;
510 sn2 += sType;
511 sn2Mirror += sType;
512 snfb += sType;
513
4b8bdb60 514 for (Int_t i = 1; i <= 4; i++) {
2869f0ae 515 if (i == 3) continue;
4b8bdb60 516
2869f0ae 517 TH1D *hSEMult = (TH1D*)fInternalList->FindObject(Form("%s%d", sMult.Data(), i));
518 TH1D *hSEMultMirror = (TH1D*)fInternalList->FindObject(Form("%s%d", sMultMirror.Data(), i));
519 /*
520 TH1D *hSEMultW = (TH1D*)fInternalList->FindObject(Form("%s%d", cMultW, cType, i));
521 TH1D *hSEMultMirrorW = (TH1D*)fInternalList->FindObject(Form("%s%d", cMultMirrorW, cType, i));
522 */
4b8bdb60 523
2869f0ae 524 TH1F *hn = (TH1F*)fOutputList->FindObject(Form("%s%d", sn.Data(), i));
525 TH1F *hnMirror = (TH1F*)fOutputList->FindObject(Form("%s%d", snMirror.Data(), i));
526 TH1F *hn2 = (TH1F*)fOutputList->FindObject(Form("%s%d", sn2.Data(), i));
527 TH1F *hn2Mirror = (TH1F*)fOutputList->FindObject(Form("%s%d", sn2Mirror.Data(), i));
528 TH1F *hnfb = (TH1F*)fOutputList->FindObject(Form("%s%d", snfb.Data(), i));
529 /*
530 TH1F *hnW = (TH1F*)fOutputList->FindObject(Form("%sW%s_binning%d", cn, cType, i));
531 TH1F *hnMirrorW = (TH1F*)fOutputList->FindObject(Form("%sW%s_binning%d", cnMirror, cType, i));
532 TH1F *hn2W = (TH1F*)fOutputList->FindObject(Form("%sW%s_binning%d", cn2, cType, i));
533 TH1F *hn2MirrorW = (TH1F*)fOutputList->FindObject(Form("%sW%s_binning%d", cn2Mirror, cType, i));
534 TH1F *hnfbW = (TH1F*)fOutputList->FindObject(Form("%sW%s_binning%d", cnfb, cType, i));
535 */
536 TH1D *hTemp = (TH1D*)fInternalList->FindObject(Form("hRebinTemp_binning%d", i));
537
538 hn->Add(hSEMult);
539
540 hnMirror->Add(hSEMultMirror);
541
542 hTemp->Reset();
543 hTemp->Add(hSEMult);
544 hTemp->Multiply(hSEMult);
545 hn2->Add(hTemp);
546
547 hTemp->Reset();
548 hTemp->Add(hSEMultMirror);
549 hTemp->Multiply(hSEMultMirror);
550 hn2Mirror->Add(hTemp);
551
552 hSEMultMirror->Multiply(hSEMult);
553 hnfb->Add(hSEMultMirror);
4b8bdb60 554 }
555}
556
557//_____________________________________________________________________
2869f0ae 558void AliFMDAnalysisTaskBFCorrelation::MultiplicityVsEta(TString sType) {
559
560 sType.ToUpper();
4b8bdb60 561
2869f0ae 562 if (!sType.Contains("ESD") && !sType.Contains("MC")) {
563 std::cout << "Wrong type specification for 'MultiplicityVsEta'" << std::endl;
564 return;
565 }
566}
1c038f67 567
4b8bdb60 568
569//_____________________________________________________________________
2869f0ae 570void AliFMDAnalysisTaskBFCorrelation::CreateResponseMatrix() {
571
572 TH2F *hResponseMatrix = (TH2F*)fOutputList->FindObject("hResponseMatrix");
573
574 TH1D *hSEMultESD = (TH1D*)fInternalList->FindObject("hSEMultESD");
575 TH1D *hSEMultMC = (TH1D*)fInternalList->FindObject("hSEMultMC");
1c038f67 576
2869f0ae 577 Float_t HitsESD = 0;
578 Float_t HitsMC = 0;
579
580 for (Int_t bin = 1; bin<= fnBinsX; bin++) {
581
582 if ((hSEMultMC->GetBinLowEdge(bin) > -3.4) &&
583 (hSEMultMC->GetBinLowEdge(bin+1) < 5.)) {
584
585 HitsESD += hSEMultESD->GetBinContent(bin);
586 HitsMC += hSEMultMC->GetBinContent(bin);
587 }
588 }
589
590 hResponseMatrix->Fill(HitsMC, HitsESD);
1c038f67 591}
2869f0ae 592
593//_____________________________________________________________________
594void AliFMDAnalysisTaskBFCorrelation::Terminate(Option_t */*option*/) {
595 /*
596 std::cout << "Terminating !" << std::endl;
597
598 TH1I *hnEvents = (TH1I*)fOutputList->FindObject("nEvents");
599 TH1I *hnMCEvents = (TH1I*)fOutputList->FindObject("nEvents");
600
601 Int_t nEvents = hnEvents->GetEntries();
602 Int_t nMCEvents = hnMCEvents->GetEntries();
603
604 for (Int_t i = 1; i <= 4; i++) {
605 if (i == 3) continue;
606
607 TH1F *hnFnB = (TH1F*)fOutputList->FindObject(Form("hnFnB_binning%d", i));
608 TH1F *hnF2 = (TH1F*)fOutputList->FindObject(Form("hnF2_binning%d", i));
609 TH1F *hnB2 = (TH1F*)fOutputList->FindObject(Form("hnB2_binning%d", i));
610 TH1F *hnF = (TH1F*)fOutputList->FindObject(Form("hnF_binning%d", i));
611 TH1F *hnB = (TH1F*)fOutputList->FindObject(Form("hnB_binning%d", i));
612
613 TH1F *hnFnB_MC = (TH1F*)fOutputList->FindObject(Form("hnFnB_MC_binning%d", i));
614 TH1F *hnF2_MC = (TH1F*)fOutputList->FindObject(Form("hnF2_MC_binning%d", i));
615 TH1F *hnB2_MC = (TH1F*)fOutputList->FindObject(Form("hnB2_MC_binning%d", i));
616 TH1F *hnF_MC = (TH1F*)fOutputList->FindObject(Form("hnF_MC_binning%d", i));
617 TH1F *hnB_MC = (TH1F*)fOutputList->FindObject(Form("hnB_MC_binning%d", i));
618
619 hnFnB->Scale(1./Float_t(nEvents));
620 hnF2->Scale(1./Float_t(nEvents));
621 hnB2->Scale(1./Float_t(nEvents));
622 hnF->Scale(1./Float_t(nEvents));
623 hnB->Scale(1./Float_t(nEvents));
624
625 hnFnB_MC->Scale(1./Float_t(nMCEvents));
626 hnF2_MC->Scale(1./Float_t(nMCEvents));
627 hnB2_MC->Scale(1./Float_t(nMCEvents));
628 hnF_MC->Scale(1./Float_t(nMCEvents));
629 hnB_MC->Scale(1./Float_t(nMCEvents));
630
631 for (Int_t bin = 1; bin <= hnFnB->GetNbinsX(); bin++) {
632
633 Double_t nFnBav = hnFnB->GetBinContent(bin);
634 Double_t nF2av = hnF2->GetBinContent(bin);
635 Double_t nB2av = hnB2->GetBinContent(bin);
636 Double_t nFav = hnF->GetBinContent(bin);
637 Double_t nBav = hnB->GetBinContent(bin);
638
639 Double_t nFnB_MCav = hnFnB_MC->GetBinContent(bin);
640 Double_t nF2_MCav = hnF2_MC->GetBinContent(bin);
641 Double_t nB2_MCav = hnB2_MC->GetBinContent(bin);
642 Double_t nF_MCav = hnF_MC->GetBinContent(bin);
643 Double_t nB_MCav = hnB_MC->GetBinContent(bin);
644
645 Double_t bF = ((nF2av-nFav*nFav) == 0 ? 0. : (nFnBav-nFav*nBav)/(nF2av-nFav*nFav));
646 Double_t bB = ((nF2av-nFav*nFav) == 0 ? 0. : (nFnBav-nFav*nBav)/(nB2av-nBav*nBav));
647
648 Double_t bF_MC = ((nF2_MCav-nF_MCav*nF_MCav) == 0 ? 0. :
649 (nFnB_MCav-nF_MCav*nB_MCav)/(nF2_MCav-nF_MCav*nF_MCav));
650 Double_t bB_MC = ((nF2_MCav-nF_MCav*nF_MCav) == 0 ? 0. :
651 (nFnB_MCav-nF_MCav*nB_MCav)/(nB2_MCav-nB_MCav*nB_MCav));
652
653 TH1D *hBFFcor = (TH1D*)fOutputList->FindObject(Form("hBFFcor_binning%d", i));
654 TH1D *hBFBcor = (TH1D*)fOutputList->FindObject(Form("hBFBcor_binning%d", i));
655 TH1D *hBFFcor_MC = (TH1D*)fOutputList->FindObject(Form("hBFFcor_MC_binning%d", i));
656 TH1D *hBFBcor_MC = (TH1D*)fOutputList->FindObject(Form("hBFBcor_MC_binning%d", i));
657
658 hBFFcor->SetBinContent(bin, bF);
659 hBFBcor->SetBinContent(bin, bB);
660 hBFFcor_MC->SetBinContent(bin, bF_MC);
661 hBFBcor_MC->SetBinContent(bin, bB_MC);
662 }
663 }*/
664}
665
1c038f67 666//_____________________________________________________________________
667void AliFMDAnalysisTaskBFCorrelation::ProcessPrimary() {
4b8bdb60 668
1b418b63 669 AliMCEventHandler* eventHandler =
670 dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()
671 ->GetMCtruthEventHandler());
672 if (!eventHandler) return;
673
1c038f67 674 AliMCEvent* mcEvent = eventHandler->MCEvent();
1b418b63 675 if(!mcEvent) return;
1c038f67 676
677 AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
678
679 AliMCParticle* particle = 0;
680 AliStack* stack = mcEvent->Stack();
681
2869f0ae 682 //TH1F* hPrimary = (TH1F*)fOutputList->FindObject("hMultvsEta");
1c038f67 683 AliHeader* header = mcEvent->Header();
684 AliGenEventHeader* genHeader = header->GenEventHeader();
2869f0ae 685 /*
4b8bdb60 686 AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
687
688 if (!pythiaGenHeader) {
689 std::cout<<" no pythia header!"<<std::endl;
690 return;
691 }
2869f0ae 692
4b8bdb60 693 Int_t pythiaType = pythiaGenHeader->ProcessType();
694
695 if(pythiaType==92||pythiaType==93){
696 std::cout<<"single diffractive"<<std::endl;
697 return;
698 }
699 if(pythiaType==94){
700 std::cout<<"double diffractive"<<std::endl;
701 return;
702 }
2869f0ae 703 */
1c038f67 704 TArrayF vertex;
705 genHeader->PrimaryVertex(vertex);
706 if(TMath::Abs(vertex.At(2)) > pars->GetVtxCutZ())
707 return;
2869f0ae 708 // Double_t delta = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
709 // Double_t vertexBinDouble = (vertex.At(2) + pars->GetVtxCutZ()) / delta;
710 // Int_t vertexBin = (Int_t)vertexBinDouble;
711
712 //if (vertexBin != 5) return;
713
4b8bdb60 714 //Bool_t firstTrack = kTRUE;
1c038f67 715
716 // we loop over the primaries only unless we need the hits (diagnostics running slowly)
717 Int_t nTracks = stack->GetNprimary();
4b8bdb60 718 // if(pars->GetProcessHits())
719 // nTracks = stack->GetNtrack();
1c038f67 720
2869f0ae 721 TH2F *hEtaPhiParticleMap = (TH2F*)fInternalList->FindObject("hEtaPhiParticleMap");
722 hEtaPhiParticleMap->Reset();
723
724 for(Int_t i= 0; i < nTracks; i++) {
725 particle = (AliMCParticle*) mcEvent->GetTrack(i);
726 if(!particle) continue;
4b8bdb60 727
2869f0ae 728 if((stack->IsPhysicalPrimary(i)) && (particle->Charge() != 0)) {
729 hEtaPhiParticleMap->Fill(particle->Eta(), particle->Phi());
730 //std::cout << "Hans er nederen" << std::endl;
1c038f67 731 }
732 }
2869f0ae 733 ProjectAndMirror("MC");
4b8bdb60 734 CalculateValues("MC");
1c038f67 735}
2869f0ae 736
1c038f67 737//_____________________________________________________________________
738//
739//
740// EOF