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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.1.1.1 by mocchiut, Fri Nov 28 14:58:50 2008 UTC revision 1.2 by mocchiut, Thu Feb 5 09:12:35 2009 UTC
# Line 72  Int_t Loop(TString ddir,TString list, TS Line 72  Int_t Loop(TString ddir,TString list, TS
72    //  read input file/list    //  read input file/list
73    //  --------------------    //  --------------------
74    pam_event = new PamLevel2(ddir,list,options);    pam_event = new PamLevel2(ddir,list,options);
75      pam_event->SetSELLI(2);
76    TTree::SetMaxTreeSize(1000*Long64_t(2000000000));    TTree::SetMaxTreeSize(1000*Long64_t(2000000000));
77   //    //
78    TString outfile_t =0;    TString outfile_t =0;
79    outfile_t = OUTF;    outfile_t = OUTF;
80    outfile_t.Append("b.root");    outfile_t.Append("b.root");
# Line 145  Int_t Loop(TString ddir,TString list, TS Line 146  Int_t Loop(TString ddir,TString list, TS
146    for(ULong64_t iev=0; iev<nevents; iev++){          for(ULong64_t iev=0; iev<nevents; iev++){      
147    
148    
149  //==================================================================      //==================================================================
150    
151    
152      //  for(ULong64_t iev=0; iev<50; iev++){            //  for(ULong64_t iev=0; iev<50; iev++){      
# Line 182  Int_t Loop(TString ddir,TString list, TS Line 183  Int_t Loop(TString ddir,TString list, TS
183        //        //
184    
185    
186  //================================================================        //================================================================
187  //==================================================================        //==================================================================
188  //---  Define absolute time        //---  Define absolute time
189         Float_t tabs=pam_event->GetOrbitalInfo()->absTime;        Float_t tabs=pam_event->GetOrbitalInfo()->absTime;
190    
191  //==================================================================        //==================================================================
192  Float_t dedx_corr_m[2000][48],dedx_corr[48];        Float_t dedx_corr_m[2000][48],dedx_corr[48];
193  Double_t mtime[2000],t1,t2,tm;        Double_t mtime[2000],t1,t2,tm;
194  Float_t yhelp1,yhelp2,slope,inter,thelp1,thelp2;        Float_t yhelp1,yhelp2,slope,inter,thelp1,thelp2;
195  Float_t xmean1,xwidth1;        Float_t xmean1,xwidth1;
196    
197  Int_t ical,ii,j,jj;        Int_t ical,ii,j,jj;
198    
199  if (iev==0) {        if (iev==0) {
200    
201  ical=0;  // counter set to zero if first-time reading          ical=0;  // counter set to zero if first-time reading
202    
203  //-----------------------------------------------------------          //-----------------------------------------------------------
204  // Here I  read the dEdx_korr parameters          // Here I  read the dEdx_korr parameters
205  //-----------------------------------------------------------          //-----------------------------------------------------------
206    
207  jj=0;          jj=0;
208    
209  ifstream fin("adcmonitor.7th.100k.cut.dat");          ifstream fin("adcmonitor.7th.100k.cut.dat");
210    
211  while (! fin.eof()) {          while (! fin.eof()) {
212  fin>>t1>>tm>>t2;            fin>>t1>>tm>>t2;
213  //cout<<jj<<" "<<tm<<endl;            //      cout<<jj<<" "<<tm<<endl;
214  cout << setiosflags(ios::fixed)  << setw(10) << setprecision(3) << tm << endl;            cout << setiosflags(ios::fixed)  << setw(10) << setprecision(3) << tm << endl;
215  mtime[jj]=tm;            mtime[jj]=tm;
216  for (ii=0; ii<48;ii++) {            for (ii=0; ii<48;ii++) {
217  fin>>j>>xmean1>>xwidth1;              fin>>j>>xmean1>>xwidth1;
218  dedx_corr_m[jj][ii]=xmean1;              dedx_corr_m[jj][ii]=xmean1;
219                         }            }
220  jj=jj+1;            jj=jj+1;
221  }            //      printf(" kk %i \n",jj);
222            }
223  fin.close();          //      printf(" 1ical %i \n",ical);
224    
225            fin.close();
226  while (tabs > mtime[ical]) {  
227  ical = ical+1;  
228            }          while (tabs<mtime[ical] || tabs > mtime[ical+1]) {
229  ical = ical-1;            //      printf(" ical %i \n",ical);
230  cout<<"abs time "<<tabs<<" limit low "<<mtime[ical]<<" limit up "<<mtime[ical+1]<<" ical "<<ical<<endl;            ical = ical+1;
231  } // if iev==0...          }
232            //      ical = ical-1;
233  //==================================================================          cout<<"abs time "<<tabs<<" limit low "<<mtime[ical]<<" limit up "<<mtime[ical+1]<<" ical "<<ical<<endl;
234  //==  End of first time reading & filling the array        } // if iev==0...
235  //==================================================================  
236          //==================================================================
237  //==================================================================        //==  End of first time reading & filling the array
238  //==  if time is outside time limits:        //==================================================================
239  //==================================================================  
240          //==================================================================
241  if (tabs<mtime[ical] || tabs>mtime[ical+1]) {        //==  if time is outside time limits:
242  cout<<"Checking Time Limits!"<<endl;        //==================================================================
243  ical=0;  
244  while (tabs > mtime[ical]) {        if (tabs<mtime[ical] || tabs>mtime[ical+1]) {
245  ical = ical+1;          cout<<"Checking Time Limits!"<<endl;
246            }          ical=0;
247  ical = ical-1;          while (tabs > mtime[ical+1] || tabs<mtime[ical]) {
248  cout<<"abs time "<<tabs<<" limit low "<<mtime[ical]<<" limit up "<<mtime[ical+1]<<" ical "<<ical<<endl;            ical = ical+1;
249                                             };          }
250            //      ical = ical-1;
251  //==================================================================          cout<<"abs time "<<tabs<<" limit low "<<mtime[ical]<<" limit up "<<mtime[ical+1]<<" ical "<<ical<<endl;
252  //== interpolate betwen time limits        };
 //==================================================================  
   
        thelp1 = mtime[ical];  
        thelp2 = mtime[ical+1];  
   
        for (ii=0; ii<48;ii++) {  
        yhelp1 = dedx_corr_m[ical][ii];  
        yhelp2 = dedx_corr_m[ical+1][ii];  
        slope  = (yhelp2-yhelp1)/(thelp2-thelp1);  
        inter  = yhelp1 - slope*thelp1;  
        dedx_corr[ii] = slope*tabs + inter;  
 //       if (ii==0) cout<<thelp1<<" "<<thelp2<<" "<<tabs<<" "<<yhelp1<<" "<<yhelp2<<" "<<dedx_corr[0]<<endl;  
                        }  
253    
254  //================================================================        //      printf(" 2ical %i \n",ical);
255  //================================================================        //==================================================================
256          //== interpolate betwen time limits
257          //==================================================================
258    
259          thelp1 = mtime[ical];
260          thelp2 = mtime[ical+1];
261    
262          for (ii=0; ii<48;ii++) {
263            yhelp1 = dedx_corr_m[ical][ii];
264            yhelp2 = dedx_corr_m[ical+1][ii];
265            slope  = (yhelp2-yhelp1)/(thelp2-thelp1);
266            inter  = yhelp1 - slope*thelp1;
267            dedx_corr[ii] = slope*tabs + inter;
268            //       if (ii==0) cout<<thelp1<<" "<<thelp2<<" "<<tabs<<" "<<yhelp1<<" "<<yhelp2<<" "<<dedx_corr[0]<<endl;
269          }
270    
271          //================================================================
272          //================================================================
273          //      printf("cpippo \n");
274    
275        // process tof data        // process tof data
276        //        //
277        for (Int_t hh=0; hh<12;hh++){        for (Int_t hh=0; hh<12;hh++){
278          for (Int_t kk=0; kk<4;kk++){          for (Int_t kk=0; kk<4;kk++){
279                  adc[kk][hh] = 4095;            adc[kk][hh] = 4095;
280                  tdc[kk][hh] = 4095;            tdc[kk][hh] = 4095;
281                  tdcc[kk][hh] = 4095.;            tdcc[kk][hh] = 4095.;
282                  tofinput_.adc[hh][kk] = 4095;            tofinput_.adc[hh][kk] = 4095;
283                  tofinput_.tdc[hh][kk] = 4095;            tofinput_.tdc[hh][kk] = 4095;
284          };          };
285        };        };
286  //      memset(adc, 0, 12*4*sizeof(Int_t));        //      memset(adc, 0, 12*4*sizeof(Int_t));
287  //      memset(tdc, 0, 12*4*sizeof(Int_t));        //      memset(tdc, 0, 12*4*sizeof(Int_t));
288        Int_t gg = 0;        Int_t gg = 0;
289        Int_t hh = 0;        Int_t hh = 0;
290        Int_t adcf[48];        Int_t adcf[48];
# Line 289  cout<<"abs time "<<tabs<<" limit low "<< Line 294  cout<<"abs time "<<tabs<<" limit low "<<
294        for (Int_t pm=0; pm < tofl2->ntrk() ; pm++){        for (Int_t pm=0; pm < tofl2->ntrk() ; pm++){
295          ToFTrkVar *ttf = tofl2->GetToFTrkVar(pm);          ToFTrkVar *ttf = tofl2->GetToFTrkVar(pm);
296          for ( Int_t nc=0; nc < ttf->npmttdc; nc++){          for ( Int_t nc=0; nc < ttf->npmttdc; nc++){
297                  if ( (ttf->tdcflag).At(nc) != 0 ) tdcf[(ttf->pmttdc).At(nc)] = 1;            if ( (ttf->tdcflag).At(nc) != 0 ) tdcf[(ttf->pmttdc).At(nc)] = 1;
298          };          };
299          for ( Int_t nc=0; nc < ttf->npmtadc; nc++){          for ( Int_t nc=0; nc < ttf->npmtadc; nc++){
300                  if ( (ttf->adcflag).At(nc) != 0 ) adcf[(ttf->pmtadc).At(nc)] = 1;            if ( (ttf->adcflag).At(nc) != 0 ) adcf[(ttf->pmtadc).At(nc)] = 1;
301          };          };
302        };        };
303        //        //
# Line 300  cout<<"abs time "<<tabs<<" limit low "<< Line 305  cout<<"abs time "<<tabs<<" limit low "<<
305          ToFPMT *pmt = tofl2->GetToFPMT(pm);          ToFPMT *pmt = tofl2->GetToFPMT(pm);
306          tofl2->GetPMTIndex(pmt->pmt_id, gg, hh);          tofl2->GetPMTIndex(pmt->pmt_id, gg, hh);
307          if ( adcf[pmt->pmt_id] == 0 ){          if ( adcf[pmt->pmt_id] == 0 ){
308                  tofinput_.adc[gg][hh] = (int)pmt->adc;            tofinput_.adc[gg][hh] = (int)pmt->adc;
309                  adc[hh][gg] = (int)pmt->adc;            adc[hh][gg] = (int)pmt->adc;
310          };          };
311          if ( tdcf[pmt->pmt_id] == 0 ){          if ( tdcf[pmt->pmt_id] == 0 ){
312                  tofinput_.tdc[gg][hh] = (int)pmt->tdc;            tofinput_.tdc[gg][hh] = (int)pmt->tdc;
313                  tdc[hh][gg] = (int)pmt->tdc;            tdc[hh][gg] = (int)pmt->tdc;
314          };          };
315          tdcc[hh][gg] = (float)pmt->tdc_tw;          tdcc[hh][gg] = (float)pmt->tdc_tw;
316          //        Int_t pppid = tof->GetPMTid(hh,gg);          //        Int_t pppid = tof->GetPMTid(hh,gg);
# Line 313  cout<<"abs time "<<tabs<<" limit low "<< Line 318  cout<<"abs time "<<tabs<<" limit low "<<
318        };        };
319    
320    
321  //      for (Int_t pm=0; pm <48 ; pm++){        //      for (Int_t pm=0; pm <48 ; pm++){
322  //      tof->GetPMTIndex(pm, gg, hh);        //        tof->GetPMTIndex(pm, gg, hh);
323  //      tofinput_.tdc[hh][gg] = (int)500.;        //        tofinput_.tdc[hh][gg] = (int)500.;
324  //      tofinput_.adc[hh][gg] = (int)500.;        //        tofinput_.adc[hh][gg] = (int)500.;
325  //      tdc[hh][gg] = (int)500.;        //        tdc[hh][gg] = (int)500.;
326  //      adc[hh][gg] = (int)500.;        //        adc[hh][gg] = (int)500.;
327  //      //      printf(" hh %i gg %i tdc %f adc %f \n",hh,gg,pmt->tdc,pmt->adc);        //        //      printf(" hh %i gg %i tdc %f adc %f \n",hh,gg,pmt->tdc,pmt->adc);
328  //    };        //    };
329        //        //
330        for (Int_t hh=0; hh<5;hh++){        for (Int_t hh=0; hh<5;hh++){
331          tofinput_.patterntrig[hh]=pam_event->GetTrigLevel2()->patterntrig[hh];          tofinput_.patterntrig[hh]=pam_event->GetTrigLevel2()->patterntrig[hh];
# Line 397  cout<<"abs time "<<tabs<<" limit low "<< Line 402  cout<<"abs time "<<tabs<<" limit low "<<
402        //        //
403        for (Int_t hh=0; hh<12;hh++){        for (Int_t hh=0; hh<12;hh++){
404          for (Int_t kk=0; kk<4;kk++){          for (Int_t kk=0; kk<4;kk++){
405           // new WM            // new WM
406  //        if ( tofoutput_.tdc_c[hh][kk] < 4095 || adc[kk][hh] < 4095  || tdc[kk][hh] < 4095 ){                //      if ( tofoutput_.tdc_c[hh][kk] < 4095 || adc[kk][hh] < 4095  || tdc[kk][hh] < 4095 ){    
407            if ( tdcc[kk][hh] < 4095. || adc[kk][hh] < 4095  || tdc[kk][hh] < 4095 ){                    if ( tdcc[kk][hh] < 4095. || adc[kk][hh] < 4095  || tdc[kk][hh] < 4095 ){        
408              //              //
409              t_pmt->pmt_id = tof->GetPMTid(kk,hh);              t_pmt->pmt_id = tof->GetPMTid(kk,hh);
# Line 435  cout<<"abs time "<<tabs<<" limit low "<< Line 440  cout<<"abs time "<<tabs<<" limit low "<<
440            // Get tracker related variables for this track            // Get tracker related variables for this track
441            //            //
442            toftrk();            toftrk();
443  //        toftrk(thelp);            //      toftrk(thelp);
444            //            //
445            // Copy values in the class from the structure (we need to use a temporary class to store variables).            // Copy values in the class from the structure (we need to use a temporary class to store variables).
446            //            //
# Line 458  cout<<"abs time "<<tabs<<" limit low "<< Line 463  cout<<"abs time "<<tabs<<" limit low "<<
463            for (Int_t hh=0; hh<12;hh++){            for (Int_t hh=0; hh<12;hh++){
464              for (Int_t kk=0; kk<4;kk++){              for (Int_t kk=0; kk<4;kk++){
465                if ( tofoutput_.adc_c[hh][kk] < 1000 ){                if ( tofoutput_.adc_c[hh][kk] < 1000 ){
466               //     t_tof->dedx.AddAt(tofoutput_.adc_c[hh][kk],t_tof->npmtadc); // EMILIANO                  //          t_tof->dedx.AddAt(tofoutput_.adc_c[hh][kk],t_tof->npmtadc); // EMILIANO
467                  pmt_id = tof->GetPMTid(kk,hh);                  pmt_id = tof->GetPMTid(kk,hh);
468                  t_tof->dedx.AddAt((tofoutput_.adc_c[hh][kk]*4./dedx_corr[pmt_id]),t_tof->npmtadc); // EMILIANO                  t_tof->dedx.AddAt((tofoutput_.adc_c[hh][kk]*4./dedx_corr[pmt_id]),t_tof->npmtadc); // EMILIANO
469                  t_tof->pmtadc.AddAt(pmt_id,t_tof->npmtadc);                  t_tof->pmtadc.AddAt(pmt_id,t_tof->npmtadc);
# Line 526  cout<<"abs time "<<tabs<<" limit low "<< Line 531  cout<<"abs time "<<tabs<<" limit low "<<
531    
532  void usage(){  void usage(){
533    
534      cout << "------------------------------------------------------------"<<endl;    cout << "------------------------------------------------------------"<<endl;
535      cout << "Loop over events (on one or more Level2-files), applying some selection cuts (defined in My-Selection.cpp), \n";    cout << "Loop over events (on one or more Level2-files), applying some selection cuts (defined in My-Selection.cpp), \n";
536      cout << "creating output histograms (defined in My-Histos.cpp) and/or trees with selected events. \n \n ";    cout << "creating output histograms (defined in My-Histos.cpp) and/or trees with selected events. \n \n ";
537      cout << "USAGE:"<<endl;    cout << "USAGE:"<<endl;
538      cout << "-processDir  DIR     -  Level2 data directory \n";    cout << "-processDir  DIR     -  Level2 data directory \n";
539      cout << "-processList LIST    -  list  of files (.txt) or single file (.root) to be analysed \n";    cout << "-processList LIST    -  list  of files (.txt) or single file (.root) to be analysed \n";
540      cout << "-outputFile  PATH    -  name of the output file \n";    cout << "-outputFile  PATH    -  name of the output file \n";
541      cout << "-NumEvents   XXX     -  number of events to be analysed \n";    cout << "-NumEvents   XXX     -  number of events to be analysed \n";
542      cout << "--debug, -g          -  debug mode \n";    cout << "--debug, -g          -  debug mode \n";
543      cout << "--help, -h           -  print this help \n";    cout << "--help, -h           -  print this help \n";
544      cout << "-options [ options ] -  options: \n";    cout << "-options [ options ] -  options: \n";
545      cout << "                        fillHistos --> create an output file with histograms  \n";    cout << "                        fillHistos --> create an output file with histograms  \n";
546      cout << "                        fillTree   --> create an output file with trees storing the selected events \n ";    cout << "                        fillTree   --> create an output file with trees storing the selected events \n ";
547      cout << "                        +(-)ALL    --> inlcude(exclude) all trees and branches \n "   ;    cout << "                        +(-)ALL    --> inlcude(exclude) all trees and branches \n "   ;
548      cout << "                        +(-)TRK1 +(-)TRK2 +(-)CAL1 +(-)CAL2 +(-)TOF +(-)TRG +(-)ND +(-)S4 +(-)ORB --> inlcude(exclude) trees and branches  \n"  ;    cout << "                        +(-)TRK1 +(-)TRK2 +(-)CAL1 +(-)CAL2 +(-)TOF +(-)TRG +(-)ND +(-)S4 +(-)ORB --> inlcude(exclude) trees and branches  \n"  ;
549      cout << "------------------------------------------------------------"<<endl;    cout << "------------------------------------------------------------"<<endl;
550  }  }
551  //  //
552  int HandleInputPar(int argc, char **argv){  int HandleInputPar(int argc, char **argv){
553    
554      if(argc>1){    if(argc>1){
555                                    
556          if(!strcmp(argv[1], "-h") || !strcmp(argv[1], "--help") ){      if(!strcmp(argv[1], "-h") || !strcmp(argv[1], "--help") ){
557              usage();        usage();
558              return(1);        return(1);
559          };      };
560  // -----------------------      // -----------------------
561  // Read input parameters      // Read input parameters
562  // -----------------------        // -----------------------  
563          DIR     = gSystem->WorkingDirectory();      DIR     = gSystem->WorkingDirectory();
564          LIST    = "";      LIST    = "";
565          OUTF    = "myfile";      OUTF    = "myfile";
566          OPTIONS = "+AUTO -TOF";      OPTIONS = "+AUTO -TOF";
567                                    
568          for (int i = 1; i < argc; i++){              for (int i = 1; i < argc; i++){            
569              // -----------------------------------------------------//        // -----------------------------------------------------//
570              if (!strcmp(argv[i], "-processDir")){        if (!strcmp(argv[i], "-processDir")){
571                  if (++i >= argc) throw -1;          if (++i >= argc) throw -1;
572                  DIR = argv[i];          DIR = argv[i];
573                  cout << "processDir "<<DIR<<endl;          cout << "processDir "<<DIR<<endl;
574                  continue;          continue;
575              }        }
576              // -----------------------------------------------------//        // -----------------------------------------------------//
577              else if (!strcmp(argv[i], "-processList")){        else if (!strcmp(argv[i], "-processList")){
578                  if (++i >= argc) throw -1;          if (++i >= argc) throw -1;
579                  LIST = argv[i];          LIST = argv[i];
580                  cout << "processList "<<LIST<<endl;          cout << "processList "<<LIST<<endl;
581                  continue;          continue;
582              }          }  
583              // -----------------------------------------------------//        // -----------------------------------------------------//
584              else if (!strcmp(argv[i], "-outputFile")){        else if (!strcmp(argv[i], "-outputFile")){
585                  if (++i >= argc) throw -1;          if (++i >= argc) throw -1;
586                  OUTF = argv[i];          OUTF = argv[i];
587                  cout << "outputFile "<<OUTF<<endl;          cout << "outputFile "<<OUTF<<endl;
588                  continue;          continue;
589              }          }  
590              // -----------------------------------------------------//        // -----------------------------------------------------//
591              else if (!strcmp(argv[i], "-options")){        else if (!strcmp(argv[i], "-options")){
592                  if (++i >= argc) throw -1;          if (++i >= argc) throw -1;
593                  OPTIONS = argv[i];          OPTIONS = argv[i];
594                  if( OPTIONS.Contains("[") ){          if( OPTIONS.Contains("[") ){
595                      do{            do{
596                          if (++i >= argc) throw -1;              if (++i >= argc) throw -1;
597                          OPTIONS.Append(argv[i]);              OPTIONS.Append(argv[i]);
598                      }while(!OPTIONS.Contains("]"));            }while(!OPTIONS.Contains("]"));
599                  }else cout << "wrong option format --> ignoring " << endl;          }else cout << "wrong option format --> ignoring " << endl;
600              }        }
601              else{        else{
602                  cout << "Unidentified input parameter. Ignored."<< endl;          cout << "Unidentified input parameter. Ignored."<< endl;
603              };        };
         };  
     }else{  
         usage();  
         return(1);  
604      };      };
605  // -----------------------    }else{
606  // Check input parameters      usage();
607  // -----------------------        return(1);
608      };
609      // -----------------------
610      // Check input parameters
611      // -----------------------  
612                    
613                                    
614      return(0);    return(0);
615                    
616  };  };
617  //  //
# Line 614  int HandleInputPar(int argc, char **argv Line 619  int HandleInputPar(int argc, char **argv
619  int main(int argc, char **argv)  int main(int argc, char **argv)
620  {  {
621                    
622      if( HandleInputPar(argc,argv) )return(1);    if( HandleInputPar(argc,argv) )return(1);
623                    
624      cout << "OPTIONS "<<OPTIONS<<endl;    cout << "OPTIONS "<<OPTIONS<<endl;
625      Loop(DIR,LIST,OPTIONS);    Loop(DIR,LIST,OPTIONS);
626            
627      cout << "Back to main - end"<<endl;    cout << "Back to main - end"<<endl;
628            
629      return 0;    return 0;
630                    
631  };  };
632    
Legend:
Removed from v.1.1.1.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.23