/[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.10 by mocchiut, Wed Jan 24 13:16:26 2007 UTC revision 1.30 by mocchiut, Wed Oct 1 15:25:44 2008 UTC
# Line 34  Line 34 
34  #include <PscuHeader.h>  #include <PscuHeader.h>
35  #include <PscuEvent.h>  #include <PscuEvent.h>
36  #include <EventHeader.h>  #include <EventHeader.h>
37    #include <mcmd/McmdEvent.h>
38    #include <mcmd/McmdRecord.h>
39  //  //
40  // This program headers  // This program headers
41  //  //
42  #include <OrbitalInfo.h>  #include <OrbitalInfo.h>
43  #include <OrbitalInfoVerl2.h>  #include <OrbitalInfoVerl2.h>
44  #include <OrbitalInfoCore.h>  #include <OrbitalInfoCore.h>
45    #include <InclinationInfo.h>
46    
47  using namespace std;  using namespace std;
48    
# Line 47  using namespace std; Line 50  using namespace std;
50  // CORE ROUTINE  // CORE ROUTINE
51  //  //
52  //  //
53  int OrbitalInfoCore(UInt_t run, TFile *file, TSQLServer *dbc, Int_t OrbitalInfoargc, char *OrbitalInfoargv[]){  int OrbitalInfoCore(UInt_t run, TFile *file, GL_TABLES *glt, Int_t OrbitalInfoargc, char *OrbitalInfoargv[]){
54      //
 //   // Temporary check to use igrf routines.  We need dat files in the  
 //   // current directory.  
 //   fstream igrfdat1("igrf05.dat");  
 //   fstream igrfdat2("igrf05s.dat");  
 //   if( (!igrfdat1) && (!igrfdat2)) {  
 //     cerr<<"\n**************************************\n"  
 //      <<"igrf05.dat or igrf05s.dat not in the current directory.  Exiting.\n"  
 //      <<"**************************************\n";  
 //     exit(EXIT_FAILURE);  
 //   }  
 //   igrfdat1.close();  
 //   igrfdat2.close();  
 //   // end of temporary code  
 //   //  
   
55    Int_t i = 0;    Int_t i = 0;
56      TString host = glt->CGetHost();
57      TString user = glt->CGetUser();
58      TString psw = glt->CGetPsw();
59      TSQLServer *dbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data());
60      //
61      stringstream myquery;
62      myquery.str("");
63      myquery << "SET time_zone='+0:00'";
64      dbc->Query(myquery.str().c_str());
65    //    //
66    TString processFolder = "OrbitalInfoFolder";    TString processFolder = Form("OrbitalInfoFolder_%u",run);
67    //    //
68    // Set these to true to have a very verbose output.    // Set these to true to have a very verbose output.
69    //    //
70    Bool_t debug = false;    Bool_t debug = false;
71    //    //
72    Bool_t verbose = false;    Bool_t verbose = false;
73      //
74      Bool_t standalone = true;
75      //
76    if ( OrbitalInfoargc > 0 ){    if ( OrbitalInfoargc > 0 ){
77      i = 0;      i = 0;
78      while ( i < OrbitalInfoargc ){      while ( i < OrbitalInfoargc ){
# Line 86  int OrbitalInfoCore(UInt_t run, TFile *f Line 85  int OrbitalInfoCore(UInt_t run, TFile *f
85        };        };
86        if ( (!strcmp(OrbitalInfoargv[i],"--debug")) || (!strcmp(OrbitalInfoargv[i],"-g")) ) {        if ( (!strcmp(OrbitalInfoargv[i],"--debug")) || (!strcmp(OrbitalInfoargv[i],"-g")) ) {
87          verbose = true;          verbose = true;
88            debug = true;
89        };        };
90        if ( (!strcmp(OrbitalInfoargv[i],"--verbose")) || (!strcmp(OrbitalInfoargv[i],"-v")) ) {        if ( (!strcmp(OrbitalInfoargv[i],"--verbose")) || (!strcmp(OrbitalInfoargv[i],"-v")) ) {
91          verbose = true;          verbose = true;
92        };        };
93          if ( (!strcmp(OrbitalInfoargv[i],"--standalone")) ) {
94            standalone = true;
95          };
96          if ( (!strcmp(OrbitalInfoargv[i],"--calculate-pitch")) ) {
97            standalone = false;
98          };
99        i++;        i++;
100      };      };
101    };    };
# Line 98  int OrbitalInfoCore(UInt_t run, TFile *f Line 104  int OrbitalInfoCore(UInt_t run, TFile *f
104    //    //
105    TTree *OrbitalInfotr = 0;    TTree *OrbitalInfotr = 0;
106    UInt_t nevents = 0;    UInt_t nevents = 0;
107      UInt_t neventsm = 0;
108    //    //
109    // variables needed to reprocess data    // variables needed to reprocess data
110    //    //
# Line 116  int OrbitalInfoCore(UInt_t run, TFile *f Line 123  int OrbitalInfoCore(UInt_t run, TFile *f
123    stringstream ftmpname;    stringstream ftmpname;
124    TString fname;    TString fname;
125    UInt_t totfileentries = 0;    UInt_t totfileentries = 0;
126    UInt_t idRun = 0;    UInt_t idRun = 0;
127      //
128      // My variables. Vitaly.
129      //
130      //  UInt_t iev = 0;
131      //  UInt_t j3 = 0;
132      UInt_t oi = 0;
133      Int_t tmpSize = 0;
134    //    //
135    // variables needed to handle error signals    // variables needed to handle error signals
136    //    //
# Line 132  int OrbitalInfoCore(UInt_t run, TFile *f Line 146  int OrbitalInfoCore(UInt_t run, TFile *f
146    //    //
147    TFile *l0File = 0;    TFile *l0File = 0;
148    TTree *l0tr = 0;    TTree *l0tr = 0;
149      TTree *l0trm = 0;
150    // EM: open also header branch    // EM: open also header branch
151    TBranch *l0head = 0;    TBranch *l0head = 0;
152    pamela::EventHeader *eh = 0;    pamela::EventHeader *eh = 0;
153    pamela::PscuHeader *ph = 0;    pamela::PscuHeader *ph = 0;
154      pamela::McmdEvent *mcmdev = 0;
155      pamela::McmdRecord *mcmdrc = 0;
156    // end EM    // end EM
157      
158      //  pamela::RunHeaderEvent *reh = new pamela::RunHeaderEvent;
159      //  pamela::EventHeader    *eH  = new pamela::EventHeader;
160      
161    //    //
162    // Define other basic variables    // Define other basic variables
163    //    //
# Line 147  int OrbitalInfoCore(UInt_t run, TFile *f Line 168  int OrbitalInfoCore(UInt_t run, TFile *f
168    Int_t totevent = 0;    Int_t totevent = 0;
169    UInt_t atime = 0;    UInt_t atime = 0;
170    UInt_t re = 0;    UInt_t re = 0;
171      UInt_t ik = 0;
172    
173    // Position    // Position
174    Float_t lon, lat, alt;    Float_t lon, lat, alt;
# Line 179  int OrbitalInfoCore(UInt_t run, TFile *f Line 201  int OrbitalInfoCore(UInt_t run, TFile *f
201    TTree *tempOrbitalInfo = 0;    TTree *tempOrbitalInfo = 0;
202    stringstream tempname;    stringstream tempname;
203    stringstream OrbitalInfofolder;    stringstream OrbitalInfofolder;
204      Bool_t myfold = false;
205    tempname.str("");    tempname.str("");
206    tempname << outDir;    tempname << outDir;
207    tempname << "/" << processFolder.Data();    tempname << "/" << processFolder.Data();
208    OrbitalInfofolder.str("");    OrbitalInfofolder.str("");
209    OrbitalInfofolder << tempname.str().c_str();    OrbitalInfofolder << tempname.str().c_str();
   gSystem->MakeDirectory(OrbitalInfofolder.str().c_str());  
