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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.7 - (hide 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 mocchiut 1.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 mocchiut 1.6 //extern "C" int rdtofcal(char [], int *);
48     extern "C" int rdtofcal(const char *, int *);
49 mocchiut 1.1
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 mocchiut 1.3 TString CALIBF;
62 mocchiut 1.1
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 mocchiut 1.2 pam_event->SetSELLI(2);
78 mocchiut 1.1 TTree::SetMaxTreeSize(1000*Long64_t(2000000000));
79 mocchiut 1.2 //
80 mocchiut 1.3 TString outfile_t = "";
81 mocchiut 1.1 outfile_t = OUTF;
82 mocchiut 1.3 outfile_t.Append(".root");
83 mocchiut 1.1 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 mocchiut 1.3 ToFdEdx *tofdedx = new ToFdEdx();
95 mocchiut 1.1 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 mocchiut 1.4 Float_t xleft=0;
105     Float_t xright=0;
106     Float_t yleft=0;
107     Float_t yright=0;
108 mocchiut 1.1 //
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 mocchiut 1.2 //==================================================================
157 mocchiut 1.1
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 mocchiut 1.3
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 mocchiut 1.6 printf(" Reading Attenuation file: %s \n",(glparam->PATH+glparam->NAME).Data());
188    
189 mocchiut 1.3 //
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 mocchiut 1.6 printf(" Reading desaturation1 file: %s \n",(glparam->PATH+glparam->NAME).Data());
197    
198 mocchiut 1.3 //
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 mocchiut 1.6 printf(" Reading BBneg file: %s \n",(glparam->PATH+glparam->NAME).Data());
206    
207 mocchiut 1.3 //
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 mocchiut 1.6 printf(" Reading BBpos file: %s \n",(glparam->PATH+glparam->NAME).Data());
215    
216 mocchiut 1.3 //
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 mocchiut 1.6 printf(" Reading desaturation2 file: %s \n",(glparam->PATH+glparam->NAME).Data());
224 mocchiut 1.3
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 mocchiut 1.1 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 mocchiut 1.6 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 mocchiut 1.1 //
243 mocchiut 1.6
244     tofdedx->Clear();
245    
246 mocchiut 1.1 };
247     //
248    
249    
250 mocchiut 1.2 //================================================================
251     //==================================================================
252     //--- Define absolute time
253 mocchiut 1.3 UInt_t tabs=pam_event->GetOrbitalInfo()->absTime;
254    
255 mocchiut 1.6 if ( !(iev%100000) ) printf(" ATIME %u re %u \n",(Int_t)tabs,(UInt_t)iev);
256 mocchiut 1.2
257     //==================================================================
258 mocchiut 1.7 Float_t dedx_corr_m[2000][48];//,dedx_corr[48];
259 mocchiut 1.2 Double_t mtime[2000],t1,t2,tm;
260 mocchiut 1.7 // Float_t yhelp1,yhelp2,slope,inter,thelp1,thelp2;
261     Float_t yhelp1,yhelp2;
262 mocchiut 1.2 Float_t xmean1,xwidth1;
263    
264 mocchiut 1.7 Int_t ical=0;
265     Int_t ii,j,jj;
266 mocchiut 1.2
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 mocchiut 1.4 printf(" READING NEW CALIBRATION FILE: %s \n",CALIBF.Data());
277 mocchiut 1.2
278 mocchiut 1.3 ifstream fin(CALIBF.Data());
279 mocchiut 1.6 //cout << "topolino" << endl;
280 mocchiut 1.2 while (! fin.eof()) {
281     fin>>t1>>tm>>t2;
282 mocchiut 1.6 //cout<<jj<<" "<<tm<<endl;
283     //cout << setiosflags(ios::fixed) << setw(10) << setprecision(3) << tm << endl;
284 mocchiut 1.2 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 mocchiut 1.6 //cout << "pippo" << endl;
294 mocchiut 1.2 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 mocchiut 1.6 //cout<<"abs time "<<tabs<<" limit low "<<mtime[ical]<<" limit up "<<mtime[ical+1]<<" ical "<<ical<<endl;
303 mocchiut 1.2 } // 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 mocchiut 1.6 //cout<<"abs time "<<tabs<<" limit low "<<mtime[ical]<<" limit up "<<mtime[ical+1]<<" ical "<<ical<<endl;
321 mocchiut 1.2 };
322 mocchiut 1.1
323 mocchiut 1.2 // printf(" 2ical %i \n",ical);
324     //==================================================================
325     //== interpolate betwen time limits
326     //==================================================================
327    
328 mocchiut 1.7 // thelp1 = mtime[ical];
329     // thelp2 = mtime[ical+1];
330 mocchiut 1.2
331     for (ii=0; ii<48;ii++) {
332 mocchiut 1.4 yhelp1 = fabs(dedx_corr_m[ical][ii]);
333 mocchiut 1.3 // yhelp1 = 6.;
334 mocchiut 1.4 if ( yhelp1 < 0.1 ) yhelp1 = 4.;
335     yhelp2 = fabs(dedx_corr_m[ical+1][ii]);
336 mocchiut 1.3 // yhelp2 = 6.;
337 mocchiut 1.4 if ( yhelp2 < 0.1 ) yhelp2 = 4.;
338 mocchiut 1.7 // slope = (yhelp2-yhelp1)/(thelp2-thelp1);
339     // inter = yhelp1 - slope*thelp1;
340     // dedx_corr[ii] = slope*tabs + inter;
341 mocchiut 1.2 // if (ii==0) cout<<thelp1<<" "<<thelp2<<" "<<tabs<<" "<<yhelp1<<" "<<yhelp2<<" "<<dedx_corr[0]<<endl;
342     }
343 mocchiut 1.1
344 mocchiut 1.2 //================================================================
345     //================================================================
346     // printf("cpippo \n");
347 mocchiut 1.1
348     // process tof data
349     //
350     for (Int_t hh=0; hh<12;hh++){
351     for (Int_t kk=0; kk<4;kk++){
352 mocchiut 1.2 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 mocchiut 1.1 };
358     };
359 mocchiut 1.2 // memset(adc, 0, 12*4*sizeof(Int_t));
360     // memset(tdc, 0, 12*4*sizeof(Int_t));
361 mocchiut 1.1 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 mocchiut 1.2 if ( (ttf->tdcflag).At(nc) != 0 ) tdcf[(ttf->pmttdc).At(nc)] = 1;
371 mocchiut 1.1 };
372     for ( Int_t nc=0; nc < ttf->npmtadc; nc++){
373 mocchiut 1.2 if ( (ttf->adcflag).At(nc) != 0 ) adcf[(ttf->pmtadc).At(nc)] = 1;
374 mocchiut 1.1 };
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 mocchiut 1.2 tofinput_.adc[gg][hh] = (int)pmt->adc;
382     adc[hh][gg] = (int)pmt->adc;
383 mocchiut 1.1 };
384     if ( tdcf[pmt->pmt_id] == 0 ){
385 mocchiut 1.2 tofinput_.tdc[gg][hh] = (int)pmt->tdc;
386     tdc[hh][gg] = (int)pmt->tdc;
387 mocchiut 1.1 };
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 mocchiut 1.3 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 mocchiut 1.2 // 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 mocchiut 1.1 //
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 mocchiut 1.3 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 mocchiut 1.4 // {
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 mocchiut 1.3
463    
464     // for (Int_t hh=0; hh<12;hh++){
465     // for (Int_t kk=0; kk<4;kk++){
466 mocchiut 1.4 // 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 mocchiut 1.3 // t_tof->pmtadc.AddAt(pmt_id,t_tof->npmtadc);
470 mocchiut 1.4 // t_tof->adcflag.AddAt(0,t_tof->npmtadc); // gf: Jan 09/07
471 mocchiut 1.3 // t_tof->npmtadc++;
472     // };
473     // };
474     // };
475 mocchiut 1.4 {
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 mocchiut 1.6 //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 mocchiut 1.4 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 mocchiut 1.3
551 mocchiut 1.4
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 mocchiut 1.1 for (Int_t hh=0; hh<12;hh++){
576     for (Int_t kk=0; kk<4;kk++){
577 mocchiut 1.2 // new WM
578     // if ( tofoutput_.tdc_c[hh][kk] < 4095 || adc[kk][hh] < 4095 || tdc[kk][hh] < 4095 ){
579 mocchiut 1.1 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 mocchiut 1.4
612     // new input for 9th reduction: tracker dEdx
613     tofinput_.trkmip = ptt->GetDEDX();
614    
615 mocchiut 1.1 //
616     // Get tracker related variables for this track
617     //
618     toftrk();
619 mocchiut 1.2 // toftrk(thelp);
620 mocchiut 1.1 //
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 mocchiut 1.3 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 mocchiut 1.1 t_tof->npmtadc = 0;
645 mocchiut 1.4
646 mocchiut 1.1 for (Int_t hh=0; hh<12;hh++){
647     for (Int_t kk=0; kk<4;kk++){
648 mocchiut 1.3 pmt_id = tof->GetPMTid(kk,hh);
649 mocchiut 1.4 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 mocchiut 1.6 //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 mocchiut 1.1 t_tof->pmtadc.AddAt(pmt_id,t_tof->npmtadc);
659 mocchiut 1.3 t_tof->adcflag.AddAt(0,t_tof->npmtadc); // gf: Jan 09/07
660 mocchiut 1.1 t_tof->npmtadc++;
661     };
662 mocchiut 1.3
663 mocchiut 1.1 };
664     };
665 mocchiut 1.4
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 mocchiut 1.3 // 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 mocchiut 1.1 //
694 mocchiut 1.3
695 mocchiut 1.1 //
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 mocchiut 1.2 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 mocchiut 1.3 cout << "-calibFile PATH+NAME - name of the calibration file \n";
757 mocchiut 1.2 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 mocchiut 1.1 }
767     //
768     int HandleInputPar(int argc, char **argv){
769    
770 mocchiut 1.2 if(argc>1){
771 mocchiut 1.1
772 mocchiut 1.2 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 mocchiut 1.1
784 mocchiut 1.2 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 mocchiut 1.3 // -----------------------------------------------------//
808     else if (!strcmp(argv[i], "-calibFile")){
809     if (++i >= argc) throw -1;
810 mocchiut 1.5 CALIBF = gSystem->ExpandPathName(argv[i]);
811 mocchiut 1.3 cout << "calibFile "<<CALIBF<<endl;
812     continue;
813     }
814     // -----------------------------------------------------//
815 mocchiut 1.2 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 mocchiut 1.1 };
829 mocchiut 1.2 }else{
830     usage();
831     return(1);
832     };
833     // -----------------------
834     // Check input parameters
835     // -----------------------
836 mocchiut 1.1
837    
838 mocchiut 1.2 return(0);
839 mocchiut 1.1
840     };
841     //
842    
843     int main(int argc, char **argv)
844     {
845    
846 mocchiut 1.2 if( HandleInputPar(argc,argv) )return(1);
847 mocchiut 1.1
848 mocchiut 1.2 cout << "OPTIONS "<<OPTIONS<<endl;
849     Loop(DIR,LIST,OPTIONS);
850 mocchiut 1.1
851 mocchiut 1.2 cout << "Back to main - end"<<endl;
852 mocchiut 1.1
853 mocchiut 1.2 return 0;
854 mocchiut 1.1
855     };
856    
857     #endif

  ViewVC Help
Powered by ViewVC 1.1.23