/[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.6 - (show annotations) (download)
Thu Apr 12 12:04:20 2012 UTC (12 years, 10 months ago) by mocchiut
Branch: MAIN
Changes since 1.5: +32 -18 lines
New ToFReprocessing

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 xmean1,xwidth1;
262
263 Int_t ical,ii,j,jj;
264
265 if (iev==0) {
266
267 ical=0; // counter set to zero if first-time reading
268
269 //-----------------------------------------------------------
270 // Here I read the dEdx_korr parameters
271 //-----------------------------------------------------------
272
273 jj=0;
274 printf(" READING NEW CALIBRATION FILE: %s \n",CALIBF.Data());
275
276 ifstream fin(CALIBF.Data());
277 //cout << "topolino" << endl;
278 while (! fin.eof()) {
279 fin>>t1>>tm>>t2;
280 //cout<<jj<<" "<<tm<<endl;
281 //cout << setiosflags(ios::fixed) << setw(10) << setprecision(3) << tm << endl;
282 mtime[jj]=tm;
283 for (ii=0; ii<48;ii++) {
284 fin>>j>>xmean1>>xwidth1;
285 dedx_corr_m[jj][ii]=xmean1;
286 }
287 jj=jj+1;
288 // printf(" kk %i \n",jj);
289 }
290 // printf(" 1ical %i \n",ical);
291 //cout << "pippo" << endl;
292 fin.close();
293
294
295 while (tabs<mtime[ical] || tabs > mtime[ical+1]) {
296 // printf(" ical %i \n",ical);
297 ical = ical+1;
298 }
299 // ical = ical-1;
300 //cout<<"abs time "<<tabs<<" limit low "<<mtime[ical]<<" limit up "<<mtime[ical+1]<<" ical "<<ical<<endl;
301 } // if iev==0...
302
303 //==================================================================
304 //== End of first time reading & filling the array
305 //==================================================================
306
307 //==================================================================
308 //== if time is outside time limits:
309 //==================================================================
310
311 if (tabs<mtime[ical] || tabs>mtime[ical+1]) {
312 cout<<"Checking Time Limits!"<<endl;
313 ical=0;
314 while (tabs > mtime[ical+1] || tabs<mtime[ical]) {
315 ical = ical+1;
316 }
317 // ical = ical-1;
318 //cout<<"abs time "<<tabs<<" limit low "<<mtime[ical]<<" limit up "<<mtime[ical+1]<<" ical "<<ical<<endl;
319 };
320
321 // printf(" 2ical %i \n",ical);
322 //==================================================================
323 //== interpolate betwen time limits
324 //==================================================================
325
326 thelp1 = mtime[ical];
327 thelp2 = mtime[ical+1];
328
329 for (ii=0; ii<48;ii++) {
330 yhelp1 = fabs(dedx_corr_m[ical][ii]);
331 // yhelp1 = 6.;
332 if ( yhelp1 < 0.1 ) yhelp1 = 4.;
333 yhelp2 = fabs(dedx_corr_m[ical+1][ii]);
334 // yhelp2 = 6.;
335 if ( yhelp2 < 0.1 ) yhelp2 = 4.;
336 slope = (yhelp2-yhelp1)/(thelp2-thelp1);
337 inter = yhelp1 - slope*thelp1;
338 dedx_corr[ii] = slope*tabs + inter;
339 // if (ii==0) cout<<thelp1<<" "<<thelp2<<" "<<tabs<<" "<<yhelp1<<" "<<yhelp2<<" "<<dedx_corr[0]<<endl;
340 }
341
342 //================================================================
343 //================================================================
344 // printf("cpippo \n");
345
346 // process tof data
347 //
348 for (Int_t hh=0; hh<12;hh++){
349 for (Int_t kk=0; kk<4;kk++){
350 adc[kk][hh] = 4095;
351 tdc[kk][hh] = 4095;
352 tdcc[kk][hh] = 4095.;
353 tofinput_.adc[hh][kk] = 4095;
354 tofinput_.tdc[hh][kk] = 4095;
355 };
356 };
357 // memset(adc, 0, 12*4*sizeof(Int_t));
358 // memset(tdc, 0, 12*4*sizeof(Int_t));
359 Int_t gg = 0;
360 Int_t hh = 0;
361 Int_t adcf[48];
362 memset(adcf, 0, 48*sizeof(Int_t));
363 Int_t tdcf[48];
364 memset(tdcf, 0, 48*sizeof(Int_t));
365 for (Int_t pm=0; pm < tofl2->ntrk() ; pm++){
366 ToFTrkVar *ttf = tofl2->GetToFTrkVar(pm);
367 for ( Int_t nc=0; nc < ttf->npmttdc; nc++){
368 if ( (ttf->tdcflag).At(nc) != 0 ) tdcf[(ttf->pmttdc).At(nc)] = 1;
369 };
370 for ( Int_t nc=0; nc < ttf->npmtadc; nc++){
371 if ( (ttf->adcflag).At(nc) != 0 ) adcf[(ttf->pmtadc).At(nc)] = 1;
372 };
373 };
374 //
375 for (Int_t pm=0; pm < tofl2->npmt() ; pm++){
376 ToFPMT *pmt = tofl2->GetToFPMT(pm);
377 tofl2->GetPMTIndex(pmt->pmt_id, gg, hh);
378 if ( adcf[pmt->pmt_id] == 0 ){
379 tofinput_.adc[gg][hh] = (int)pmt->adc;
380 adc[hh][gg] = (int)pmt->adc;
381 };
382 if ( tdcf[pmt->pmt_id] == 0 ){
383 tofinput_.tdc[gg][hh] = (int)pmt->tdc;
384 tdc[hh][gg] = (int)pmt->tdc;
385 };
386 tdcc[hh][gg] = (float)pmt->tdc_tw;
387 // Int_t pppid = tof->GetPMTid(hh,gg);
388 // 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);
389 };
390 for (Int_t hh=0; hh<12;hh++){
391 for (Int_t kk=0; kk<4;kk++){
392 tofdedx->Init(kk,hh,adc[kk][hh]);
393 };
394 };
395 // for (Int_t pm=0; pm <48 ; pm++){
396 // tof->GetPMTIndex(pm, gg, hh);
397 // tofinput_.tdc[hh][gg] = (int)500.;
398 // tofinput_.adc[hh][gg] = (int)500.;
399 // tdc[hh][gg] = (int)500.;
400 // adc[hh][gg] = (int)500.;
401 // // printf(" hh %i gg %i tdc %f adc %f \n",hh,gg,pmt->tdc,pmt->adc);
402 // };
403 //
404 for (Int_t hh=0; hh<5;hh++){
405 tofinput_.patterntrig[hh]=pam_event->GetTrigLevel2()->patterntrig[hh];
406 };
407 //
408 tof->Clear();
409 Int_t pmt_id = 0;
410 ToFPMT *t_pmt = new ToFPMT();
411 if(!(tof->PMT))tof->PMT = new TClonesArray("ToFPMT",12); //ELENA
412 TClonesArray &tpmt = *tof->PMT;
413 ToFTrkVar *t_tof = new ToFTrkVar();
414 if(!(tof->ToFTrk))tof->ToFTrk = new TClonesArray("ToFTrkVar",2); //ELENA
415 TClonesArray &t = *tof->ToFTrk;
416 //
417 //
418 // Here we have calibrated data, ready to be passed to the FORTRAN routine which will extract common and track-related variables.
419 //
420 npmtentry = 0;
421 //
422 ntrkentry = 0;
423 //
424 // Calculate tracks informations from ToF alone
425 //
426 tofl2com();
427 //
428 memcpy(tof->tof_j_flag,tofoutput_.tof_j_flag,6*sizeof(Int_t));
429 //
430 t_tof->trkseqno = -1;
431 //
432 // and now we must copy from the output structure to the level2 class:
433 //
434 t_tof->npmttdc = 0;
435 //
436 for (Int_t hh=0; hh<12;hh++){
437 for (Int_t kk=0; kk<4;kk++){
438 if ( tofoutput_.tofmask[hh][kk] != 0 ){
439 pmt_id = tof->GetPMTid(kk,hh);
440 t_tof->pmttdc.AddAt(pmt_id,t_tof->npmttdc);
441 t_tof->tdcflag.AddAt(tofoutput_.tdcflagtof[hh][kk],t_tof->npmttdc); // gf: Jan 09/07
442 t_tof->npmttdc++;
443 };
444 };
445 };
446 for (Int_t kk=0; kk<13;kk++){
447 t_tof->beta[kk] = tofoutput_.betatof_a[kk];
448 }
449 //
450 memcpy(t_tof->xtofpos,tofoutput_.xtofpos,sizeof(t_tof->xtofpos));
451 memcpy(t_tof->ytofpos,tofoutput_.ytofpos,sizeof(t_tof->ytofpos));
452 memcpy(t_tof->xtr_tof,tofoutput_.xtr_tof,sizeof(t_tof->xtr_tof));
453 memcpy(t_tof->ytr_tof,tofoutput_.ytr_tof,sizeof(t_tof->ytr_tof));
454 // {
455 // Float_t xtof_temp[6]={0.,t_tof->xtofpos[0],t_tof->xtofpos[1],0.,0.,t_tof->xtofpos[2]};
456 // Float_t ytof_temp[6]={t_tof->ytofpos[0],0.,0.,t_tof->ytofpos[1],t_tof->ytofpos[2],0.};
457 // tofdedx->Process(pam_event->GetOrbitalInfo()->absTime,t_tof->beta[12], (Float_t *)xtof_temp,(Float_t *)ytof_temp);
458 // }
459 // t_tof->npmtadc = 0;
460
461
462 // for (Int_t hh=0; hh<12;hh++){
463 // for (Int_t kk=0; kk<4;kk++){
464 // pmt_id = tof->GetPMTid(kk,hh);
465 // if ( tofdedx->GetdEdx_pmt(pmt_id)>-1. ){
466 // t_tof->dedx.AddAt((tofdedx->GetdEdx_pmt(pmt_id)*36./pow(dedx_corr[pmt_id],2)),t_tof->npmtadc);
467 // t_tof->pmtadc.AddAt(pmt_id,t_tof->npmtadc);
468 // t_tof->adcflag.AddAt(0,t_tof->npmtadc); // gf: Jan 09/07
469 // t_tof->npmtadc++;
470 // };
471 // };
472 // };
473 {
474
475 Float_t xtof_temp[6]={100.,100.,100.,100.,100.,100.};
476 Float_t ytof_temp[6]={100.,100.,100.,100.,100.,100.};
477
478 if(t_tof->xtofpos[0]<100. && t_tof->ytofpos[0]<100.){
479 xtof_temp[1]=t_tof->xtofpos[0];
480 ytof_temp[0]=t_tof->ytofpos[0];
481 }else if(t_tof->xtofpos[0]>=100. && t_tof->ytofpos[0]<100.){
482 ytof_temp[0]=t_tof->ytofpos[0];
483 tof->GetPaddleGeometry(0,(Int_t)log2(tof->tof_j_flag[0]),xleft, xright, yleft, yright);
484 xtof_temp[1]=xleft+2.55;
485 }else if(t_tof->ytofpos[0]>=100. && t_tof->xtofpos[0]<100.){
486 xtof_temp[1]=t_tof->xtofpos[0];
487 tof->GetPaddleGeometry(1,(Int_t)log2(tof->tof_j_flag[1]),xleft, xright, yleft, yright);
488 ytof_temp[0]=yleft+2.75;
489 }
490
491 if(t_tof->xtofpos[1]<100. && t_tof->ytofpos[1]<100.){
492 xtof_temp[2]=t_tof->xtofpos[1];
493 ytof_temp[3]=t_tof->ytofpos[1];
494 }else if(t_tof->xtofpos[1]>=100. && t_tof->ytofpos[1]<100.){
495 ytof_temp[3]=t_tof->ytofpos[1];
496 tof->GetPaddleGeometry(3,(Int_t)log2(tof->tof_j_flag[3]),xleft, xright, yleft, yright);
497 xtof_temp[2]=xleft+4.5;
498 }else if(t_tof->ytofpos[1]>=100. && t_tof->xtofpos[1]<100.){
499 xtof_temp[2]=t_tof->xtofpos[1];
500 tof->GetPaddleGeometry(2,(Int_t)log2(tof->tof_j_flag[2]),xleft, xright, yleft, yright);
501 ytof_temp[3]=yleft+3.75;
502 }
503
504 if(t_tof->xtofpos[2]<100. && t_tof->ytofpos[2]<100.){
505 xtof_temp[5]=t_tof->xtofpos[2];
506 ytof_temp[4]=t_tof->ytofpos[2];
507 }else if(t_tof->xtofpos[2]>=100. && t_tof->ytofpos[2]<100.){
508 ytof_temp[4]=t_tof->ytofpos[2];
509 tof->GetPaddleGeometry(4,(Int_t)log2(tof->tof_j_flag[4]),xleft, xright, yleft, yright);
510 xtof_temp[5]=xleft+3;
511 }else if(t_tof->ytofpos[2]>=100. && t_tof->xtofpos[2]<100.){
512 xtof_temp[5]=t_tof->xtofpos[2];
513 tof->GetPaddleGeometry(5,(Int_t)log2(tof->tof_j_flag[5]),xleft, xright, yleft, yright);
514 ytof_temp[4]=yleft+2.5;
515 }
516 // Float_t xtof_temp[6]={0.,t_tof->xtofpos[0],t_tof->xtofpos[1],0.,0.,t_tof->xtofpos[2]};
517 // Float_t ytof_temp[6]={t_tof->ytofpos[0],0.,0.,t_tof->ytofpos[1],t_tof->ytofpos[2],0.};
518 // tofdedx->Process(atime,t_tof->beta[12], (Float_t *)xtof_temp,(Float_t *)ytof_temp);
519 tofdedx->Process(pam_event->GetOrbitalInfo()->absTime,t_tof->beta[12], (Float_t *)xtof_temp,(Float_t *)ytof_temp);
520 t_tof->npmtadc = 0;
521 for (Int_t hh=0; hh<12;hh++){
522 for (Int_t kk=0; kk<4;kk++){
523 pmt_id = tof->GetPMTid(kk,hh);
524 Int_t Iplane=-1;
525 Int_t Ipaddle=-1;
526 // Int_t IpaddleT=-1;
527 tof->GetPMTPaddle(pmt_id, Iplane, Ipaddle);
528 tof->GetPaddleGeometry(Iplane,Ipaddle,xleft,xright,yleft,yright);
529 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 !!!!
530 //t_tof->dedx.AddAt((tofdedx->GetdEdx_pmt(pmt_id)*4./dedx_corr[pmt_id]),t_tof->npmtadc);
531 //t_tof->dedx.AddAt((tofdedx->GetdEdx_pmt(pmt_id)*dedx_corr[pmt_id]/4.),t_tof->npmtadc);//annullo wolfrizzazione
532 t_tof->dedx.AddAt((tofdedx->GetdEdx_pmt(pmt_id)),t_tof->npmtadc);//annullo wolfrizzazione
533 t_tof->pmtadc.AddAt(pmt_id,t_tof->npmtadc);
534 t_tof->adcflag.AddAt(0,t_tof->npmtadc); // gf: Jan 09/07
535 t_tof->npmtadc++;
536 };
537 };
538 };
539 };
540
541
542
543
544
545
546
547
548
549
550 // for (Int_t hh=0; hh<12;hh++){
551 // for (Int_t kk=0; kk<4;kk++){
552 // if ( tofoutput_.adctof_c[hh][kk] < 1000 ){
553 // // t_tof->dedx.AddAt(tofoutput_.adctof_c[hh][kk],t_tof->npmtadc); // EMILIANO
554 // pmt_id = tof->GetPMTid(kk,hh);
555 // t_tof->dedx.AddAt((tofoutput_.adctof_c[hh][kk]*4./dedx_corr[pmt_id]),t_tof->npmtadc); // EMILIANO
556 // t_tof->pmtadc.AddAt(pmt_id,t_tof->npmtadc);
557 // t_tof->adcflag.AddAt(tofoutput_.adcflagtof[hh][kk],t_tof->npmtadc); // gf: Jan 09/07
558 // t_tof->npmtadc++;
559 // };
560 // };
561 // };
562 //
563
564 //
565 new(t[ntrkentry]) ToFTrkVar(*t_tof);
566 ntrkentry++;
567 t_tof->Clear();
568 //
569 //
570 //
571 t_pmt->Clear();
572 //
573 for (Int_t hh=0; hh<12;hh++){
574 for (Int_t kk=0; kk<4;kk++){
575 // new WM
576 // if ( tofoutput_.tdc_c[hh][kk] < 4095 || adc[kk][hh] < 4095 || tdc[kk][hh] < 4095 ){
577 if ( tdcc[kk][hh] < 4095. || adc[kk][hh] < 4095 || tdc[kk][hh] < 4095 ){
578 //
579 t_pmt->pmt_id = tof->GetPMTid(kk,hh);
580 t_pmt->tdc_tw = tofoutput_.tdc_c[hh][kk];
581 t_pmt->adc = (Float_t)adc[kk][hh];
582 t_pmt->tdc = (Float_t)tdc[kk][hh];
583 //
584 new(tpmt[npmtentry]) ToFPMT(*t_pmt);
585 npmtentry++;
586 t_pmt->Clear();
587 };
588 };
589 };
590 //
591 // Calculate track-related variables
592 //
593 if ( pam_event->GetTrkLevel2()->ntrk() > 0 ){
594 //
595 // We have at least one track
596 //
597 //
598 // Run over tracks
599 //
600 for(Int_t nt=0; nt < pam_event->GetTrkLevel2()->ntrk(); nt++){
601 //
602 TrkTrack *ptt = pam_event->GetTrkLevel2()->GetStoredTrack(nt);
603 //
604 // Copy the alpha vector in the input structure
605 //
606 for (Int_t e = 0; e < 5 ; e++){
607 tofinput_.al_pp[e] = ptt->al[e];
608 };
609
610 // new input for 9th reduction: tracker dEdx
611 tofinput_.trkmip = ptt->GetDEDX();
612
613 //
614 // Get tracker related variables for this track
615 //
616 toftrk();
617 // toftrk(thelp);
618 //
619 // Copy values in the class from the structure (we need to use a temporary class to store variables).
620 //
621 t_tof->npmttdc = 0;
622 for (Int_t hh=0; hh<12;hh++){
623 for (Int_t kk=0; kk<4;kk++){
624 if ( tofoutput_.tofmask[hh][kk] != 0 ){
625 pmt_id = tof->GetPMTid(kk,hh);
626 t_tof->pmttdc.AddAt(pmt_id,t_tof->npmttdc);
627 t_tof->tdcflag.AddAt(tofoutput_.tdcflag[hh][kk],t_tof->npmttdc); // gf: Jan 09/07
628 t_tof->npmttdc++;
629 };
630 };
631 };
632 for (Int_t kk=0; kk<13;kk++){
633 t_tof->beta[kk] = tofoutput_.beta_a[kk];
634 };
635 //
636 memcpy(t_tof->xtofpos,tofoutput_.xtofpos,sizeof(t_tof->xtofpos));
637 memcpy(t_tof->ytofpos,tofoutput_.ytofpos,sizeof(t_tof->ytofpos));
638 memcpy(t_tof->xtr_tof,tofoutput_.xtr_tof,sizeof(t_tof->xtr_tof));
639 memcpy(t_tof->ytr_tof,tofoutput_.ytr_tof,sizeof(t_tof->ytr_tof));
640 //
641 tofdedx->Process(pam_event->GetOrbitalInfo()->absTime,t_tof->beta[12], (Float_t *)t_tof->xtr_tof,(Float_t *)t_tof->ytr_tof);
642 t_tof->npmtadc = 0;
643
644 for (Int_t hh=0; hh<12;hh++){
645 for (Int_t kk=0; kk<4;kk++){
646 pmt_id = tof->GetPMTid(kk,hh);
647 Int_t Iplane=-1;
648 Int_t Ipaddle=-1;
649 Int_t IpaddleT=-1;
650 tof->GetPMTPaddle(pmt_id, Iplane, Ipaddle);
651 IpaddleT=tof->GetPaddleIdOfTrack(t_tof->xtr_tof[Iplane],t_tof->ytr_tof[Iplane], Iplane,0.0);
652 if ( tofdedx->GetdEdx_pmt(pmt_id) > -1. && Ipaddle==IpaddleT ){
653 //t_tof->dedx.AddAt((tofdedx->GetdEdx_pmt(pmt_id)*4./dedx_corr[pmt_id]),t_tof->npmtadc);
654 //t_tof->dedx.AddAt((tofdedx->GetdEdx_pmt(pmt_id)*dedx_corr[pmt_id]/4.),t_tof->npmtadc);//annullo wolfrizzazione
655 t_tof->dedx.AddAt((tofdedx->GetdEdx_pmt(pmt_id)),t_tof->npmtadc);//annullo wolfrizzazione
656 t_tof->pmtadc.AddAt(pmt_id,t_tof->npmtadc);
657 t_tof->adcflag.AddAt(0,t_tof->npmtadc); // gf: Jan 09/07
658 t_tof->npmtadc++;
659 };
660
661 };
662 };
663
664
665 // for (Int_t hh=0; hh<12;hh++){
666 // for (Int_t kk=0; kk<4;kk++){
667 // pmt_id = tof->GetPMTid(kk,hh);
668 // if ( tofdedx->GetdEdx_pmt(pmt_id) > -1. ){
669 // t_tof->dedx.AddAt((tofdedx->GetdEdx_pmt(pmt_id)*36./pow(dedx_corr[pmt_id],2)),t_tof->npmtadc);
670 // t_tof->pmtadc.AddAt(pmt_id,t_tof->npmtadc);
671 // t_tof->adcflag.AddAt(0,t_tof->npmtadc); // gf: Jan 09/07
672 // 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]);
673 // t_tof->npmtadc++;
674 // };
675
676 // };
677 // };
678 // t_tof->npmtadc = 0;
679 // for (Int_t hh=0; hh<12;hh++){
680 // for (Int_t kk=0; kk<4;kk++){
681 // if ( tofoutput_.adc_c[hh][kk] < 1000 ){
682 // // t_tof->dedx.AddAt(tofoutput_.adc_c[hh][kk],t_tof->npmtadc); // EMILIANO
683 // pmt_id = tof->GetPMTid(kk,hh);
684 // t_tof->dedx.AddAt((tofoutput_.adc_c[hh][kk]*4./dedx_corr[pmt_id]),t_tof->npmtadc); // EMILIANO
685 // t_tof->pmtadc.AddAt(pmt_id,t_tof->npmtadc);
686 // t_tof->adcflag.AddAt(tofoutput_.adcflag[hh][kk],t_tof->npmtadc); // gf: Jan 09/07
687 // t_tof->npmtadc++;
688 // };
689 // };
690 // };
691 //
692
693 //
694 // Store the tracker track number in order to be sure to have shyncronized data during analysis
695 //
696 t_tof->trkseqno = nt;
697 //
698 // create a new object for this event with track-related variables
699 //
700 new(t[ntrkentry]) ToFTrkVar(*t_tof);
701 ntrkentry++;
702 t_tof->Clear();
703 //
704 }; // loop on all the tracks
705 };
706 //
707 tof->unpackError = tofl2->unpackError;
708 if ( defcal ){
709 tof->default_calib = 1;
710 } else {
711 tof->default_calib = 0;
712 };
713 //
714 // Fill the rootple
715 //
716 toft->Fill();
717 //
718 // toft->Show();
719 // toftr->Show();
720 //
721 delete t_tof;
722 //
723 //
724 };
725 //
726 };
727 //
728 pam_event->Clear();
729 //
730 outft->cd();
731 pam_event->WriteCloneTrees();
732 toftr->Delete();
733
734 toft->Write();
735 outft->Close();
736 //
737 return 0;
738 //
739 }
740
741 /////////////////////////////////////////////////////////////////
742
743 #if !defined(__CINT__)
744
745 void usage(){
746
747 cout << "------------------------------------------------------------"<<endl;
748 cout << "Loop over events (on one or more Level2-files), applying some selection cuts (defined in My-Selection.cpp), \n";
749 cout << "creating output histograms (defined in My-Histos.cpp) and/or trees with selected events. \n \n ";
750 cout << "USAGE:"<<endl;
751 cout << "-processDir DIR - Level2 data directory \n";
752 cout << "-processList LIST - list of files (.txt) or single file (.root) to be analysed \n";
753 cout << "-outputFile PATH - name of the output file \n";
754 cout << "-calibFile PATH+NAME - name of the calibration file \n";
755 cout << "-NumEvents XXX - number of events to be analysed \n";
756 cout << "--debug, -g - debug mode \n";
757 cout << "--help, -h - print this help \n";
758 cout << "-options [ options ] - options: \n";
759 cout << " fillHistos --> create an output file with histograms \n";
760 cout << " fillTree --> create an output file with trees storing the selected events \n ";
761 cout << " +(-)ALL --> inlcude(exclude) all trees and branches \n " ;
762 cout << " +(-)TRK1 +(-)TRK2 +(-)CAL1 +(-)CAL2 +(-)TOF +(-)TRG +(-)ND +(-)S4 +(-)ORB --> inlcude(exclude) trees and branches \n" ;
763 cout << "------------------------------------------------------------"<<endl;
764 }
765 //
766 int HandleInputPar(int argc, char **argv){
767
768 if(argc>1){
769
770 if(!strcmp(argv[1], "-h") || !strcmp(argv[1], "--help") ){
771 usage();
772 return(1);
773 };
774 // -----------------------
775 // Read input parameters
776 // -----------------------
777 DIR = gSystem->WorkingDirectory();
778 LIST = "";
779 OUTF = "myfile";
780 OPTIONS = "+AUTO -TOF";
781
782 for (int i = 1; i < argc; i++){
783 // -----------------------------------------------------//
784 if (!strcmp(argv[i], "-processDir")){
785 if (++i >= argc) throw -1;
786 DIR = argv[i];
787 cout << "processDir "<<DIR<<endl;
788 continue;
789 }
790 // -----------------------------------------------------//
791 else if (!strcmp(argv[i], "-processList")){
792 if (++i >= argc) throw -1;
793 LIST = argv[i];
794 cout << "processList "<<LIST<<endl;
795 continue;
796 }
797 // -----------------------------------------------------//
798 else if (!strcmp(argv[i], "-outputFile")){
799 if (++i >= argc) throw -1;
800 OUTF = argv[i];
801 cout << "outputFile "<<OUTF<<endl;
802 continue;
803 }
804 // -----------------------------------------------------//
805 // -----------------------------------------------------//
806 else if (!strcmp(argv[i], "-calibFile")){
807 if (++i >= argc) throw -1;
808 CALIBF = gSystem->ExpandPathName(argv[i]);
809 cout << "calibFile "<<CALIBF<<endl;
810 continue;
811 }
812 // -----------------------------------------------------//
813 else if (!strcmp(argv[i], "-options")){
814 if (++i >= argc) throw -1;
815 OPTIONS = argv[i];
816 if( OPTIONS.Contains("[") ){
817 do{
818 if (++i >= argc) throw -1;
819 OPTIONS.Append(argv[i]);
820 }while(!OPTIONS.Contains("]"));
821 }else cout << "wrong option format --> ignoring " << endl;
822 }
823 else{
824 cout << "Unidentified input parameter. Ignored."<< endl;
825 };
826 };
827 }else{
828 usage();
829 return(1);
830 };
831 // -----------------------
832 // Check input parameters
833 // -----------------------
834
835
836 return(0);
837
838 };
839 //
840
841 int main(int argc, char **argv)
842 {
843
844 if( HandleInputPar(argc,argv) )return(1);
845
846 cout << "OPTIONS "<<OPTIONS<<endl;
847 Loop(DIR,LIST,OPTIONS);
848
849 cout << "Back to main - end"<<endl;
850
851 return 0;
852
853 };
854
855 #endif

  ViewVC Help
Powered by ViewVC 1.1.23