210    tempname << "/OrbitalInfotree_run";    tempname << "/OrbitalInfotree_run";
211    tempname << run << ".root";      tempname << run << ".root";  
212    //    //
# Line 193  int OrbitalInfoCore(UInt_t run, TFile *f Line 215  int OrbitalInfoCore(UInt_t run, TFile *f
215    GL_ROOT *glroot = new GL_ROOT();    GL_ROOT *glroot = new GL_ROOT();
216    GL_TIMESYNC *dbtime = 0;    GL_TIMESYNC *dbtime = 0;
217    GL_TLE *gltle = new GL_TLE();    GL_TLE *gltle = new GL_TLE();
218      //
219      //Quaternions classes
220      //
221      Quaternions *L_QQ_Q_l_lower = new Quaternions();
222      InclinationInfo *RYPang_lower = new InclinationInfo();
223      Quaternions *L_QQ_Q_l_upper = new Quaternions();
224      InclinationInfo *RYPang_upper = new InclinationInfo();
225      
226      cEci eCi;
227      
228    // Initialize fortran routines!!!    // Initialize fortran routines!!!
229    Int_t ltp2 = 0;    Int_t ltp2 = 0;
230    Int_t ltp3 = 0;    Int_t ltp3 = 0;
# Line 201  int OrbitalInfoCore(UInt_t run, TFile *f Line 233  int OrbitalInfoCore(UInt_t run, TFile *f
233    GL_PARAM *glparam = new GL_PARAM();    GL_PARAM *glparam = new GL_PARAM();
234    GL_PARAM *glparam2 = new GL_PARAM();    GL_PARAM *glparam2 = new GL_PARAM();
235    Int_t parerror=glparam->Query_GL_PARAM(1,301,dbc); // parameters stored in DB in GL_PRAM table    Int_t parerror=glparam->Query_GL_PARAM(1,301,dbc); // parameters stored in DB in GL_PRAM table
236      //
237      // Orientation variables
238      //
239      UInt_t evfrom = 0;
240      UInt_t jumped = 0;
241      Int_t itr = -1;    
242      Double_t A1;
243      Double_t A2;
244      Double_t A3;
245      Double_t Px = 0;
246      Double_t Py = 0;      
247      Double_t Pz = 0;  
248      TTree *ttof = 0;
249      ToFLevel2 *tof = new ToFLevel2();
250      OrientationInfo *PO = new OrientationInfo();
251      Int_t nz = 6;
252      Float_t zin[6];
253      Int_t nevtofl2 = 0;
254      //  
255    if ( parerror<0 ) {    if ( parerror<0 ) {
256      code = parerror;      code = parerror;
257      goto closeandexit;      goto closeandexit;
# Line 217  int OrbitalInfoCore(UInt_t run, TFile *f Line 268  int OrbitalInfoCore(UInt_t run, TFile *f
268    if ( verbose ) printf(" Reading Earth's Magnetic Field parameter file: %s \n",(glparam2->PATH+glparam2->NAME).Data());    if ( verbose ) printf(" Reading Earth's Magnetic Field parameter file: %s \n",(glparam2->PATH+glparam2->NAME).Data());
269    //    //
270    initize_((char *)niente,&uno,(char *)(glparam->PATH+glparam->NAME).Data(),&ltp2,(char *)(glparam2->PATH+glparam2->NAME).Data(),&ltp3);    initize_((char *)niente,&uno,(char *)(glparam->PATH+glparam->NAME).Data(),&ltp2,(char *)(glparam2->PATH+glparam2->NAME).Data(),&ltp3);
   
271    //    //
272    // End IGRF stuff//    // End IGRF stuff//
273    //    //
274      for (Int_t ip=0;ip<nz;ip++){
275        zin[ip] = tof->GetZTOF(tof->GetToFPlaneID(ip));
276      };
277      //
278      if ( !standalone ){
279        //
280        // Does it contain the Tracker tree?
281        //
282        ttof = (TTree*)file->Get("ToF");
283        if ( !ttof ) {
284          if ( verbose ) printf(" OrbitalInfo - ERROR: no tof tree\n");
285          code = -900;
286          goto closeandexit;
287        };
288        ttof->SetBranchAddress("ToFLevel2",&tof);  
289        nevtofl2 = ttof->GetEntries();
290      };
291    //    //
292    // Let's start!    // Let's start!
293    //    //
# Line 276  int OrbitalInfoCore(UInt_t run, TFile *f Line 342  int OrbitalInfoCore(UInt_t run, TFile *f
342    // number of run to be processed    // number of run to be processed
343    //    //
344    numbofrun = runinfo->GetNoRun();    numbofrun = runinfo->GetNoRun();
345      UInt_t totnorun = runinfo->GetRunEntries();
346    //    //
347    // Try to access the OrbitalInfo tree in the file, if it exists we are reprocessing data if not we are processing a new run    // Try to access the OrbitalInfo tree in the file, if it exists we are reprocessing data if not we are processing a new run
348    //    //
# Line 302  int OrbitalInfoCore(UInt_t run, TFile *f Line 369  int OrbitalInfoCore(UInt_t run, TFile *f
369      //      //
370      if (verbose) printf("\n Preparing the pre-processing...\n");      if (verbose) printf("\n Preparing the pre-processing...\n");
371      //      //
372      if ( run == 0 ){      if ( run == 0 || totnorun == 1 ){
373        //        //
374        // we are reprocessing all the file        // we are reprocessing all the file
375        // 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 321  int OrbitalInfoCore(UInt_t run, TFile *f Line 388  int OrbitalInfoCore(UInt_t run, TFile *f
388        //        //
389        // copying old tree to a new file        // copying old tree to a new file
390        //        //
391          gSystem->MakeDirectory(OrbitalInfofolder.str().c_str());
392          myfold = true;
393        tempfile = new TFile(tempname.str().c_str(),"RECREATE");        tempfile = new TFile(tempname.str().c_str(),"RECREATE");
394        tempOrbitalInfo = OrbitalInfotrclone->CloneTree(-1,"fast");        tempOrbitalInfo = OrbitalInfotrclone->CloneTree(-1,"fast");
395        tempOrbitalInfo->SetName("OrbitalInfo-old");        tempOrbitalInfo->SetName("OrbitalInfo-old");
# Line 341  int OrbitalInfoCore(UInt_t run, TFile *f Line 410  int OrbitalInfoCore(UInt_t run, TFile *f
410    file->cd();    file->cd();
411    OrbitalInfotr = new TTree("OrbitalInfo-new","PAMELA OrbitalInfo data");    OrbitalInfotr = new TTree("OrbitalInfo-new","PAMELA OrbitalInfo data");
412    OrbitalInfotr->SetAutoSave(900000000000000LL);    OrbitalInfotr->SetAutoSave(900000000000000LL);
413      orbitalinfo->Set();//ELENA **TEMPORANEO?**
414    OrbitalInfotr->Branch("OrbitalInfo","OrbitalInfo",&orbitalinfo);    OrbitalInfotr->Branch("OrbitalInfo","OrbitalInfo",&orbitalinfo);
415    //    //
416    if ( reproc && !reprocall ){    if ( reproc && !reprocall ){
# Line 387  int OrbitalInfoCore(UInt_t run, TFile *f Line 457  int OrbitalInfoCore(UInt_t run, TFile *f
457      //      //
458      // retrieve the first run ID to be processed using the RunInfo list      // retrieve the first run ID to be processed using the RunInfo list
459      //      //
460        
461      idRun = runlist->At(irun);      idRun = runlist->At(irun);
462      if (verbose){      if (verbose){
463        printf("\n\n\n ####################################################################### \n");        printf("\n\n\n ####################################################################### \n");
# Line 419  int OrbitalInfoCore(UInt_t run, TFile *f Line 490  int OrbitalInfoCore(UInt_t run, TFile *f
490      // prepare the timesync for the db      // prepare the timesync for the db
491      //      //
492      dbtime = new GL_TIMESYNC(runinfo->ID_ROOT_L0,"ID",dbc);      dbtime = new GL_TIMESYNC(runinfo->ID_ROOT_L0,"ID",dbc);
493      
494      //      //
495      // 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.
496      //      //
# Line 428  int OrbitalInfoCore(UInt_t run, TFile *f Line 500  int OrbitalInfoCore(UInt_t run, TFile *f
500      ftmpname << glroot->PATH.Data() << "/";      ftmpname << glroot->PATH.Data() << "/";
501      ftmpname << glroot->NAME.Data();      ftmpname << glroot->NAME.Data();
502      fname = ftmpname.str().c_str();      fname = ftmpname.str().c_str();
503        ftmpname.str("");
504      //      //
505      // print out informations      // print out informations
506      //      //
507      totevent = runinfo->NEVENTS;      totevent = runinfo->NEVENTS;
508        evfrom = runinfo->EV_FROM;
509        //cout<<"totevents = "<<totevent<<"\n";
510      if (verbose){      if (verbose){
511        printf("\n LEVEL0 data file: %s \n",fname.Data());        printf("\n LEVEL0 data file: %s \n",fname.Data());
512        printf(" RUN HEADER absolute time is:  %u \n",runinfo->RUNHEADER_TIME);        printf(" RUN HEADER absolute time is:  %u \n",runinfo->RUNHEADER_TIME);
513        printf(" RUN TRAILER absolute time is: %u \n",runinfo->RUNTRAILER_TIME);        printf(" RUN TRAILER absolute time is: %u \n",runinfo->RUNTRAILER_TIME);
514        printf(" %i events to be processed for run %u: from %i to %i \n\n",totevent,idRun,runinfo->EV_FROM,runinfo->EV_FROM+totevent);        printf(" %i events to be processed for run %u: from %i to %i \n\n",totevent,idRun,runinfo->EV_FROM+1,runinfo->EV_FROM+totevent);
515      }//      }//
516        //
517        //    if ( !totevent ) goto closeandexit;
518      // Open Level0 file      // Open Level0 file
519      l0File = new TFile(fname.Data());      l0File = new TFile(fname.Data());
520      if ( !l0File ) {      if ( !l0File ) {
# Line 464  int OrbitalInfoCore(UInt_t run, TFile *f Line 541  int OrbitalInfoCore(UInt_t run, TFile *f
541      // end EM      // end EM
542      nevents = l0head->GetEntries();      nevents = l0head->GetEntries();
543      //      //
544      if ( nevents < 1 ) {      if ( nevents < 1 && totevent ) {
545        if ( debug ) printf(" OrbitalInfo - ERROR: Level0 file is empty\n\n");        if ( debug ) printf(" OrbitalInfo - ERROR: Level0 file is empty\n\n");
546        l0File->Close();        l0File->Close();
547        code = -11;        code = -11;
548        goto closeandexit;        goto closeandexit;
549      };      };
550      //      //
551      if ( runinfo->EV_TO > nevents-1 ) {      if ( runinfo->EV_TO > nevents-1 && totevent ) {
552        if ( debug ) printf(" OrbitalInfo - ERROR: too few entries in the registry tree\n");        if ( debug ) printf(" OrbitalInfo - ERROR: too few entries in the registry tree\n");
553        l0File->Close();        l0File->Close();
554        code = -12;        code = -12;
555        goto closeandexit;        goto closeandexit;
556      };      };
557      //      //
558    //     TTree *tp = (TTree*)l0File->Get("RunHeader");
559    //     tp->SetBranchAddress("Header", &eH);
560    //     tp->SetBranchAddress("RunHeader", &reh);
561    //     tp->GetEntry(0);
562    //     ph = eH->GetPscuHeader();
563    //     ULong_t TimeSync = reh->LAST_TIME_SYNC_INFO;
564    //     ULong_t ObtSync = reh->OBT_TIME_SYNC;    
565    //     if ( debug ) printf(" 1 TimeSync %lu ObtSync %lu DeltaOBT %lu\n",TimeSync,ObtSync,TimeSync-ObtSync);
566    //
567        ULong_t TimeSync = (ULong_t)dbtime->GetTimesync();
568        ULong_t ObtSync = (ULong_t)(dbtime->GetObt0()/1000);
569        ULong_t DeltaOBT = TimeSync - ObtSync;
570    
571        if ( debug ) printf(" 2 TimeSync %lu ObtSync %lu DeltaOBT %lu\n",(ULong_t)(dbtime->GetTimesync()/1000),(ULong_t)dbtime->GetObt0(),TimeSync-ObtSync);
572        
573        l0trm = (TTree*)l0File->Get("Mcmd");
574        neventsm = l0trm->GetEntries();
575        //    neventsm = 0;
576        //
577        if (neventsm == 0){
578          if ( debug ) printf("InclinationInfo - WARNING: No quaternions in this File");
579          //      l0File->Close();
580          code = 900;
581          //      goto closeandexit;
582        }
583        //
584        
585        l0trm->SetBranchAddress("Mcmd", &mcmdev);
586        //    l0trm->SetBranchAddress("Header", &eh);
587        //
588        //
589        //
590        UInt_t mctren = 0;    
591        UInt_t mcreen = 0;  
592        UInt_t numrec = 0;
593        //
594        Double_t upperqtime = 0;
595        Double_t lowerqtime = 0;
596        
597        Double_t incli = 0;
598        oi = 0;
599        UInt_t ooi = 0;
600        //
601        // init quaternions sync
602        //
603        Bool_t isf = true;
604        Int_t fgh = 0;
605        //
606      // run over all the events of the run      // run over all the events of the run
607      //      //
608      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");
609      //      //
610         //
611      for ( re = runinfo->EV_FROM; re < (runinfo->EV_FROM+runinfo->NEVENTS); re++){      for ( re = runinfo->EV_FROM; re < (runinfo->EV_FROM+runinfo->NEVENTS); re++){
612        
613        //        //
614        if ( procev%1000 == 0 && procev > 0 && verbose ) printf(" %iK \n",procev/1000);          if ( procev%1000 == 0 && procev > 0 && verbose ) printf(" %iK \n",procev/1000);  
615          if ( debug ) printf(" %i \n",procev);      
616        //        //
617        l0head->GetEntry(re);        l0head->GetEntry(re);
618        //        //
619        // absolute time of this event        // absolute time of this event
620        //        //
621        ph = eh->GetPscuHeader();        ph = eh->GetPscuHeader();
622        atime = dbtime->DBabsTime(ph->GetOrbitalTime());          atime = dbtime->DBabsTime(ph->GetOrbitalTime());
623        //        //
624        // paranoid check        // paranoid check
625        //        //
626        if ( (atime > runinfo->RUNTRAILER_TIME) || (atime < runinfo->RUNHEADER_TIME)  ) {        if ( (atime > runinfo->RUNTRAILER_TIME) || (atime < runinfo->RUNHEADER_TIME)  ) {
627          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");
628          debug = true;          jumped++;
629    //      debug = true;
630          continue;          continue;
631        }        }
632    
633          //
634          // retrieve tof informations
635          //
636          if ( !reprocall ){
637            itr = nobefrun + (re - evfrom - jumped);
638            //itr = re-(46438+200241);
639          } else {
640            itr = runinfo->GetFirstEntry() + (re - evfrom - jumped);
641          };
642          //
643          if ( !standalone ){
644            if ( itr > nevtofl2 ){  
645              if ( verbose ) printf(" OrbitalInfo - ERROR: no tof events with entry = %i in Level2 file\n",itr);
646              if ( debug ) printf(" nobefrun %u re %u evfrom %u jumped %u reprocall %i \n",nobefrun,re,evfrom,jumped,reprocall);
647              l0File->Close();
648              code = -901;
649              goto closeandexit;
650            };
651            //
652            tof->Clear();
653            //
654            ttof->GetEntry(itr);
655            //
656          };
657        //        //
658        procev++;        procev++;
659        //        //
# Line 507  int OrbitalInfoCore(UInt_t run, TFile *f Line 661  int OrbitalInfoCore(UInt_t run, TFile *f
661        //        //
662        orbitalinfo->Clear();        orbitalinfo->Clear();
663        //        //
664        // CHANGE HERE!!!!        OrbitalInfoTrkVar *t_orb = new OrbitalInfoTrkVar();
665          if( !(orbitalinfo->OrbitalInfoTrk) ) orbitalinfo->OrbitalInfoTrk = new TClonesArray("OrbitalInfoTrkVar",2);
666          TClonesArray &tor = *orbitalinfo->OrbitalInfoTrk;
667        //        //
668        orbitalinfo->absTime = atime;        // Fill OBT, pkt_num and absTime
669        // EM: add OBT and plt_num infos from the header        //      
       ph = eh->GetPscuHeader();  
670        orbitalinfo->pkt_num = ph->GetCounter();        orbitalinfo->pkt_num = ph->GetCounter();
671        orbitalinfo->OBT = ph->GetOrbitalTime();        orbitalinfo->OBT = ph->GetOrbitalTime();
672          orbitalinfo->absTime = atime;
673        // If the absolute time of the event overpass the time of the        //
674        // tle, get a new tle.  GL_TLE::GetToTime() default to zero.        // Propagate the orbit from the tle time to atime, using SGP(D)4.
675          //
676          cCoordGeo coo;
677          float jyear=0;    
678        //        //
       // I also use this condition to compute the dipole moment dimo.  
       // It's really redundant to compute it so often because  
       // probably it will not change at all.  But the overhead is  
       // minimum.  
       float jyear=0;  
   
679        if(atime >= gltle->GetToTime()) {        if(atime >= gltle->GetToTime()) {
680          if ( !gltle->Query(atime, dbc) ){          if ( !gltle->Query(atime, dbc) ){
681                        //      
682            // Compute the magnetic dipole moment.            // Compute the magnetic dipole moment.
683              //
684            UInt_t year, month, day, hour, min, sec;            UInt_t year, month, day, hour, min, sec;
685                        //
686            TTimeStamp t = TTimeStamp(atime, kTRUE);            TTimeStamp t = TTimeStamp(atime, kTRUE);
687            t.GetDate(kTRUE, 0, &year, &month, &day);            t.GetDate(kTRUE, 0, &year, &month, &day);
688            t.GetTime(kTRUE, 0, &hour, &min, &sec);            t.GetTime(kTRUE, 0, &hour, &min, &sec);
689            jyear = (float) year            jyear = (float) year
690              + (month*31.+ (float) day)/365.              + (month*31.+ (float) day)/365.
691              + (hour*3600.+min*60.+(float)sec)/(24*3600*365.);              + (hour*3600.+min*60.+(float)sec)/(24*3600*365.);
692                        //
693            feldcof_(&jyear, &dimo); // get dipole moment for year            feldcof_(&jyear, &dimo); // get dipole moment for year
694          } else {          } else {
695            code = -56;            code = -56;
696            goto closeandexit;            goto closeandexit;
697          };          };
698        }        }
699          coo = getCoo(atime, gltle->GetFromTime(), gltle->GetTle());
700        // Propagate the orbit from the tle time to atime, using SGP(D)4.        //
701        cCoordGeo coo = getCoo(atime, gltle->GetFromTime(), gltle->GetTle());        cOrbit orbits(*gltle->GetTle());
702          //
703          if ( debug ) printf(" I am Here \n");
704          //
705          // synchronize with quaternions data
706          //
707          if ( isf && neventsm>0 ){
708            if ( debug ) printf(" I am here \n");
709            //
710            // First event
711            //
712            isf = false;
713            upperqtime = atime;
714            lowerqtime = runinfo->RUNHEADER_TIME;
715            for ( ik = 0; ik < neventsm; ik++){
716              l0trm->GetEntry(ik);
717              tmpSize = mcmdev->Records->GetEntries();
718              numrec = tmpSize;
719              for (Int_t j3 = 0;j3<tmpSize;j3++){
720                if ( debug ) printf(" eh eh eh \n");
721                mcmdrc = (pamela::McmdRecord*)mcmdev->Records->At(j3);
722                if ((int)mcmdrc->ID1 == 226){
723                  L_QQ_Q_l_upper->fill(mcmdrc->McmdData);
724                  for (UInt_t ui = 0; ui < 6; ui++){
725                    if (ui>0){
726                      if (L_QQ_Q_l_upper->time[ui]>L_QQ_Q_l_upper->time[0]){
727                        if (dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[ui]*1000-DeltaOBT*1000))<atime){
728                          upperqtime = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[ui]*1000-DeltaOBT*1000));
729                          orbits.getPosition((double) (upperqtime - gltle->GetFromTime())/60., &eCi);
730                          RYPang_upper->TransAngle(eCi.getPos().m_x,eCi.getPos().m_y,eCi.getPos().m_z,eCi.getVel().m_x,eCi.getVel().m_y,eCi.getVel().m_z,L_QQ_Q_l_upper->quat[ui][0],L_QQ_Q_l_upper->quat[ui][1],L_QQ_Q_l_upper->quat[ui][2],L_QQ_Q_l_upper->quat[ui][3]);
731                        }else {
732                          lowerqtime = upperqtime;
733                          upperqtime = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[ui]*1000-DeltaOBT*1000));
734                          orbits.getPosition((double) (upperqtime - gltle->GetFromTime())/60., &eCi);
735                          RYPang_upper->TransAngle(eCi.getPos().m_x,eCi.getPos().m_y,eCi.getPos().m_z,eCi.getVel().m_x,eCi.getVel().m_y,eCi.getVel().m_z,L_QQ_Q_l_upper->quat[ui][0],L_QQ_Q_l_upper->quat[ui][1],L_QQ_Q_l_upper->quat[ui][2],L_QQ_Q_l_upper->quat[ui][3]);
736                          mcreen = j3;
737                          mctren = ik;
738                          if(fgh==0){
739                            CopyQ(L_QQ_Q_l_lower,L_QQ_Q_l_upper);
740                            CopyAng(RYPang_lower,RYPang_upper);
741                          }
742                          oi=ui;
743                          goto closethisloop;
744                        }
745                        fgh++;
746                        CopyQ(L_QQ_Q_l_lower,L_QQ_Q_l_upper);
747                        CopyAng(RYPang_lower,RYPang_upper);
748                      }
749                    }else{
750                      if (dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000))<atime){
751                        upperqtime = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
752                        orbits.getPosition((double) (upperqtime - gltle->GetFromTime())/60., &eCi);
753                        RYPang_upper->TransAngle(eCi.getPos().m_x,eCi.getPos().m_y,eCi.getPos().m_z,eCi.getVel().m_x,eCi.getVel().m_y,eCi.getVel().m_z,L_QQ_Q_l_upper->quat[0][0],L_QQ_Q_l_upper->quat[0][1],L_QQ_Q_l_upper->quat[0][2],L_QQ_Q_l_upper->quat[0][3]);
754                      }
755                      else {
756                        lowerqtime = upperqtime;
757                        upperqtime = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
758                        orbits.getPosition((double) (upperqtime - gltle->GetFromTime())/60., &eCi);
759                        RYPang_upper->TransAngle(eCi.getPos().m_x,eCi.getPos().m_y,eCi.getPos().m_z,eCi.getVel().m_x,eCi.getVel().m_y,eCi.getVel().m_z,L_QQ_Q_l_upper->quat[0][0],L_QQ_Q_l_upper->quat[0][1],L_QQ_Q_l_upper->quat[0][2],L_QQ_Q_l_upper->quat[0][3]);
760                        mcreen = j3;
761                        mctren = ik;
762                        if(fgh==0){
763                          CopyQ(L_QQ_Q_l_lower,L_QQ_Q_l_upper);
764                          CopyAng(RYPang_lower,RYPang_upper);
765                          lowerqtime = atime-1;
766                        }
767                        oi=ui;
768                        goto closethisloop;
769                        //_0 = true;
770                      }
771                      fgh++;
772                      CopyQ(L_QQ_Q_l_lower,L_QQ_Q_l_upper);
773                      CopyAng(RYPang_lower,RYPang_upper);
774                      //_0 = true;
775                    };
776                    //cin>>grib;
777                  };
778                };
779              };
780            };
781          };
782        closethisloop:
783          //
784          if ( debug ) printf(" I am There \n");
785          //
786          if (((atime>(UInt_t)upperqtime)||(atime<(UInt_t)lowerqtime)) && neventsm>0 ){
787            if ( debug ) printf(" I am there \n");
788            //
789            lowerqtime = upperqtime;
790            Long64_t maxloop = 100000000LL;
791            Long64_t mn = 0;
792            bool gh=false;
793            ooi=oi;
794            if ( verbose ) printf(" OrbitalInfoCore: sync with quaternions data upperqtime %u lowerqtime %u atime %u \n",(UInt_t)upperqtime,(UInt_t)lowerqtime,atime);
795            while (!gh){      
796              if ( mn > maxloop ){
797                if ( verbose ) printf(" OrbitalInfoCore: quaternions sync out of range! exiting\n");
798                gh = true;
799                neventsm = 0;
800              };
801              mn++;
802              if (oi<5) oi++;
803              else oi=0;
804              if (oi==0 && numrec > 0){
805                if ( debug ) printf(" mumble \n");
806                mcreen++;
807                if (mcreen == numrec){
808                  mctren++;
809                  mcreen = 0;
810                  l0trm->GetEntry(mctren);
811                  numrec = mcmdev->Records->GetEntries();
812                }
813                CopyQ(L_QQ_Q_l_lower,L_QQ_Q_l_upper);
814                CopyAng(RYPang_lower,RYPang_upper);
815                mcmdrc = (pamela::McmdRecord*)mcmdev->Records->At(mcreen);
816                if ((int)mcmdrc->ID1 == 226){
817                  L_QQ_Q_l_upper->fill(mcmdrc->McmdData);
818                  upperqtime = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
819                  if (upperqtime<lowerqtime){
820                    upperqtime=runinfo->RUNTRAILER_TIME;
821                    CopyQ(L_QQ_Q_l_upper,L_QQ_Q_l_lower);
822                    CopyAng(RYPang_upper,RYPang_lower);
823                  }else{
824                    orbits.getPosition((double) (upperqtime - gltle->GetFromTime())/60., &eCi);
825                    RYPang_upper->TransAngle(eCi.getPos().m_x,eCi.getPos().m_y,eCi.getPos().m_z,eCi.getVel().m_x,eCi.getVel().m_y,eCi.getVel().m_z,L_QQ_Q_l_upper->quat[0][0],L_QQ_Q_l_upper->quat[0][1],L_QQ_Q_l_upper->quat[0][2],L_QQ_Q_l_upper->quat[0][3]);
826                  }
827                  //              re--;
828                  gh=true;
829                }
830              }else{
831                if ((Int_t)L_QQ_Q_l_upper->time[oi]>(Int_t)L_QQ_Q_l_upper->time[0]){
832                  upperqtime = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000));
833                  orbits.getPosition((double) (upperqtime - gltle->GetFromTime())/60., &eCi);
834                  RYPang_upper->TransAngle(eCi.getPos().m_x,eCi.getPos().m_y,eCi.getPos().m_z,eCi.getVel().m_x,eCi.getVel().m_y,eCi.getVel().m_z,L_QQ_Q_l_upper->quat[oi][0],L_QQ_Q_l_upper->quat[oi][1],L_QQ_Q_l_upper->quat[oi][2],L_QQ_Q_l_upper->quat[oi][3]);
835                  orbits.getPosition((double) (lowerqtime - gltle->GetFromTime())/60., &eCi);
836                  RYPang_lower->TransAngle(eCi.getPos().m_x,eCi.getPos().m_y,eCi.getPos().m_z,eCi.getVel().m_x,eCi.getVel().m_y,eCi.getVel().m_z,L_QQ_Q_l_upper->quat[oi-1][0],L_QQ_Q_l_upper->quat[oi-1][1],L_QQ_Q_l_upper->quat[oi-1][2],L_QQ_Q_l_upper->quat[oi-1][3]);
837                  //              re--;
838                  gh=true;
839                };
840              };
841            };
842            if ( verbose ) printf(" OrbitalInfoCore: sync with quaternions data now we have upperqtime %u lowerqtime %u atime %u \n",(UInt_t)upperqtime,(UInt_t)lowerqtime,atime);
843          };
844          //
845          if ( debug ) printf(" I am THIS \n");
846          //
847          // Fill in quaternions and angles
848          //
849          if ((atime<=(UInt_t)upperqtime)&&(atime>=(UInt_t)lowerqtime)&& neventsm>0){      
850            if ( debug ) printf(" I am this \n");
851            UInt_t tut = holeq(lowerqtime, upperqtime, L_QQ_Q_l_lower, L_QQ_Q_l_upper, oi);
852            if (oi == 0){
853              if ((tut!=5)||(tut!=6)){
854                incli = (L_QQ_Q_l_upper->quat[0][0]-L_QQ_Q_l_lower->quat[ooi][0])/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_lower->time[ooi]*1000-DeltaOBT*1000)));
855                orbitalinfo->q0 =  incli*atime+L_QQ_Q_l_upper->quat[0][0]-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
856                incli =     (L_QQ_Q_l_upper->quat[0][1]-L_QQ_Q_l_lower->quat[ooi][1])/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_lower->time[ooi]*1000-DeltaOBT*1000)));
857                orbitalinfo->q1 =  incli*atime+L_QQ_Q_l_upper->quat[0][1]-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
858                incli = (L_QQ_Q_l_upper->quat[0][2]-L_QQ_Q_l_lower->quat[ooi][2])/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_lower->time[ooi]*1000-DeltaOBT*1000)));
859                orbitalinfo->q2 =  incli*atime+L_QQ_Q_l_upper->quat[0][2]-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
860                incli = (L_QQ_Q_l_upper->quat[0][3]-L_QQ_Q_l_lower->quat[ooi][3])/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_lower->time[ooi]*1000-DeltaOBT*1000)));
861                orbitalinfo->q3 =  incli*atime+L_QQ_Q_l_upper->quat[0][3]-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
862            
863                incli = (RYPang_upper->Tangazh-RYPang_lower->Tangazh)/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_lower->time[ooi]*1000-DeltaOBT*1000)));
864                orbitalinfo->theta =  incli*atime+RYPang_upper->Tangazh-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
865                incli = (RYPang_upper->Ryskanie-RYPang_lower->Ryskanie)/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_lower->time[ooi]*1000-DeltaOBT*1000)));
866                orbitalinfo->phi =  incli*atime+RYPang_upper->Ryskanie-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
867                incli = (RYPang_upper->Kren-RYPang_lower->Kren)/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_lower->time[ooi]*1000)));
868                orbitalinfo->etha =  incli*atime+RYPang_upper->Kren-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
869              }
870              if (tut==6){
871                if (fabs(RYPang_lower->Kren-RYPang_upper->Kren)<0.1){
872                  incli = (L_QQ_Q_l_upper->quat[0][0]-L_QQ_Q_l_lower->quat[ooi][0])/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_lower->time[ooi]*1000-DeltaOBT*1000)));
873                  orbitalinfo->q0 =  incli*atime+L_QQ_Q_l_upper->quat[0][0]-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
874                  incli =           (L_QQ_Q_l_upper->quat[0][1]-L_QQ_Q_l_lower->quat[ooi][1])/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_lower->time[ooi]*1000-DeltaOBT*1000)));
875                  orbitalinfo->q1 =  incli*atime+L_QQ_Q_l_upper->quat[0][1]-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
876                  incli = (L_QQ_Q_l_upper->quat[0][2]-L_QQ_Q_l_lower->quat[ooi][2])/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_lower->time[ooi]*1000-DeltaOBT*1000)));
877                  orbitalinfo->q2 =  incli*atime+L_QQ_Q_l_upper->quat[0][2]-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
878                  incli = (L_QQ_Q_l_upper->quat[0][3]-L_QQ_Q_l_lower->quat[ooi][3])/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_lower->time[ooi]*1000-DeltaOBT*1000)));
879                  orbitalinfo->q3 =  incli*atime+L_QQ_Q_l_upper->quat[0][3]-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
880            
881                  incli = (RYPang_upper->Tangazh-RYPang_lower->Tangazh)/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_lower->time[ooi]*1000-DeltaOBT*1000)));
882                  orbitalinfo->theta =  incli*atime+RYPang_upper->Tangazh-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
883                  incli = (RYPang_upper->Ryskanie-RYPang_lower->Ryskanie)/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_lower->time[ooi]*1000-DeltaOBT*1000)));
884                  orbitalinfo->phi =  incli*atime+RYPang_upper->Ryskanie-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
885                  //cout<<"upper = "<<RYPang_upper->Ryskanie<<" lower = "<<RYPang_lower->Ryskanie<<" timeupper[0] = "<<L_QQ_Q_l_upper->time[0]-5500000<<" timelower["<<ooi<<"] = "<<L_QQ_Q_l_lower->time[ooi]-5500000<<" Ryscanie = "<<orbitalinfo->phi<<" incli = "<<incli<<" upper-lower = "<<RYPang_upper->Ryskanie-RYPang_lower->Ryskanie<<" Dtime = "<<dbtime->DBabsTime((UInt_t)L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000)<<"-"<<dbtime->DBabsTime((UInt_t)L_QQ_Q_l_lower->time[ooi]*1000-DeltaOBT*1000)<<" atime = "<<atime<<"\n";
886                  //cin>>grib;
887                  incli = (RYPang_upper->Kren-RYPang_lower->Kren)/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_lower->time[ooi]*1000-DeltaOBT*1000)));
888                  orbitalinfo->etha =  incli*atime+RYPang_upper->Kren-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
889                }
890              }
891            } else {
892              if((tut!=6)||(tut!=7)||(tut!=9)){
893                incli = (L_QQ_Q_l_upper->quat[oi][0]-L_QQ_Q_l_upper->quat[oi-1][0])/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi-1]*1000-DeltaOBT*1000)));
894                orbitalinfo->q0 =  incli*atime+L_QQ_Q_l_upper->quat[oi][0]-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000));
895                incli = (L_QQ_Q_l_upper->quat[oi][1]-L_QQ_Q_l_upper->quat[oi-1][1])/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi-1]*1000-DeltaOBT*1000)));
896                orbitalinfo->q1 =  incli*atime+L_QQ_Q_l_upper->quat[oi][1]-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000));
897                incli = (L_QQ_Q_l_upper->quat[oi][2]-L_QQ_Q_l_upper->quat[oi-1][2])/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi-1]*1000-DeltaOBT*1000)));
898                orbitalinfo->q2 =  incli*atime+L_QQ_Q_l_upper->quat[oi][2]-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000));
899                incli = (L_QQ_Q_l_upper->quat[oi][3]-L_QQ_Q_l_upper->quat[oi-1][3])/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi-1]*1000-DeltaOBT*1000)));
900                orbitalinfo->q3 =  incli*atime+L_QQ_Q_l_upper->quat[oi][3]-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000));
901            
902                incli = (RYPang_upper->Tangazh-RYPang_lower->Tangazh)/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi-1]*1000-DeltaOBT*1000)));
903                orbitalinfo->theta =  incli*atime+RYPang_upper->Tangazh-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000));
904                incli = (RYPang_upper->Ryskanie-RYPang_lower->Ryskanie)/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi-1]*1000-DeltaOBT*1000)));
905                orbitalinfo->phi =  incli*atime+RYPang_upper->Ryskanie-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000));
906                //cout<<"upper = "<<RYPang_upper->Ryskanie<<" lower = "<<RYPang_lower->Ryskanie<<" timeupper["<<oi<<"] = "<<L_QQ_Q_l_upper->time[oi]-5500000<<" timelower["<<oi-1<<"] = "<<L_QQ_Q_l_lower->time[oi-1]-5500000<<" Ryscanie = "<<orbitalinfo->phi<<" incli = "<<incli<<" upper-lower = "<<RYPang_upper->Ryskanie-RYPang_lower->Ryskanie<<" Dtime = "<<dbtime->DBabsTime((UInt_t)L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000)<<"-"<<dbtime->DBabsTime((UInt_t)L_QQ_Q_l_lower->time[oi-1]*1000-DeltaOBT*1000)<<" atime = "<<atime<<"\n";
907                //cin>>grib;
908                incli = (RYPang_upper->Kren-RYPang_lower->Kren)/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi-1]*1000-DeltaOBT*1000)));
909                orbitalinfo->etha =  incli*atime+RYPang_upper->Kren-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000));
910              }
911              if (tut==6){
912                if (fabs(RYPang_lower->Kren-RYPang_upper->Kren)<0.1){
913                  incli = (L_QQ_Q_l_upper->quat[oi][0]-L_QQ_Q_l_upper->quat[oi-1][0])/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi-1]*1000-DeltaOBT*1000)));
914                  orbitalinfo->q0 =  incli*atime+L_QQ_Q_l_upper->quat[oi][0]-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000));
915                  incli = (L_QQ_Q_l_upper->quat[oi][1]-L_QQ_Q_l_upper->quat[oi-1][1])/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi-1]*1000-DeltaOBT*1000)));
916                  orbitalinfo->q1 =  incli*atime+L_QQ_Q_l_upper->quat[oi][1]-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000));
917                  incli = (L_QQ_Q_l_upper->quat[oi][2]-L_QQ_Q_l_upper->quat[oi-1][2])/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi-1]*1000-DeltaOBT*1000)));
918                  orbitalinfo->q2 =  incli*atime+L_QQ_Q_l_upper->quat[oi][2]-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000));
919                  incli = (L_QQ_Q_l_upper->quat[oi][3]-L_QQ_Q_l_upper->quat[oi-1][3])/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi-1]*1000-DeltaOBT*1000)));
920                  orbitalinfo->q3 =  incli*atime+L_QQ_Q_l_upper->quat[oi][3]-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000));
921            
922                  incli = (RYPang_upper->Tangazh-RYPang_lower->Tangazh)/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi-1]*1000-DeltaOBT*1000)));
923                  orbitalinfo->theta =  incli*atime+RYPang_upper->Tangazh-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000));
924                  incli = (RYPang_upper->Ryskanie-RYPang_lower->Ryskanie)/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi-1]*1000-DeltaOBT*1000)));
925                  orbitalinfo->phi =  incli*atime+RYPang_upper->Ryskanie-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000));
926                  //cout<<"upper = "<<RYPang_upper->Ryskanie<<" lower = "<<RYPang_lower->Ryskanie<<" timeupper["<<oi<<"] = "<<L_QQ_Q_l_upper->time[oi]-5500000<<" timelower["<<oi-1<<"] = "<<L_QQ_Q_l_lower->time[oi-1]-5500000<<" Ryscanie = "<<orbitalinfo->phi<<" incli = "<<incli<<" upper-lower = "<<RYPang_upper->Ryskanie-RYPang_lower->Ryskanie<<" Dtime = "<<dbtime->DBabsTime((UInt_t)L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000)<<"-"<<dbtime->DBabsTime((UInt_t)L_QQ_Q_l_lower->time[oi-1]*1000-DeltaOBT*1000)<<" atime = "<<atime<<"\n";
927                  //cin>>grib;
928                  incli = (RYPang_upper->Kren-RYPang_lower->Kren)/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi-1]*1000-DeltaOBT*1000)));
929                  orbitalinfo->etha =  incli*atime+RYPang_upper->Kren-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000));
930                }
931              }            
932            }
933            //
934            orbitalinfo->mode = holeq(lowerqtime, upperqtime, L_QQ_Q_l_lower, L_QQ_Q_l_upper, oi);
935            //
936          } else {
937            if ( debug ) printf(" ops no incl! \n");
938            orbitalinfo->mode = 10;
939          };
940          //
941          // ops no inclination information
942          //
943          if ( orbitalinfo->q0< -999 || orbitalinfo->q1 < -999 || orbitalinfo->q2 < -999 || orbitalinfo->q3 < -999 || orbitalinfo->q0 != orbitalinfo->q0 || orbitalinfo->q1 != orbitalinfo->q1 || orbitalinfo->q2 != orbitalinfo->q2 || orbitalinfo->q3 != orbitalinfo->q3 ){
944            orbitalinfo->mode = 10;
945            orbitalinfo->q0 = -1000.;
946            orbitalinfo->q1 = -1000.;
947            orbitalinfo->q2 = -1000.;
948            orbitalinfo->q3 = -1000.;
949            orbitalinfo->etha = -1000.;
950            orbitalinfo->phi = -1000.;
951            orbitalinfo->theta = -1000.;
952          };
953          //
954          // #########################################################################################################################  
955          //
956          // fill orbital positions
957          //        
958        // Build coordinates in the right range.  We want to convert,        // Build coordinates in the right range.  We want to convert,
959        // longitude from (0, 2*pi) to (-180deg, 180deg).  Altitude is        // longitude from (0, 2*pi) to (-180deg, 180deg).  Altitude is
960        // in meters.        // in meters.
961        lon = (coo.m_Lon > M_PI) ? rad2deg(coo.m_Lon - 2*M_PI) : rad2deg(coo.m_Lon);        lon = (coo.m_Lon > M_PI) ? rad2deg(coo.m_Lon - 2*M_PI) : rad2deg(coo.m_Lon);
962        lat = rad2deg(coo.m_Lat);        lat = rad2deg(coo.m_Lat);
963        alt = coo.m_Alt;        alt = coo.m_Alt;
964          //
   
       //      if((lon>180) || (lon<-180) || (lat>90) || (lat<-90) || (alt<0))  
       //        continue;  
