BmnRoot 26.10.0-18-g4e3414f9
Software framework of the BM@N experiment
Loading...
Searching...
No Matches
BmnMCEfficiencyTOF700.C
Go to the documentation of this file.
1#include "TParticle.h"
2
3#include <Rtypes.h>
4#include <fstream>
5#ifndef ROOT_TParticlePDG
6#include "TParticlePDG.h"
7#endif
8#ifndef ROOT_TDatabasePDG
9#include "TDatabasePDG.h"
10#endif
11R__ADD_INCLUDE_PATH($VMCWORKDIR)
12#include "bmndata/BmnEnums.h"
13// #include "macro/run/bmnloadlibs.C"
14
15// #define TOFGEOMETRY "TOF700_geometry_run7_1st.txt"
16#define TOFGEOMETRY "TOF700_geometry_run7.txt"
17
18#define USETRUEHIT 0
19#define SORT 0
20#define AARGON 40
21#define MCMATCH 0
22
23#define MINTRUE 0.5f
24#define GEMCUTS 0
25#define PRIMARY_VERTEX_CUTS 0
26#define MINGEMNDF1 6
27#define MINGEMNDF 4
28#define MINSIHITS 2
29#define GEMDEBUG false
30#define MCDEBUG false
31#define DXGEMDCHMAX 5.
32#define DYGEMDCHMAX 3.
33// #define DXGEMDCHMAX 2.5
34// #define DYGEMDCHMAX 1.5
35
36// CFACTOR - beta factor for main time peak particles
37#define CFACTOR 0.97f
38
39// BETAMAX max velocity for proton selection
40#define BETAMAX 1.00f
41
42// PMIN min momentum for track selection
43#define PMIN 0.5f
44
45// PMAX max momentum for track selection
46#define PMAX 5.0f
47
48// PMAX1 max momentum for proton mass limits for calculate corrections
49#define PMAX1 3.0f
50
51// PMAX2 max momentum for proton mass limits for resolution hists
52#define PMAX2 2.5f
53
54// PRES min momentum for resolution hists
55#define PRES 2.0f
56
57// PRES1 another min momentum for resolution hists
58#define PRES1 3.0f
59
60// PMAXPI max momentum for pion track selection for calculate corrections
61#define PMAXPI 3.0f
62
63// MPIMAX max limit pion mass
64#define MPIMAX 0.2f
65
66// MP proton mass
67#define MP 0.938f
68
69// MP2 proton mass sqare
70#define MP2 (MP * MP)
71
72// pion mass
73#define MPI 0.139
74
75// pion mass square
76#define MPI2 (MPI * MPI)
77
78// scale factor for min mass proton region selection
79#define DMPSCALE1 1.
80
81// scale factor for max mass proton region selection
82#define DMPSCALE2 1.
83
84// scale factor for min mass proton region selection for time hist
85#define DMPSCALET1 1.
86
87// scale factor for max mass proton region selection for time hist
88#define DMPSCALET2 1.5
89
90// scale factor for pion region selection
91#define DMPISCALE 1.
92
93#define XYPANIN 1
94#define XFACTOR 1.0f
95#define DISTMAX 20.0f
96#define DXMAX 5.0f
97#define DYMAX 3.0f
98#define DXMAXPAN 5.0f
99#define DYMAXPAN 5.0f
100
101#define ACCEPTANCE false
102#define DXBORDER 5.
103#define DYBORDER 5.
104
105#define TOFMAX 60.0f
106
107#define NUMBER_CHAMBERS 59
108
109Float_t halfwidthx[NUMBER_CHAMBERS] = {8.0};
110Float_t halfwidthy[NUMBER_CHAMBERS] = {0.5};
111
112int nhi = 0;
113int nhit[NUMBER_CHAMBERS] = {0};
114int nhits[NUMBER_CHAMBERS][32] = {{0}};
115float tof[NUMBER_CHAMBERS][32] = {{0.}};
116double xhi[1000], yhi[1000], zhi[1000];
117double xhic[1000], yhic[1000], zhic[1000];
118using namespace std;
119using namespace TMath;
120
121BmnTof2Raw2DigitNew* TOF2 = NULL;
122
123float idchambers[59] = {27.1, 28.1, 3.1, 1.1, 29.1, 4.1, 33.1, 30.1, 5.1, 19.3, 31.1, 6.1,
124 2.1, 32.1, 15.2, 16.2, 17.2, 18.2, 19.2, 20.2, 7.1, 115.2, 113.1, 117.1,
125 35.1, 9.1, 37.1, 11.1, 39.1, 13.1, 34.1, 8.1, 36.1, 10.1, 38.1, 12.1,
126 21.2, 23.2, 25.2, 22.2, 24.2, 26.2, 107.2, 108.2, 109.2, 110.2, 111.2, 112.2,
127 114.1, 116.2, 118.1, 14.1, 40.1, 119.2, 120.2, 121.2, 122.2, 123.2, 124.2};
128
129// float idchambers[59] = {0.};
130int order[59] = {-1};
131
132Int_t readMCTimeFile(const char* MCTimeFile = "TOF700_MC_argon_qgsm_time_run7.txt");
133
134Double_t tofcalmc[NUMBER_CHAMBERS] = {0.f};
135
136// #define TRACKFILE "/eos/nica/bmn/users/kovachev/efficiency/TracksForYuri/MC_From_Gleb/recoMC_GEM_TOF.root"
137// #define MCFILE "/eos/nica/bmn/users/kovachev/efficiency/TracksForYuri/MC_From_Gleb/mc_GEM_CSC_TOF1.root"
138
139// #define TRACKFILE "/eos/nica/bmn/users/petukhov/MC/recoMC_GEM_TOF.root"
140// #define MCFILE "/eos/nica/bmn/users/petukhov/MC/mc_GEM_CSC_TOF1.root"
141
142// #define TRACKFILE "/eos/nica/bmn/users/petukhov/MC/reco_GEM_CSC_DCH_Tof_Ar_all_targets.root"
143
144// #define TRACKFILE "/eos/nica/bmn/users/petukhov/MC/EXTR510_GEM_CSC_DCH_Tof_Ar_all_targets.root"
145// #define MCFILE "/eos/nica/bmn/users/petukhov/MC/mc_GEM_CSC_DCH_Tof_Ar.root"
146
147// #define TRACKFILE "/nica/mpd15/kovachev/MC_FILES/EXTR510_DCMQGSM_Ar_All_3.2AGeV_mb_250k.root"
148// #define MCFILE "/nica/mpd15/kovachev/MC_FILES/mc_DCMQGSM_Ar_All_3.2AGeV_mb_250k.root"
149
150/*
151/nica/mpd16/kovachev/MC_run7_efficiency/MC_Files_run7_For_Yuri/EXTR510_DQGSM_Ar_old_target.root
152/nica/mpd16/kovachev/MC_run7_efficiency/MC_Files_run7_For_Yuri/reco_DCMQGSM_Ar_Old_target_1-23.root
153*/
154
155// #define MCFILE "/nica/mpd16/kovachev/MC_run7_efficiency/MC_Files_run7_For_Yuri/sim_DCMQGSM_Ar_1-23_Old_target.root"
156// #define TRACKFILE
157// "/nica/mpd16/kovachev/MC_run7_efficiency/MC_Files_run7_For_Yuri/GemDch_Link_DCMQGSM_old_target.root"
158
159// #define MCFILE "/nica/mpd16/kovachev/MC_run7_efficiency/MC_Files_run7_For_Yuri/sim_DCMQGSM_Ar_1-23_Old_target.root"
160// #define TRACKFILE
161// "/nica/mpd16/kovachev/MC_run7_efficiency/MC_Files_run7_For_Yuri/GemDch_Link_DCMQGSM_old_target_efficiency.root"
162
163// #define TRACKFILE
164// "/nica/mpd16/kovachev/MC_run7_efficiency/MC_Files_run7_For_Yuri/GemDch_Link_DCMQGSM_Ar_New_target_without_SetBeam_AR.root"
165// #define TRACKFILE "mk1_argon_files.txt"
166// #define TRACKFILE "mk3_argon_files.txt"
167#define TRACKFILE "/eos/nica/bmn/users/kovachev/test_Urqmd/Gem_Dch_Tracks_Urqmd2-80.root"
168// #define TRACKFILE "/eos/nica/bmn/users/kovachev/test_Urqmd/Gem_Dch_Tracks_Urqmd.root"
169#define MCFILE \
170 "/nica/mpd16/kovachev/MC_run7_efficiency/MC_Files_run7_For_Yuri/" \
171 "sim_DCMQGSM_New_target_without_SetBeam_Ar_1-25_.root"
172
173void BmnMCEfficiencyTOF700(const char* TracksGemDch = TRACKFILE, const char* DataMC = MCFILE)
174{
175
176 /* Load basic libraries */
177 // #if ROOT_VERSION_CODE < ROOT_VERSION(5,99,99)
178 // gROOT->LoadMacro("$VMCWORKDIR/macro/run/bmnloadlibs.C");
179 // #endif
180 // bmnloadlibs(); // load BmnRoot libraries
181
182 int RunPeriod = 7;
183 char mapping[256];
184 if (RunPeriod >= 1 && RunPeriod <= 7) {
185 sprintf(mapping, "TOF700_map_period_%d.txt", RunPeriod);
186 } else {
187 printf("Non-existing run period number %d!", RunPeriod);
188 return;
189 }
190
191 TClonesArray* eventHeaderMC = NULL;
192 TClonesArray *tof700Digits = NULL, *eventHeader = NULL, *mcmatch = NULL, *stspoints = NULL;
193 TClonesArray *tracks = NULL, *mtracks = NULL, *MCtracks = NULL, *tofpoints = NULL;
194 CbmVertex* primaryVertex = NULL;
195 TString filename;
196
197 char fnameMC[256];
198 // TChain *bmnTree = new TChain("bmndata");
199 TChain* bmnTree = new TChain("BMN_DIGIT");
200 if (strstr(TracksGemDch, ".root")) {
201 bmnTree->Add(TracksGemDch);
202 strcpy(fnameMC, TracksGemDch);
203 } else {
204 ifstream fin(TracksGemDch);
205 if (fin.is_open()) {
206 for (int idx = 0; idx < 1000; idx++) {
207 if (fin.eof())
208 break;
209 fin >> filename;
210 if (strlen(filename.Data()) == 0)
211 break;
212 if ((filename.Data())[0] == '#')
213 continue;
214 cout << " GEM+DCH tracks file #" << idx + 1 << " " << filename.Data() << endl;
215 bmnTree->Add(filename);
216 strcpy(fnameMC, filename.Data());
217 }
218 fin.close();
219 } else {
220 printf("Non-existing run list file: %s!", TracksGemDch);
221 return;
222 }
223 }
224 // bmnTree->SetBranchAddress("EventHeader", &eventHeader);
225 bmnTree->SetBranchAddress("GemDchTracks", &mtracks);
226 // bmnTree->SetBranchAddress("BmnTofHit", &tof700Digits);
227 bmnTree->SetBranchAddress("BmnTof700Hit", &tof700Digits);
228 if (MCMATCH)
229 bmnTree->SetBranchAddress("StsTrackMatch", &mcmatch);
230 bmnTree->SetBranchAddress("MCTrack", &MCtracks);
231 bmnTree->SetBranchAddress("StsPoint", &stspoints);
232 bmnTree->SetBranchAddress("TOF700Point", &tofpoints);
233 bmnTree->SetBranchAddress("PrimaryVertex.", &primaryVertex);
234 /*
235 char fnameMC[256];
236 TChain *bmnTreeMC = new TChain("bmndata");
237 if (strstr(DataMC,".root"))
238 {
239 bmnTreeMC->Add(DataMC); //
240 strcpy(fnameMC, DataMC);
241 }
242 else
243 {
244 ifstream fin(DataMC);
245 if (fin.is_open())
246 {
247 for(int idx=0; idx<1000; idx++) {
248 if (fin.eof()) break;
249 fin >> filename;
250 if (strlen(filename.Data())==0) break;
251 if ((filename.Data())[0]=='#') continue;
252 cout << " TOF Data file #" << idx+1 << " " << filename.Data() << endl;
253 bmnTreeMC->Add(filename);
254 strcpy(fnameMC, filename.Data());
255 }
256 fin.close();
257 }
258 else
259 {
260 printf("Non-existing run list file: %s!", DataMC);
261 return;
262 }
263 }
264 // bmnTreeMC->SetBranchAddress("EventHeader", &eventHeaderMC);
265 bmnTreeMC->SetBranchAddress("MCTrack", &MCtracks);
266 bmnTreeMC->SetBranchAddress("StsPoint", &stspoints);
267 bmnTreeMC->SetBranchAddress("TOF700Point", &tofpoints);
268 */
269 TOF2 = new BmnTof2Raw2DigitNew(mapping, "bmn_run4649_raw.root", 0, 0, TOFGEOMETRY);
270
271 readMCTimeFile();
272
273 Double_t xmintof700 = +1000., xmaxtof700 = -1000., ymintof700 = +1000., ymaxtof700 = -1000.;
274 Float_t xs = 0., ys = 0., zs = 0.;
275 Double_t xsmin = 0., ysmin = 0., xsmax = 0., ysmax = 0.;
276 for (int c = 0; c < NUMBER_CHAMBERS; c++) {
277 halfwidthy[c] = 0.5f;
278 if (idchambers[c] < 100.f) {
279 halfwidthx[c] = 8.0f;
280 } else {
281 halfwidthx[c] = 28.0f;
282 }
283 int nstr = 32;
284 if (idchambers[c] > 100.f)
285 nstr = 16;
286 TOF2->get_strip_xyz(c, 0, &xs, &ys, &zs);
287 xsmin = xs - halfwidthx[c];
288 xsmax = xs + halfwidthx[c];
289 ysmin = ys - halfwidthy[c];
290 ysmax = ys + halfwidthy[c];
291 if (xsmin < xmintof700)
292 xmintof700 = xsmin;
293 if (xsmax > xmaxtof700)
294 xmaxtof700 = xsmax;
295 if (ysmin < ymintof700)
296 ymintof700 = ysmin;
297 if (ysmax > ymaxtof700)
298 ymaxtof700 = ysmax;
299 TOF2->get_strip_xyz(c, nstr - 1, &xs, &ys, &zs);
300 xsmin = xs - halfwidthx[c];
301 xsmax = xs + halfwidthx[c];
302 ysmin = ys - halfwidthy[c];
303 ysmax = ys + halfwidthy[c];
304 if (xsmin < xmintof700)
305 xmintof700 = xsmin;
306 if (xsmax > xmaxtof700)
307 xmaxtof700 = xsmax;
308 if (ysmin < ymintof700)
309 ymintof700 = ysmin;
310 if (ysmax > ymaxtof700)
311 ymaxtof700 = ysmax;
312 }
313 xmintof700 -= DXBORDER;
314 xmaxtof700 += DXBORDER;
315 ymintof700 -= DYBORDER;
316 ymaxtof700 += DYBORDER;
317
318 char name[128], title[128];
319
320 gStyle->SetOptFit(111);
321 gStyle->SetOptStat(111111);
322
323 TFile* f = NULL;
324 if (strlen(fnameMC)) {
325 char fname_root[128], *delim = 0;
326 strcpy(fname_root, gSystem->BaseName(fnameMC));
327 if ((delim = strrchr(fname_root, (int)'.')))
328 *delim = '\0';
329 strcat(fname_root, "_MC_efficiency_results.root");
330 f = new TFile(fname_root, "RECREATE", "Results of BmnTOF700 MC");
331 }
332 TH1F* hMCmass0 = new TH1F("hMCmass0", "Mass/Z, Monte-Carlo, nominal", 200, -0.5, 9.5);
333 TH1F* hMCrmass0 = new TH1F("hMCrmass0", "Mass/Z, Monte-Carlo, nominal", 200, -0.5, 9.5);
334 TH1F* hMCamass = new TH1F("hMCamass", "Mass/Z, Monte-Carlo, nominal", 200, -0.5, 9.5);
335 TH1F* h1MCmass0 = new TH1F("h1MCmass0", "Mass/Z, Monte-Carlo, reconstructable", 200, -0.5, 9.5);
336 TH1F* h1MCrmass0 = new TH1F("h1MCrmass0", "Mass/Z, Monte-Carlo, reconstructable", 200, -0.5, 9.5);
337 TH1F* h1MCamass = new TH1F("h1MCamass", "Mass/Z, Monte-Carlo, reconstructable", 200, -0.5, 9.5);
338 TH1F* h1MCeffmass = new TH1F("h1MCeffmass", "Eff Mass/Z, Monte-Carlo, reconstructable", 200, -0.5, 9.5);
339 TH1F* h2MCamass = new TH1F("h2MCamass", "Mass/Z, Monte-Carlo, reconstructable, p<2", 200, -0.5, 9.5);
340 TH1F* h2MCeffmass = new TH1F("h2MCeffmass", "Eff Mass/Z, Monte-Carlo, reconstructable, p<2", 200, -0.5, 9.5);
341 TH1F* h3MCamass = new TH1F("h3MCamass", "Mass/Z, Monte-Carlo, NDF>5", 200, -0.5, 9.5);
342 TH1F* h3MCeffmass = new TH1F("h3MCeffmass", "Eff Mass/Z, Monte-Carlo, NDF>5", 200, -0.5, 9.5);
343 TH1F* h4MCamass = new TH1F("h4MCamass", "Mass/Z, Monte-Carlo, TOF700 only", 200, -0.5, 9.5);
344 TH1F* h4MCeffmass = new TH1F("h4MCeffmass", "Eff Mass/Z, Monte-Carlo, TOF700 only", 200, -0.5, 9.5);
345 TH1F* hMCmass2 = new TH1F("hMCmass2", "Mass2/Z2, Monte-Carlo, simple reconstion", 200, -0.5, 9.5);
346 TH1F* hMCmass = new TH1F("hMCmass", "Mass/Z, Monte-Carlo, simple reconstion", 200, -0.5, 9.5);
347 TH1F* hMCmassr = new TH1F("hMCmassr", "Mass/Z, Monte-Carlo, simple reconstion", 200, -0.5, 9.5);
348 TH1F* hMCmassrp = new TH1F("hMCmassrp", "Mass/Z, Monte-Carlo, simple reconstion, p<2", 200, -0.5, 9.5);
349 TH1F* hMCeffmass = new TH1F("hMCeffmass", "Eff Mass/Z, Monte-Carlo, simple reconstruction", 200, -0.5, 9.5);
350 TH1F* hMCeffmassp = new TH1F("hMCeffmassp", "Eff Mass/Z, Monte-Carlo, simple reconstruction, p<2", 200, -0.5, 9.5);
351 TH2F* hMCavsz = new TH2F("hMCavsz", "A vs Z, Monte-Carlo", 30, 0, 30, 30, 0, 30);
352 hMCavsz->SetOption("COLZ");
353 TH2F* hMCravsz = new TH2F("hMCravsz", "A vs Z, Monte-Carlo, reconstructable", 30, 0, 30, 30, 0, 30);
354 hMCravsz->SetOption("COLZ");
355 TH2F* h1MCavsz = new TH2F("h1MCavsz", "A vs Z, Monte-Carlo", 30, 0, 30, 30, 0, 30);
356 h1MCavsz->SetOption("COLZ");
357 TH2F* h1MCravsz = new TH2F("h1MCravsz", "A vs Z, Monte-Carlo, reconstructable", 30, 0, 30, 30, 0, 30);
358 h1MCravsz->SetOption("COLZ");
359 TH2F* hpbetamc = new TH2F("hpbetamc", "Beta vs momentum, Monte-Carlo", 200, 0, 10, 220, 0, 1.1f);
360 hpbetamc->SetOption("COLZ");
361 TProfile* hpanglep =
362 new TProfile("hpanglep", "Average Angle vs momentum, fragments, Monte-Carlo", 150, 0, 15, 0, 70.f, "e");
363 TH2F* hpangle = new TH2F("hpangle", "Angle vs momentum, fragments, Monte-Carlo", 200, -5, 15, 50, 0, 50.f);
364 hpangle->SetOption("COLZ");
365 TH1F* hangle = new TH1F("hangle", "Theta angle of tracks at vertex", 100, 0, 50);
366 TProfile* hpanglep1 = new TProfile(
367 "hpanglep1", "Average Angle vs momentum, fragments, Monte-Carlo, reconstructable", 150, 0, 15, 0, 70.f, "e");
368 TH2F* hpangle1 =
369 new TH2F("hpangle1", "Angle vs momentum, fragments, Monte-Carlo, reconstructable", 200, -5, 15, 50, 0, 50.f);
370 hpangle1->SetOption("COLZ");
371 TH1F* hangle1 = new TH1F("hangle1", "Theta angle of tracks at vertex, reconstructable", 100, 0, 50);
372 TProfile* hpanglep2 = new TProfile("hpanglep2", "Average Angle vs momentum, fragments, Monte-Carlo, TOF700 hit",
373 150, 0, 15, 0, 70.f, "e");
374 TH2F* hpangle2 =
375 new TH2F("hpangle2", "Angle vs momentum, fragments, Monte-Carlo, TOF700 hit", 200, -5, 15, 50, 0, 50.f);
376 hpangle2->SetOption("COLZ");
377 TH1F* hangle2 = new TH1F("hangle2", "Theta angle of tracks at vertex, TOF700 hit", 100, 0, 50);
378 TH1F* hangleff = new TH1F("hangleff", "Efficiency vs Theta angle of tracks at vertex", 100, 0, 50);
379 TH1F* hangleff2 = new TH1F("hangleff2", "Efficiency vs Theta angle of tracks at vertex, TOF700 only", 100, 0, 50);
380 TH1F* hMCpndf = new TH1F("hMCpndf", "Momentum of tracks", 1000, -5, 15);
381 TH1F* hMCpndf0 = new TH1F("hMCpndf0", "Momentum of tracks, NDF<6", 1000, -5, 15);
382 TH1F* hMCpndf1 = new TH1F("hMCpndf1", "Momentum of tracks, NDF>=6", 1000, -5, 15);
383 TH1F* hMCpndf2 = new TH1F("hMCpndf2", "Momentum of tracks, NDF>=6, NTOF>0", 1000, -5, 15);
384 // TH1F *hnhgem = new TH1F("hnhgem","MC GEM hits for fragments", 10,0,10);
385 // TH1F *hnhgem0 = new TH1F("hnhgem0","MC GEM hits for fragments (no TOF700 hits)", 10,0,10);
386 // TH1F *hnhgem1 = new TH1F("hnhgem1","MC GEM hits for fragments (with TOF700 hits)", 10,0,10);
387
388 TH1F* hnsts = new TH1F("hnsts", "Number of STS points", 20, 0, 20);
389 TH1F* hntr = new TH1F("hntr", "Number of all tracks", 200, 0, 200);
390 TH1F* hntr0 = new TH1F("hntr0", "Number of selected tracks", 200, 0, 200);
391 TH1F* hntr1 = new TH1F("hntr1", "Number of tracks with FLAG conditions", 200, 0, 200);
392 TH2F* hxydch = new TH2F("hxydch", "XY DCH tracks", 70, -140, 140, 60, -90, 90);
393 TH2F* hxygem = new TH2F("hxygem", "XY GEM tracks", 70, -140, 140, 60, -90, 90);
394 TH2F* hpxysts = new TH2F("hpxysts", "XY GEM points (P/Z > 7)", 140, -70, 70, 120, -60, 60);
395 TH2F* hpxygem = new TH2F("hpxygem", "XY GEM REC tracks (P/Z > 7)", 70, -140, 140, 60, -90, 90);
396 TH2F* hp0xygem = new TH2F("hp0xygem", "XY GEM MC tracks (P/Z > 7)", 70, -140, 140, 60, -90, 90);
397 TH1F* hptr = new TH1F("hptr", "Momentum of GEM tracks", 1000, -10, 10);
398 TH1F* hpgemdch = new TH1F("hpgemdch", "Momentum of GEM+DCH tracks", 1000, -10, 10);
399 TH1F* hpgoodhit = new TH1F("hpgoodhit", "Momentum of GEM+DCH+TOF700 tracks", 1000, -10, 10);
400 TH1F* hpbadhit = new TH1F("hpbadhit", "Momentum of GEM+DCH+NoTOF700 tracks", 1000, -10, 10);
401 TH1F* htruehits = new TH1F("htruehits", "True hits fraction in reconstructed GEM track", 110, 0, 1.1);
402 TH1F* hallhits = new TH1F("hallhits", "Number of hits in reconstructed GEM track", 20, 0, 20);
403 TH1F* hp0 = new TH1F("hp0", "Momentum of tracks (all)", 1000, -10, 10);
404 TH1F* hltr = new TH1F("hltr", "Lenght of GEM tracks", 100, 500, 700);
405 TH1F* hltra = new TH1F("hltra", "Lenght of GEM+DCH tracks", 100, 500, 700);
406 TH2F* hxytrtof = new TH2F("hxytrtof", "XY tracks with TOF700 hits", 1500, -150, 150, 1000, -100, 100);
407 hxytrtof->SetOption("COLZ");
408 TH2F* hxyhit = new TH2F("hxyhit", "XY track hit distance from TOF700 hit", 100, -50, 50, 40, -10, 10);
409 TH1F* hxdiff = new TH1F("hxdiff", "X difference for TOF700 hits with track", 100, -50, +50);
410 TH2F* hxxhits = new TH2F("hxxhits", "X track hit vs TOF700 hit (small chambers)", 200, -100, 100, 200, -100, 100);
411 TH2F* hxxhitb = new TH2F("hxxhitb", "X track hit vs TOF700 hit (big chambers)", 200, -100, 100, 200, -100, 100);
412
413 TH1F* htoftr = new TH1F("htoftr", "TOF of tracks with TOF700 hits", 600, -20, 40);
414 TH1F* htof = new TH1F("htof", "TOF for hits", 600, -20, 40);
415 TH1F* hntof = new TH1F("hntof", "Number of hits", 100, 0, 100);
416 TH1F* hnhit = new TH1F("hnhit", "Number of hits in width limits", 100, 0, 100);
417 TH2F* hxytof = new TH2F("hxytof", "XY TOF700 hits", 1500, -150, 150, 1000, -100, 100);
418 hxytof->SetOption("COLZ");
419 TH2F* hpbeta = new TH2F("hpbeta", "Beta vs momentum", 200, 0, 10, 220, 0, 1.1f);
420 hpbeta->SetOption("COLZ");
421 TH2F* hpbetac = new TH2F("hpbetac", "Beta vs momentum, chamber corrections", 200, 0, 10, 220, 0, 1.1f);
422 hpbetac->SetOption("COLZ");
423 TH2F* hpbetas = new TH2F("hpbetas", "Beta vs momentum, strip corrections", 200, 0, 10, 220, 0, 1.1f);
424 hpbetas->SetOption("COLZ");
425 TH2F* hpbetacm = new TH2F("hpbetacm", "Beta vs momentum, Z<0, chamber corrections", 200, 0, 10, 220, 0, 1.1f);
426 hpbetacm->SetOption("COLZ");
427 TH2F* hpbetasm = new TH2F("hpbetasm", "Beta vs momentum, Z<0, strip corrections", 200, 0, 10, 220, 0, 1.1f);
428 hpbetasm->SetOption("COLZ");
429 TH1F* hmass = new TH1F("hmass", "Mass, p<2", 200, -0.5, 9.5);
430 TH1F* hmassc = new TH1F("hmassc", "Mass, chamber corrections", 100, -0.5, 4.5);
431 TH1F* hmassc15 = new TH1F("hmassc15", "Mass, chamber corrections, p<1.5", 100, -0.5, 4.5);
432 TH1F* hmassc20 = new TH1F("hmassc20", "Mass, chamber corrections, p<2.0", 100, -0.5, 4.5);
433 TH1F* hmasss = new TH1F("hmasss", "Mass, strip corrections", 100, -0.5, 4.5);
434 TH1F* hmasss15 = new TH1F("hmasss15", "Mass, strip corrections, p<1.5", 100, -0.5, 4.5);
435 TH1F* hmasss20 = new TH1F("hmasss20", "Mass, strip corrections, p<2.0", 100, -0.5, 4.5);
436 TH1F* hmasscm = new TH1F("hmasscm", "Mass, Z<0, chamber corrections", 200, -1, 4);
437 TH1F* hmasssm = new TH1F("hmasssm", "Mass, Z<0, strip corrections", 200, -1, 4);
438
439 TH1F* hdtps = new TH1F("hdtps", "Time resolution, protons, strip corrections", 200, -1., +1.);
440 TH2F* hdtpsvsp =
441 new TH2F("hdtpsvsp", "Time resolution vs momentum, protons, strip corrections", 240, 0., 6., 200, -1., +1.);
442 hdtpsvsp->SetOption("COLZ");
443
444 TH1F* hdtpc = new TH1F("hdtpc", "Time resolution, protons, chamber corrections", 200, -1., +1.);
445 TH2F* hdtpcvsp =
446 new TH2F("hdtpcvsp", "Time resolution vs momentum, protons, chamber corrections", 240, 0., 6., 200, -1., +1.);
447 hdtpcvsp->SetOption("COLZ");
448 TH1F* hdtpcplim = new TH1F("hdtpcplim", "Time resolution, protons, p>3", 200, -1., +1.);
449 TH1F* hdtpcp1lim = new TH1F("hdtpcp1lim", "Time resolution, protons, p>4", 200, -1., +1.);
450 TH1F* hdtpcp2lim = new TH1F("hdtpcp2lim", "Time resolution, protons, p<2", 200, -1., +1.);
451 TH1F* hdtpcp3lim = new TH1F("hdtpcp3lim", "Time resolution, protons, 2<p<5", 200, -1., +1.);
452
453 TH2F* hmassvsp = new TH2F("hmassvsp", "Mass vs momentum", 240, 0., 6., 200, -0.2, 4.8);
454 hmassvsp->SetOption("COLZ");
455 TH2F* hmasssvsp = new TH2F("hmasssvsp", "Mass vs momentum, strip correction", 240, 0., 6., 200, -0.2, 4.8);
456 hmasssvsp->SetOption("COLZ");
457 TH2F* hmasscvsp = new TH2F("hmasscvsp", "Mass vs momentum, chamber corrections", 240, 0., 6., 200, -0.2, 4.8);
458 hmasscvsp->SetOption("COLZ");
459 TH2F* hmasssmvsp = new TH2F("hmasssmvsp", "Mass vs momentum, Z<0, strip correction", 240, 0., 6., 200, -0.2, 4.8);
460 hmasssmvsp->SetOption("COLZ");
461 TH2F* hmasscmvsp =
462 new TH2F("hmasscmvsp", "Mass vs momentum, Z<0, chamber corrections", 240, 0., 6., 200, -0.2, 4.8);
463 hmasscmvsp->SetOption("COLZ");
464
465 TH1F* htofp = new TH1F("htofp", "TOF for proton hits, p<2", 600, -20, 40);
466
467 TH1F* hstr = new TH1F("hstr", "Strips", 32, 0, 32);
468 TH1F* hcha = new TH1F("hcha", "Chambers", 60, 0, 60);
469 TH1F* hstrtr = new TH1F("hstrtr", "Strips with tracks", 32, 0, 32);
470 TH1F* hchatr = new TH1F("hchatr", "Chambers with tracks", 60, 0, 60);
471
472 TH2F* hdxy0 = new TH2F("hdxy0", "DxDy all", 40, -5, +5, 40, -5, +5);
473 hdxy0->SetOption("COLZ");
474 TH1F* hdx = new TH1F("hdx", "dx all", 40, -5, +5);
475 TH1F* hdy = new TH1F("hdy", "dy all", 40, -5, +5);
476 TH2F* hdxy = new TH2F("hdxy", "DxDy selected", 40, -5, +5, 40, -5, +5);
477 hdxy->SetOption("COLZ");
478 TH2F* hdxytof = new TH2F("hdxytof", "DxDy with TOF700 hits", 40, -5, +5, 40, -5, +5);
479 hdxytof->SetOption("COLZ");
480 TH2F* hdxyeff = new TH2F("hdxyeff", "Efficiency vs DxDy", 40, -5, +5, 40, -5, +5);
481 hdxyeff->SetOption("COLZ");
482 TH1F* hdxeff = new TH1F("hdxeff", "Efficiency vs Dx", 40, -5, +5);
483 TH1F* hdyeff = new TH1F("hdyeff", "Efficiency vs Dy", 40, -5, +5);
484 TH2F* hxygoodtr = new TH2F("hxygoodtr", "XY(dch) GEM+DCH good tracks", 70, -140, 140, 60, -90, 90);
485 hxygoodtr->SetOption("COLZ");
486 TH2F* hxygoodtrgem = new TH2F("hxygoodtrgem", "XY(gem) GEM+DCH good tracks", 70, -140, 140, 60, -90, 90);
487 hxygoodtrgem->SetOption("COLZ");
488 TH2F* hxygoodhit = new TH2F("hxygoodhit", "XY good tracks (GEM+DCH) with good hits", 70, -140, 140, 60, -90, 90);
489 hxygoodhit->SetOption("COLZ");
490 TH2F* hxygoodhita = new TH2F("hxygoodhita", "XY good tracks (GEM+DCH) with good hits", 70, -140, 140, 60, -90, 90);
491 hxygoodhita->SetOption("COLZ");
492 TH2F* hxybadhit = new TH2F("hxybadhit", "XY good tracks (GEM+DCH) with no hits", 70, -140, 140, 60, -90, 90);
493 hxybadhit->SetOption("COLZ");
494 TH2F* hxybadhita = new TH2F("hxybadhita", "XY good tracks (GEM+DCH) with no hits", 70, -140, 140, 60, -90, 90);
495 hxybadhita->SetOption("COLZ");
496 TH2F* hxygoodhit1 = new TH2F("hxygoodhit1", "XY good tracks (DCH) with good hits", 70, -140, 140, 60, -90, 90);
497 hxygoodhit1->SetOption("COLZ");
498 TH2F* hxygoodhit1a = new TH2F("hxygoodhit1a", "XY good tracks (DCH) with good hits", 70, -140, 140, 60, -90, 90);
499 hxygoodhit1a->SetOption("COLZ");
500 TH2F* hxyeff =
501 new TH2F("hxyeff", "Efficiency GEM+DCH+TOF700 vs XY good GEM+DCH tracks", 70, -140, 140, 60, -90, 90);
502 hxyeff->SetOption("COLZ");
503 TH2F* hxyeffa =
504 new TH2F("hxyeffa", "Efficiency GEM+DCH+TOF700 vs XY good GEM+DCH tracks", 70, -140, 140, 60, -90, 90);
505 hxyeffa->SetOption("COLZ");
506 TH1F* hxeff = new TH1F("hxeff", "Efficiency GEM+DCH+TOF700 vs X good GEM+DCH tracks", 70, -140, 140);
507 TH1F* hxeffa = new TH1F("hxeffa", "Efficiency GEM+DCH+TOF700 vs X good GEM+DCH tracks", 70, -140, 140);
508 TH1F* hyeff = new TH1F("hyeff", "Efficiency GEM+DCH+TOF700 vs Y good GEM+DCH tracks", 60, -90, 90);
509 TH1F* hyeffa = new TH1F("hyeffa", "Efficiency GEM+DCH+TOF700 vs Y good GEM+DCH tracks", 60, -90, 90);
510 TH2F* hxyeffmg = new TH2F("hxyeffmg", "Efficiency GEM+DCH vs XY good GEM tracks", 70, -140, 140, 60, -90, 90);
511 hxyeffmg->SetOption("COLZ");
512 TH2F* hxyeffmd = new TH2F("hxyeffmd", "Efficiency GEM+DCH vs XY good DCH tracks", 70, -140, 140, 60, -90, 90);
513 hxyeffmd->SetOption("COLZ");
514 TH2F* hxyeffm1d = new TH2F("hxyeffm1d", "Efficiency DCH+TOF700 vs XY good DCH tracks", 70, -140, 140, 60, -90, 90);
515 hxyeffm1d->SetOption("COLZ");
516 TH2F* hxyeffm1da =
517 new TH2F("hxyeffm1da", "Efficiency DCH+TOF700 vs XY good DCH tracks", 70, -140, 140, 60, -90, 90);
518 hxyeffm1da->SetOption("COLZ");
519 TH1F* hpeff = new TH1F("hpeff", "Efficiency GEM+DCH+TOF700 vs P good GEM+DCH tracks", 1000, -10, 10);
520 TH1F* hpeffmg = new TH1F("hpeffmg", "Efficiency GEM+DCH vs P good GEM tracks", 1000, -10, 10);
521 TH2F* hxyeffdg = new TH2F("hxyeffdg", "Efficiency DCH vs XY good GEM tracks", 70, -140, 140, 60, -90, 90);
522 hxyeffdg->SetOption("COLZ");
523 TH2F* hxyeffgd = new TH2F("hxyeffgd", "Efficiency GEM vs XY good DCH tracks", 70, -140, 140, 60, -90, 90);
524 hxyeffgd->SetOption("COLZ");
525
526 TH1F* hdtx = new TH1F("hdtx", "dtx", 100, -0.25, +0.25);
527 TH1F* hdty = new TH1F("hdty", "dty", 100, -0.25, +0.25);
528 TH1F* hDdx = new TH1F("hDdx", "Err dx", 100, 0, 10);
529 TH1F* hDdy = new TH1F("hDdy", "Err dy", 100, 0, 10);
530 TH1F* hDdtx = new TH1F("hDdtx", "Err dtx", 100, 0, 0.05);
531 TH1F* hDdty = new TH1F("hDdty", "Err dty", 100, 0, 0.05);
532 TH1F* hpdx = new TH1F("hpdx", "Pull dx", 500, -2.5, +2.5);
533 TH1F* hpdy = new TH1F("hpdy", "Pull dy", 500, -2.5, +2.5);
534 TH1F* hpdtx = new TH1F("hpdtx", "Pull dtx", 500, -25., +25.);
535 TH1F* hpdty = new TH1F("hpdty", "Pull dty", 500, -25., +25.);
536 TH1F* hchi2m = new TH1F("hChi2m", "Chi2 GemDch track", 1000, 0, 100);
537 TH1F* hchi2 = new TH1F("hchi2", "Chi2 tracks matching", 1000, 0, 100);
538
539 TH1F* hrej1 = new TH1F("hrej1", "Reject because not own hit", 5, 0, 5);
540 TH1F* hrej2 = new TH1F("hrej2", "Reject because DX", 5, 0, 5);
541 TH1F* hrej3 = new TH1F("hrej3", "Reject because DY", 5, 0, 5);
542 TH1F* hrej5 = new TH1F("hrej5", "Reject because DIST", 5, 0, 5);
543 TH1F* hrej6 = new TH1F("hrej6", "Reject because TIME", 5, 0, 5);
544
545 TProfile2D* tcalibrs = new TProfile2D("tcalibr", "Time calibration based on proton mass, strips", NUMBER_CHAMBERS,
546 0, NUMBER_CHAMBERS, 32, 0, 32, -1., +1., "e");
547 TProfile* tcalibrc = new TProfile("tcalibrc", "Time calibration based on proton mass, chambers", NUMBER_CHAMBERS, 0,
548 NUMBER_CHAMBERS, -1., +1., "e");
549
550 TProfile* tpeak = new TProfile("tpeak", "Peak time vs chamber", NUMBER_CHAMBERS, 0, NUMBER_CHAMBERS, 0., +30., "e");
551 TProfile* bpeak =
552 new TProfile("bpeak", "Peak velosity vs chamber", NUMBER_CHAMBERS, 0, NUMBER_CHAMBERS, 0., +1.5, "e");
553
554 char nam[64], tit[64];
555
556 float cvel = 29.97925f;
557 float ekin = 3.2f, ebeam = ekin + MP, pbeam = TMath::Sqrt(ebeam * ebeam + MP * MP);
558 float cbeam = pbeam / ebeam;
559 float cpeak = cbeam * CFACTOR * cvel;
560 // float ztarg = -2.2f; // EXP
561 float ztarg = -2.3f; // MC
562 float zc[NUMBER_CHAMBERS] = {0.f};
563
564 for (int i = 0; i < NUMBER_CHAMBERS; i++) {
565 TOF2->get_chamber_z(i, &zc[i]);
566 }
567
568 Int_t startEvent = 0;
569 Int_t nEvents = bmnTree->GetEntries();
570 Int_t nEventsMC = nEvents; // bmnTreeMC->GetEntries();
571 printf("\nNumber of events in tracks chain is %d (MC chain %d)\n", nEvents, nEventsMC);
572
573 int Ntr = 0, Ntr1 = 0, Ntr1g = 0, Ntr0 = 0, Ntr700 = 0;
574 Int_t iEv = 0, iEvMC = 0, fCurrent = -1;
575 UInt_t trackEvId = 0, MCEvId = 0;
576 UInt_t trackRunId = 0, MCRunId = 0;
577 // nEvents = 10;
578 TH1D *projtmp = 0, *projtmp1 = 0;
579 for (Int_t iEv0 = startEvent; iEv0 < startEvent + nEvents; iEv0++) {
580 if (iEv % 10000 == 0) {
581 cout << "Event: " << iEv << "/" << startEvent + nEvents << " ( Run " << trackRunId << " EvId " << trackEvId
582 << " )" << endl;
583 }
584 BmnGemDchTrack* gemdchtrack0 = NULL;
585 // do
586 // {
587 if (iEv >= nEvents)
588 break;
589 bmnTree->GetEntry(iEv);
590 // if (mtracks == NULL) { iEv++; continue; };
591 gemdchtrack0 = (BmnGemDchTrack*)(mtracks->UncheckedAt(0));
592 // gemdchtrack0 = (BmnGemDchTrack*)(mtracks->At(0));
593 // printf("First attempt, iEv = %d, iEv0 = %d\n",iEv,iEv0);
594 iEv++;
595 // }
596 // while (gemdchtrack0 == NULL);
597 // printf("Point 0, iEv = %d\n",iEv);
598 if (iEv >= nEvents)
599 break;
600 // trackEvId = gemdchtrack0->GetEventId();
601 /*
602 trackEvId = ((BmnEventHeader*)eventHeader)->GetEventId();
603 // printf("Point 1, iEv = %d, trackEvId = %d\n",iEv,trackEvId);
604 if (eventHeader == NULL) continue;
605 // trackRunId = ((BmnEventHeader*)eventHeader->UncheckedAt(0))->GetRunId();
606 trackRunId = ((BmnEventHeader*)eventHeader)->GetRunId();
607 // printf("Point 2, iEv = %d, trackRunId = %d\n",iEv,trackRunId);
608 if (MCRunId == trackRunId && MCEvId > trackEvId) continue;
609 if (MCRunId > trackRunId) continue;
610 // printf("Point 3, iEv = %d\n",iEv);
611 if (!(MCEvId==trackEvId && MCRunId==trackRunId))
612 {
613 do
614 { */
615 // bmnTreeMC->GetEntry(iEvMC);
616 // if (eventHeaderMC == NULL) break;
617 // MCEvId = ((BmnEventHeader*) eventHeaderMC->UncheckedAt(0))->GetEventId();
618 // MCRunId = ((BmnEventHeader*) eventHeaderMC->UncheckedAt(0))->GetRunId();
619 // cout << "Tracks Event: " << iEv << ", MC Event: " << iEvMC << endl;
620 // cout << "Tracks EventId: " << trackEvId << ", MC EventId: " << MCEvId << endl;
621 // cout << "Tracks RunId: " << trackRunId << ", MC RunId: " << MCRunId << endl;
622 // iEvMC++;
623 // if (iEvMC >= nEventsMC) break;
624 // if (MCRunId > trackRunId) break;
625 // if (MCRunId == trackRunId && MCEvId > trackEvId) break;
626 // }
627 // while (!(MCEvId==trackEvId && MCRunId==trackRunId));
628 // }
629 // if (eventHeaderMC == NULL) continue;
630 // if (!(MCEvId==trackEvId && MCRunId==trackRunId))
631 // {
632 // cout << "Not found EventId " << trackEvId << " in MC Data file" << endl;
633 // cout << "Tracks Event: " << iEv << ", MC Event: " << iEvMC << endl;
634 // cout << "Tracks RunId: " << trackRunId << ", MC RunId: " << MCRunId << endl;
635 // cout << "Tracks EventId: " << trackEvId << ", MC EventId: " << MCEvId << endl;
636 // if (iEvMC >= nEventsMC) break;
637 // if (MCRunId > trackRunId) continue;
638 // if (MCRunId == trackRunId && MCEvId > trackEvId) continue;
639 // break;
640 // continue;
641 // }
642 int istr[NUMBER_CHAMBERS] = {-1};
643
644 for (int i = 0; i < NUMBER_CHAMBERS; i++) {
645 nhit[i] = 0;
646 for (int j = 0; j < 32; j++)
647 nhits[i][j] = 0;
648 for (int j = 0; j < 32; j++)
649 tof[i][j] = 0.;
650 }
651 nhi = 0;
652
653 int nsts = stspoints->GetEntriesFast();
654 hnsts->Fill(nsts);
655 int ntof = tof700Digits->GetEntriesFast();
656 hntof->Fill(ntof);
657 int ntrmc = MCtracks->GetEntriesFast();
658 CbmMCTrack* mctrack = NULL;
659 float x = 0., y = 0., z = 0.;
660 int hchamb[1000] = {0}, hstrip[1000] = {0};
661 int hmctr[1000] = {0};
662 int clchamb[1000] = {0}, clstrip[1000] = {0};
663 Bool_t usedstrip[1000] = {false};
664 Bool_t usedstrip1[1000] = {false};
665 Int_t numhit[NUMBER_CHAMBERS][32] = {{0}};
666 BmnTOFPoint *tofpoint[1000] = {NULL}, *tofpoin = NULL;
667 Int_t toftrackid[1000] = {-1}, toftrid = -1;
668 Int_t indpoint = -1;
669 if (MCDEBUG)
670 printf("\n*************** TOF hits %d %d **********************\n", iEv0, ntof);
671 for (Int_t iDig = 0; iDig < ntof; ++iDig) {
672 BmnTofHit* digit = (BmnTofHit*)tof700Digits->At(iDig);
673 if (digit == NULL)
674 continue;
675 UInt_t DetId = digit->GetDetectorID();
676 Short_t plane = ((DetId & 0xFFFF) >> 8) - 1;
677 Short_t planeexp = order[plane];
678 Short_t strip = (DetId & 0xFF) - 1;
679 Float_t time = digit->GetTimeStamp();
680 // if (iEv % 10 == 0 && iEv < 100) cout << "plane = " << plane << " plane(exp) = " << planeexp << "
681 // chamberID = " << idchambers[planeexp] << " strip = " << strip << " time = " << time << endl;
682 if (planeexp >= NUMBER_CHAMBERS)
683 continue;
684 // TOF2->get_hit_xyzng(planeexp, strip, 0, &x, &y, &z);
685 // if (z == 0.) continue;
686 plane = planeexp;
687 htof->Fill(time);
688 hchamb[nhi] = plane;
689 hstrip[nhi] = strip;
690 nhit[plane]++;
691 nhits[plane][strip]++;
692 tof[plane][strip] = time;
693 istr[plane] = strip;
694 x = digit->GetX();
695 y = digit->GetY();
696 z = digit->GetZ();
697 xhi[nhi] = x;
698 yhi[nhi] = y;
699 zhi[nhi] = z;
700 hxytof->Fill(x, y);
701 TOF2->get_strip_xyz(plane, strip, &x, &y, &z);
702 xhic[nhi] = x;
703 yhic[nhi] = y;
704 zhic[nhi] = z;
705 hstr->Fill(strip);
706 hcha->Fill(plane);
707 numhit[plane][strip] = nhi;
708 indpoint = digit->GetRefIndex();
709 tofpoin = (BmnTOFPoint*)tofpoints->At(indpoint);
710 if (!tofpoin)
711 continue;
712 toftrid = tofpoin->GetTrackID();
713 mctrack = NULL;
714 if (MCDEBUG)
715 printf(" *********\n");
716 if (toftrid >= 0)
717 mctrack = (CbmMCTrack*)MCtracks->At(toftrid);
718 if (mctrack) {
719 if (mctrack->GetMotherId() < 0 && mctrack->GetStartZ() < 1.f) {
720 if (MCDEBUG) {
721 printf(" hit %d\n", iDig);
722 digit->Print("");
723 printf(" point %d\n", indpoint);
724 tofpoin->Print("");
725 printf(" track %d\n", toftrid);
726 LOG(info) << "Track " << toftrid << ", mother : " << mctrack->GetMotherId() << ", Type "
727 << mctrack->GetPdgCode() << ", momentum (" << mctrack->GetPx() << ", "
728 << mctrack->GetPy() << ", " << mctrack->GetPz() << ") GeV";
729 LOG(info) << ", vertex (" << mctrack->GetStartX() << ", " << mctrack->GetStartY() << ", "
730 << mctrack->GetStartZ() << ") cm";
731 LOG(info) << " Ref " << mctrack->GetNPoints(kREF) << ", MVD " << mctrack->GetNPoints(kMVD)
732 << ", GEM " << mctrack->GetNPoints(kGEM) << ", TOF1 " << mctrack->GetNPoints(kTOF1)
733 << ", DCH1 " << mctrack->GetNPoints(kDCH) << ", TOF2 " << mctrack->GetNPoints(kTOF);
734 }
735 tofpoint[nhi] = (BmnTOFPoint*)tofpoints->At(indpoint);
736 toftrackid[nhi] = tofpoint[nhi]->GetTrackID();
737 } else {
738 if (MCDEBUG) {
739 printf(" reject hit %d\n", iDig);
740 digit->Print("");
741 printf(" reject point %d\n", indpoint);
742 tofpoin->Print("");
743 printf(" reject track %d\n", toftrid);
744 LOG(info) << "Track " << toftrid << ", mother : " << mctrack->GetMotherId() << ", Type "
745 << mctrack->GetPdgCode() << ", momentum (" << mctrack->GetPx() << ", "
746 << mctrack->GetPy() << ", " << mctrack->GetPz() << ") GeV";
747 LOG(info) << ", vertex (" << mctrack->GetStartX() << ", " << mctrack->GetStartY() << ", "
748 << mctrack->GetStartZ() << ") cm";
749 LOG(info) << " Ref " << mctrack->GetNPoints(kREF) << ", MVD " << mctrack->GetNPoints(kMVD)
750 << ", GEM " << mctrack->GetNPoints(kGEM) << ", TOF1 " << mctrack->GetNPoints(kTOF1)
751 << ", DCH1 " << mctrack->GetNPoints(kDCH) << ", TOF2 " << mctrack->GetNPoints(kTOF);
752 }
753 tofpoint[nhi] = NULL;
754 toftrackid[nhi] = -10;
755 }
756 } else {
757 if (MCDEBUG) {
758 printf(" hit %d\n", iDig);
759 digit->Print("");
760 printf(" point %d\n", indpoint);
761 tofpoin->Print("");
762 printf(" track %d null\n", toftrid);
763 }
764 }
765 nhi++;
766 } // tof700Digits loop
767 hnhit->Fill(nhi);
768
769 float xyz[3] = {0.}, cxyz[3] = {0.}, xrpc = 0., yrpc = 0., dz = 100., tx = 0., ty = 0., zcha = 0.;
770 float beta = 1.0f, mass = 1.0f, mass2 = 1.0f, tcor = 0., tcorc = 0., tcors = 0., tcorspi = 0., tcorcpi = 0.;
771 int nh = 0;
772 int ntr = mtracks->GetEntriesFast();
773 int nmatch = 0;
774 if (MCMATCH)
775 if (mcmatch)
776 nmatch = mcmatch->GetEntriesFast();
777 hntr->Fill(ntr);
778 Ntr += ntr;
779 int ntr0 = 0, ntr700 = 0;
780 int ntr1 = 0, ntr1g = 0;
781 int iflag = 0, igem = 0, idch = 0;
782 float chi2 = 0.f;
783 BmnGemDchTrack* gemdchtrack = NULL;
784 BmnTrack* track = NULL;
785 CbmStsTrack* gtrack = NULL; // Gleb GEM tracks
786 CbmTrackMatch* matchinfo = NULL; // Gleb GEM track match to MC
787 Bool_t gooddch = false, goodgem = false, goodgemdch = false, gooddchtof = false, goodgemdchtof = false;
788 Double_t xGem = 0.;
789 Double_t yGem = 0.;
790 Double_t zGem = 0.;
791 Double_t txGem = 0.;
792 Double_t tyGem = 0.;
793 Double_t QpGem = 0.;
794 Double_t QpDch = 0.;
795 Int_t nGemHits = 0;
796 Int_t nSiHits = 0;
797 Double_t Xpv = 0.;
798 Double_t Ypv = 0.;
799 Double_t dxm = 0.;
800 Double_t dym = 0.;
801 Double_t dtx = 0.;
802 Double_t dty = 0.;
803 Double_t Ddx = 0.;
804 Double_t Ddy = 0.;
805 Double_t Ddtx = 0.;
806 Double_t Ddty = 0.;
807 Double_t Chi2m = 0.;
808 Double_t xgoodtr = 0.;
809 Double_t ygoodtr = 0.;
810 Double_t le = 0., p = 0., ps = 0., pmod = 0., lenr = 0.;
811 int himin = -1, himin1 = -1;
812 double dx, dy, dl;
813 double dist = 0., distmin = 1000000000.0f, distmin1 = 1000000000.0f, tofmin = 1000000000.0f;
814 Float_t xc, yc, zcc;
815 int ic;
816 int is;
817 double cosa = 0.;
818 float hidist = 0.;
819 Double_t tofr0 = 0.f;
820 Double_t tofr = 0.f;
821 Double_t tofp = 0.f;
822 Double_t tofpi = 0.f;
823 float massmin1 = 0.;
824 float massmax1 = 0.;
825 float massmin2 = 0.;
826 float massmax2 = 0.;
827
828 float massmaxpi = 0.;
829 float massminpi = 0.;
830 float tcorpi = 0.;
831 float tofrc = 0.;
832 float tofrs = 0.;
833 float betac = 0.;
834 float betas = 0.;
835 float masss2 = 0.;
836 float masss = 0.;
837 float massc2 = 0.;
838 float massc = 0.;
839 double timepeak = 0., betapeak = 0.;
840 Int_t mctrid = -1, ntrue = 0, nwrong = 0, nfake = 0, nhall = 0, nmctr = 0;
841 Float_t ptrue = 0.f;
842
843 Double_t MCm = 0., MCp = 0., MCe = 0., MCmz = 0., MCpz = 0., MCl = 0., MCt = 0., sMCpz = 0., sMCl = 0.,
844 sMCt = 0., sMCb = 0., rMCmz2 = 0.;
845 Int_t MCim = 0, MCrim = 0, ipart = 0;
846 Int_t pdg = 0, MCn = 0, MCz = 0, MCnt = 0;
847 // sort by M/Z
848 Double_t pfra[15000], pzfra[15000] = {0.}, mfra[15000] = {0.}, gfra[15000] = {0.};
849 Int_t index[15000] = {0};
850 Double_t MCpmold = 100.;
851 for (Int_t iTrmc = 0; iTrmc < ntrmc; ++iTrmc) {
852 mctrack = (CbmMCTrack*)MCtracks->At(iTrmc);
853 if (mctrack == NULL)
854 continue;
855 pfra[iTrmc] = mctrack->GetP();
856 pdg = mctrack->GetPdgCode();
857 if (pdg < 1000000000) {
858 MCim = 0;
859 if (TDatabasePDG::Instance()) {
860 TParticlePDG* particle = TDatabasePDG::Instance()->GetParticle(pdg);
861 if (particle) {
862 MCm = particle->Mass(); /* printf("pdg %d ch %d ma %f\n",pdg,MCz,MCm);*/
863 } else
864 continue;
865 MCim = (int)(MCm / 0.93);
866 } else {
867 continue;
868 }
869 } else {
870 pdg -= 1000000000;
871 pdg = pdg % 10000000;
872 pdg = pdg % 10000;
873 MCim = (int)(pdg / 10);
874 }
875 mfra[iTrmc] = MCim;
876 gfra[iTrmc] = 0.;
877 if (MCim >= 1)
878 gfra[iTrmc] = pfra[iTrmc] / MCim;
879 }
880 TMath::Sort(ntrmc, gfra, index, true);
881 //
882 if (MCDEBUG)
883 printf("\n*************** MC tracks %d %d **********************\n", iEv0, ntrmc);
884 Int_t iTrmc = 0;
885 Int_t SumA = 0, SumA0 = 0;
886 for (Int_t iTrmc0 = 0; iTrmc0 < ntrmc; ++iTrmc0) {
887 iTrmc = iTrmc0;
888 if (SORT)
889 iTrmc = index[iTrmc0];
890 mctrack = (CbmMCTrack*)MCtracks->At(iTrmc);
891 if (mctrack == NULL)
892 continue;
893 if (SORT) {
894 SumA0 += (int)mfra[iTrmc];
895 }
896 if (SORT)
897 printf("iTrmc0 %d iTrmc %d gamma %f A %d SumA %d\n", iTrmc0, iTrmc, gfra[iTrmc], (int)mfra[iTrmc],
898 SumA0);
899 if (MCDEBUG) {
900 LOG(info) << "Track " << iTrmc << ", mother : " << mctrack->GetMotherId() << ", Type "
901 << mctrack->GetPdgCode() << ", momentum (" << mctrack->GetPx() << ", " << mctrack->GetPy()
902 << ", " << mctrack->GetPz() << ") GeV";
903 LOG(info) << ", vertex (" << mctrack->GetStartX() << ", " << mctrack->GetStartY() << ", "
904 << mctrack->GetStartZ() << ") cm";
905 LOG(info) << " Ref " << mctrack->GetNPoints(kREF) << ", MVD " << mctrack->GetNPoints(kMVD)
906 << ", GEM " << mctrack->GetNPoints(kGEM) << ", TOF1 " << mctrack->GetNPoints(kTOF1)
907 << ", DCH1 " << mctrack->GetNPoints(kDCH) << ", TOF2 " << mctrack->GetNPoints(kTOF);
908 }
909 if (mctrack->GetMotherId() >= 0)
910 continue;
911 if (mctrack->GetStartZ() > 1.f)
912 continue;
913 pdg = mctrack->GetPdgCode();
914 if (pdg < 1000000000) {
915 ipart = 1;
916 MCim = 0;
917 // if (true) continue;
918 if (TDatabasePDG::Instance()) {
919 TParticlePDG* particle = TDatabasePDG::Instance()->GetParticle(pdg);
920 if (particle) {
921 MCz = (int)(particle->Charge() / 3.);
922 MCm = particle->Mass(); /* printf("pdg %d ch %d ma %f\n",pdg,MCz,MCm);*/
923 } else
924 continue;
925 if (MCz <= 0)
926 continue;
927 MCim = (int)(MCm / 0.93);
928 MCmz = ((double)MCim) / MCz;
929 } else {
930 continue;
931 }
932 if (MCm > 0.93) {
933 hMCamass->Fill(MCmz);
934 hMCmass0->Fill(MCmz);
935 hMCavsz->Fill(MCz, MCim);
936 }
937 } else {
938 ipart = 0;
939 pdg -= 1000000000;
940 pdg = pdg % 10000000;
941 MCz = pdg / 10000;
942 if (MCz <= 0)
943 continue;
944 pdg = pdg % 10000;
945 MCim = (int)(pdg / 10);
946 MCm = MCim * MP;
947 MCmz = ((double)MCim) / MCz;
948 if (MCm > 0.93) {
949 hMCamass->Fill(MCmz);
950 hMCrmass0->Fill(MCmz);
951 hMCravsz->Fill(MCz, MCim);
952 }
953 }
954 Bool_t goodfra = true;
955 if (SORT && MCim >= 1) {
956 SumA += MCim;
957 if (SumA > AARGON)
958 goodfra = false;
959 }
960 // if ((MCm = mctrack->GetMass()) < 0.93) continue;
961 MCp = mctrack->GetP();
962 if (((MCp / MCim) > MCpmold) && SORT && (MCim >= 1.)) {
963 // printf("Previous gamma = %f, current one = %f, track #%d(%d), SumA = %d, IEv = %d !\n", MCpmold,
964 // MCp/MCim, iTrmc0, iTrmc, SumA, iEv0);
965 }
966 MCpmold = MCp / MCim;
967 if (MCp < 0.01)
968 continue;
969 MCnt = mctrack->GetNPoints(kTOF);
970 MCn = mctrack->GetNPoints(kGEM);
971 // if (MCm > 0.93 && MCm < 1.07)
972 // {
973 // hnhgem->Fill(MCn);
974 // if (MCnt) hnhgem1->Fill(MCn);
975 // else hnhgem0->Fill(MCn);
976 // }
977 if (goodfra) {
978 if (MCm > 0.93) {
979 pzfra[iTrmc0] = MCp / MCz;
980 hMCpndf->Fill(MCp / MCz);
981 hpangle->Fill(MCp / MCz, TMath::RadToDeg() * TMath::ACos(mctrack->GetPz() / MCp));
982 hpanglep->Fill(MCp / MCz, TMath::RadToDeg() * TMath::ACos(mctrack->GetPz() / MCp));
983 hangle->Fill(TMath::RadToDeg() * TMath::ACos(mctrack->GetPz() / MCp));
984 }
985 if (MCm > 0.93 && MCn < MINGEMNDF1) {
986 hMCpndf0->Fill(MCp / MCz);
987 }
988 if (MCm > 0.93 && MCn >= MINGEMNDF1) {
989 h3MCamass->Fill(MCmz);
990 hMCpndf1->Fill(MCp / MCz);
991 if (MCnt > 0)
992 hMCpndf2->Fill(MCp / MCz);
993 }
994 if (MCnt < 1)
995 continue;
996 if (MCm > 0.93 && ipart) {
997 h4MCamass->Fill(MCmz);
998 hpangle2->Fill(MCp / MCz, TMath::RadToDeg() * TMath::ACos(mctrack->GetPz() / MCp));
999 hpanglep2->Fill(MCp / MCz, TMath::RadToDeg() * TMath::ACos(mctrack->GetPz() / MCp));
1000 hangle2->Fill(TMath::RadToDeg() * TMath::ACos(mctrack->GetPz() / MCp));
1001 }
1002 if (MCm > 0.93 && !ipart) {
1003 h4MCamass->Fill(MCmz);
1004 hpangle2->Fill(MCp / MCz, TMath::RadToDeg() * TMath::ACos(mctrack->GetPz() / MCp));
1005 hpanglep2->Fill(MCp / MCz, TMath::RadToDeg() * TMath::ACos(mctrack->GetPz() / MCp));
1006 hangle2->Fill(TMath::RadToDeg() * TMath::ACos(mctrack->GetPz() / MCp));
1007 }
1008 if (MCn < MINGEMNDF1)
1009 continue;
1010 if (MCm > 0.93 && ipart) {
1011 h1MCamass->Fill(MCmz);
1012 h1MCmass0->Fill(MCmz);
1013 h1MCavsz->Fill(MCz, MCim);
1014 hpangle1->Fill(MCp / MCz, TMath::RadToDeg() * TMath::ACos(mctrack->GetPz() / MCp));
1015 hpanglep1->Fill(MCp / MCz, TMath::RadToDeg() * TMath::ACos(mctrack->GetPz() / MCp));
1016 hangle1->Fill(TMath::RadToDeg() * TMath::ACos(mctrack->GetPz() / MCp));
1017 }
1018 if (MCm > 0.93 && !ipart) {
1019 h1MCamass->Fill(MCmz);
1020 h1MCrmass0->Fill(MCmz);
1021 h1MCravsz->Fill(MCz, MCim);
1022 hpangle1->Fill(MCp / MCz, TMath::RadToDeg() * TMath::ACos(mctrack->GetPz() / MCp));
1023 hpanglep1->Fill(MCp / MCz, TMath::RadToDeg() * TMath::ACos(mctrack->GetPz() / MCp));
1024 hangle1->Fill(TMath::RadToDeg() * TMath::ACos(mctrack->GetPz() / MCp));
1025 }
1026 if (MCm > 0.93 && MCp < 2.0) {
1027 h2MCamass->Fill(MCmz);
1028 }
1029 }
1030 MCe = mctrack->GetEnergy();
1031 MCpz = MCp / MCz;
1032 sMCpz = MCpz; // + gRandom->Gaus()*(MCpz*0.04);
1033 Int_t hused[10000] = {0};
1034 for (int ih = 0; ih < nhi; ih++) {
1035 if (iTrmc == toftrackid[ih]) {
1036 MCl = tofpoint[ih]->GetLength();
1037 if (MCl < 540.)
1038 continue;
1039 MCt = tofpoint[ih]->GetTime();
1040 sMCl = MCl; // + gRandom->Gaus()*5.;
1041 sMCt = MCt; // + gRandom->Gaus()*0.115;
1042 sMCb = sMCl / sMCt / cvel;
1043 rMCmz2 = sMCpz * sMCpz * (1.f / sMCb / sMCb - 1.f);
1044 hMCmass2->Fill(rMCmz2);
1045 hMCmass->Fill(TMath::Sqrt(rMCmz2));
1046 if (hused[iTrmc] == 0) {
1047 hMCmassr->Fill(MCmz);
1048 if (MCp < 2.0)
1049 hMCmassrp->Fill(MCmz);
1050 hused[iTrmc] = 1;
1051 }
1052 hpbetamc->Fill(sMCpz, sMCb);
1053 if (pzfra[iTrmc0] > 7.)
1054 hp0xygem->Fill(tofpoint[ih]->GetX(), tofpoint[ih]->GetY());
1055 }
1056 }
1057 CbmStsPoint* stspoint = 0;
1058 Double_t x_sts = 0., y_sts = 0., z_sts = 0.;
1059 for (int ih = 0; ih < nsts; ih++) {
1060 stspoint = (CbmStsPoint*)stspoints->At(ih);
1061 if (stspoint == NULL)
1062 continue;
1063 if (iTrmc == stspoint->GetTrackID()) {
1064 x_sts = stspoint->GetXIn();
1065 y_sts = stspoint->GetYIn();
1066 z_sts = stspoint->GetZIn();
1067 if (pzfra[iTrmc0] > 7.)
1068 hpxysts->Fill(x_sts, y_sts);
1069 }
1070 }
1071 }
1072
1073 for (Int_t iTr = 0; iTr < ntr; ++iTr) {
1074 gemdchtrack = (BmnGemDchTrack*)mtracks->At(iTr);
1075 if (gemdchtrack == NULL)
1076 continue;
1077 track = (BmnTrack*)gemdchtrack->GetDchTrack();
1078 gooddch = false;
1079 if (track != NULL) {
1080 QpDch = track->GetParamLast()->GetQp();
1081 if (!TMath::IsNaN(QpDch) && (track->GetParamLast()->GetZ() > 500.f)) {
1082 // get xyz and cxyz here
1083 xyz[0] = track->GetParamLast()->GetX();
1084 xyz[1] = track->GetParamLast()->GetY();
1085 xyz[2] = track->GetParamLast()->GetZ();
1086 tx = track->GetParamLast()->GetTx();
1087 ty = track->GetParamLast()->GetTy();
1088 cxyz[0] = tx / TMath::Sqrt(1. + tx * tx + ty * ty);
1089 cxyz[1] = ty / TMath::Sqrt(1. + tx * tx + ty * ty);
1090 cxyz[2] = 1. / TMath::Sqrt(1. + tx * tx + ty * ty);
1091 xgoodtr = xyz[0] + tx * 100.;
1092 ygoodtr = xyz[1] + ty * 100.;
1093 hxydch->Fill(xgoodtr, ygoodtr);
1094 gooddch = true;
1095 }
1096 }
1097 Double_t mcp = 0.;
1098 if (MCMATCH) {
1099 if (nmatch == 0)
1100 continue;
1101 if (iTr >= nmatch)
1102 continue;
1103 matchinfo = (CbmTrackMatch*)mcmatch->At(iTr);
1104 if (matchinfo == NULL)
1105 continue;
1106 mctrid = matchinfo->GetMCTrackId();
1107 /** Number of true hits on track (from matched MC track) **/
1108 ntrue = matchinfo->GetNofTrueHits();
1109 /** Number of wrong hits on track (from other MC tracks) **/
1110 nwrong = matchinfo->GetNofWrongHits();
1111 /** Number of fake hits on track (from no MC track) **/
1112 nfake = matchinfo->GetNofFakeHits();
1113 /** Number of MCTracks with common hits **/
1114 nmctr = matchinfo->GetNofMCTracks();
1115 nhall = matchinfo->GetNofHits();
1116 ptrue = matchinfo->GetTrueOverAllHitsRatio();
1117 htruehits->Fill(ptrue);
1118 hallhits->Fill(nhall);
1119 // if (iEv < 10) printf("mcid %d true %d wrong %d fake %d Nmctr %d quality %f\n", mctrid, ntrue, nwrong,
1120 // nfake, nmctr, ptrue);
1121 mctrack = NULL;
1122 if (mctrid >= 0 && mctrid < ntrmc)
1123 mctrack = (CbmMCTrack*)MCtracks->At(mctrid);
1124 if (mctrack)
1125 mcp = mctrack->GetP();
1126 }
1127 gtrack = (CbmStsTrack*)gemdchtrack->GetGemTrack(); // Merts GEM tracks
1128 goodgem = false;
1129 if (gtrack != NULL && (ptrue >= MINTRUE || !MCMATCH)) {
1130 // Selection criteria
1131 QpGem = gtrack->GetParamLast()->GetQp();
1132 if (!TMath::IsNaN(QpGem)) {
1133 if (QpGem != 0.) {
1134 hp0->Fill(1. / QpGem);
1135 } else {
1136 // hp0->Fill(0.);
1137 }
1138 // ntr1++;
1139 xGem = gtrack->GetParamLast()->GetX();
1140 yGem = gtrack->GetParamLast()->GetY();
1141 zGem = gtrack->GetParamLast()->GetZ();
1142 txGem = gtrack->GetParamLast()->GetTx();
1143 tyGem = gtrack->GetParamLast()->GetTy();
1144 // nGemHits = gtrack->GetNDF(); // Merts GEM tracks
1145 nGemHits = gemdchtrack->GetGemHits(); // Gleb GEM tracks
1146 nSiHits = gemdchtrack->GetSiHits(); // Gleb GEM tracks
1147 Xpv = gemdchtrack->GetXpv();
1148 Ypv = gemdchtrack->GetYpv();
1149
1150 if (GEMCUTS) {
1151 // Xpv < -3.5f || Xpv > 4.0f || Ypv < -1.0f || Ypv > 6.0f ||
1152 Bool_t goodpv = true;
1153 if (PRIMARY_VERTEX_CUTS) {
1154 if (ntr < 2)
1155 goodpv = false;
1156 if (primaryVertex->GetX() < -1 || primaryVertex->GetX() > 1 || primaryVertex->GetY() < -1
1157 || primaryVertex->GetY() > 1 || primaryVertex->GetZ() < -5 || primaryVertex->GetZ() > 3)
1158 goodpv = false;
1159 };
1160 if ((!goodpv) || nGemHits < MINGEMNDF || nSiHits < MINSIHITS || xGem < -140.f || xGem > 240.f
1161 || yGem < -40.f || yGem > 200.f || txGem < -0.5f || txGem > 1.0f || tyGem < -0.1f
1162 || tyGem > 0.5f || Xpv < -1.0f || Xpv > 1.0f || Ypv < -1.0f || Ypv > 1.0f || QpGem == 0.0f
1163 || 1.0f / TMath::Abs(QpGem) > 15.0f || 1.0f / TMath::Abs(QpGem) < 0.15f || zGem < 500.f)
1164 {
1165 if (GEMDEBUG)
1166 printf("Nhits %d xGem %f yGem %f zGem %f txGem %f tyGem %f Xpv %f Ypv %f QpGem %f\n",
1167 nGemHits, xGem, yGem, zGem, txGem, tyGem, Xpv, Ypv, QpGem);
1168 } else {
1169 if (GEMDEBUG)
1170 printf("Good GEM track!\n");
1171 if (ACCEPTANCE) {
1172 Double_t xgema = xGem + txGem * 100.;
1173 Double_t ygema = yGem + tyGem * 100.;
1174 if (xgema >= xmintof700 && xgema <= xmaxtof700 && ygema >= ymintof700
1175 && ygema <= ymaxtof700)
1176 {
1177 goodgem = true;
1178 le = gemdchtrack->GetLength() - ztarg;
1179 p = QpGem;
1180 ps = p;
1181 if (p != 0.)
1182 p = 1. / p;
1183 hptr->Fill(p); // GeV/c
1184 hltr->Fill(le); // cm
1185 pmod = TMath::Abs(p);
1186 hxygem->Fill(xGem, yGem);
1187 if (!MCMATCH)
1188 hpxygem->Fill(xGem, yGem);
1189 else if (pzfra[mctrid] > 7.)
1190 hpxygem->Fill(xGem, yGem);
1191 }
1192 } else {
1193 goodgem = true;
1194 le = gemdchtrack->GetLength() - ztarg;
1195 p = QpGem;
1196 ps = p;
1197 if (p != 0.)
1198 p = 1. / p;
1199 hptr->Fill(p); // GeV/c
1200 hltr->Fill(le); // cm
1201 pmod = TMath::Abs(p);
1202 hxygem->Fill(xGem, yGem);
1203 if (!MCMATCH)
1204 hpxygem->Fill(xGem, yGem);
1205 else if (pzfra[mctrid] > 7.)
1206 hpxygem->Fill(xGem, yGem);
1207 }
1208 }
1209 } else {
1210 if (nGemHits >= MINGEMNDF && nSiHits >= MINSIHITS && QpGem != 0.0f && zGem >= 500.f) {
1211 if (ACCEPTANCE) {
1212 Double_t xgema = xGem + txGem * 100.;
1213 Double_t ygema = yGem + tyGem * 100.;
1214 if (xgema >= xmintof700 && xgema <= xmaxtof700 && ygema >= ymintof700
1215 && ygema <= ymaxtof700)
1216 {
1217 goodgem = true;
1218 le = gemdchtrack->GetLength() - ztarg;
1219 p = QpGem;
1220 ps = p;
1221 if (p != 0.)
1222 p = 1. / p;
1223 hptr->Fill(p); // GeV/c
1224 hltr->Fill(le); // cm
1225 pmod = TMath::Abs(p);
1226 hxygem->Fill(xGem, yGem);
1227 if (!MCMATCH)
1228 hpxygem->Fill(xGem, yGem);
1229 else if (pzfra[mctrid] > 7.)
1230 hpxygem->Fill(xGem, yGem);
1231 }
1232 } else {
1233 goodgem = true;
1234 le = gemdchtrack->GetLength() - ztarg;
1235 p = QpGem;
1236 ps = p;
1237 if (p != 0.)
1238 p = 1. / p;
1239 hptr->Fill(p); // GeV/c
1240 hltr->Fill(le); // cm
1241 pmod = TMath::Abs(p);
1242 hxygem->Fill(xGem, yGem);
1243 if (!MCMATCH)
1244 hpxygem->Fill(xGem, yGem);
1245 else if (pzfra[mctrid] > 7.)
1246 hpxygem->Fill(xGem, yGem);
1247 }
1248 }
1249 }
1250 // if (iEv < 10 && goodgem && mctrack != NULL) printf("MC momentum %f REC momentum %f\n",
1251 // mcp,p);
1252 if (track == NULL) {
1253 // Xpv < -3.5f || Xpv > 4.0f || Ypv < -1.0f || Ypv > 6.0f ||
1254 Bool_t goodpv = true;
1255 if (PRIMARY_VERTEX_CUTS) {
1256 if (ntr < 2)
1257 goodpv = false;
1258 if (primaryVertex->GetX() < -1 || primaryVertex->GetX() > 1 || primaryVertex->GetY() < -1
1259 || primaryVertex->GetY() > 1 || primaryVertex->GetZ() < -5 || primaryVertex->GetZ() > 3)
1260 goodpv = false;
1261 };
1262 if ((!goodpv) || nGemHits < MINGEMNDF || nSiHits < MINSIHITS || xGem < -140.f || xGem > 240.f
1263 || yGem < -40.f || yGem > 200.f || txGem < -0.5f || txGem > 1.0f || tyGem < -0.1f
1264 || tyGem > 0.5f || Xpv < -1.0f || Xpv > 1.0f || Ypv < -1.0f || Ypv > 1.0f || QpGem == 0.0f
1265 || 1.0f / TMath::Abs(QpGem) > 15.0f || 1.0f / TMath::Abs(QpGem) < 0.15f || zGem < 500.f)
1266 {
1267 QpDch = QpGem;
1268 if (!TMath::IsNaN(QpDch) && (zGem > 500.f)) {
1269 // get xyz and cxyz here
1270 xyz[0] = xGem;
1271 xyz[1] = yGem;
1272 xyz[2] = zGem;
1273 tx = txGem;
1274 ty = tyGem;
1275 cxyz[0] = tx / TMath::Sqrt(1. + tx * tx + ty * ty);
1276 cxyz[1] = ty / TMath::Sqrt(1. + tx * tx + ty * ty);
1277 cxyz[2] = 1. / TMath::Sqrt(1. + tx * tx + ty * ty);
1278 xgoodtr = xyz[0] + tx * 100.;
1279 ygoodtr = xyz[1] + ty * 100.;
1280 if (ACCEPTANCE) {
1281 if (xgoodtr >= xmintof700 && xgoodtr <= xmaxtof700 && ygoodtr >= ymintof700
1282 && ygoodtr <= ymaxtof700)
1283 {
1284 hxydch->Fill(xgoodtr, ygoodtr);
1285 gooddch = true;
1286 }
1287 }
1288 }
1289 }
1290 }
1291 }
1292 }
1293 goodgemdch = false;
1294 int nrej1 = 0, nrej2 = 0, nrej3 = 0, nrej4 = 0, nrej5 = 0, nrej6 = 0;
1295 if (goodgem && gooddch) {
1296 ntr1++;
1297 // fill matching parameters
1298 dxm = gemdchtrack->GetDx();
1299 dym = gemdchtrack->GetDy();
1300 dtx = gemdchtrack->GetDtx();
1301 dty = gemdchtrack->GetDty();
1302 Ddx = gemdchtrack->GetDdx();
1303 Ddy = gemdchtrack->GetDdy();
1304 Ddtx = gemdchtrack->GetDdtx();
1305 Ddty = gemdchtrack->GetDdty();
1306 Chi2m = gemdchtrack->GetChi2();
1307 hdx->Fill(dxm);
1308 hdy->Fill(dym);
1309 hdxy0->Fill(dxm, dym);
1310 hdtx->Fill(dtx);
1311 hdty->Fill(dty);
1312 hDdx->Fill(sqrt(Ddx));
1313 hDdy->Fill(sqrt(Ddy));
1314 hDdtx->Fill(sqrt(Ddtx));
1315 hDdty->Fill(sqrt(Ddty));
1316 hpdx->Fill(dxm / sqrt(Ddx));
1317 hpdy->Fill(dym / sqrt(Ddy));
1318 hpdtx->Fill(dtx / sqrt(Ddtx));
1319 hpdty->Fill(dty / sqrt(Ddty));
1320 hchi2m->Fill(Chi2m);
1321 chi2 = dxm * dxm / Ddx + dym * dym / Ddy;
1322 hchi2->Fill(chi2);
1323 // Selection criteria - x and y tracks matching
1324 if (dxm >= -DXGEMDCHMAX && dxm <= +DXGEMDCHMAX && dym >= -DYGEMDCHMAX
1325 && dym <= +DYGEMDCHMAX /* && Chi2m <= 30.f && chi2 > 30.f*/)
1326 {
1327 ntr1g++;
1328 hdxy->Fill(dxm, dym);
1329 hxygoodtr->Fill(xgoodtr, ygoodtr);
1330 hxygoodtrgem->Fill(xGem, yGem);
1331 hpgemdch->Fill(p); // GeV/c
1332 goodgemdch = true;
1333 }
1334 }
1335
1336 // find corresponding strip (if any)
1337 /* if (USETRUEHIT)
1338 {
1339 nh = nhi;
1340 himin = -1;
1341 distmin = 1000000000.0f;
1342 tofmin = 1000000000.0f;
1343 for (int ih = 0; (ih<nh) && goodgemdch; ih++)
1344 {
1345 if (usedstrip[ih]) continue;
1346 if (mctrid == toftrackid[ih])
1347 {
1348 if (tof[hchamb[ih]][hstrip[ih]] < tofmin)
1349 {
1350 himin = ih;
1351 distmin = 0.;
1352 tofmin = tof[hchamb[ih]][hstrip[ih]];
1353 }
1354 }
1355 }
1356
1357 nh = nhi;
1358 himin1 = -1;
1359 distmin1 = 1000000000.0f;
1360 tofmin = 1000000000.0f;
1361 for (int ih = 0; (ih<nh) && gooddch; ih++)
1362 {
1363 if (usedstrip1[ih]) continue;
1364 if (mctrid == toftrackid[ih])
1365 {
1366 if (tof[hchamb[ih]][hstrip[ih]] < tofmin)
1367 {
1368 himin1 = ih;
1369 distmin1 = 0.;
1370 tofmin = tof[hchamb[ih]][hstrip[ih]];
1371 }
1372 }
1373 }
1374 }
1375 else
1376 { */
1377 nh = nhi;
1378 himin = -1;
1379 dist = 0.;
1380 distmin = 1000000000.0f;
1381 for (int ih = 0; (ih < nh) && goodgemdch; ih++) {
1382 if (usedstrip[ih])
1383 continue;
1384 if (mctrid != toftrackid[ih] && USETRUEHIT) {
1385 nrej1++;
1386 continue;
1387 }
1388 ic = hchamb[ih];
1389 is = hstrip[ih];
1390 dz = zhi[ih] - xyz[2];
1391 xrpc = xyz[0] + dz * tx;
1392 yrpc = xyz[1] + dz * ty;
1393 if (!XYPANIN) {
1394 if ((xrpc - xhic[ih]) < -(halfwidthx[ic] + DXMAX) || (xrpc - xhic[ih]) > +(halfwidthx[ic] + DXMAX))
1395 continue;
1396 if ((yrpc - yhic[ih]) < -(halfwidthy[ic] + DYMAX) || (yrpc - yhic[ih]) > +(halfwidthy[ic] + DYMAX))
1397 continue;
1398 }
1399 dx = xrpc - xhi[ih];
1400 dy = yrpc - yhi[ih];
1401 if (XYPANIN) {
1402 if (TMath::Abs(dx) > DXMAXPAN) {
1403 nrej2++;
1404 continue;
1405 }
1406 if (TMath::Abs(dy) > DYMAXPAN) {
1407 nrej3++;
1408 continue;
1409 }
1410 dist = TMath::Sqrt(dx * dx + dy * dy);
1411 } else {
1412 dist = TMath::Sqrt(dx * dx / XFACTOR + dy * dy);
1413 }
1414 if (dist < distmin) {
1415 distmin = dist;
1416 himin = ih;
1417 }
1418 }
1419
1420 nh = nhi;
1421 himin1 = -1;
1422 dist = 0.;
1423 distmin1 = 1000000000.0f;
1424 for (int ih = 0; (ih < nh) && gooddch; ih++) {
1425 if (usedstrip1[ih])
1426 continue;
1427 if (mctrid != toftrackid[ih] && USETRUEHIT)
1428 continue;
1429 ic = hchamb[ih];
1430 is = hstrip[ih];
1431 dz = zhi[ih] - xyz[2];
1432 xrpc = xyz[0] + dz * tx;
1433 yrpc = xyz[1] + dz * ty;
1434 if (!XYPANIN) {
1435 if ((xrpc - xhic[ih]) < -(halfwidthx[ic] + DXMAX) || (xrpc - xhic[ih]) > +(halfwidthx[ic] + DXMAX))
1436 continue;
1437 if ((yrpc - yhic[ih]) < -(halfwidthy[ic] + DYMAX) || (yrpc - yhic[ih]) > +(halfwidthy[ic] + DYMAX))
1438 continue;
1439 }
1440 dx = xrpc - xhi[ih];
1441 dy = yrpc - yhi[ih];
1442 if (XYPANIN) {
1443 if (TMath::Abs(dx) > DXMAXPAN)
1444 continue;
1445 if (TMath::Abs(dy) > DYMAXPAN)
1446 continue;
1447 dist = TMath::Sqrt(dx * dx + dy * dy);
1448 } else {
1449 dist = TMath::Sqrt(dx * dx / XFACTOR + dy * dy);
1450 }
1451 if (dist < distmin1) {
1452 distmin1 = dist;
1453 himin1 = ih;
1454 }
1455 }
1456 // }
1457 // DCH+TOF700 matching
1458 // Selection criteria on minimal distance track-hit
1459 // Selection criteria time-of-flight
1460 gooddchtof = false;
1461 if (himin1 >= 0) {
1462 ic = hchamb[himin1];
1463 is = hstrip[himin1];
1464 if (distmin1 <= DISTMAX && tof[ic][is] <= TOFMAX) {
1465 gooddchtof = true;
1466 usedstrip1[himin1] = true;
1467 hxygoodhit1->Fill(xgoodtr, ygoodtr);
1468 hxygoodhit1a->Fill(xhi[himin1], yhi[himin1]);
1469 }
1470 }
1471 // GEM+DCH+TOF700 matching
1472 // Selection criteria on minimal distance track-hit
1473 // Selection criteria time-of-flight
1474 goodgemdchtof = false;
1475 if (himin >= 0) {
1476 ic = hchamb[himin];
1477 is = hstrip[himin];
1478 if (distmin > DISTMAX)
1479 nrej5++;
1480 if (tof[ic][is] > TOFMAX)
1481 nrej6++;
1482 if (distmin <= DISTMAX && tof[ic][is] <= TOFMAX) {
1483 ntr0++;
1484 goodgemdchtof = true;
1485 usedstrip[himin] = true;
1486 ntr700++;
1487 hdxytof->Fill(dxm, dym);
1488 hpgoodhit->Fill(p); // GeV/c
1489 xc = xhic[himin];
1490 yc = yhic[himin];
1491 zcc = zhic[himin];
1492 dz = zhi[himin] - xyz[2];
1493 xrpc = xyz[0] + dz * tx;
1494 yrpc = xyz[1] + dz * ty;
1495 hxygoodhit->Fill(xgoodtr, ygoodtr);
1496 hxygoodhita->Fill(xhi[himin], yhi[himin]);
1497
1498 hxyhit->Fill(xrpc - xhi[himin], yrpc - yhi[himin]);
1499 hxdiff->Fill(xhi[himin] - xc);
1500 if (idchambers[ic] < 100.f)
1501 hxxhits->Fill(xhi[himin] - xc, xrpc - xc);
1502 if (idchambers[ic] > 100.f)
1503 hxxhitb->Fill(xhi[himin] - xc, xrpc - xc);
1504
1505 cosa = sqrt(1. + tx * tx + ty * ty);
1506 dl = dz * cosa;
1507 lenr = le + dl;
1508 hltra->Fill(lenr);
1509 hidist = TMath::Sqrt(xc * xc + yc * yc + (zcc - ztarg) * (zcc - ztarg));
1510 hstrtr->Fill(is);
1511 hchatr->Fill(ic);
1512 hxytrtof->Fill(xrpc, yrpc);
1513
1514 htoftr->Fill(tof[ic][is]);
1515 //
1516 tofr = tof[ic][is];
1517 tofr0 = tof[ic][is];
1518 tofp = 0.f;
1519 tofpi = 0.f;
1520 beta = lenr / tofr / cvel;
1521 // Selection criteria
1522 if (beta <= 0.)
1523 continue;
1524 mass2 = p * p * (1.f / beta / beta - 1.f);
1525 mass = TMath::Sqrt(mass2);
1526 hpbeta->Fill(fabs(p), beta);
1527 hmass->Fill(mass2);
1528 hmassvsp->Fill(fabs(p), mass2);
1529
1530 massmin1 = MP2 - MP2 / (DMPSCALE1 * PMAX1) * (p - 0.5);
1531 massmax1 = MP2 + MP2 / (DMPSCALE2 * PMAX1) * (p - 0.5);
1532 if (p > 0.f && mass2 > massmin1 && mass2 < massmax1 && pmod > PMIN && pmod < PMAX) {
1533 tofp = (lenr / 29.97925f) * TMath::Sqrt(MP2 + p * p) / p;
1534 tcor = tofp - tofr;
1535 tcorc = tcor;
1536 tcors = tcor;
1537 tcalibrc->Fill(ic, tcorc);
1538 tcalibrs->Fill(ic, is, tcors);
1539 htofp->Fill(tofr0);
1540 massmin2 = MP2 - MP2 / (DMPSCALET1 * PMAX2) * (p - 0.5);
1541 massmax2 = MP2 + MP2 / (DMPSCALET2 * PMAX2) * (p - 0.5);
1542 if (mass2 > massmin2 && mass2 < massmax2) {
1543 hdtpc->Fill(tcorc);
1544 hdtpcvsp->Fill(p, tcorc);
1545 hdtps->Fill(tcors);
1546 hdtpsvsp->Fill(p, tcors);
1547 if (p > PRES)
1548 hdtpcplim->Fill(tcorc);
1549 if (p > PRES1)
1550 hdtpcp1lim->Fill(tcorc);
1551 if (p < 2.0f)
1552 hdtpcp2lim->Fill(tcorc);
1553 if (p >= 2.0f && p < 5.0f)
1554 hdtpcp3lim->Fill(tcorc);
1555 }
1556 }
1557 // massmaxpi = 0.1 + MPIMAX*TMath::Abs(p-0.5)/PMAXPI/DMPISCALE; // Kapishin
1558 massmaxpi = 0.1 + MPIMAX * TMath::Abs(pmod - 0.5) / PMAXPI;
1559 if (p > 0 && massmaxpi > massmin1)
1560 massmaxpi = massmin1;
1561 // massminpi =-0.1 - MPIMAX*TMath::Abs(p-0.5)/PMAXPI/DMPISCALE; // Kapishin
1562 massminpi = -0.1 - MPIMAX * TMath::Abs(pmod - 0.5) / PMAXPI;
1563 if (mass2 < massmaxpi && mass2 > massminpi && pmod > PMIN && pmod < PMAXPI) {
1564 tofpi = (lenr / 29.97925f) * TMath::Sqrt(MPI2 + p * p) / pmod;
1565 tcorpi = tofpi - tofr;
1566 tcorcpi = tcorpi;
1567 tcorspi = tcorpi;
1568 tcalibrc->Fill(ic, tcorcpi);
1569 tcalibrs->Fill(ic, is, tcorspi);
1570 }
1571 tofrc = tofr;
1572 tofrs = tofr;
1573 betac = lenr / tofrc / cvel;
1574 betas = lenr / tofrs / cvel;
1575 masss2 = p * p * (1.f / betas / betas - 1.f);
1576 masss = TMath::Sqrt(masss2);
1577 massc2 = p * p * (1.f / betac / betac - 1.f);
1578 massc = TMath::Sqrt(massc2);
1579 if (p > 0) {
1580 hpbetac->Fill(p, betac);
1581 hmassc->Fill(massc2);
1582 hmasscvsp->Fill(p, massc2);
1583 hpbetas->Fill(p, betas);
1584 hmasss->Fill(masss2);
1585 hmasssvsp->Fill(p, masss2);
1586 if (p < 1.5) {
1587 hmasss15->Fill(masss2);
1588 hmassc15->Fill(massc2);
1589 }
1590 if (p < 2.0) {
1591 hmasss20->Fill(masss2);
1592 hmassc20->Fill(massc2);
1593 }
1594 } else {
1595 hpbetacm->Fill(fabs(p), betac);
1596 hmasscm->Fill(massc2);
1597 hmasscmvsp->Fill(fabs(p), massc2);
1598 hpbetasm->Fill(fabs(p), betas);
1599 hmasssm->Fill(masss2);
1600 hmasssmvsp->Fill(fabs(p), masss2);
1601 }
1602 } // hit is good!
1603 hrej5->Fill(nrej5);
1604 hrej6->Fill(nrej6);
1605 } // hit exist!
1606 Double_t rej1 = (Double_t)nrej1 / nhi;
1607 Double_t rej2 = (Double_t)nrej2 / nhi;
1608 Double_t rej3 = (Double_t)nrej3 / nhi;
1609 hrej1->Fill(rej1);
1610 hrej2->Fill(rej2);
1611 hrej3->Fill(rej3);
1612 } // tracks loop
1613 hntr1->Fill(ntr1);
1614 Ntr1 += ntr1;
1615 Ntr1g += ntr1g;
1616 hntr0->Fill(ntr0);
1617 Ntr0 += ntr0;
1618 Ntr700 += ntr700;
1619 // t_out->Fill();
1620 if (iEv % 50000 == 0) {
1621 hdxyeff->Divide(hdxytof, hdxy);
1622 hdxeff->Divide((projtmp = hdxytof->ProjectionX()), (projtmp1 = hdxy->ProjectionX()));
1623 if (projtmp)
1624 delete projtmp;
1625 if (projtmp1)
1626 delete projtmp1;
1627 hdyeff->Divide((projtmp = hdxytof->ProjectionY()), (projtmp1 = hdxy->ProjectionY()));
1628 if (projtmp)
1629 delete projtmp;
1630 if (projtmp1)
1631 delete projtmp1;
1632 hxyeff->Divide(hxygoodhit, hxygoodtr);
1633 hxeff->Divide((projtmp = hxygoodhit->ProjectionX()), (projtmp1 = hxygoodtr->ProjectionX()));
1634 if (projtmp)
1635 delete projtmp;
1636 if (projtmp1)
1637 delete projtmp1;
1638 hyeff->Divide((projtmp = hxygoodhit->ProjectionY()), (projtmp1 = hxygoodtr->ProjectionY()));
1639 if (projtmp)
1640 delete projtmp;
1641 if (projtmp1)
1642 delete projtmp1;
1643
1644 hxyeffa->Divide(hxygoodhita, hxygoodtr);
1645 hxeffa->Divide((projtmp = hxygoodhita->ProjectionX()), (projtmp1 = hxygoodtr->ProjectionX()));
1646 if (projtmp)
1647 delete projtmp;
1648 if (projtmp1)
1649 delete projtmp1;
1650 hyeffa->Divide((projtmp = hxygoodhita->ProjectionY()), (projtmp1 = hxygoodtr->ProjectionY()));
1651 if (projtmp)
1652 delete projtmp;
1653 if (projtmp1)
1654 delete projtmp1;
1655
1656 hpeff->Divide(hpgoodhit, hpgemdch);
1657 hpeffmg->Divide(hpgemdch, hptr);
1658 hxyeffmg->Divide(hxygoodtrgem, hxygem);
1659 hxyeffmd->Divide(hxygoodtr, hxydch);
1660 hxyeffdg->Divide(hxydch, hxygem);
1661 hxyeffgd->Divide(hxygem, hxydch);
1662 hxyeffm1d->Divide(hxygoodhit1, hxydch);
1663 hxyeffm1da->Divide(hxygoodhit1a, hxydch);
1664 hxybadhit->Add(hxygoodtr, hxygoodhit, 1., -1.);
1665 hxybadhita->Add(hxygoodtr, hxygoodhita, 1., -1.);
1666 hpbadhit->Add(hpgemdch, hpgoodhit, 1., -1.);
1667
1668 h1MCeffmass->Divide(h1MCamass, hMCamass);
1669 h2MCeffmass->Divide(h2MCamass, hMCamass);
1670 h3MCeffmass->Divide(h3MCamass, hMCamass);
1671 h4MCeffmass->Divide(h4MCamass, hMCamass);
1672 hMCeffmass->Divide(hMCmassr, hMCamass);
1673 hMCeffmassp->Divide(hMCmassrp, hMCamass);
1674 hangleff->Divide(hangle1, hangle);
1675 hangleff2->Divide(hangle2, hangle);
1676
1677 // TH1F *hpnohit = new TH1F("hpnohit","Momentum of GEM+DCH+NoTOF700 tracks", 1000,-10,10);
1678
1679 printf("Number of tracks %d, per event %f\n", Ntr, (float)Ntr / nEvents);
1680 printf("GEM+DCH tracks %d, per event %f (%f%% per track)\n", Ntr1, (float)Ntr1 / nEvents,
1681 Ntr > 0 ? (float)Ntr1 / Ntr * 100.f : 0.f);
1682 printf("Good GEM+DCH tracks %d, per event %f (%f%% per valid track)\n", Ntr1g, (float)Ntr1g / nEvents,
1683 Ntr1 > 0 ? (float)Ntr1g / Ntr1 * 100.f : 0.f);
1684 printf("Matched GEM+DCH tracks %d, per event %f (%f%% per good track)\n", Ntr0, (float)Ntr0 / nEvents,
1685 Ntr1g > 0 ? (float)Ntr0 / Ntr1g * 100.f : 0.f);
1686 // printf("Matched GEM+DCH+TOF700 tracks %d, per event %f (%f%% per matched track)\n", Ntr700,
1687 //(float)Ntr700/nEvents, Ntr0>0?(float)Ntr700/Ntr0*100.f:0.f);
1688
1689 if (f) {
1690 f->Write(NULL, TObject::kOverwrite);
1691 }
1692 }
1693 } // event loop
1694
1695 // TH1D *projtmp = 0, *projtmp1 = 0;
1696 hdxyeff->Divide(hdxytof, hdxy);
1697 hdxeff->Divide((projtmp = hdxytof->ProjectionX()), (projtmp1 = hdxy->ProjectionX()));
1698 if (projtmp)
1699 delete projtmp;
1700 if (projtmp1)
1701 delete projtmp1;
1702 hdyeff->Divide((projtmp = hdxytof->ProjectionY()), (projtmp1 = hdxy->ProjectionY()));
1703 if (projtmp)
1704 delete projtmp;
1705 if (projtmp1)
1706 delete projtmp1;
1707 hxyeff->Divide(hxygoodhit, hxygoodtr);
1708 hxeff->Divide((projtmp = hxygoodhit->ProjectionX()), (projtmp1 = hxygoodtr->ProjectionX()));
1709 if (projtmp)
1710 delete projtmp;
1711 if (projtmp1)
1712 delete projtmp1;
1713 hyeff->Divide((projtmp = hxygoodhit->ProjectionY()), (projtmp1 = hxygoodtr->ProjectionY()));
1714 if (projtmp)
1715 delete projtmp;
1716 if (projtmp1)
1717 delete projtmp1;
1718
1719 hxyeffa->Divide(hxygoodhita, hxygoodtr);
1720 hxeffa->Divide((projtmp = hxygoodhita->ProjectionX()), (projtmp1 = hxygoodtr->ProjectionX()));
1721 if (projtmp)
1722 delete projtmp;
1723 if (projtmp1)
1724 delete projtmp1;
1725 hyeffa->Divide((projtmp = hxygoodhita->ProjectionY()), (projtmp1 = hxygoodtr->ProjectionY()));
1726 if (projtmp)
1727 delete projtmp;
1728 if (projtmp1)
1729 delete projtmp1;
1730
1731 hpeff->Divide(hpgoodhit, hpgemdch);
1732 hpeffmg->Divide(hpgemdch, hptr);
1733 hxyeffmg->Divide(hxygoodtrgem, hxygem);
1734 hxyeffmd->Divide(hxygoodtr, hxydch);
1735 hxyeffdg->Divide(hxydch, hxygem);
1736 hxyeffgd->Divide(hxygem, hxydch);
1737 hxyeffm1d->Divide(hxygoodhit1, hxydch);
1738 hxyeffm1da->Divide(hxygoodhit1a, hxydch);
1739 hxybadhit->Add(hxygoodtr, hxygoodhit, 1., -1.);
1740 hxybadhita->Add(hxygoodtr, hxygoodhita, 1., -1.);
1741 hpbadhit->Add(hpgemdch, hpgoodhit, 1., -1.);
1742
1743 h1MCeffmass->Divide(h1MCamass, hMCamass);
1744 h2MCeffmass->Divide(h2MCamass, hMCamass);
1745 h3MCeffmass->Divide(h3MCamass, hMCamass);
1746 h4MCeffmass->Divide(h4MCamass, hMCamass);
1747 hMCeffmass->Divide(hMCmassr, hMCamass);
1748 hMCeffmassp->Divide(hMCmassrp, hMCamass);
1749 hangleff->Divide(hangle1, hangle);
1750 hangleff2->Divide(hangle2, hangle);
1751
1752 // TH1F *hpnohit = new TH1F("hpnohit","Momentum of GEM+DCH+NoTOF700 tracks", 1000,-10,10);
1753
1754 printf("Number of tracks %d, per event %f\n", Ntr, (float)Ntr / nEvents);
1755 printf("GEM+DCH tracks %d, per event %f (%f%% per track)\n", Ntr1, (float)Ntr1 / nEvents,
1756 Ntr > 0 ? (float)Ntr1 / Ntr * 100.f : 0.f);
1757 printf("Good GEM+DCH tracks %d, per event %f (%f%% per valid track)\n", Ntr1g, (float)Ntr1g / nEvents,
1758 Ntr1 > 0 ? (float)Ntr1g / Ntr1 * 100.f : 0.f);
1759 printf("Matched GEM+DCH tracks %d, per event %f (%f%% per good track)\n", Ntr0, (float)Ntr0 / nEvents,
1760 Ntr1g > 0 ? (float)Ntr0 / Ntr1g * 100.f : 0.f);
1761 // printf("Matched GEM+DCH+TOF700 tracks %d, per event %f (%f%% per matched track)\n", Ntr700,
1762 //(float)Ntr700/nEvents, Ntr0>0?(float)Ntr700/Ntr0*100.f:0.f);
1763
1764 if (f) {
1765 f->Write(NULL, TObject::kOverwrite);
1766 }
1767
1768 return;
1769}
1770
1771void select_hist()
1772{
1773 // example of macro called when a mouse event occurs in a pad.
1774 // Example:
1775 // Root > TFile f("hsimple.root");
1776 // Root > hpxpy.Draw();
1777 // Root > c1.AddExec("exsel","select_hist()");
1778 // When moving the mouse in the canvas, a second canvas shows the
1779 // projection along X of the bin corresponding to the Y position
1780 // of the mouse. The resulting histogram is fitted with a gaussian.
1781 // A "dynamic" line shows the current bin position in Y.
1782 // This more elaborated example can be used as a starting point
1783 // to develop more powerful interactive applications exploiting CINT
1784 // as a development engine.
1785 // Author: Rene Brun
1786
1787 int event = gPad->GetEvent();
1788 if (event != 11)
1789 return;
1790
1791 TObject* select = gPad->GetSelected();
1792 if (!select)
1793 return;
1794 if (!(select->InheritsFrom("TH1")))
1795 return;
1796
1797 // create or set the new canvas c2
1798 TCanvas* csingle = 0;
1799 if ((csingle = (TCanvas*)gROOT->GetListOfCanvases()->FindObject("csingle")))
1800 ;
1801 else
1802 csingle = new TCanvas("csingle");
1803 csingle->cd();
1804
1805 // draw slice corresponding to mouse position
1806 if (!(select->InheritsFrom("TH1F"))) {
1807 TH1F* h = (TH1F*)select;
1808 h->Draw();
1809 } else if (!(select->InheritsFrom("TH2F"))) {
1810 TH2F* h = (TH2F*)select;
1811 h->Draw();
1812 } else if (!(select->InheritsFrom("TProfile"))) {
1813 TProfile* h = (TProfile*)select;
1814 h->Draw();
1815 } else if (!(select->InheritsFrom("TH1D"))) {
1816 TH1D* h = (TH1D*)select;
1817 h->Draw();
1818 } else if (!(select->InheritsFrom("TH2D"))) {
1819 TH2D* h = (TH2D*)select;
1820 h->Draw();
1821 } else if (!(select->InheritsFrom("TProfile2D"))) {
1822 TProfile2D* h = (TProfile2D*)select;
1823 h->Draw();
1824 } else
1825 return;
1826 csingle->Update();
1827}
1828
1829int readMCTimeFile(const char* MCTimeFile)
1830{
1831 // cout << "readMCTimeFile: " << MCTimeFile << endl;
1832 // float idchambers[59] = { 27.1,28.1,3.1,1.1,29.1,4.1,33.1,30.1,5.1,19.3,31.1,6.1,2.1,32.1,15.2,16.2,17.2,
1833 // 18.2,19.2,20.2,7.1,115.2,113.1,117.1,35.1,9.1,37.1,11.1,39.1,13.1,34.1,8.1,36.1,10.1,38.1,12.1,21.2,
1834 // 23.2,25.2,22.2,24.2,26.2,107.2,108.2,109.2,110.2,111.2,112.2,114.1,116.2,118.1,14.1,40.1,119.2,120.2,
1835 // 121.2,122.2,123.2,124.2 };
1836 bool notused[59] = {true};
1837 int c = 0, cmin = -1, c0 = 0;
1838 float idmin = 200.f;
1839 for (c0 = 0; c0 < 59; c0++) {
1840 notused[c0] = true;
1841 order[c0] = -1;
1842 }
1843 for (c0 = 0; c0 < 59; c0++) {
1844 cmin = -1;
1845 idmin = 200.f;
1846 for (c = 0; c < 59; c++) {
1847 if (notused[c] && (idchambers[c] < idmin)) {
1848 idmin = idchambers[c];
1849 cmin = c;
1850 }
1851 }
1852 // printf("c0 %d cmin %d notused %d\n",c0,cmin,(int)notused[cmin]);
1853 if (cmin >= 0) {
1854 notused[cmin] = false;
1855 order[c0] = cmin;
1856 }
1857 }
1858 char fname[128];
1859 FILE* ft = 0;
1860 float time = 0.f, timesigma = 0.f;
1861 if (MCTimeFile == NULL) {
1862 printf("TOF700 MC time-of-flight file name not defined!\n");
1863 return 0;
1864 }
1865 if (strlen(MCTimeFile) == 0) {
1866 printf("TOF700 MC time-of-flight file name not defined!\n");
1867 return 0;
1868 }
1869 TString dir = getenv("VMCWORKDIR");
1870 sprintf(fname, "%s/geometry/%s", dir.Data(), MCTimeFile);
1871 ft = fopen(fname, "r");
1872 if (ft == NULL) {
1873 printf("TOF700 MC time-of-flight file %s open error!\n", fname);
1874 return 0;
1875 }
1876 int cexp = 0, ic = -1;
1877 while (fscanf(ft, "%d %f %f\n", &ic, &time, &timesigma) == 3) {
1878 if (ic > 0 && ic <= 59) {
1879 cexp = order[ic - 1];
1880 tofcalmc[cexp] = time;
1881 if (1)
1882 printf("Chamber %.1f (MC %d, EXP %d) average time-of-flight %f with sigma %f\n", idchambers[cexp], ic,
1883 cexp, time, timesigma);
1884 }
1885 }
1886 fclose(ft);
1887 return 1;
1888}
dist
vector< Double_t > dist(vector< Double_t > qp, Double_t mu)
Definition BmnMath.cxx:869
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition P4_F32vec4.h:52
i
int i
Definition P4_F32vec4.h:22
Bool_t
bool Bool_t
Int_t
int Int_t
UInt_t
unsigned int UInt_t
Float_t
float Float_t
Short_t
short Short_t
Double_t
double Double_t
R__ADD_INCLUDE_PATH
#define R__ADD_INCLUDE_PATH(PATH)
gROOT
#define gROOT
gStyle
R__EXTERN TStyle * gStyle
gSystem
R__EXTERN TSystem * gSystem
gPad
#define gPad
kGEM
@ kGEM
Definition BmnDetectorList.h:15
kTOF
@ kTOF
Definition BmnDetectorList.h:18
kREF
@ kREF
Definition BmnDetectorList.h:13
kTOF1
@ kTOF1
Definition BmnDetectorList.h:16
kDCH
@ kDCH
Definition BmnDetectorList.h:17
BmnGemDchTrack::GetDdtx
Double_t GetDdtx() const
Definition BmnGemDchTrack.h:46
BmnGemDchTrack::GetDtx
Double_t GetDtx() const
Definition BmnGemDchTrack.h:42
BmnGemDchTrack::GetYpv
Double_t GetYpv() const
Definition BmnGemDchTrack.h:50
BmnGemDchTrack::GetChi2
Double_t GetChi2() const
Definition BmnGemDchTrack.h:48
BmnGemDchTrack::GetDdx
Double_t GetDdx() const
Definition BmnGemDchTrack.h:44
BmnGemDchTrack::GetSiHits
Int_t GetSiHits() const
Definition BmnGemDchTrack.h:38
BmnGemDchTrack::GetGemTrack
CbmStsTrack * GetGemTrack()
Definition BmnGemDchTrack.h:31
BmnGemDchTrack::GetDdty
Double_t GetDdty() const
Definition BmnGemDchTrack.h:47
BmnGemDchTrack::GetLength
Double_t GetLength() const
Definition BmnGemDchTrack.h:34
BmnGemDchTrack::GetXpv
Double_t GetXpv() const
Definition BmnGemDchTrack.h:49
BmnGemDchTrack::GetDdy
Double_t GetDdy() const
Definition BmnGemDchTrack.h:45
BmnGemDchTrack::GetDx
Double_t GetDx() const
Definition BmnGemDchTrack.h:40
BmnGemDchTrack::GetDty
Double_t GetDty() const
Definition BmnGemDchTrack.h:43
BmnGemDchTrack::GetDy
Double_t GetDy() const
Definition BmnGemDchTrack.h:41
BmnGemDchTrack::GetGemHits
Int_t GetGemHits() const
Definition BmnGemDchTrack.h:37
BmnGemDchTrack::GetDchTrack
BmnTrack * GetDchTrack()
Definition BmnGemDchTrack.h:32
BmnTrack::GetParamLast
FairTrackParam * GetParamLast()
Definition BmnTrack.h:65
CbmMCTrack
CbmMCTrack.h.
Definition CbmMCTrack.h:26
CbmMCTrack::GetPy
Double_t GetPy() const
Definition CbmMCTrack.h:59
CbmMCTrack::GetEnergy
Double_t GetEnergy() const
Definition CbmMCTrack.h:133
CbmMCTrack::GetNPoints
Long64_t GetNPoints(DetectorId detId) const
Accessors to the number of MCPoints in the detectors.
Definition CbmMCTrack.cxx:151
CbmMCTrack::GetPz
Double_t GetPz() const
Definition CbmMCTrack.h:60
CbmMCTrack::GetPx
Double_t GetPx() const
Definition CbmMCTrack.h:58
CbmMCTrack::GetStartZ
Double_t GetStartZ() const
Definition CbmMCTrack.h:63
CbmMCTrack::GetStartY
Double_t GetStartY() const
Definition CbmMCTrack.h:62
CbmMCTrack::GetMotherId
Int_t GetMotherId() const
Definition CbmMCTrack.h:57
CbmMCTrack::GetP
Double_t GetP() const
Definition CbmMCTrack.h:68
CbmMCTrack::GetStartX
Double_t GetStartX() const
Definition CbmMCTrack.h:61
CbmMCTrack::GetPdgCode
Int_t GetPdgCode() const
Accessors.
Definition CbmMCTrack.h:56
CbmStsPoint::GetXIn
Double_t GetXIn() const
Accessors.
Definition CbmStsPoint.h:64
CbmStsPoint::GetZIn
Double_t GetZIn() const
Definition CbmStsPoint.h:66
CbmStsPoint::GetYIn
Double_t GetYIn() const
Definition CbmStsPoint.h:65
CbmStsTrack::GetParamLast
FairTrackParam * GetParamLast()
Definition CbmStsTrack.h:70
CbmTrackMatch
CbmStsTrackMatch.h.
Definition CbmTrackMatch.h:16
CbmTrackMatch::GetNofTrueHits
Int_t GetNofTrueHits() const
Number of true hits on track (from matched MC track).
Definition CbmTrackMatch.h:38
CbmTrackMatch::GetMCTrackId
Int_t GetMCTrackId() const
Index of matched MC track.
Definition CbmTrackMatch.h:35
CbmTrackMatch::GetNofMCTracks
Int_t GetNofMCTracks() const
Number of MCTracks with common hits.
Definition CbmTrackMatch.h:47
CbmTrackMatch::GetNofFakeHits
Int_t GetNofFakeHits() const
Number of fake hits on track (from no MC track).
Definition CbmTrackMatch.h:44
CbmTrackMatch::GetTrueOverAllHitsRatio
Double_t GetTrueOverAllHitsRatio() const
Definition CbmTrackMatch.h:50
CbmTrackMatch::GetNofHits
Int_t GetNofHits() const
Definition CbmTrackMatch.h:49
CbmTrackMatch::GetNofWrongHits
Int_t GetNofWrongHits() const
Number of wrong hits on track (from other MC tracks).
Definition CbmTrackMatch.h:41
CbmVertex
CbmVertex.h.
Definition CbmVertex.h:21
CbmVertex::GetZ
Double_t GetZ() const
Definition CbmVertex.h:60
CbmVertex::GetX
Double_t GetX() const
Accessors.
Definition CbmVertex.h:58
CbmVertex::GetY
Double_t GetY() const
Definition CbmVertex.h:59
TCanvas::cd
TVirtualPad * cd(Int_t subpadnumber=0) override
TCanvas::Update
void Update() override
TChain::SetBranchAddress
Int_t SetBranchAddress(const char *bname, T **add, TBranch **ptr=nullptr)
TChain::Add
virtual Int_t Add(const char *name, Long64_t nentries=TTree::kMaxEntries)
TChain::GetEntries
Long64_t GetEntries() const override
TChain::GetEntry
Int_t GetEntry(Long64_t entry=0, Int_t getall=0) override
TDatabasePDG::Instance
static TDatabasePDG * Instance()
TDatabasePDG::GetParticle
TParticlePDG * GetParticle(const char *name) const
TH1::Add
virtual Bool_t Add(const TH1 *h, const TH1 *h2, Double_t c1=1, Double_t c2=1)
TH1::Draw
void Draw(Option_t *option="") override
TH1::Fill
virtual Int_t Fill(const char *name, Double_t w)
TH1::Divide
virtual Bool_t Divide(const TH1 *h1)
TH1::SetOption
virtual void SetOption(Option_t *option=" ")
TH2::ProjectionY
TH1D * ProjectionY(const char *name="_py", Int_t firstxbin=0, Int_t lastxbin=-1, Option_t *option="") const
TH2::Fill
virtual Int_t Fill(const char *namex, const char *namey, Double_t w)
TH2::ProjectionX
TH1D * ProjectionX(const char *name="_px", Int_t firstybin=0, Int_t lastybin=-1, Option_t *option="") const
TObjArray::GetEntriesFast
Int_t GetEntriesFast() const
TObjArray::At
TObject * At(Int_t idx) const override
TParticlePDG::Charge
Double_t Charge() const
TParticlePDG::Mass
Double_t Mass() const
TProfile2D::Fill
Int_t Fill(const char *namex, const char *namey, Double_t z) override
TProfile::Fill
Int_t Fill(const char *namex, Double_t y) override
TString::Data
const char * Data() const
beta
double beta(double x, double y)
sqrt
RVec< PromoteType< T > > sqrt(const RVec< T > &v)
NUMBER_CHAMBERS
#define NUMBER_CHAMBERS
Definition BmnMCEfficiencyTOF700.C:107
idchambers
float idchambers[59]
Definition BmnMCEfficiencyTOF700.C:123
PMAX2
#define PMAX2
Definition BmnMCEfficiencyTOF700.C:52
nhits
int nhits[NUMBER_CHAMBERS][32]
Definition BmnMCEfficiencyTOF700.C:114
MPIMAX
#define MPIMAX
Definition BmnMCEfficiencyTOF700.C:64
DMPSCALET2
#define DMPSCALET2
Definition BmnMCEfficiencyTOF700.C:88
yhi
double yhi[1000]
Definition BmnMCEfficiencyTOF700.C:116
GEMDEBUG
#define GEMDEBUG
Definition BmnMCEfficiencyTOF700.C:29
XYPANIN
#define XYPANIN
Definition BmnMCEfficiencyTOF700.C:93
DYMAXPAN
#define DYMAXPAN
Definition BmnMCEfficiencyTOF700.C:99
DISTMAX
#define DISTMAX
Definition BmnMCEfficiencyTOF700.C:95
PMAX
#define PMAX
Definition BmnMCEfficiencyTOF700.C:46
halfwidthy
Float_t halfwidthy[NUMBER_CHAMBERS]
Definition BmnMCEfficiencyTOF700.C:110
nhi
int nhi
Definition BmnMCEfficiencyTOF700.C:112
XFACTOR
#define XFACTOR
Definition BmnMCEfficiencyTOF700.C:94
halfwidthx
Float_t halfwidthx[NUMBER_CHAMBERS]
Definition BmnMCEfficiencyTOF700.C:109
zhic
double zhic[1000]
Definition BmnMCEfficiencyTOF700.C:117
MP
#define MP
Definition BmnMCEfficiencyTOF700.C:67
DXMAXPAN
#define DXMAXPAN
Definition BmnMCEfficiencyTOF700.C:98
PRES
#define PRES
Definition BmnMCEfficiencyTOF700.C:55
tof
float tof[NUMBER_CHAMBERS][32]
Definition BmnMCEfficiencyTOF700.C:115
xhi
double xhi[1000]
Definition BmnMCEfficiencyTOF700.C:116
PMIN
#define PMIN
Definition BmnMCEfficiencyTOF700.C:43
DXGEMDCHMAX
#define DXGEMDCHMAX
Definition BmnMCEfficiencyTOF700.C:31
MPI2
#define MPI2
Definition BmnMCEfficiencyTOF700.C:76
yhic
double yhic[1000]
Definition BmnMCEfficiencyTOF700.C:117
DMPSCALE2
#define DMPSCALE2
Definition BmnMCEfficiencyTOF700.C:82
DXMAX
#define DXMAX
Definition BmnMCEfficiencyTOF700.C:96
GEMCUTS
#define GEMCUTS
Definition BmnMCEfficiencyTOF700.C:24
PMAXPI
#define PMAXPI
Definition BmnMCEfficiencyTOF700.C:61
tofcalmc
Double_t tofcalmc[NUMBER_CHAMBERS]
Definition BmnMCEfficiencyTOF700.C:134
zhi
double zhi[1000]
Definition BmnMCEfficiencyTOF700.C:116
PRES1
#define PRES1
Definition BmnMCEfficiencyTOF700.C:58
DMPSCALE1
#define DMPSCALE1
Definition BmnMCEfficiencyTOF700.C:79
MINGEMNDF
#define MINGEMNDF
Definition BmnMCEfficiencyTOF700.C:27
DYMAX
#define DYMAX
Definition BmnMCEfficiencyTOF700.C:97
xhic
double xhic[1000]
Definition BmnMCEfficiencyTOF700.C:117
CFACTOR
#define CFACTOR
Definition BmnMCEfficiencyTOF700.C:37
MP2
#define MP2
Definition BmnMCEfficiencyTOF700.C:70
nhit
int nhit[NUMBER_CHAMBERS]
Definition BmnMCEfficiencyTOF700.C:113
PMAX1
#define PMAX1
Definition BmnMCEfficiencyTOF700.C:49
DYGEMDCHMAX
#define DYGEMDCHMAX
Definition BmnMCEfficiencyTOF700.C:32
TOFMAX
#define TOFMAX
Definition BmnMCEfficiencyTOF700.C:105
DMPSCALET1
#define DMPSCALET1
Definition BmnMCEfficiencyTOF700.C:85
DXBORDER
#define DXBORDER
Definition BmnMCEfficiencyTOF700.C:102
DYBORDER
#define DYBORDER
Definition BmnMCEfficiencyTOF700.C:103
PRIMARY_VERTEX_CUTS
#define PRIMARY_VERTEX_CUTS
Definition BmnMCEfficiencyTOF700.C:25
ACCEPTANCE
#define ACCEPTANCE
Definition BmnMCEfficiencyTOF700.C:101
MINSIHITS
#define MINSIHITS
Definition BmnMCEfficiencyTOF700.C:28
TRACKFILE
#define TRACKFILE
Definition BmnMCEfficiencyTOF700.C:167
MCDEBUG
#define MCDEBUG
Definition BmnMCEfficiencyTOF700.C:30
MCMATCH
#define MCMATCH
Definition BmnMCEfficiencyTOF700.C:21
order
int order[59]
Definition BmnMCEfficiencyTOF700.C:130
MCFILE
#define MCFILE
Definition BmnMCEfficiencyTOF700.C:169
USETRUEHIT
#define USETRUEHIT
Definition BmnMCEfficiencyTOF700.C:18
SORT
#define SORT
Definition BmnMCEfficiencyTOF700.C:19
TOFGEOMETRY
#define TOFGEOMETRY
Definition BmnMCEfficiencyTOF700.C:16
MINTRUE
#define MINTRUE
Definition BmnMCEfficiencyTOF700.C:23
AARGON
#define AARGON
Definition BmnMCEfficiencyTOF700.C:20
MINGEMNDF1
#define MINGEMNDF1
Definition BmnMCEfficiencyTOF700.C:26
readMCTimeFile
Int_t readMCTimeFile(const char *MCTimeFile="TOF700_MC_argon_qgsm_time_run7.txt")
Definition BmnMCEfficiencyTOF700.C:1829
BmnMCEfficiencyTOF700
void BmnMCEfficiencyTOF700(const char *TracksGemDch=TRACKFILE, const char *DataMC=MCFILE)
Definition BmnMCEfficiencyTOF700.C:173
select_hist
void select_hist()
Definition BmnMCEfficiencyTOF700.C:1771
TMath::ACos
Double_t ACos(Double_t)
TMath::IsNaN
Bool_t IsNaN(Double_t x)
TMath::Sqrt
Double_t Sqrt(Double_t x)
TMath::Abs
Double_t Abs(Double_t d)
TMath::Sort
void Sort(Index n, const Element *a, Index *index, Bool_t down=kTRUE)
TMath::RadToDeg
constexpr Double_t RadToDeg()
std
Definition json.hpp:4446