CbmRoot
Loading...
Searching...
No Matches
CbmHadronAnalysis.cxx
Go to the documentation of this file.
1/* Copyright (C) 2012-2021 PI-UHd, GSI
2 SPDX-License-Identifier: GPL-3.0-only
3 Authors: Norbert Herrmann [committer], Florian Uhlig */
4
5// ------------------------------------------------------------------
6// ----- CbmHadronAnalysis -----
7// ----- Created 14/12/2012 by nh -----
8// ------------------------------------------------------------------
9#include <iostream>
10using namespace std;
11
12#include "CbmCluster.h"
13#include "CbmDigiManager.h"
14#include "CbmEvent.h"
15#include "CbmGlobalTrack.h"
16#include "CbmHadron.h"
17#include "CbmHadronAnalysis.h"
18#include "CbmLink.h"
19#include "CbmMCDataArray.h"
20#include "CbmMCDataManager.h"
21#include "CbmMCTrack.h"
22#include "CbmMatch.h"
23#include "CbmStsCluster.h"
24#include "CbmStsDigi.h"
25#include "CbmStsHit.h"
26#include "CbmStsKFTrackFitter.h"
27#include "CbmStsPoint.h"
28#include "CbmStsTrack.h"
29#include "CbmTofCell.h"
31#include "CbmTofGeoHandler.h"
32#include "CbmTofHit.h"
33#include "CbmTofPoint.h"
34#include "CbmTofTrack.h"
35#include "CbmTrdHit.h"
36#include "CbmTrdPoint.h"
37#include "CbmVertex.h"
38
39#include <FairMCEventHeader.h>
40#include <FairMCPoint.h>
41#include <FairRootFileSink.h>
42#include <FairRootManager.h>
43#include <FairRunAna.h>
44#include <Logger.h>
45
46#include <TClonesArray.h>
47#include <TFile.h>
48#include <TH1.h>
49#include <TH1F.h>
50#include <TH2F.h>
51#include <TH3.h>
52#include <TH3F.h>
53#include <TMath.h>
54#include <TROOT.h>
55#include <TRandom.h>
56#include <TString.h>
57
58CbmDigiManager* fDigiMan; // TOF Input Digis
59TClonesArray* fEventsColl; // CBMEvents (time based)
60TClonesArray* fTofTracks; // CbmTofTrack array
61TClonesArray* fTofHitsColl;
62TClonesArray* fStsHitsColl;
63
64static Int_t iNbTs = 0;
65
66#define M2PI 0.019479835
67#define M2KA 0.24371698
68#define M2PROT 0.88035435
69#define clight 29.9792
70#define MinWallDist 550.
71
72Int_t nMCTracks = 0, nTofPoints = 0, nTofHits = 0;
73Int_t nTofTracks = 0;
74Int_t nStsHits = 0;
75Int_t nGlobTracks = 0;
76Int_t NMASS = 3;
77Float_t refMass[3] = {0.139, 0.494, 0.938};
78
79static TH1F* fhTofHitMul;
80static TH1F* fhStsHitMul;
81const Int_t fiTofHitMulMax = 200;
82static TFile* fHist;
83
84//___________________________________________________________________
85//
86// CbmHadronAnalysis
87//
88// Task for analysis of hadron spectra
89//
90// ------------------------------------------------------------------
92{
93 cout << "CbmHadronAnalysis: Task started with defaults" << endl;
94}
95// ------------------------------------------------------------------
96
97// ------------------------------------------------------------------
98CbmHadronAnalysis::CbmHadronAnalysis(const char* name, Int_t verbose)
99 : FairTask(name, verbose)
100 , fEvents(0)
101 , fBeamMomentum(10)
102 , fMidY(0.)
103 , fDY(0.)
104 , fFlowMinPtm(0.)
105 , fBSelMin(0.)
106 , fBSelMax(0.)
107 , fwxy2(0.)
108 , fWMax(0.)
109 , fVtxBMax(0.)
110 , fPdfFileName("")
111 , fFlowFileName("")
112 , fflowFile(NULL)
113 , fMCEventHeader(NULL)
114 , fGeoHandler(NULL)
115 , fCellInfo(NULL)
116 , fMCTracks(NULL)
117 , fStsPoints(NULL)
118 , fMCTracksColl(NULL)
119 , fStsPointsColl(NULL)
120 , fStsHits(NULL)
121 , fStsClusters(NULL)
122 , fStsTracks(NULL)
123 , fStsDigis(NULL)
124 , fStsDigiMatchColl(NULL)
125 , fTrdPoints(NULL)
126 , fTrdHits(NULL)
127 , fTofPoints(NULL)
128 , fTofHits(NULL)
129 , fTofDigis(NULL)
130 , fTofDigiMatchColl(NULL)
131 , fTofDigiMatchPointsColl(NULL)
132 , fGlobalTracks(NULL)
133 , fHadrons(NULL)
134 , fPrimVertex(NULL)
135 , fTrackFitter()
136 , fa_ptm_rap_gen_pip(NULL)
137 , fa_ptm_rap_gen_pim(NULL)
138 , fa_ptm_rap_gen_kp(NULL)
139 , fa_ptm_rap_gen_km(NULL)
140 , fa_ptm_rap_gen_p(NULL)
141 , fa_ptm_rap_gen_pbar(NULL)
142 , fa_ptm_rap_gen_d(NULL)
143 , fa_ptm_rap_gen_t(NULL)
144 , fa_ptm_rap_gen_h(NULL)
145 , fa_ptm_rap_gen_a(NULL)
146 , fa_ptm_rap_gen_imf(NULL)
147 , fa_plab_sts_pip(NULL)
148 , fa_plab_sts_pim(NULL)
149 , fa_plab_sts_kp(NULL)
150 , fa_plab_sts_km(NULL)
151 , fa_plab_sts_p(NULL)
152 , fa_plab_sts_pbar(NULL)
153 , fa_ptm_rap_sts_pip(NULL)
154 , fa_ptm_rap_sts_pim(NULL)
155 , fa_ptm_rap_sts_kp(NULL)
156 , fa_ptm_rap_sts_km(NULL)
157 , fa_ptm_rap_sts_p(NULL)
158 , fa_ptm_rap_sts_pbar(NULL)
159 , fa_ptm_rap_sts_d(NULL)
160 , fa_ptm_rap_sts_t(NULL)
161 , fa_ptm_rap_sts_h(NULL)
162 , fa_ptm_rap_sts_a(NULL)
163 , fa_ptm_rap_sts_imf(NULL)
164 , fa_ptm_rap_poi_pip(NULL)
165 , fa_ptm_rap_poi_pim(NULL)
166 , fa_ptm_rap_poi_kp(NULL)
167 , fa_ptm_rap_poi_km(NULL)
168 , fa_ptm_rap_poi_p(NULL)
169 , fa_ptm_rap_poi_pbar(NULL)
170 , fa_ptm_rap_poi_d(NULL)
171 , fa_ptm_rap_poi_t(NULL)
172 , fa_ptm_rap_poi_h(NULL)
173 , fa_ptm_rap_poi_a(NULL)
174 , fa_ptm_rap_poi_imf(NULL)
175 , fa_ptm_rap_hit_pip(NULL)
176 , fa_ptm_rap_hit_pim(NULL)
177 , fa_ptm_rap_hit_kp(NULL)
178 , fa_ptm_rap_hit_km(NULL)
179 , fa_ptm_rap_hit_p(NULL)
180 , fa_ptm_rap_hit_pbar(NULL)
181 , fa_ptm_rap_hit_d(NULL)
182 , fa_ptm_rap_hit_t(NULL)
183 , fa_ptm_rap_hit_h(NULL)
184 , fa_ptm_rap_hit_a(NULL)
185 , fa_ptm_rap_hit_imf(NULL)
186 , fa_ptm_rap_glo_pip(NULL)
187 , fa_ptm_rap_glo_pim(NULL)
188 , fa_ptm_rap_glo_kp(NULL)
189 , fa_ptm_rap_glo_km(NULL)
190 , fa_ptm_rap_glo_p(NULL)
191 , fa_ptm_rap_glo_pbar(NULL)
192 , fa_ptm_rap_glo_d(NULL)
193 , fa_ptm_rap_glo_t(NULL)
194 , fa_ptm_rap_glo_h(NULL)
195 , fa_ptm_rap_glo_a(NULL)
196 , fa_ptm_rap_glo_imf(NULL)
197 , fa_mul_b_gen(NULL)
198 , fa_mul_b_poi(NULL)
199 , fa_mul_b_hit(NULL)
200 , fa_mul_b_glo(NULL)
201 , fa_mul_b_had(NULL)
202 , fa_phirp_b_gen(NULL)
203 , fa_phgrp_b_gen(NULL)
204 , fa_phphrp_gen(NULL)
205 , fa_delrp_b_gen(NULL)
206 , fa_delrp_b_poi(NULL)
207 , fa_delrp_b_hit(NULL)
208 , fa_delrp_b_glo(NULL)
209 , fa_drp_b_gen(NULL)
210 , fa_cdrp_b_gen(NULL)
211 , fa_drp_b_poi(NULL)
212 , fa_cdrp_b_poi(NULL)
213 , fa_drp_b_hit(NULL)
214 , fa_cdrp_b_hit(NULL)
215 , fa_drp_b_glo(NULL)
216 , fa_cdrp_b_glo(NULL)
217 , fa_drp_b_had(NULL)
218 , fa_cdrp_b_had(NULL)
219 , fa_cdelrp_b_gen(NULL)
220 , fa_cdelrp_b_poi(NULL)
221 , fa_cdelrp_b_hit(NULL)
222 , fa_cdelrp_b_glo(NULL)
223 , fa_cdelrp_b_had(NULL)
224 , fa_phirp_gen(NULL)
225 , fa_phirp_poi(NULL)
226 , fa_phirp_hit(NULL)
227 , fa_phirp_glo(NULL)
228 , fa_phirp_had(NULL)
229 , fa_phirps_gen(NULL)
230 , fa_phirps_poi(NULL)
231 , fa_phirps_hit(NULL)
232 , fa_phirps_glo(NULL)
233 , fa_phirps_had(NULL)
234 , fa_v1_rap_gen_pip(NULL)
235 , fa_v1_rap_gen_pim(NULL)
236 , fa_v1_rap_gen_kp(NULL)
237 , fa_v1_rap_gen_km(NULL)
238 , fa_v1_rap_gen_p(NULL)
239 , fa_v1_rap_gen_pbar(NULL)
240 , fa_v1_rap_gen_d(NULL)
241 , fa_v1_rap_gen_t(NULL)
242 , fa_v1_rap_gen_h(NULL)
243 , fa_v1_rap_gen_a(NULL)
244 , fa_v1_rap_gen_imf(NULL)
245 , fa_v2_rap_gen_pip(NULL)
246 , fa_v2_rap_gen_pim(NULL)
247 , fa_v2_rap_gen_kp(NULL)
248 , fa_v2_rap_gen_km(NULL)
249 , fa_v2_rap_gen_p(NULL)
250 , fa_v2_rap_gen_pbar(NULL)
251 , fa_v2_rap_gen_d(NULL)
252 , fa_v2_rap_gen_t(NULL)
253 , fa_v2_rap_gen_h(NULL)
254 , fa_v2_rap_gen_a(NULL)
255 , fa_v2_rap_gen_imf(NULL)
256 , fa_v1_rap_poi_pip(NULL)
257 , fa_v1_rap_poi_pim(NULL)
258 , fa_v1_rap_poi_kp(NULL)
259 , fa_v1_rap_poi_km(NULL)
260 , fa_v1_rap_poi_p(NULL)
261 , fa_v1_rap_poi_pbar(NULL)
262 , fa_v1_rap_poi_d(NULL)
263 , fa_v1_rap_poi_t(NULL)
264 , fa_v1_rap_poi_h(NULL)
265 , fa_v1_rap_poi_a(NULL)
266 , fa_v1_rap_poi_imf(NULL)
267 , fa_v2_rap_poi_pip(NULL)
268 , fa_v2_rap_poi_pim(NULL)
269 , fa_v2_rap_poi_kp(NULL)
270 , fa_v2_rap_poi_km(NULL)
271 , fa_v2_rap_poi_p(NULL)
272 , fa_v2_rap_poi_pbar(NULL)
273 , fa_v2_rap_poi_d(NULL)
274 , fa_v2_rap_poi_t(NULL)
275 , fa_v2_rap_poi_h(NULL)
276 , fa_v2_rap_poi_a(NULL)
277 , fa_v2_rap_poi_imf(NULL)
278 , fa_v1_rap_hit_pip(NULL)
279 , fa_v1_rap_hit_pim(NULL)
280 , fa_v1_rap_hit_kp(NULL)
281 , fa_v1_rap_hit_km(NULL)
282 , fa_v1_rap_hit_p(NULL)
283 , fa_v1_rap_hit_pbar(NULL)
284 , fa_v1_rap_hit_d(NULL)
285 , fa_v1_rap_hit_t(NULL)
286 , fa_v1_rap_hit_h(NULL)
287 , fa_v1_rap_hit_a(NULL)
288 , fa_v1_rap_hit_imf(NULL)
289 , fa_v2_rap_hit_pip(NULL)
290 , fa_v2_rap_hit_pim(NULL)
291 , fa_v2_rap_hit_kp(NULL)
292 , fa_v2_rap_hit_km(NULL)
293 , fa_v2_rap_hit_p(NULL)
294 , fa_v2_rap_hit_pbar(NULL)
295 , fa_v2_rap_hit_d(NULL)
296 , fa_v2_rap_hit_t(NULL)
297 , fa_v2_rap_hit_h(NULL)
298 , fa_v2_rap_hit_a(NULL)
299 , fa_v2_rap_hit_imf(NULL)
300 , fa_v1_rap_glo_pip(NULL)
301 , fa_v1_rap_glo_pim(NULL)
302 , fa_v1_rap_glo_kp(NULL)
303 , fa_v1_rap_glo_km(NULL)
304 , fa_v1_rap_glo_p(NULL)
305 , fa_v1_rap_glo_pbar(NULL)
306 , fa_v1_rap_glo_d(NULL)
307 , fa_v1_rap_glo_t(NULL)
308 , fa_v1_rap_glo_h(NULL)
309 , fa_v1_rap_glo_a(NULL)
310 , fa_v1_rap_glo_imf(NULL)
311 , fa_v2_rap_glo_pip(NULL)
312 , fa_v2_rap_glo_pim(NULL)
313 , fa_v2_rap_glo_kp(NULL)
314 , fa_v2_rap_glo_km(NULL)
315 , fa_v2_rap_glo_p(NULL)
316 , fa_v2_rap_glo_pbar(NULL)
317 , fa_v2_rap_glo_d(NULL)
318 , fa_v2_rap_glo_t(NULL)
319 , fa_v2_rap_glo_h(NULL)
320 , fa_v2_rap_glo_a(NULL)
321 , fa_v2_rap_glo_imf(NULL)
322 , fa_xy_poi1(NULL)
323 , fa_xy_poi2(NULL)
324 , fa_xy_poi3(NULL)
325 , fa_xy_hit1(NULL)
326 , fa_xy_hit2(NULL)
327 , fa_xy_hit3(NULL)
328 , fa_xy_glo1(NULL)
329 , fa_xy_glo_pip(NULL)
330 , fa_xy_glo_pim(NULL)
331 , fa_xy_glo_kp(NULL)
332 , fa_xy_glo_km(NULL)
333 , fa_xy_glo_p(NULL)
334 , fa_xy_glo_pbar(NULL)
335 , fa_xy_glo_d(NULL)
336 , fa_xy_glo_t(NULL)
337 , fa_xy_glo_h(NULL)
338 , fa_xy_glo_a(NULL)
339 , fa_pv_poi(NULL)
340 , fa_tm_poi(NULL)
341 , fa_tm_poiprim(NULL)
342 , fa_dxx(NULL)
343 , fa_dxy(NULL)
344 , fa_dxz(NULL)
345 , fa_dyx(NULL)
346 , fa_dyy(NULL)
347 , fa_dyz(NULL)
348 , fa_dzx(NULL)
349 , fa_dzy(NULL)
350 , fa_dzz(NULL)
351 , fa_hit_ch(NULL)
352 , fa_dhit_ch(NULL)
353 , fa_tof_hit(NULL)
354 , fa_dtof_hit(NULL)
355 , fa_tof_hitprim(NULL)
356 , fa_pv_hit(NULL)
357 , fa_tm_hit(NULL)
358 , fa_tm_hitprim(NULL)
359 , fa_tn_hit(NULL)
360 , fa_t0_hit(NULL)
361 , fa_t0m_hit(NULL)
362 , fa_t0mn_hit(NULL)
363 , fa_t0m_b_hit(NULL)
364 , fa_t0mn_b_hit(NULL)
365 , fa_t0m_f_hit(NULL)
366 , fa_t0mn_f_hit(NULL)
367 , fa_t0m_f_b_hit(NULL)
368 , fa_t0mn_f_b_hit(NULL)
369 , fa_t0m_nf_hit(NULL)
370 , fa_t0mn_nf_hit(NULL)
371 , fa_t0m_nf_b_hit(NULL)
372 , fa_t0mn_nf_b_hit(NULL)
373 , fa_TofTrackMul(NULL)
374 , fa_VtxB(NULL)
375 , fa_chi2_mom_glo(NULL)
376 , fa_chi2_mom_gloprim(NULL)
377 , fa_len_mom_glo(NULL)
378 , fa_pv_glo(NULL)
379 , fa_tm_glo(NULL)
380 , fa_tm_glo_pip(NULL)
381 , fa_tm_glo_pim(NULL)
382 , fa_tm_glo_kp(NULL)
383 , fa_tm_glo_km(NULL)
384 , fa_tm_glo_p(NULL)
385 , fa_tm_glo_pbar(NULL)
386 , fa_tm_glo_d(NULL)
387 , fa_tm_glo_t(NULL)
388 , fa_tm_glo_h(NULL)
389 , fa_tm_glo_a(NULL)
390 , fa_tm_gloprim(NULL)
391 , fa_tm_glomis(NULL)
392 , fa_tm_glovtxb(NULL)
393 , fa_tm_gloprimvtxb(NULL)
394 , fa_m2mom_glo(NULL)
395 , fa_m2mom_glovtxb(NULL)
396 , fa_m2mom_gloprim(NULL)
397 , fa_m2mom_gloprimvtxb(NULL)
398 , fa_m2mom_glo_pip(NULL)
399 , fa_m2mom_glo_pim(NULL)
400 , fa_m2mom_glo_kp(NULL)
401 , fa_m2mom_glo_km(NULL)
402 , fa_m2mom_glo_p(NULL)
403 , fa_m2mom_glo_pbar(NULL)
404 , fa_m2mom_glo_d(NULL)
405 , fa_m2mom_glo_t(NULL)
406 , fa_m2mom_glo_h(NULL)
407 , fa_m2mom_glo_a(NULL)
408 , fa_pMCmom_glo(NULL)
409 , fa_pMCmom_glo_pip(NULL)
410 , fa_pMCmom_glo_pim(NULL)
411 , fa_pMCmom_glo_kp(NULL)
412 , fa_pMCmom_glo_km(NULL)
413 , fa_pMCmom_glo_p(NULL)
414 , fa_pMCmom_glo_pbar(NULL)
415 , fa_pMCmom_glo_d(NULL)
416 , fa_pMCmom_glo_t(NULL)
417 , fa_pMCmom_glo_h(NULL)
418 , fa_pMCmom_glo_a(NULL)
419 , fa_w_mom_glo(NULL)
420 , fa_w_mom_glo_pip(NULL)
421 , fa_w_mom_glo_pim(NULL)
422 , fa_w_mom_glo_kp(NULL)
423 , fa_w_mom_glo_km(NULL)
424 , fa_w_mom_glo_p(NULL)
425 , fa_w_mom_glo_pbar(NULL)
426 , fa_w_mom_glo_d(NULL)
427 , fa_w_mom_glo_t(NULL)
428 , fa_w_mom_glo_h(NULL)
429 , fa_w_mom_glo_a(NULL)
430 , fa_w_mom_gloprim(NULL)
431 , fa_w_mom_glomis(NULL)
432 , fa_LenDismom_glo(NULL)
433 , fa_LenDismom_glo_pip(NULL)
434 , fa_LenDismom_glo_pim(NULL)
435 , fa_LenDismom_glo_kp(NULL)
436 , fa_LenDismom_glo_km(NULL)
437 , fa_LenDismom_glo_p(NULL)
438 , fa_LenDismom_glo_pbar(NULL)
439 , fa_LenDismom_glo_d(NULL)
440 , fa_LenDismom_glo_t(NULL)
441 , fa_LenDismom_glo_h(NULL)
442 , fa_LenDismom_glo_a(NULL)
443 , fa_LenMcLenGlomom_glo(NULL)
444 , fa_LenMcDismom_glo(NULL)
445 , fhwdist(NULL)
446 , fhwmindelmass(NULL)
447 , fhwminlen(NULL)
448 , fhwdelp(NULL)
449 , fhTofTrkDx(NULL)
450 , fhTofTrkDy(NULL)
451 , fhTofTrkDxsel(NULL)
452 , fhTofTrkDysel(NULL)
453 , bRecSec(kFALSE)
454 , fdDistPrimLim(1.5)
455 , // Ext Parameter: Max Tof-Sts trans distance for primaries
456 fdDistPrimLim2(0.3)
457 , // Ext Parameter: Max Sts-Sts trans distance for primaries
458 fdDistSecLim2(0.5)
459 , // Ext Parameter: Max Sts-Sts trans distance from TOF direction for secondaries
460 fdD0ProtLim(0.4)
461 , // Ext Parameter: Min impact parameter for secondary proton
462 fdOpAngMin(0.01)
463 , // Ext Parameter: Min opening angle for accepting pair
464 fdDCALim(0.2)
465 , // Ext Parameter: Max DCA for accepting pair
466 fdVLenMin(5.)
467 , // Ext Parameter: Min Lambda flight path length for accepting pair
468 fdVLenMax(25.)
469 , // Ext Parameter: Max Lambda flight path length for accepting pair
470 fdDistTRD(10.)
471 , // max accepted distance of Trd Hit from STS-TOF line
472 fdTRDHmulMin(0.)
473 , // min associated Trd Hits to Track candidates
474 fNMixedEvents(1)
475{
476}
477// ------------------------------------------------------------------
478
479
480// ------------------------------------------------------------------
482{
483 // Destructor
484 fPrimVertex = NULL;
485}
486// ------------------------------------------------------------------
487
488
489// ------------------------------------------------------------------
491{
492 // Create histogramms
493 // gROOT->cd();
494 auto sink = FairRunAna::Instance()->GetSink();
495 assert(sink->GetSinkType() == kFILESINK);
496 auto rootFileSink = static_cast<FairRootFileSink*>(sink);
497 fHist = rootFileSink->GetRootFile();
498
499 TString hname = fHist->GetName();
500 hname.Insert(hname.Length() - 5, ".HadAna");
501 fHist = new TFile(hname, "recreate");
502 LOG(info) << "CreateHistograms in " << fHist->GetName();
503 // gSystem->cd(fHist->GetName());
504 fHist->ReOpen("Update");
505
506 fhTofHitMul = new TH1F(Form("hTofHitMul"), Form("TofHit Multiplicity; M_{TofHit} "), fiTofHitMulMax, 0.,
507 (Double_t) fiTofHitMulMax);
508 fhStsHitMul = new TH1F(Form("hStsHitMul"), Form("StsHit Multiplicity; M_{StsHit} "), fiTofHitMulMax * 10, 0.,
509 (Double_t) fiTofHitMulMax * 10);
510
511 Float_t ymin = -1.;
512 Float_t ymax = 4.;
513 Float_t ptmmax = 2.5;
514 Int_t ptm_nbx = 30;
515 Int_t ptm_nby = 30;
516
517 // generator level
518
520 new TH2F("ptm_rap_gen_pip", "MCTrack-gen pi-plus; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
522 new TH2F("ptm_rap_gen_pim", "MCTrack-gen pi-minus;y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
524 new TH2F("ptm_rap_gen_kp", "MCTrack-gen k-plus; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
526 new TH2F("ptm_rap_gen_km", "MCTrack-gen k-minus; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
528 new TH2F("ptm_rap_gen_p", "MCTrack-gen proton; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
530 new TH2F("ptm_rap_gen_pbar", "MCTrack-gen antiproton;y;p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
532 new TH2F("ptm_rap_gen_d", "MCTrack-gen deuteron;y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
534 new TH2F("ptm_rap_gen_t", "MCTrack-gen triton; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
536 new TH2F("ptm_rap_gen_h", "MCTrack-gen 3he; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
538 new TH2F("ptm_rap_gen_a", "MCTrack-gen alpha; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
540 new TH2F("ptm_rap_gen_imf", "MCTrack-gen imf; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
541
542 Float_t v1_nbx = 20.;
543 Float_t v1_nby = 20.;
544 Float_t yvmax = 1.3;
545
547 new TH2F("v1_rap_gen_pip", "MCTrack-gen pi-plus; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
549 new TH2F("v1_rap_gen_pim", "MCTrack-gen pi-minus;y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
550 fa_v1_rap_gen_kp = new TH2F("v1_rap_gen_kp", "MCTrack-gen k-plus; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
551 fa_v1_rap_gen_km = new TH2F("v1_rap_gen_km", "MCTrack-gen k-minus; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
552 fa_v1_rap_gen_p = new TH2F("v1_rap_gen_p", "MCTrack-gen proton; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
554 new TH2F("v1_rap_gen_pbar", "MCTrack-gen antiproton;y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
555 fa_v1_rap_gen_d = new TH2F("v1_rap_gen_d", "MCTrack-gen deuteron;y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
556 fa_v1_rap_gen_t = new TH2F("v1_rap_gen_t", "MCTrack-gen triton; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
557 fa_v1_rap_gen_h = new TH2F("v1_rap_gen_h", "MCTrack-gen 3he; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
558 fa_v1_rap_gen_a = new TH2F("v1_rap_gen_a", "MCTrack-gen alpha; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
560 new TH2F("v1_rap_gen_imf", "MCTrack-gen imf; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
561
563 new TH2F("v2_rap_gen_pip", "MCTrack-gen pi-plus; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
565 new TH2F("v2_rap_gen_pim", "MCTrack-gen pi-minus;y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
566 fa_v2_rap_gen_kp = new TH2F("v2_rap_gen_kp", "MCTrack-gen k-plus; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
567 fa_v2_rap_gen_km = new TH2F("v2_rap_gen_km", "MCTrack-gen k-minus; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
568 fa_v2_rap_gen_p = new TH2F("v2_rap_gen_p", "MCTrack-gen proton; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
570 new TH2F("v2_rap_gen_pbar", "MCTrack-gen antiproton;y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
571 fa_v2_rap_gen_d = new TH2F("v2_rap_gen_d", "MCTrack-gen deuteron;y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
572 fa_v2_rap_gen_t = new TH2F("v2_rap_gen_t", "MCTrack-gen triton; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
573 fa_v2_rap_gen_h = new TH2F("v2_rap_gen_h", "MCTrack-gen 3he; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
574 fa_v2_rap_gen_a = new TH2F("v2_rap_gen_a", "MCTrack-gen alpha; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
576 new TH2F("v2_rap_gen_imf", "MCTrack-gen imf; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
577
578 // TofPoint level
579
581 new TH2F("ptm_rap_poi_pip", "MCTrack-poi pi-plus; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
583 new TH2F("ptm_rap_poi_pim", "MCTrack-poi pi-minus;y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
585 new TH2F("ptm_rap_poi_kp", "MCTrack-poi k-plus; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
587 new TH2F("ptm_rap_poi_km", "MCTrack-poi k-minus; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
589 new TH2F("ptm_rap_poi_p", "MCTrack-poi proton; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
591 new TH2F("ptm_rap_poi_pbar", "MCTrack-poi antiproton;y;p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
593 new TH2F("ptm_rap_poi_d", "MCTrack-poi deuteron;y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
595 new TH2F("ptm_rap_poi_t", "MCTrack-poi triton; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
597 new TH2F("ptm_rap_poi_h", "MCTrack-poi 3he; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
599 new TH2F("ptm_rap_poi_a", "MCTrack-poi alpha; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
601 new TH2F("ptm_rap_poi_imf", "MCTrack-poi imf; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
602
604 new TH2F("v1_rap_poi_pip", "MCTrack-poi pi-plus; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
606 new TH2F("v1_rap_poi_pim", "MCTrack-poi pi-minus;y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
607 fa_v1_rap_poi_kp = new TH2F("v1_rap_poi_kp", "MCTrack-poi k-plus; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
608 fa_v1_rap_poi_km = new TH2F("v1_rap_poi_km", "MCTrack-poi k-minus; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
609 fa_v1_rap_poi_p = new TH2F("v1_rap_poi_p", "MCTrack-poi proton; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
611 new TH2F("v1_rap_poi_pbar", "MCTrack-poi antiproton;y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
612 fa_v1_rap_poi_d = new TH2F("v1_rap_poi_d", "MCTrack-poi deuteron;y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
613 fa_v1_rap_poi_t = new TH2F("v1_rap_poi_t", "MCTrack-poi triton; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
614 fa_v1_rap_poi_h = new TH2F("v1_rap_poi_h", "MCTrack-poi 3he; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
615 fa_v1_rap_poi_a = new TH2F("v1_rap_poi_a", "MCTrack-poi alpha; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
617 new TH2F("v1_rap_poi_imf", "MCTrack-poi imf; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
618
620 new TH2F("v2_rap_poi_pip", "MCTrack-poi pi-plus; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
622 new TH2F("v2_rap_poi_pim", "MCTrack-poi pi-minus;y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
623 fa_v2_rap_poi_kp = new TH2F("v2_rap_poi_kp", "MCTrack-poi k-plus; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
624 fa_v2_rap_poi_km = new TH2F("v2_rap_poi_km", "MCTrack-poi k-minus; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
625 fa_v2_rap_poi_p = new TH2F("v2_rap_poi_p", "MCTrack-poi proton; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
627 new TH2F("v2_rap_poi_pbar", "MCTrack-poi antiproton;y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
628 fa_v2_rap_poi_d = new TH2F("v2_rap_poi_d", "MCTrack-poi deuteron;y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
629 fa_v2_rap_poi_t = new TH2F("v2_rap_poi_t", "MCTrack-poi triton; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
630 fa_v2_rap_poi_h = new TH2F("v2_rap_poi_h", "MCTrack-poi 3he; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
631 fa_v2_rap_poi_a = new TH2F("v2_rap_poi_a", "MCTrack-poi alpha; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
633 new TH2F("v2_rap_poi_imf", "MCTrack-poi imf; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
634
635 // TofHit level
636
638 new TH2F("ptm_rap_hit_pip", "MCTrack-hit pi-plus; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
640 new TH2F("ptm_rap_hit_pim", "MCTrack-hit pi-minus;y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
642 new TH2F("ptm_rap_hit_kp", "MCTrack-hit k-plus; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
644 new TH2F("ptm_rap_hit_km", "MCTrack-hit k-minus; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
646 new TH2F("ptm_rap_hit_p", "MCTrack-hit proton; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
648 new TH2F("ptm_rap_hit_pbar", "MCTrack-hit antiproton;y;p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
650 new TH2F("ptm_rap_hit_d", "MCTrack-hit deuteron;y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
652 new TH2F("ptm_rap_hit_t", "MCTrack-hit triton; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
654 new TH2F("ptm_rap_hit_h", "MCTrack-hit 3he; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
656 new TH2F("ptm_rap_hit_a", "MCTrack-hit alpha; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
658 new TH2F("ptm_rap_hit_imf", "MCTrack-hit imf; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
659
661 new TH2F("v1_rap_hit_pip", "MCTrack-hit pi-plus; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
663 new TH2F("v1_rap_hit_pim", "MCTrack-hit pi-minus;y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
664 fa_v1_rap_hit_kp = new TH2F("v1_rap_hit_kp", "MCTrack-hit k-plus; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
665 fa_v1_rap_hit_km = new TH2F("v1_rap_hit_km", "MCTrack-hit k-minus; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
666 fa_v1_rap_hit_p = new TH2F("v1_rap_hit_p", "MCTrack-hit proton; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
668 new TH2F("v1_rap_hit_pbar", "MCTrack-hit antiproton;y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
669 fa_v1_rap_hit_d = new TH2F("v1_rap_hit_d", "MCTrack-hit deuteron;y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
670 fa_v1_rap_hit_t = new TH2F("v1_rap_hit_t", "MCTrack-hit triton; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
671 fa_v1_rap_hit_h = new TH2F("v1_rap_hit_h", "MCTrack-hit 3he; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
672 fa_v1_rap_hit_a = new TH2F("v1_rap_hit_a", "MCTrack-hit alpha; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
674 new TH2F("v1_rap_hit_imf", "MCTrack-hit imf; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
675
677 new TH2F("v2_rap_hit_pip", "MCTrack-hit pi-plus; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
679 new TH2F("v2_rap_hit_pim", "MCTrack-hit pi-minus;y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
680 fa_v2_rap_hit_kp = new TH2F("v2_rap_hit_kp", "MCTrack-hit k-plus; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
681 fa_v2_rap_hit_km = new TH2F("v2_rap_hit_km", "MCTrack-hit k-minus; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
682 fa_v2_rap_hit_p = new TH2F("v2_rap_hit_p", "MCTrack-hit proton; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
684 new TH2F("v2_rap_hit_pbar", "MCTrack-hit antiproton;y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
685 fa_v2_rap_hit_d = new TH2F("v2_rap_hit_d", "MCTrack-hit deuteron;y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
686 fa_v2_rap_hit_t = new TH2F("v2_rap_hit_t", "MCTrack-hit triton; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
687 fa_v2_rap_hit_h = new TH2F("v2_rap_hit_h", "MCTrack-hit 3he; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
688 fa_v2_rap_hit_a = new TH2F("v2_rap_hit_a", "MCTrack-hit alpha; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
690 new TH2F("v2_rap_hit_imf", "MCTrack-hit imf; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
691
692 // GlobalTrack level
693
694 fa_plab_sts_pip = new TH1F("plab_sts_pip", "MCTrack-sts pi-plus; p_{Lab}(GeV/c)", 100, 0., 10.);
695 fa_plab_sts_pim = new TH1F("plab_sts_pim", "MCTrack-sts pi-minus; p_{Lab}(GeV/c)", 100, 0., 10.);
696 fa_plab_sts_kp = new TH1F("plab_sts_kp", "MCTrack-sts k-plus; p_{Lab}(GeV/c)", 100, 0., 10.);
697 fa_plab_sts_km = new TH1F("plab_sts_km", "MCTrack-sts k-minus; p_{Lab}(GeV/c)", 100, 0., 10.);
698 fa_plab_sts_p = new TH1F("plab_sts_p", "MCTrack-sts proton; p_{Lab}(GeV/c)", 100, 0., 10.);
699 fa_plab_sts_pbar = new TH1F("plab_sts_pbar", "MCTrack-sts pbar; p_{Lab}(GeV/c)", 100, 0., 10.);
700
702 new TH2F("ptm_rap_sts_pip", "MCTrack-sts pi-plus; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
704 new TH2F("ptm_rap_sts_pim", "MCTrack-sts pi-minus;y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
706 new TH2F("ptm_rap_sts_kp", "MCTrack-sts k-plus; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
708 new TH2F("ptm_rap_sts_km", "MCTrack-sts k-minus; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
710 new TH2F("ptm_rap_sts_p", "MCTrack-sts proton; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
712 new TH2F("ptm_rap_sts_pbar", "MCTrack-sts antiproton;y;p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
714 new TH2F("ptm_rap_sts_d", "MCTrack-sts deuteron;y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
716 new TH2F("ptm_rap_sts_t", "MCTrack-sts triton; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
718 new TH2F("ptm_rap_sts_h", "MCTrack-sts 3he; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
720 new TH2F("ptm_rap_sts_a", "MCTrack-sts alpha; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
722 new TH2F("ptm_rap_sts_imf", "MCTrack-sts imf; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
723
725 new TH2F("ptm_rap_glo_pip", "MCTrack-glo pi-plus; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
727 new TH2F("ptm_rap_glo_pim", "MCTrack-glo pi-minus;y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
729 new TH2F("ptm_rap_glo_kp", "MCTrack-glo k-plus; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
731 new TH2F("ptm_rap_glo_km", "MCTrack-glo k-minus; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
733 new TH2F("ptm_rap_glo_p", "MCTrack-glo proton; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
735 new TH2F("ptm_rap_glo_pbar", "MCTrack-glo antiproton;y;p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
737 new TH2F("ptm_rap_glo_d", "MCTrack-glo deuteron;y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
739 new TH2F("ptm_rap_glo_t", "MCTrack-glo triton; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
741 new TH2F("ptm_rap_glo_h", "MCTrack-glo 3he; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
743 new TH2F("ptm_rap_glo_a", "MCTrack-glo alpha; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
745 new TH2F("ptm_rap_glo_imf", "MCTrack-glo imf; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
746
748 new TH2F("v1_rap_glo_pip", "MCTrack-glo pi-plus; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
750 new TH2F("v1_rap_glo_pim", "MCTrack-glo pi-minus;y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
751 fa_v1_rap_glo_kp = new TH2F("v1_rap_glo_kp", "MCTrack-glo k-plus; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
752 fa_v1_rap_glo_km = new TH2F("v1_rap_glo_km", "MCTrack-glo k-minus; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
753 fa_v1_rap_glo_p = new TH2F("v1_rap_glo_p", "MCTrack-glo proton; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
755 new TH2F("v1_rap_glo_pbar", "MCTrack-glo antiproton;y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
756 fa_v1_rap_glo_d = new TH2F("v1_rap_glo_d", "MCTrack-glo deuteron;y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
757 fa_v1_rap_glo_t = new TH2F("v1_rap_glo_t", "MCTrack-glo triton; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
758 fa_v1_rap_glo_h = new TH2F("v1_rap_glo_h", "MCTrack-glo 3he; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
759 fa_v1_rap_glo_a = new TH2F("v1_rap_glo_a", "MCTrack-glo alpha; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
761 new TH2F("v1_rap_glo_imf", "MCTrack-glo imf; y; v_{1}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
762
764 new TH2F("v2_rap_glo_pip", "MCTrack-glo pi-plus; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
766 new TH2F("v2_rap_glo_pim", "MCTrack-glo pi-minus;y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
767 fa_v2_rap_glo_kp = new TH2F("v2_rap_glo_kp", "MCTrack-glo k-plus; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
768 fa_v2_rap_glo_km = new TH2F("v2_rap_glo_km", "MCTrack-glo k-minus; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
769 fa_v2_rap_glo_p = new TH2F("v2_rap_glo_p", "MCTrack-glo proton; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
771 new TH2F("v2_rap_glo_pbar", "MCTrack-glo antiproton;y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
772 fa_v2_rap_glo_d = new TH2F("v2_rap_glo_d", "MCTrack-glo deuteron;y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
773 fa_v2_rap_glo_t = new TH2F("v2_rap_glo_t", "MCTrack-glo triton; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
774 fa_v2_rap_glo_h = new TH2F("v2_rap_glo_h", "MCTrack-glo 3he; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
775 fa_v2_rap_glo_a = new TH2F("v2_rap_glo_a", "MCTrack-glo alpha; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
777 new TH2F("v2_rap_glo_imf", "MCTrack-glo imf; y; v_{2}", v1_nbx, -yvmax, yvmax, v1_nby, -1., 1.);
778
779 // xy - hit densities and rates
780 Int_t nbinx = 200;
781 Int_t nbiny = 200;
782 Float_t xrange = 750.;
783 Float_t yrange = 500.;
784 fwxy2 = nbinx * nbiny / 4. / xrange / yrange;
785
786 fa_xy_poi1 = new TH2F("xy_poi1", "TofPoint; ;", nbinx, -xrange, xrange, nbiny, -yrange, yrange);
787 fa_xy_poi2 = new TH2F("xy_poi2", "TofPoint /cm^{2}; ;", nbinx, -xrange, xrange, nbiny, -yrange, yrange);
788 fa_xy_poi3 = new TH2F("xy_poi3", "TofPoint /cm^{2}/s; ;", nbinx, -xrange, xrange, nbiny, -yrange, yrange);
789
790 fa_pv_poi = new TH2F("pv_poi", "TofPoint(all); velocity;momentum;", 100, 0., 32., 100., 0., 5.);
791 fa_tm_poi =
792 new TH2F("tm_poi", "Tofpoi(all); momentum (GeV/c); M_{ToF}*sign(Z) (GeV/c^{2});", 100, 0., 10., 200., -1.5, 2.5);
793 fa_tm_poiprim = new TH2F("tm_poiprim", "Tofpoi(primary); momentum; Tofmass", 100, 0., 10., 200., -1.5, 2.5);
794
795 fa_xy_hit1 = new TH2F("xy_hit1", "TofHit; ;", nbinx, -xrange, xrange, nbiny, -yrange, yrange);
796 fa_xy_hit2 = new TH2F("xy_hit2", "TofHit /cm^{2}; ;", nbinx, -xrange, xrange, nbiny, -yrange, yrange);
797 fa_xy_hit3 = new TH2F("xy_hit3", "TofHit /cm^{2}/s; ;", nbinx, -xrange, xrange, nbiny, -yrange, yrange);
798
799
800 fa_tof_hit = new TH1F("tof_hit", "TofHit(all); t (ns); counts;", 100, 0., 100.);
801 fa_dtof_hit = new TH1F("dtof_hit", "TofHit(all); #Deltat (ns); counts;", 100, -1., 1.);
802 fa_tof_hitprim = new TH1F("tof_hitprim", "TofHit(prim); t (ns); counts;", 100, 0., 100.);
803 fa_pv_hit = new TH2F("pv_hit", "TofHit(all); velocity; momentum;", 100, 0., 32., 100., 0., 5.);
804 fa_tm_hit =
805 new TH2F("tm_hit", "TofHit(all); momentum (GeV/c); M_{ToF}*sign(Z) (GeV/c^{2});", 100, 0., 10., 200., -1.5, 2.5);
806 fa_tm_hitprim = new TH2F("tm_hitprim", "TofHit(primary); momentum; Tofmass", 100, 0., 10., 200., -1.5, 2.5);
807 fa_tn_hit = new TH1F("tn_hit", "fastest TofHit(all); t (ns); counts;", 100, -1.0, 1.);
808 fa_t0_hit = new TH1F("t0_hit", "time zero; t0 (ns) ; counts;", 100, -0.5, 0.5);
809 fa_t0m_hit = new TH1F("t0m_hit", "average time zero; t0 (ns) ; counts;", 100, -0.1, 0.1);
810 fa_t0mn_hit = new TH1F("t0mn_hit", "average time zero hits; number of hits ; counts;", 100, 0., 500.);
811 fa_t0m_b_hit = new TH2F("t0m_b_hit", "average time zero; b (fm); t0 (ns) ; counts;", 28, 0., 14., 500, -0.3, 0.7);
813 new TH2F("t0mn_b_hit", "average time zero hits; b(fm); number of hits ; counts;", 28, 0., 14., 100, 0., 500.);
814
815 fa_t0m_f_hit = new TH1F("t0m_f_hit", "average time zero forward; t0 (ns) ; counts;", 100, -0.1, 0.1);
816 fa_t0mn_f_hit = new TH1F("t0mn_f_hit", "average time zero hits forward; number of hits ; counts;", 100, 0., 500.);
818 new TH2F("t0m_f_b_hit", "average time zero forward; b (fm); t0 (ns) ; counts;", 28, 0., 14., 500, -0.3, 0.7);
819 fa_t0mn_f_b_hit = new TH2F("t0mn_f_b_hit", "average time zero hits forward ; b(fm); number of hits ; counts;", 28, 0.,
820 14., 50, 0., 50.);
821
822 fa_t0m_nf_hit = new TH1F("t0m_nf_hit", "average time zero noforward; t0 (ns) ; counts;", 100, -0.1, 0.1);
823 fa_t0mn_nf_hit = new TH1F("t0mn_nf_hit", "average time zero hits noforward; number of hits ; counts;", 100, 0., 500.);
825 new TH2F("t0m_nf_b_hit", "average time zero noforward; b (fm); t0 (ns) ; counts;", 28, 0., 14., 500, -0.3, 0.7);
826 fa_t0mn_nf_b_hit = new TH2F("t0mn_nf_b_hit", "average time zero hits noforward ; b(fm); number of hits ; counts;", 28,
827 0., 14., 100, 0., 500.);
828
829 fa_dxx = new TH2F("dxx", "TofHit - TofPoint; x (cm); #Delta x (cm);", 100, -xrange, xrange, 50., -10., 10.);
830 fa_dxy = new TH2F("dxy", "TofHit - TofPoint; y (cm); #Delta x (cm);", 100, -yrange, yrange, 50., -10., 10.);
831 fa_dxz = new TH2F("dxz", "TofHit - TofPoint; z (cm); #Delta x (cm);", 100, 400., 650., 50., -10., 10.);
832 fa_dyx = new TH2F("dyx", "TofHit - TofPoint; x (cm); #Delta y (cm);", 100, -xrange, xrange, 50., -10., 10.);
833 fa_dyy = new TH2F("dyy", "TofHit - TofPoint; y (cm); #Delta y (cm);", 100, -yrange, yrange, 50., -10., 10.);
834 fa_dyz = new TH2F("dyz", "TofHit - TofPoint; z (cm); #Delta y (cm);", 100, 400., 650., 50., -10., 10.);
835 fa_dzx = new TH2F("dzx", "TofHit - TofPoint; x (cm); #Delta z (cm);", 100, -xrange, xrange, 50, -20., 20.);
836 fa_dzy = new TH2F("dzy", "TofHit - TofPoint; y (cm); #Delta z (cm);", 100, -yrange, yrange, 50, -20., 20.);
837 fa_dzz = new TH2F("dzz", "TofHit - TofPoint; z (cm); #Delta z (cm);", 100, 400., 650., 50, -20., 20.);
838
839 fa_hit_ch = new TH1F("hit_ch", "TofHits; channel; rate (Hz/s);", 50000, 0., 50000.);
840 fa_dhit_ch = new TH2F("dhit_ch", "Tof Double Hits; channel; counts;", 50000, 0., 50000., 2, 0.5, 2.5);
841
842 fa_xy_glo1 = new TH2F("xy_glo1", "GlobalTrack(all); x (cm); y(cm);", nbinx, -xrange, xrange, nbiny, -yrange, yrange);
844 new TH2F("xy_glo_pip", "GlobalTrack(pip); x (cm); y(cm);", nbinx, -xrange, xrange, nbiny, -yrange, yrange);
846 new TH2F("xy_glo_pim", "GlobalTrack(pim); x (cm); y(cm);", nbinx, -xrange, xrange, nbiny, -yrange, yrange);
847 fa_xy_glo_p = new TH2F("xy_glo_p", "GlobalTrack(p); x (cm); y(cm);", nbinx, -xrange, xrange, nbiny, -yrange, yrange);
849 new TH2F("xy_glo_pbar", "GlobalTrack(pbar); x (cm); y(cm);", nbinx, -xrange, xrange, nbiny, -yrange, yrange);
851 new TH2F("xy_glo_kp", "GlobalTrack(kp); x (cm); y(cm);", nbinx, -xrange, xrange, nbiny, -yrange, yrange);
853 new TH2F("xy_glo_km", "GlobalTrack(km); x (cm); y(cm);", nbinx, -xrange, xrange, nbiny, -yrange, yrange);
854 fa_xy_glo_d = new TH2F("xy_glo_d", "GlobalTrack(d); x (cm); y(cm);", nbinx, -xrange, xrange, nbiny, -yrange, yrange);
855 fa_xy_glo_t = new TH2F("xy_glo_t", "GlobalTrack(t); x (cm); y(cm);", nbinx, -xrange, xrange, nbiny, -yrange, yrange);
856 fa_xy_glo_h = new TH2F("xy_glo_h", "GlobalTrack(h); x (cm); y(cm);", nbinx, -xrange, xrange, nbiny, -yrange, yrange);
857 fa_xy_glo_a = new TH2F("xy_glo_a", "GlobalTrack(a); x (cm); y(cm);", nbinx, -xrange, xrange, nbiny, -yrange, yrange);
858
859 fa_TofTrackMul = new TH1F("TofTrackMul", "number of assigned TofTrack / global track; multiplicity; ", 50, 0., 50.);
860 fa_VtxB = new TH1F("VtxB", "Chi2 to primary vertex; ", 100, 0., 20.);
861
862 Float_t TMMIN = -1.5;
863 Float_t TMMAX = 2.5;
864 Int_t TMYBIN = 400;
865 fa_pv_glo = new TH2F("pv_glo", "GlobalTrack(all); velocity; momentum;", 100, 0., 32., 100., 0., 5.);
866 fa_tm_glo = new TH2F("tm_glo", "GlobalTrack(all); momentum (GeV/c); M_{ToF}*sign(Z) (GeV/c^{2});", 100, 0., 10.,
867 TMYBIN, TMMIN, TMMAX);
868 fa_chi2_mom_glo = new TH2F("chi2_mom_glo", "GlobalTrack(all); momentum; chi2;", 100, 0., 10., 100., 0., 50.);
870 new TH2F("chi2_mom_gloprim", "GlobalTrack(primaries); momentum; chi2;", 100, 0., 10., 100., 0., 50.);
871 fa_len_mom_glo = new TH2F("len_mom_glo", "GlobalTrack(all); momentum; len;", 100, 0., 10., 100., 0., 1500.);
872 fa_tm_gloprim = new TH2F("tm_gloprim", "GlobalTrack(primary); momentum; Tofmass", 100, 0., 10., TMYBIN, TMMIN, TMMAX);
873 fa_tm_glovtxb = new TH2F("tm_glovtxb", "GlobalTrack(vtxb); momentum (GeV/c); M_{ToF}*sign(Z) (GeV/c^{2});", 100, 0.,
874 10., TMYBIN, TMMIN, TMMAX);
876 new TH2F("tm_gloprimvtxb", "GlobalTrack(prim,vtxb); momentum; Tofmass", 100, 0., 10., TMYBIN, TMMIN, TMMAX);
877 fa_tm_glomis = new TH2F("tm_glomis", "GlobalTrack(mis); momentum (GeV/c); M_{ToF}*sign(Z) (GeV/c^{2});", 100, 0., 10.,
878 TMYBIN, TMMIN, TMMAX);
879 fa_tm_glo_pip = new TH2F("tm_glo_pip", "GlobalTrack(pip); momentum (GeV/c); M_{ToF}*sign(Z) (GeV/c^{2});", 100, 0.,
880 10., TMYBIN, TMMIN, TMMAX);
881 fa_tm_glo_pim = new TH2F("tm_glo_pim", "GlobalTrack(pim); momentum (GeV/c); M_{ToF}*sign(Z) (GeV/c^{2});", 100, 0.,
882 10., TMYBIN, TMMIN, TMMAX);
883 fa_tm_glo_kp = new TH2F("tm_glo_kp", "GlobalTrack(kp); momentum (GeV/c); M_{ToF}*sign(Z) (GeV/c^{2});", 100, 0., 10.,
884 TMYBIN, TMMIN, TMMAX);
885 fa_tm_glo_km = new TH2F("tm_glo_km", "GlobalTrack(km); momentum (GeV/c); M_{ToF}*sign(Z) (GeV/c^{2});", 100, 0., 10.,
886 TMYBIN, TMMIN, TMMAX);
887 fa_tm_glo_p = new TH2F("tm_glo_p", "GlobalTrack(p); momentum (GeV/c); M_{ToF}*sign(Z) (GeV/c^{2});", 100, 0., 10.,
888 TMYBIN, TMMIN, TMMAX);
889 fa_tm_glo_pbar = new TH2F("tm_glo_pbar", "GlobalTrack(pbar); momentum (GeV/c); M_{ToF}*sign(Z) (GeV/c^{2});", 100, 0.,
890 10., TMYBIN, TMMIN, TMMAX);
891 fa_tm_glo_d = new TH2F("tm_glo_d", "GlobalTrack(d); momentum (GeV/c); M_{ToF}*sign(Z) (GeV/c^{2});", 100, 0., 10.,
892 TMYBIN, TMMIN, TMMAX);
893 fa_tm_glo_t = new TH2F("tm_glo_t", "GlobalTrack(t); momentum (GeV/c); M_{ToF}*sign(Z) (GeV/c^{2});", 100, 0., 10.,
894 TMYBIN, TMMIN, TMMAX);
895 fa_tm_glo_h = new TH2F("tm_glo_h", "GlobalTrack(h); momentum (GeV/c); M_{ToF}*sign(Z) (GeV/c^{2});", 100, 0., 10.,
896 TMYBIN, TMMIN, TMMAX);
897 fa_tm_glo_a = new TH2F("tm_glo_a", "GlobalTrack(a); momentum (GeV/c); M_{ToF}*sign(Z) (GeV/c^{2});", 100, 0., 10.,
898 TMYBIN, TMMIN, TMMAX);
899
900 Double_t M2MIN = -0.4;
901 Double_t M2MAX = 1.4;
902 Int_t M2YBIN = 360;
903 fa_m2mom_glo = new TH2F("m2mom_glo", "GlobalTrack(all); p ^{.} sign(Z)(GeV); M_{ToF}^{2} (GeV^{2});", 200, -10., 10.,
904 M2YBIN, M2MIN, M2MAX);
905 fa_m2mom_glovtxb = new TH2F("m2mom_glovtxb", "GlobalTrack(vtxb); p ^{.} sign(Z)(GeV); M_{ToF}^{2} (GeV^{2});", 200,
906 -10., 10., M2YBIN, M2MIN, M2MAX);
907 fa_m2mom_gloprim = new TH2F("m2mom_gloprim", "GlobalTrack(prim); p ^{.} sign(Z)(GeV); M_{ToF}^{2} (GeV^{2});", 200,
908 -10., 10., M2YBIN, M2MIN, M2MAX);
910 new TH2F("m2mom_gloprimvtxb", "GlobalTrack(primvtxb); p ^{.} sign(Z)(GeV); M_{ToF}^{2} (GeV^{2});", 200, -10., 10.,
911 M2YBIN, M2MIN, M2MAX);
912 fa_m2mom_glo_pip = new TH2F("m2mom_glo_pip", "GlobalTrack(pip); p ^{.} sign(Z)(GeV); M_{ToF}^{2} (GeV^{2});", 200,
913 -10., 10., M2YBIN, M2MIN, M2MAX);
914 fa_m2mom_glo_pim = new TH2F("m2mom_glo_pim", "GlobalTrack(pim); p ^{.} sign(Z)(GeV); M_{ToF}^{2} (GeV^{2});", 200,
915 -10., 10., M2YBIN, M2MIN, M2MAX);
916 fa_m2mom_glo_kp = new TH2F("m2mom_glo_kp", "GlobalTrack(kp); p ^{.} sign(Z)(GeV); M_{ToF}^{2} (GeV^{2});", 200, -10.,
917 10., M2YBIN, M2MIN, M2MAX);
918 fa_m2mom_glo_km = new TH2F("m2mom_glo_km", "GlobalTrack(km); p ^{.} sign(Z)(GeV); M_{ToF}^{2} (GeV^{2});", 200, -10.,
919 10., M2YBIN, M2MIN, M2MAX);
920 fa_m2mom_glo_p = new TH2F("m2mom_glo_p", "GlobalTrack(p); p ^{.} sign(Z)(GeV); M_{ToF}^{2} (GeV^{2});", 200, -10.,
921 10., M2YBIN, M2MIN, M2MAX);
922 fa_m2mom_glo_pbar = new TH2F("m2mom_glo_pbar", "GlobalTrack(pbar); p ^{.} sign(Z)(GeV); M_{ToF}^{2} (GeV^{2});", 200,
923 -10., 10., M2YBIN, M2MIN, M2MAX);
924
925 fa_pMCmom_glo = new TH2F("pMCmom_glo", "GlobalTrack(all); momentum; p_{MC};", 100, 0., 10., 100, 0., 10.);
926 fa_pMCmom_glo_pip = new TH2F("pMCmom_glo_pip", "GlobalTrack(pip); momentum; p_{MC};", 100, 0., 10., 100, 0., 10.);
927 fa_pMCmom_glo_pim = new TH2F("pMCmom_glo_pim", "GlobalTrack(pim); momentum; p_{MC};", 100, 0., 10., 100, 0., 10.);
928 fa_pMCmom_glo_kp = new TH2F("pMCmom_glo_kp", "GlobalTrack(kp); momentum; p_{MC};", 100, 0., 10., 100, 0., 10.);
929 fa_pMCmom_glo_km = new TH2F("pMCmom_glo_km", "GlobalTrack(km); momentum; p_{MC};", 100, 0., 10., 100, 0., 10.);
930 fa_pMCmom_glo_p = new TH2F("pMCmom_glo_p", "GlobalTrack(p); momentum; p_{MC};", 100, 0., 10., 100, 0., 10.);
931 fa_pMCmom_glo_pbar = new TH2F("pMCmom_glo_pbar", "GlobalTrack(pbar); momentum; p_{MC};", 100, 0., 10., 100, 0., 10.);
932 fa_pMCmom_glo_d = new TH2F("pMCmom_glo_d", "GlobalTrack(d); momentum; p_{MC};", 100, 0., 10., 100, 0., 10.);
933 fa_pMCmom_glo_t = new TH2F("pMCmom_glo_t", "GlobalTrack(t); momentum; p_{MC};", 100, 0., 10., 100, 0., 10.);
934 fa_pMCmom_glo_h = new TH2F("pMCmom_glo_h", "GlobalTrack(h); momentum; p_{MC};", 100, 0., 10., 100, 0., 10.);
935 fa_pMCmom_glo_a = new TH2F("pMCmom_glo_a", "GlobalTrack(a); momentum; p_{MC};", 100, 0., 10., 100, 0., 10.);
936
937 Double_t LenDisMax = 20.;
939 new TH2F("LenDismom_glo", "GlobalTrack(all); momentum; len-dis;", 100, 0., 10., 100., -LenDisMax, LenDisMax);
941 new TH2F("LenDismom_glo_pip", "GlobalTrack(pip); momentum; len-dis;", 100, 0., 10., 100., -LenDisMax, LenDisMax);
943 new TH2F("LenDismom_glo_pim", "GlobalTrack(pim); momentum; len-dis;", 100, 0., 10., 100., -LenDisMax, LenDisMax);
945 new TH2F("LenDismom_glo_kp", "GlobalTrack(kp); momentum; len-dis;", 100, 0., 10., 100., -LenDisMax, LenDisMax);
947 new TH2F("LenDismom_glo_km", "GlobalTrack(km); momentum; len-dis;", 100, 0., 10., 100., -LenDisMax, LenDisMax);
949 new TH2F("LenDismom_glo_p", "GlobalTrack(p); momentum; len-dis;", 100, 0., 10., 100., -LenDisMax, LenDisMax);
951 new TH2F("LenDismom_glo_pbar", "GlobalTrack(pbar); momentum; len-dis;", 100, 0., 10., 100., -LenDisMax, LenDisMax);
953 new TH2F("LenDismom_glo_d", "GlobalTrack(d); momentum; len-dis;", 100, 0., 10., 100., -LenDisMax, LenDisMax);
955 new TH2F("LenDismom_glo_t", "GlobalTrack(t); momentum; len-dis;", 100, 0., 10., 100., -LenDisMax, LenDisMax);
957 new TH2F("LenDismom_glo_h", "GlobalTrack(h); momentum; len-dis;", 100, 0., 10., 100., -LenDisMax, LenDisMax);
959 new TH2F("LenDismom_glo_a", "GlobalTrack(a); momentum; len-dis;", 100, 0., 10., 100., -LenDisMax, LenDisMax);
960
961 fa_LenMcLenGlomom_glo = new TH2F("LenMcLenGlomom_glo", "GlobalTrack(all); momentum [GeV]; len glo - len MC pnt[cm];",
962 100, 0., 10., 400., -LenDisMax, LenDisMax);
963 fa_LenMcDismom_glo = new TH2F("LenMcDismom_glo", "GlobalTrack(all); momentum [GeV]; len MC pnt - dis [cm];", 100, 0.,
964 10., 400., -LenDisMax, LenDisMax);
965
966 Int_t WYBIN = 100;
967 Float_t WYMAX = 20.;
968 fa_w_mom_glo = new TH2F("w_mom_glo", "GlobalTrack(all); momentum; weight;", 100, 0., 10., WYBIN, 0., WYMAX);
969 fa_w_mom_glo_pip = new TH2F("w_mom_glo_pip", "GlobalTrack(pip); momentum; weight;", 100, 0., 10., WYBIN, 0., WYMAX);
970 fa_w_mom_glo_pim = new TH2F("w_mom_glo_pim", "GlobalTrack(pim); momentum; weight;", 100, 0., 10., WYBIN, 0., WYMAX);
971 fa_w_mom_glo_kp = new TH2F("w_mom_glo_kp", "GlobalTrack(kp); momentum; weight;", 100, 0., 10., WYBIN, 0., WYMAX);
972 fa_w_mom_glo_km = new TH2F("w_mom_glo_km", "GlobalTrack(km); momentum; weight;", 100, 0., 10., WYBIN, 0., WYMAX);
973 fa_w_mom_glo_p = new TH2F("w_mom_glo_p", "GlobalTrack(p); momentum; weight;", 100, 0., 10., WYBIN, 0., WYMAX);
975 new TH2F("w_mom_glo_pbar", "GlobalTrack(pbar); momentum; weight;", 100, 0., 10., WYBIN, 0., WYMAX);
976 fa_w_mom_glo_d = new TH2F("w_mom_glo_d", "GlobalTrack(d); momentum; weight;", 100, 0., 10., WYBIN, 0., WYMAX);
977 fa_w_mom_glo_t = new TH2F("w_mom_glo_t", "GlobalTrack(t); momentum; weight;", 100, 0., 10., WYBIN, 0., WYMAX);
978 fa_w_mom_glo_h = new TH2F("w_mom_glo_h", "GlobalTrack(h); momentum; weight;", 100, 0., 10., WYBIN, 0., WYMAX);
979 fa_w_mom_glo_a = new TH2F("w_mom_glo_a", "GlobalTrack(a); momentum; weight;", 100, 0., 10., WYBIN, 0., WYMAX);
980 fa_w_mom_glomis = new TH2F("w_mom_glomis", "GlobalTrack(mis); momentum; weight;", 100, 0., 10., WYBIN, 0., WYMAX);
981
982 //centrality
984 new TH2F("mul_b_gen", "Centrality - gen;impact parameter b; multiplicity;", 15, 0., 15., 100, 0., 10000.);
986 new TH2F("mul_b_poi", "Centrality - poi;impact parameter b; multiplicity;", 15, 0., 15., 100, 0., 10000.);
988 new TH2F("mul_b_hit", "Centrality - hit;impact parameter b; multiplicity;", 15, 0., 15., 100, 0., 2000.);
990 new TH2F("mul_b_glo", "Centrality - glo;impact parameter b; multiplicity;", 15, 0., 15., 100, 0., 2000.);
992 new TH2F("mul_b_had", "Centrality - had;impact parameter b; multiplicity;", 15, 0., 15., 100, 0., 2000.);
993
994 //reaction plane
995 fa_cdrp_b_gen = new TH2F("cdrp_b_gen", "reaction plane resolution - gen;impact parameter b; cos#Delta#phi_{rp};", 15,
996 0., 15., 20, -1., 1.);
998 new TH2F("drp_b_gen", "#Delta#phi-reaction plane - gen ;impact parameter b;", 15, 0., 15., 36, -180., 180.);
1000 new TH2F("phirp_b_gen", "#phi_{reaction plane} - gen ;impact parameter b;", 15, 0., 15., 36, -180., 180.);
1002 new TH2F("phgrp_b_gen", "#phi_{G reaction plane} - gen ;impact parameter b;", 15, 0., 15., 36, -180., 180.);
1003 fa_phphrp_gen = new TH2F("phphrp_gen", "#phi_#phi - gen ;#phi_{rec}; #phi_{gen} ;", 36, -180., 180., 36, -180., 180.);
1004 fa_delrp_b_gen = new TH2F("delrp_b_gen",
1005 "#Delta#phi_{G}-reaction plane - gen ;impact "
1006 "parameter b;#phi_{rec}-#phi_{gen}",
1007 15, 0., 15., 36, -180., 180.);
1008 fa_delrp_b_poi = new TH2F("delrp_b_poi",
1009 "#Delta#phi_{G}-reaction plane - poi ;impact "
1010 "parameter b;#phi_{rec}-#phi_{gen}",
1011 15, 0., 15., 36, -180., 180.);
1012 fa_delrp_b_hit = new TH2F("delrp_b_hit",
1013 "#Delta#phi_{G}-reaction plane - hit ;impact "
1014 "parameter b;#phi_{rec}-#phi_{gen}",
1015 15, 0., 15., 36, -180., 180.);
1016 fa_delrp_b_glo = new TH2F("delrp_b_glo",
1017 "#Delta#phi_{G}-reaction plane - glo ;impact "
1018 "parameter b;#phi_{rec}-#phi_{gen}",
1019 15, 0., 15., 36, -180., 180.);
1020
1021 fa_cdelrp_b_gen = new TH2F("cdelrp_b_gen",
1022 "reaction plane resolution - gen;impact parameter "
1023 "b;cos(#phi_{rec}-#phi_{gen})",
1024 15, 0., 15., 20, -1., 1.);
1025 fa_cdelrp_b_poi = new TH2F("cdelrp_b_poi",
1026 "reaction plane resolution - poi;impact parameter "
1027 "b;cos(#phi_{rec}-#phi_{gen})",
1028 15, 0., 15., 20, -1., 1.);
1029 fa_cdelrp_b_hit = new TH2F("cdelrp_b_hit",
1030 "reaction plane resolution - hit;impact parameter "
1031 "b;cos(#phi_{rec}-#phi_{gen})",
1032 15, 0., 15., 20, -1., 1.);
1033 fa_cdelrp_b_glo = new TH2F("cdelrp_b_glo",
1034 "reaction plane resolution - glo;impact parameter "
1035 "b;cos(#phi_{rec}-#phi_{gen})",
1036 15, 0., 15., 20, -1., 1.);
1037 fa_cdelrp_b_had = new TH2F("cdelrp_b_had",
1038 "reaction plane resolution - had;impact parameter "
1039 "b;cos(#phi_{rec}-#phi_{gen})",
1040 15, 0., 15., 20, -1., 1.);
1041
1042 fa_cdrp_b_poi = new TH2F("cdrp_b_poi", "reaction plane resolution - poi;impact parameter b; cos#Delta#phi_{rp};", 15,
1043 0., 15., 20, -1., 1.);
1044 fa_drp_b_poi =
1045 new TH2F("drp_b_poi", "#Delta#phi-reaction plane - poi ;impact parameter b;", 15, 0., 15., 36, -180., 180.);
1046
1047 fa_cdrp_b_hit = new TH2F("cdrp_b_hit", "reaction plane resolution - hit;impact parameter b; cos#Delta#phi_{rp};", 15,
1048 0., 15., 20, -1., 1.);
1049 fa_drp_b_hit =
1050 new TH2F("drp_b_hit", "#Delta#phi-reaction plane - hit ;impact parameter b;", 15, 0., 15., 36, -180., 180.);
1051
1052 fa_cdrp_b_glo = new TH2F("cdrp_b_glo", "reaction plane resolution - glo;impact parameter b; cos#Delta#phi_{rp};", 15,
1053 0., 15., 20, -1., 1.);
1054 fa_drp_b_glo =
1055 new TH2F("drp_b_glo", "#Delta#phi-reaction plane - glo ;impact parameter b;", 15, 0., 15., 36, -180., 180.);
1056
1057 fa_cdrp_b_had = new TH2F("cdrp_b_had", "reaction plane resolution - had;impact parameter b; cos#Delta#phi_{rp};", 15,
1058 0., 15., 20, -1., 1.);
1059 fa_drp_b_had =
1060 new TH2F("drp_b_had", "#Delta#phi-reaction plane - had ;impact parameter b;", 15, 0., 15., 36, -180., 180.);
1061
1062 fa_phirp_gen = new TH1F("phirp_gen", "#phi_{reaction plane} - gen ;#phi_{RPgen};", 72, -180., 180.);
1063 fa_phirp_poi = new TH1F("phirp_poi", "#phi_{reaction plane} - poi ;#phi_{RPpoi};", 72, -180., 180.);
1064 fa_phirp_hit = new TH1F("phirp_hit", "#phi_{reaction plane} - hit ;#phi_{RPhit};", 72, -180., 180.);
1065 fa_phirp_glo = new TH1F("phirp_glo", "#phi_{reaction plane} - glo ;#phi_{RPglo};", 72, -180., 180.);
1066 fa_phirp_had = new TH1F("phirp_had", "#phi_{reaction plane} - had ;#phi_{RPhad};", 72, -180., 180.);
1067
1068 fa_phirps_gen = new TH1F("phirps_gen", "#phi_{reaction plane} - gen ;#phi_{sRPgen};", 72, -180., 180.);
1069 fa_phirps_poi = new TH1F("phirps_poi", "#phi_{reaction plane} - poi ;#phi_{sRPpoi};", 72, -180., 180.);
1070 fa_phirps_hit = new TH1F("phirps_hit", "#phi_{reaction plane} - hit ;#phi_{sRPhit};", 72, -180., 180.);
1071 fa_phirps_glo = new TH1F("phirps_glo", "#phi_{reaction plane} - glo ;#phi_{sRPglo};", 72, -180., 180.);
1072
1073 fhwdist = new TH2F("hwdist", "matching wdist; p (GeV/c); dist;", 100, 0., 10., 50, 0., 10.);
1074 fhwmindelmass = new TH2F("hwmindelmass", "matching wmindelmass ; p (GeV/c); #Delta M;", 100, 0., 10., 50, 0., 1.);
1075 fhwminlen = new TH2F("hwminlen", "matching wminlen ; p (GeV/c); MinPathlength;", 100, 0., 10., 50, 0., 2.);
1076 fhwdelp = new TH2F("hwdelp", "matching wdelp ; p (GeV/c); #Delta p/p;", 100, 0., 10., 50, 0., 1.);
1077
1078 fhTofTrkDx = new TH1F("hTofTrkDx", " x - position resolution ; #deltax;", 50, -1., 1.);
1079 fhTofTrkDy = new TH1F("hTofTrkDy", " y - position resolution ; #deltay;", 50, -1., 1.);
1080 fhTofTrkDxsel = new TH1F("hTofTrkDxsel", " x - position resolution ; #deltax;", 50, -1., 1.);
1081 fhTofTrkDysel = new TH1F("hTofTrkDysel", " y - position resolution ; #deltay;", 50, -1., 1.);
1082
1083 cout << "CbmHadronAnalysis::CreateHistogramms: histograms booked in directory " << gDirectory->GetName() << endl;
1084}
1085
1086// ------------------------------------------------------------------
1088{ // Open PDF file and get histogramms for PID
1089 if (fPdfFileName.IsNull()) return kSUCCESS;
1090
1092 TFile* oldFile = gFile;
1093 TDirectory* oldDir = gDirectory;
1094
1095 TFile* pdfFile = new TFile(fPdfFileName);
1096
1098 gFile = oldFile;
1099 gDirectory = oldDir;
1100
1101 if (NULL == pdfFile) {
1102 cout << "-E- CbmHadronAnalysis::ReadPdfFile : "
1103 << "file " << fPdfFileName << " does not exist!" << endl;
1104 return kSUCCESS;
1105 }
1106 return kSUCCESS;
1107}
1108
1109// ------------------------------------------------------------------
1111{ // Open file and get histogramms for RP corrections
1112 if (fFlowFileName.IsNull()) return kSUCCESS;
1113
1115 TFile* oldFile = gFile;
1116 TDirectory* oldDir = gDirectory;
1117
1118 fflowFile = new TFile(fFlowFileName);
1119
1121 gFile = oldFile;
1122 gDirectory = oldDir;
1123
1124 if (NULL == fflowFile) {
1125 cout << "-E- CbmHadronAnalysis::ReadFlowFile : "
1126 << "file " << fFlowFileName << " does not exist!" << endl;
1127 return kSUCCESS;
1128 }
1129 cout << "-I- CbmHadronAnalysis::ReadFlowFile : RP corrections from " << fFlowFileName << endl;
1130
1131 return kSUCCESS;
1132}
1133
1134// ------------------------------------------------------------------
1136{
1137 // Task initialisation
1138 FairRootManager* rootMgr = FairRootManager::Instance();
1139 if (NULL == rootMgr) {
1140 cout << "-E- CbmHadronAnalysis::Init : ROOT manager is not instantiated." << endl;
1141 return kFATAL;
1142 }
1143
1144 fEventsColl = dynamic_cast<TClonesArray*>(rootMgr->GetObject("Event"));
1145 if (NULL == fEventsColl) fEventsColl = dynamic_cast<TClonesArray*>(rootMgr->GetObject("CbmEvent"));
1146
1147 if (NULL == fEventsColl) LOG(info) << "CbmEvent not found in input file, assume eventwise input";
1148 else
1149 LOG(info) << "CbmEvent found in input file, timebased analysis";
1150
1151 fMCEventHeader = (FairMCEventHeader*) rootMgr->GetObject("MCEventHeader.");
1152 if (NULL == fMCEventHeader) { cout << "-W- CbmHadronAnalysis::Init : no MC Header Info" << endl; }
1153
1154 fTofPoints = (TClonesArray*) rootMgr->GetObject("TofPoint");
1155 if (NULL == fTofPoints) { cout << "-W- CbmHadronAnalysis::Init : no TOF point array!" << endl; }
1156
1158 fDigiMan->Init();
1160 LOG(error) << GetName() << ": No input digis!";
1161 return kFATAL;
1162 }
1163 else
1164 LOG(info) << "DigiManager has Tof Digis";
1165
1166 if (!fDigiMan->IsMatchPresent(ECbmModuleId::kTof)) { LOG(warn) << GetName() << ": No TofDigiMatches!"; }
1167 else
1168 LOG(info) << "DigiManager has TofDigiMatches";
1169
1171 LOG(error) << GetName() << ": No STS input digis!";
1172 return kFATAL;
1173 }
1174 else
1175 LOG(info) << "DigiManager has Sts Digis";
1176
1177 if (!fDigiMan->IsMatchPresent(ECbmModuleId::kSts)) { LOG(warn) << GetName() << ": No StsDigiMatches!"; }
1178 else
1179 LOG(info) << "DigiManager has StsDigiMatches";
1180
1181 if (!fDigiMan->IsMatchPresent(ECbmModuleId::kTof)) { LOG(warn) << GetName() << ": No TofDigiMatches!"; }
1182 else
1183 LOG(info) << "DigiManager has TofDigiMatches";
1184
1186 LOG(error) << GetName() << ": No STS input digis!";
1187 return kFATAL;
1188 }
1189 else
1190 LOG(info) << "DigiManager has Sts Digis";
1191
1192 if (!fDigiMan->IsMatchPresent(ECbmModuleId::kSts)) { LOG(warn) << GetName() << ": No StsDigiMatches!"; }
1193 else
1194 LOG(info) << "DigiManager has StsDigiMatches";
1195
1196 /*
1197 fTofDigis = (TClonesArray *) rootMgr->GetObject("TofDigi");
1198
1199 if( NULL == fTofDigis)
1200 fTofDigis = (TClonesArray *) rootMgr->GetObject("CbmTofDigiExp");
1201
1202 if( NULL == fTofDigis)
1203 fTofDigis = (TClonesArray *) rootMgr->GetObject("CbmTofDigi");
1204 if( NULL == fTofDigis)
1205 {
1206 LOG(error)<<"CbmHadronAnalysis::RegisterInputs => Could not get the TofDigi TClonesArray!!!";
1207 } // if( NULL == fTofDigis)
1208
1209 fTofDigiMatchPointsColl = (TClonesArray*) rootMgr->GetObject("TofDigiMatch");
1210 if (NULL == fTofDigiMatchPointsColl) { cout << "-I- CbmHadronAnalysis::Init : no TofDigiMatch array!" << endl; }
1211
1212 fStsDigis = (TClonesArray*) rootMgr->GetObject("StsDigi");
1213 if (NULL == fStsDigis) { cout << "-W- CbmHadronAnalysis::Init : no STS Digi array!" << endl; }
1214 fStsDigiMatchColl = (TClonesArray*) rootMgr->GetObject("StsDigiMatch");
1215 if (NULL == fStsDigiMatchColl) {
1216 cout << "-W- CbmHadronAnalysis::Init : no STS DigiMatch array!" << endl;
1217 }
1218 */
1219
1220 fTofDigiMatchColl = (TClonesArray*) rootMgr->GetObject("TofHitCalDigiMatch");
1221 if (NULL == fTofDigiMatchColl) {
1222 cout << "-I- CbmHadronAnalysis::Init : no TOF HitCalDigiMatch array!" << endl;
1223 }
1224
1225 fTofHitsColl = (TClonesArray*) rootMgr->GetObject("TofCalHit");
1226 if (NULL == fTofHitsColl) {
1227 cout << "-W- CbmHadronAnalysis::Init : no TOFCalHit array!" << endl;
1228 fTofHitsColl = (TClonesArray*) rootMgr->GetObject("TofHit");
1229 if (NULL == fTofHitsColl) { cout << "-W- CbmHadronAnalysis::Init : no TOF Hit array!" << endl; }
1230 }
1231
1232 fTofTracks = (TClonesArray*) rootMgr->GetObject("TofTrack");
1233 if (NULL == fTofTracks) { cout << "-W- CbmHadronAnalysis::Init : no Tof Track array!" << endl; }
1234
1235 fTrdPoints = (TClonesArray*) rootMgr->GetObject("TrdPoint");
1236 if (NULL == fTrdPoints) { cout << "-W- CbmHadronAnalysis::Init : no TRD point array!" << endl; }
1237
1238 fTrdHits = (TClonesArray*) rootMgr->GetObject("TrdHit");
1239 if (NULL == fTrdHits) { cout << "-W- CbmHadronAnalysis::Init : no TRD Hit array!" << endl; }
1240
1241 fGlobalTracks = (TClonesArray*) rootMgr->GetObject("GlobalTrack");
1242 if (NULL == fGlobalTracks) { cout << "-W- CbmHadronAnalysis::Init : no Global Track array!" << endl; }
1243
1244 fHadrons = (TClonesArray*) rootMgr->GetObject("Hadron");
1245 if (NULL == fHadrons) { cout << "-W- CbmHadronAnalysis::Init : no Hadron array!" << endl; }
1246
1247 fStsTracks = (TClonesArray*) rootMgr->GetObject("StsTrack");
1248 if (NULL == fStsTracks) { cout << "-W- CbmHadronAnalysis::Init : no STS Track array!" << endl; }
1249 fStsHitsColl = (TClonesArray*) rootMgr->GetObject("StsHit");
1250 if (NULL == fStsHitsColl) { cout << "-W- CbmHadronAnalysis::Init : no STS Hit array!" << endl; }
1251 fStsClusters = (TClonesArray*) rootMgr->GetObject("StsCluster");
1252 if (NULL == fStsClusters) { cout << "-W- CbmHadronAnalysis::Init : no STS Cluster array!" << endl; }
1253
1254 // Get pointer to PrimaryVertex object from IOManager if it exists
1255 // The old name for the object is "PrimaryVertex" the new one
1256 // "PrimaryVertex." Check first for the new name
1257 fPrimVertex = dynamic_cast<CbmVertex*>(rootMgr->GetObject("PrimaryVertex."));
1258 if (nullptr == fPrimVertex) { fPrimVertex = dynamic_cast<CbmVertex*>(rootMgr->GetObject("PrimaryVertex")); }
1259 if (nullptr == fPrimVertex) { LOG(warn) << "No primary vertex"; }
1260
1261 if (NULL == fEventsColl) {
1264 }
1265 else { // for time based analysis generate temporary Event TClonesArrays
1266 fTofHits = new TClonesArray("CbmTofHit", 100);
1267 fStsHits = new TClonesArray("CbmStsHit", 100);
1268 }
1269
1270 // --- MC data manager
1271 CbmMCDataManager* mcManager = (CbmMCDataManager*) rootMgr->GetObject("MCDataManager");
1272
1273 if (NULL != mcManager) {
1274 // --- Data arrays
1275 fMCTracks = mcManager->InitBranch("MCTrack");
1276 fStsPoints = mcManager->InitBranch("StsPoint");
1277 }
1278
1279 fMCTracksColl = (TClonesArray*) rootMgr->GetObject("MCTrack");
1280 if (NULL == fMCTracksColl) { cout << "-W- CbmHadronAnalysis::Init : no MC track array!" << endl; }
1281
1282 fStsPointsColl = (TClonesArray*) rootMgr->GetObject("StsPoint");
1283 if (NULL == fStsPointsColl) { cout << "-W- CbmHadronAnalysis::Init : no STS Point array!" << endl; }
1284
1286
1287 InitStatus status = ReadPdfFile();
1288 if (kSUCCESS != status) { return status; }
1289
1290 InitStatus statf = ReadFlowFile();
1291 if (kSUCCESS != statf) { return statf; }
1292
1293 Float_t pbeam = GetBeamMomentum();
1294 {
1295 Float_t ep = TMath::Sqrt(pbeam * pbeam + M2PROT);
1296 Float_t gp = ep / TMath::Sqrt(M2PROT);
1297 Float_t bp = TMath::Sqrt(1. - 1. / gp / gp);
1298 Float_t yp = 0.5 * TMath::Log((1. + bp) / (1. - bp));
1299 SetMidY(yp * 0.5);
1300 cout << "-I- CbmHadronAnalysis: Initialize for pbeam = " << pbeam << " rap:" << GetMidY() << endl;
1301 }
1302
1304
1305 cout << "-I- CbmHadronAnalysis::Init : "
1306 << "initialisation completed." << endl;
1307
1308 return kSUCCESS;
1309}
1310// ------------------------------------------------------------------
1311// ------------------------------------------------------------------
1312void CbmHadronAnalysis::Exec(Option_t* option)
1313{
1314 // Task execution on TS or Events
1315 if (NULL != fEventsColl) { // working for STS and TOF only
1316 iNbTs++;
1317 LOG(debug) << "process TS " << iNbTs << " with " << fEventsColl->GetEntriesFast() << " events";
1318 if (NULL == fTofHits || NULL == fTofHitsColl) { assert("Invalid pointer to Tof TClonesArray"); }
1319 if (NULL == fStsHits || NULL == fStsHitsColl) { assert("Invalid pointer to Sts TClonesArray"); }
1320 LOG(debug) << "TS contains " << fTofHitsColl->GetEntriesFast() << "TOF and " << fStsHitsColl->GetEntriesFast()
1321 << " Sts hits";
1322
1323 for (Int_t iEvent = 0; iEvent < fEventsColl->GetEntriesFast(); iEvent++) {
1324 CbmEvent* tEvent = dynamic_cast<CbmEvent*>(fEventsColl->At(iEvent));
1325 // Inspect CbmEvent
1326 LOG(debug) << "CbmEvent with " << tEvent->GetNofData() << " total data";
1327 LOG(debug) << " " << tEvent->ToString();
1328 // copy TOF hits
1329 fTofHits->Clear();
1330 Int_t iNbHits = 0;
1331 LOG(debug) << "Fill Tof array with mul " << tEvent->GetNofData(ECbmDataType::kTofHit);
1332 for (size_t iHit = 0; iHit < tEvent->GetNofData(ECbmDataType::kTofHit); iHit++) {
1333 Int_t iHitIndex = static_cast<Int_t>(tEvent->GetIndex(ECbmDataType::kTofHit, iHit));
1334 CbmTofHit* pHit = (CbmTofHit*) fTofHitsColl->At(iHitIndex);
1335 new ((*fTofHits)[iNbHits]) CbmTofHit(*pHit); // fill temporary working TClonesArray
1336 iNbHits++;
1337 }
1338 // copy STS hits
1339 iNbHits = 0;
1340 fStsHits->Clear();
1341 LOG(debug) << "Fill Sts array with mul " << tEvent->GetNofData(ECbmDataType::kStsHit) << " out of "
1342 << fStsHitsColl->GetEntriesFast();
1343 for (size_t iHit = 0; iHit < tEvent->GetNofData(ECbmDataType::kStsHit); iHit++) {
1344 Int_t iHitIndex = static_cast<Int_t>(tEvent->GetIndex(ECbmDataType::kStsHit, iHit));
1345 const CbmStsHit* pHit = (CbmStsHit*) fStsHitsColl->At(iHitIndex);
1346 new ((*fStsHits)[iNbHits]) CbmStsHit(*pHit); // fill temporary working TClonesArray
1347 iNbHits++;
1348 }
1349 LOG(debug) << Form("analyze TS %d, Ev %d with %d STS, %d TOF hits", iNbTs, iEvent, fStsHits->GetEntriesFast(),
1350 fTofHits->GetEntriesFast());
1351 ExecEvent(option);
1352 }
1353 }
1354 else { // event based analysis
1355 ExecEvent(option);
1356 }
1357}
1358
1359// ------------------------------------------------------------------
1361{
1362 // Task execution
1363
1364 // Declare variables outside the loop
1365 CbmMCTrack* MCTrack = NULL;
1366 CbmStsTrack* StsTrack = NULL;
1367 CbmTofPoint* TofPoint = NULL;
1368 CbmTofHit* TofHit = NULL;
1369 CbmTofTrack* TofTrack = NULL;
1370 // CbmTofTrack *TofTrackh;
1371 CbmTofTrack* BestTofTrack = NULL;
1372 CbmGlobalTrack* GlobTrack = NULL;
1373 // CbmMatch *tofHitMatch;
1374
1375
1376 Int_t pdgCode, Np1, Np2;
1377 Float_t Qx1, Qy1, Qx2, Qy2, phirp1, phirp2, phirp, delrp, rp_weight;
1378 Float_t RADDEG = 57.29577951;
1379 Float_t p_MC, px_MC, py_MC, pz_MC;
1380 Float_t mfrag = 0.;
1381 Float_t MaxT0 = 0.1;
1382
1383 Int_t TrackP[100000];
1384 Float_t t_hit;
1385 Bool_t use_pions_for_flow = kTRUE;
1386
1387 Int_t verbose = 1;
1388
1389 if (GetBSelMax() > 0.) {
1390 if (fMCEventHeader->GetB() > GetBSelMax()) { return; };
1391 }
1392
1393 if (GetBSelMin() > 0.) {
1394 if (fMCEventHeader->GetB() < GetBSelMin()) { return; };
1395 }
1396
1397 Float_t yrp_mid = GetMidY(); // midrapidity -> update from simulation!
1398 if (fMCTracksColl != NULL) nMCTracks = fMCTracksColl->GetEntriesFast();
1399 if (fTofPoints != NULL) nTofPoints = fTofPoints->GetEntriesFast();
1400 if (fTofHits != NULL) nTofHits = fTofHits->GetEntriesFast();
1401 if (fTofTracks != NULL) nTofTracks = fTofTracks->GetEntriesFast();
1402 if (fGlobalTracks != NULL) nGlobTracks = fGlobalTracks->GetEntriesFast();
1403 if (fStsHits != NULL) nStsHits = fStsHits->GetEntriesFast();
1404
1405 fhTofHitMul->Fill((Double_t) nTofHits);
1406 fhStsHitMul->Fill((Double_t) nStsHits);
1407
1408 if (verbose > 0) { //nh-debug
1409 LOG(debug) << "<D> HadronAnalysis::Exec starting with MCtrks " << nMCTracks << ", TofPoi " << nTofPoints
1410 << ", TofHit " << nTofHits << ", TofTrk " << nTofTracks << ", GlbTrk " << nGlobTracks << ", StsHit "
1411 << nStsHits;
1412 if (fMCEventHeader != NULL)
1413 LOG(debug) << "-D- b = " << fMCEventHeader->GetB() << ", phi = " << fMCEventHeader->GetRotZ();
1414 }
1415 if (fair::Logger::Logging(fair::Severity::debug)) {
1416 for (Int_t j = 0; j < nTofHits; j++) {
1417 TofHit = (CbmTofHit*) fTofHits->At(j);
1418 if (NULL == TofHit) continue;
1419 LOG(debug) << Form("TofHit %d, addr 0x%08x, x %6.1f, y %6.1f, z %6.1f, t %6.1f ", j, TofHit->GetAddress(),
1420 TofHit->GetX(), TofHit->GetY(), TofHit->GetZ(), TofHit->GetTime());
1421 }
1422 }
1423
1424 Double_t dT0 = 0.;
1425 for (Int_t j = 0; j < nTofHits; j++) {
1426 TofHit = (CbmTofHit*) fTofHits->At(j);
1427 if (NULL == TofHit) continue;
1428 if (TofHit->GetZ() == 0.) dT0 = TofHit->GetTime();
1429 }
1430
1431 if (dT0 != 0.) {
1432 for (Int_t j = 0; j < nTofHits; j++) {
1433 TofHit = (CbmTofHit*) fTofHits->At(j);
1434 if (NULL == TofHit) continue;
1435 TofHit->SetTime(TofHit->GetTime() - dT0);
1436 }
1437 }
1438
1439 if (fair::Logger::Logging(fair::Severity::debug)) {
1440 for (Int_t j = 0; j < nTofHits; j++) {
1441 TofHit = (CbmTofHit*) fTofHits->At(j);
1442 if (NULL == TofHit) continue;
1443 LOG(debug) << Form("TofHit %d, addr 0x%08x, x %6.1f, y %6.1f, z %6.1f, t %6.1f ", j, TofHit->GetAddress(),
1444 TofHit->GetX(), TofHit->GetY(), TofHit->GetZ(), TofHit->GetTime());
1445 }
1446 }
1447
1448 if (bRecSec && nTofHits > 1) ReconstructSecondaries(); // independent method
1449
1450 // some local arrays
1451 Int_t MAXNHT = 50;
1452 Int_t NTHMUL[nGlobTracks]; // number of candidate TofTracks
1453 Int_t IndTHMUL[nGlobTracks][MAXNHT]; // ordered array of candidate number of TofTrack indices
1454 Float_t Weight_THMUL[nGlobTracks][MAXNHT]; // weights for ordered array of candidate number of TofTrack indices
1455 // Int_t IndTofTracks[nGlobTracks][MAXNHT][MAXNHT]; // array of TofTrack Indices that selected TofHit is assigned to
1456 Int_t NTofHitTMul[nTofHits]; // number of GlobalTracks assigned to a specific TofHit
1457 Int_t IndTofTrack_TofHit[nTofHits][MAXNHT]; // index of TofTracks assigned to specific TofHit
1458
1459 // generator level
1460 if (NULL != fMCEventHeader) fa_mul_b_gen->Fill(fMCEventHeader->GetB(), nMCTracks);
1461
1462 Qx1 = 0.;
1463 Qy1 = 0.;
1464 Np1 = 0;
1465 Qx2 = 0.;
1466 Qy2 = 0.;
1467 Np2 = 0;
1468 for (Int_t k = 0; k < nMCTracks; k++) { // inspect MCTracks
1469
1470 MCTrack = (CbmMCTrack*) fMCTracksColl->At(k);
1471
1472 if (MCTrack->GetMotherId() != -1) continue; // primary particles only
1473
1474 // Process only hadrons
1475 pdgCode = MCTrack->GetPdgCode();
1476 // cout << " Track k="<<k<<", pdgCode = "<<pdgCode<<endl;
1477
1478 if (pdgCode < 100000000) {
1479 if (211 != TMath::Abs(pdgCode) && // pions
1480 321 != TMath::Abs(pdgCode) && // kaons
1481 2212 != TMath::Abs(pdgCode)) // protons
1482 {
1483 LOG(info) << "skip pdgCode " << pdgCode;
1484 continue;
1485 }
1486 }
1487 else {
1488 mfrag = (pdgCode % 1000) / 10 * .931494028; // ignoring binding energies ...
1489 // where is the proper mass stored ?
1490 }
1491
1492 Float_t Phip = RADDEG * atan2(MCTrack->GetPy(), MCTrack->GetPx());
1493 Float_t dphi = Phip - RADDEG * fMCEventHeader->GetRotZ();
1494 if (dphi < -180.) { dphi += 360.; };
1495 if (dphi > 180.) { dphi -= 360.; };
1496 dphi = dphi / RADDEG;
1497 rp_weight = 0.;
1498
1499 LOG(info) << " Track k=" << k << ", pdgCode = " << pdgCode << " Mass " << MCTrack->GetMass() << "," << mfrag
1500 << " Y " << MCTrack->GetRapidity() << " Pt " << MCTrack->GetPt();
1501
1502 switch (pdgCode) {
1503 case 211: {
1504 fa_ptm_rap_gen_pip->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / MCTrack->GetMass());
1505 fa_v1_rap_gen_pip->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(dphi));
1506 fa_v2_rap_gen_pip->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(2 * dphi));
1507
1508 if (use_pions_for_flow && TMath::Abs((MCTrack->GetRapidity() - yrp_mid) / yrp_mid) > GetDY()
1509 && MCTrack->GetPt() / MCTrack->GetMass() > GetFlowMinPtm()) {
1510 if (MCTrack->GetRapidity() > yrp_mid) { // set weights for reaction plane
1511 rp_weight = -1.;
1512 }
1513 else {
1514 rp_weight = +1.;
1515 }
1516 }
1517 else {
1518 rp_weight = 0.;
1519 }
1520 break;
1521 };
1522 case -211: {
1523 fa_ptm_rap_gen_pim->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / MCTrack->GetMass());
1524 fa_v1_rap_gen_pim->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(dphi));
1525 fa_v2_rap_gen_pim->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(2 * dphi));
1526
1527 if (use_pions_for_flow && TMath::Abs((MCTrack->GetRapidity() - yrp_mid) / yrp_mid) > GetDY()
1528 && MCTrack->GetPt() / MCTrack->GetMass() > GetFlowMinPtm()) {
1529 if (MCTrack->GetRapidity() > yrp_mid) { // set weights for reaction plane
1530 rp_weight = -1.;
1531 }
1532 else {
1533 rp_weight = +1.;
1534 }
1535 }
1536 else {
1537 rp_weight = 0.;
1538 }
1539 break;
1540 };
1541 case 321: {
1542 fa_ptm_rap_gen_kp->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / MCTrack->GetMass());
1543 fa_v1_rap_gen_kp->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(dphi));
1544 fa_v2_rap_gen_kp->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(2 * dphi));
1545 break;
1546 };
1547 case -321: {
1548 fa_ptm_rap_gen_km->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / MCTrack->GetMass());
1549 fa_v1_rap_gen_km->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(dphi));
1550 fa_v2_rap_gen_km->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(2 * dphi));
1551 break;
1552 };
1553 case 2212: { //protons
1554 fa_ptm_rap_gen_p->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / MCTrack->GetMass());
1555 fa_v1_rap_gen_p->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(dphi));
1556 fa_v2_rap_gen_p->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(2 * dphi));
1557
1558 // reaction plane determination
1559 if (TMath::Abs((MCTrack->GetRapidity() - yrp_mid) / yrp_mid) > GetDY()
1560 && MCTrack->GetPt() / MCTrack->GetMass() > GetFlowMinPtm()) {
1561 if (MCTrack->GetRapidity() > yrp_mid) { // set weights for reaction plane
1562 rp_weight = 1.;
1563 }
1564 else {
1565 rp_weight = -1.;
1566 }
1567 }
1568 break;
1569 };
1570 case -2212: {
1571 fa_ptm_rap_gen_pbar->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / MCTrack->GetMass());
1572 fa_v1_rap_gen_pbar->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(dphi));
1573 fa_v2_rap_gen_pbar->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(2 * dphi));
1574 break;
1575 };
1576
1577 case 1000010020: { // deuteron
1578 fa_ptm_rap_gen_d->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / mfrag);
1579 fa_v1_rap_gen_d->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(dphi));
1580 fa_v2_rap_gen_d->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(2 * dphi));
1581 // reaction plane determination
1582 if (TMath::Abs((MCTrack->GetRapidity() - yrp_mid) / yrp_mid) > GetDY()
1583 && MCTrack->GetPt() / MCTrack->GetMass() > GetFlowMinPtm()) {
1584 if (MCTrack->GetRapidity() > yrp_mid) { // set weights for reaction plane
1585 rp_weight = 1.;
1586 }
1587 else {
1588 rp_weight = -1.;
1589 }
1590 }
1591 break;
1592 };
1593
1594 case 1000010030: { // triton
1595 fa_ptm_rap_gen_t->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / mfrag);
1596 fa_v1_rap_gen_t->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(dphi));
1597 fa_v2_rap_gen_t->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(2 * dphi));
1598 // reaction plane determination
1599 if (TMath::Abs((MCTrack->GetRapidity() - yrp_mid) / yrp_mid) > GetDY()
1600 && MCTrack->GetPt() / MCTrack->GetMass() > GetFlowMinPtm()) {
1601 if (MCTrack->GetRapidity() > yrp_mid) { // set weights for reaction plane
1602 rp_weight = 1.;
1603 }
1604 else {
1605 rp_weight = -1.;
1606 }
1607 }
1608 break;
1609 };
1610
1611 case 1000020030: { // 3he
1612 fa_ptm_rap_gen_h->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / mfrag);
1613 fa_v1_rap_gen_h->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(dphi));
1614 fa_v2_rap_gen_h->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(2 * dphi));
1615 // reaction plane determination
1616 if (TMath::Abs((MCTrack->GetRapidity() - yrp_mid) / yrp_mid) > GetDY()
1617 && MCTrack->GetPt() / MCTrack->GetMass() > GetFlowMinPtm()) {
1618 if (MCTrack->GetRapidity() > yrp_mid) { // set weights for reaction plane
1619 rp_weight = 1.;
1620 }
1621 else {
1622 rp_weight = -1.;
1623 }
1624 }
1625 break;
1626 };
1627
1628 case 1000020040: { // alpha
1629 fa_ptm_rap_gen_a->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / mfrag);
1630 fa_v1_rap_gen_a->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(dphi));
1631 fa_v2_rap_gen_a->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(2 * dphi));
1632 // reaction plane determination
1633 if (TMath::Abs((MCTrack->GetRapidity() - yrp_mid) / yrp_mid) > GetDY()
1634 && MCTrack->GetPt() / MCTrack->GetMass() > GetFlowMinPtm()) {
1635 if (MCTrack->GetRapidity() > yrp_mid) { // set weights for reaction plane
1636 rp_weight = 1.;
1637 }
1638 else {
1639 rp_weight = -1.;
1640 }
1641 }
1642 break;
1643 };
1644
1645 default:
1646 { // intermediate mass fragments
1647 //cout << " Track k="<<k<<", pdgCode = "<<pdgCode<<
1648 //" Mass " << MCTrack->GetMass()<<","<<mfrag<<" Y " << MCTrack->GetRapidity() <<
1649 //" Pt " << MCTrack->GetPt() <<endl;
1650 fa_ptm_rap_gen_imf->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / mfrag);
1651 fa_v1_rap_gen_imf->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(dphi));
1652 fa_v2_rap_gen_imf->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(2 * dphi));
1653 // reaction plane determination (optimistic)
1654 if (TMath::Abs((MCTrack->GetRapidity() - yrp_mid) / yrp_mid) > GetDY()
1655 && MCTrack->GetPt() / MCTrack->GetMass() > GetFlowMinPtm()) {
1656 if (MCTrack->GetRapidity() > yrp_mid) { // set weights for reaction plane
1657 rp_weight = 1.;
1658 }
1659 else {
1660 rp_weight = -1.;
1661 }
1662 }
1663 break;
1664 };
1665 }
1666 if (rp_weight != 0.) {
1667 if (gRandom->Uniform(1) > 0.5) { //subdivide events into 2 random subevents
1668 Np1++;
1669 Qx1 = Qx1 + rp_weight * MCTrack->GetPx();
1670 Qy1 = Qy1 + rp_weight * MCTrack->GetPy();
1671 }
1672 else {
1673 Np2++;
1674 Qx2 = Qx2 + rp_weight * MCTrack->GetPx();
1675 Qy2 = Qy2 + rp_weight * MCTrack->GetPy();
1676 }
1677 }
1678 }
1679
1680 if (Np1 > 0 && Np2 > 0) {
1681 phirp1 = atan2(Qy1, Qx1);
1682 phirp2 = atan2(Qy2, Qx2);
1683 if (fflowFile != NULL) { // subevent RP flattening
1684 TH1F* phirp_gen_fpar = fflowFile->Get<TH1F>("phirps_gen_fpar");
1685 Float_t dphir = 0.;
1686 for (int j = 0; j < 4; j++) {
1687 Float_t i = (float) (j + 1);
1688 dphir += (-phirp_gen_fpar->GetBinContent(j) * TMath::Cos(i * phirp1)
1689 + phirp_gen_fpar->GetBinContent(j + 4) * TMath::Sin(i * phirp1))
1690 / i;
1691 }
1692 phirp1 += dphir;
1693
1694 dphir = 0.;
1695 for (int j = 0; j < 4; j++) {
1696 Float_t i = (float) (j + 1);
1697 dphir += (-phirp_gen_fpar->GetBinContent(j) * TMath::Cos(i * phirp2)
1698 + phirp_gen_fpar->GetBinContent(j + 4) * TMath::Sin(i * phirp2))
1699 / i;
1700 }
1701 phirp2 += dphir;
1702 } // subevent RP flattening end
1703
1704 delrp = (phirp1 - phirp2);
1705 fa_phirps_gen->Fill(phirp1 * RADDEG); // 1D histo
1706 fa_phirps_gen->Fill(phirp2 * RADDEG); // 1D histo
1707 if (0) { //nh-debug
1708 cout << "<D-gen> Impact parameter " << fMCEventHeader->GetB() << ", delrp = " << delrp << endl;
1709 }
1710 fa_cdrp_b_gen->Fill(fMCEventHeader->GetB(), TMath::Cos(delrp));
1711 delrp = delrp * RADDEG;
1712 while (delrp < -180.) {
1713 delrp += 360.;
1714 }
1715 while (delrp > 180.) {
1716 delrp -= 360.;
1717 }
1718 fa_drp_b_gen->Fill(fMCEventHeader->GetB(), delrp);
1719
1720 phirp = RADDEG * atan2(Qy1 + Qy2, Qx1 + Qx2); // full reaction plane
1721 while (phirp < -180.) {
1722 phirp += 360.;
1723 }
1724 while (phirp > 180.) {
1725 phirp -= 360.;
1726 }
1727 if (fflowFile != NULL) { // RP flattening
1728 TH1F* phirp_gen_fpar = fflowFile->Get<TH1F>("phirp_gen_fpar");
1729 Float_t dphir = 0.;
1730 for (int j = 0; j < 4; j++) {
1731 Float_t i = (float) (j + 1);
1732 //cout << " RP flat par "<< i << ","<<j<< " par " << phirp_gen_fpar->GetBinContent(j)
1733 // << ","<< phirp_gen_fpar->GetBinContent(j+4) << " phirp "<<phirp<<" dphir "<< dphir << endl;
1734 dphir += ((-phirp_gen_fpar->GetBinContent(j) * TMath::Cos(i * phirp / RADDEG)
1735 + phirp_gen_fpar->GetBinContent(j + 4) * TMath::Sin(i * phirp / RADDEG))
1736 / i);
1737 }
1738 //cout << " phirp " << phirp << " dphir " << dphir*RADDEG << endl;
1739
1740 phirp += dphir * RADDEG;
1741 while (phirp < -180.) {
1742 phirp += 360.;
1743 }
1744 while (phirp > 180.) {
1745 phirp -= 360.;
1746 }
1747 } // RP flattening end
1748 delrp = phirp - RADDEG * fMCEventHeader->GetRotZ();
1749 while (delrp < -180.) {
1750 delrp += 360.;
1751 }
1752 while (delrp > 180.) {
1753 delrp -= 360.;
1754 }
1755
1756 fa_phirp_gen->Fill(phirp);
1757 fa_delrp_b_gen->Fill(fMCEventHeader->GetB(), delrp);
1758 fa_cdelrp_b_gen->Fill(fMCEventHeader->GetB(), TMath::Cos(delrp / RADDEG));
1759 fa_phirp_b_gen->Fill(fMCEventHeader->GetB(), phirp);
1760 fa_phgrp_b_gen->Fill(fMCEventHeader->GetB(), RADDEG * fMCEventHeader->GetRotZ());
1761 fa_phphrp_gen->Fill(phirp, RADDEG * fMCEventHeader->GetRotZ());
1762 } // Np1 && Np2 end
1763
1764 // TofPoint level
1765
1766 for (Int_t k = 0; k < nMCTracks; k++)
1767 TrackP[k] = 0; // reset track detected flag
1768
1769 Qx1 = 0.;
1770 Qy1 = 0.;
1771 Np1 = 0;
1772 Qx2 = 0.;
1773 Qy2 = 0.;
1774 Np2 = 0;
1775 if (NULL != fMCEventHeader) fa_mul_b_poi->Fill(fMCEventHeader->GetB(), nTofPoints);
1776
1777 for (Int_t l = 0; l < nTofPoints; l++) {
1778 TofPoint = (CbmTofPoint*) fTofPoints->At(l);
1779 if (NULL == TofPoint) {
1780 LOG(warning) << " Missing TofPoint at " << l << ", mul " << nTofPoints;
1781 continue;
1782 }
1783 Int_t k = TofPoint->GetTrackID();
1784 MCTrack = (CbmMCTrack*) fMCTracksColl->At(k);
1785 pdgCode = MCTrack->GetPdgCode();
1786 if (pdgCode > 100000000) {
1787 mfrag = (pdgCode % 1000) / 10 * .931494028; // ignoring binding energies ...
1788 }
1789 px_MC = MCTrack->GetPx();
1790 py_MC = MCTrack->GetPy();
1791 pz_MC = MCTrack->GetPz();
1792 p_MC = sqrt(px_MC * px_MC + py_MC * py_MC + pz_MC * pz_MC);
1793
1794 if (TrackP[k] == 0) { // for efficiency
1795 TrackP[k] = 1;
1796
1797 fa_xy_poi1->Fill(TofPoint->GetX(), TofPoint->GetY());
1798 fa_xy_poi2->Fill(TofPoint->GetX(), TofPoint->GetY(), fwxy2);
1799
1800 Float_t vel = TofPoint->GetLength() / TofPoint->GetTime();
1801 Float_t bet = vel / clight;
1802 if (bet > 0.99999) { bet = 0.99999; }
1803 Float_t tofmass = p_MC / bet * TMath::Sqrt(1. - bet * bet) * TMath::Sign(1, pdgCode);
1804
1805 fa_pv_poi->Fill(vel, p_MC);
1806 fa_tm_poi->Fill(p_MC, tofmass);
1807
1808 if (MCTrack->GetMotherId() != -1) continue; // primary particles only
1809 fa_tm_poiprim->Fill(p_MC, tofmass);
1810
1811 Float_t Phip = RADDEG * atan2(MCTrack->GetPy(), MCTrack->GetPx());
1812 Float_t dphi = Phip - RADDEG * fMCEventHeader->GetRotZ();
1813 if (dphi < -180.) { dphi += 360.; };
1814 if (dphi > 180.) { dphi -= 360.; };
1815 dphi = dphi / RADDEG;
1816 rp_weight = 0.;
1817
1818 switch (pdgCode) {
1819 case 211: {
1820 fa_ptm_rap_poi_pip->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / MCTrack->GetMass());
1821 if (use_pions_for_flow && TMath::Abs((MCTrack->GetRapidity() - yrp_mid) / yrp_mid) > GetDY()
1822 && MCTrack->GetPt() / MCTrack->GetMass() > GetFlowMinPtm()) {
1823 if (MCTrack->GetRapidity() > yrp_mid) { // set weights for reaction plane
1824 rp_weight = -1.;
1825 }
1826 else {
1827 rp_weight = +1.;
1828 }
1829 }
1830 else {
1831 rp_weight = 0.;
1832 }
1833 break;
1834 };
1835 case -211: {
1836 fa_ptm_rap_poi_pim->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / MCTrack->GetMass());
1837 if (use_pions_for_flow && TMath::Abs((MCTrack->GetRapidity() - yrp_mid) / yrp_mid) > GetDY()
1838 && MCTrack->GetPt() / MCTrack->GetMass() > GetFlowMinPtm()) {
1839 if (MCTrack->GetRapidity() > yrp_mid) { // set weights for reaction plane
1840 rp_weight = -1.;
1841 }
1842 else {
1843 rp_weight = +1.;
1844 }
1845 }
1846 else {
1847 rp_weight = 0.;
1848 }
1849 break;
1850 };
1851 case 321: {
1852 fa_ptm_rap_poi_kp->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / MCTrack->GetMass());
1853 break;
1854 };
1855 case -321: {
1856 fa_ptm_rap_poi_km->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / MCTrack->GetMass());
1857 break;
1858 };
1859 case 2212: {
1860 fa_ptm_rap_poi_p->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / MCTrack->GetMass());
1861 // reaction plane determination
1862 if (TMath::Abs((MCTrack->GetRapidity() - yrp_mid) / yrp_mid) > GetDY()
1863 && MCTrack->GetPt() / MCTrack->GetMass() > GetFlowMinPtm()) {
1864 if (MCTrack->GetRapidity() > yrp_mid) { // set weights for reaction plane
1865 rp_weight = +1.;
1866 }
1867 else {
1868 rp_weight = -1.;
1869 }
1870 }
1871 break;
1872 };
1873 case -2212: {
1874 fa_ptm_rap_poi_pbar->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / MCTrack->GetMass());
1875 break;
1876 };
1877
1878 case 1000010020: { // deuteron
1879 fa_ptm_rap_poi_d->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / mfrag);
1880 fa_v1_rap_poi_d->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(dphi));
1881 fa_v2_rap_poi_d->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(2 * dphi));
1882 // reaction plane determination
1883 if (TMath::Abs((MCTrack->GetRapidity() - yrp_mid) / yrp_mid) > GetDY()
1884 && MCTrack->GetPt() / MCTrack->GetMass() > GetFlowMinPtm()) {
1885 if (MCTrack->GetRapidity() > yrp_mid) { // set weights for reaction plane
1886 rp_weight = 1.;
1887 }
1888 else {
1889 rp_weight = -1.;
1890 }
1891 }
1892 break;
1893 };
1894
1895 case 1000010030: { // triton
1896 fa_ptm_rap_poi_t->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / mfrag);
1897 fa_v1_rap_poi_t->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(dphi));
1898 fa_v2_rap_poi_t->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(2 * dphi));
1899 // reaction plane determination
1900 if (TMath::Abs((MCTrack->GetRapidity() - yrp_mid) / yrp_mid) > GetDY()
1901 && MCTrack->GetPt() / MCTrack->GetMass() > GetFlowMinPtm()) {
1902 if (MCTrack->GetRapidity() > yrp_mid) { // set weights for reaction plane
1903 rp_weight = 1.;
1904 }
1905 else {
1906 rp_weight = -1.;
1907 }
1908 }
1909 break;
1910 };
1911
1912 case 1000020030: { // 3he
1913 fa_ptm_rap_poi_h->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / mfrag);
1914 fa_v1_rap_poi_h->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(dphi));
1915 fa_v2_rap_poi_h->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(2 * dphi));
1916 // reaction plane determination
1917 if (TMath::Abs((MCTrack->GetRapidity() - yrp_mid) / yrp_mid) > GetDY()
1918 && MCTrack->GetPt() / MCTrack->GetMass() > GetFlowMinPtm()) {
1919 if (MCTrack->GetRapidity() > yrp_mid) { // set weights for reaction plane
1920 rp_weight = 1.;
1921 }
1922 else {
1923 rp_weight = -1.;
1924 }
1925 }
1926 break;
1927 };
1928
1929 case 1000020040: { // alpha
1930 fa_ptm_rap_poi_a->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / mfrag);
1931 fa_v1_rap_poi_a->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(dphi));
1932 fa_v2_rap_poi_a->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(2 * dphi));
1933 // reaction plane determination
1934 if (TMath::Abs((MCTrack->GetRapidity() - yrp_mid) / yrp_mid) > GetDY()
1935 && MCTrack->GetPt() / MCTrack->GetMass() > GetFlowMinPtm()) {
1936 if (MCTrack->GetRapidity() > yrp_mid) { // set weights for reaction plane
1937 rp_weight = 1.;
1938 }
1939 else {
1940 rp_weight = -1.;
1941 }
1942 }
1943 break;
1944 };
1945 default:
1946 { // intermediate mass fragments
1947 //cout << " Track k="<<k<<", pdgCode = "<<pdgCode<<
1948 //" Mass " << MCTrack->GetMass()<<","<<mfrag<<" Y " << MCTrack->GetRapidity() <<
1949 //" Pt " << MCTrack->GetPt() <<endl;
1950 fa_ptm_rap_poi_imf->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / mfrag);
1951 fa_v1_rap_poi_imf->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(dphi));
1952 fa_v2_rap_poi_imf->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(2 * dphi));
1953 // reaction plane determination (optimistic)
1954 if (TMath::Abs((MCTrack->GetRapidity() - yrp_mid) / yrp_mid) > GetDY()
1955 && MCTrack->GetPt() / MCTrack->GetMass() > GetFlowMinPtm()) {
1956 if (MCTrack->GetRapidity() > yrp_mid) { // set weights for reaction plane
1957 rp_weight = 1.;
1958 }
1959 else {
1960 rp_weight = -1.;
1961 }
1962 }
1963 break;
1964 };
1965 }
1966 if (rp_weight != 0.) {
1967 if (gRandom->Uniform(1) > 0.5) { //subdivide events into 2 random subevents
1968 Np1++;
1969 Qx1 = Qx1 + rp_weight * MCTrack->GetPx();
1970 Qy1 = Qy1 + rp_weight * MCTrack->GetPy();
1971 }
1972 else {
1973 Np2++;
1974 Qx2 = Qx2 + rp_weight * MCTrack->GetPx();
1975 Qy2 = Qy2 + rp_weight * MCTrack->GetPy();
1976 }
1977 }
1978 }
1979 }
1980 if (Np1 > 0 && Np2 > 0) {
1981 phirp1 = atan2(Qy1, Qx1);
1982 phirp2 = atan2(Qy2, Qx2);
1983 if (fflowFile != NULL) { // subevent RP flattening
1984 TH1F* phirp_poi_fpar = fflowFile->Get<TH1F>("phirps_poi_fpar");
1985 Float_t dphir = 0.;
1986 for (int j = 0; j < 4; j++) {
1987 Float_t i = (float) (j + 1);
1988 dphir += (-phirp_poi_fpar->GetBinContent(j) * TMath::Cos(i * phirp1)
1989 + phirp_poi_fpar->GetBinContent(j + 4) * TMath::Sin(i * phirp1))
1990 / i;
1991 }
1992 phirp1 += dphir;
1993
1994 dphir = 0.;
1995 for (int j = 0; j < 4; j++) {
1996 Float_t i = (float) (j + 1);
1997 dphir += (-phirp_poi_fpar->GetBinContent(j) * TMath::Cos(i * phirp2)
1998 + phirp_poi_fpar->GetBinContent(j + 4) * TMath::Sin(i * phirp2))
1999 / i;
2000 }
2001 phirp2 += dphir;
2002 } // subevent RP flattening end
2003 delrp = (phirp1 - phirp2);
2004 fa_phirps_poi->Fill(phirp1 * RADDEG); // 1D histo
2005 fa_phirps_poi->Fill(phirp2 * RADDEG); // 1D histo
2006 if (0) { //nh-debug
2007 cout << "<D-poi> Impact parameter " << fMCEventHeader->GetB() << ", delrp = " << delrp << endl;
2008 }
2009 fa_cdrp_b_poi->Fill(fMCEventHeader->GetB(), TMath::Cos(delrp));
2010 delrp = delrp * RADDEG;
2011 if (delrp < -180.) delrp += 360.;
2012 if (delrp > 180.) delrp -= 360.;
2013 fa_drp_b_poi->Fill(fMCEventHeader->GetB(), delrp);
2014 }
2015
2016 phirp = RADDEG * atan2(Qy1 + Qy2, Qx1 + Qx2); // full reaction plane
2017 while (phirp < -180.) {
2018 phirp += 360.;
2019 }
2020 while (phirp > 180.) {
2021 phirp -= 360.;
2022 }
2023 if (fflowFile != NULL) { // RP flattening
2024 TH1F* phirp_poi_fpar = fflowFile->Get<TH1F>("phirp_poi_fpar");
2025 Float_t dphir = 0.;
2026 for (int j = 0; j < 4; j++) {
2027 Float_t i = (float) (j + 1);
2028 //cout << " RP flat par "<< i << ","<<j<< " par " << phirp_gen_fpar->GetBinContent(j)
2029 // << ","<< phirp_gen_fpar->GetBinContent(j+4) << " phirp "<<phirp<<" dphir "<< dphir << endl;
2030 dphir += ((-phirp_poi_fpar->GetBinContent(j) * TMath::Cos(i * phirp / RADDEG)
2031 + phirp_poi_fpar->GetBinContent(j + 4) * TMath::Sin(i * phirp / RADDEG))
2032 / i);
2033 }
2034 //cout << " phirp " << phirp << " dphir " << dphir*RADDEG << endl;
2035
2036 phirp += dphir * RADDEG;
2037 while (phirp < -180.) {
2038 phirp += 360.;
2039 }
2040 while (phirp > 180.) {
2041 phirp -= 360.;
2042 }
2043 } // RP flattening end
2044 if (NULL != fMCEventHeader) {
2045 delrp = phirp - RADDEG * fMCEventHeader->GetRotZ();
2046 while (delrp < -180.) {
2047 delrp += 360.;
2048 }
2049 while (delrp > 180.) {
2050 delrp -= 360.;
2051 }
2052 fa_phirp_poi->Fill(phirp); // 1D histo
2053 fa_delrp_b_poi->Fill(fMCEventHeader->GetB(), delrp);
2054 fa_cdelrp_b_poi->Fill(fMCEventHeader->GetB(), TMath::Cos(delrp / RADDEG));
2055 }
2056 // TofHit level
2057 for (Int_t k = 0; k < nMCTracks; k++)
2058 TrackP[k] = 0; // reset track detected flag
2059
2060 Qx1 = 0.;
2061 Qy1 = 0.;
2062 Np1 = 0;
2063 Qx2 = 0.;
2064 Qy2 = 0.;
2065 Np2 = 0;
2066 Int_t NT0 = 0;
2067 Float_t t0m_hit = 0.;
2068 Int_t NT0F = 0;
2069 Float_t t0mf_hit = 0.;
2070 Int_t NT0NF = 0;
2071 Float_t t0mnf_hit = 0.;
2072 if (NULL != fMCEventHeader) fa_mul_b_hit->Fill(fMCEventHeader->GetB(), nTofHits);
2073
2074 Float_t T0MIN = 0.;
2075 Int_t NFHITS = 10;
2076 Int_t nFTofHits = 0;
2077 Float_t T0FAST[NFHITS];
2078 for (Int_t l = 0; l < NFHITS; l++) {
2079 T0FAST[l] = 100;
2080 };
2081
2082 for (Int_t j = 0; j < nTofHits; j++) {
2083 // reset track assignment vector
2084 NTofHitTMul[j] = 0;
2085 //cout << "<D-hit> j= " << j << endl;
2086 TofHit = (CbmTofHit*) fTofHits->At(j);
2087 if (NULL == TofHit) continue;
2088 t_hit = TofHit->GetTime();
2089 if (1) { // use everything
2090 //if ((TMath::Abs(TofHit->GetX())<55.&&TMath::Abs(TofHit->GetY())<55.)) { // use region E
2091 //if (!(TMath::Abs(TofHit->GetX())<55.&&TMath::Abs(TofHit->GetY())<55.)) { // exclude region E
2092 nFTofHits++;
2093 Float_t dist = TMath::Sqrt(TMath::Power(TofHit->GetX(), 2) + TMath::Power(TofHit->GetY(), 2)
2094 + TMath::Power(TofHit->GetZ(), 2));
2095 Float_t t0_hit = t_hit - dist / clight;
2096 for (Int_t l = 0; l < NFHITS; l++) {
2097 if (t0_hit < T0FAST[l]) {
2098 //cout << "Sort fasted hits: j="<<j<<", l="<<l<<","<<t0_hit<<","<<T0FAST[l]<<endl;
2099 if (T0FAST[l] < 100.) {
2100 for (Int_t ll = NFHITS - 1; ll >= l; ll--) {
2101 T0FAST[ll + 1] = T0FAST[ll];
2102 //cout <<"Move: "<<ll<<" -> "<< ll+1 <<": "<< T0FAST[ll] << endl;
2103 }
2104 }
2105 T0FAST[l] = t0_hit;
2106 //cout<<"Sort:"; for(Int_t ll=0; ll<NFHITS; ll++){ cout << T0FAST[ll] <<",";}cout << endl;
2107 break;
2108 }
2109 else {
2110 }
2111 }
2112 }
2113 }
2114 // cout << "Interpreting result: ";
2115 // T0MIN=average
2116 Int_t nfh = 0;
2117 Int_t nfhmax = nFTofHits;
2118 Int_t nT0It = 0;
2119 Int_t lmin = 0;
2120 Int_t lmax = NFHITS;
2121 Float_t T02 = 0.;
2122 // nfhmax=TMath::Max(nfhmax,3);
2123 for (Int_t l = 0; l < NFHITS; l++) {
2124 if (T0FAST[l] < 100. && nfh < nfhmax) {
2125 lmax = l;
2126 // cout << T0FAST[l] << ", ";
2127 T0MIN += T0FAST[l];
2128 T02 += TMath::Power(T0FAST[l], 2.);
2129 nfh++;
2130 }
2131 }
2132 T0MIN = T0MIN / (Float_t) nfh;
2133 Float_t T0RMS = TMath::Sqrt(T02 - T0MIN * T0MIN);
2134 // cout << " T0MIN = "<<T0MIN << ", T0RMS = "<< T0RMS << ", nfh = "<<nfh
2135 // << " lmin = "<<lmin << ", lmax = "<<lmax <<endl;
2136
2137 while (T0RMS > 1. && nT0It < 10 && nfh > 2) {
2138 nT0It++;
2139 Float_t T0AV = T0MIN;
2140 if (TMath::Abs(T0FAST[lmin] - T0AV) > TMath::Abs(T0FAST[lmax] - T0AV)) {
2141 //remove first entry
2142 lmin++;
2143 }
2144 else { // remove last entry
2145 lmax--;
2146 }
2147 T0MIN = 0.;
2148 nfh = 0;
2149 T02 = 0.;
2150 for (Int_t l = lmin; l < lmax; l++) {
2151 if (T0FAST[l] < 100. && nfh < nfhmax) {
2152 //cout << T0FAST[l] << ", ";
2153 T0MIN += T0FAST[l];
2154 T02 += TMath::Power(T0FAST[l], 2.);
2155 nfh++;
2156 }
2157 }
2158 T0MIN = T0MIN / (Float_t) nfh;
2159 T0RMS = TMath::Sqrt(T02 - T0MIN * T0MIN);
2160 //cout << "Redo("<<nT0It<<"): T0MIN = "<<T0MIN << ", T0RMS = "<< T0RMS << ", nfh = "
2161 // << nfh << " lmin = "<<lmin << ", lmax = "<<lmax <<endl;
2162 }
2163
2164 Int_t lp = -1;
2165 for (Int_t j = 0; j < nTofHits; j++) {
2166 TofHit = (CbmTofHit*) fTofHits->At(j);
2167 if (NULL == TofHit) continue;
2168 if (TofHit->GetTime() < 0.2) continue; //start counter has no poi matches
2169
2170 if (fTofDigiMatchColl != NULL) {
2171 CbmMatch* digiMatch = (CbmMatch*) fTofDigiMatchColl->At(j);
2172 // take first digi's point link
2173 CbmLink L0 = digiMatch->GetLink(0);
2174 Int_t iDigInd0 = L0.GetIndex();
2175 if (fDigiMan->GetMatch(ECbmModuleId::kTof, iDigInd0) != NULL) {
2176 CbmMatch* poiMatch = (CbmMatch*) fDigiMan->GetMatch(ECbmModuleId::kTof, iDigInd0);
2177 if (NULL == poiMatch) {
2178 LOG(warn) << "No MC point found for hit " << j << ", digi " << iDigInd0;
2179 continue;
2180 }
2181 LOG(debug) << "Got PoiMatch for TofHit " << j << ", digi " << iDigInd0 << ": " << poiMatch;
2182 CbmLink LP;
2183 try {
2184 LP = poiMatch->GetMatchedLink();
2185 }
2186 catch (...) {
2187 LOG(info) << "Got invalid PoiMatch for TofHit " << j << ", digi " << iDigInd0 << ": " << poiMatch;
2188 continue;
2189 }
2190 lp = LP.GetIndex();
2191 }
2192 }
2193 else {
2194 lp = -1;
2195 LOG(debug) << "No Link to MCTofPoint found for hit " << j;
2196 continue;
2197 }
2198
2199 if (NULL != fTofPoints) {
2200 TofPoint = (CbmTofPoint*) fTofPoints->At(lp);
2201 Int_t k = TofPoint->GetTrackID();
2202 MCTrack = (CbmMCTrack*) fMCTracksColl->At(k);
2203 pdgCode = MCTrack->GetPdgCode();
2204 px_MC = MCTrack->GetPx();
2205 py_MC = MCTrack->GetPy();
2206 pz_MC = MCTrack->GetPz();
2207 p_MC = sqrt(px_MC * px_MC + py_MC * py_MC + pz_MC * pz_MC);
2208 if (k > 100000) {
2209 cout << "<E-hit> Too many MCTracks " << k << " from " << nMCTracks << endl;
2210 continue;
2211 }
2212 LOG(debug) << Form("HadronAnalysis:: hit %d, poi %d, MCt %d, Eff %d, mom %7.2f, pdg %d ", j, lp, k, TrackP[k],
2213 p_MC, pdgCode);
2214 /*
2215 if(TofHit->GetX()<-90.) { //nh-debug
2216 LOG(info) << Form(" Invalid Hit in ev %d: %6.1f, %6.1f, %6.1f, poi: %6.1f, %6.1f, %6.1f, pdg: %d",
2217 fEvents,TofHit->GetX(),TofHit->GetY(),TofHit->GetZ(),
2218 TofPoint->GetX(),TofPoint->GetY(),TofPoint->GetZ(),
2219 pdgCode);
2220 LOG(fatal) << "<D-hit> j " << j << ", l " << l << ", k " << k << ", lp ";
2221 }
2222 */
2223 if (TrackP[k] == 0) { // for efficiency
2224 // if (1) {
2225 TrackP[k]++;
2226 t_hit = TofHit->GetTime();
2227 Float_t delta_tof = TofPoint->GetTime() - t_hit;
2228 Float_t delta_x = TofPoint->GetX() - TofHit->GetX();
2229 Float_t delta_y = TofPoint->GetY() - TofHit->GetY();
2230 Float_t delta_z = TofPoint->GetZ() - TofHit->GetZ();
2231
2232 fa_dxx->Fill(TofPoint->GetX(), delta_x);
2233 fa_dxy->Fill(TofPoint->GetY(), delta_x);
2234 fa_dxz->Fill(TofPoint->GetZ(), delta_x);
2235 fa_dyx->Fill(TofPoint->GetX(), delta_y);
2236 fa_dyy->Fill(TofPoint->GetY(), delta_y);
2237 fa_dyz->Fill(TofPoint->GetZ(), delta_y);
2238 fa_dzx->Fill(TofPoint->GetX(), delta_z);
2239 fa_dzy->Fill(TofPoint->GetY(), delta_z);
2240 fa_dzz->Fill(TofPoint->GetZ(), delta_z);
2241
2242 fa_xy_hit1->Fill(TofHit->GetX(), TofHit->GetY());
2243 fa_xy_hit2->Fill(TofHit->GetX(), TofHit->GetY(), fwxy2);
2244 fa_hit_ch->Fill(TofHit->GetCh());
2245 fa_dhit_ch->Fill(TofHit->GetCh(), TofHit->GetFlag());
2246
2247 Float_t vel = TofPoint->GetLength() / t_hit;
2248 Float_t bet = vel / clight;
2249 if (bet > 0.99999) { bet = 0.99999; }
2250 Float_t tofmass = p_MC / bet * TMath::Sqrt(1. - bet * bet) * TMath::Sign(1, pdgCode);
2251 Float_t dist = TMath::Sqrt(TMath::Power(TofHit->GetX(), 2) + TMath::Power(TofHit->GetY(), 2)
2252 + TMath::Power(TofHit->GetZ(), 2));
2253 fa_tof_hit->Fill(t_hit);
2254 fa_pv_hit->Fill(vel, p_MC);
2255 fa_tm_hit->Fill(p_MC, tofmass);
2256 fa_dtof_hit->Fill(delta_tof);
2257 Float_t t0_hit = t_hit - dist / clight;
2258 fa_t0_hit->Fill(t0_hit);
2259 //if(t0_hit<MaxT0) {
2260 if (t0_hit < T0MIN + 2.4 * MaxT0) {
2261 NT0++;
2262 t0m_hit = ((Float_t)(NT0 - 1) * t0m_hit + t0_hit) / (Float_t) NT0;
2263 if ((TMath::Abs(TofHit->GetX()) < 55. && TMath::Abs(TofHit->GetY()) < 55.)) { // region E
2264 NT0F++;
2265 t0mf_hit = ((Float_t)(NT0F - 1) * t0mf_hit + t0_hit) / (Float_t) NT0F;
2266 }
2267 else {
2268 NT0NF++;
2269 t0mnf_hit = ((Float_t)(NT0NF - 1) * t0mnf_hit + t0_hit) / (Float_t) NT0NF;
2270 }
2271 }
2272
2273 if (MCTrack->GetMotherId() != -1) continue; // primary particles only
2274 if (TrackP[k] > 1) continue; // for efficiency consider only first hit of track
2275 fa_tm_hitprim->Fill(p_MC, tofmass);
2276 fa_tof_hitprim->Fill(t_hit);
2277
2278 Float_t Phip = RADDEG * atan2(MCTrack->GetPy(), MCTrack->GetPx());
2279 Float_t dphi = Phip - RADDEG * fMCEventHeader->GetRotZ();
2280 if (dphi < -180.) { dphi += 360.; };
2281 if (dphi > 180.) { dphi -= 360.; };
2282 dphi = dphi / RADDEG;
2283
2284 rp_weight = 0.;
2285
2286 switch (pdgCode) {
2287 case 211: {
2288 fa_ptm_rap_hit_pip->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / MCTrack->GetMass());
2289 fa_v1_rap_hit_pip->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(dphi));
2290 fa_v2_rap_hit_pip->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(2 * dphi));
2291
2292 if (use_pions_for_flow && TMath::Abs((MCTrack->GetRapidity() - yrp_mid) / yrp_mid) > GetDY()
2293 && MCTrack->GetPt() / MCTrack->GetMass() > GetFlowMinPtm()) {
2294 if (MCTrack->GetRapidity() > yrp_mid) { // set weights for reaction plane
2295 rp_weight = -1.;
2296 }
2297 else {
2298 rp_weight = +1.;
2299 }
2300 }
2301 else {
2302 rp_weight = 0.;
2303 }
2304 break;
2305 };
2306 case -211: {
2307 fa_ptm_rap_hit_pim->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / MCTrack->GetMass());
2308 fa_v1_rap_hit_pim->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(dphi));
2309 fa_v2_rap_hit_pim->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(2 * dphi));
2310
2311 if (use_pions_for_flow && TMath::Abs((MCTrack->GetRapidity() - yrp_mid) / yrp_mid) > GetDY()
2312 && MCTrack->GetPt() / MCTrack->GetMass() > GetFlowMinPtm()) {
2313 if (MCTrack->GetRapidity() > yrp_mid) { // set weights for reaction plane
2314 rp_weight = -1.;
2315 }
2316 else {
2317 rp_weight = +1.;
2318 }
2319 }
2320 else {
2321 rp_weight = 0.;
2322 }
2323 break;
2324 };
2325 case 321: {
2326 fa_ptm_rap_hit_kp->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / MCTrack->GetMass());
2327 fa_v1_rap_hit_kp->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(dphi));
2328 fa_v2_rap_hit_kp->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(2 * dphi));
2329 break;
2330 };
2331 case -321: {
2332 fa_ptm_rap_hit_km->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / MCTrack->GetMass());
2333 fa_v1_rap_hit_km->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(dphi));
2334 fa_v2_rap_hit_km->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(2 * dphi));
2335 break;
2336 };
2337 case 2212: {
2338 fa_ptm_rap_hit_p->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / MCTrack->GetMass());
2339 fa_v1_rap_hit_p->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(dphi));
2340 fa_v2_rap_hit_p->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(2 * dphi));
2341 // reaction plane determination
2342 if (TMath::Abs((MCTrack->GetRapidity() - yrp_mid) / yrp_mid) > GetDY()
2343 && MCTrack->GetPt() / MCTrack->GetMass() > GetFlowMinPtm()) {
2344 if (MCTrack->GetRapidity() > yrp_mid) { // set weights for reaction plane
2345 rp_weight = +1.;
2346 }
2347 else {
2348 rp_weight = -1.;
2349 }
2350 }
2351 else {
2352 rp_weight = 0.;
2353 }
2354 break;
2355 };
2356 case -2212: {
2357 fa_ptm_rap_hit_pbar->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / MCTrack->GetMass());
2358 fa_v1_rap_hit_pbar->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(dphi));
2359 fa_v2_rap_hit_pbar->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(2 * dphi));
2360 break;
2361 };
2362
2363 case 1000010020: { // deuteron
2364 fa_ptm_rap_hit_d->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / MCTrack->GetMass());
2365 fa_v1_rap_hit_d->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(dphi));
2366 fa_v2_rap_hit_d->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(2 * dphi));
2367 // reaction plane determination
2368 if (TMath::Abs((MCTrack->GetRapidity() - yrp_mid) / yrp_mid) > GetDY()
2369 && MCTrack->GetPt() / MCTrack->GetMass() > GetFlowMinPtm()) {
2370 if (MCTrack->GetRapidity() > yrp_mid) { // set weights for reaction plane
2371 rp_weight = 1.;
2372 }
2373 else {
2374 rp_weight = -1.;
2375 }
2376 }
2377 break;
2378 };
2379
2380 case 1000010030: { // triton
2381 fa_ptm_rap_hit_t->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / MCTrack->GetMass());
2382 fa_v1_rap_hit_t->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(dphi));
2383 fa_v2_rap_hit_t->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(2 * dphi));
2384 // reaction plane determination
2385 if (TMath::Abs((MCTrack->GetRapidity() - yrp_mid) / yrp_mid) > GetDY()
2386 && MCTrack->GetPt() / MCTrack->GetMass() > GetFlowMinPtm()) {
2387 if (MCTrack->GetRapidity() > yrp_mid) { // set weights for reaction plane
2388 rp_weight = 1.;
2389 }
2390 else {
2391 rp_weight = -1.;
2392 }
2393 }
2394 break;
2395 };
2396
2397 case 1000020030: { // 3he
2398 fa_ptm_rap_hit_h->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / MCTrack->GetMass());
2399 fa_v1_rap_hit_h->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(dphi));
2400 fa_v2_rap_hit_h->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(2 * dphi));
2401 // reaction plane determination
2402 if (TMath::Abs((MCTrack->GetRapidity() - yrp_mid) / yrp_mid) > GetDY()
2403 && MCTrack->GetPt() / MCTrack->GetMass() > GetFlowMinPtm()) {
2404 if (MCTrack->GetRapidity() > yrp_mid) { // set weights for reaction plane
2405 rp_weight = 1.;
2406 }
2407 else {
2408 rp_weight = -1.;
2409 }
2410 }
2411 break;
2412 };
2413 case 1000020040: { // alpha
2414 fa_ptm_rap_hit_a->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / MCTrack->GetMass());
2415 fa_v1_rap_hit_a->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(dphi));
2416 fa_v2_rap_hit_a->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(2 * dphi));
2417 // reaction plane determination
2418 if (TMath::Abs((MCTrack->GetRapidity() - yrp_mid) / yrp_mid) > GetDY()
2419 && MCTrack->GetPt() / MCTrack->GetMass() > GetFlowMinPtm()) {
2420 if (MCTrack->GetRapidity() > yrp_mid) { // set weights for reaction plane
2421 rp_weight = 1.;
2422 }
2423 else {
2424 rp_weight = -1.;
2425 }
2426 }
2427 break;
2428 };
2429
2430 default: { // intermediate mass fragments
2431 //cout << " Track k="<<k<<", pdgCode = "<<pdgCode<<
2432 //" Mass " << MCTrack->GetMass()<<","<<MCTrack->GetMass()<<" Y " << MCTrack->GetRapidity() <<
2433 //" Pt " << MCTrack->GetPt() <<endl;
2434 fa_ptm_rap_hit_imf->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / MCTrack->GetMass());
2435 fa_v1_rap_hit_imf->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(dphi));
2436 fa_v2_rap_hit_imf->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(2 * dphi));
2437 // reaction plane determination (optimistic)
2438 if (TMath::Abs((MCTrack->GetRapidity() - yrp_mid) / yrp_mid) > GetDY()
2439 && MCTrack->GetPt() / MCTrack->GetMass() > GetFlowMinPtm()) {
2440 if (MCTrack->GetRapidity() > yrp_mid) { // set weights for reaction plane
2441 rp_weight = 1.;
2442 }
2443 else {
2444 rp_weight = -1.;
2445 }
2446 }
2447 break;
2448 };
2449 }
2450
2451 if (rp_weight != 0.) {
2452 if (gRandom->Uniform(1) > 0.5) { //subdivide events into 2 random subevents
2453 Np1++;
2454 Qx1 = Qx1 + rp_weight * MCTrack->GetPx();
2455 Qy1 = Qy1 + rp_weight * MCTrack->GetPy();
2456 }
2457 else {
2458 Np2++;
2459 Qx2 = Qx2 + rp_weight * MCTrack->GetPx();
2460 Qy2 = Qy2 + rp_weight * MCTrack->GetPy();
2461 }
2462 }
2463 }
2464 else {
2465 if (0) {
2466 cout << "<W> CHA: >=2. hit from track k =" << k << " Hit# "
2467 << j
2468 // << " Str(Cell) "<< TofHit->GetCell() << "," << TofPoint->GetCell()
2469 // << " Module " << TofHit->GetModule() << "," << TofPoint->GetModule()
2470 // << " SM " << TofHit->GetSModule() << "," << TofPoint->GetSModule()
2471 // << " SMType " << TofHit->GetSMtype() << "," << TofPoint->GetSMtype()
2472 << endl;
2473 }
2474 }
2475 }
2476 if (Np1 > 0 && Np2 > 0) {
2477 phirp1 = atan2(Qy1, Qx1);
2478 phirp2 = atan2(Qy2, Qx2);
2479 if (fflowFile != NULL) { // subevent RP flattening
2480 TH1F* phirp_hit_fpar = fflowFile->Get<TH1F>("phirps_hit_fpar");
2481 Float_t dphir = 0.;
2482 for (int jj = 0; jj < 4; jj++) {
2483 Float_t i = (float) (jj + 1);
2484 dphir += (-phirp_hit_fpar->GetBinContent(jj) * TMath::Cos(i * phirp1)
2485 + phirp_hit_fpar->GetBinContent(jj + 4) * TMath::Sin(i * phirp1))
2486 / i;
2487 }
2488 phirp1 += dphir;
2489
2490 dphir = 0.;
2491 for (int jj = 0; jj < 4; jj++) {
2492 Float_t i = (float) (jj + 1);
2493 dphir += (-phirp_hit_fpar->GetBinContent(jj) * TMath::Cos(i * phirp2)
2494 + phirp_hit_fpar->GetBinContent(jj + 4) * TMath::Sin(i * phirp2))
2495 / i;
2496 }
2497 phirp2 += dphir;
2498 } // subevent RP flattening end
2499 fa_phirps_hit->Fill(phirp1 * RADDEG); // 1D histo
2500 fa_phirps_hit->Fill(phirp2 * RADDEG); // 1D histo
2501 delrp = (phirp1 - phirp2);
2502 if (NULL != fMCEventHeader) {
2503 if (0) { //nh-debug
2504 cout << "<D-hit> Impact parameter " << fMCEventHeader->GetB() << ", delrp = " << delrp << endl;
2505 }
2506 fa_cdrp_b_hit->Fill(fMCEventHeader->GetB(), TMath::Cos(delrp));
2507 delrp = delrp * RADDEG;
2508 if (delrp < -180.) delrp += 360.;
2509 if (delrp > 180.) delrp -= 360.;
2510 fa_drp_b_hit->Fill(fMCEventHeader->GetB(), delrp);
2511 }
2512 }
2513 phirp = RADDEG * atan2(Qy1 + Qy2, Qx1 + Qx2); // full reaction plane
2514 while (phirp < -180.) {
2515 phirp += 360.;
2516 }
2517 while (phirp > 180.) {
2518 phirp -= 360.;
2519 }
2520 if (fflowFile != NULL) { // RP flattening
2521 TH1F* phirp_hit_fpar = fflowFile->Get<TH1F>("phirps_hit_fpar");
2522 Float_t dphir = 0.;
2523 for (int jj = 0; jj < 4; jj++) {
2524 Float_t i = (float) (jj + 1);
2525 //cout << " RP flat par "<< i << ","<<jj<< " par " << phirp_gen_fpar->GetBinContent(j)
2526 // << ","<< phirp_gen_fpar->GetBinContent(j+4) << " phirp "<<phirp<<" dphir "<< dphir << endl;
2527 dphir += ((-phirp_hit_fpar->GetBinContent(jj) * TMath::Cos(i * phirp / RADDEG)
2528 + phirp_hit_fpar->GetBinContent(jj + 4) * TMath::Sin(i * phirp / RADDEG))
2529 / i);
2530 }
2531 //cout << " phirp " << phirp << " dphir " << dphir*RADDEG << endl;
2532
2533 phirp += dphir * RADDEG;
2534 while (phirp < -180.) {
2535 phirp += 360.;
2536 }
2537 while (phirp > 180.) {
2538 phirp -= 360.;
2539 }
2540 } // RP flattening end
2541 if (NULL != fMCEventHeader) {
2542 delrp = phirp - RADDEG * fMCEventHeader->GetRotZ();
2543 while (delrp < -180.) {
2544 delrp += 360.;
2545 }
2546 while (delrp > 180.) {
2547 delrp -= 360.;
2548 }
2549 fa_phirp_hit->Fill(phirp); // 1D histo
2550 fa_delrp_b_hit->Fill(fMCEventHeader->GetB(), delrp);
2551 fa_cdelrp_b_hit->Fill(fMCEventHeader->GetB(), TMath::Cos(delrp / RADDEG));
2552
2553 fa_tn_hit->Fill(T0MIN);
2554 fa_t0mn_hit->Fill((Float_t) NT0);
2555 fa_t0mn_b_hit->Fill(fMCEventHeader->GetB(), (Float_t) NT0);
2556 if (NT0 > 0) {
2557 fa_t0m_hit->Fill(t0m_hit);
2558 fa_t0m_b_hit->Fill(fMCEventHeader->GetB(), t0m_hit);
2559 }
2560
2561 fa_t0mn_f_hit->Fill((Float_t) NT0F);
2562 fa_t0mn_f_b_hit->Fill(fMCEventHeader->GetB(), (Float_t) NT0F);
2563 if (NT0F > 0) {
2564 fa_t0m_f_hit->Fill(t0mf_hit);
2565 fa_t0m_f_b_hit->Fill(fMCEventHeader->GetB(), t0mf_hit);
2566 }
2567 fa_t0mn_nf_hit->Fill((Float_t) NT0NF);
2568 fa_t0mn_nf_b_hit->Fill(fMCEventHeader->GetB(), (Float_t) NT0NF);
2569 if (NT0NF > 0) {
2570 fa_t0m_nf_hit->Fill(t0mnf_hit);
2571 fa_t0m_nf_b_hit->Fill(fMCEventHeader->GetB(), t0mnf_hit);
2572 }
2573 }
2574 // cout << "<I> CbmHadronAnalysis: Number of T0 particles: NTO = " << NT0 << endl;
2576 } // TofPoints check end
2577
2578 // GlobalTrack level analysis
2579
2580 const Double_t DISTMAX = 100.;
2581 Double_t WMAX = 10.; //0.02;
2582 if (fWMax != 0.) WMAX = fWMax;
2583
2584 Int_t NGTofTrack = 0;
2585 Qx1 = 0.;
2586 Qy1 = 0.;
2587 Np1 = 0;
2588 Qx2 = 0.;
2589 Qy2 = 0.;
2590 Np2 = 0;
2591 if (NULL != fMCEventHeader) fa_mul_b_glo->Fill(fMCEventHeader->GetB(), nGlobTracks);
2592
2593 Int_t NReas = 0; //100; // activate reassignment of hits to global tracks
2594 Int_t NRIt = 0;
2595 while (NReas > 0) {
2596 NRIt++;
2597 if (verbose > 5) { // nh-debug readability
2598 cout << endl;
2599 cout << Form(" TofTrack selection: %d. iteration, reassigned hits: %d, "
2600 "GlobTracks %d, TofTracks %d",
2601 NRIt, NReas, nGlobTracks, nTofTracks)
2602 << endl;
2603 }
2604 NReas = 0;
2605 for (Int_t i = 0; i < nGlobTracks; i++) { // inspect global tracks
2606 GlobTrack = (CbmGlobalTrack*) fGlobalTracks->At(i);
2607 if (NRIt == 1) NTHMUL[i] = 0; // number of TofTracks assigned to GlobTrack
2608 Int_t Btt = -1; // best TofTrack index
2609 Int_t Bthi = -1; // best TofHit index
2610
2611 Int_t s = GlobTrack->GetStsTrackIndex();
2612 Int_t j = GlobTrack->GetTofHitIndex();
2613
2614 if (verbose > 10) { // nh-debug
2615 cout << "<Di> NRIt " << NRIt << ": Global Track " << i << ", TofHit " << j << " StsTrk " << s << endl;
2616 }
2617
2618 const FairTrackParam* tparf = GlobTrack->GetParamFirst();
2619 if (0 == tparf->GetQp()) {
2620 if (verbose > 10) { // nh-debug
2621 cout << "<W> Global Track " << i << " without Qp!, take from Sts " << s << endl;
2622 }
2623 StsTrack = (CbmStsTrack*) fStsTracks->At(s);
2624 if (NULL == StsTrack) {
2625 cout << "<E> Invalid StsTrack pointer at location " << s << endl;
2626 continue;
2627 }
2628 GlobTrack->SetParamFirst(StsTrack->GetParamFirst());
2629 }
2630 tparf = GlobTrack->GetParamFirst();
2631 if (0 == tparf->GetQp()) {
2632 cout << "<E> Global Track " << i << " without Qp! " << endl;
2633 continue;
2634 }
2635
2636 if (211 != StsTrack->GetPidHypo()) {
2637 cout << "<E> Invalid StsTrack PID " << StsTrack->GetPidHypo() << " at " << s << endl;
2638 continue;
2639 }
2640
2641 FairTrackParam paramExtr;
2642 fTrackFitter.FitToVertex(StsTrack, fPrimVertex, &paramExtr);
2643 // GlobTrack -> SetParamFirst(&paramExtr); // nh: attach track parameter to global track
2644 // cout <<Form("<D> Extrapolate StsTrack %d with PidHypo %6.0f to vertex (%5.2f,%5.2f,%5.2f)",
2645 // s,StsTrack->GetPidHypo(),fPrimVertex->GetX(),fPrimVertex->GetY(),fPrimVertex->GetZ()) <<endl;
2646 Double_t vtxb = fTrackFitter.GetChiToVertex(StsTrack, fPrimVertex); //impact paramter ???
2647 if (verbose > 10) { // nh-debug
2648 cout << Form("<D> Extrapolate Glob Track %d to prim. vertex %6.2f with "
2649 "chi2 %6.2f",
2650 i, fPrimVertex->GetZ(), vtxb)
2651 << endl;
2652 //GlobTrack->GetParamFirst()->Print();
2653 }
2654
2655
2656 Float_t momf = 1. / tparf->GetQp();
2657 if (momf < 0.) momf = -momf; // positive momentum at vertex
2658
2659 Float_t dist = 0.;
2660
2661 if (j > -1) { // TOF Track analysis
2662
2663 if (NRIt == 1) {
2664 Float_t DISTMIN = WMAX;
2665 Int_t nth = -1; // number of TofHits for global track i
2666 Weight_THMUL[i][0] = WMAX; // initialize limit
2667 for (Int_t tt = 0; tt < nTofTracks; tt++) { // loop over all TofTracks
2668 TofTrack = (CbmTofTrack*) fTofTracks->At(tt);
2669 if (i == TofTrack->GetTrackIndex()) { // select TofTrack belonging to global track i
2670 if (verbose > 10) { // nh-debug
2671 cout << "<Dt> Global Track " << i << ", TofHit " << j << ", StsTrk " << s << ", TofTrack " << tt
2672 << endl;
2673 }
2674
2675 Int_t thi = TofTrack->GetTofHitIndex();
2676 TofHit = (CbmTofHit*) fTofHits->At(thi);
2677 if (NULL == TofHit) continue;
2678
2679 dist = TofTrack->GetDistance();
2680
2681 if (isinf(dist)) {
2682 cout << "<E> invalid dist for gt " << i << ", tt " << tt << ", d:" << dist << endl;
2683 break;
2684 }
2685 /*
2686 fhTofTrkDx->Fill(TofTrack->GetTrackDx());
2687 fhTofTrkDy->Fill(TofTrack->GetTrackDy());
2688 */
2689
2690 //dist=TMath::Abs(TMath::Abs(dist)-0.5);
2691
2692 const FairTrackParam* tpar = TofTrack->GetTrackParameter();
2693 //cout << " Inspect TrackParameter "; tpar->Print();
2694 Float_t moml = 1. / tpar->GetQp();
2695 if (moml < 0.) moml = -moml;
2696 Float_t bet = TofHit->GetR() / TofHit->GetTime() / clight;
2697 //nh-inconsistent (?), TrackLength needs to be determined experimentally
2698 if (bet > 0.99999) { bet = 0.99999; }
2699 Float_t tofmass = momf / bet * sqrt(1. - bet * bet) * TMath::Sign(1., tpar->GetQp());
2700 TofTrack->SetMass(tofmass); // store tofmass as part of PID info
2701 if (TofTrack->GetMass() != tofmass) { cout << "<E> did not store tofmass properly " << tofmass << endl; }
2702 // calculate attachment weight, to be refined ... (nh, 03.01.2014)
2703 Float_t mindelmass = 1.E6;
2704 Float_t minlen;
2705 Int_t immin; // minlen = (TofTrack->GetTrackLength()-MinWallDist)/MinWallDist;
2706 minlen = (TofHit->GetR() - MinWallDist) / MinWallDist;
2707 mindelmass = 1.E6;
2708 immin = 0;
2709 for (Int_t im = 0; im < NMASS; im++) {
2710 Float_t delmass = TMath::Abs(TMath::Abs(TofTrack->GetMass()) - refMass[im]);
2711 if (delmass < mindelmass) {
2712 mindelmass = delmass;
2713 immin = im;
2714 }
2715 }
2716 mindelmass /= refMass[immin];
2717 Float_t delp = TMath::Abs((momf - moml) / momf);
2718 // Float_t w = dist * mindelmass * minlen * delp * 10000.;
2719 Float_t w = dist;
2720 //Float_t w = minlen;
2721 //Float_t w = mindelmass ;
2722 //Float_t w = delp;
2723 //Float_t w = TMath::Sqrt((dist*dist + delp*delp)/2.);
2724 //Float_t w = TMath::Sqrt((dist*dist + mindelmass*mindelmass)/2.);
2725 //Float_t w = TMath::Sqrt((dist*dist + delp*delp + mindelmass*mindelmass + minlen*minlen)/4.);
2726 //Float_t w = dist * mindelmass * minlen ;
2727 //Float_t w = dist * mindelmass * 15.;
2728 //Float_t w = dist * mindelmass *delp;
2729 //Float_t w = dist * delp;
2730
2731 fhwdist->Fill(momf, dist);
2732 fhwmindelmass->Fill(momf, mindelmass);
2733 fhwminlen->Fill(momf, minlen);
2734 fhwdelp->Fill(momf, delp);
2735 if (verbose > 5) {
2736 cout << Form("<D> w for gt %3d, tt %3d, w: %9.5f, d: %7.2f, m: "
2737 "%7.3f, l: %7.2f, dp: %7.3f, p: %7.2f ",
2738 i, tt, w, dist, mindelmass, minlen, delp, momf)
2739 << endl;
2740 }
2741 if (w < WMAX) {
2742 nth++;
2743 if (nth == MAXNHT) {
2744 if (verbose > 1) {
2745 cout << "<W> Too many TofTrack candidates for track " << i << ", limit to " << nth << endl;
2746 }
2747 nth = MAXNHT - 1;
2748 }
2749 // sort TofTracks according to weight into array IndTHMUL
2750 Int_t jthmin = nth;
2751 for (Int_t jth = 0; jth < nth; jth++) { //determine position in array
2752 if (verbose > 10) {
2753 cout << " Compare for position " << jth << " w " << w << " - " << Weight_THMUL[i][jth] << endl;
2754 }
2755 if (w < Weight_THMUL[i][jth]) {
2756 jthmin = jth;
2757 break;
2758 }
2759 }
2760 if (verbose > 10) {
2761 cout << " Put Track " << tt << " with w = " << w << " to position " << jthmin << " of " << nth
2762 << endl;
2763 }
2764 for (Int_t jth = nth; jth > jthmin; jth--) { // save old entries
2765 if (verbose > 10) {
2766 cout << " Save Track " << IndTHMUL[i][jth - 1] << " with w " << Weight_THMUL[i][jth - 1]
2767 << " to position " << jth << endl;
2768 }
2769 IndTHMUL[i][jth] = IndTHMUL[i][jth - 1];
2770 Weight_THMUL[i][jth] = Weight_THMUL[i][jth - 1];
2771 }
2772 IndTHMUL[i][jthmin] = tt; // store index of TofTrack
2773 Weight_THMUL[i][jthmin] = w; // store weight of TofTrack
2774
2775 if (w < DISTMIN) {
2776 DISTMIN = w;
2777 BestTofTrack = TofTrack;
2778 Bthi = thi; // best TofHit index
2779 Btt = tt; // best TofTrack index
2780 if (verbose > 5) { cout << Form("<DMin> gt %d, hit %d, tt %d, w: %6.2f", i, Bthi, Btt, w) << endl; }
2781 }
2782 } //w < WMAX end
2783 if (verbose > 10) { cout << Form("<D> tt-loop: gt %d, tt %d, w: %6.2f", i, tt, w) << endl; }
2784 } // (i==TofTrack->GetTrackIndex())
2785 } // inspection of all TofTracks finished
2786 NTHMUL[i] = nth + 1; // number of TofHit candidates
2787 fa_TofTrackMul->Fill(NTHMUL[i]);
2788
2789 // report summary:
2790 if (verbose > 5) {
2791 for (Int_t k = 0; k < NTHMUL[i]; k++) {
2792 if (verbose > 3) {
2793 cout << Form("<Ddeb> i %d, k %d, M %d, Ind %d ", i, k, NTHMUL[i], IndTHMUL[i][k]) << endl;
2794 }
2795
2796 TofTrack = (CbmTofTrack*) fTofTracks->At(IndTHMUL[i][k]);
2797 cout << "<DSum> GlobTrack " << i << ", TMul: " << NTHMUL[i] << ", TofTrack " << IndTHMUL[i][k]
2798 << ", TofHit " << TofTrack->GetTofHitIndex() << ", TMul_hit "
2799 << NTofHitTMul[TofTrack->GetTofHitIndex()] << ", dist " << TofTrack->GetDistance() << ", len "
2800 << TofTrack->GetTrackLength() << ", R " << ((CbmTofHit*) fTofHits->At(j))->GetR() << ", mass "
2801 << TofTrack->GetMass() << ", mom " << momf << ", w " << Weight_THMUL[i][k] << endl;
2802 }
2803 }
2804 }
2805 else { // NRIt>1; initialize from array
2806 if (NTHMUL[i] > 0) {
2807 Btt = IndTHMUL[i][0];
2808 if (Btt < 0 || Btt > fTofTracks->GetEntriesFast()) {
2809 cout << "<E> invalid TofTrackIndex " << Btt << ", gt " << i << ", NRIt " << NRIt << endl;
2810 Btt = -1;
2811 continue;
2812 }
2813
2814 BestTofTrack = (CbmTofTrack*) fTofTracks->At(Btt);
2815 Bthi = BestTofTrack->GetTofHitIndex();
2816 if (verbose > 5) {
2817 cout << "<DBest> GloBTrack " << i << ", TMul: " << NTHMUL[i] << ", TofTrack " << IndTHMUL[i][0]
2818 << ", TofHit " << BestTofTrack->GetTofHitIndex() << ", TMul_hit "
2819 << NTofHitTMul[BestTofTrack->GetTofHitIndex()] << ", dist " << BestTofTrack->GetDistance()
2820 << ", len " << BestTofTrack->GetTrackLength() << ", mass " << BestTofTrack->GetMass() << ", w "
2821 << Weight_THMUL[i][0] << endl;
2822 }
2823 }
2824 }
2825 } //(j > -1) end
2826
2827 // now do global distribution of TofHits to GlobalTracks
2828 // attach BestTofTrack candidate to GlobalTrack
2829 if (verbose > 10) {
2830 cout << Form("<Ddis> NRIt %d, gt %d, BestTofTrack j=%d, best 0x%p, %d, "
2831 "%d, w: %7.2f ",
2832 NRIt, i, j, BestTofTrack, Btt, Bthi, Weight_THMUL[i][0])
2833 << endl;
2834 }
2835
2836 if (NRIt == 1) GlobTrack->SetLength(0.);
2837 if (Btt > -1)
2838 while (j > -1 && GlobTrack->GetLength() != BestTofTrack->GetTrackLength()) {
2839 if (verbose > 10) {
2840 cout << Form("<Ddeb> BestTofTrack j=%d, best 0x%p, %d", j, BestTofTrack, BestTofTrack->GetTofHitIndex())
2841 << endl;
2842 }
2843 dist = BestTofTrack->GetDistance();
2844 if (isinf(dist)) {
2845 cout << "<E2> invalid dist for gt " << i << ", Btt " << Btt << ", d:" << dist << endl;
2846 break;
2847 }
2848 if (dist < DISTMAX && Weight_THMUL[i][0] < WMAX) {
2849 Int_t jh = NTofHitTMul[Bthi]++;
2850 Int_t nht = NTofHitTMul[Bthi];
2851 if (nht == MAXNHT) {
2852 if (verbose > -1) { cout << "<E> Too many TofTrack candidates for hit " << Bthi << ", break!" << endl; }
2853 break;
2854 }
2855 IndTofTrack_TofHit[Bthi][jh] = Btt; // index of TofTrack assigned to specific TofHit
2856 if (verbose > 3) {
2857 cout << "<Ias> GlobTrack " << i << " -> TofTrack " << Btt << ", TofHitIndex " << Bthi << ", TMul_hit "
2858 << nht << endl;
2859 }
2860
2861 Int_t io = -1;
2862 if (NTofHitTMul[BestTofTrack->GetTofHitIndex()] > 1) {
2863 CbmTofTrack* TofTracko =
2864 (CbmTofTrack*) fTofTracks->At(IndTofTrack_TofHit[BestTofTrack->GetTofHitIndex()][0]);
2865 io = TofTracko->GetTrackIndex(); // Global Track index
2866 if (verbose > 2) { // nh-debug
2867 cout << "<D> GlobTrack " << i << ": update TofHitIndex from " << j << " (Mul " << NTofHitTMul[j] << ") "
2868 << " to " << BestTofTrack->GetTofHitIndex() << " (Mul "
2869 << NTofHitTMul[BestTofTrack->GetTofHitIndex()] << ")"
2870 << ", m " << BestTofTrack->GetMass() << ", w " << Weight_THMUL[i][0] << ", cur: tt "
2871 << IndTofTrack_TofHit[Bthi][0] << ", gt " << io << ", w " << Weight_THMUL[io][0] << " ? " << endl;
2872 }
2873
2874 // decide now!
2875 if (Weight_THMUL[i][0] < Weight_THMUL[io][0]) { // new assignment better than old one -> change
2876 if (verbose > 1) { //nh-debug
2877 cout << "<D> New cand. is better, invalidate entry for gt " << io << endl;
2878 }
2879 NReas++;
2880 NTofHitTMul[Bthi]--; // deregister old toftrack
2881 IndTofTrack_TofHit[Bthi][0] = Btt; // update
2882 CbmGlobalTrack* GlobTrack2 = (CbmGlobalTrack*) fGlobalTracks->At(io);
2883 GlobTrack2->SetLength(0.); // signal entry invalid
2884 }
2885 else { // old assignment better than current candidate
2886 if (verbose > 0) { //nh-debug
2887 cout << Form("<D> Stick to old assignment, Bthi %d, TM %d, THM %d", Bthi, NTofHitTMul[Bthi],
2888 NTHMUL[i])
2889 << endl;
2890 }
2891 NTofHitTMul[Bthi]--; // deregister toftrack
2892 if (NTHMUL[i] > 1) { // take next one from list
2893 NTHMUL[i]--;
2894 for (Int_t jth = 0; jth < NTHMUL[i]; jth++) { // shift old entries
2895 IndTHMUL[i][jth] = IndTHMUL[i][jth + 1];
2896 Weight_THMUL[i][jth] = Weight_THMUL[i][jth + 1];
2897 }
2898 Btt = IndTHMUL[i][0]; // next best TofTrack index
2899 BestTofTrack = (CbmTofTrack*) fTofTracks->At(Btt);
2900 Bthi = BestTofTrack->GetTofHitIndex(); // next best TofHit index
2901 }
2902 else { // no other candidate available
2903 if (verbose > 0) { //nh-debug
2904 cout << "<I> no TofTrack candidate for Global Track " << i << endl;
2905 }
2906 // BestTofTrack->Delete();
2907 Bthi = -1;
2908 Btt = -1;
2909 GlobTrack->SetTofHitIndex(-1);
2910 j = -1;
2911 continue;
2912 }
2913 if (verbose > 0) { //nh-debug
2914 cout << "<D> Old choice better, current options: NTHMUL " << NTHMUL[i]
2915 << ", take btt = " << IndTHMUL[i][0] << ", bthi " << Bthi << endl;
2916 }
2917 }
2918
2919 if (NTofHitTMul[BestTofTrack->GetTofHitIndex()] > 1) {
2920 if (verbose > -1) {
2921 cout << "<E> GlobTrack " << i
2922 << ": double assignment of hit, check all possibilities "
2923 "... "
2924 << endl;
2925 continue;
2926 }
2927 }
2928 j = BestTofTrack->GetTofHitIndex();
2929 GlobTrack->SetTofHitIndex(j); // update Global Track info
2930 }
2931 GlobTrack->SetParamLast(BestTofTrack->GetTrackParameter());
2932 GlobTrack->SetLength(BestTofTrack->GetTrackLength());
2933 }
2934 else {
2935 if (verbose > 3) {
2936 cout << "<D> GlobTrack " << i << ", dist = " << dist << ", w = " << Weight_THMUL[i][0]
2937 << " -> remove TofTrack" << endl;
2938 }
2939 GlobTrack->SetTofHitIndex(-1);
2940 j = -1;
2941 }
2942 } // GetTrackLength matching while end
2943 else {
2944 if (verbose > 3) {
2945 cout << "<D> GlobTrack " << i << ", dist = " << dist << ", w = " << Weight_THMUL[i][0] << " -> no TofTrack"
2946 << endl;
2947 }
2948 GlobTrack->SetTofHitIndex(-1);
2949 j = -1;
2950 }
2951 if (verbose > 10) { cout << "<Dch> GlobTrack " << i << "(" << nGlobTracks << "), Btt " << Btt << endl; }
2952 } // endfor of GlobalTrack inspection and TofHit assignment
2953 if (verbose > 1) {
2954 cout << "<Q> Reassignment iteration for b= " << fMCEventHeader->GetB() << ": " << NReas << endl;
2955 }
2956 } //end of reassignment while
2957 // Analysis of GlobalTracks
2958 for (Int_t i = 0; i < nGlobTracks; i++) { // loop over global tracks
2959 GlobTrack = (CbmGlobalTrack*) fGlobalTracks->At(i);
2960 Int_t s = GlobTrack->GetStsTrackIndex();
2961 Int_t j = GlobTrack->GetTofHitIndex();
2962
2963 if (j > -1 && Weight_THMUL[i][0] >= WMAX) {
2964 cout << "<E> TofHit assigned beyond w-limit, Track " << i << " w= " << Weight_THMUL[i][0] << endl;
2965 break; // less drastic response
2966 }
2967
2968 if (verbose > 10) { cout << "<Da> gt " << i << ", th " << j << ", s " << s << endl; }
2969 const FairTrackParam* tparf = GlobTrack->GetParamFirst();
2970 if (0 == tparf->GetQp()) {
2971 if (verbose > 10) cout << "<W2> Global Track " << i << " without Qp!, take from Sts " << s << endl;
2972 StsTrack = (CbmStsTrack*) fStsTracks->At(s);
2973 GlobTrack->SetParamFirst(StsTrack->GetParamFirst());
2974 }
2975 tparf = GlobTrack->GetParamFirst();
2976 Float_t qpf = tparf->GetQp(); //
2977 if (qpf == 0.) {
2978 cout << "<E2> GlobTrack " << i << ", STS " << s << ", TofHit " << j << " without momentum " << endl;
2979 break;
2980 }
2981
2982 // STS
2983 Double_t vtxb = 100.;
2984 Int_t smc = -1;
2985 Int_t StsMCt[100]; // array of MC track indices for current StsTrack
2986 Int_t NStsMCc[100] = {100 * 0}; // number of contributions
2987 Int_t NStsMCt = 0; // number of MC tracks contributing to this Ststrack
2988 if (s > -1) { // STS Track analysis, disable, bad referencing to StsHits
2989 StsTrack = (CbmStsTrack*) fStsTracks->At(s);
2990
2991 FairTrackParam paramExtr;
2992 fTrackFitter.FitToVertex(StsTrack, fPrimVertex, &paramExtr);
2993 vtxb = fTrackFitter.GetChiToVertex(StsTrack,
2994 fPrimVertex); //impact parameter ???
2995 fa_VtxB->Fill(vtxb);
2996
2997 Int_t NStsHits = StsTrack->GetNofStsHits();
2998 //if(NStsHits<8) continue; // nh-debugging
3000 for (Int_t ih = 0; ih < NStsHits; ih++) {
3001 Int_t iHind = StsTrack->GetStsHitIndex(ih);
3002
3003 LOG(debug1) << " inspect STS track " << s << ", hit " << ih << ", hitindex " << iHind;
3004 if (NULL == fStsHits) LOG(fatal) << " No STS Hits available ";
3005 //CbmStsHit* hit = (CbmStsHit*) fStsHits->At(iHind); // still valid ? - ok?
3006 CbmStsHit* hit = dynamic_cast<CbmStsHit*>(fStsHits->At(iHind));
3007 if (NULL == hit) continue;
3008 LOG(debug1) << " valid hit " << ih << ", hitindex " << iHind
3009 << " cluster index f: " << hit->GetFrontClusterId() << ", b: " << hit->GetBackClusterId();
3010
3013 std::stringstream ss;
3014 ss << " Mul f: " << fclu->GetNofDigis() << " (";
3015 for (Int_t iDigi = 0; iDigi < fclu->GetNofDigis(); iDigi++) {
3016 ss << fclu->GetDigi(iDigi) << " ";
3017 //CbmStsDigi* stsdigi = (CbmStsDigi*) fStsDigis->At( fclu->GetDigi(iDigi) );
3018 //CbmMatch* stsdigiMatch = (CbmMatch*) fStsDigiMatchColl->At(fclu->GetDigi(iDigi));
3019 CbmMatch* stsdigiMatch = (CbmMatch*) fDigiMan->GetMatch(ECbmModuleId::kSts, fclu->GetDigi(iDigi));
3020 ss << stsdigiMatch->GetNofLinks() << " ";
3021 for (Int_t iL = 0; iL < stsdigiMatch->GetNofLinks(); iL++) {
3022 const CbmLink& link = stsdigiMatch->GetLink(iL);
3023 CbmStsPoint* poi = (CbmStsPoint*) fStsPointsColl->At(link.GetIndex());
3024 if (NULL == poi) continue;
3025 Int_t MCInd = poi->GetTrackID();
3026 ss << " MCInd " << poi->GetTrackID() << " ";
3027 Int_t iMCt = 0;
3028 for (; iMCt < NStsMCt; iMCt++) {
3029 if (MCInd == StsMCt[iMCt]) {
3030 NStsMCc[iMCt]++;
3031 break;
3032 }
3033 }
3034 if (iMCt == NStsMCt) {
3035 LOG(debug) << "contribution by new MC track: " << MCInd;
3036 StsMCt[iMCt] = MCInd;
3037 NStsMCc[iMCt] = 1;
3038 NStsMCt++;
3039 }
3040 }
3041 }
3042
3043 for (Int_t iDigi = 0; iDigi < bclu->GetNofDigis(); iDigi++) {
3044 ss << bclu->GetDigi(iDigi) << " ";
3045 //CbmStsDigi* stsdigi = (CbmStsDigi*) fStsDigis->At( bclu->GetDigi(iDigi) );
3046 //CbmMatch* stsdigiMatch = (CbmMatch*) fStsDigiMatchColl->At(bclu->GetDigi(iDigi));
3047 CbmMatch* stsdigiMatch = (CbmMatch*) fDigiMan->GetMatch(ECbmModuleId::kSts, bclu->GetDigi(iDigi));
3048 ss << stsdigiMatch->GetNofLinks() << " ";
3049 for (Int_t iL = 0; iL < stsdigiMatch->GetNofLinks(); iL++) {
3050 const CbmLink& link = stsdigiMatch->GetLink(iL);
3051 CbmStsPoint* poi = (CbmStsPoint*) fStsPointsColl->At(link.GetIndex());
3052 Int_t MCInd = poi->GetTrackID();
3053 ss << " MCInd " << poi->GetTrackID() << " ";
3054 Int_t iMCt = 0;
3055 for (; iMCt < NStsMCt; iMCt++) {
3056 if (MCInd == StsMCt[iMCt]) {
3057 NStsMCc[iMCt]++;
3058 break;
3059 }
3060 }
3061 if (iMCt == NStsMCt) {
3062 LOG(debug) << "contribution by new back MC track: " << MCInd;
3063 StsMCt[iMCt] = MCInd;
3064 NStsMCc[iMCt] = 1;
3065 NStsMCt++;
3066 }
3067 }
3068 }
3069
3070 ss << "), mul b: " << bclu->GetNofDigis();
3071 LOG(debug1) << ss.str();
3072 } // loop over STS hits finished
3073
3074 std::stringstream ss;
3075 ss << "STS summary: NStsMCt =" << NStsMCt;
3076 for (Int_t iT = 0; iT < NStsMCt; iT++) {
3077 ss << " iT " << iT << " NMCc " << NStsMCc[iT] << " MCt " << StsMCt[iT];
3078 }
3079 LOG(debug) << ss.str();
3080
3081 if (NStsMCt == 0) {
3082 smc = -1;
3083 LOG(debug) << "StsTrack " << s << " with " << NStsHits << " Hits without StsPoints ??? from Global Track " << i
3084 << ", TofHit " << j;
3085 }
3086 else { // find most probable MCtrack
3087 smc = -1;
3088 Int_t iMaxCount = 0;
3089 for (Int_t k = 0; k < NStsMCt; k++) {
3090 if (NStsMCc[k] > iMaxCount) {
3091 smc = StsMCt[k];
3092 iMaxCount = NStsMCc[k];
3093 // cout << "-D- STS Track "<<smc<<" with "<<NStsMCc[k]<<"("<< NStsMCt<<") matches "
3094 // <<" in "<<k<<". position"<<endl;
3095 // continue;
3096 }
3097 }
3098 }
3099 // analysis of STS tracks
3100 if (-1 < smc) {
3101 MCTrack = (CbmMCTrack*) fMCTracksColl->At(smc);
3102 pdgCode = MCTrack->GetPdgCode();
3103 px_MC = MCTrack->GetPx();
3104 py_MC = MCTrack->GetPy();
3105 pz_MC = MCTrack->GetPz();
3106 p_MC = sqrt(px_MC * px_MC + py_MC * py_MC + pz_MC * pz_MC);
3107 if (MCTrack->GetMotherId() == -1) { // select primaries
3108 // if (0 == MCTrack->GetMass()) continue;
3109 switch (pdgCode) {
3110 case 211: {
3111 fa_ptm_rap_sts_pip->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / MCTrack->GetMass());
3112 fa_plab_sts_pip->Fill(p_MC);
3113 break;
3114 };
3115 case -211: {
3116 fa_ptm_rap_sts_pim->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / MCTrack->GetMass());
3117 fa_plab_sts_pim->Fill(p_MC);
3118 break;
3119 };
3120 case 321: {
3121 fa_ptm_rap_sts_kp->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / MCTrack->GetMass());
3122 fa_plab_sts_kp->Fill(p_MC);
3123 break;
3124 };
3125 case -321: {
3126 fa_ptm_rap_sts_km->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / MCTrack->GetMass());
3127 fa_plab_sts_km->Fill(p_MC);
3128 break;
3129 };
3130 case 2212: {
3131 fa_ptm_rap_sts_p->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / MCTrack->GetMass());
3132 fa_plab_sts_p->Fill(p_MC);
3133 break;
3134 };
3135 case -2212: {
3136 fa_ptm_rap_sts_pbar->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / MCTrack->GetMass());
3137 fa_plab_sts_pbar->Fill(p_MC);
3138 break;
3139 };
3140 case 1000010020: { // deuteron
3141 fa_ptm_rap_sts_d->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / MCTrack->GetMass());
3142 break;
3143 };
3144 case 1000010030: { // triton
3145 fa_ptm_rap_sts_t->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / MCTrack->GetMass());
3146 break;
3147 };
3148 case 1000020030: { // 3he
3149 fa_ptm_rap_sts_h->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / MCTrack->GetMass());
3150 break;
3151 };
3152 case 1000020040: { // alpha
3153 fa_ptm_rap_sts_a->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / MCTrack->GetMass());
3154 break;
3155 };
3156 default:
3157 { // intermediate mass fragments
3158 //cout << " Track k="<<k<<", pdgCode = "<<pdgCode<<
3159 //" Mass " << MCTrack->GetMass()<<","<<MCTrack->GetMass()<<" Y " << MCTrack->GetRapidity() <<
3160 //" Pt " << MCTrack->GetPt() <<endl;
3161 fa_ptm_rap_sts_imf->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / MCTrack->GetMass());
3162 break;
3163 };
3164 }
3165 }
3166 }
3167 }
3168
3169 if (j > -1) { // TofHit available for global track
3170 NGTofTrack++;
3171 TofHit = (CbmTofHit*) fTofHits->At(j);
3172 if (NULL == TofHit) continue;
3173 if (TofHit->GetTime() < 0.2) continue; // skip start counter
3174
3175 Int_t k = -1;
3176
3177 if (NULL == fTofDigiMatchColl) {
3178 LOG(fatal) << "No Digi Info available for TofHit !?? ";
3179 //TofPoint = (CbmTofPoint*) fTofPoints->At( TofHit->GetRefId() );
3180 //k = TofPoint->GetTrackID();
3181 }
3182 else {
3183 CbmMatch* digiMatch = (CbmMatch*) fTofDigiMatchColl->At(j);
3184 // take first digi's point link
3185 Int_t iDigiMul = digiMatch->GetNofLinks();
3186 Int_t iPoiMul = 0;
3187 Int_t iTrkMul = 0;
3188 Int_t iPoiArr[iDigiMul];
3189 Int_t iTrkArr[iDigiMul];
3190 iPoiArr[0] = -1;
3191 iTrkArr[0] = -1;
3192 for (Int_t iLink = 0; iLink < digiMatch->GetNofLinks(); iLink++) { // loop over digis
3193 CbmLink L = digiMatch->GetLink(iLink);
3194 Int_t iDigInd = L.GetIndex();
3196 LOG(error) << "MC-Points not available for DigInd " << iDigInd;
3197 continue;
3198 }
3199 CbmMatch* poiMatch = (CbmMatch*) fDigiMan->GetMatch(ECbmModuleId::kTof, iDigInd);
3200 //CbmMatch* poiMatch = (CbmMatch*) fTofDigiMatchPointsColl->At(iDigInd);
3201 LOG(info) << "Got PoiMatch for TofHit " << j << ", t " << TofHit->GetTime() << ", digi " << iDigInd << ": "
3202 << poiMatch;
3203 if (NULL == poiMatch) continue;
3204
3205 CbmLink LP;
3206 try {
3207 LP = poiMatch->GetMatchedLink();
3208 }
3209 catch (...) {
3210 LOG(info) << "Got invalid PoiMatch for TofHit " << j << ", digi " << iDigInd << ": " << poiMatch;
3211 continue;
3212 }
3213 lp = LP.GetIndex();
3214
3215 if (lp != iPoiArr[iPoiMul]) {
3216 // cout << Form("<D> HadronAnalysis: gt %d, Hit %d, Link %d, poi %d, lpoi %d, PoiMul %d",
3217 // i,j,iLink,lp, iPoiArr[iPoiMul], iPoiMul)<<endl;
3218 iPoiArr[iPoiMul] = lp;
3219 iPoiMul++;
3220 iPoiArr[iPoiMul] = lp;
3221 }
3222
3223 TofPoint = (CbmTofPoint*) fTofPoints->At(lp);
3224 k = TofPoint->GetTrackID();
3225 if (k != iTrkArr[iTrkMul]) {
3226 iTrkArr[iTrkMul] = k;
3227 iTrkMul++;
3228 iTrkArr[iTrkMul] = k;
3229 }
3230 }
3231 /*
3232 if(iTrkMul>1 || iPoiMul>1) {
3233 // cout << Form("HadronAnalysis: McTrkMul %d for TofHit %d, PoiMul %d,",iTrkMul,j,iPoiMul)<<endl;
3234 continue; // for debugging response
3235 }
3236 */
3237 TofPoint = (CbmTofPoint*) fTofPoints->At(iPoiArr[0]);
3238 k = iTrkArr[0];
3239 }
3240 if (NULL == fMCTracksColl) {
3241 LOG(error) << "MC-Tracks not available for k " << k;
3242 continue;
3243 }
3244 MCTrack = (CbmMCTrack*) fMCTracksColl->At(k);
3245 if (NULL == MCTrack) {
3246 LOG(warn) << "No Track for TofHit " << j << ", Tofpoint " << lp << ", TrackId " << k;
3247 continue;
3248 }
3249 pdgCode = MCTrack->GetPdgCode();
3250 px_MC = MCTrack->GetPx();
3251 py_MC = MCTrack->GetPy();
3252 pz_MC = MCTrack->GetPz();
3253 p_MC = sqrt(px_MC * px_MC + py_MC * py_MC + pz_MC * pz_MC);
3254
3255 Double_t len = GlobTrack->GetLength();
3256 const FairTrackParam* tpar = GlobTrack->GetParamFirst();
3257 // FairTrackParam *tpar = GlobTrack->GetParamLast();
3258 if (tpar->GetQp() == 0.) {
3259 cout << "Invalid momentum for global track " << i << " TofHit " << j << endl;
3260 break;
3261 }
3262 Double_t mom = 1. / tpar->GetQp();
3263 if (mom < 0.) mom = -mom;
3264 Float_t vel = TofHit->GetR() / TofHit->GetTime(); // GetR() instead of len
3265 Float_t bet = vel / clight;
3266 Double_t m2 = mom * mom * (1. / bet / bet - 1.);
3267
3268 if (bet > 0.99999) { bet = 0.99999; }
3269 Float_t tofmass = mom / bet * sqrt(1. - bet * bet) * TMath::Sign(1., tpar->GetQp());
3270 // Double_t chi2=0.;//(Double_t)(GlobTrack->GetChi2())/(GlobTrack->GetNDF());
3271 //cout << "GlobTrack-Chi2 "<<GlobTrack->GetChi2()<<", "<<GlobTrack->GetNDF()<<", "<<chi2<<endl;
3272
3273 if (k != smc) {
3274 // cout << " Ana GlobalTrack: MCTrack TOF - STS mismatch: "<< k <<" - "<<smc<<endl;
3275 fa_tm_glomis->Fill(mom, tofmass);
3276 fa_w_mom_glomis->Fill(mom, Weight_THMUL[i][0]);
3277 // continue; // for debugging
3278 }
3279 // if(TofHit->GetRt()<150.) continue; // nh-debugging
3280
3281 fa_xy_glo1->Fill(TofHit->GetX(), TofHit->GetY());
3282 fa_pv_glo->Fill(vel, mom);
3283 fa_tm_glo->Fill(mom, tofmass);
3284 fa_m2mom_glo->Fill(mom * TMath::Sign(1., tpar->GetQp()), m2);
3285 fa_pMCmom_glo->Fill(mom, p_MC);
3286 fa_chi2_mom_glo->Fill(mom, vtxb);
3287 fa_w_mom_glo->Fill(mom, Weight_THMUL[i][0]);
3288 fa_len_mom_glo->Fill(mom, len);
3289 fa_LenDismom_glo->Fill(mom, len - TofHit->GetR());
3290
3291 if (NULL != TofPoint) {
3292 fa_LenMcLenGlomom_glo->Fill(mom, len - TofPoint->GetLength());
3293 fa_LenMcDismom_glo->Fill(mom, TofPoint->GetLength() - TofHit->GetR());
3294 }
3295 /*
3296 fhTofTrkDxsel->Fill(TofTrack->GetTrackDx());
3297 fhTofTrkDysel->Fill(TofTrack->GetTrackDy());
3298 */
3299 if (vtxb < fVtxBMax) {
3300 fa_tm_glovtxb->Fill(mom, tofmass);
3301 fa_m2mom_glovtxb->Fill(mom * TMath::Sign(1., tpar->GetQp()), m2);
3302 }
3303
3304 if (MCTrack->GetMotherId() == -1) { // select primaries
3305 fa_tm_gloprim->Fill(mom, tofmass);
3306 fa_m2mom_gloprim->Fill(mom * TMath::Sign(1., tpar->GetQp()), m2);
3307 fa_chi2_mom_gloprim->Fill(mom, vtxb);
3308
3309 if (vtxb < fVtxBMax) {
3310 fa_tm_gloprimvtxb->Fill(mom, tofmass);
3311 fa_m2mom_gloprimvtxb->Fill(mom * TMath::Sign(1., tpar->GetQp()), m2);
3312 }
3313 Float_t Phip = RADDEG * atan2(MCTrack->GetPy(), MCTrack->GetPx());
3314 Float_t dphi = Phip - RADDEG * fMCEventHeader->GetRotZ();
3315 if (dphi < -180.) { dphi += 360.; };
3316 if (dphi > 180.) { dphi -= 360.; };
3317 dphi = dphi / RADDEG;
3318 rp_weight = 0.;
3319
3320 switch (pdgCode) {
3321 case 211: {
3322 fa_ptm_rap_glo_pip->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / MCTrack->GetMass());
3323 fa_v1_rap_glo_pip->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(dphi));
3324 fa_v2_rap_glo_pip->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(2 * dphi));
3325 fa_xy_glo_pip->Fill(TofHit->GetX(), TofHit->GetY());
3326 fa_tm_glo_pip->Fill(mom, tofmass);
3327 fa_m2mom_glo_pip->Fill(mom * TMath::Sign(1., tpar->GetQp()), m2);
3328 fa_pMCmom_glo_pip->Fill(mom, p_MC);
3329 fa_w_mom_glo_pip->Fill(mom, Weight_THMUL[i][0]);
3330 fa_LenDismom_glo_pip->Fill(mom, len - TofHit->GetR());
3331
3332 if (use_pions_for_flow && TMath::Abs((MCTrack->GetRapidity() - yrp_mid) / yrp_mid) > GetDY()
3333 && MCTrack->GetPt() / MCTrack->GetMass() > GetFlowMinPtm()) {
3334 if (MCTrack->GetRapidity() > yrp_mid) { // set weights for reaction plane
3335 rp_weight = -1.;
3336 }
3337 else {
3338 rp_weight = +1.;
3339 }
3340 }
3341 else {
3342 rp_weight = 0.;
3343 }
3344 break;
3345 };
3346 case -211: {
3347 fa_ptm_rap_glo_pim->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / MCTrack->GetMass());
3348 fa_v1_rap_glo_pim->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(dphi));
3349 fa_v2_rap_glo_pim->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(2 * dphi));
3350 fa_xy_glo_pim->Fill(TofHit->GetX(), TofHit->GetY());
3351 fa_tm_glo_pim->Fill(mom, tofmass);
3352 fa_m2mom_glo_pim->Fill(mom * TMath::Sign(1., tpar->GetQp()), m2);
3353 fa_pMCmom_glo_pim->Fill(mom, p_MC);
3354 fa_w_mom_glo_pim->Fill(mom, Weight_THMUL[i][0]);
3355 fa_LenDismom_glo_pim->Fill(mom, len - TofHit->GetR());
3356
3357 if (use_pions_for_flow && TMath::Abs((MCTrack->GetRapidity() - yrp_mid) / yrp_mid) > GetDY()
3358 && MCTrack->GetPt() / MCTrack->GetMass() > GetFlowMinPtm()) {
3359 if (MCTrack->GetRapidity() > yrp_mid) { // set weights for reaction plane
3360 rp_weight = -1.;
3361 }
3362 else {
3363 rp_weight = +1.;
3364 }
3365 }
3366 else {
3367 rp_weight = 0.;
3368 }
3369 break;
3370 };
3371 case 321: {
3372 fa_ptm_rap_glo_kp->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / MCTrack->GetMass());
3373 fa_v1_rap_glo_kp->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(dphi));
3374 fa_v2_rap_glo_kp->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(2 * dphi));
3375 fa_xy_glo_kp->Fill(TofHit->GetX(), TofHit->GetY());
3376 fa_tm_glo_kp->Fill(mom, tofmass);
3377 fa_m2mom_glo_kp->Fill(mom * TMath::Sign(1., tpar->GetQp()), m2);
3378 fa_pMCmom_glo_kp->Fill(mom, p_MC);
3379 fa_w_mom_glo_kp->Fill(mom, Weight_THMUL[i][0]);
3380 fa_LenDismom_glo_kp->Fill(mom, len - TofHit->GetR());
3381
3382 break;
3383 };
3384 case -321: {
3385 fa_ptm_rap_glo_km->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / MCTrack->GetMass());
3386 fa_v1_rap_glo_km->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(dphi));
3387 fa_v2_rap_glo_km->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(2 * dphi));
3388 fa_xy_glo_km->Fill(TofHit->GetX(), TofHit->GetY());
3389 fa_tm_glo_km->Fill(mom, tofmass);
3390 fa_m2mom_glo_km->Fill(mom * TMath::Sign(1., tpar->GetQp()), m2);
3391 fa_pMCmom_glo_km->Fill(mom, p_MC);
3392 fa_w_mom_glo_km->Fill(mom, Weight_THMUL[i][0]);
3393 fa_LenDismom_glo_km->Fill(mom, len - TofHit->GetR());
3394
3395 break;
3396 };
3397 case 2212: {
3398 fa_ptm_rap_glo_p->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / MCTrack->GetMass());
3399 fa_v1_rap_glo_p->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(dphi));
3400 fa_v2_rap_glo_p->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(2 * dphi));
3401 fa_xy_glo_p->Fill(TofHit->GetX(), TofHit->GetY());
3402 fa_tm_glo_p->Fill(mom, tofmass);
3403 fa_m2mom_glo_p->Fill(mom * TMath::Sign(1., tpar->GetQp()), m2);
3404 fa_pMCmom_glo_p->Fill(mom, p_MC);
3405 fa_w_mom_glo_p->Fill(mom, Weight_THMUL[i][0]);
3406 fa_LenDismom_glo_p->Fill(mom, len - TofHit->GetR());
3407
3408 // reaction plane determination
3409 if (TMath::Abs((MCTrack->GetRapidity() - yrp_mid) / yrp_mid) > GetDY()
3410 && MCTrack->GetPt() / MCTrack->GetMass() > GetFlowMinPtm()) {
3411 if (MCTrack->GetRapidity() > yrp_mid) { // set weights for reaction plane
3412 rp_weight = 1.;
3413 }
3414 else {
3415 rp_weight = -1.;
3416 }
3417 }
3418 else {
3419 rp_weight = 0.;
3420 }
3421 break;
3422 };
3423 case -2212: {
3424 fa_ptm_rap_glo_pbar->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / MCTrack->GetMass());
3425 fa_v1_rap_glo_pbar->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(dphi));
3426 fa_v2_rap_glo_pbar->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(2 * dphi));
3427 fa_xy_glo_pbar->Fill(TofHit->GetX(), TofHit->GetY());
3428 fa_tm_glo_pbar->Fill(mom, tofmass);
3429 fa_m2mom_glo_pbar->Fill(mom * TMath::Sign(1., tpar->GetQp()), m2);
3430 fa_pMCmom_glo_pbar->Fill(mom, p_MC);
3431 fa_w_mom_glo_pbar->Fill(mom, Weight_THMUL[i][0]);
3432 fa_LenDismom_glo_pbar->Fill(mom, len - TofHit->GetR());
3433
3434 break;
3435 };
3436 case 1000010020: { // deuteron
3437 fa_ptm_rap_glo_d->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / MCTrack->GetMass());
3438 fa_v1_rap_glo_d->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(dphi));
3439 fa_v2_rap_glo_d->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(2 * dphi));
3440 fa_xy_glo_d->Fill(TofHit->GetX(), TofHit->GetY());
3441 fa_tm_glo_d->Fill(mom, tofmass);
3442 fa_pMCmom_glo_d->Fill(mom, p_MC);
3443 fa_w_mom_glo_d->Fill(mom, Weight_THMUL[i][0]);
3444
3445 // reaction plane determination
3446 if (TMath::Abs((MCTrack->GetRapidity() - yrp_mid) / yrp_mid) > GetDY()
3447 && MCTrack->GetPt() / MCTrack->GetMass() > GetFlowMinPtm()) {
3448 if (MCTrack->GetRapidity() > yrp_mid) { // set weights for reaction plane
3449 rp_weight = 1.;
3450 }
3451 else {
3452 rp_weight = -1.;
3453 }
3454 }
3455 break;
3456 };
3457 case 1000010030: { // triton
3458 fa_ptm_rap_glo_t->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / MCTrack->GetMass());
3459 fa_v1_rap_glo_t->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(dphi));
3460 fa_v2_rap_glo_t->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(2 * dphi));
3461 fa_xy_glo_t->Fill(TofHit->GetX(), TofHit->GetY());
3462 fa_tm_glo_t->Fill(mom, tofmass);
3463 fa_pMCmom_glo_t->Fill(mom, p_MC);
3464 fa_w_mom_glo_t->Fill(mom, Weight_THMUL[i][0]);
3465
3466 // reaction plane determination
3467 if (TMath::Abs((MCTrack->GetRapidity() - yrp_mid) / yrp_mid) > GetDY()
3468 && MCTrack->GetPt() / MCTrack->GetMass() > GetFlowMinPtm()) {
3469 if (MCTrack->GetRapidity() > yrp_mid) { // set weights for reaction plane
3470 rp_weight = 1.;
3471 }
3472 else {
3473 rp_weight = -1.;
3474 }
3475 }
3476 break;
3477 };
3478 case 1000020030: { // 3he
3479 fa_ptm_rap_glo_h->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / MCTrack->GetMass());
3480 fa_v1_rap_glo_h->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(dphi));
3481 fa_v2_rap_glo_h->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(2 * dphi));
3482 fa_xy_glo_h->Fill(TofHit->GetX(), TofHit->GetY());
3483 fa_tm_glo_h->Fill(mom, tofmass);
3484 fa_pMCmom_glo_h->Fill(mom, p_MC);
3485 fa_w_mom_glo_h->Fill(mom, Weight_THMUL[i][0]);
3486
3487 // reaction plane determination
3488 if (TMath::Abs((MCTrack->GetRapidity() - yrp_mid) / yrp_mid) > GetDY()
3489 && MCTrack->GetPt() / MCTrack->GetMass() > GetFlowMinPtm()) {
3490 if (MCTrack->GetRapidity() > yrp_mid) { // set weights for reaction plane
3491 rp_weight = 1.;
3492 }
3493 else {
3494 rp_weight = -1.;
3495 }
3496 }
3497 break;
3498 };
3499 case 1000020040: { // alpha
3500 fa_ptm_rap_glo_a->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / MCTrack->GetMass());
3501 fa_v1_rap_glo_a->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(dphi));
3502 fa_v2_rap_glo_a->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(2 * dphi));
3503 fa_xy_glo_a->Fill(TofHit->GetX(), TofHit->GetY());
3504 fa_tm_glo_a->Fill(mom, tofmass);
3505 fa_pMCmom_glo_a->Fill(mom, p_MC);
3506 fa_w_mom_glo_a->Fill(mom, Weight_THMUL[i][0]);
3507
3508 // reaction plane determination
3509 if (TMath::Abs((MCTrack->GetRapidity() - yrp_mid) / yrp_mid) > GetDY()
3510 && MCTrack->GetPt() / MCTrack->GetMass() > GetFlowMinPtm()) {
3511 if (MCTrack->GetRapidity() > yrp_mid) { // set weights for reaction plane
3512 rp_weight = 1.;
3513 }
3514 else {
3515 rp_weight = -1.;
3516 }
3517 }
3518 break;
3519 };
3520 default:
3521 { // intermediate mass fragments
3522 //cout << " Track k="<<k<<", pdgCode = "<<pdgCode<<
3523 //" Mass " << MCTrack->GetMass()<<","<<MCTrack->GetMass()<<" Y " << MCTrack->GetRapidity() <<
3524 //" Pt " << MCTrack->GetPt() <<endl;
3525 fa_ptm_rap_glo_imf->Fill(MCTrack->GetRapidity(), MCTrack->GetPt() / MCTrack->GetMass());
3526 fa_v1_rap_glo_imf->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(dphi));
3527 fa_v2_rap_glo_imf->Fill((MCTrack->GetRapidity() - GetMidY()) / GetMidY(), TMath::Cos(2 * dphi));
3528 // reaction plane determination (optimistic)
3529 if (TMath::Abs((MCTrack->GetRapidity() - yrp_mid) / yrp_mid) > GetDY()
3530 && MCTrack->GetPt() / MCTrack->GetMass() > GetFlowMinPtm()) {
3531 if (MCTrack->GetRapidity() > yrp_mid) { // set weights for reaction plane
3532 rp_weight = 1.;
3533 }
3534 else {
3535 rp_weight = -1.;
3536 }
3537 }
3538 break;
3539 };
3540 }
3541 if (rp_weight != 0.) {
3542 if (gRandom->Uniform(1) > 0.5) { //subdivide events into 2 random subevents
3543 Np1++;
3544 Qx1 = Qx1 + rp_weight * MCTrack->GetPx();
3545 Qy1 = Qy1 + rp_weight * MCTrack->GetPy();
3546 }
3547 else {
3548 Np2++;
3549 Qx2 = Qx2 + rp_weight * MCTrack->GetPx();
3550 Qy2 = Qy2 + rp_weight * MCTrack->GetPy();
3551 }
3552 }
3553 }
3554 }
3555 }
3556 if (verbose > 10) { cout << "<D> RP analysis " << Np1 << ", " << Np2 << endl; }
3557 if (Np1 > 0 && Np2 > 0) {
3558 phirp1 = atan2(Qy1, Qx1);
3559 phirp2 = atan2(Qy2, Qx2);
3560 if (fflowFile != NULL) { // subevent RP flattening
3561 TH1F* phirp_glo_fpar = fflowFile->Get<TH1F>("phirps_glo_fpar");
3562 Float_t dphir = 0.;
3563 for (int j = 0; j < 4; j++) {
3564 Float_t i = (float) (j + 1);
3565 dphir += (-phirp_glo_fpar->GetBinContent(j) * TMath::Cos(i * phirp1)
3566 + phirp_glo_fpar->GetBinContent(j + 4) * TMath::Sin(i * phirp1))
3567 / i;
3568 }
3569 phirp1 += dphir;
3570
3571 dphir = 0.;
3572 for (int j = 0; j < 4; j++) {
3573 Float_t i = (float) (j + 1);
3574 dphir += (-phirp_glo_fpar->GetBinContent(j) * TMath::Cos(i * phirp2)
3575 + phirp_glo_fpar->GetBinContent(j + 4) * TMath::Sin(i * phirp2))
3576 / i;
3577 }
3578 phirp2 += dphir;
3579 } // subevent RP flattening end
3580 delrp = (phirp1 - phirp2);
3581 fa_phirps_glo->Fill(phirp1 * RADDEG); // 1D histo
3582 fa_phirps_glo->Fill(phirp2 * RADDEG); // 1D histo
3583 if (NULL != fMCEventHeader) {
3584 // cout << " Impact parameter "<<fMCEventHeader->GetB()<< ", delrp = "<< delrp << endl;
3585 fa_cdrp_b_glo->Fill(fMCEventHeader->GetB(), TMath::Cos(delrp));
3586 delrp = delrp * RADDEG;
3587 if (delrp < -180.) delrp += 360.;
3588 if (delrp > 180.) delrp -= 360.;
3589 fa_drp_b_glo->Fill(fMCEventHeader->GetB(), delrp);
3590 }
3591 phirp = RADDEG * atan2(Qy1 + Qy2, Qx1 + Qx2); // full reaction plane
3592 while (phirp < -180.) {
3593 phirp += 360.;
3594 }
3595 while (phirp > 180.) {
3596 phirp -= 360.;
3597 }
3598 if (fflowFile != NULL) { // RP flattening
3599 TH1F* phirp_glo_fpar = fflowFile->Get<TH1F>("phirp_glo_fpar");
3600 Float_t dphir = 0.;
3601 for (int j = 0; j < 4; j++) {
3602 Float_t i = (float) (j + 1);
3603 //cout << " RP flat par "<< i << ","<<j<< " par " << phirp_glo_fpar->GetBinContent(j)
3604 // << ","<< phirp_glo_fpar->GetBinContent(j+4) << " phirp "<<phirp<<" dphir "<< dphir << endl;
3605 dphir += ((-phirp_glo_fpar->GetBinContent(j) * TMath::Cos(i * phirp / RADDEG)
3606 + phirp_glo_fpar->GetBinContent(j + 4) * TMath::Sin(i * phirp / RADDEG))
3607 / i);
3608 }
3609 //cout << " phirp " << phirp << " dphir " << dphir*RADDEG << endl;
3610
3611 phirp += dphir * RADDEG;
3612 while (phirp < -180.) {
3613 phirp += 360.;
3614 }
3615 while (phirp > 180.) {
3616 phirp -= 360.;
3617 }
3618 } // RP flattening end
3619 if (NULL != fMCEventHeader) {
3620 delrp = phirp - RADDEG * fMCEventHeader->GetRotZ();
3621 while (delrp < -180.) {
3622 delrp += 360.;
3623 }
3624 while (delrp > 180.) {
3625 delrp -= 360.;
3626 }
3627 fa_phirp_glo->Fill(phirp); // 1D histo
3628 fa_delrp_b_glo->Fill(fMCEventHeader->GetB(), delrp);
3629 fa_cdelrp_b_glo->Fill(fMCEventHeader->GetB(), TMath::Cos(delrp / RADDEG));
3630
3631 fa_mul_b_had->Fill(fMCEventHeader->GetB(), NGTofTrack);
3632 }
3633 }
3634
3635 // Hadron level
3636
3637 if (0 == (fEvents % 1000)) {
3638 LOG(info) << "-I- CbmHadronAnalysis::Exec : "
3639 << "event " << fEvents << " processed.";
3640 }
3641 fEvents += 1;
3642}
3643// ------------------------------------------------------------------
3644
3645
3646// ------------------------------------------------------------------
3648{
3649 // Normalisations
3650 LOG(info) << "CbmHadronAnalysis::Finish up with " << fEvents << " analyzed events ";
3651
3652 Double_t sfe = 1. / fEvents;
3653 Double_t sfac = 1.E7;
3654
3655 LOG(info) << "<I> Normalisation factors " << sfe << "," << sfac;
3656
3657 fa_hit_ch->Scale(sfe * sfac);
3658
3659 fa_xy_poi2->Scale(sfe);
3660 fa_xy_poi3->Add(fa_xy_poi2, fa_xy_poi2, sfac, 0.);
3661
3662 fa_xy_hit2->Scale(sfe);
3663 fa_xy_hit3->Add(fa_xy_hit2, fa_xy_hit2, sfac, 0.);
3664
3665 // Finish of the task execution
3666
3668}
3669// ------------------------------------------------------------------
3670
3671
3672// ------------------------------------------------------------------
3674{
3675 // Write histogramms to the file
3676 if (NULL != fHist) {
3677 fHist->cd();
3678 TIter next(gDirectory->GetList());
3679 TH1* h;
3680 TObject* obj;
3681 while ((obj = (TObject*) next())) {
3682 if (obj->InheritsFrom(TH1::Class())) {
3683 h = (TH1*) obj;
3684 //cout << "Write histo " << h->GetTitle() << endl;
3685 h->Write();
3686 }
3687 }
3688 }
3689 LOG(info) << "Histos in " << fHist->GetName();
3690 //fHist->ls();
3691 fHist->Close();
3692}
3693// ------------------------------------------------------------------
3694static Int_t iCandEv = 0;
3696{
3697#include "TLorentzVector.h"
3698#include "TVector3.h"
3699
3700 static TH1F* fhTofChi;
3701 static TH1F* fhDperp;
3702 static TH2F* fhdEdxMul;
3703 static TH2F* fhdEdxMulsec;
3704 static TH1F* fhDTRDprim;
3705 static TH1F* fhDTRDsec;
3706 static TH1F* fhDperp2;
3707 static TH1F* fhDperpS;
3708 static TH1F* fhD0prim;
3709 static TH1F* fhOpAng;
3710 static TH1F* fhDCA;
3711 static TH1F* fhMinv;
3712 static TH1F* fhPathLen;
3713 static TH1F* fhMMom;
3714 static TH1F* fhMIXOpAng;
3715 static TH1F* fhMIXDCA;
3716 static TH1F* fhMIXMinv;
3717 static TH1F* fhMIXPathLen;
3718 static TH1F* fhMIXMMom;
3719 static TH1F* fhMCLamMom;
3720 static TH1F* fhMCPathLen;
3721
3722 static TH2F* fa_ptm_rap_gen_lam;
3723 static TH2F* fa_ptm_rap_rec_lam;
3724 static TH2F* fa_ptm_rap_mix_lam;
3725 //static TH2F* fa_ptm_rap_sub_lam;
3726
3727 static TH2F* fhTofStsXX;
3728 static TH2F* fhTofStsYY;
3729 static TH2F* fhTofStsdXdY;
3730 static TH2F* fhTofStsdXX;
3731 static TH2F* fhTofStsdXY;
3732 static TH2F* fhTofStsdYX;
3733 static TH2F* fhTofStsdYY;
3734 static TH2F* fhTofVelNhit;
3735 static TH1F* fhNT0;
3736 static TH2F* fhStsStsX0Y0;
3737 static TH2F* fhStsStsX0X1;
3738 static TH2F* fhStsStsY0X1;
3739 static TH2F* fhStsStsX0Y1;
3740 static TH2F* fhStsStsY0Y1;
3741 static TH2F* fhStsStsX0X2;
3742 static TH2F* fhStsStsY0X2;
3743 static TH2F* fhStsStsX0Y2;
3744 static TH2F* fhStsStsY0Y2;
3745
3746 static TH3F* fhTofSts1dXXY;
3747 static TH3F* fhTofSts1dYXY;
3748 static TH3F* fhTofSts2dXXY;
3749 static TH3F* fhTofSts2dYXY;
3750
3751 static std::vector<std::list<std::vector<TLorentzVector>>> fvP;
3752 static std::vector<std::list<std::vector<TLorentzVector>>> fvX;
3753 static std::vector<std::list<std::vector<TVector3>>> fvX0;
3754 static std::vector<std::list<std::vector<TVector3>>> fvDX;
3755
3756 /*
3757 Double_t fdDistPrimLim =1.5; // Ext Parameter: Max Tof-Sts trans distance for primaries
3758 Double_t fdDistPrimLim2=0.3; // Ext Parameter: Max Sts-Sts trans distance for primaries
3759 Double_t fdDistSecLim2=0.5; // Ext Parameter: Max Sts-Sts trans distance from TOF direction for secondaries
3760 Double_t fdD0ProtLim=0.4; // Ext Parameter: Min impact parameter for secondary proton
3761 Double_t fdOpAngMin=0.01; // Ext Parameter: Min opening angle for accepting pair
3762 Double_t fdDCALim=0.2; // Ext Parameter: Max DCA for accepting pair
3763 Double_t fdVLenMax=25.; // Ext Parameter: Max Lambda flight path length for accepting pair
3764 Double_t fdDistTRD = 10.; // max accepted distance of Trd Hit from STS-TOF line
3765 Double_t fdTRDHmulMin = 0.; // min associated Trd Hits to Track candidates
3766 */
3767
3768
3769 const Int_t fiNMixClasses = 10;
3770 const Double_t beamRotY = -25.;
3771 const Double_t MLAM = 1.1156;
3772 const Double_t DMLAM = 0.015;
3773 const Int_t NSECMASS = 2; // pi-minus, proton, he3, alpha
3774 const Int_t iMode = 0;
3775 Float_t secMass[NSECMASS] = {0.139, 0.938};
3776 switch (iMode) {
3777 case 0: // Lambda
3778 break;
3779 case 1: // hypertriton
3780 secMass[1] = 2.808381; //3he
3781 break;
3782 }
3783 const Double_t dTofSigX = 0.5;
3784 const Double_t dTofSigY = 0.8;
3785 const Double_t dTofSigT = 0.08;
3786 const Double_t dChiTofLim = 3.;
3787
3788 nStsHits = 0;
3789 if (NULL != fStsHits) nStsHits = fStsHits->GetEntriesFast();
3790 Int_t nTrdHits = 0;
3791 if (NULL != fTrdHits) nTrdHits = fTrdHits->GetEntriesFast();
3792
3793 LOG(debug) << "Secondaries from " << nTofHits << " TofHits, " << nStsHits << " StsHits and " << nTrdHits
3794 << " TrdHits in event " << iCandEv;
3795
3796 if (iCandEv == 0) { //initialize
3797 // define some histograms
3798 TDirectory* oldDir = gDirectory;
3799 fHist->ReOpen("Update");
3800
3801 Double_t MinvMin = secMass[0] + secMass[1];
3802 cout << "Add secondary histos to " << fHist->GetName() << endl;
3803 fhTofChi = new TH1F(Form("hTofChi"), Form("TofHit Merger; #chi "), 100, 0., dChiTofLim * 2.);
3804 fhDperp = new TH1F(Form("hDperp"), Form("transverse matching distance; d [cm]"), 100, 0., 5.);
3805 fhdEdxMul =
3806 new TH2F(Form("hdEdxMul"), Form("average energy loss vs TrdHitMul; TrdHitMul; dE []"), 4, 1, 5, 100, 0., 5.E-5);
3807 fhdEdxMulsec = new TH2F(Form("hdEdxMulsec"), Form("av. energy loss vs TrdHitMul for secondaries; TrdHitMul; dE []"),
3808 4, 1, 5, 100, 0., 5.E-5);
3809 fhDTRDprim =
3810 new TH1F(Form("hDTRDprim"), Form("TRD transverse matching distance (prim); d [cm]"), 100, 0., 2. * fdDistTRD);
3811 fhDTRDsec =
3812 new TH1F(Form("hDTRDsec"), Form("TRD transverse matching distance (sec); d [cm]"), 100, 0., 2. * fdDistTRD);
3813
3814 fhDperp2 = new TH1F(Form("hDperp2"), Form("transverse matching distance (prim); d [cm]"), 100, 0., 1.);
3815 fhDperpS = new TH1F(Form("hDperpS"), Form("transverse matching distance (sec); d [cm]"), 100, 0., 1.);
3816 fhD0prim = new TH1F(Form("hD0prim"), Form("transverse distance to primary vertex; d [cm]"), 100, 0., 2.);
3817 fhOpAng = new TH1F(Form("hOpAng"), Form("opening angle; #alpha [rad]"), 100, 0., 0.5);
3818 fhDCA = new TH1F(Form("hDCA"), Form("distance of closest approach; d [cm]"), 100, 0., 2.);
3819 fhMinv = new TH1F(Form("hMinv"), Form("invariant mass; M_{inv} [GeV]"), 100, MinvMin, MinvMin + 0.2);
3820 fhPathLen = new TH1F(Form("hPathLen"), Form("path length; L [cm]"), 100, 0., 30.);
3821 fhMMom = new TH1F(Form("hMMom"), Form("momentum of mother ; p [GeV]"), 100, 0., 5.);
3822 fhMIXOpAng = new TH1F(Form("hMIXOpAng"), Form("opening angle; #alpha [rad]"), 100, 0., 0.5);
3823 fhMIXDCA = new TH1F(Form("hMIXDCA"), Form("distance of closest approach; d [cm]"), 100, 0., 2.);
3824 fhMIXMinv = new TH1F(Form("hMIXMinv"), Form("invariant mass; M_{inv} [GeV]"), 100, MinvMin, MinvMin + 0.2);
3825 fhMIXPathLen = new TH1F(Form("hMIXPathLen"), Form("path length; L [cm]"), 100, 0., 30.);
3826 fhMIXMMom = new TH1F(Form("hMIXMMom"), Form("momentum of mother ; p [GeV]"), 100, 0., 5.);
3827 fhMCPathLen = new TH1F(Form("hMCPathLen"), Form("MC hyperon path length; L [cm]"), 100, 0., 30.);
3828 fhMCLamMom = new TH1F(Form("hMCLamMom"), Form("MC hyperon momentum; p [GeV]"), 100, 0., 5.);
3829
3830 fhTofStsXX = new TH2F(Form("hTofStsXX"), Form("Position correlation XX; x_{Sts} [cm]; x_{TofExt} [cm]"), 200, -10.,
3831 10., 200, -15., 15.);
3832 fhTofStsYY = new TH2F(Form("hTofStsYY"), Form("Position correlation YY; y_{Sts} [cm]; y_{TofExt} [cm]"), 200, -10.,
3833 10., 200, -15., 15.);
3834 fhTofStsdXdY =
3835 new TH2F(Form("hTofStsdXdY"), Form("Position correlation dXX; #Deltax_{Tof-Sts} [cm]; #Deltay_{Tof-Sts} [cm]"),
3836 100, -2., 2., 100, -2., 2.);
3837 fhTofStsdXX = new TH2F(Form("hTofStsdXX"), Form("Position correlation dXX; x_{Tof} [cm]; #Deltax_{Tof-Sts} [cm]"),
3838 200, -15., 15., 100, -2., 2.);
3839 fhTofStsdXY = new TH2F(Form("hTofStsdXY"), Form("Position correlation dXY; y_{Tof} [cm]; #Deltax_{Tof-Sts} [cm]"),
3840 200, -15., 15., 100, -2., 2.);
3841 fhTofStsdYX = new TH2F(Form("hTofStsdYX"), Form("Position correlation dYX; x_{Tof} [cm]; #Deltay_{Tof-Sts} [cm]"),
3842 200, -15., 15., 100, -2., 2.);
3843 fhTofStsdYY = new TH2F(Form("hTofStsdYY"), Form("Position correlation dYY; y_{Tof} [cm]; #Deltay_{Tof-Sts} [cm]"),
3844 200, -15., 15., 100, -2., 2.);
3845 fhTofVelNhit = new TH2F(Form("hTofVelNhit"), Form("Velocity vs Nhit; Nhit []; v [cm/ns]"), 10, 0, 10, 100, 0, 50.);
3846 fhNT0 = new TH1F(Form("hNT0"), Form("Number of T0 signal per event; N []; cts []"), 10, 0, 10);
3847 fhStsStsX0Y0 = new TH2F(Form("hStsStsX0Y0"), Form("vertex location; x0_{StsSts} [cm]; y0_{StsSts} [cm]"), 100, -2.,
3848 2., 100, -2., 2.);
3849
3850 fhStsStsX0X1 = new TH2F(Form("hStsStsX0X1"), Form("vertex position dependence; x1_{Sts} [cm]; x0_{StsSts} [cm]"),
3851 100, -10., 10., 100, -2., 2.);
3852 fhStsStsY0X1 = new TH2F(Form("hStsStsY0X1"), Form("vertex position dependence; x1_{Sts} [cm]; y0_{StsSts} [cm]"),
3853 100, -10., 10., 100, -2., 2.);
3854 fhStsStsX0Y1 = new TH2F(Form("hStsStsX0Y1"), Form("vertex position dependence; y1_{Sts} [cm]; x0_{StsSts} [cm]"),
3855 100, -10., 10., 100, -2., 2.);
3856 fhStsStsY0Y1 = new TH2F(Form("hStsStsY0Y1"), Form("vertex position dependence; y1_{Sts} [cm]; y0_{StsSts} [cm]"),
3857 100, -10., 10., 100, -2., 2.);
3858
3859 fhStsStsX0X2 = new TH2F(Form("hStsStsX0X2"), Form("vertex position dependence; x2_{Sts} [cm]; x0_{StsSts} [cm]"),
3860 100, -10., 10., 100, -2., 2.);
3861 fhStsStsY0X2 = new TH2F(Form("hStsStsY0X2"), Form("vertex position dependence; x2_{Sts} [cm]; y0_{StsSts} [cm]"),
3862 100, -10., 10., 100, -2., 2.);
3863 fhStsStsX0Y2 = new TH2F(Form("hStsStsX0Y2"), Form("vertex position dependence; y2_{Sts} [cm]; x0_{StsSts} [cm]"),
3864 100, -10., 10., 100, -2., 2.);
3865 fhStsStsY0Y2 = new TH2F(Form("hStsStsY0Y2"), Form("vertex position dependence; y2_{Sts} [cm]; y0_{StsSts} [cm]"),
3866 100, -10., 10., 100, -2., 2.);
3867
3868 fhTofSts1dXXY =
3869 new TH3F(Form("hTofSts1dXXY"), Form("Position deviation dXXY; x_{Sts1} [cm]; y_{Sts1}; #Deltax_{Tof-Sts} [cm]"),
3870 100, -10., 10., 100, 10, 10, 50, -2., 2.);
3871 fhTofSts1dYXY =
3872 new TH3F(Form("hTofSts1dYXY"), Form("Position deviation dYXY; x_{Sts1} [cm]; y_{Sts1}; #Deltay_{Tof-Sts} [cm]"),
3873 100, -10., 10., 100, 10, 10, 50, -2., 2.);
3874 fhTofSts2dXXY =
3875 new TH3F(Form("hTofSts2dXXY"), Form("Position deviation dXXY; x_{Sts1} [cm]; y_{Sts1}; #Deltax_{Tof-Sts} [cm]"),
3876 100, -10., 10., 100, 10, 10, 50, -2., 2.);
3877 fhTofSts2dYXY =
3878 new TH3F(Form("hTofSts2dYXY"), Form("Position deviation dYXY; x_{Sts1} [cm]; y_{Sts1}; #Deltay_{Tof-Sts} [cm]"),
3879 100, -10., 10., 100, 10, 10, 50, -2., 2.);
3880
3881 // physics observables
3882 Float_t ymin = -1.;
3883 Float_t ymax = 4.;
3884 Float_t ptmmax = 2.5;
3885 Int_t ptm_nbx = 30;
3886 Int_t ptm_nby = 30;
3887
3888 fa_ptm_rap_gen_lam =
3889 new TH2F("ptm_rap_gen_lam", "MCTrack-gen lam; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
3890 fa_ptm_rap_rec_lam = new TH2F("ptm_rap_rec_lam", "rec lam; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
3891 fa_ptm_rap_mix_lam = new TH2F("ptm_rap_mix_lam", "mix lam; y; p_{T}/m", ptm_nbx, ymin, ymax, ptm_nby, 0., ptmmax);
3892 // fa_ptm_rap_sub_lam = new TH2F("ptm_rap_sub_lam","sub lam; y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
3893 gDirectory = oldDir;
3894 }
3895 iCandEv++; // count events locally
3896
3897 // fill generator distributions for reference
3898 for (Int_t k = 0; k < nMCTracks; k++) { // inspect MCTracks
3899 CbmMCTrack* MCTrack = (CbmMCTrack*) fMCTracksColl->At(k);
3900 Int_t pdgCode = MCTrack->GetPdgCode();
3901 //cout<<"MCTrack pdg "<< pdgCode << ", Mother "<<MCTrack->GetMotherId()<<endl;
3902 if (MCTrack->GetMotherId() == -1 && pdgCode == 3122) {
3903 fhMCLamMom->Fill(MCTrack->GetP());
3904 TLorentzVector PLAM(MCTrack->GetPx(), MCTrack->GetPy(), MCTrack->GetPz(), MCTrack->GetEnergy());
3905 PLAM.RotateY(beamRotY * TMath::Pi() / 180.);
3906 fa_ptm_rap_gen_lam->Fill(PLAM.Rapidity(), PLAM.Pt() / MCTrack->GetMass());
3907 }
3908 if (MCTrack->GetMotherId() > -1 && pdgCode == -211) { // decay pion
3909 CbmMCTrack* MCTrackm = (CbmMCTrack*) fMCTracksColl->At(MCTrack->GetMotherId());
3910 if (MCTrackm->GetPdgCode() == 3122) {
3911 TVector3 MCV;
3912 MCTrack->GetStartVertex(MCV);
3913 fhMCPathLen->Fill(MCV.Mag());
3914 LOG(debug) << "MC vertex at Pathlen = " << MCV.Mag() << ", pi-mom " << MCTrack->GetP();
3915 }
3916 }
3917 }
3918
3919 fvP.resize(fiNMixClasses);
3920 fvX.resize(fiNMixClasses);
3921 fvX0.resize(fiNMixClasses);
3922 fvDX.resize(fiNMixClasses);
3923
3924 Int_t iMixClass = nTofHits * fiNMixClasses / fiTofHitMulMax;
3925 if (iMixClass >= fiNMixClasses) iMixClass = fiNMixClasses - 1;
3926
3927 std::vector<TLorentzVector> P;
3928 P.resize(nTofHits); // define array of momentum Lorentzvectors
3929 std::vector<TLorentzVector> X;
3930 X.resize(nTofHits); // define array of space Lorentzvectors
3931 std::vector<TVector3> X0;
3932 X0.resize(nTofHits); // first measured point of line
3933 std::vector<TVector3> DX;
3934 DX.resize(nTofHits); // direction of line
3935 const Int_t NTrdStations = 4;
3936 std::vector<std::vector<Int_t>> iTRD;
3937 iTRD.resize(nTofHits);
3938
3939 Double_t dStsDistMin[nTofHits];
3940 Double_t dSts2DistMin[nTofHits];
3941 Double_t dTofDistMin[nStsHits];
3942 Double_t dTofDist2Min[nStsHits];
3943 Int_t iStsMin[nTofHits][2]; // storage of track candidate STS hits
3944 Int_t iTofMin[nStsHits];
3945
3946 Double_t dTrdDistMin[nTofHits][NTrdStations];
3947
3948 Int_t proton_cand = 0;
3949 Int_t pion_cand = 0;
3950
3951 for (Int_t j = 0; j < nStsHits; j++) {
3952 iTofMin[j] = -1; //initialize
3953 dTofDistMin[j] = 100.; //initialize
3954 dTofDist2Min[j] = 100.; //initialize
3955 }
3956
3957 //-2. find T0
3958 Double_t dT0 = 0.;
3959 Double_t dNT0 = 0;
3960 for (Int_t i = 0; i < nTofHits; i++) {
3961 CbmTofHit* pTofHit = (CbmTofHit*) fTofHits->At(i);
3962 if (NULL == pTofHit) continue;
3963 if (pTofHit->GetZ() == 0.) {
3964 dT0 += pTofHit->GetTime();
3965 dNT0++;
3966 }
3967 }
3968 fhNT0->Fill(dNT0);
3969 if (dNT0 > 0) {
3970 dT0 /= dNT0;
3971 LOG(debug) << "Found T0 " << dT0 << " from " << dNT0 << " signals ";
3972 }
3973 else
3974 return;
3975
3976 //-1. prepare TOF hits for counter overlaps and substract T0
3977 for (Int_t i = 0; i < nTofHits; i++) {
3978 CbmTofHit* pTofHit = (CbmTofHit*) fTofHits->At(i);
3979 if (NULL == pTofHit) continue;
3980 if (pTofHit->GetZ() == 0.) continue; // don't merge with fake beam counter
3981 pTofHit->SetFlag(1);
3982 pTofHit->SetTime(pTofHit->GetTime() - dT0);
3983 }
3984
3985 //0. merge TOF hits due to counter overlaps
3986 for (Int_t i = 0; i < nTofHits; i++) {
3987 CbmTofHit* pTofHit = (CbmTofHit*) fTofHits->At(i);
3988 if (NULL == pTofHit) continue;
3989 if (pTofHit->GetZ() == 0.) continue; // don't merge with fake beam counter
3990 for (Int_t i2 = 0; i2 < nTofHits; i2++) {
3991 if (i2 != i) {
3992 CbmTofHit* pTofHit2 = (CbmTofHit*) fTofHits->At(i2);
3993 if (NULL == pTofHit2) continue;
3994 if (pTofHit2->GetZ() == 0.) continue; // don't merge with fake beam counter
3995 // Project to plane with smallest z coordinate
3996 if (pTofHit2->GetZ() < pTofHit->GetZ()) { //invert order
3997 CbmTofHit* pTofHittmp = pTofHit;
3998 pTofHit = pTofHit2;
3999 pTofHit2 = pTofHittmp;
4000 }
4001 Double_t dPosZExp = pTofHit->GetZ() / pTofHit2->GetZ();
4002 Double_t dPosXExp = pTofHit2->GetX() * dPosZExp;
4003 Double_t dPosYExp = pTofHit2->GetY() * dPosZExp;
4004 Double_t dTimeExp = pTofHit2->GetTime() * dPosZExp;
4005 Double_t dChi2 = TMath::Power((dPosXExp - pTofHit->GetX()) / dTofSigX, 2)
4006 + TMath::Power((dPosYExp - pTofHit->GetY()) / dTofSigY, 2)
4007 + TMath::Power((dTimeExp - pTofHit->GetTime()) / dTofSigT, 2);
4008 Double_t dChi = TMath::Sqrt(dChi2) / 3.;
4009 fhTofChi->Fill(dChi);
4010 if (dChi < dChiTofLim) { // merge info
4011 pTofHit->SetTime((pTofHit->GetTime() * pTofHit->GetFlag() + dTimeExp)
4012 / (pTofHit->GetFlag() + 1)); // update time
4013 pTofHit->SetFlag(pTofHit->GetFlag() + 1);
4014 pTofHit2->SetFlag(0);
4015 fTofHits->Remove(pTofHit2);
4016 //pTofHit2->Delete(); // remove from TClonesArray
4017 LOG(debug) << "Tof Hits " << i << " and " << i2 << " merged ";
4018 LOG(debug) << "Tof " << i << ", xyz " << pTofHit->GetX() << ", " << pTofHit->GetY() << ", "
4019 << pTofHit->GetZ();
4020 }
4021 }
4022 }
4023 fhTofVelNhit->Fill(pTofHit->GetFlag(), pTofHit->GetR() / pTofHit->GetTime());
4024 }
4025
4026
4027 //1. find best tof silicon match for primary track hypothesis
4028 for (Int_t i = 0; i < nTofHits; i++) {
4029 CbmTofHit* pTofHit = (CbmTofHit*) fTofHits->At(i);
4030 if (NULL == pTofHit) {
4031 P[i].SetPxPyPzE(0., 0., 0., 0.); // avoid nan
4032 continue;
4033 }
4034 if (pTofHit->GetZ() == 0.) continue; // skip fake beam counter
4035 dStsDistMin[i] = 1.E3;
4036 dSts2DistMin[i] = 1.E3;
4037 Double_t dTofX = 0, dTofY = 0, dStsX = 0, dStsY = 0;
4038 for (Int_t l = 0; l < NTrdStations; l++)
4039 dTrdDistMin[i][l] = 1.E3;
4040 iStsMin[i][0] = -1;
4041 iStsMin[i][1] = -1;
4042 for (Int_t j = 0; j < nStsHits; j++) {
4043 CbmStsHit* pStsHit = (CbmStsHit*) fStsHits->At(j);
4044 // Check for primary track
4045 Double_t sPosZ = pStsHit->GetZ();
4046 Double_t sPosXext = pTofHit->GetX() * sPosZ / pTofHit->GetZ();
4047 Double_t sPosYext = pTofHit->GetY() * sPosZ / pTofHit->GetZ();
4048 Double_t dDist2 = TMath::Power(pStsHit->GetX() - sPosXext, 2) + TMath::Power(pStsHit->GetY() - sPosYext, 2);
4049 Double_t dDist = TMath::Sqrt(dDist2);
4050
4051 LOG(debug) << "Sts " << j << ", xyz " << pStsHit->GetX() << ", " << pStsHit->GetY() << ", " << sPosZ;
4052 LOG(debug) << "Tof " << i << ", xyz " << pTofHit->GetX() << ", " << pTofHit->GetY() << ", " << pTofHit->GetZ();
4053 LOG(debug) << "Tof " << i << ", Sts " << j << Form(" -> dist %6.3f, Min %6.3f ", dDist, dStsDistMin[i]);
4054
4055 if (dDist < fdDistPrimLim && dDist < dStsDistMin[i]) { // primary or proton candidate
4056 if (iTofMin[j] > -1) {
4057 LOG(debug) << Form("Sts hit %d already assigned to tof hit %d with "
4058 "dist= %6.3f, prev %6.3f",
4059 j, iTofMin[j], dDist, dTofDistMin[j]);
4060 if (dDist > dTofDistMin[j]) continue; // previous assignment was better
4061 }
4062 dStsDistMin[i] = dDist;
4063 iTofMin[j] = i;
4064 dTofDistMin[j] = dDist;
4065 LOG(debug) << "Prim Track cand started for Tof " << i << ", Sts " << j
4066 << Form(": dist %6.3f, Min %6.3f at z = %4.1f", dDist, dStsDistMin[i], sPosZ);
4067 dTofX = sPosXext;
4068 dTofY = sPosYext;
4069 dStsX = pStsHit->GetX();
4070 dStsY = pStsHit->GetY();
4071 }
4072 } // for (Int_t j=0; j<nStsHits; j++) {
4073 // diagnose closest pair
4074 if (dStsDistMin[i] < 100.) {
4075 fhDperp->Fill(dStsDistMin[i]);
4076 fhTofStsXX->Fill(dStsX, dTofX);
4077 fhTofStsYY->Fill(dStsY, dTofY);
4078 fhTofStsdXdY->Fill(dTofX - dStsX, dTofY - dStsY);
4079 fhTofStsdXX->Fill(dTofX, dTofX - dStsX);
4080 fhTofStsdXY->Fill(dTofY, dTofX - dStsX);
4081 fhTofStsdYX->Fill(dTofX, dTofY - dStsY);
4082 fhTofStsdYY->Fill(dTofY, dTofY - dStsY);
4083 }
4084 } //for (Int_t i=0; i<nTofHits; i++) {
4085
4086 //2.: find second silicon hit for primary tracks
4087 for (Int_t j = 0; j < nStsHits; j++) {
4088 if (iTofMin[j] < 0) continue;
4089 Int_t i = iTofMin[j]; // index of Tof Hit
4090 CbmStsHit* pStsHit = (CbmStsHit*) fStsHits->At(j);
4091 // Check for confirmation of primary track
4092 for (Int_t k = 0; k < nStsHits; k++) {
4093 if (j == k) continue;
4094 CbmStsHit* pSts2Hit = (CbmStsHit*) fStsHits->At(k);
4095 Double_t sPos2Z = pSts2Hit->GetZ();
4096 Double_t sPos2Xext = pStsHit->GetX() * sPos2Z / pStsHit->GetZ();
4097 Double_t sPos2Yext = pStsHit->GetY() * sPos2Z / pStsHit->GetZ();
4098 Double_t dDist2 = TMath::Power(pSts2Hit->GetX() - sPos2Xext, 2) + TMath::Power(pSts2Hit->GetY() - sPos2Yext, 2);
4099 Double_t dDist = TMath::Sqrt(dDist2);
4100 LOG(debug) << "Tof " << i << ", Sts " << j
4101 << Form(" Sts2 %d -> dist %6.3f, Min %6.3f at z = %4.1f", k, dDist, dSts2DistMin[i], sPos2Z);
4102
4103 if (dDist < fdDistPrimLim2 && dDist < dSts2DistMin[i]) { // look for primary or proton candidate
4104 if (iTofMin[k] > -1) {
4105 LOG(debug) << Form("Sts2hit %d already assigned to tof hit %d with "
4106 "dist= %6.3f, prev %6.3f",
4107 k, iTofMin[k], dDist, dTofDistMin[k]);
4108 if (dDist > dTofDist2Min[k]) continue; // previous assignment was better
4109 }
4110 dSts2DistMin[i] = dDist;
4111 iTofMin[k] = i;
4112 if (pStsHit->GetZ() < pSts2Hit->GetZ()) {
4113 iStsMin[i][0] = j;
4114 iStsMin[i][1] = k;
4115 }
4116 else {
4117 iStsMin[i][0] = k;
4118 iStsMin[i][1] = j;
4119 }
4120 dTofDistMin[k] = dDist;
4121 LOG(debug) << "Prim Track cand found for Tof " << i << ", Sts " << j << ", Sts2 " << k
4122 << Form(": dist %6.3f, Min %6.3f at z = %4.1f", dDist, dSts2DistMin[i], sPos2Z);
4123 } // primary or proton candidate
4124 } // for (Int_t k=0; k<nStsHits; k++) {
4125 if (dSts2DistMin[i] < fdDistPrimLim2) {
4126 fhDperp2->Fill(dSts2DistMin[i]);
4127 // diagnose Sts
4128 CbmStsHit* pSts0 = (CbmStsHit*) fStsHits->At(iStsMin[i][0]);
4129 CbmStsHit* pSts1 = (CbmStsHit*) fStsHits->At(iStsMin[i][1]);
4130 Double_t dX0 = pSts0->GetX() - (pSts1->GetX() - pSts0->GetX()) / (pSts1->GetZ() - pSts0->GetZ()) * pSts0->GetZ();
4131 Double_t dY0 = pSts0->GetY() - (pSts1->GetY() - pSts0->GetY()) / (pSts1->GetZ() - pSts0->GetZ()) * pSts0->GetZ();
4132 fhStsStsX0Y0->Fill(dX0, dY0);
4133 fhStsStsX0X1->Fill(pSts0->GetX(), dX0);
4134 fhStsStsX0Y1->Fill(pSts0->GetY(), dX0);
4135 fhStsStsY0X1->Fill(pSts0->GetX(), dY0);
4136 fhStsStsY0Y1->Fill(pSts0->GetY(), dY0);
4137 fhStsStsX0X2->Fill(pSts1->GetX(), dX0);
4138 fhStsStsX0Y2->Fill(pSts1->GetY(), dX0);
4139 fhStsStsY0X2->Fill(pSts1->GetX(), dY0);
4140 fhStsStsY0Y2->Fill(pSts1->GetY(), dY0);
4141
4142 CbmTofHit* pTofHit = (CbmTofHit*) fTofHits->At(i);
4143 fhTofSts1dXXY->Fill(pSts0->GetX(), pSts0->GetY(),
4144 pTofHit->GetX() * pSts0->GetZ() / pTofHit->GetZ() - pSts0->GetX());
4145 fhTofSts1dYXY->Fill(pSts0->GetX(), pSts0->GetY(),
4146 pTofHit->GetY() * pSts0->GetZ() / pTofHit->GetZ() - pSts0->GetY());
4147 fhTofSts2dXXY->Fill(pSts1->GetX(), pSts1->GetY(),
4148 pTofHit->GetX() * pSts1->GetZ() / pTofHit->GetZ() - pSts1->GetX());
4149 fhTofSts2dYXY->Fill(pSts1->GetX(), pSts1->GetY(),
4150 pTofHit->GetY() * pSts1->GetZ() / pTofHit->GetZ() - pSts1->GetY());
4151 }
4152 } // for (Int_t j=0; j<nStsHits; j++) {
4153
4154 //2.a - confirm primary track hypothesis by TRD
4155 if (nTrdHits > 0) {
4156 for (Int_t i = 0; i < nTofHits; i++) {
4157 Int_t j = iStsMin[i][1]; // index of STS hit in second plane
4158 if (j < 0) continue; // no STS assigned
4159 CbmTofHit* pTofHit = (CbmTofHit*) fTofHits->At(i);
4160 if (NULL == pTofHit) continue;
4161 if (pTofHit->GetZ() == 0) continue; // skip fake beam counter
4162 CbmStsHit* pStsHit = (CbmStsHit*) fStsHits->At(j);
4163 Double_t dDx = pTofHit->GetX() - pStsHit->GetX();
4164 Double_t dDy = pTofHit->GetY() - pStsHit->GetY();
4165 Double_t dDz = pTofHit->GetZ() - pStsHit->GetZ();
4166 LOG(debug) << "Check for TRD hits between STS " << j << " and TOF " << i;
4167
4168 for (Int_t l = 0; l < nTrdHits; l++) {
4169 CbmTrdHit* pTrdHit = (CbmTrdHit*) fTrdHits->At(l);
4170 // calculate expected position in Trd layer
4171 Double_t dXexp = pStsHit->GetX() + dDx * (pTrdHit->GetZ() - pStsHit->GetZ()) / dDz;
4172 Double_t dYexp = pStsHit->GetY() + dDy * (pTrdHit->GetZ() - pStsHit->GetZ()) / dDz;
4173 Double_t dDtrans =
4174 TMath::Sqrt(TMath::Power(dXexp - pTrdHit->GetX(), 2) + TMath::Power(dYexp - pTrdHit->GetY(), 2));
4175 UInt_t iTrdLayer = CbmTrdAddress::GetLayerId(pTrdHit->GetAddress());
4176 LOG(debug) << "Inspect TRD hit " << l << " in "
4177 << Form("Module 0x%08x, layer %d", pTrdHit->GetAddress(),
4179 << " at z= " << pTrdHit->GetZ() << " dD = " << dDtrans << " < " << fdDistTRD;
4180 fhDTRDprim->Fill(dDtrans);
4181 if (dDtrans < fdDistTRD && dDtrans < dTrdDistMin[i][iTrdLayer]) { // check if acceptable and take best match
4182 Int_t iMul = iTRD[i].size();
4183 if (dTrdDistMin[i][iTrdLayer] < 1.E3) { // modify previous entry
4184 //find old entry in vector
4185 Int_t ll = 0;
4186 for (; ll < iMul; ll++)
4187 if (CbmTrdAddress::GetLayerId(((CbmTrdHit*) fTrdHits->At(iTRD[i][ll]))->GetAddress()) == iTrdLayer) break;
4188 iTRD[i][ll] = l;
4189 }
4190 else { //add hit
4191 dTrdDistMin[i][iTrdLayer] = dDtrans;
4192 iTRD[i].resize(iMul + 1);
4193 iTRD[i][iMul] = l;
4194 }
4195 LOG(debug) << "assign TrdHit " << l << " to TofHit " << i << " in layer " << iTrdLayer
4196 << " with d = " << dDtrans << ", TrdMul" << iMul << ", dEdx = " << pTrdHit->GetELoss();
4197 }
4198 }
4199 }
4200 //2.b - monitor TRD dEdx
4201 for (Int_t i = 0; i < nTofHits; i++) {
4202 Int_t iMul = iTRD[i].size();
4203 if (iMul > 0) {
4204 Double_t ddEdx = 0.;
4205 for (Int_t l = 0; l < iMul; l++) {
4206 CbmTrdHit* pTrdHit = (CbmTrdHit*) fTrdHits->At(iTRD[i][l]);
4207 ddEdx += pTrdHit->GetELoss();
4208 }
4209 ddEdx /= (Double_t) iMul;
4210 fhdEdxMul->Fill((Double_t) iMul, ddEdx);
4211 }
4212 }
4213 }
4214
4215 //3. find secondary proton candidate
4216 for (Int_t i = 0; i < nTofHits; i++) {
4217 CbmTofHit* pTofHit = (CbmTofHit*) fTofHits->At(i);
4218 if (NULL == pTofHit) continue;
4219 if (pTofHit->GetZ() == 0.) continue; // skip fake beam counter
4220 if (iStsMin[i][0] > -1 && iStsMin[i][1] > -1) {
4221 CbmStsHit* pStsHit = (CbmStsHit*) fStsHits->At(iStsMin[i][0]);
4222 CbmStsHit* pSts2Hit = (CbmStsHit*) fStsHits->At(iStsMin[i][1]);
4223
4224 Double_t dDx = pStsHit->GetX() - pSts2Hit->GetX();
4225 Double_t dDy = pStsHit->GetY() - pSts2Hit->GetY();
4226 Double_t dDz = pStsHit->GetZ() - pSts2Hit->GetZ();
4227 // Extrapolate to z=0 (from Si information)
4228 Double_t dX0 = pSts2Hit->GetX() - dDx / dDz * pSts2Hit->GetZ();
4229 Double_t dY0 = pSts2Hit->GetY() - dDy / dDz * pSts2Hit->GetZ();
4230 Double_t dD0 = TMath::Sqrt(dX0 * dX0 + dY0 * dY0);
4231 fhD0prim->Fill(dD0);
4232
4233 if (dD0 > fdD0ProtLim
4234 && (Double_t) iTRD[i].size() >= fdTRDHmulMin) { // secondary proton canditate, memorize relevant quantities
4235 Double_t dDd = TMath::Sqrt(dDx * dDx + dDy * dDy + dDz * dDz);
4236 Double_t vel = pTofHit->GetR() / pTofHit->GetTime();
4237 Double_t bet = vel / clight;
4238 if (bet > 0.9999) continue; // bet=0.9999;
4239 Double_t m = secMass[1]; // assume proton
4240 Double_t pmag = m * bet / TMath::Sqrt(1. - bet * bet); // natural units
4241 Double_t pz = pmag * dDz / dDd;
4242 Double_t px = pmag * dDx / dDd;
4243 Double_t py = pmag * dDy / dDd;
4244 Double_t E = TMath::Sqrt(pmag * pmag + m * m);
4245 P[i].SetPxPyPzE(px, py, pz, E);
4246 X[i].SetXYZT(pTofHit->GetX(), pTofHit->GetY(), pTofHit->GetZ(), pTofHit->GetTime());
4247 LOG(debug) << "Init proton LV at ind " << Form("%d %d %d", i, iStsMin[i][0], iStsMin[i][1])
4248 << " with beta = " << bet << ", minv = " << P[i].M() << ", tof " << X[i].T() << ", TRDHmul "
4249 << iTRD[i].size();
4250 X0[i].SetXYZ(pSts2Hit->GetX(), pSts2Hit->GetY(), pSts2Hit->GetZ());
4251 DX[i].SetXYZ(dDx, dDy, dDz);
4252 proton_cand++;
4253 } // secondary proton canditate,
4254 }
4255 }
4256
4257 //4. secondary pion candidate
4258 for (Int_t i = 0; i < nTofHits; i++) {
4259 CbmTofHit* pTofHit = (CbmTofHit*) fTofHits->At(i);
4260 if (NULL == pTofHit) continue;
4261 if (pTofHit->GetZ() == 0.) continue; // skip fake beam counter
4262 LOG(debug) << "Tof " << i << Form(" sec cand Min %6.3f > %6.3f ?", dStsDistMin[i], fdDistPrimLim);
4263 if (dStsDistMin[i] > fdDistPrimLim) { // Tof hit not in the primary class
4264 Double_t dDistMin = 100.;
4265 Int_t jbest = -1;
4266 Int_t kbest = -1;
4267 for (Int_t j = 0; j < nStsHits; j++) {
4268 LOG(debug) << "Tof " << i << ", Sts " << j
4269 << Form(" ? sec cand %6.3f Min %6.3f ", dTofDistMin[j], fdDistPrimLim);
4270 if (dTofDistMin[j] > fdDistPrimLim) { // Sts hit not in the primary class
4271 CbmStsHit* pStsHit = (CbmStsHit*) fStsHits->At(j);
4272 // check for extension of pair to 2nd silicon plane
4273 Double_t dDx = pTofHit->GetX() - pStsHit->GetX();
4274 Double_t dDy = pTofHit->GetY() - pStsHit->GetY();
4275 Double_t dDz = pTofHit->GetZ() - pStsHit->GetZ();
4276
4277 for (Int_t k = 0; k < nStsHits; k++) {
4278 if (j == k) continue;
4279 CbmStsHit* pSts2Hit = (CbmStsHit*) fStsHits->At(k);
4280 Double_t sPos2Z = pSts2Hit->GetZ();
4281 Double_t sPos2Xext = pStsHit->GetX() + dDx / dDz * (sPos2Z - pStsHit->GetZ());
4282 Double_t sPos2Yext = pStsHit->GetY() + dDy / dDz * (sPos2Z - pStsHit->GetZ());
4283 Double_t dDist2 =
4284 TMath::Power(pSts2Hit->GetX() - sPos2Xext, 2) + TMath::Power(pSts2Hit->GetY() - sPos2Yext, 2);
4285 Double_t dDist = TMath::Sqrt(dDist2);
4286 LOG(debug) << "Sec Tof " << i << ", Sts " << j
4287 << Form(" Sts2 %d -> dist %6.3f < %6.3f ? at z = %4.1f", k, dDist,
4288 TMath::Min(dDistMin, fdDistSecLim2), sPos2Z);
4289 if (dDist < fdDistSecLim2 && dDist < dDistMin) { // secondary or pion candidate
4290 dDistMin = dDist;
4291 jbest = j;
4292 kbest = k;
4293 }
4294 } // for (Int_t k=0; k<nStsHits; k++) {
4295 } //if( dTofDistMin[j] > dDistPrimLim) { // Sts hit not in the primary class
4296 } // for (Int_t j=0; j<nStsHits; j++) {
4297 fhDperpS->Fill(dDistMin);
4298 LOG(debug) << "Sec Dist for TofHit " << i << ": " << dDistMin << ", j " << jbest << ", k " << kbest;
4299
4300 if (dDistMin < 100.) { // secondary candidate found, store vectors
4301 CbmStsHit* pStsHit = (CbmStsHit*) fStsHits->At(jbest);
4302 CbmStsHit* pSts2Hit = (CbmStsHit*) fStsHits->At(kbest);
4303 if (pSts2Hit->GetZ() > pStsHit->GetZ()) { // swap order
4304 pStsHit = (CbmStsHit*) fStsHits->At(kbest);
4305 pSts2Hit = (CbmStsHit*) fStsHits->At(jbest);
4306 }
4307 //4.a - confirm secondary track hypothesis by TRD
4308 if (nTrdHits > 0) {
4309 Double_t dDx = pTofHit->GetX() - pStsHit->GetX();
4310 Double_t dDy = pTofHit->GetY() - pStsHit->GetY();
4311 Double_t dDz = pTofHit->GetZ() - pStsHit->GetZ();
4312 for (Int_t l = 0; l < nTrdHits; l++) {
4313 CbmTrdHit* pTrdHit = (CbmTrdHit*) fTrdHits->At(l);
4314 // calculate expected position in Trd layer
4315 Double_t dXexp = pStsHit->GetX() + dDx * (pTrdHit->GetZ() - pStsHit->GetZ()) / dDz;
4316 Double_t dYexp = pStsHit->GetY() + dDy * (pTrdHit->GetZ() - pStsHit->GetZ()) / dDz;
4317 Double_t dDtrans =
4318 TMath::Sqrt(TMath::Power(dXexp - pTrdHit->GetX(), 2) + TMath::Power(dYexp - pTrdHit->GetY(), 2));
4319 UInt_t iTrdLayer = CbmTrdAddress::GetLayerId(pTrdHit->GetAddress());
4320 fhDTRDsec->Fill(dDtrans);
4321 LOG(debug) << "Inspect sec. TRD hit " << l << " in "
4322 << Form("Module 0x%08x, layer %d", pTrdHit->GetAddress(),
4324 << " at z= " << pTrdHit->GetZ() << " dD = " << dDtrans << " < " << fdDistTRD;
4325 if (dDtrans < fdDistTRD
4326 && dDtrans < dTrdDistMin[i][iTrdLayer]) { // check if acceptable and take best match
4327 Int_t iMul = iTRD[i].size();
4328 if (dTrdDistMin[i][iTrdLayer] < 1.E3) { // modify previous entry
4329 //find old entry in vector
4330 Int_t ll = 0;
4331 for (; ll < iMul; ll++)
4332 if (CbmTrdAddress::GetLayerId(((CbmTrdHit*) fTrdHits->At(iTRD[i][ll]))->GetAddress()) == iTrdLayer)
4333 break;
4334 iTRD[i][ll] = l;
4335 }
4336 else { //add hit
4337 dTrdDistMin[i][iTrdLayer] = dDtrans;
4338 iTRD[i].resize(iMul + 1);
4339 iTRD[i][iMul] = l;
4340 }
4341 LOG(debug) << "assign TrdHit " << l << " to TofHit " << i << " in layer " << iTrdLayer
4342 << " with d = " << dDtrans << ", TrdMul" << iMul << ", dEdx = " << pTrdHit->GetELoss();
4343 }
4344 }
4345 //4.b - monitor TRD dEdx
4346 Int_t iMul = iTRD[i].size();
4347 if (iMul > 0) {
4348 Double_t ddEdx = 0.;
4349 for (Int_t l = 0; l < iMul; l++) {
4350 CbmTrdHit* pTrdHit = (CbmTrdHit*) fTrdHits->At(iTRD[i][l]);
4351 ddEdx += pTrdHit->GetELoss();
4352 }
4353 ddEdx /= (Double_t) iMul;
4354 fhdEdxMulsec->Fill((Double_t) iMul, ddEdx);
4355 }
4356 } // (nTrdHits > 0) end
4357 if ((Double_t) iTRD[i].size() >= fdTRDHmulMin) {
4358 Double_t dDx = pStsHit->GetX() - pSts2Hit->GetX();
4359 Double_t dDy = pStsHit->GetY() - pSts2Hit->GetY();
4360 Double_t dDz = pStsHit->GetZ() - pSts2Hit->GetZ();
4361 Double_t dDd = TMath::Sqrt(dDx * dDx + dDy * dDy + dDz * dDz);
4362 Double_t vel = pTofHit->GetR() / pTofHit->GetTime(); // approximation, ignoring decay kinematics
4363 Double_t bet = vel / clight;
4364 if (bet > 0.9999) continue; //bet=0.9999;
4365 Double_t m = secMass[0]; // assume pion
4366 Double_t pmag = m * bet / TMath::Sqrt(1. - bet * bet); // natural units
4367 Double_t pz = pmag * dDz / dDd;
4368 Double_t px = pmag * dDx / dDd;
4369 Double_t py = pmag * dDy / dDd;
4370 Double_t E = TMath::Sqrt(pmag * pmag + m * m);
4371 P[i].SetPxPyPzE(px, py, pz, E);
4372 X[i].SetXYZT(pTofHit->GetX(), pTofHit->GetY(), pTofHit->GetZ(), pTofHit->GetTime());
4373 LOG(debug) << "Init pion LV at ind " << i << " with beta = " << bet << ", minv = " << P[i].M() << ", tof "
4374 << X[i].T() << ", TRDHmul " << iTRD[i].size();
4375 X0[i].SetXYZ(pSts2Hit->GetX(), pSts2Hit->GetY(), pSts2Hit->GetZ());
4376 DX[i].SetXYZ(dDx, dDy, dDz);
4377 pion_cand++;
4378 }
4379 }
4380 } //if( dStsDistMin[i] > dDistPrimLim) { // Sts hit not in the primary class
4381 } //for (Int_t i=0; i<nTofHits; i++)
4382 LOG(debug) << " Ev " << iCandEv << " has " << proton_cand << " protons and " << pion_cand << " pion candidates";
4383 if (proton_cand > 0 && pion_cand > 0) {
4384 LOG(debug) << "add event " << iCandEv << " to mixing class " << iMixClass << " of size " << fvP[iMixClass].size();
4385
4386 fvP[iMixClass].push_front(P); //insert to mixed event vector
4387 fvX[iMixClass].push_front(X); //insert to mixed event vector
4388 fvX0[iMixClass].push_front(X0); //insert to mixed event vector
4389 fvDX[iMixClass].push_front(DX); //insert to mixed event vector
4390
4391 if (fvP[iMixClass].size() > fNMixedEvents + 1) {
4392 fvP[iMixClass].pop_back();
4393 fvX[iMixClass].pop_back();
4394 fvX0[iMixClass].pop_back();
4395 fvDX[iMixClass].pop_back();
4396 }
4397 }
4398 else
4399 return; // nothing to be done
4400
4401 //5. do combinatorics
4402 for (Int_t i = 0; i < nTofHits; i++) {
4403 if (TMath::Abs(P[i].M() - secMass[1]) < 0.01) { // proton candidate
4404 std::list<std::vector<TLorentzVector>>::iterator itX = fvX[iMixClass].begin();
4405 std::list<std::vector<TVector3>>::iterator itX0 = fvX0[iMixClass].begin();
4406 std::list<std::vector<TVector3>>::iterator itDX = fvDX[iMixClass].begin();
4407 Int_t iMixEv = 0;
4408 itX0--;
4409 itDX--;
4410 itX--;
4411 LOG(debug) << "LV P has size " << P.size() << ", fvP size " << fvP[iMixClass].size() << " in mix class "
4412 << iMixClass;
4413 for (std::list<std::vector<TLorentzVector>>::iterator itP = fvP[iMixClass].begin(); itP != fvP[iMixClass].end();
4414 ++itP) {
4415 iMixEv++;
4416 ++itX;
4417 ++itX0;
4418 ++itDX;
4419 std::vector<TLorentzVector> PE = *itP; //fvP[iEv];
4420 std::vector<TLorentzVector> XE = *itX; //fvX[iEv];
4421 std::vector<TVector3> X0E = *itX0; //fvX0[iEv];
4422 std::vector<TVector3> DXE = *itDX; //fvDX[iEv];
4423 LOG(debug) << "iMixEv " << iMixEv << ": PE has size " << PE.size() << ", X0E: " << X0E.size();
4424 if (iMixEv == 1) {
4425 if (PE != P) {
4426 for (UInt_t j = 0; j < PE.size(); j++) {
4427 LOG(debug) << "P comp " << j << ": " << P[j].M() << " " << PE[j].M();
4428 }
4429 LOG(debug) << "P not properly restored from list ";
4430 }
4431 }
4432
4433 for (UInt_t j = 0; j < PE.size(); j++) {
4434 if (TMath::Abs(PE[j].M() - secMass[0]) < 0.01) { // pion candidate
4435 //request minimum opening angle
4436 Double_t dOpAngle = DX[i].Angle(DXE[j]);
4437 if (iMixEv == 1) fhOpAng->Fill(dOpAngle);
4438 else
4439 fhMIXOpAng->Fill(dOpAngle);
4440 if (dOpAngle < fdOpAngMin) continue;
4441 // calculate decay vertex
4442 TVector3 N = DX[i].Cross(DXE[j]);
4443 if (N.Mag() == 0.) continue;
4444 N.SetMag(1.);
4445 Double_t dDCA = TMath::Abs((X0[i] - X0E[j]) * N);
4446 if (iMixEv == 1) fhDCA->Fill(dDCA);
4447 else
4448 fhMIXDCA->Fill(dDCA);
4449 LOG(debug) << "DCA for iMixEv " << iMixEv << " at ind i " << i << ", j " << j << ": " << dDCA;
4450 if (dDCA == 0.) continue;
4451 if (dDCA < fdDCALim) {
4452 // vertex position
4453 TVector3 D = dDCA * N;
4454 TVector3 Ni = DX[i].Cross(N);
4455 Double_t cj = -(X0E[j] - X0[i] - D) * Ni / (DXE[j] * Ni);
4456 TVector3 V = X0E[j] + cj * DXE[j] - 0.5 * D; // Vertex vector
4457 Double_t dVLen = V.Mag(); // Lambda flight pathlength
4458 if (dVLen > fdVLenMax) continue;
4459 if (dVLen < fdVLenMin) continue;
4460 TLorentzVector PM = P[i] + PE[j];
4461 TVector3 PV = TVector3(PM.Px(), PM.Py(), PM.Pz());
4462 Double_t TofM = dVLen / PM.Beta() / clight; // flight time of M in ns
4463
4464 Double_t TofMLast = 100.;
4465 Int_t Niter = 0;
4466 std::vector<TLorentzVector> Ptmp = P;
4467 std::vector<TLorentzVector> PEtmp = PE;
4468 std::vector<TLorentzVector> Xtmp = X;
4469 std::vector<TLorentzVector> XEtmp = XE;
4470
4471 while (TMath::Abs(TofM - TofMLast) > 0.001 && Niter++ < 3) {
4472 LOG(debug) << "MinvI at ind i " << i << ", j " << j << ": " << PM.M() << ", vertex: " << V[0] << " "
4473 << V[1] << " " << V[2] << ", Len " << dVLen << ", mom = " << PV.Mag() << ", TofM " << TofM;
4474
4475 // update momentum calculation
4476 for (Int_t ii = 0; ii < 2; ii++) {
4477 Int_t k;
4478 TVector3 TofX;
4479 Double_t m;
4480 Double_t tof;
4481 if (ii == 0) {
4482 k = i;
4483 TofX = Xtmp[k].Vect();
4484 m = Ptmp[k].M();
4485 tof = Xtmp[k].T();
4486 }
4487 else {
4488 k = j;
4489 TofX = XEtmp[k].Vect();
4490 m = PEtmp[k].M();
4491 tof = XEtmp[k].T();
4492 }
4493 TVector3 vDTofV = TofX - V;
4494 Double_t vel = vDTofV.Mag() / (tof - TofM);
4495 Double_t bet = vel / clight;
4496 if (bet > 0.9999) bet = 0.9999;
4497 Double_t pmag = m * bet / TMath::Sqrt(1. - bet * bet); // natural units
4498 TVector3 vPsec = vDTofV;
4499 vPsec.SetMag(pmag);
4500 Double_t E = TMath::Sqrt(pmag * pmag + m * m);
4501 if (ii == 0) {
4502 Ptmp[k].SetVect(vPsec);
4503 Ptmp[k].SetE(E);
4504 }
4505 else {
4506 PEtmp[k].SetVect(vPsec);
4507 PEtmp[k].SetE(E);
4508 }
4509 }
4510 PM = Ptmp[i] + PEtmp[j];
4511 PV = TVector3(PM.Px(), PM.Py(), PM.Pz());
4512 TofMLast = TofM;
4513 TofM = dVLen / PM.Beta() / clight; // flight time of M in ns
4514 }
4515 Double_t minv = PM.M();
4516 if (iMixEv == 1) { // fill correlated event distributions
4517 fhMinv->Fill(minv);
4518 fhPathLen->Fill(dVLen);
4519 fhMMom->Fill(PV.Mag());
4520
4521 LOG(debug) << "MinvII in event " << fEvents << " at ind i " << i << ", j " << j << ": " << minv
4522 << ", vertex: " << V[0] << " " << V[1] << " " << V[2] << ", Len " << dVLen
4523 << ", mom = " << PV.Mag() << ", tof " << TofM;
4524
4525 if (TMath::Abs(MLAM - minv) < DMLAM) {
4526 PM.RotateY(beamRotY * TMath::Pi() / 180.);
4527 fa_ptm_rap_rec_lam->Fill(PM.Rapidity(), PM.Pt() / MLAM);
4528 }
4529 }
4530 else { // fill mixed event distributions
4531 fhMIXMinv->Fill(minv);
4532 fhMIXPathLen->Fill(dVLen);
4533 fhMIXMMom->Fill(PV.Mag());
4534 if (TMath::Abs(MLAM - minv) < DMLAM) {
4535 PM.RotateY(beamRotY * TMath::Pi() / 180.);
4536 fa_ptm_rap_mix_lam->Fill(PM.Rapidity(), PM.Pt() / MLAM);
4537 }
4538 }
4539 }
4540 }
4541 }
4542 }
4543 }
4544 }
4545} //void CbmHadronAnalysis::ReconstructSecondaries()
4546
Base class for cluster objects.
@ kTof
Time-of-flight Detector.
@ kSts
Silicon Tracking System.
static FairRootManager * rootMgr
static Int_t iNbTs
Int_t nGlobTracks
static TFile * fHist
const Int_t fiTofHitMulMax
TClonesArray * fTofHitsColl
Int_t nTofHits
#define MinWallDist
CbmDigiManager * fDigiMan
ClassImp(CbmHadronAnalysis)
Float_t refMass[3]
static TH1F * fhStsHitMul
static TH1F * fhTofHitMul
TClonesArray * fStsHitsColl
Int_t nStsHits
TClonesArray * fEventsColl
static Int_t iCandEv
Int_t nTofPoints
Int_t NMASS
#define M2PROT
Int_t nMCTracks
#define clight
Int_t nTofTracks
TClonesArray * fTofTracks
Data class for STS clusters.
Data class for a reconstructed hit in the STS.
Data class for STS tracks.
CbmDigiManager * fDigiMan
Class for hits in TRD detector.
friend fvec sqrt(const fvec &a)
static constexpr size_t size()
Definition KfSimdPseudo.h:2
Generates beam ions for transport simulation.
int32_t GetDigi(int32_t index) const
Get digi at position index.
Definition CbmCluster.h:76
int32_t GetNofDigis() const
Number of digis in cluster.
Definition CbmCluster.h:69
CbmDigiManager.
static Bool_t IsPresent(ECbmModuleId systemId)
Presence of a digi branch.
static Bool_t IsMatchPresent(ECbmModuleId systemId)
Presence of a digi match branch.
InitStatus Init()
Initialisation.
static CbmDigiManager * Instance()
Static instance.
const CbmMatch * GetMatch(ECbmModuleId systemId, UInt_t index) const
Get a match object.
Class characterising one event by a collection of links (indices) to data objects,...
Definition CbmEvent.h:34
size_t GetNofData() const
Definition CbmEvent.cxx:53
uint32_t GetIndex(ECbmDataType type, uint32_t iData)
Definition CbmEvent.cxx:42
std::string ToString() const
Definition CbmEvent.cxx:96
void SetParamLast(const FairTrackParam *parLast)
int32_t GetStsTrackIndex() const
void SetTofHitIndex(int32_t iTofHit)
int32_t GetTofHitIndex() const
void SetParamFirst(const FairTrackParam *parFirst)
void SetLength(double length)
const FairTrackParam * GetParamFirst() const
double GetLength() const
virtual InitStatus Init()
void SetMidY(Float_t midY)
TClonesArray * fTofDigiMatchColl
TClonesArray * fTrdPoints
TClonesArray * fTofHits
Float_t GetDY() const
TClonesArray * fGlobalTracks
TClonesArray * fHadrons
Float_t GetBSelMax() const
TClonesArray * fMCTracksColl
Float_t GetBSelMin() const
Float_t GetBeamMomentum() const
Float_t GetFlowMinPtm() const
CbmMCDataArray * fStsPoints
CbmStsKFTrackFitter fTrackFitter
FairMCEventHeader * fMCEventHeader
TClonesArray * fStsClusters
TClonesArray * fStsPointsColl
TClonesArray * fStsHits
TClonesArray * fTrdHits
CbmMCDataArray * fMCTracks
virtual void Exec(Option_t *option)
Float_t GetMidY() const
virtual void ExecEvent(Option_t *option)
TClonesArray * fTofPoints
TClonesArray * fStsTracks
double GetTime() const
Definition CbmHit.h:76
int32_t GetAddress() const
Definition CbmHit.h:74
double GetZ() const
Definition CbmHit.h:71
void SetTime(double time)
Definition CbmHit.h:85
Task class creating and managing CbmMCDataArray objects.
CbmMCDataObject * GetObject(const char *name)
CbmMCDataArray * InitBranch(const char *name)
double GetPz() const
Definition CbmMCTrack.h:72
double GetPt() const
Definition CbmMCTrack.h:97
double GetP() const
Definition CbmMCTrack.h:98
double GetPx() const
Definition CbmMCTrack.h:70
int32_t GetMotherId() const
Definition CbmMCTrack.h:69
double GetMass() const
Mass of the associated particle.
int32_t GetPdgCode() const
Definition CbmMCTrack.h:68
double GetPy() const
Definition CbmMCTrack.h:71
double GetRapidity() const
void GetStartVertex(TVector3 &vertex) const
Definition CbmMCTrack.h:176
double GetEnergy() const
Definition CbmMCTrack.h:162
const CbmLink & GetLink(int32_t i) const
Definition CbmMatch.h:39
int32_t GetNofLinks() const
Definition CbmMatch.h:42
const CbmLink & GetMatchedLink() const
Definition CbmMatch.h:41
double GetY() const
Definition CbmPixelHit.h:74
double GetX() const
Definition CbmPixelHit.h:73
Data class for STS clusters.
data class for a reconstructed 3-d hit in the STS
Definition CbmStsHit.h:35
int32_t GetFrontClusterId() const
Definition CbmStsHit.h:105
int32_t GetBackClusterId() const
Definition CbmStsHit.h:65
Double_t FitToVertex(CbmStsTrack *track, CbmVertex *vtx, FairTrackParam *v_track)
Double_t GetChiToVertex(CbmStsTrack *track, CbmVertex *vtx=nullptr)
int32_t GetStsHitIndex(int32_t iHit) const
int32_t GetNofStsHits() const
Definition CbmStsTrack.h:93
void SetFlag(int32_t flag)
Definition CbmTofHit.h:92
double GetR() const
Definition CbmTofHit.h:78
int32_t GetFlag() const
Definition CbmTofHit.h:75
Geometric intersection of a MC track with a TOFb detector.
Definition CbmTofPoint.h:44
float GetMass() const
Definition CbmTofTrack.h:52
int32_t GetTofHitIndex() const
Definition CbmTofTrack.h:55
double GetTrackLength() const
Definition CbmTofTrack.h:67
double GetDistance() const
Definition CbmTofTrack.h:79
void SetMass(float mass)
int32_t GetTrackIndex() const
Definition CbmTofTrack.h:64
const FairTrackParam * GetTrackParameter() const
Definition CbmTofTrack.h:70
int32_t GetPidHypo() const
Definition CbmTrack.h:61
const FairTrackParam * GetParamFirst() const
Definition CbmTrack.h:68
static uint32_t GetLayerId(uint32_t address)
Return layer ID from address.
data class for a reconstructed Energy-4D measurement in the TRD
Definition CbmTrdHit.h:40
double GetELoss() const
Definition CbmTrdHit.h:79
double GetZ() const
Definition CbmVertex.h:69
Hash for CbmL1LinkKey.