965        if( lon<180 && lon>-180 && lat<90 && lat>-90 && alt>0 ){          if( lon<180 && lon>-180 && lat<90 && lat>-90 && alt>0 ){  
966            //      
967          orbitalinfo->lon = lon;          orbitalinfo->lon = lon;
968          orbitalinfo->lat = lat;          orbitalinfo->lat = lat;
969          orbitalinfo->alt = alt ;          orbitalinfo->alt = alt ;
970                    //
971          // compute mag field components and L shell.          // compute mag field components and L shell.
972            //
973          feldg_(&lat, &lon, &alt, &bnorth, &beast, &bdown, &babs);          feldg_(&lat, &lon, &alt, &bnorth, &beast, &bdown, &babs);
974          shellg_(&lat, &lon, &alt, &dimo, &xl, &icode, &bab1);          shellg_(&lat, &lon, &alt, &dimo, &xl, &icode, &bab1);
975          findb0_(&stps, &bdel, &value, &bequ, &rr0);          findb0_(&stps, &bdel, &value, &bequ, &rr0);
976                    //
977          orbitalinfo->Bnorth = bnorth;          orbitalinfo->Bnorth = bnorth;
978          orbitalinfo->Beast = beast;          orbitalinfo->Beast = beast;
979          orbitalinfo->Bdown = bdown;          orbitalinfo->Bdown = bdown;
980          orbitalinfo->Babs = babs;          orbitalinfo->Babs = babs;
981          orbitalinfo->BB0 = babs/bequ;          orbitalinfo->BB0 = babs/bequ;
982          orbitalinfo->L = xl;          orbitalinfo->L = xl;      
           
983          // Set Stormer vertical cutoff using L shell.          // Set Stormer vertical cutoff using L shell.
984          orbitalinfo->cutoff[0] = 14.9/(xl*xl);          orbitalinfo->cutoff[0] = 14.9/(xl*xl);
985            //
986          };      
987          //
988          // pitch angles
989          //
990          if ( orbitalinfo->mode != 10 && orbitalinfo->mode != 5 && orbitalinfo->mode !=7 && orbitalinfo->mode != 9 ){
991            //
992            Float_t Bx = -orbitalinfo->Bdown;                       //don't need for PamExp ExpOnly for all geography areas
993            Float_t By = orbitalinfo->Beast;                        //don't need for PamExp ExpOnly for all geography areas
994            Float_t Bz = orbitalinfo->Bnorth;                       //don't need for PamExp ExpOnly for all geography areas
995            //
996            TMatrixD Fij = PO->ECItoGreenwich(PO->QuatoECI(orbitalinfo->q0,orbitalinfo->q1,orbitalinfo->q2,orbitalinfo->q3),orbitalinfo->absTime);
997            TMatrixD Dij = PO->GreenwichtoGEO(orbitalinfo->lat,orbitalinfo->lon,Fij);
998            TMatrixD Iij = PO->ColPermutation(Dij);
999            //
1000            A1 = Iij(0,2);
1001            A2 = Iij(1,2);
1002            A3 = Iij(2,2);
1003            //      
1004            orbitalinfo->pamzenitangle = (Float_t)PO->GetPitchAngle(1,0,0,A1,A2,A3);                        // Angle between zenit and Pamela's main axiz
1005            orbitalinfo->pamBangle = (Float_t)PO->GetPitchAngle(A1,A2,A3,Bx,By,Bz);                 // Angle between Pamela's main axiz and B
1006            //
1007            if ( !standalone && tof->ntrk() > 0 ){
1008              //
1009              Int_t nn = 0;
1010              for(Int_t nt=0; nt < tof->ntrk(); nt++){  
1011                //
1012                ToFTrkVar *ptt = tof->GetToFTrkVar(nt);
1013                Double_t E11x = ptt->xtr_tof[0]; // tr->x[0];
1014                Double_t E11y = ptt->ytr_tof[0]; //tr->y[0];
1015                Double_t E11z = zin[0];
1016                Double_t E22x = ptt->xtr_tof[3];//tr->x[3];
1017                Double_t E22y = ptt->ytr_tof[3];//tr->y[3];
1018                Double_t E22z = zin[3];
1019                if ( E11x < 100. && E11y < 100. && E22x < 100. && E22y < 100.  ){
1020                  Double_t norm = sqrt(pow(E22x-E11x,2)+pow(E22y-E11y,2)+pow(E22z-E11z,2));
1021                  Double_t MyAzim = TMath::RadToDeg()*atan(TMath::Abs(E22y-E11y)/TMath::Abs(E22x-E11x));
1022                  if(E22x-E11x>=0 && E22y-E11y <0) MyAzim =  360. - MyAzim;
1023                  if(E22x-E11x>=0 && E22y-E11y >=0) MyAzim = MyAzim;
1024                  if(E22x-E11x<0 && E22y-E11y >0) MyAzim = 180. - MyAzim;
1025                  if(E22x-E11x<0 && E22y-E11y <0) MyAzim = 180. + MyAzim;
1026                  Px = (E22x-E11x)/norm;
1027                  Py = (E22y-E11y)/norm;
1028                  Pz = (E22z-E11z)/norm;
1029                  //
1030                  TMatrixD Eij = PO->PamelatoGEO(Iij,Px,Py,Pz);            
1031                  //            
1032                  t_orb->trkseqno = ptt->trkseqno;
1033                  t_orb->pitch = (Float_t)PO->GetPitchAngle(Eij(0,0),Eij(1,0),Eij(2,0),Bx,By,Bz);
1034                  //
1035                  new(tor[nn]) OrbitalInfoTrkVar(*t_orb);
1036                  nn++;
1037                  //
1038                  t_orb->Clear();
1039                  //
1040                };
1041                //
1042              };
1043            } else {
1044              if ( debug ) printf(" mmm... mode %u standalone %i ntrk %i \n",orbitalinfo->mode,standalone,tof->ntrk());
1045            };
1046            //
1047        };        };
1048                //
1049        // end EM        // Fill the class
1050          //
1051        OrbitalInfotr->Fill();        OrbitalInfotr->Fill();
1052      }        //
1053              delete t_orb;
1054          //
1055        }; // loop over the events in the run
1056      //      //
1057      // Here you may want to clear some variables before processing another run        // Here you may want to clear some variables before processing another run  
1058      //      //
1059      delete dbtime;      delete dbtime;
1060        if ( L_QQ_Q_l_upper ) delete L_QQ_Q_l_upper;
1061        if ( L_QQ_Q_l_lower ) delete L_QQ_Q_l_lower;
1062        if ( RYPang_upper ) delete RYPang_upper;
1063        if ( RYPang_lower ) delete RYPang_lower;
1064    }; // process all the runs    }; // process all the runs
1065    //    
1066    if (verbose) printf("\n Finished processing data \n");    if (verbose) printf("\n Finished processing data \n");
1067    //    //
1068   closeandexit:   closeandexit:
# Line 626  int OrbitalInfoCore(UInt_t run, TFile *f Line 1100  int OrbitalInfoCore(UInt_t run, TFile *f
1100    //    //
1101    if ( l0File ) l0File->Close();    if ( l0File ) l0File->Close();
1102    if ( tempfile ) tempfile->Close();                if ( tempfile ) tempfile->Close();            
1103    gSystem->Unlink(tempname.str().c_str());    if ( myfold ) gSystem->Unlink(tempname.str().c_str());
1104    //    //
1105    if ( runinfo ) runinfo->Close();        if ( runinfo ) runinfo->Close();    
1106    if ( OrbitalInfotr ) OrbitalInfotr->SetName("OrbitalInfo");        if ( OrbitalInfotr ) OrbitalInfotr->SetName("OrbitalInfo");    
1107      if ( tof ) tof->Delete();
1108      if ( ttof ) ttof->Delete();
1109      //
1110    if ( file ){    if ( file ){
1111      file->cd();      file->cd();
1112      file->Write();      file->Write();
1113    };    };
1114    //    //
1115    gSystem->Unlink(OrbitalInfofolder.str().c_str());    if ( myfold ) gSystem->Unlink(OrbitalInfofolder.str().c_str());
1116    //    //
1117    // the end    // the end
1118    //    //
1119      if ( dbc ){
1120        dbc->Close();
1121        delete dbc;
1122      };
1123    if (verbose) printf("\n Exiting...\n");    if (verbose) printf("\n Exiting...\n");
1124    if(OrbitalInfotr)OrbitalInfotr->Delete();    if(OrbitalInfotr)OrbitalInfotr->Delete();
1125    //    //
1126      if ( PO ) delete PO;
1127    if ( orbitalinfo ) delete orbitalinfo;    if ( orbitalinfo ) delete orbitalinfo;
1128    if ( orbitalinfoclone ) delete orbitalinfoclone;    if ( orbitalinfoclone ) delete orbitalinfoclone;
1129    if ( glroot ) delete glroot;    if ( glroot ) delete glroot;
# Line 663  cCoordGeo getCoo(UInt_t atime, UInt_t tl Line 1145  cCoordGeo getCoo(UInt_t atime, UInt_t tl
1145  {  {
1146    cEci eci;    cEci eci;
1147    cOrbit orbit(*tle);    cOrbit orbit(*tle);
     
1148    orbit.getPosition((double) (atime - tletime)/60., &eci);    orbit.getPosition((double) (atime - tletime)/60., &eci);
1149        
1150    return eci.toGeo();    return eci.toGeo();
1151  }  }
1152    
1153    // function of copyng of quatrnions classes
1154    
1155    void CopyQ(Quaternions *Q1, Quaternions *Q2){
1156      for(UInt_t i = 0; i < 6; i++){
1157        Q1->time[i]=Q2->time[i];
1158        for (UInt_t j = 0; j < 4; j++)Q1->quat[i][j]=Q2->quat[i][j];
1159      }
1160      return;
1161    }
1162    
1163    // functions of copyng InclinationInfo classes
1164    
1165    void CopyAng(InclinationInfo *A1, InclinationInfo *A2){
1166      A1->Tangazh = A2->Tangazh;
1167      A1->Ryskanie = A2->Ryskanie;
1168      A1->Kren = A2->Kren;
1169      return;
1170    }
1171    
1172    UInt_t holeq(Double_t lower,Double_t upper,Quaternions *Qlower, Quaternions *Qupper, UInt_t f){
1173      
1174      UInt_t hole = 10;
1175      bool R10l = false;     // Sign of R10 mode in lower quaternions array
1176      bool R10u = false;     // Sign of R10 mode in upper quaternions array
1177      bool insm = false;     // Sign that we inside quaternions array
1178      bool mxtml = false;    // Sign of mixt mode in lower quaternions array
1179      bool mxtmu = false;    // Sign of mixt mode in upper quaternions array
1180      bool npasm = false;     // Sign of normall pass between R10 and non R10 or between non R10 and R10
1181      UInt_t NCQl = 6;       // Number of correct quaternions in lower array
1182      UInt_t NCQu = 6;       // Number of correct quaternions in upper array
1183      if (f>0){
1184        insm = true;
1185        if(Qupper->time[f]-Qupper->time[f-1]==30) R10u = false;
1186        if(Qupper->time[f]-Qupper->time[f-1]<1) R10u = true;
1187      }else{
1188        insm = false;
1189        if((Qlower->time[5]-Qlower->time[0]<2)&&(Qlower->time[1]-Qlower->time[0]<2)) R10l = true;
1190        if((Qupper->time[5]-Qupper->time[0]<2)&&(Qupper->time[1]-Qupper->time[0]<2)) R10u = true;
1191        if((Qlower->time[5]-Qlower->time[0]==150)&&(Qlower->time[1]-Qlower->time[0]==30)) R10l = false;
1192        if((Qupper->time[5]-Qupper->time[0]==150)&&(Qupper->time[1]-Qupper->time[0]==30)) R10u = false;
1193        if((Qlower->time[5]-Qlower->time[0]<2)&&(Qlower->time[1]-Qlower->time[0]==30)){
1194          mxtml = true;
1195          for(UInt_t i = 1; i < 6; i++){
1196            if(Qlower->time[i]-Qlower->time[0]==30*i) NCQl=i;
1197          }
1198        }
1199        if((Qupper->time[5]-Qupper->time[0]<2)&&(Qupper->time[1]-Qupper->time[0]==30)){
1200          mxtmu = true;
1201          for(UInt_t i = 1; i < 6; i++){
1202            if(Qupper->time[i]-Qupper->time[0]==30*i) NCQu=i;
1203          }
1204        }
1205      }
1206      
1207      if(((upper-lower==1.5)||(upper-lower==3.)||(upper-lower==30.)||(upper-lower==31.5)||(upper-lower==33.)||(upper-lower==181.5)||(upper-lower==210.)||(upper-lower==211.5))&&!insm) npasm = true;
1208      
1209      
1210      if (R10u&&insm) hole=0; // best event R10
1211      if ((upper-lower<=5)&&(!insm)&&R10l&&R10u) hole = 1; // when first of 6 quaternions in array is correct
1212      if (((!R10u)&&insm)||((!insm)&&(!R10u)&&(!R10l)&&((upper-lower==210+(6-NCQl)*30)||(upper-lower==30)))) hole = 2; //non R10
1213      if (npasm&&(!insm)&&((R10l&&!R10u)||(R10u&&!R10l))) hole = 3; //normall pass from R10 to non R10 or from non R10 to R10
1214      if ((!npasm)&&(upper-lower<=300)&&(!insm)&&((R10l&&!R10u)||(R10u&&!R10l))) hole = 4; // eliminable hole between R10 and non R10 or between non R10 and R10
1215      if ((upper-lower>=300)&&(!insm)&&((R10l&&!R10u)||(R10u&&!R10l))) hole = 5; //uneliminable hole between R10 and non R10 or between non R10 and R10
1216      if ((upper-lower>5)&&(upper-lower<=300)&&R10u&&R10l) hole = 6; // eliminable hole inside R10
1217      if ((upper-lower>300)&&R10u&&R10l) hole = 7; //uneliminable hole inside R10
1218      if ((upper-lower>210)&&(upper-lower<=1200)&&(!R10u)&&(!R10l)) hole = 8; //eliminable hole inside non R10
1219      if ((upper-lower>1200)&&!R10u&&!R10l) hole = 9; // uneliminable hole inside non R10
1220      return hole;
1221    }
1222    
Legend:
Removed from v.1.10  
changed lines
  Added in v.1.30
  ViewVC Help
Powered by ViewVC 1.1.23