/[PAMELA software]/tof/flight/ToFdEdx_patch/macros/monitor_daily_10th_tri.cxx
ViewVC logotype

Contents of /tof/flight/ToFdEdx_patch/macros/monitor_daily_10th_tri.cxx

Parent Directory Parent Directory | Revision Log Revision Log

Revision 1.1 - (show annotations) (download)
Mon Dec 10 18:47:01 2018 UTC (7 years ago) by mayorov
Branch: MAIN
CVS Tags: HEAD 
ToFdEdx_patch from Wolfgang. Follows Rita Carbone's advanced method to
derive the ToF dEdx (including PMT saturation, Bethe-Bloch effect, etc.)
1 #include <fstream>
2 #include <sstream>
3 #include <string>
4 #include <iostream>
5 #include <TTree.h>
6 #include <TClassEdit.h>
7 #include <TObject.h>
8 #include <TList.h>
9 #include <TSystem.h>
10 #include <TSystemDirectory.h>
11 #include <TString.h>
12 #include <TFile.h>
13 #include <TClass.h>
14 #include <TCanvas.h>
15 #include <TGraphErrors.h>
16 #include <TH1.h>
17 #include <TH1F.h>
18 #include <TH2D.h>
19 #include <TLatex.h>
20 #include <TPad.h>
21 #include <TPaveLabel.h>
22 #include <TChain.h>
23 #include <TSpline.h>
24 #include <TGraph.h>
25 #include <TMultiGraph.h>
26 #include <TLine.h>
27 #include <TStyle.h>
28 #include <TF1.h>
29 #include <TH1F.h>
30 #include <TTree.h>
31 #include <TCanvas.h>
32 #include <PamLevel2.h>
33 #include <RunInfo.h>
34 #include <OrbitalInfo.h>
35 #include <TStopwatch.h>
36 #include <iomanip>
37
38 #include <CaloNuclei.h>
39 #include <ToFNuclei.h>
40 //#include <ToFNaNuclei.h>
41
42 #include </gpfs/wizard/flight/analysis/menn/devel/10th/calo/flight/CaloNuclei/inc/CaloNuclei.h>
43 //#include </gpfs/wizard/flight/analysis/menn/devel/10th/tof/flight/ToFNuclei/inc/ToFNuclei.h>
44 //#include </gpfs/wizard/flight/analysis/menn/devel/10th/tof/flight/ToFNaNuclei/inc/ToFNaNuclei.h>
45 #include </gpfs/wizard/flight/analysis/menn/devel/10th/tof/flight/ToFPatch_v5/inc/ToFPatch.h>
46
47 #include </gpfs/wizard/flight/analysis/menn/devel/10th/tof/flight/ToFdEdx_patch/inc/ToFdEdx_patch.h>
48
49 using namespace std;
50 //
51 //
52
53
54 void monitor(TString dir, TString list, TString outputname, Int_t iyear, TString ioptions, Int_t nmax=500000000){
55
56 const char* alg = "STD"; // or "EXT" "STD" "EXTF" "NUC" "NUCEXT" "NUCEXTF"
57 const char* alg1 = "NUC"; // or "EXT" "STD" "EXTF" "NUC" "NUCEXT" "NUCEXTF"
58
59 const char* tri_or_bi = "TRI"; // gives tri-scale calibartion, "BI" gives bi-scale
60
61 char htitle[50],info[50];
62 char *s = new char[1];
63
64 Double_t tmin,tmax,tint;
65 ULong64_t tmin1,tmax1;
66
67
68 // 2006
69 if (iyear==2006) {
70 tmin = 1152250000 ;
71 tmax = 1167606000 ;
72 }
73
74 // 2007
75 if (iyear==2007) {
76 tmin = 1167606000 ;
77 tmax = 1199142000 ;
78 }
79
80 // 2008
81 if (iyear==2008) {
82 tmin = 1199142000;
83 tmax = 1230764400;
84 }
85
86 // 2009
87 if (iyear==2009) {
88 tmin = 1230764400;
89 tmax = 1262300400;
90 }
91
92 // 2010
93 if (iyear==2010) {
94 tmin = 1262300400;
95 tmax = 1293836400;
96 }
97
98 // 2011
99 if (iyear==2011) {
100 tmin = 1293836400;
101 tmax = 1325372400;
102 }
103
104 // 2012
105 if (iyear==2012) {
106 tmin = 1325372400;
107 tmax = 1356994800;
108 }
109
110 // 2013
111 if (iyear==2013) {
112 tmin = 1356994800;
113 tmax = 1388530800;
114 }
115
116 // 2014
117 if (iyear==2014) {
118 tmin = 1388530800;
119 tmax = 1420066800;
120 }
121
122 // 2015
123 if (iyear==2015) {
124 tmin = 1420066800;
125 tmax = 1451602800;
126 }
127
128 // Andrey definition:
129 if (iyear==0) {
130 // 1.1.2006: 1136070000
131 // 31.12.2015: 1451602800
132 tmin = 1136070000;
133 tmax = 1451602800;
134 }
135
136
137 //tmin = tmin - 10;
138 tmin1=tmin;
139
140 //tmax = tmax + 10;
141 tmax1=tmax;
142
143
144 cout<<"year "<<iyear<<" tmin "<<tmin1<<" tmax "<<tmax1<<endl;
145
146
147 Int_t icount = 0;
148
149 //------------------------------------------------------------------
150
151 gROOT->SetStyle("Plain");
152 gStyle->SetOptStat(000);
153 gStyle->SetOptFit(111);
154 gStyle->SetFuncColor( kRed );
155
156
157 TH2F *sc_adc_he[48];
158
159 TH2F *sc_dedx_p1[48];
160 TH2F *sc_dedx_he1[48];
161 TH2F *sc_dedx_c1[48];
162 TH2F *sc_dedx_p2[48];
163 TH2F *sc_dedx_he2[48];
164 TH2F *sc_dedx_c2[48];
165
166 TH1F *dedx_he_trk = new TH1F("dedx_he_trk","",100.,0.,15.);
167
168 TH1F *adc_he[48];
169 for(int i=0;i<48;i++) {
170 sprintf(htitle,"adc_%d",(i+1));
171 adc_he[i] = new TH1F(htitle,htitle,100,0.,1500.);
172 }
173
174
175 TH1F *dedx_p1[48];
176 for(int i=0;i<48;i++) {
177 sprintf(htitle,"dedx_p1_%d",(i+1));
178 dedx_p1[i] = new TH1F(htitle,htitle,100,0.,5.);
179 }
180
181 TH1F *dedx_he1[48];
182 for(int i=0;i<48;i++) {
183 sprintf(htitle,"dedx_he1_%d",(i+1));
184 dedx_he1[i] = new TH1F(htitle,htitle,100,0.,15.);
185 }
186
187 TH1F *dedx_c1[48];
188 for(int i=0;i<48;i++) {
189 sprintf(htitle,"dedx_c1_%d",(i+1));
190 dedx_c1[i] = new TH1F(htitle,htitle,100,10.,100.);
191 }
192
193 TH1F *dedx_p2[48];
194 for(int i=0;i<48;i++) {
195 sprintf(htitle,"dedx_p2_%d",(i+1));
196 dedx_p2[i] = new TH1F(htitle,htitle,100,0.,5.);
197 }
198
199 TH1F *dedx_he2[48];
200 for(int i=0;i<48;i++) {
201 sprintf(htitle,"dedx_he2_%d",(i+1));
202 dedx_he2[i] = new TH1F(htitle,htitle,100,0.,15.);
203 }
204
205 TH1F *dedx_c2[48];
206 for(int i=0;i<48;i++) {
207 sprintf(htitle,"dedx_c2_%d",(i+1));
208 dedx_c2[i] = new TH1F(htitle,htitle,100,10.,100.);
209 }
210
211
212 TH1F *invbeta[164];
213 for(int i=0;i<164;i++) {
214 sprintf(htitle, "invbeta_%d",(i+1));
215 invbeta[i] = new TH1F(htitle,htitle,100,0.2,1.8);
216 }
217
218 TH1F *pos[24];
219 for(int i=0;i<24;i++) {
220 sprintf(htitle, "pos_%d",(i+1));
221 pos[i] = new TH1F(htitle,htitle,100,-10.,10.);
222 }
223
224
225
226 //------------------------------------------------------------------
227
228
229 // Create a timer object to benchmark this loop
230 TStopwatch timer;
231 timer.Start();
232
233 Float_t fRealTime=0.;
234 Double_t cputime,rtime,rate;
235
236
237 TString wd=gSystem->WorkingDirectory();
238 //
239 // Elena Ischia talk
240
241 PamLevel2* pam_event = new PamLevel2(dir,list,"+AUTO"); // << create pamela event
242
243 // CaloNuclei *cn = new CaloNuclei(pam_event);
244 CaloNuclei *cn = new CaloNuclei(pam_event,alg);
245
246 // ToFNaNuclei *tna = new ToFNaNuclei();
247 // tna->InitPar("/gpfs/wizard/flight/analysis/menn/devel/9th/tof/flight/ToFNaNuclei/GoodParam/");
248
249
250 ToFdEdx_patch *tdedx = new ToFdEdx_patch();
251 tdedx->InitPar("/gpfs/wizard/flight/analysis/menn/devel/10th/tof/flight/ToFdEdx_patch/GoodParam/","flight");
252
253 ToFdEdx_patch *tdedx1 = new ToFdEdx_patch();
254 tdedx1->InitPar("/gpfs/wizard/flight/analysis/menn/devel/10th/tof/flight/ToFdEdx_patch/GoodParam/","flight");
255
256 // ToFPatch *tofpatch = new ToFPatch();
257 // tofpatch->InitPar("/gpfs/wizard/flight/analysis/menn/devel/9th/tof/flight/ToFPatch_v5/GoodParam/");
258
259
260 Int_t ipad_a[6] = {0,8,14,16,18,21};
261 Int_t ipmt_a[6] = {0,16,28,32,36,42};
262
263
264 Int_t nevent = pam_event->GetEntries();
265
266 //--------------------------------------------------------------------
267
268 Int_t nbins = 1500;
269
270 Double_t thelp;
271
272 for(int i=0;i<48;i++) {
273 sprintf(htitle,"sc_adc_%d",(i+1));
274 sc_adc_he[i] = new TH2F(htitle,htitle,nbins,tmin,tmax,100,0.,2000.);
275
276 sprintf(htitle,"sc_dedx_p1_%d",(i+1));
277 sc_dedx_p1[i] = new TH2F(htitle,htitle,nbins,tmin,tmax,100,0.,5.);
278
279 sprintf(htitle,"sc_dedx_he1_%d",(i+1));
280 sc_dedx_he1[i] = new TH2F(htitle,htitle,nbins,tmin,tmax,100,0.,15.);
281
282 sprintf(htitle,"sc_dedx_c1_%d",(i+1));
283 sc_dedx_c1[i] = new TH2F(htitle,htitle,nbins,tmin,tmax,100,10.,80.);
284
285 sprintf(htitle,"sc_dedx_p2_%d",(i+1));
286 sc_dedx_p2[i] = new TH2F(htitle,htitle,nbins,tmin,tmax,100,0.,5.);
287
288 sprintf(htitle,"sc_dedx_he2_%d",(i+1));
289 sc_dedx_he2[i] = new TH2F(htitle,htitle,nbins,tmin,tmax,100,0.,15.);
290
291 sprintf(htitle,"sc_dedx_c2_%d",(i+1));
292 sc_dedx_c2[i] = new TH2F(htitle,htitle,nbins,tmin,tmax,100,10.,80.);
293 }
294
295
296 TH2F *trkdedx_rig1 = new TH2F("trkdedx_rig1","",100,0.1,10,100,0.,50);
297 TH2F *trkdedx_rig2 = new TH2F("trkdedx_rig2","",100,0.1,10,100,0.,50);
298 TH2F *zcalo_rig = new TH2F("zcalo_rig","",100,0.5,10,100,0.,10);
299 TH2F *zcalo_rig1 = new TH2F("zcalo_rig1","",100,0.5,10,100,0.,10);
300
301 //--------------------------------------------------------------------
302
303 Int_t totevent = -1;
304 Int_t trkevent = -1;
305 Int_t p_event = -1;
306 Int_t he_event = -1;
307 Int_t c_event = -1;
308 Int_t invb_event = -1;
309
310 Int_t paddlestrack[6]={0};
311
312 Int_t trigconfsaved = 0;
313
314 Double_t absolute_time;
315 UInt_t atime;
316
317 Double_t tsum = 0.;
318 Double_t nsum = 0.;
319
320 Double_t yhelp,yhelp0,yhelp1,yhelp2,yhelp3;
321 Double_t adcl,adcr,adclo,adcro,am,theta,xkorr,dx,dy,dr;
322 Int_t ihelp,jj,kk,hh;
323
324 Float_t na_dedx[48]; // ToFNaNuclei dEdx
325 Float_t patch_dedx[48]; // ToFPatch dEdx
326 Float_t dedx_wm[48]; // new ToFdEdx_patch
327 Float_t dedx_wm1[48]; // new ToFdEdx_patch
328
329
330 //
331
332 Int_t ntrk1=0;
333 Int_t ntrk2=0;
334 if(nevent < nmax)nmax=nevent;
335 cout << endl<<" Start loop over "<<nmax<<" events "<<endl;
336
337 for (Int_t iev=0; iev<nmax;iev++){
338 //
339
340 pam_event->Clear();
341 pam_event->GetEntry(iev);
342
343 //----- Define time
344 thelp=pam_event->GetOrbitalInfo()->absTime;
345 atime = thelp;
346
347 if (iev==0) tmin = thelp; // first event
348
349 if ( iev == 0) cout<<"First event time : "<<thelp<<endl;
350 // if ( ntrk1%50000 == 0) cout<<iev<<" "<<thelp<<" ntrk "<<ntrk1<<" H "<<p_event<<" He "<<he_event<<" C "<<c_event<<" 1/beta "<<invb_event<<endl;
351 if ( iev%50000 == 0) {
352 timer.Stop();
353 fRealTime = fRealTime + timer.RealTime();
354 timer.Start();
355 rate=(1.*iev)/fRealTime;
356 cout<<" RealTime= "<<fRealTime<<" Rate "<<rate<<" events/sec"<<endl;
357 cout<<iev<<" "<<thelp<<" ntrk "<<ntrk1<<" H "<<p_event<<" He "<<he_event<<" C "<<c_event<<" 1/beta "<<invb_event<<endl;
358 }
359
360
361
362 Int_t ip = 0;
363 Int_t ihe = 0;
364 Int_t ic=0;
365
366
367 Int_t istandard = 1;
368
369 if (istandard == 1) {
370
371
372 //////////////////////////////////////////////////////////
373 ///
374 /// standard
375 ///
376 //////////////////////////////////////////////////////////
377
378
379 // old
380 // int npt = pam_event->GetTrkLevel2()->GetNTracks();
381 // new
382 int npt = pam_event->GetNTracks(alg);
383
384 if( npt>0 ) ntrk1++; //if n.physical tracks>0
385
386 //======================= only single tracks ========================
387
388 if( npt>0 ) { //if n.physical tracks>0
389
390 /*
391 PamTrack *pamtrack = pam_event->GetTrack(0); //this method selects the 0th physical track, discarding the image
392 TrkTrack* trktrack = pamtrack->GetTrkTrack();
393 ToFTrkVar* toftrack = pamtrack->GetToFTrack();
394 CaloTrkVar* calotrack = pamtrack->GetCaloTrack();
395 */
396
397 PamTrack *pamtrack = pam_event->GetTrack(0,alg); //the method selects the 0th physical track, discarding the image
398 ExtTrack* trktrack = pamtrack->GetExtTrack();
399 ToFTrkVar* toftrack = pamtrack->GetToFTrack();
400 CaloTrkVar *calotrack = pam_event->GetTrack(0,alg)->GetCaloTrack();
401
402
403
404
405
406 if(
407 npt ==1 && //the tracks is single
408 trktrack->chi2 > 0 && //the fit converged
409 true){
410
411
412 Float_t chi=trktrack->chi2;
413 Float_t deflection=trktrack->al[4];
414 Float_t rig = 1./deflection;
415 Float_t beta = toftrack->beta[12] ;
416
417 //cout<<"beta "<<beta<<endl;
418
419 //cout<<atime<<" "<<rig<<" "<<"beta "<<beta<<" chi "<<endl;
420
421 if(
422 chi > 0 &&
423 pow((double)chi,0.25) < (3.5+1.85*fabs(deflection)) &&
424 trktrack->GetNX() >= 4 &&
425 trktrack->GetNY() >= 3 &&
426 deflection > 0 &&
427 deflection < 10 &&
428 // beta > 0.1 &&
429 // beta < 2.0 &&
430 // beta > 0.1 &&
431 // beta < 0.6 &&
432 true
433 ){
434
435
436 //cout<<"passing trk cuts "<<endl;
437 //gets(s);
438
439 //cout<<iev<<" STD rig "<<1/deflection<<" chi "<<chi<<" beta "<<beta<<endl;
440
441 //-----------------------------------------------------------------------
442
443
444
445 Double_t def = trktrack->al[4];
446 yhelp0 = 0.5 + 0.6*def*def ;
447 yhelp1 = 4. + 0.8*def*def ;
448 yhelp2 = 4. + 8.*def*def;
449 yhelp3 = 10. + 12.5*def*def;
450
451 yhelp = trktrack->GetDEDX();
452
453 trkdedx_rig1->Fill((1./def),(yhelp/1.3));
454
455 if ((yhelp>yhelp0) && (yhelp<yhelp1)) ip=1; //protons
456 if ((yhelp>yhelp2) && (yhelp<yhelp3)) ihe=1; //helium
457 if (yhelp>yhelp3) ic=1; // Z>2
458
459 //cout<<" rig "<<rig<<" trktrack->GetDEDX() "<<yhelp<<" ip= "<<ip<<endl;
460
461 if (ip==1) p_event = p_event + 1 ;
462 if (ihe==1) he_event = he_event + 1 ;
463
464 if (ihe==1) dedx_he_trk->Fill(yhelp/1.3);
465
466
467 Float_t xv[6],yv[6];
468 for (jj=0; jj<6; jj++){
469 xv[jj]=toftrack->xtr_tof[jj];
470 yv[jj]=toftrack->ytr_tof[jj];
471 }
472
473 Float_t xt,yt,xl,xh,yl,yh;
474 //Float_t beta;
475 Int_t npad,i1,i2;
476
477 Int_t paddleidoftrack[6]= {0} ;
478
479 for (Int_t ilay=0; ilay<6; ilay++) {
480 paddleidoftrack[ilay] = pam_event->GetToFLevel2()->GetPaddleIdOfTrack(xv[ilay], yv[ilay], ilay, 0.0) ;
481 }
482
483 Int_t tof11_i,tof12_i,tof21_i,tof22_i,tof31_i,tof32_i;
484
485 tof11_i = paddleidoftrack[0];
486 tof12_i = paddleidoftrack[1];
487 tof21_i = paddleidoftrack[2];
488 tof22_i = paddleidoftrack[3];
489 tof31_i = paddleidoftrack[4];
490 tof32_i = paddleidoftrack[5];
491
492 //cout<<tof11_i<<" "<<tof12_i<<" "<<tof21_i<<" "<<tof22_i<<" "<<tof31_i<<" "<<tof32_i<<endl;
493
494 //----------------------------------------------------------------
495 //----------------------------------------------------------------
496
497 Float_t xhelp;
498
499 jj = tof11_i ;
500 if (jj != -1) {
501 yhelp = toftrack->ytofpos[0];
502 if (fabs(yhelp-yv[0]) > 0.01) pos[0+jj]->Fill(yhelp-yv[0]);
503 }
504 jj = tof12_i ;
505 if (jj != -1) {
506 xhelp = toftrack->xtofpos[0];
507 if (fabs(xhelp-xv[1]) > 0.01) pos[8+jj]->Fill(xhelp-xv[1]);
508 }
509
510 jj = tof21_i ;
511 if (jj != -1) {
512 xhelp = toftrack->xtofpos[1];
513 if (fabs(xhelp-xv[2]) > 0.01) pos[14+jj]->Fill(xhelp-xv[2]);
514 }
515 jj = tof22_i ;
516 if (jj != -1) {
517 yhelp = toftrack->ytofpos[1];
518 if (fabs(yhelp-yv[3]) > 0.01) pos[16+jj]->Fill(yhelp-yv[3]);
519 }
520
521 jj = tof31_i ;
522 if (jj != -1) {
523 yhelp = toftrack->ytofpos[2];
524 if (fabs(yhelp-yv[4]) > 0.01) pos[18+jj]->Fill(yhelp-yv[4]);
525 }
526 jj = tof32_i ;
527 if (jj != -1) {
528 xhelp = toftrack->xtofpos[2];
529 if (fabs(xhelp-xv[5]) > 0.01) pos[21+jj]->Fill(xhelp-xv[5]);
530 }
531
532 //======================================================================
533 //======== select high rigidity protons to check inverse beta ========
534 //======================================================================
535
536 if ((ip == 1) && (trktrack->al[4] > 0.) && (trktrack->al[4] < 0.08)) {
537
538 invb_event++;
539
540 //----------------------------------------------------------------
541 //--------- S1 - S3 ---------------------------------------------
542 //----------------------------------------------------------------
543
544 //--------- S11 - S31 -------------------------------------------
545
546 i1 = tof11_i;
547 i2 = tof31_i;
548 if ((i1 != -1) && (i2 != -1)) {
549 ihelp=i1*3+i2;
550 beta = toftrack->beta[0] ;
551 invbeta[ihelp]->Fill(1./beta);
552 }
553
554 //--------- S11 - S32 -------------------------------------------
555
556 i1 = tof11_i;
557 i2 = tof32_i;
558 if ((i1 != -1) && (i2 != -1)) {
559 ihelp=24+i1*3+i2;
560 beta = toftrack->beta[1] ;
561 invbeta[ihelp]->Fill(1./beta);
562 }
563
564
565 //--------- S12 - S31 -------------------------------------------
566
567 i1 = tof12_i;
568 i2 = tof31_i;
569 if ((i1 != -1) && (i2 != -1)) {
570 ihelp=48+i1*3+i2;
571 beta = toftrack->beta[2] ;
572 invbeta[ihelp]->Fill(1./beta);
573 }
574
575 //--------- S12 - S32 -------------------------------------------
576
577 i1 = tof12_i;
578 i2 = tof32_i;
579 if ((i1 != -1) && (i2 != -1)) {
580 ihelp=66+i1*3+i2;
581 beta = toftrack->beta[3] ;
582 invbeta[ihelp]->Fill(1./beta);
583 }
584
585
586 //----------------------------------------------------------------
587 //--------- S2 - S3 ---------------------------------------------
588 //----------------------------------------------------------------
589
590 //--------- S21 - S31 -------------------------------------------
591
592 i1 = tof21_i;
593 i2 = tof31_i;
594 if ((i1 != -1) && (i2 != -1)) {
595 ihelp=84+i1*3+i2;
596 beta = toftrack->beta[4] ;
597 invbeta[ihelp]->Fill(1./beta);
598 }
599
600 //--------- S21 - S32 -------------------------------------------
601
602 i1 = tof21_i;
603 i2 = tof32_i;
604 if ((i1 != -1) && (i2 != -1)) {
605 ihelp=90+i1*3+i2;
606 beta = toftrack->beta[5] ;
607 invbeta[ihelp]->Fill(1./beta);
608 }
609
610
611 //--------- S22 - S31 -------------------------------------------
612
613 i1 = tof22_i;
614 i2 = tof31_i;
615 if ((i1 != -1) && (i2 != -1)) {
616 ihelp=96+i1*3+i2;
617 beta = toftrack->beta[6] ;
618 invbeta[ihelp]->Fill(1./beta);
619 }
620
621 //--------- S22 - S32 -------------------------------------------
622
623 i1 = tof22_i;
624 i2 = tof32_i;
625 if ((i1 != -1) && (i2 != -1)) {
626 ihelp=102+i1*3+i2;
627 beta = toftrack->beta[7] ;
628 invbeta[ihelp]->Fill(1./beta);
629 }
630
631
632 //----------------------------------------------------------------
633 //--------- S1 - S2 ---------------------------------------------
634 //----------------------------------------------------------------
635
636 //--------- S11 - S21 -------------------------------------------
637
638 i1 = tof11_i;
639 i2 = tof21_i;
640 if ((i1 != -1) && (i2 != -1)) {
641 ihelp=108+i1*2+i2;
642 beta = toftrack->beta[8] ;
643 invbeta[ihelp]->Fill(1./beta);
644 }
645
646
647 //--------- S11 - S22 -------------------------------------------
648
649 i1 = tof11_i;
650 i2 = tof22_i;
651 if ((i1 != -1) && (i2 != -1)) {
652 ihelp=124+i1*2+i2;
653 beta = toftrack->beta[9] ;
654 invbeta[ihelp]->Fill(1./beta);
655 }
656
657
658 //--------- S12 - S21 -------------------------------------------
659
660 i1 = tof12_i;
661 i2 = tof21_i;
662 if ((i1 != -1) && (i2 != -1)) {
663 ihelp=140+i1*2+i2;
664 beta = toftrack->beta[10] ;
665 invbeta[ihelp]->Fill(1./beta);
666 }
667
668
669 //--------- S12 - S22 -------------------------------------------
670
671 i1 = tof12_i;
672 i2 = tof22_i;
673 if ((i1 != -1) && (i2 != -1)) {
674 ihelp=152+i1*2+i2;
675 beta = toftrack->beta[11] ;
676 invbeta[ihelp]->Fill(1./beta);
677 }
678
679 //-------------------------------------------------------------------------------
680
681 }; // charge and R>12
682
683
684 //=======================================================================
685 //===================== ToFNaNuclei ===============================
686 //=======================================================================
687 /*
688 tna->Process(pam_event,0);
689 for (Int_t ipmt=0; ipmt<48; ipmt++) na_dedx[ipmt] = tna->GetdEdx_pmt(ipmt);
690 */
691 //=======================================================================
692 //===================== ToFdEdx_patch ===============================
693 //=======================================================================
694
695 // tdedx->Process(pam_event,0);
696 tdedx->Process(pam_event,alg,tri_or_bi);
697 for (Int_t ipmt=0; ipmt<48; ipmt++) dedx_wm[ipmt] = tdedx->GetdEdx_pmt(ipmt);
698
699 //=======================================================================
700 //===================== ToFPatch ===============================
701 //=======================================================================
702 // tofpatch->Process(pam_event,0);
703 // for (Int_t ipmt=0; ipmt<48; ipmt++) patch_dedx[ipmt] = tofpatch->GetdEdx_pmt(ipmt);
704 //======================================================================
705 //=========================== select helium ===========================
706 //======================================================================
707
708 if ((ihe==1) && (deflection<0.5)) {
709 //if ((ihe==1)) {
710
711
712 // fill adc "raw" values
713 // we have raw ADC values in then "PMT" class (track independent)
714
715 Float_t adc[48]; //,dedx[48];
716 for (jj=0;jj<48;jj++) adc[jj]=4095.;
717 // for (jj=0;jj<48;jj++) dedx[jj]=1000.;
718
719
720 Float_t aa;
721 Int_t ipm,ich,ihb;
722
723
724 for (Int_t ipmt=0; ipmt<pam_event->GetToFLevel2()->npmt() ; ipmt++){
725 ToFPMT *tofpmt = pam_event->GetToFLevel2()->GetToFPMT(ipmt);
726 ipm = tofpmt->pmt_id;
727 aa = (tofpmt->adc) ;
728 adc_he[ipm]->Fill(aa);
729 sc_adc_he[ipm]->Fill(thelp,aa);
730 adc[ipm] = aa; // raw ADC values
731 //cout<<ipm<<" "<<adc[ipm]<<endl;
732 }
733
734
735 // if (dedx_wm[9]<1000) cout<<adc[9]<<" "<<dedx_wm[9]<<endl;
736
737 // normal loop for 10th red data
738
739 for (Int_t ipmt=0; ipmt<toftrack->npmtadc; ipmt++){
740 // Float_t dEdx = toftrack->dedx[ipmt];
741 Int_t pmtadc = toftrack->pmtadc[ipmt];
742 Int_t adcflag = toftrack->adcflag[ipmt];
743
744 Int_t iverbose=0;
745 // if (pmtadc==14) iverbose=1;
746 // if (pmtadc==27) iverbose=1;
747 // if (pmtadc==18) iverbose=1;
748
749 Float_t dEdx1 = toftrack->dedx[ipmt];
750 // Float_t dEdx1 = na_dedx[pmtadc]; // ToFNaNuclei dEdx
751 if (adcflag==0) dedx_he1[pmtadc]->Fill(dEdx1);
752 if (adcflag==0) sc_dedx_he1[pmtadc]->Fill(thelp,dEdx1);
753
754 }
755
756
757
758 // new loop over hitted paddles for tofdedx_patch data
759
760
761 for (Int_t ilay=0; ilay<6; ilay++) {
762 Int_t i1 = paddleidoftrack[ilay];
763 if (i1 != -1) {
764 Int_t pmtadc = ipmt_a[ilay]+(i1*2);
765
766 Int_t iverbose=0;
767 // if (pmtadc==14) iverbose=1;
768 // if (pmtadc==27) iverbose=1;
769 // if (pmtadc==17) iverbose=1;
770
771 Float_t dEdx2 = dedx_wm[pmtadc]; // ToFdEdx_patch wm
772
773 if (iverbose==1) cout<<atime<<" "<<rig<<" adc "<<adc[pmtadc]<<" tofdedxpatch "<<dEdx2<<endl;
774 if (adc[pmtadc] < 4095) dedx_he2[pmtadc]->Fill(dEdx2);
775 if (adc[pmtadc] < 4095) sc_dedx_he2[pmtadc]->Fill(thelp,dEdx2);
776
777 pmtadc = pmtadc + 1;
778 dEdx2 = dedx_wm[pmtadc]; // ToFdEdx_patch wm
779
780 iverbose=0;
781 // if (pmtadc==14) iverbose=1;
782 // if (pmtadc==27) iverbose=1;
783 // if (pmtadc==17) iverbose=1;
784
785 if (iverbose==1) cout<<atime<<" "<<rig<<" adc "<<adc[pmtadc]<<" tofdedxpatch "<<dEdx2<<endl;
786 if (adc[pmtadc] < 4095) dedx_he2[pmtadc]->Fill(dEdx2);
787 if (adc[pmtadc] < 4095) sc_dedx_he2[pmtadc]->Fill(thelp,dEdx2);
788
789
790 } // i1>-1
791 } // ilay
792
793
794
795 // End fill adc
796
797 }; //ihe==1
798
799
800 //======================================================================
801 //======================================================================
802 //=========== select Protons with Tracker ====================
803 //======================================================================
804
805 //if ( (deflection<0.5) && (ip == 1) && (fabs(charge-1.)<0.15) ) { // iz = 1
806 if ( (deflection<0.5) && (ip == 1) ) { // iz = 1
807
808 //cout<<"Proton detected !!!! "<<endl;
809
810
811 Float_t adc[48];
812 for (jj=0;jj<48;jj++) adc[jj]=4095.;
813
814 Float_t aa;
815 Int_t ipm,ich,ihb;
816
817
818 for (Int_t ipmt=0; ipmt<pam_event->GetToFLevel2()->npmt() ; ipmt++){
819 ToFPMT *tofpmt = pam_event->GetToFLevel2()->GetToFPMT(ipmt);
820 ipm = tofpmt->pmt_id;
821 aa = (tofpmt->adc) ;
822 adc[ipm] = aa; // raw ADC values
823 }
824
825
826 /*
827 for (Int_t ipmt=0; ipmt<toftrack->npmtadc; ipmt++){
828 Int_t ii = toftrack->pmtadc[ipmt];
829 cout<<"In monitor400 "<<ipmt<<" "<<ii<<" "<<adc[ii]<<endl;
830 }
831 */
832
833 // normal loop for 10th reduction data
834
835 for (Int_t ipmt=0; ipmt<toftrack->npmtadc; ipmt++){
836 Int_t pmtadc = toftrack->pmtadc[ipmt];
837
838 Int_t iverbose=0;
839 // if (pmtadc==27) iverbose=1;
840 // if (pmtadc==17) iverbose=1;
841
842 Int_t adcflag = toftrack->adcflag[ipmt];
843
844 Float_t dEdx1 = toftrack->dedx[ipmt];
845 // Float_t dEdx1 = na_dedx[pmtadc]; // ToFNaNuclei dEdx
846
847 if (iverbose==1) cout<<"10th "<<atime<<" "<<rig<<" adc "<<adc[pmtadc]<<" 10th dedx "<<dEdx1<<endl;
848
849 if (adcflag==0) dedx_p1[pmtadc]->Fill(dEdx1);
850 if (adcflag==0) sc_dedx_p1[pmtadc]->Fill(thelp,dEdx1);
851
852
853 if (iverbose==1) gets(s);
854 } // ipmt
855
856
857 // new loop over hitted paddles for tofdedx_patch data
858
859
860 for (Int_t ilay=0; ilay<6; ilay++) {
861 Int_t i1 = paddleidoftrack[ilay];
862 if (i1 != -1) {
863 Int_t pmtadc = ipmt_a[ilay]+(i1*2);
864
865 Int_t iverbose=0;
866 // if (pmtadc==14) iverbose=1;
867 // if (pmtadc==27) iverbose=1;
868 // if (pmtadc==17) iverbose=1;
869
870 Float_t dEdx2 = dedx_wm[pmtadc]; // ToFdEdx_patch wm
871
872 if (iverbose==1) cout<<"tofdedx_patch "<<atime<<" "<<rig<<" adc "<<adc[pmtadc]<<" tofdedxpatch "<<dEdx2<<endl;
873 if (adc[pmtadc] < 4095) dedx_p2[pmtadc]->Fill(dEdx2);
874 if (adc[pmtadc] < 4095) sc_dedx_p2[pmtadc]->Fill(thelp,dEdx2);
875
876 if (iverbose==1) gets(s);
877
878 pmtadc = pmtadc + 1;
879 dEdx2 = dedx_wm[pmtadc]; // ToFdEdx_patch wm
880
881 iverbose=0;
882 // if (pmtadc==14) iverbose=1;
883 // if (pmtadc==27) iverbose=1;
884 // if (pmtadc==17) iverbose=1;
885
886 if (iverbose==1) cout<<"tofdedx_patch "<<atime<<" "<<rig<<" adc "<<adc[pmtadc]<<" tofdedxpatch "<<dEdx2<<endl;
887 if (adc[pmtadc] < 4095) dedx_p2[pmtadc]->Fill(dEdx2);
888 if (adc[pmtadc] < 4095) sc_dedx_p2[pmtadc]->Fill(thelp,dEdx2);
889
890 if (iverbose==1) gets(s);
891
892 } // i1>-1
893 } // ilay
894
895
896
897
898 } // protons selected
899
900 //======================================================================
901
902 }; // tracking cut
903 }; // fit converged
904 }; // ntrk = 1 or 2
905
906 //======================================================================
907
908
909 } // istandard
910
911
912
913 Int_t inuc = 1;
914
915 if (inuc == 1) {
916
917 //////////////////////////////////////////////////////////
918 ///
919 /// NUC
920 ///
921 //////////////////////////////////////////////////////////
922
923 int npt = pam_event->GetNTracks(alg1);
924
925 if( npt>0 ) ntrk2++; //if n.physical tracks>0
926
927 //======================= only single tracks ========================
928
929 if( npt>0 ) { //if n.physical tracks>0
930
931 PamTrack *pamtrack1 = pam_event->GetTrack(0,alg1); //the method selects the 0th physical track, discarding the image
932
933 ExtTrack* trktrack1 = pamtrack1->GetExtTrack();
934 ToFTrkVar* toftrack1 = pamtrack1->GetToFTrack();
935 CaloTrkVar *calotrack1 = pam_event->GetTrack(0,alg1)->GetCaloTrack();
936
937
938 if(
939 npt ==1 && //the tracks is single
940 trktrack1->chi2 > 0 && //the fit converged
941 true){
942
943
944 //
945
946 Float_t chi=trktrack1->chi2;
947 Float_t deflection=trktrack1->al[4];
948 Float_t rig = 1./deflection;
949 Float_t beta = toftrack1->beta[12] ;
950
951 if(
952 chi > 0 &&
953 // pow((double)chi,0.25) < (3.5+1.85*fabs(deflection)) &&
954 trktrack1->GetNX() >= 3 &&
955 trktrack1->GetNY() >= 2 &&
956 trktrack1->GetLeverArmX() >= 4 &&
957 deflection > 0 &&
958 deflection < 10 &&
959 // beta > 0.1 &&
960 // beta < 0.6 &&
961 true
962 ){
963
964 //cout<<iev<<" NUC rig "<<1/deflection<<" chi "<<chi<<" beta "<<beta<<endl;
965
966
967 Double_t def = trktrack1->al[4];
968 yhelp0 = 0.5 + 0.6*def*def ;
969 yhelp1 = 4. + 0.8*def*def ;
970 yhelp2 = 4. + 8.*def*def;
971 yhelp3 = 10. + 12.5*def*def;
972
973 yhelp = trktrack1->GetDEDX();
974
975 // trkdedx_rig1->Fill((1./def),(yhelp/1.3));
976
977 if ((yhelp>yhelp0) && (yhelp<yhelp1)) ip=1; //protons
978 if ((yhelp>yhelp2) && (yhelp<yhelp3)) ihe=1; //helium
979 if (yhelp>yhelp3) ic=1; // Z>2
980
981 if (ip==1) p_event = p_event + 1 ;
982 if (ihe==1) he_event = he_event + 1 ;
983
984 // if (ihe==1) dedx_he_trk->Fill(yhelp/1.3);
985
986
987
988
989 Float_t xv[6],yv[6];
990 for (jj=0; jj<6; jj++){
991 xv[jj]=toftrack1->xtr_tof[jj];
992 yv[jj]=toftrack1->ytr_tof[jj];
993 }
994
995 Int_t paddleidoftrack[6]= {0} ;
996
997 for (Int_t ilay=0; ilay<6; ilay++) {
998 paddleidoftrack[ilay] = pam_event->GetToFLevel2()->GetPaddleIdOfTrack(xv[ilay], yv[ilay], ilay, 0.0) ;
999 }
1000
1001
1002 //======================================================================
1003 //================= select Carbon with Calo ========================
1004 //======================================================================
1005
1006
1007 Float_t charge = 1000.;
1008 Float_t mip=0;
1009
1010 //Float_t charge1 = 1000;
1011 //charge1 = cn->Get_charge_siegen1();
1012
1013 // mip=cn->Get_StdEdx1();
1014
1015
1016 charge=1000;
1017
1018 // Data from file Calo_Bands_New_7.dat
1019 Float_t C0[9] = {0 , 1 , 2 , 3 , 4 , 5 , 6 , 8 , 90 };
1020 Float_t B0[9] = {0 , -2.03769 , 7.61781 , 19.7098 , 60.5598 , 57.9226 , 14.8368 , -1358.83 , 8200 };
1021 Float_t B1[9] = {0 , 0.0211274 , 9.32057e-010 , 4.47241e-07 , 1.44826e-06 , 2.6189e-05 , 0.00278178 , 55.5445 , 0 };
1022 Float_t B2[9] = {0 , -3.91742 , -20.0359 , -16.3043 , -16.9471 , -14.4622 , -10.9594 , -2.38014 , 0 };
1023 Float_t B3[9] = {0 , 11.1469 , -6.63105 , -27.8834 , -132.044 , -55.341 , 173.25 , 4115 , 0 };
1024 Float_t B4[9] = {0 , -14.3465 , -0.485215 , 18.8122 , 117.533 , -14.0898 , -325.269 , -4388.89 , 0 };
1025 Float_t B5[9] = {0 , 6.24281 , 3.96018 , 0 , -26.1881 , 42.9731 , 182.697 , 1661.01 , 0 };
1026
1027
1028 Float_t x1[9],y1[9];
1029 Int_t n1 = 9;
1030
1031 // if ( !calotrack ) printf(" ERROR: No CaloTrack!\n");
1032
1033 Int_t view = 0;
1034 Int_t plane = 0;
1035 Int_t strip = 0;
1036
1037 Int_t indx = 0;
1038 Float_t vfpl[96];
1039 Int_t stfpl[96];
1040 memset(vfpl, 0, 96*sizeof(Float_t));
1041 memset(stfpl, 0, 96*sizeof(Int_t));
1042
1043 Float_t dedx1 = 0.;
1044 Float_t stdedx1 = 0.;
1045 Float_t maxrel = 0.;
1046
1047 mip=0.;
1048
1049 for ( Int_t is=0; is<pam_event->GetCaloLevel1()->istrip; is++ ){
1050 //
1051 mip = pam_event->GetCaloLevel1()->DecodeEstrip(is,view,plane,strip);
1052
1053 if ( strip != -1 &&
1054 view == 1 &&
1055 plane == 0 &&
1056 ( strip == (calotrack1->tibar[0][1]-1) || strip == (calotrack1->tibar[0][1]-2) || strip == (calotrack1->tibar[0][1]) )
1057 && true ){
1058 dedx1 += mip;
1059 };
1060
1061 if ( strip != -1 && view == 1 && plane == 0 ) {
1062 stfpl[indx] = strip;
1063 vfpl[indx] = mip;
1064 indx++;
1065 };
1066
1067
1068 } // is
1069
1070
1071 // find energy released along the strip of maximum on the first plane and on the two neighbour strips
1072 //
1073 if ( indx > 0 ){
1074 Int_t mindx = (Int_t)TMath::LocMax(indx,vfpl);
1075 for (Int_t ii=0; ii<indx; ii++){
1076 if ( stfpl[ii] == stfpl[mindx] ) stdedx1 += vfpl[ii];
1077 if ( (mindx-1)>=0 && stfpl[ii] == (stfpl[mindx]-1) ) stdedx1 += vfpl[ii];
1078 if ( (mindx+1)<96 && stfpl[ii] == (stfpl[mindx]+1) ) stdedx1 += vfpl[ii];
1079 // if ( (mindx-1)>=0 && stfpl[ii] == stfpl[mindx-1] ) stdedx1 += vfpl[ii];
1080 // if ( (mindx+1)<96 && stfpl[ii] == stfpl[mindx+1] ) stdedx1 += vfpl[ii];
1081 };
1082 maxrel = vfpl[mindx];
1083 } else {
1084 stdedx1 = 0.;
1085 maxrel = 0.;
1086 };
1087
1088
1089 // cout<<"cn->Get_stdedx1 "<<mip1<<" selfmade dedx1 "<<dedx1<<" selfmade stdedx1 "<<stdedx1<<"\n";
1090 //
1091 //cout<<"Calo dedx1 "<<dedx1<<endl;
1092
1093 //-----------------------------------------------------------------------
1094
1095
1096 beta = toftrack1->beta[12] ;
1097
1098 // if (beta<2.) { // it makes no sense to allow beta=5 or so...
1099 if ((beta>0.1)&&(beta<2.)) { // it makes no sense to allow beta=5 or so... reject negative beta
1100
1101 //------ Calo function
1102
1103 Float_t mip;
1104 // mip=dedx1;
1105 // mip=cn->Get_StdEdx1();
1106 mip = stdedx1 ;
1107
1108
1109 if (mip>0) {
1110
1111 Float_t betahelp = pow(beta, 1.8);
1112 Float_t ym = mip*betahelp;
1113 Float_t xb = beta;
1114
1115 for ( Int_t jj=0; jj<9; jj++ ){
1116 x1[jj] = B0[jj]+B1[jj]*pow(xb,B2[jj])+B3[jj]*xb+B4[jj]*xb*xb+B5[jj]*xb*xb*xb;
1117 y1[jj] = C0[jj]*C0[jj] ;
1118 }
1119
1120 TGraph *gr1 = new TGraph(n1,x1,y1);
1121 TSpline3 *spl1 = new TSpline3("grs",gr1); // use a cubic spline
1122 Float_t chelp = spl1->Eval(ym);
1123 charge = TMath::Sqrt(chelp);
1124 gr1->Delete();
1125 spl1->Delete();
1126
1127 } // if (mip1>0)
1128 } // beta < 2
1129
1130 //cout<<"Handmade charge "<<charge<<endl;
1131 //cout<<"Handmade : stdedx1 "<<stdedx1<<" charge "<<charge<<endl;
1132
1133 //cout<<"Handmade charge "<<charge<<" new CaloNuclei "<<charge1<<endl;
1134
1135
1136 //cout<<iev<<" rig "<<(1./def)<<" beta "<<beta<<" stdedx1 "<<stdedx1<<" charge "<<charge<<endl;
1137
1138 //=================================================================
1139
1140 //cout<<"CaloNuclei : stdedx1 "<<mip<<" charge "<<charge<<endl;
1141
1142 zcalo_rig->Fill((1./def),charge);
1143
1144 if (ic==1) { // otherwise the count is swamped by background of protons...
1145
1146 zcalo_rig1->Fill((1./def),charge);
1147
1148 ic=0;
1149 //if (fabs(charge-6.)<0.40) ic = 1;
1150 if ((5.6<charge)&&(charge<6.4)) ic = 1;
1151
1152 if (ic==1) c_event = c_event + 1 ;
1153
1154 //=======================================================================
1155 //=======================================================================
1156 //=======================================================================
1157
1158 if (fabs(charge-6.)<0.40) { // iz = 6
1159 if (deflection<0.5) {
1160
1161
1162 // cout<<" main program: npt = pam_event->GetNTracks(alg1) "<<npt<<endl;
1163
1164 //===================== ToFdEdx_patch ==========================
1165
1166 tdedx1->Process(pam_event,alg1);
1167 for (Int_t ipmt=0; ipmt<48; ipmt++) dedx_wm1[ipmt] = tdedx1->GetdEdx_pmt(ipmt);
1168
1169 //=======================================================================
1170
1171 Float_t adc[48];
1172 for (jj=0;jj<48;jj++) adc[jj]=4095.;
1173
1174 Float_t aa;
1175 Int_t ipm,ich,ihb;
1176
1177
1178 for (Int_t ipmt=0; ipmt<pam_event->GetToFLevel2()->npmt() ; ipmt++){
1179 ToFPMT *tofpmt = pam_event->GetToFLevel2()->GetToFPMT(ipmt);
1180 ipm = tofpmt->pmt_id;
1181 aa = (tofpmt->adc) ;
1182 adc[ipm] = aa; // raw ADC values
1183 }
1184
1185
1186
1187 for (Int_t ipmt=0; ipmt<toftrack1->npmtadc; ipmt++){
1188 Int_t pmtadc = toftrack1->pmtadc[ipmt];
1189 // Float_t dEdx = toftrack1->dedx[ipmt];
1190 // Float_t dEdx10th = toftrack->dedx[ipmt];
1191 // cout<<ipmt<<" 10th "<<dEdx10th<<" TNa "<<dEdx<<endl;
1192 Int_t adcflag = toftrack1->adcflag[ipmt];
1193
1194 Float_t dEdx1 = toftrack1->dedx[ipmt];
1195 // Float_t dEdx1 = na_dedx[pmtadc]; // ToFNaNuclei dEdx
1196
1197 Int_t iverbose=0;
1198 // if (pmtadc==33) iverbose=1;
1199 if (iverbose==1) cout<<ipmt<<" "<<pmtadc<<" adc "<<adc[pmtadc]<<" 10th "<<dEdx1<<endl;
1200
1201 if (adcflag==0) dedx_c1[pmtadc]->Fill(dEdx1);
1202 if (adcflag==0) sc_dedx_c1[pmtadc]->Fill(thelp,dEdx1);
1203
1204 }
1205
1206
1207
1208 // new loop over hitted paddles for tofdedx_patch data
1209
1210
1211 for (Int_t ilay=0; ilay<6; ilay++) {
1212 Int_t i1 = paddleidoftrack[ilay];
1213 if (i1 != -1) {
1214 Int_t pmtadc = ipmt_a[ilay]+(i1*2);
1215
1216 Int_t iverbose=0;
1217 // if (pmtadc==14) iverbose=1;
1218 // if (pmtadc==27) iverbose=1;
1219 // if (pmtadc==33) iverbose=1;
1220
1221 Float_t dEdx2 = dedx_wm1[pmtadc]; // ToFdEdx_patch wm
1222
1223 // if (iverbose==1) cout<<atime<<" "<<rig<<" adc "<<adc[pmtadc]<<" tofdedxpatch "<<dEdx2<<endl;
1224 if (iverbose==1) cout<<pmtadc<<" adc "<<adc[pmtadc]<<" tofdedxpatch "<<dEdx2<<endl;
1225 if (adc[pmtadc] < 4095) dedx_c2[pmtadc]->Fill(dEdx2);
1226 if (adc[pmtadc] < 4095) sc_dedx_c2[pmtadc]->Fill(thelp,dEdx2);
1227
1228 pmtadc = pmtadc + 1;
1229 dEdx2 = dedx_wm1[pmtadc]; // ToFdEdx_patch wm
1230
1231 iverbose=0;
1232 // if (pmtadc==14) iverbose=1;
1233 // if (pmtadc==27) iverbose=1;
1234 // if (pmtadc==33) iverbose=1;
1235
1236 // cout<<pmtadc<<" adc "<<adc[pmtadc]<<" tofdedxpatch "<<dEdx2<<endl;
1237 if (iverbose==1) cout<<pmtadc<<" adc "<<adc[pmtadc]<<" tofdedxpatch "<<dEdx2<<endl;
1238
1239 // if (iverbose==1) cout<<atime<<" "<<rig<<" adc "<<adc[pmtadc]<<" tofdedxpatch "<<dEdx2<<endl;
1240 if (adc[pmtadc] < 4095) dedx_c2[pmtadc]->Fill(dEdx2);
1241 if (adc[pmtadc] < 4095) sc_dedx_c2[pmtadc]->Fill(thelp,dEdx2);
1242
1243
1244 } // i1>-1
1245 } // ilay
1246
1247
1248 } // deflection < 0.5
1249 } // charge-6 < 0.4
1250
1251 //=======================================================================
1252 //=======================================================================
1253 //=======================================================================
1254
1255 }
1256
1257
1258
1259 //======================================================================
1260
1261 }; // tracking cut
1262 }; // fit converged
1263 }; // ntrk = 1 or 2
1264
1265 //======================================================================
1266
1267 } // inuc
1268
1269 //======================================================================
1270
1271
1272
1273
1274 thelp=pam_event->GetOrbitalInfo()->absTime;
1275
1276 tsum = tsum + thelp ;
1277 nsum = nsum + 1. ;
1278
1279
1280 //---------------------------------------------------------------------
1281 //---------------------------------------------------------------------
1282
1283 }; // end loop over the events
1284
1285 //---------------------------------------------------------------------
1286 //---------------------------------------------------------------------
1287
1288 cout << endl << endl << " Done "<< endl<<endl;
1289 cout <<"Std " << ntrk1 <<" tracks over "<<nmax<<" events ("<< 100*ntrk1/nmax<<"%)"<<endl;
1290 cout <<"NUC " << ntrk2 <<" tracks over "<<nmax<<" events ("<< 100*ntrk2/nmax<<"%)"<<endl;
1291 cout<<" H "<<p_event<<" He "<<he_event<<" C "<<c_event<<" 1/beta "<<invb_event<<endl;
1292
1293
1294 // Stop timer and print results
1295 timer.Stop();
1296 rtime = timer.RealTime();
1297 cputime = timer.CpuTime();
1298 printf("RealTime=%f seconds, CpuTime=%f seconds\n",rtime,cputime);
1299
1300
1301 //================================================================================
1302 //================================================================================
1303
1304 Double_t xhelp = tsum/nsum ; // mean time in interval
1305
1306
1307 cout<<" mean time in interval "<<xhelp<<endl;
1308
1309
1310
1311 Float_t XX[1],YY[1],eXX[1],eYY[1];
1312
1313 XX[0] = xhelp;
1314 eXX[0] = 0.00001;
1315
1316 // TGraphErrors *gr_ibeta[164];
1317 // TGraphErrors *gr_pos[24];
1318
1319 TGraphErrors *gr_adc[48];
1320
1321 TGraphErrors *gr_dedx_p1[48];
1322 TGraphErrors *gr_dedx_he1[48];
1323 TGraphErrors *gr_dedx_c1[48];
1324
1325 TGraphErrors *gr_dedx_p2[48];
1326 TGraphErrors *gr_dedx_he2[48];
1327 TGraphErrors *gr_dedx_c2[48];
1328
1329
1330
1331 TF1 *f1 = new TF1("f1","landau",0,100);
1332 //TF1 *f1 = new TF1("f1","gaus",0.,100.);
1333
1334 for (jj=0;jj<48;jj++) {
1335
1336 Float_t xhelp1 = adc_he[jj]->GetMean();
1337 Float_t xhelp2 = 0.01;
1338 YY[0] = xhelp1;
1339 eYY[0] = xhelp2;
1340 gr_adc[jj] = new TGraphErrors(1,XX,YY,eXX,eYY);
1341
1342 f1->SetParameters(0.2*dedx_p1[jj]->GetEntries(), dedx_p1[jj]->GetMean(), 0.5*dedx_p1[jj]->GetRMS());
1343 f1->SetRange(0.5*dedx_p1[jj]->GetMean(), 2.0*dedx_p1[jj]->GetMean());
1344 dedx_p1[jj]->Fit("f1", "R");
1345 Float_t xhelp3 = f1->GetParameter(1);
1346 Float_t xhelp4 = f1->GetParError(1);
1347 YY[0] = xhelp3;
1348 eYY[0] = xhelp4;
1349 gr_dedx_p1[jj] = new TGraphErrors(1,XX,YY,eXX,eYY);
1350
1351 f1->SetParameters(0.2*dedx_he1[jj]->GetEntries(), dedx_he1[jj]->GetMean(), 0.5*dedx_he1[jj]->GetRMS());
1352 f1->SetRange(0.5*dedx_he1[jj]->GetMean(), 2.0*dedx_he1[jj]->GetMean());
1353 dedx_he1[jj]->Fit("f1", "R");
1354 Float_t xhelp5 = f1->GetParameter(1);
1355 Float_t xhelp6 = f1->GetParError(1);
1356 YY[0] = xhelp5;
1357 eYY[0] = xhelp6;
1358 gr_dedx_he1[jj] = new TGraphErrors(1,XX,YY,eXX,eYY);
1359
1360 f1->SetParameters(0.2*dedx_c1[jj]->GetEntries(), dedx_c1[jj]->GetMean(), 0.5*dedx_c1[jj]->GetRMS());
1361 f1->SetRange(0.5*dedx_c1[jj]->GetMean(), 2.0*dedx_c1[jj]->GetMean());
1362 dedx_c1[jj]->Fit("f1", "R");
1363 Float_t xhelp7 = f1->GetParameter(1);
1364 Float_t xhelp8 = f1->GetParError(1);
1365 YY[0] = xhelp7;
1366 eYY[0] = xhelp8;
1367 gr_dedx_c1[jj] = new TGraphErrors(1,XX,YY,eXX,eYY);
1368
1369
1370 f1->SetParameters(0.2*dedx_p2[jj]->GetEntries(), dedx_p2[jj]->GetMean(), 0.5*dedx_p2[jj]->GetRMS());
1371 f1->SetRange(0.5*dedx_p2[jj]->GetMean(), 2.0*dedx_p2[jj]->GetMean());
1372 dedx_p2[jj]->Fit("f1", "R");
1373 xhelp3 = f1->GetParameter(1);
1374 xhelp4 = f1->GetParError(1);
1375 YY[0] = xhelp3;
1376 eYY[0] = xhelp4;
1377 gr_dedx_p2[jj] = new TGraphErrors(1,XX,YY,eXX,eYY);
1378
1379 f1->SetParameters(0.2*dedx_he2[jj]->GetEntries(), dedx_he2[jj]->GetMean(), 0.5*dedx_he2[jj]->GetRMS());
1380 f1->SetRange(0.5*dedx_he2[jj]->GetMean(), 2.0*dedx_he2[jj]->GetMean());
1381 dedx_he2[jj]->Fit("f1", "R");
1382 xhelp5 = f1->GetParameter(1);
1383 xhelp6 = f1->GetParError(1);
1384 YY[0] = xhelp5;
1385 eYY[0] = xhelp6;
1386 gr_dedx_he2[jj] = new TGraphErrors(1,XX,YY,eXX,eYY);
1387
1388 f1->SetParameters(0.2*dedx_c2[jj]->GetEntries(), dedx_c2[jj]->GetMean(), 0.5*dedx_c2[jj]->GetRMS());
1389 f1->SetRange(0.5*dedx_c2[jj]->GetMean(), 2.0*dedx_c2[jj]->GetMean());
1390 dedx_c2[jj]->Fit("f1", "R");
1391 xhelp7 = f1->GetParameter(1);
1392 xhelp8 = f1->GetParError(1);
1393 YY[0] = xhelp7;
1394 eYY[0] = xhelp8;
1395 gr_dedx_c2[jj] = new TGraphErrors(1,XX,YY,eXX,eYY);
1396
1397
1398 // cout<<jj<<" "<<xhelp1<<" "<<xhelp2<<" "<<xhelp3<<" "<<xhelp4<<" "<<xhelp5<<" "<<xhelp6<<" "<<xhelp7<<" "<<xhelp8<<"\n";
1399 // fout1<<jj<<" "<<xhelp1<<" "<<xhelp2<<" "<<xhelp3<<" "<<xhelp4<<" "<<xhelp5<<" "<<xhelp6<<" "<<xhelp7<<" "<<xhelp8<<"\n";
1400 }
1401
1402 /*
1403 TF1 *f2 = new TF1("f2","gaus",0.3,1.6);
1404
1405
1406 for (jj=0;jj<164;jj++) {
1407 invbeta[jj]->Fit("f2");
1408 Float_t xhelp5 = f2->GetParameter(1);
1409 Float_t xhelp6 = f2->GetParError(1);
1410 Float_t xhelp7 = f2->GetParameter(2);
1411 cout<<jj<<" "<<xhelp5<<" "<<xhelp6<<" "<<xhelp7<<"\n";
1412 YY[0] = xhelp5;
1413 eYY[0] = xhelp6;
1414 gr_ibeta[jj] = new TGraphErrors(1,XX,YY,eXX,eYY);
1415 }
1416
1417 for (jj=0;jj<24;jj++) {
1418 pos[jj]->Fit("f2");
1419 Float_t xhelp8 = f2->GetParameter(1);
1420 Float_t xhelp9 = f2->GetParError(1);
1421 Float_t xhelp10 = f2->GetParameter(2);
1422 cout<<jj<<" "<<xhelp8<<" "<<xhelp9<<" "<<xhelp10<<"\n";
1423 YY[0] = xhelp8;
1424 eYY[0] = xhelp9;
1425 gr_pos[jj] = new TGraphErrors(1,XX,YY,eXX,eYY);
1426 }
1427 */
1428
1429 //================================================================================
1430 //================================================================================
1431
1432 cout<<"Finished gr-calculation"<<endl;
1433 //cout<<" mean time in interval "<<xhelp<<endl;
1434
1435 // -------------------
1436 // create output files
1437 // -------------------
1438
1439 TFile *outfh = new TFile(outputname,"RECREATE");
1440
1441 for(int i=0;i<48;i++) {
1442 sc_adc_he[i]->Write();
1443 sc_dedx_p1[i]->Write();
1444 sc_dedx_he1[i]->Write();
1445 sc_dedx_c1[i]->Write();
1446 sc_dedx_p2[i]->Write();
1447 sc_dedx_he2[i]->Write();
1448 sc_dedx_c2[i]->Write();
1449 }
1450
1451 /*
1452 for(int i=0;i<164;i++) {
1453 sprintf(info,"gr_ibeta_%d",i);
1454 gr_ibeta[i]->Write(info);
1455 }
1456
1457 for(int i=0;i<24;i++) {
1458 sprintf(info,"gr_pos_%d",i);
1459 gr_pos[i]->Write(info);
1460 }
1461 */
1462
1463 for(int i=0;i<48;i++) {
1464 sprintf(info,"gr_adc_%d",i+1);
1465 // cout<<info<<endl;
1466 gr_adc[i]->Write(info);
1467
1468 sprintf(info,"gr_dedx_p1_%d",i+1);
1469 // cout<<info<<endl;
1470 gr_dedx_p1[i]->Write(info);
1471 sprintf(info,"gr_dedx_he1_%d",i+1);
1472 // cout<<info<<endl;
1473 gr_dedx_he1[i]->Write(info);
1474 sprintf(info,"gr_dedx_c1_%d",i+1);
1475 // cout<<info<<endl;
1476 gr_dedx_c1[i]->Write(info);
1477 sprintf(info,"gr_dedx_p2_%d",i+1);
1478 // cout<<info<<endl;
1479 gr_dedx_p2[i]->Write(info);
1480 sprintf(info,"gr_dedx_he2_%d",i+1);
1481 // cout<<info<<endl;
1482 gr_dedx_he2[i]->Write(info);
1483 sprintf(info,"gr_dedx_c2_%d",i+1);
1484 // cout<<info<<endl;
1485 gr_dedx_c2[i]->Write(info);
1486
1487 }
1488
1489
1490 trkdedx_rig1->Write();
1491 trkdedx_rig2->Write();
1492 zcalo_rig->Write();
1493 zcalo_rig1->Write();
1494
1495 outfh->Close();
1496
1497
1498 // end WM
1499
1500
1501
1502 gSystem->ChangeDirectory(wd);
1503
1504 };
1505
1506
  ViewVC Help
Powered by ViewVC 1.1.23