/[PAMELA software]/DarthVader/OrbitalInfo/src/OrbitalInfoCore.cpp
ViewVC logotype

Diff of /DarthVader/OrbitalInfo/src/OrbitalInfoCore.cpp

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

revision 1.4 by mocchiut, Fri Aug 4 10:31:29 2006 UTC revision 1.10 by mocchiut, Wed Jan 24 13:16:26 2007 UTC
# Line 13  Line 13 
13  #include <TClassEdit.h>  #include <TClassEdit.h>
14  #include <TObject.h>  #include <TObject.h>
15  #include <TList.h>  #include <TList.h>
16  #include <TArrayL.h>  #include <TArrayI.h>
17  #include <TSystem.h>  #include <TSystem.h>
18  #include <TSystemDirectory.h>  #include <TSystemDirectory.h>
19  #include <TString.h>  #include <TString.h>
# Line 23  Line 23 
23  #include <TSQLRow.h>  #include <TSQLRow.h>
24  #include <TSQLResult.h>  #include <TSQLResult.h>
25  //  //
26    // RunInfo header
27    //
28    #include <RunInfo.h>
29    #include <GLTables.h>
30    //
31  // YODA headers  // YODA headers
32  //  //
33  #include <PamelaRun.h>  #include <PamelaRun.h>
 #include <RegistryEvent.h>  
34  #include <PscuHeader.h>  #include <PscuHeader.h>
35  #include <PscuEvent.h>  #include <PscuEvent.h>
36  #include <EventHeader.h>  #include <EventHeader.h>
 #include <RegistryEvent.h>  
 //  
 // RunInfo header  
 //  
 #include <RunInfo.h>  
 #include <GLTables.h>  
37  //  //
38  // This program headers  // This program headers
39  //  //
40  #include <OrbitalInfo.h>  #include <OrbitalInfo.h>
 #include <OrbitalInfoCore.h>  
