]> 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
Fixed warnings [-Wunused-but-set-variable] from GCC 4.6 -
[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())));
ca47c865 402 if(!hMult) {
403 AliWarning("no hist - returning");
404 return;
405 }
2869f0ae 406 hMult->Reset();
407
408 // Create generic names for retrieving histograms
4b8bdb60 409
2869f0ae 410 TList *list = 0; // List for getting either ESD or MC "hit map"
4b8bdb60 411
2869f0ae 412 TString sEtaPhiMap;
4b8bdb60 413
2869f0ae 414 if (sType.Contains("ESD")) {
4b8bdb60 415
2869f0ae 416 list = fInputList;
4b8bdb60 417
2869f0ae 418 sEtaPhiMap = "dNdetadphiHistogramSPDTrVtx";
419 }
420
421 if (sType.Contains("MC")) {
422
423 list = fInternalList;
424
425 sEtaPhiMap = "hEtaPhiParticleMap";
4b8bdb60 426 }
427
2869f0ae 428 TString sMult("hSERebinMult");
429 TString sMultMirror("hSERebinMultMirror");
430 TString sMultW("hSERebinMultW");
431 TString sMultMirrorW("hSERebinMultMirrorW");
432
433 sType += "_binning";
434
435 sMult += sType;
436 sMultMirror += sType;
437 sMultW += sType;
438 sMultMirrorW += sType;
1c038f67 439
2869f0ae 440 // Get the 2D histogram containing the particles of the event being analyzed
441
442 TH2F *hEtaPhiMap = dynamic_cast<TH2F*>(list->FindObject(sEtaPhiMap));
443
444 TH1D *hProjection = hEtaPhiMap->ProjectionX("hTemporary");
445 hMult->Add(hProjection);
446
447 // Loop over the 3 binnings
448
449 for (Int_t i = 1; i<=4; i++) {
450 if (i == 3) continue;
451
452 // Project the 2D "hit map" onto the eta-axis and rebin
453
454 TH1D *hProjRebin = hEtaPhiMap->ProjectionX("hProjRebin");
455 hProjRebin->Rebin(i*5);
456
457 // Retrieve the histograms to store the Singe Event information and reset
458
1b418b63 459 TH1D *hSEMult = static_cast<TH1D*>(fInternalList->FindObject(Form("%s%d", sMult.Data(), i)));
2869f0ae 460 hSEMult->Reset();
461
1b418b63 462 TH1D *hSEMultMirror = static_cast<TH1D*>(fInternalList->FindObject(Form("%s%d", sMultMirror.Data(), i)));
2869f0ae 463 hSEMultMirror->Reset();
464
1b418b63 465 TH1D *hSEMultW = static_cast<TH1D*>(fInternalList->FindObject(Form("%s%d", sMultW.Data(), i)));
2869f0ae 466 hSEMultW->Reset();
467
1b418b63 468 TH1D *hSEMultMirrorW = static_cast<TH1D*>(fInternalList->FindObject(Form("%s%d", sMultMirrorW.Data(), i)));
2869f0ae 469 hSEMultMirrorW->Reset();
470
471 // Fill the histograms with the Single Event information
472
473 hSEMult->Add(hProjRebin);
474
475 for (Int_t bin = 1; bin <= hSEMult->GetNbinsX(); bin++) {
476
477 hSEMultMirror->SetBinContent(hSEMultMirror->FindBin(-hProjRebin->GetBinCenter(bin)), hProjRebin->GetBinContent(bin));
478 hSEMultMirror->SetBinError(hSEMultMirror->FindBin(-hProjRebin->GetBinCenter(bin)), hProjRebin->GetBinError(bin));
479
480 // hMultDist->Fill(bin, hMultTemp->GetBinContent(bin));
481 }
482 hProjRebin->Delete();
483 }
1c038f67 484}
4b8bdb60 485
2869f0ae 486//_____________________________________________________________________
487void AliFMDAnalysisTaskBFCorrelation::CalculateValues(TString sType) {
488
489 sType.ToUpper();
490
491 if (!sType.Contains("ESD") && !sType.Contains("MC")) {
492 std::cout << "Wrong type specification for 'CalculateValues'" << std::endl;
493 return;
494 }
495
496 TString sMult("hSERebinMult");
497 TString sMultMirror("hSERebinMultMirror");
498 // TString_t *sMultW("hSEMultW");
499 // TString_t *sMultMirrorW("hSEMultMirrorW");
500
501 TString sn("hRebinn");
502 TString snMirror("hRebinnMirror");
503 TString sn2("hRebinn2");
504 TString sn2Mirror("hRebinn2Mirror");
505 TString snfb("hRebinnfb");
4b8bdb60 506
2869f0ae 507 sType += "_binning";
4b8bdb60 508
2869f0ae 509 sMult += sType;
510 sMultMirror += sType;
4b8bdb60 511
2869f0ae 512 sn += sType;
513 snMirror += sType;
514 sn2 += sType;
515 sn2Mirror += sType;
516 snfb += sType;
517
4b8bdb60 518 for (Int_t i = 1; i <= 4; i++) {
2869f0ae 519 if (i == 3) continue;
4b8bdb60 520
2869f0ae 521 TH1D *hSEMult = (TH1D*)fInternalList->FindObject(Form("%s%d", sMult.Data(), i));
522 TH1D *hSEMultMirror = (TH1D*)fInternalList->FindObject(Form("%s%d", sMultMirror.Data(), i));
523 /*
524 TH1D *hSEMultW = (TH1D*)fInternalList->FindObject(Form("%s%d", cMultW, cType, i));
525 TH1D *hSEMultMirrorW = (TH1D*)fInternalList->FindObject(Form("%s%d", cMultMirrorW, cType, i));
526 */
4b8bdb60 527
2869f0ae 528 TH1F *hn = (TH1F*)fOutputList->FindObject(Form("%s%d", sn.Data(), i));
529 TH1F *hnMirror = (TH1F*)fOutputList->FindObject(Form("%s%d", snMirror.Data(), i));
530 TH1F *hn2 = (TH1F*)fOutputList->FindObject(Form("%s%d", sn2.Data(), i));
531 TH1F *hn2Mirror = (TH1F*)fOutputList->FindObject(Form("%s%d", sn2Mirror.Data(), i));
532 TH1F *hnfb = (TH1F*)fOutputList->FindObject(Form("%s%d", snfb.Data(), i));
533 /*
534 TH1F *hnW = (TH1F*)fOutputList->FindObject(Form("%sW%s_binning%d", cn, cType, i));
535 TH1F *hnMirrorW = (TH1F*)fOutputList->FindObject(Form("%sW%s_binning%d", cnMirror, cType, i));
536 TH1F *hn2W = (TH1F*)fOutputList->FindObject(Form("%sW%s_binning%d", cn2, cType, i));
537 TH1F *hn2MirrorW = (TH1F*)fOutputList->FindObject(Form("%sW%s_binning%d", cn2Mirror, cType, i));
538 TH1F *hnfbW = (TH1F*)fOutputList->FindObject(Form("%sW%s_binning%d", cnfb, cType, i));
539 */
540 TH1D *hTemp = (TH1D*)fInternalList->FindObject(Form("hRebinTemp_binning%d", i));
541
542 hn->Add(hSEMult);
543
544 hnMirror->Add(hSEMultMirror);
545
546 hTemp->Reset();
547 hTemp->Add(hSEMult);
548 hTemp->Multiply(hSEMult);
549 hn2->Add(hTemp);
550
551 hTemp->Reset();
552 hTemp->Add(hSEMultMirror);
553 hTemp->Multiply(hSEMultMirror);
554 hn2Mirror->Add(hTemp);
555
556 hSEMultMirror->Multiply(hSEMult);
557 hnfb->Add(hSEMultMirror);
4b8bdb60 558 }
559}
560
561//_____________________________________________________________________
2869f0ae 562void AliFMDAnalysisTaskBFCorrelation::MultiplicityVsEta(TString sType) {
563
564 sType.ToUpper();
4b8bdb60 565
2869f0ae 566 if (!sType.Contains("ESD") && !sType.Contains("MC")) {
567 std::cout << "Wrong type specification for 'MultiplicityVsEta'" << std::endl;
568 return;
569 }
570}
1c038f67 571
4b8bdb60 572
573//_____________________________________________________________________
2869f0ae 574void AliFMDAnalysisTaskBFCorrelation::CreateResponseMatrix() {
575
576 TH2F *hResponseMatrix = (TH2F*)fOutputList->FindObject("hResponseMatrix");
577
578 TH1D *hSEMultESD = (TH1D*)fInternalList->FindObject("hSEMultESD");
579 TH1D *hSEMultMC = (TH1D*)fInternalList->FindObject("hSEMultMC");
1c038f67 580
2869f0ae 581 Float_t HitsESD = 0;
582 Float_t HitsMC = 0;
583
584 for (Int_t bin = 1; bin<= fnBinsX; bin++) {
585
586 if ((hSEMultMC->GetBinLowEdge(bin) > -3.4) &&
587 (hSEMultMC->GetBinLowEdge(bin+1) < 5.)) {
588
589 HitsESD += hSEMultESD->GetBinContent(bin);
590 HitsMC += hSEMultMC->GetBinContent(bin);
591 }
592 }
593
594 hResponseMatrix->Fill(HitsMC, HitsESD);
1c038f67 595}
2869f0ae 596
597//_____________________________________________________________________
598void AliFMDAnalysisTaskBFCorrelation::Terminate(Option_t */*option*/) {
599 /*
600 std::cout << "Terminating !" << std::endl;
601
602 TH1I *hnEvents = (TH1I*)fOutputList->FindObject("nEvents");
603 TH1I *hnMCEvents = (TH1I*)fOutputList->FindObject("nEvents");
604
605 Int_t nEvents = hnEvents->GetEntries();
606 Int_t nMCEvents = hnMCEvents->GetEntries();
607
608 for (Int_t i = 1; i <= 4; i++) {
609 if (i == 3) continue;
610
611 TH1F *hnFnB = (TH1F*)fOutputList->FindObject(Form("hnFnB_binning%d", i));
612 TH1F *hnF2 = (TH1F*)fOutputList->FindObject(Form("hnF2_binning%d", i));
613 TH1F *hnB2 = (TH1F*)fOutputList->FindObject(Form("hnB2_binning%d", i));
614 TH1F *hnF = (TH1F*)fOutputList->FindObject(Form("hnF_binning%d", i));
615 TH1F *hnB = (TH1F*)fOutputList->FindObject(Form("hnB_binning%d", i));
616
617 TH1F *hnFnB_MC = (TH1F*)fOutputList->FindObject(Form("hnFnB_MC_binning%d", i));
618 TH1F *hnF2_MC = (TH1F*)fOutputList->FindObject(Form("hnF2_MC_binning%d", i));
619 TH1F *hnB2_MC = (TH1F*)fOutputList->FindObject(Form("hnB2_MC_binning%d", i));
620 TH1F *hnF_MC = (TH1F*)fOutputList->FindObject(Form("hnF_MC_binning%d", i));
621 TH1F *hnB_MC = (TH1F*)fOutputList->FindObject(Form("hnB_MC_binning%d", i));
622
623 hnFnB->Scale(1./Float_t(nEvents));
624 hnF2->Scale(1./Float_t(nEvents));
625 hnB2->Scale(1./Float_t(nEvents));
626 hnF->Scale(1./Float_t(nEvents));
627 hnB->Scale(1./Float_t(nEvents));
628
629 hnFnB_MC->Scale(1./Float_t(nMCEvents));
630 hnF2_MC->Scale(1./Float_t(nMCEvents));
631 hnB2_MC->Scale(1./Float_t(nMCEvents));
632 hnF_MC->Scale(1./Float_t(nMCEvents));
633 hnB_MC->Scale(1./Float_t(nMCEvents));
634
635 for (Int_t bin = 1; bin <= hnFnB->GetNbinsX(); bin++) {
636
637 Double_t nFnBav = hnFnB->GetBinContent(bin);
638 Double_t nF2av = hnF2->GetBinContent(bin);
639 Double_t nB2av = hnB2->GetBinContent(bin);
640 Double_t nFav = hnF->GetBinContent(bin);
641 Double_t nBav = hnB->GetBinContent(bin);
642
643 Double_t nFnB_MCav = hnFnB_MC->GetBinContent(bin);
644 Double_t nF2_MCav = hnF2_MC->GetBinContent(bin);
645 Double_t nB2_MCav = hnB2_MC->GetBinContent(bin);
646 Double_t nF_MCav = hnF_MC->GetBinContent(bin);
647 Double_t nB_MCav = hnB_MC->GetBinContent(bin);
648
649 Double_t bF = ((nF2av-nFav*nFav) == 0 ? 0. : (nFnBav-nFav*nBav)/(nF2av-nFav*nFav));
650 Double_t bB = ((nF2av-nFav*nFav) == 0 ? 0. : (nFnBav-nFav*nBav)/(nB2av-nBav*nBav));
651
652 Double_t bF_MC = ((nF2_MCav-nF_MCav*nF_MCav) == 0 ? 0. :
653 (nFnB_MCav-nF_MCav*nB_MCav)/(nF2_MCav-nF_MCav*nF_MCav));
654 Double_t bB_MC = ((nF2_MCav-nF_MCav*nF_MCav) == 0 ? 0. :
655 (nFnB_MCav-nF_MCav*nB_MCav)/(nB2_MCav-nB_MCav*nB_MCav));
656
657 TH1D *hBFFcor = (TH1D*)fOutputList->FindObject(Form("hBFFcor_binning%d", i));
658 TH1D *hBFBcor = (TH1D*)fOutputList->FindObject(Form("hBFBcor_binning%d", i));
659 TH1D *hBFFcor_MC = (TH1D*)fOutputList->FindObject(Form("hBFFcor_MC_binning%d", i));
660 TH1D *hBFBcor_MC = (TH1D*)fOutputList->FindObject(Form("hBFBcor_MC_binning%d", i));
661
662 hBFFcor->SetBinContent(bin, bF);
663 hBFBcor->SetBinContent(bin, bB);
664 hBFFcor_MC->SetBinContent(bin, bF_MC);
665 hBFBcor_MC->SetBinContent(bin, bB_MC);
666 }
667 }*/
668}
669
1c038f67 670//_____________________________________________________________________
671void AliFMDAnalysisTaskBFCorrelation::ProcessPrimary() {
4b8bdb60 672
1b418b63 673 AliMCEventHandler* eventHandler =
674 dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()
675 ->GetMCtruthEventHandler());
676 if (!eventHandler) return;
677
1c038f67 678 AliMCEvent* mcEvent = eventHandler->MCEvent();
1b418b63 679 if(!mcEvent) return;
1c038f67 680
681 AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
682
683 AliMCParticle* particle = 0;
684 AliStack* stack = mcEvent->Stack();
685
2869f0ae 686 //TH1F* hPrimary = (TH1F*)fOutputList->FindObject("hMultvsEta");
1c038f67 687 AliHeader* header = mcEvent->Header();
688 AliGenEventHeader* genHeader = header->GenEventHeader();
2869f0ae 689 /*
4b8bdb60 690 AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
691
692 if (!pythiaGenHeader) {
693 std::cout<<" no pythia header!"<<std::endl;
694 return;
695 }
2869f0ae 696
4b8bdb60 697 Int_t pythiaType = pythiaGenHeader->ProcessType();
698
699 if(pythiaType==92||pythiaType==93){
700 std::cout<<"single diffractive"<<std::endl;
701 return;
702 }
703 if(pythiaType==94){
704 std::cout<<"double diffractive"<<std::endl;
705 return;
706 }
2869f0ae 707 */
1c038f67 708 TArrayF vertex;
709 genHeader->PrimaryVertex(vertex);
710 if(TMath::Abs(vertex.At(2)) > pars->GetVtxCutZ())
711 return;
2869f0ae 712 // Double_t delta = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
713 // Double_t vertexBinDouble = (vertex.At(2) + pars->GetVtxCutZ()) / delta;
714 // Int_t vertexBin = (Int_t)vertexBinDouble;
715
716 //if (vertexBin != 5) return;
717
4b8bdb60 718 //Bool_t firstTrack = kTRUE;
1c038f67 719
720 // we loop over the primaries only unless we need the hits (diagnostics running slowly)
721 Int_t nTracks = stack->GetNprimary();
4b8bdb60 722 // if(pars->GetProcessHits())
723 // nTracks = stack->GetNtrack();
1c038f67 724
2869f0ae 725 TH2F *hEtaPhiParticleMap = (TH2F*)fInternalList->FindObject("hEtaPhiParticleMap");
726 hEtaPhiParticleMap->Reset();
727
728 for(Int_t i= 0; i < nTracks; i++) {
729 particle = (AliMCParticle*) mcEvent->GetTrack(i);
730 if(!particle) continue;
4b8bdb60 731
2869f0ae 732 if((stack->IsPhysicalPrimary(i)) && (particle->Charge() != 0)) {
733 hEtaPhiParticleMap->Fill(particle->Eta(), particle->Phi());
734 //std::cout << "Hans er nederen" << std::endl;
1c038f67 735 }
736 }
2869f0ae 737 ProjectAndMirror("MC");
4b8bdb60 738 CalculateValues("MC");
1c038f67 739}
2869f0ae 740
1c038f67 741//_____________________________________________________________________
742//
743//
744// EOF