/[PAMELA software]/tof/flight/ToFReprocessing/ToFreproc.cpp
ViewVC logotype

Contents of /tof/flight/ToFReprocessing/ToFreproc.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.7 - (show annotations) (download)
Fri Jan 17 15:16:24 2014 UTC (10 years, 10 months ago) by mocchiut
Branch: MAIN
CVS Tags: HEAD
Changes since 1.6: +10 -8 lines
Compilation warnings using GCC4.7 fixed

1 #if !defined(__CINT__) || defined(__MAKECINT__)
2
3
4 #include <fstream>
5 #include <sstream>
6 #include <string>
7 #include <math.h>
8 #include <iostream>
9 #include <TTree.h>
10 #include <TClassEdit.h>
11 #include <TObject.h>
12 #include <TList.h>
13 #include <TSystem.h>
14 #include <TSystemDirectory.h>
15 #include <TString.h>
16 #include <TFile.h>
17 #include <TClass.h>
18 #include <TCanvas.h>
19 #include <TH1.h>
20 #include <TH1F.h>
21 #include <TH2D.h>
22 #include <TLatex.h>
23 #include <TPad.h>
24 #include <TPaveLabel.h>
25 #include <TChain.h>
26 #include <TGraph.h>
27 #include <TMultiGraph.h>
28 #include <TLine.h>
29 #include <TStyle.h>
30 #include <TF1.h>
31 #include <TH1F.h>
32 #include <TTree.h>
33 #include <TCanvas.h>
34 #include <PamLevel2.h>
35 #include <iomanip>
36 //
37 #include <ToFCore.h>
38 //
39 // Declaration of the core fortran routines
40 //
41 #define tofl2com tofl2com_
42 extern "C" int tofl2com();
43 #define toftrk toftrk_
44 extern "C" int toftrk();
45 ///extern "C" int toftrk(float);
46 #define rdtofcal rdtofcal_
47 //extern "C" int rdtofcal(char [], int *);
48 extern "C" int rdtofcal(const char *, int *);
49
50 using namespace std;
51
52 #endif
53
54 //===================================
55 // global variables
56 //===================================
57 TString DIR;
58 TString OUTF;
59 TString LIST;
60 TString OPTIONS;
61 TString CALIBF;
62
63 PamLevel2 *pam_event = NULL;
64
65 //==========================================
66 //000000000000000000000000000000000000000000
67 //==========================================
68 Int_t Loop(TString ddir,TString list, TString options){
69 //
70 //
71 extern struct ToFInput tofinput_;
72 extern struct ToFOutput tofoutput_;
73 // --------------------
74 // read input file/list
75 // --------------------
76 pam_event = new PamLevel2(ddir,list,options);
77 pam_event->SetSELLI(2);
78 TTree::SetMaxTreeSize(1000*Long64_t(2000000000));
79 //
80 TString outfile_t = "";
81 outfile_t = OUTF;
82 outfile_t.Append(".root");
83 TFile *outft = (TFile*)gROOT->FindObject(outfile_t);
84 if (outft) outft->Close();
85 outft = new TFile(outfile_t,"RECREATE");
86 if(outft->IsZombie()){
87 cout << "Output file could not be created\n";
88 return 1;
89 };
90 cout << "Created output file: "<<outft->GetName()<<endl;
91 //
92 outft->cd();
93 ToFLevel2 *tof = new ToFLevel2();
94 ToFdEdx *tofdedx = new ToFdEdx();
95 TTree *toft = new TTree("ToF","PAMELA Level2 ToF data");
96 toft->SetAutoSave(900000000000000LL);
97 tof->Set();//ELENA **TEMPORANEO?**
98 toft->Branch("ToFLevel2","ToFLevel2",&tof);
99 //
100 pam_event->CreateCloneTrees(outft);
101 //
102 Int_t ntrkentry = 0;
103 Int_t npmtentry = 0;
104 Float_t xleft=0;
105 Float_t xright=0;
106 Float_t yleft=0;
107 Float_t yright=0;
108 //
109 ULong64_t nevents = pam_event->GetEntries();
110 printf("\n\n Running on %llu events \n\n",nevents);
111 //
112 TFile *tfile = TFile::Open(list.Data(),"READ");
113 TTree *toftr = (TTree*)tfile->Get("ToF");
114 ToFLevel2 *tofl2 = new ToFLevel2();
115 toftr->SetBranchAddress("ToFLevel2", &tofl2);
116 ULong64_t ntofev = toftr->GetEntries();
117 printf(" ToF events are %llu \n",ntofev);
118 //
119 //
120 GL_PARAM *glparam = new GL_PARAM();
121 TrkLevel2 *trk = new TrkLevel2();
122 //
123 TString host;
124 TString user;
125 TString psw;
126 const char *pamdbhost=gSystem->Getenv("PAM_DBHOST");
127 const char *pamdbuser=gSystem->Getenv("PAM_DBUSER");
128 const char *pamdbpsw=gSystem->Getenv("PAM_DBPSW");
129 if ( !pamdbhost ) pamdbhost = "";
130 if ( !pamdbuser ) pamdbuser = "";
131 if ( !pamdbpsw ) pamdbpsw = "";
132 if ( strcmp(pamdbhost,"") ) host = pamdbhost;
133 if ( strcmp(pamdbuser,"") ) user = pamdbuser;
134 if ( strcmp(pamdbpsw,"") ) psw = pamdbpsw;
135 //
136 //
137 TSQLServer *dbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data());
138 if ( !dbc->IsConnected() ) return 1;
139 stringstream myquery;
140 myquery.str("");
141 myquery << "SET time_zone='+0:00'";
142 dbc->Query(myquery.str().c_str());
143 glparam->Query_GL_PARAM(1,1,dbc); // parameters stored in DB in GL_PRAM table
144 trk->LoadField(glparam->PATH+glparam->NAME);
145 //
146 Long64_t myrun = 0;
147 Long64_t mysrun = -1;
148 Int_t adc[4][12];
149 Int_t tdc[4][12];
150 Float_t tdcc[4][12];
151 //
152 Bool_t defcal = true;
153 for(ULong64_t iev=0; iev<nevents; iev++){
154
155
156 //==================================================================
157
158
159 // for(ULong64_t iev=0; iev<50; iev++){
160 //
161 if ( !(iev%1000) ) printf(" %iK \n",(int)iev/1000);
162 pam_event->Clear();
163 //
164 if( pam_event->GetEntry(iev) ){
165 //
166 pam_event->FillCloneTrees();
167 toftr->GetEntry(iev);
168 //
169 // variable to save information about the tof calibration used
170 //
171 //
172 myrun=pam_event->GetRunInfo()->ID;
173 if ( myrun != mysrun ){
174 // printf(" rhtime %u myrun %i mysrun %i \n",pam_event->GetRunInfo()->RUNHEADER_TIME,(int)myrun,(int)mysrun);
175 mysrun = myrun;
176 //
177
178
179
180
181 Int_t error=glparam->Query_GL_PARAM(pam_event->GetRunInfo()->RUNHEADER_TIME,204,dbc); // parameters stored in DB in GL_PRAM table
182 if ( error<0 ) {
183 return(1);
184 };
185 //
186 tofdedx->ReadParAtt((glparam->PATH+glparam->NAME).Data());
187 printf(" Reading Attenuation file: %s \n",(glparam->PATH+glparam->NAME).Data());
188
189 //
190 error=glparam->Query_GL_PARAM(pam_event->GetRunInfo()->RUNHEADER_TIME,205,dbc); // parameters stored in DB in GL_PRAM table
191 if ( error<0 ) {
192 return(1);
193 };
194 //
195 tofdedx->ReadParPos((glparam->PATH+glparam->NAME).Data());
196 printf(" Reading desaturation1 file: %s \n",(glparam->PATH+glparam->NAME).Data());
197
198 //
199 error=glparam->Query_GL_PARAM(pam_event->GetRunInfo()->RUNHEADER_TIME,206,dbc); // parameters stored in DB in GL_PRAM table
200 if ( error<0 ) {
201 return(1);
202 };
203 //
204 tofdedx->ReadParBBneg((glparam->PATH+glparam->NAME).Data());
205 printf(" Reading BBneg file: %s \n",(glparam->PATH+glparam->NAME).Data());
206
207 //
208 error=glparam->Query_GL_PARAM(pam_event->GetRunInfo()->RUNHEADER_TIME,207,dbc); // parameters stored in DB in GL_PRAM table
209 if ( error<0 ) {
210 return(1);
211 };
212 //
213 tofdedx->ReadParBBpos((glparam->PATH+glparam->NAME).Data());
214 printf(" Reading BBpos file: %s \n",(glparam->PATH+glparam->NAME).Data());
215
216 //
217 error=glparam->Query_GL_PARAM(pam_event->GetRunInfo()->RUNHEADER_TIME,208,dbc); // parameters stored in DB in GL_PRAM table
218 if ( error<0 ) {
219 return(1);
220 };
221 //
222 tofdedx->ReadParDesatBB((glparam->PATH+glparam->NAME).Data());
223 printf(" Reading desaturation2 file: %s \n",(glparam->PATH+glparam->NAME).Data());
224
225 tofdedx->CheckConnectors(pam_event->GetRunInfo()->RUNHEADER_TIME,glparam,dbc);
226
227
228 error=glparam->Query_GL_PARAM(pam_event->GetRunInfo()->RUNHEADER_TIME,201,dbc); // parameters stored in DB in GL_PRAM table
229 if ( error<0 ) {
230 return(1);
231 };
232 //
233 printf(" Reading ToF parameter file: %s \n",(glparam->PATH+glparam->NAME).Data());
234 //
235 if ( (UInt_t)glparam->TO_TIME != (UInt_t)4294967295UL ) defcal = false;
236 //
237 TString pippo=(glparam->PATH+glparam->NAME).Data();
238 //Int_t nlen = (Int_t)(glparam->PATH+glparam->NAME).Length();
239 //rdtofcal((char *)(glparam->PATH+glparam->NAME).Data(),&nlen);
240 Int_t nlen = (Int_t)pippo.Length();
241 rdtofcal(pippo.Data(),&nlen);
242 //
243
244 tofdedx->Clear();
245
246 };
247 //
248
249
250 //================================================================
251 //==================================================================
252 //--- Define absolute time
253 UInt_t tabs=pam_event->GetOrbitalInfo()->absTime;
254
255 if ( !(iev%100000) ) printf(" ATIME %u re %u \n",(Int_t)tabs,(UInt_t)iev);
256
257 //==================================================================
258 Float_t dedx_corr_m[2000][48];//,dedx_corr[48];
259 Double_t mtime[2000],t1,t2,tm;
260 // Float_t yhelp1,yhelp2,slope,inter,thelp1,thelp2;
261 Float_t yhelp1,yhelp2;
262 Float_t xmean1,xwidth1;
263
264 Int_t ical=0;
265 Int_t ii,j,jj;
266
267 if (iev==0) {
268
269 ical=0; // counter set to zero if first-time reading
270
271 //-----------------------------------------------------------
272 // Here I read the dEdx_korr parameters
273 //-----------------------------------------------------------
274
275 jj=0;
276 printf(" READING NEW CALIBRATION FILE: %s \n",CALIBF.Data());
277
278 ifstream fin(CALIBF.Data());
279 //cout << "topolino" << endl;
280 while (! fin.eof()) {
281 fin>>t1>>tm>>t2;
282 //cout<<jj<<" "<<tm<<endl;
283 //cout << setiosflags(ios::fixed) << setw(10) << setprecision(3) << tm << endl;
284 mtime[jj]=tm;
285 for (ii=0; ii<48;ii++) {
286 fin>>j>>xmean1>>xwidth1;
287 dedx_corr_m[jj][ii]=xmean1;
288 }
289 jj=jj+1;
290 // printf(" kk %i \n",jj);
291 }
292 // printf(" 1ical %i \n",ical);
293 //cout << "pippo" << endl;
294 fin.close();
295
296
297 while (tabs<mtime[ical] || tabs > mtime[ical+1]) {
298 // printf(" ical %i \n",ical);
299 ical = ical+1;
300 }
301 // ical = ical-1;
302 //cout<<"abs time "<<tabs<<" limit low "<<mtime[ical]<<" limit up "<<mtime[ical+1]<<" ical "<<ical<<endl;
303 } // if iev==0...
304
305 //==================================================================
306 //== End of first time reading & filling the array
307 //==================================================================
308
309 //==================================================================
310 //== if time is outside time limits:
311 //==================================================================
312
313 if (tabs<mtime[ical] || tabs>mtime[ical+1]) {
314 cout<<"Checking Time Limits!"<<endl;
315 ical=0;
316 while (tabs > mtime[ical+1] || tabs<mtime[ical]) {
317 ical = ical+1;
318 }
319 // ical = ical-1;
320 //cout<<"abs time "<<tabs<<" limit low "<<mtime[ical]<<" limit up "<<mtime[ical+1]<<" ical "<<ical<<endl;
321 };
322
323 // printf(" 2ical %i \n",ical);
324 //==================================================================
325 //== interpolate betwen time limits
326 //==================================================================
327
328 // thelp1 = mtime[ical];
329 // thelp2 = mtime[ical+1];
330
331 for (ii=0; ii<48;ii++) {
332 yhelp1 = fabs(dedx_corr_m[ical][ii]);
333 // yhelp1 = 6.;
334 if ( yhelp1 < 0.1 ) yhelp1 = 4.;
335 yhelp2 = fabs(dedx_corr_m[ical+1][ii]);
336 // yhelp2 = 6.;
337 if ( yhelp2 < 0.1 ) yhelp2 = 4.;
338 // slope = (yhelp2-yhelp1)/(thelp2-thelp1);
339 // inter = yhelp1 - slope*thelp1;
340 // dedx_corr[ii] = slope*tabs + inter;
341 // if (ii==0) cout<<thelp1<<" "<<thelp2<<" "<<tabs<<" "<<yhelp1<<" "<<yhelp2<<" "<<dedx_corr[0]<<endl;
342 }
343
344 //================================================================
345 //================================================================
346 // printf("cpippo \n");
347
348 // process tof data
349 //
350 for (Int_t hh=0; hh<12;hh++){
351 for (Int_t kk=0; kk<4;kk++){
352 adc[kk][hh] = 4095;
353 tdc[kk][hh] = 4095;
354 tdcc[kk][hh] = 4095.;
355 tofinput_.adc[hh][kk] = 4095;
356 tofinput_.tdc[hh][kk] = 4095;
357 };
358 };
359 // memset(adc, 0, 12*4*sizeof(Int_t));
360 // memset(tdc, 0, 12*4*sizeof(Int_t));
361 Int_t gg = 0;
362 Int_t hh = 0;
363 Int_t adcf[48];
364 memset(adcf, 0, 48*sizeof(Int_t));
365 Int_t tdcf[48];
366 memset(tdcf, 0, 48*sizeof(Int_t));
367 for (Int_t pm=0; pm < tofl2->ntrk() ; pm++){
368 ToFTrkVar *ttf = tofl2->GetToFTrkVar(pm);
369 for ( Int_t nc=0; nc < ttf->npmttdc; nc++){
370 if ( (ttf->tdcflag).At(nc) != 0 ) tdcf[(ttf->pmttdc).At(nc)] = 1;
371 };
372 for ( Int_t nc=0; nc < ttf->npmtadc; nc++){
373 if ( (ttf->adcflag).At(nc) != 0 ) adcf[(ttf->pmtadc).At(nc)] = 1;
374 };
375 };
376 //
377 for (Int_t pm=0; pm < tofl2->npmt() ; pm++){
378 ToFPMT *pmt = tofl2->GetToFPMT(pm);
379 tofl2->GetPMTIndex(pmt->pmt_id, gg, hh);
380 if ( adcf[pmt->pmt_id] == 0 ){
381 tofinput_.adc[gg][hh] = (int)pmt->adc;
382 adc[hh][gg] = (int)pmt->adc;
383 };
384 if ( tdcf[pmt->pmt_id] == 0 ){
385 tofinput_.tdc[gg][hh] = (int)pmt->tdc;
386 tdc[hh][gg] = (int)pmt->tdc;
387 };
388 tdcc[hh][gg] = (float)pmt->tdc_tw;
389 // Int_t pppid = tof->GetPMTid(hh,gg);
390 // printf(" pm %i pmt_id %i pppid %i hh %i gg %i tdcc %f tdc %f adc %f \n",pm,pmt->pmt_id,pppid,hh,gg,pmt->tdc_tw,pmt->tdc,pmt->adc);
391 };
392 for (Int_t hh=0; hh<12;hh++){
393 for (Int_t kk=0; kk<4;kk++){
394 tofdedx->Init(kk,hh,adc[kk][hh]);
395 };
396 };
397 // for (Int_t pm=0; pm <48 ; pm++){
398 // tof->GetPMTIndex(pm, gg, hh);
399 // tofinput_.tdc[hh][gg] = (int)500.;
400 // tofinput_.adc[hh][gg] = (int)500.;
401 // tdc[hh][gg] = (int)500.;
402 // adc[hh][gg] = (int)500.;
403 // // printf(" hh %i gg %i tdc %f adc %f \n",hh,gg,pmt->tdc,pmt->adc);
404 // };
405 //
406 for (Int_t hh=0; hh<5;hh++){
407 tofinput_.patterntrig[hh]=pam_event->GetTrigLevel2()->patterntrig[hh];
408 };
409 //
410 tof->Clear();
411 Int_t pmt_id = 0;
412 ToFPMT *t_pmt = new ToFPMT();
413 if(!(tof->PMT))tof->PMT = new TClonesArray("ToFPMT",12); //ELENA
414 TClonesArray &tpmt = *tof->PMT;
415 ToFTrkVar *t_tof = new ToFTrkVar();
416 if(!(tof->ToFTrk))tof->ToFTrk = new TClonesArray("ToFTrkVar",2); //ELENA
417 TClonesArray &t = *tof->ToFTrk;
418 //
419 //
420 // Here we have calibrated data, ready to be passed to the FORTRAN routine which will extract common and track-related variables.
421 //
422 npmtentry = 0;
423 //
424 ntrkentry = 0;
425 //
426 // Calculate tracks informations from ToF alone
427 //
428 tofl2com();
429 //
430 memcpy(tof->tof_j_flag,tofoutput_.tof_j_flag,6*sizeof(Int_t));
431 //
432 t_tof->trkseqno = -1;
433 //
434 // and now we must copy from the output structure to the level2 class:
435 //
436 t_tof->npmttdc = 0;
437 //
438 for (Int_t hh=0; hh<12;hh++){
439 for (Int_t kk=0; kk<4;kk++){
440 if ( tofoutput_.tofmask[hh][kk] != 0 ){
441 pmt_id = tof->GetPMTid(kk,hh);
442 t_tof->pmttdc.AddAt(pmt_id,t_tof->npmttdc);
443 t_tof->tdcflag.AddAt(tofoutput_.tdcflagtof[hh][kk],t_tof->npmttdc); // gf: Jan 09/07
444 t_tof->npmttdc++;
445 };
446 };
447 };
448 for (Int_t kk=0; kk<13;kk++){
449 t_tof->beta[kk] = tofoutput_.betatof_a[kk];
450 }
451 //
452 memcpy(t_tof->xtofpos,tofoutput_.xtofpos,sizeof(t_tof->xtofpos));
453 memcpy(t_tof->ytofpos,tofoutput_.ytofpos,sizeof(t_tof->ytofpos));
454 memcpy(t_tof->xtr_tof,tofoutput_.xtr_tof,sizeof(t_tof->xtr_tof));
455 memcpy(t_tof->ytr_tof,tofoutput_.ytr_tof,sizeof(t_tof->ytr_tof));
456 // {
457 // Float_t xtof_temp[6]={0.,t_tof->xtofpos[0],t_tof->xtofpos[1],0.,0.,t_tof->xtofpos[2]};
458 // Float_t ytof_temp[6]={t_tof->ytofpos[0],0.,0.,t_tof->ytofpos[1],t_tof->ytofpos[2],0.};
459 // tofdedx->Process(pam_event->GetOrbitalInfo()->absTime,t_tof->beta[12], (Float_t *)xtof_temp,(Float_t *)ytof_temp);
460 // }
461 // t_tof->npmtadc = 0;
462
463
464 // for (Int_t hh=0; hh<12;hh++){
465 // for (Int_t kk=0; kk<4;kk++){
466 // pmt_id = tof->GetPMTid(kk,hh);
467 // if ( tofdedx->GetdEdx_pmt(pmt_id)>-1. ){
468 // t_tof->dedx.AddAt((tofdedx->GetdEdx_pmt(pmt_id)*36./pow(dedx_corr[pmt_id],2)),t_tof->npmtadc);
469 // t_tof->pmtadc.AddAt(pmt_id,t_tof->npmtadc);
470 // t_tof->adcflag.AddAt(0,t_tof->npmtadc); // gf: Jan 09/07
471 // t_tof->npmtadc++;
472 // };
473 // };
474 // };
475 {
476
477 Float_t xtof_temp[6]={100.,100.,100.,100.,100.,100.};
478 Float_t ytof_temp[6]={100.,100.,100.,100.,100.,100.};
479
480 if(t_tof->xtofpos[0]<100. && t_tof->ytofpos[0]<100.){
481 xtof_temp[1]=t_tof->xtofpos[0];
482 ytof_temp[0]=t_tof->ytofpos[0];
483 }else if(t_tof->xtofpos[0]>=100. && t_tof->ytofpos[0]<100.){
484 ytof_temp[0]=t_tof->ytofpos[0];
485 tof->GetPaddleGeometry(0,(Int_t)log2(tof->tof_j_flag[0]),xleft, xright, yleft, yright);
486 xtof_temp[1]=xleft+2.55;
487 }else if(t_tof->ytofpos[0]>=100. && t_tof->xtofpos[0]<100.){
488 xtof_temp[1]=t_tof->xtofpos[0];
489 tof->GetPaddleGeometry(1,(Int_t)log2(tof->tof_j_flag[1]),xleft, xright, yleft, yright);
490 ytof_temp[0]=yleft+2.75;
491 }
492
493 if(t_tof->xtofpos[1]<100. && t_tof->ytofpos[1]<100.){
494 xtof_temp[2]=t_tof->xtofpos[1];
495 ytof_temp[3]=t_tof->ytofpos[1];
496 }else if(t_tof->xtofpos[1]>=100. && t_tof->ytofpos[1]<100.){
497 ytof_temp[3]=t_tof->ytofpos[1];
498 tof->GetPaddleGeometry(3,(Int_t)log2(tof->tof_j_flag[3]),xleft, xright, yleft, yright);
499 xtof_temp[2]=xleft+4.5;
500 }else if(t_tof->ytofpos[1]>=100. && t_tof->xtofpos[1]<100.){
501 xtof_temp[2]=t_tof->xtofpos[1];
502 tof->GetPaddleGeometry(2,(Int_t)log2(tof->tof_j_flag[2]),xleft, xright, yleft, yright);
503 ytof_temp[3]=yleft+3.75;
504 }
505
506 if(t_tof->xtofpos[2]<100. && t_tof->ytofpos[2]<100.){
507 xtof_temp[5]=t_tof->xtofpos[2];
508 ytof_temp[4]=t_tof->ytofpos[2];
509 }else if(t_tof->xtofpos[2]>=100. && t_tof->ytofpos[2]<100.){
510 ytof_temp[4]=t_tof->ytofpos[2];
511 tof->GetPaddleGeometry(4,(Int_t)log2(tof->tof_j_flag[4]),xleft, xright, yleft, yright);
512 xtof_temp[5]=xleft+3;
513 }else if(t_tof->ytofpos[2]>=100. && t_tof->xtofpos[2]<100.){
514 xtof_temp[5]=t_tof->xtofpos[2];
515 tof->GetPaddleGeometry(5,(Int_t)log2(tof->tof_j_flag[5]),xleft, xright, yleft, yright);
516 ytof_temp[4]=yleft+2.5;
517 }
518 // Float_t xtof_temp[6]={0.,t_tof->xtofpos[0],t_tof->xtofpos[1],0.,0.,t_tof->xtofpos[2]};
519 // Float_t ytof_temp[6]={t_tof->ytofpos[0],0.,0.,t_tof->ytofpos[1],t_tof->ytofpos[2],0.};
520 // tofdedx->Process(atime,t_tof->beta[12], (Float_t *)xtof_temp,(Float_t *)ytof_temp);
521 tofdedx->Process(pam_event->GetOrbitalInfo()->absTime,t_tof->beta[12], (Float_t *)xtof_temp,(Float_t *)ytof_temp);
522 t_tof->npmtadc = 0;
523 for (Int_t hh=0; hh<12;hh++){
524 for (Int_t kk=0; kk<4;kk++){
525 pmt_id = tof->GetPMTid(kk,hh);
526 Int_t Iplane=-1;
527 Int_t Ipaddle=-1;
528 // Int_t IpaddleT=-1;
529 tof->GetPMTPaddle(pmt_id, Iplane, Ipaddle);
530 tof->GetPaddleGeometry(Iplane,Ipaddle,xleft,xright,yleft,yright);
531 if ( tofdedx->GetdEdx_pmt(pmt_id)>-1. &&((xtof_temp[Iplane]>=xleft&&xtof_temp[Iplane]<=xright) || (ytof_temp[Iplane]>=yleft&&ytof_temp[Iplane]<=yright)) ){ //attenzione:qui va inserito un controllo sulla traccia tof o sulle variabili di posizione !!!!
532 //t_tof->dedx.AddAt((tofdedx->GetdEdx_pmt(pmt_id)*4./dedx_corr[pmt_id]),t_tof->npmtadc);
533 //t_tof->dedx.AddAt((tofdedx->GetdEdx_pmt(pmt_id)*dedx_corr[pmt_id]/4.),t_tof->npmtadc);//annullo wolfrizzazione
534 t_tof->dedx.AddAt((tofdedx->GetdEdx_pmt(pmt_id)),t_tof->npmtadc);//annullo wolfrizzazione
535 t_tof->pmtadc.AddAt(pmt_id,t_tof->npmtadc);
536 t_tof->adcflag.AddAt(0,t_tof->npmtadc); // gf: Jan 09/07
537 t_tof->npmtadc++;
538 };
539 };
540 };
541 };
542
543
544
545
546
547
548
549
550
551
552 // for (Int_t hh=0; hh<12;hh++){
553 // for (Int_t kk=0; kk<4;kk++){
554 // if ( tofoutput_.adctof_c[hh][kk] < 1000 ){
555 // // t_tof->dedx.AddAt(tofoutput_.adctof_c[hh][kk],t_tof->npmtadc); // EMILIANO
556 // pmt_id = tof->GetPMTid(kk,hh);
557 // t_tof->dedx.AddAt((tofoutput_.adctof_c[hh][kk]*4./dedx_corr[pmt_id]),t_tof->npmtadc); // EMILIANO
558 // t_tof->pmtadc.AddAt(pmt_id,t_tof->npmtadc);
559 // t_tof->adcflag.AddAt(tofoutput_.adcflagtof[hh][kk],t_tof->npmtadc); // gf: Jan 09/07
560 // t_tof->npmtadc++;
561 // };
562 // };
563 // };
564 //
565
566 //
567 new(t[ntrkentry]) ToFTrkVar(*t_tof);
568 ntrkentry++;
569 t_tof->Clear();
570 //
571 //
572 //
573 t_pmt->Clear();
574 //
575 for (Int_t hh=0; hh<12;hh++){
576 for (Int_t kk=0; kk<4;kk++){
577 // new WM
578 // if ( tofoutput_.tdc_c[hh][kk] < 4095 || adc[kk][hh] < 4095 || tdc[kk][hh] < 4095 ){
579 if ( tdcc[kk][hh] < 4095. || adc[kk][hh] < 4095 || tdc[kk][hh] < 4095 ){
580 //
581 t_pmt->pmt_id = tof->GetPMTid(kk,hh);
582 t_pmt->tdc_tw = tofoutput_.tdc_c[hh][kk];
583 t_pmt->adc = (Float_t)adc[kk][hh];
584 t_pmt->tdc = (Float_t)tdc[kk][hh];
585 //
586 new(tpmt[npmtentry]) ToFPMT(*t_pmt);
587 npmtentry++;
588 t_pmt->Clear();
589 };
590 };
591 };
592 //
593 // Calculate track-related variables
594 //
595 if ( pam_event->GetTrkLevel2()->ntrk() > 0 ){
596 //
597 // We have at least one track
598 //
599 //
600 // Run over tracks
601 //
602 for(Int_t nt=0; nt < pam_event->GetTrkLevel2()->ntrk(); nt++){
603 //
604 TrkTrack *ptt = pam_event->GetTrkLevel2()->GetStoredTrack(nt);
605 //
606 // Copy the alpha vector in the input structure
607 //
608 for (Int_t e = 0; e < 5 ; e++){
609 tofinput_.al_pp[e] = ptt->al[e];
610 };
611
612 // new input for 9th reduction: tracker dEdx
613 tofinput_.trkmip = ptt->GetDEDX();
614
615 //
616 // Get tracker related variables for this track
617 //
618 toftrk();
619 // toftrk(thelp);
620 //
621 // Copy values in the class from the structure (we need to use a temporary class to store variables).
622 //
623 t_tof->npmttdc = 0;
624 for (Int_t hh=0; hh<12;hh++){
625 for (Int_t kk=0; kk<4;kk++){
626 if ( tofoutput_.tofmask[hh][kk] != 0 ){
627 pmt_id = tof->GetPMTid(kk,hh);
628 t_tof->pmttdc.AddAt(pmt_id,t_tof->npmttdc);
629 t_tof->tdcflag.AddAt(tofoutput_.tdcflag[hh][kk],t_tof->npmttdc); // gf: Jan 09/07
630 t_tof->npmttdc++;
631 };
632 };
633 };
634 for (Int_t kk=0; kk<13;kk++){
635 t_tof->beta[kk] = tofoutput_.beta_a[kk];
636 };
637 //
638 memcpy(t_tof->xtofpos,tofoutput_.xtofpos,sizeof(t_tof->xtofpos));
639 memcpy(t_tof->ytofpos,tofoutput_.ytofpos,sizeof(t_tof->ytofpos));
640 memcpy(t_tof->xtr_tof,tofoutput_.xtr_tof,sizeof(t_tof->xtr_tof));
641 memcpy(t_tof->ytr_tof,tofoutput_.ytr_tof,sizeof(t_tof->ytr_tof));
642 //
643 tofdedx->Process(pam_event->GetOrbitalInfo()->absTime,t_tof->beta[12], (Float_t *)t_tof->xtr_tof,(Float_t *)t_tof->ytr_tof);
644 t_tof->npmtadc = 0;
645
646 for (Int_t hh=0; hh<12;hh++){
647 for (Int_t kk=0; kk<4;kk++){
648 pmt_id = tof->GetPMTid(kk,hh);
649 Int_t Iplane=-1;
650 Int_t Ipaddle=-1;
651 Int_t IpaddleT=-1;
652 tof->GetPMTPaddle(pmt_id, Iplane, Ipaddle);
653 IpaddleT=tof->GetPaddleIdOfTrack(t_tof->xtr_tof[Iplane],t_tof->ytr_tof[Iplane], Iplane,0.0);
654 if ( tofdedx->GetdEdx_pmt(pmt_id) > -1. && Ipaddle==IpaddleT ){
655 //t_tof->dedx.AddAt((tofdedx->GetdEdx_pmt(pmt_id)*4./dedx_corr[pmt_id]),t_tof->npmtadc);
656 //t_tof->dedx.AddAt((tofdedx->GetdEdx_pmt(pmt_id)*dedx_corr[pmt_id]/4.),t_tof->npmtadc);//annullo wolfrizzazione
657 t_tof->dedx.AddAt((tofdedx->GetdEdx_pmt(pmt_id)),t_tof->npmtadc);//annullo wolfrizzazione
658 t_tof->pmtadc.AddAt(pmt_id,t_tof->npmtadc);
659 t_tof->adcflag.AddAt(0,t_tof->npmtadc); // gf: Jan 09/07
660 t_tof->npmtadc++;
661 };
662
663 };
664 };
665
666
667 // for (Int_t hh=0; hh<12;hh++){
668 // for (Int_t kk=0; kk<4;kk++){
669 // pmt_id = tof->GetPMTid(kk,hh);
670 // if ( tofdedx->GetdEdx_pmt(pmt_id) > -1. ){
671 // t_tof->dedx.AddAt((tofdedx->GetdEdx_pmt(pmt_id)*36./pow(dedx_corr[pmt_id],2)),t_tof->npmtadc);
672 // t_tof->pmtadc.AddAt(pmt_id,t_tof->npmtadc);
673 // t_tof->adcflag.AddAt(0,t_tof->npmtadc); // gf: Jan 09/07
674 // printf(" nt %i npmtadc %i dedx %f dedx corr %f\n",nt,t_tof->npmtadc,(tofdedx->GetdEdx_pmt(pmt_id)*36./pow(dedx_corr[pmt_id],2)),dedx_corr[pmt_id]);
675 // t_tof->npmtadc++;
676 // };
677
678 // };
679 // };
680 // t_tof->npmtadc = 0;
681 // for (Int_t hh=0; hh<12;hh++){
682 // for (Int_t kk=0; kk<4;kk++){
683 // if ( tofoutput_.adc_c[hh][kk] < 1000 ){
684 // // t_tof->dedx.AddAt(tofoutput_.adc_c[hh][kk],t_tof->npmtadc); // EMILIANO
685 // pmt_id = tof->GetPMTid(kk,hh);
686 // t_tof->dedx.AddAt((tofoutput_.adc_c[hh][kk]*4./dedx_corr[pmt_id]),t_tof->npmtadc); // EMILIANO
687 // t_tof->pmtadc.AddAt(pmt_id,t_tof->npmtadc);
688 // t_tof->adcflag.AddAt(tofoutput_.adcflag[hh][kk],t_tof->npmtadc); // gf: Jan 09/07
689 // t_tof->npmtadc++;
690 // };
691 // };
692 // };
693 //
694
695 //
696 // Store the tracker track number in order to be sure to have shyncronized data during analysis
697 //
698 t_tof->trkseqno = nt;
699 //
700 // create a new object for this event with track-related variables
701 //
702 new(t[ntrkentry]) ToFTrkVar(*t_tof);
703 ntrkentry++;
704 t_tof->Clear();
705 //
706 }; // loop on all the tracks
707 };
708 //
709 tof->unpackError = tofl2->unpackError;
710 if ( defcal ){
711 tof->default_calib = 1;
712 } else {
713 tof->default_calib = 0;
714 };
715 //
716 // Fill the rootple
717 //
718 toft->Fill();
719 //
720 // toft->Show();
721 // toftr->Show();
722 //
723 delete t_tof;
724 //
725 //
726 };
727 //
728 };
729 //
730 pam_event->Clear();
731 //
732 outft->cd();
733 pam_event->WriteCloneTrees();
734 toftr->Delete();
735
736 toft->Write();
737 outft->Close();
738 //
739 return 0;
740 //
741 }
742
743 /////////////////////////////////////////////////////////////////
744
745 #if !defined(__CINT__)
746
747 void usage(){
748
749 cout << "------------------------------------------------------------"<<endl;
750 cout << "Loop over events (on one or more Level2-files), applying some selection cuts (defined in My-Selection.cpp), \n";
751 cout << "creating output histograms (defined in My-Histos.cpp) and/or trees with selected events. \n \n ";
752 cout << "USAGE:"<<endl;
753 cout << "-processDir DIR - Level2 data directory \n";
754 cout << "-processList LIST - list of files (.txt) or single file (.root) to be analysed \n";
755 cout << "-outputFile PATH - name of the output file \n";
756 cout << "-calibFile PATH+NAME - name of the calibration file \n";
757 cout << "-NumEvents XXX - number of events to be analysed \n";
758 cout << "--debug, -g - debug mode \n";
759 cout << "--help, -h - print this help \n";
760 cout << "-options [ options ] - options: \n";
761 cout << " fillHistos --> create an output file with histograms \n";
762 cout << " fillTree --> create an output file with trees storing the selected events \n ";
763 cout << " +(-)ALL --> inlcude(exclude) all trees and branches \n " ;
764 cout << " +(-)TRK1 +(-)TRK2 +(-)CAL1 +(-)CAL2 +(-)TOF +(-)TRG +(-)ND +(-)S4 +(-)ORB --> inlcude(exclude) trees and branches \n" ;
765 cout << "------------------------------------------------------------"<<endl;
766 }
767 //
768 int HandleInputPar(int argc, char **argv){
769
770 if(argc>1){
771
772 if(!strcmp(argv[1], "-h") || !strcmp(argv[1], "--help") ){
773 usage();
774 return(1);
775 };
776 // -----------------------
777 // Read input parameters
778 // -----------------------
779 DIR = gSystem->WorkingDirectory();
780 LIST = "";
781 OUTF = "myfile";
782 OPTIONS = "+AUTO -TOF";
783
784 for (int i = 1; i < argc; i++){
785 // -----------------------------------------------------//
786 if (!strcmp(argv[i], "-processDir")){
787 if (++i >= argc) throw -1;
788 DIR = argv[i];
789 cout << "processDir "<<DIR<<endl;
790 continue;
791 }
792 // -----------------------------------------------------//
793 else if (!strcmp(argv[i], "-processList")){
794 if (++i >= argc) throw -1;
795 LIST = argv[i];
796 cout << "processList "<<LIST<<endl;
797 continue;
798 }
799 // -----------------------------------------------------//
800 else if (!strcmp(argv[i], "-outputFile")){
801 if (++i >= argc) throw -1;
802 OUTF = argv[i];
803 cout << "outputFile "<<OUTF<<endl;
804 continue;
805 }
806 // -----------------------------------------------------//
807 // -----------------------------------------------------//
808 else if (!strcmp(argv[i], "-calibFile")){
809 if (++i >= argc) throw -1;
810 CALIBF = gSystem->ExpandPathName(argv[i]);
811 cout << "calibFile "<<CALIBF<<endl;
812 continue;
813 }
814 // -----------------------------------------------------//
815 else if (!strcmp(argv[i], "-options")){
816 if (++i >= argc) throw -1;
817 OPTIONS = argv[i];
818 if( OPTIONS.Contains("[") ){
819 do{
820 if (++i >= argc) throw -1;
821 OPTIONS.Append(argv[i]);
822 }while(!OPTIONS.Contains("]"));
823 }else cout << "wrong option format --> ignoring " << endl;
824 }
825 else{
826 cout << "Unidentified input parameter. Ignored."<< endl;
827 };
828 };
829 }else{
830 usage();
831 return(1);
832 };
833 // -----------------------
834 // Check input parameters
835 // -----------------------
836
837
838 return(0);
839
840 };
841 //
842
843 int main(int argc, char **argv)
844 {
845
846 if( HandleInputPar(argc,argv) )return(1);
847
848 cout << "OPTIONS "<<OPTIONS<<endl;
849 Loop(DIR,LIST,OPTIONS);
850
851 cout << "Back to main - end"<<endl;
852
853 return 0;
854
855 };
856
857 #endif

  ViewVC Help
Powered by ViewVC 1.1.23