41  #include <OrbitalInfoVerl2.h>  #include <OrbitalInfoVerl2.h>
42    #include <OrbitalInfoCore.h>
43    
44  using namespace std;  using namespace std;
45    
# Line 49  using namespace std; Line 47  using namespace std;
47  // CORE ROUTINE  // CORE ROUTINE
48  //  //
49  //  //
50    int OrbitalInfoCore(UInt_t run, TFile *file, TSQLServer *dbc, Int_t OrbitalInfoargc, char *OrbitalInfoargv[]){
51    
52    //   // Temporary check to use igrf routines.  We need dat files in the
53    //   // current directory.
54    //   fstream igrfdat1("igrf05.dat");
55    //   fstream igrfdat2("igrf05s.dat");
56    //   if( (!igrfdat1) && (!igrfdat2)) {
57    //     cerr<<"\n**************************************\n"
58    //      <<"igrf05.dat or igrf05s.dat not in the current directory.  Exiting.\n"
59    //      <<"**************************************\n";
60    //     exit(EXIT_FAILURE);
61    //   }
62    //   igrfdat1.close();
63    //   igrfdat2.close();
64    //   // end of temporary code
65    //   //
66    
 int OrbitalInfoCore(ULong64_t run, TFile *file, TSQLServer *dbc, Int_t OrbitalInfoargc, char *OrbitalInfoargv[]){  
67    Int_t i = 0;    Int_t i = 0;
68    //    //
69    TString processFolder = "OrbitalInfoFolder";    TString processFolder = "OrbitalInfoFolder";
# Line 84  int OrbitalInfoCore(ULong64_t run, TFile Line 97  int OrbitalInfoCore(ULong64_t run, TFile
97    const char* outDir = gSystem->DirName(gSystem->DirName(file->GetPath()));    const char* outDir = gSystem->DirName(gSystem->DirName(file->GetPath()));
98    //    //
99    TTree *OrbitalInfotr = 0;    TTree *OrbitalInfotr = 0;
100    Long64_t nevents = 0LL;    UInt_t nevents = 0;
101    //    //
102    // variables needed to reprocess data    // variables needed to reprocess data
103    //    //
104      Long64_t maxsize = 10000000000LL;  
105      TTree::SetMaxTreeSize(maxsize);
106      //
107    TString OrbitalInfoversion;    TString OrbitalInfoversion;
108    ItoRunInfo *runinfo = 0;    ItoRunInfo *runinfo = 0;
109    TArrayL *runlist = 0;    TArrayI *runlist = 0;
110    TTree *OrbitalInfotrclone = 0;    TTree *OrbitalInfotrclone = 0;
111    Bool_t reproc = false;    Bool_t reproc = false;
112    Bool_t reprocall = false;    Bool_t reprocall = false;
# Line 99  int OrbitalInfoCore(ULong64_t run, TFile Line 115  int OrbitalInfoCore(ULong64_t run, TFile
115    UInt_t numbofrun = 0;    UInt_t numbofrun = 0;
116    stringstream ftmpname;    stringstream ftmpname;
117    TString fname;    TString fname;
118    Long64_t totfileentries = 0ULL;    UInt_t totfileentries = 0;
119    Long64_t idRun = 0LL;    UInt_t idRun = 0;
120    //    //
121    // variables needed to handle error signals    // variables needed to handle error signals
122    //    //
# Line 116  int OrbitalInfoCore(ULong64_t run, TFile Line 132  int OrbitalInfoCore(ULong64_t run, TFile
132    //    //
133    TFile *l0File = 0;    TFile *l0File = 0;
134    TTree *l0tr = 0;    TTree *l0tr = 0;
   TBranch *l0registry = 0;  
   pamela::RegistryEvent *l0reg=0;  
135    // EM: open also header branch    // EM: open also header branch
136    TBranch *l0head = 0;    TBranch *l0head = 0;
137    pamela::EventHeader *eh = 0;    pamela::EventHeader *eh = 0;
# Line 131  int OrbitalInfoCore(ULong64_t run, TFile Line 145  int OrbitalInfoCore(ULong64_t run, TFile
145    stringstream file3;    stringstream file3;
146    stringstream qy;    stringstream qy;
147    Int_t totevent = 0;    Int_t totevent = 0;
148    ULong64_t atime = 0ULL;    UInt_t atime = 0;
149    Int_t ei = 0;    UInt_t re = 0;
150    Int_t re = 0;  
151      // Position
152      Float_t lon, lat, alt;
153    
154      //
155      // IGRF stuff
156      //
157      float dimo = 0.0; // dipole moment (computed from dat files)
158      float bnorth, beast, bdown, babs;
159      float xl; // L value
160      float icode; // code value for L accuracy (see fortran code)
161      float bab1; // What's  the difference with babs?
162      float stps = 0.005; // step size for field line tracing
163      float bdel = 0.01; // required accuracy
164      float bequ;  // equatorial b value (also called b_0)
165      bool value = 0; // false if bequ is not the minimum b value
166      float rr0; // equatorial radius normalized to earth radius
167    
168    //    //
169    // Working filename    // Working filename
170    //    //
# Line 160  int OrbitalInfoCore(ULong64_t run, TFile Line 191  int OrbitalInfoCore(ULong64_t run, TFile
191    // DB classes    // DB classes
192    //    //
193    GL_ROOT *glroot = new GL_ROOT();    GL_ROOT *glroot = new GL_ROOT();
194      GL_TIMESYNC *dbtime = 0;
195      GL_TLE *gltle = new GL_TLE();
196      // Initialize fortran routines!!!
197      Int_t ltp2 = 0;
198      Int_t ltp3 = 0;
199      Int_t uno = 1;
200      char *niente = " ";
201      GL_PARAM *glparam = new GL_PARAM();
202      GL_PARAM *glparam2 = new GL_PARAM();
203      Int_t parerror=glparam->Query_GL_PARAM(1,301,dbc); // parameters stored in DB in GL_PRAM table
204      if ( parerror<0 ) {
205        code = parerror;
206        goto closeandexit;
207      };
208      ltp2 = (Int_t)(glparam->PATH+glparam->NAME).Length();
209      if ( verbose ) printf(" Reading Earth's Magnetic Field parameter file: %s \n",(glparam->PATH+glparam->NAME).Data());
210      //
211      parerror=glparam2->Query_GL_PARAM(1,302,dbc); // parameters stored in DB in GL_PRAM table
212      if ( parerror<0 ) {
213        code = parerror;
214        goto closeandexit;
215      };
216      ltp3 = (Int_t)(glparam2->PATH+glparam2->NAME).Length();
217      if ( verbose ) printf(" Reading Earth's Magnetic Field parameter file: %s \n",(glparam2->PATH+glparam2->NAME).Data());
218      //
219      initize_((char *)niente,&uno,(char *)(glparam->PATH+glparam->NAME).Data(),&ltp2,(char *)(glparam2->PATH+glparam2->NAME).Data(),&ltp3);
220    
221      //
222      // End IGRF stuff//
223      //
224    
225    //    //
226    // Let's start!    // Let's start!
227    //    //
# Line 167  int OrbitalInfoCore(ULong64_t run, TFile Line 229  int OrbitalInfoCore(ULong64_t run, TFile
229    // if run != 0 we must process only that run but first we have to check if the tree MyDetector2 already exist in the file    // if run != 0 we must process only that run but first we have to check if the tree MyDetector2 already exist in the file
230    // if it exists we are reprocessing data and we must delete that entries, if not we must create it.    // if it exists we are reprocessing data and we must delete that entries, if not we must create it.
231    //    //
232    if ( run == 0ULL )  reproc = true;    if ( run == 0 )  reproc = true;
233    //    //
234    //    //
235    // Output file is "outputfile"    // Output file is "outputfile"
# Line 224  int OrbitalInfoCore(ULong64_t run, TFile Line 286  int OrbitalInfoCore(ULong64_t run, TFile
286      // tree does not exist, we are not reprocessing      // tree does not exist, we are not reprocessing
287      //      //
288      reproc = false;      reproc = false;
289      if ( run == 0ULL ){      if ( run == 0 ){
290        if (verbose) printf(" OrbitalInfo - WARNING: you are reprocessing data but OrbitalInfo tree does not exist!\n");        if (verbose) printf(" OrbitalInfo - WARNING: you are reprocessing data but OrbitalInfo tree does not exist!\n");
291      }      }
292      if ( runinfo->IsReprocessing() && run != 0ULL ) {      if ( runinfo->IsReprocessing() && run != 0 ) {
293        if (verbose) printf(" OrbitalInfo - WARNING: it seems you are not reprocessing data but OrbitalInfo\n versioning information already exists in RunInfo.\n");        if (verbose) printf(" OrbitalInfo - WARNING: it seems you are not reprocessing data but OrbitalInfo\n versioning information already exists in RunInfo.\n");
294      }      }
295    } else {    } else {
296      //      //
297      // tree exists, we are reprocessing data. Are we reprocessing a single run or all the file?      // tree exists, we are reprocessing data. Are we reprocessing a single run or all the file?
298      //      //
299        OrbitalInfotrclone->SetAutoSave(900000000000000LL);
300      reproc = true;      reproc = true;
301      //      //
302      //      //
303      if (verbose) printf("\n Preparing the pre-processing...\n");      if (verbose) printf("\n Preparing the pre-processing...\n");
304      //      //
305      if ( run == 0ULL ){      if ( run == 0 ){
306        //        //
307        // we are reprocessing all the file        // we are reprocessing all the file
308        // if we are reprocessing everything we don't need to copy any old event and we can just work with the new tree and delete the old one immediately        // if we are reprocessing everything we don't need to copy any old event and we can just work with the new tree and delete the old one immediately
# Line 254  int OrbitalInfoCore(ULong64_t run, TFile Line 317  int OrbitalInfoCore(ULong64_t run, TFile
317        //        //
318        reprocall = false;        reprocall = false;
319        //        //
320        if (verbose) printf("\n OrbitalInfo - WARNING: Reprocessing run number %llu \n",run);        if (verbose) printf("\n OrbitalInfo - WARNING: Reprocessing run number %u \n",run);
321        //        //
322        // copying old tree to a new file        // copying old tree to a new file
323        //        //
# Line 277  int OrbitalInfoCore(ULong64_t run, TFile Line 340  int OrbitalInfoCore(ULong64_t run, TFile
340    //    //
341    file->cd();    file->cd();
342    OrbitalInfotr = new TTree("OrbitalInfo-new","PAMELA OrbitalInfo data");    OrbitalInfotr = new TTree("OrbitalInfo-new","PAMELA OrbitalInfo data");
343      OrbitalInfotr->SetAutoSave(900000000000000LL);
344    OrbitalInfotr->Branch("OrbitalInfo","OrbitalInfo",&orbitalinfo);    OrbitalInfotr->Branch("OrbitalInfo","OrbitalInfo",&orbitalinfo);
345    //    //
346    if ( reproc && !reprocall ){    if ( reproc && !reprocall ){
# Line 285  int OrbitalInfoCore(ULong64_t run, TFile Line 349  int OrbitalInfoCore(ULong64_t run, TFile
349      //      //
350      tempfile = new TFile(tempname.str().c_str(),"READ");      tempfile = new TFile(tempname.str().c_str(),"READ");
351      OrbitalInfotrclone = (TTree*)tempfile->Get("OrbitalInfo-old");      OrbitalInfotrclone = (TTree*)tempfile->Get("OrbitalInfo-old");
352        OrbitalInfotrclone->SetAutoSave(900000000000000LL);
353      OrbitalInfotrclone->SetBranchAddress("OrbitalInfo",&orbitalinfoclone);      OrbitalInfotrclone->SetBranchAddress("OrbitalInfo",&orbitalinfoclone);
354      //            //      
355      if ( nobefrun > 0 ){      if ( nobefrun > 0 ){
356        if (verbose){        if (verbose){
357        printf("\n Pre-processing: copying events from the old tree before the processed run\n");            printf("\n Pre-processing: copying events from the old tree before the processed run\n");  
358        printf(" Copying %u events in the file which are before the beginning of the run %llu \n",nobefrun,run);          printf(" Copying %u events in the file which are before the beginning of the run %u \n",nobefrun,run);
359        printf(" Start copying at event number 0, end copying at event number %u \n",nobefrun);          printf(" Start copying at event number 0, end copying at event number %u \n",nobefrun);
360        }        }
361        for (UInt_t j = 0; j < nobefrun; j++){        for (UInt_t j = 0; j < nobefrun; j++){
362          //          //
# Line 299  int OrbitalInfoCore(ULong64_t run, TFile Line 364  int OrbitalInfoCore(ULong64_t run, TFile
364          //          //
365          // copy orbitalinfoclone to mydec          // copy orbitalinfoclone to mydec
366          //          //
         //      orbitalinfo = new OrbitalInfo();  
367          orbitalinfo->Clear();          orbitalinfo->Clear();
368            //
369          memcpy(&orbitalinfo,&orbitalinfoclone,sizeof(orbitalinfoclone));          memcpy(&orbitalinfo,&orbitalinfoclone,sizeof(orbitalinfoclone));
370          //          //
371          // Fill entry in the new tree          // Fill entry in the new tree
# Line 329  int OrbitalInfoCore(ULong64_t run, TFile Line 394  int OrbitalInfoCore(ULong64_t run, TFile
394        printf(" ####################################################################### \n\n\n");        printf(" ####################################################################### \n\n\n");
395      }      }
396      //      //
397      runinfo->ID_REG_RUN = 0ULL;      runinfo->ID_ROOT_L0 = 0;
398      //      //
399      // store in the runinfo class the GL_RUN variables for our run      // store in the runinfo class the GL_RUN variables for our run
400      //      //
401      sgnl = 0;      sgnl = 0;
402      sgnl = runinfo->GetRunInfo(idRun);      sgnl = runinfo->GetRunInfo(idRun);
403      if ( sgnl ){      if ( sgnl ){
404        //printf("\n OrbitalInfo - ERROR: RunInfo exited with non-zero status\n");        if ( debug ) printf("\n OrbitalInfo - ERROR: RunInfo exited with non-zero status\n");
405        code = sgnl;        code = sgnl;
406        goto closeandexit;        goto closeandexit;
407      } else {      } else {
# Line 345  int OrbitalInfoCore(ULong64_t run, TFile Line 410  int OrbitalInfoCore(ULong64_t run, TFile
410      //      //
411      // now you can access that variables using the RunInfo class this way runinfo->ID_REG_RUN      // now you can access that variables using the RunInfo class this way runinfo->ID_REG_RUN
412      //      //
413      if ( runinfo->ID_REG_RUN == 0 ){      if ( runinfo->ID_ROOT_L0 == 0 ){
414        //printf("\n OrbitalInfo - ERROR: no run with ID_RUN = %i \n\n Exiting... \n\n",(int)idRun);        if ( debug ) printf("\n OrbitalInfo - ERROR: no run with ID_RUN = %u \n\n Exiting... \n\n",idRun);
415        code = -5;        code = -5;
416        goto closeandexit;            goto closeandexit;    
417      };      };
418      //      //
419        // prepare the timesync for the db
420        //
421        dbtime = new GL_TIMESYNC(runinfo->ID_ROOT_L0,"ID",dbc);
422        //
423      // Search in the DB the path and name of the LEVEL0 file to be processed.      // Search in the DB the path and name of the LEVEL0 file to be processed.
424      //      //
425      glroot->Query_GL_ROOT(runinfo->ID_REG_RUN,dbc);      glroot->Query_GL_ROOT(runinfo->ID_ROOT_L0,dbc);
426      //      //
427      ftmpname.str("");      ftmpname.str("");
428      ftmpname << glroot->PATH.Data() << "/";      ftmpname << glroot->PATH.Data() << "/";
# Line 362  int OrbitalInfoCore(ULong64_t run, TFile Line 431  int OrbitalInfoCore(ULong64_t run, TFile
431      //      //
432      // print out informations      // print out informations
433      //      //
434      totevent = runinfo->EV_REG_PHYS_TO - runinfo->EV_REG_PHYS_FROM + 1;      totevent = runinfo->NEVENTS;
435      if (verbose){      if (verbose){
436        printf("\n LEVEL0 data file: %s \n",fname.Data());        printf("\n LEVEL0 data file: %s \n",fname.Data());
437        printf(" RUN HEADER absolute time is:  %llu \n",runinfo->RUNHEADER_TIME);        printf(" RUN HEADER absolute time is:  %u \n",runinfo->RUNHEADER_TIME);
438        printf(" RUN TRAILER absolute time is: %llu \n",runinfo->RUNTRAILER_TIME);        printf(" RUN TRAILER absolute time is: %u \n",runinfo->RUNTRAILER_TIME);
439        printf(" %i events to be processed for run %llu: from %i to %i (reg entries)\n\n",totevent,idRun,runinfo->EV_REG_PHYS_FROM,runinfo->EV_REG_PHYS_TO);        printf(" %i events to be processed for run %u: from %i to %i \n\n",totevent,idRun,runinfo->EV_FROM,runinfo->EV_FROM+totevent);
440      }//      }//
441      // Open Level0 file      // Open Level0 file
442      l0File = new TFile(fname.Data());      l0File = new TFile(fname.Data());
443      if ( !l0File ) {      if ( !l0File ) {
444        //printf(" OrbitalInfo - ERROR: problems opening Level0 file\n");        if ( debug ) printf(" OrbitalInfo - ERROR: problems opening Level0 file\n");
445        code = -6;        code = -6;
446        goto closeandexit;        goto closeandexit;
447      };      };
448      l0tr = (TTree*)l0File->Get("Physics");      l0tr = (TTree*)l0File->Get("Physics");
449      if ( !l0tr ) {      if ( !l0tr ) {
450        //printf(" OrbitalInfo - ERROR: no Physics tree in Level0 file\n");        if ( debug ) printf(" OrbitalInfo - ERROR: no Physics tree in Level0 file\n");
451        l0File->Close();        l0File->Close();
452        code = -7;        code = -7;
453        goto closeandexit;        goto closeandexit;
# Line 386  int OrbitalInfoCore(ULong64_t run, TFile Line 455  int OrbitalInfoCore(ULong64_t run, TFile
455      // EM: open header branch as well      // EM: open header branch as well
456      l0head = l0tr->GetBranch("Header");      l0head = l0tr->GetBranch("Header");
457      if ( !l0head ) {      if ( !l0head ) {
458        //if ( verbose ) printf(" CALORIMETER - ERROR: no Header branch in Level0 tree\n");        if ( debug ) printf(" OrbitalInfo - ERROR: no Header branch in Level0 tree\n");
459        l0File->Close();        l0File->Close();
460        code = -8;        code = -8;
461        goto closeandexit;            goto closeandexit;    
462      };      };
463      l0tr->SetBranchAddress("Header", &eh);      l0tr->SetBranchAddress("Header", &eh);
464      // end EM      // end EM
465      l0registry = l0tr->GetBranch("Registry");      nevents = l0head->GetEntries();
     if ( !l0registry ) {  
       //printf(" OrbitalInfo - ERROR: no Registry branch in Level0 tree\n");  
       l0File->Close();  
       code = -9;  
       goto closeandexit;      
     };  
     //  
     l0tr->SetBranchAddress("Registry", &l0reg);  
     //  
     nevents = l0registry->GetEntries();  
466      //      //
467      if ( nevents < 1 ) {      if ( nevents < 1 ) {
468        //printf(" OrbitalInfo - ERROR: Level0 file is empty\n\n");        if ( debug ) printf(" OrbitalInfo - ERROR: Level0 file is empty\n\n");
469        l0File->Close();        l0File->Close();
470        code = -11;        code = -11;
471        goto closeandexit;        goto closeandexit;
472      };      };
473      //      //
474      if ( runinfo->EV_REG_PHYS_TO > nevents-1 ) {      if ( runinfo->EV_TO > nevents-1 ) {
475        //printf(" OrbitalInfo - ERROR: too few entries in the registry tree\n");        if ( debug ) printf(" OrbitalInfo - ERROR: too few entries in the registry tree\n");
476        l0File->Close();        l0File->Close();
477        code = -12;        code = -12;
478        goto closeandexit;        goto closeandexit;
# Line 423  int OrbitalInfoCore(ULong64_t run, TFile Line 482  int OrbitalInfoCore(ULong64_t run, TFile
482      //      //
483      if (verbose) printf("\n Ready to start! \n\n Processed events: \n\n");      if (verbose) printf("\n Ready to start! \n\n Processed events: \n\n");
484      //      //
485      for ( re = runinfo->EV_REG_PHYS_FROM; re <= runinfo->EV_REG_PHYS_TO; re++){      for ( re = runinfo->EV_FROM; re < (runinfo->EV_FROM+runinfo->NEVENTS); re++){
486        //        //
487        if ( procev%1000 == 0 && procev > 0 && verbose) printf(" %iK \n",procev/1000);            if ( procev%1000 == 0 && procev > 0 && verbose ) printf(" %iK \n",procev/1000);  
488        //        //
489        l0registry->GetEntry(re);        l0head->GetEntry(re);
490        //        //
491        // absolute time of this event        // absolute time of this event
492        //        //
493        atime = l0reg->absTime;        ph = eh->GetPscuHeader();
494        //        atime = dbtime->DBabsTime(ph->GetOrbitalTime());  
       // physics events is at entry number ei where  
       //  
       ei = l0reg->event;  
495        //        //
496        // paranoid check        // paranoid check
497        //        //
498        if ( (atime > runinfo->RUNTRAILER_TIME) || (atime < runinfo->RUNHEADER_TIME)  ) {        if ( (atime > runinfo->RUNTRAILER_TIME) || (atime < runinfo->RUNHEADER_TIME)  ) {
499          if (verbose) printf(" OrbitalInfo - WARNING: event at time outside the run time window, skipping it\n");          if (verbose) printf(" OrbitalInfo - WARNING: event at time outside the run time window, skipping it\n");
500          goto jumpev;          debug = true;
501        };          continue;
502          }
503        //        //
504        procev++;        procev++;
505        //        //
506        // start processing        // start processing
507        //        //
508        orbitalinfo->Clear();        orbitalinfo->Clear();
509        //orbitalinfo = new OrbitalInfo();        //
510        orbitalinfo->absTime = l0reg->absTime;        // CHANGE HERE!!!!
511          //
512          orbitalinfo->absTime = atime;
513        // EM: add OBT and plt_num infos from the header        // EM: add OBT and plt_num infos from the header
       l0head->GetEntry(ei);  
514        ph = eh->GetPscuHeader();        ph = eh->GetPscuHeader();
515        orbitalinfo->pkt_num = ph->GetCounter();        orbitalinfo->pkt_num = ph->GetCounter();
516        orbitalinfo->OBT = ph->GetOrbitalTime();        orbitalinfo->OBT = ph->GetOrbitalTime();
517    
518          // If the absolute time of the event overpass the time of the
519          // tle, get a new tle.  GL_TLE::GetToTime() default to zero.
520          //
521          // I also use this condition to compute the dipole moment dimo.
522          // It's really redundant to compute it so often because
523          // probably it will not change at all.  But the overhead is
524          // minimum.
525          float jyear=0;
526    
527          if(atime >= gltle->GetToTime()) {
528            if ( !gltle->Query(atime, dbc) ){
529              
530              // Compute the magnetic dipole moment.
531              UInt_t year, month, day, hour, min, sec;
532              
533              TTimeStamp t = TTimeStamp(atime, kTRUE);
534              t.GetDate(kTRUE, 0, &year, &month, &day);
535              t.GetTime(kTRUE, 0, &hour, &min, &sec);
536              jyear = (float) year
537                + (month*31.+ (float) day)/365.
538                + (hour*3600.+min*60.+(float)sec)/(24*3600*365.);
539              
540              feldcof_(&jyear, &dimo); // get dipole moment for year
541            } else {
542              code = -56;
543              goto closeandexit;
544            };
545          }
546    
547          // Propagate the orbit from the tle time to atime, using SGP(D)4.
548          cCoordGeo coo = getCoo(atime, gltle->GetFromTime(), gltle->GetTle());
549    
550          // Build coordinates in the right range.  We want to convert,
551          // longitude from (0, 2*pi) to (-180deg, 180deg).  Altitude is
552          // in meters.
553          lon = (coo.m_Lon > M_PI) ? rad2deg(coo.m_Lon - 2*M_PI) : rad2deg(coo.m_Lon);
554          lat = rad2deg(coo.m_Lat);
555          alt = coo.m_Alt;
556    
557    
558          //      if((lon>180) || (lon<-180) || (lat>90) || (lat<-90) || (alt<0))
559          //        continue;
560          if( lon<180 && lon>-180 && lat<90 && lat>-90 && alt>0 ){  
561    
562            orbitalinfo->lon = lon;
563            orbitalinfo->lat = lat;
564            orbitalinfo->alt = alt ;
565            
566            // compute mag field components and L shell.
567            feldg_(&lat, &lon, &alt, &bnorth, &beast, &bdown, &babs);
568            shellg_(&lat, &lon, &alt, &dimo, &xl, &icode, &bab1);
569            findb0_(&stps, &bdel, &value, &bequ, &rr0);
570            
571            orbitalinfo->Bnorth = bnorth;
572            orbitalinfo->Beast = beast;
573            orbitalinfo->Bdown = bdown;
574            orbitalinfo->Babs = babs;
575            orbitalinfo->BB0 = babs/bequ;
576            orbitalinfo->L = xl;
577            
578            // Set Stormer vertical cutoff using L shell.
579            orbitalinfo->cutoff[0] = 14.9/(xl*xl);
580          };
581          
582        // end EM        // end EM
583        OrbitalInfotr->Fill();        OrbitalInfotr->Fill();
584        //          }
585        //      
     jumpev:  
       debug = false;  
       //  
     };  
586      //      //
587      // Here you may want to clear some variables before processing another run        // Here you may want to clear some variables before processing another run  
588      //      //
589      ei = 0;      delete dbtime;
590    }; // process all the runs    }; // process all the runs
591    //    //
592    if (verbose) printf("\n Finished processing data \n");    if (verbose) printf("\n Finished processing data \n");
# Line 492  int OrbitalInfoCore(ULong64_t run, TFile Line 611  int OrbitalInfoCore(ULong64_t run, TFile
611          // copy orbitalinfoclone to OrbitalInfo          // copy orbitalinfoclone to OrbitalInfo
612          //          //
613          orbitalinfo->Clear();          orbitalinfo->Clear();
614          //      orbitalinfo = new OrbitalInfo();          //
615          memcpy(&orbitalinfo,&orbitalinfoclone,sizeof(orbitalinfoclone));          memcpy(&orbitalinfo,&orbitalinfoclone,sizeof(orbitalinfoclone));
616          //          //
617          // Fill entry in the new tree          // Fill entry in the new tree
# Line 508  int OrbitalInfoCore(ULong64_t run, TFile Line 627  int OrbitalInfoCore(ULong64_t run, TFile
627    if ( l0File ) l0File->Close();    if ( l0File ) l0File->Close();
628    if ( tempfile ) tempfile->Close();                if ( tempfile ) tempfile->Close();            
629    gSystem->Unlink(tempname.str().c_str());    gSystem->Unlink(tempname.str().c_str());
630      //
   //if ( code < 0 ) printf("\n OrbitalInfo - ERROR: an error occurred, try to save anyway...\n");  
   //printf("\n Writing and closing rootple\n");  
631    if ( runinfo ) runinfo->Close();        if ( runinfo ) runinfo->Close();    
632    if ( OrbitalInfotr ) OrbitalInfotr->SetName("OrbitalInfo");        if ( OrbitalInfotr ) OrbitalInfotr->SetName("OrbitalInfo");    
633    if ( file ){    if ( file ){
# Line 534  int OrbitalInfoCore(ULong64_t run, TFile Line 651  int OrbitalInfoCore(ULong64_t run, TFile
651    return(code);    return(code);
652  }  }
653    
654    
655    //
656    // Returns the cCoordGeo structure holding the geographical
657    // coordinates for the event (see sgp4.h).
658    //
659    // atime is the abstime of the event in UTC unix time.
660    // tletime is the time of the tle in UTC unix time.
661    // tle is the previous and nearest tle (compared to atime).
662    cCoordGeo getCoo(UInt_t atime, UInt_t tletime, cTle *tle)
663    {
664      cEci eci;
665      cOrbit orbit(*tle);
666      
667      orbit.getPosition((double) (atime - tletime)/60., &eci);
668      
669      return eci.toGeo();
670    }
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.10
  ViewVC Help
Powered by ViewVC 1.1.23