/[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.18 by mocchiut, Thu May 3 20:09:29 2007 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 48  using namespace std; Line 51  using namespace std;
51  //  //
52  //  //
53  int OrbitalInfoCore(UInt_t run, TFile *file, TSQLServer *dbc, Int_t OrbitalInfoargc, char *OrbitalInfoargv[]){  int OrbitalInfoCore(UInt_t run, TFile *file, TSQLServer *dbc, 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    //    //
57    TString processFolder = "OrbitalInfoFolder";    TString processFolder = Form("OrbitalInfoFolder_%u",run);
58    //    //
59    // Set these to true to have a very verbose output.    // Set these to true to have a very verbose output.
60    //    //
# Line 86  int OrbitalInfoCore(UInt_t run, TFile *f Line 74  int OrbitalInfoCore(UInt_t run, TFile *f
74        };        };
75        if ( (!strcmp(OrbitalInfoargv[i],"--debug")) || (!strcmp(OrbitalInfoargv[i],"-g")) ) {        if ( (!strcmp(OrbitalInfoargv[i],"--debug")) || (!strcmp(OrbitalInfoargv[i],"-g")) ) {
76          verbose = true;          verbose = true;
77            debug = true;
78        };        };
79        if ( (!strcmp(OrbitalInfoargv[i],"--verbose")) || (!strcmp(OrbitalInfoargv[i],"-v")) ) {        if ( (!strcmp(OrbitalInfoargv[i],"--verbose")) || (!strcmp(OrbitalInfoargv[i],"-v")) ) {
80          verbose = true;          verbose = true;
# Line 98  int OrbitalInfoCore(UInt_t run, TFile *f Line 87  int OrbitalInfoCore(UInt_t run, TFile *f
87    //    //
88    TTree *OrbitalInfotr = 0;    TTree *OrbitalInfotr = 0;
89    UInt_t nevents = 0;    UInt_t nevents = 0;
90      UInt_t neventsm = 0;
91    //    //
92    // variables needed to reprocess data    // variables needed to reprocess data
93    //    //
# Line 116  int OrbitalInfoCore(UInt_t run, TFile *f Line 106  int OrbitalInfoCore(UInt_t run, TFile *f
106    stringstream ftmpname;    stringstream ftmpname;
107    TString fname;    TString fname;
108    UInt_t totfileentries = 0;    UInt_t totfileentries = 0;
109    UInt_t idRun = 0;    UInt_t idRun = 0;
110      //
111      // My variables. Vitaly.
112      //
113      //  UInt_t iev = 0;
114      //  UInt_t j3 = 0;
115      UInt_t oi = 0;
116      Int_t tmpSize = 0;
117    //    //
118    // variables needed to handle error signals    // variables needed to handle error signals
119    //    //
# Line 132  int OrbitalInfoCore(UInt_t run, TFile *f Line 129  int OrbitalInfoCore(UInt_t run, TFile *f
129    //    //
130    TFile *l0File = 0;    TFile *l0File = 0;
131    TTree *l0tr = 0;    TTree *l0tr = 0;
132      TTree *l0trm = 0;
133    // EM: open also header branch    // EM: open also header branch
134    TBranch *l0head = 0;    TBranch *l0head = 0;
135    pamela::EventHeader *eh = 0;    pamela::EventHeader *eh = 0;
136    pamela::PscuHeader *ph = 0;    pamela::PscuHeader *ph = 0;
137      pamela::McmdEvent *mcmdev = 0;
138      pamela::McmdRecord *mcmdrc = 0;
139    // end EM    // end EM
140      
141      //  pamela::RunHeaderEvent *reh = new pamela::RunHeaderEvent;
142      //  pamela::EventHeader    *eH  = new pamela::EventHeader;
143      
144    //    //
145    // Define other basic variables    // Define other basic variables
146    //    //
# Line 147  int OrbitalInfoCore(UInt_t run, TFile *f Line 151  int OrbitalInfoCore(UInt_t run, TFile *f
151    Int_t totevent = 0;    Int_t totevent = 0;
152    UInt_t atime = 0;    UInt_t atime = 0;
153    UInt_t re = 0;    UInt_t re = 0;
154      UInt_t ik = 0;
155    
156    // Position    // Position
157    Float_t lon, lat, alt;    Float_t lon, lat, alt;
# Line 193  int OrbitalInfoCore(UInt_t run, TFile *f Line 198  int OrbitalInfoCore(UInt_t run, TFile *f
198    GL_ROOT *glroot = new GL_ROOT();    GL_ROOT *glroot = new GL_ROOT();
199    GL_TIMESYNC *dbtime = 0;    GL_TIMESYNC *dbtime = 0;
200    GL_TLE *gltle = new GL_TLE();    GL_TLE *gltle = new GL_TLE();
201      //
202      //Quaternions classes
203      //
204      Quaternions *L_QQ_Q_l_lower = new Quaternions();
205      InclinationInfo *RYPang_lower = new InclinationInfo();
206      Quaternions *L_QQ_Q_l_upper = new Quaternions();
207      InclinationInfo *RYPang_upper = new InclinationInfo();
208      
209      cEci eCi;
210      
211    // Initialize fortran routines!!!    // Initialize fortran routines!!!
212    Int_t ltp2 = 0;    Int_t ltp2 = 0;
213    Int_t ltp3 = 0;    Int_t ltp3 = 0;
# Line 201  int OrbitalInfoCore(UInt_t run, TFile *f Line 216  int OrbitalInfoCore(UInt_t run, TFile *f
216    GL_PARAM *glparam = new GL_PARAM();    GL_PARAM *glparam = new GL_PARAM();
217    GL_PARAM *glparam2 = new GL_PARAM();    GL_PARAM *glparam2 = new GL_PARAM();
218    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
219      
220    if ( parerror<0 ) {    if ( parerror<0 ) {
221      code = parerror;      code = parerror;
222      goto closeandexit;      goto closeandexit;
# Line 217  int OrbitalInfoCore(UInt_t run, TFile *f Line 233  int OrbitalInfoCore(UInt_t run, TFile *f
233    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());
234    //    //
235    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);
   
236    //    //
237    // End IGRF stuff//    // End IGRF stuff//
238    //    //
# Line 383  int OrbitalInfoCore(UInt_t run, TFile *f Line 398  int OrbitalInfoCore(UInt_t run, TFile *f
398    //    //
399    // Loop over the run to be processed    // Loop over the run to be processed
400    //    //
401      
402    for (UInt_t irun=0; irun < numbofrun; irun++){    for (UInt_t irun=0; irun < numbofrun; irun++){
403      //      //
404      // retrieve the first run ID to be processed using the RunInfo list      // retrieve the first run ID to be processed using the RunInfo list
405      //      //
406        
407      idRun = runlist->At(irun);      idRun = runlist->At(irun);
408      if (verbose){      if (verbose){
409        printf("\n\n\n ####################################################################### \n");        printf("\n\n\n ####################################################################### \n");
# Line 419  int OrbitalInfoCore(UInt_t run, TFile *f Line 436  int OrbitalInfoCore(UInt_t run, TFile *f
436      // prepare the timesync for the db      // prepare the timesync for the db
437      //      //
438      dbtime = new GL_TIMESYNC(runinfo->ID_ROOT_L0,"ID",dbc);      dbtime = new GL_TIMESYNC(runinfo->ID_ROOT_L0,"ID",dbc);
439      
440      //      //
441      // 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.
442      //      //
# Line 428  int OrbitalInfoCore(UInt_t run, TFile *f Line 446  int OrbitalInfoCore(UInt_t run, TFile *f
446      ftmpname << glroot->PATH.Data() << "/";      ftmpname << glroot->PATH.Data() << "/";
447      ftmpname << glroot->NAME.Data();      ftmpname << glroot->NAME.Data();
448      fname = ftmpname.str().c_str();      fname = ftmpname.str().c_str();
449        ftmpname.str("");
450      //      //
451      // print out informations      // print out informations
452      //      //
453      totevent = runinfo->NEVENTS;      totevent = runinfo->NEVENTS;
454        //cout<<"totevents = "<<totevent<<"\n";
455      if (verbose){      if (verbose){
456        printf("\n LEVEL0 data file: %s \n",fname.Data());        printf("\n LEVEL0 data file: %s \n",fname.Data());
457        printf(" RUN HEADER absolute time is:  %u \n",runinfo->RUNHEADER_TIME);        printf(" RUN HEADER absolute time is:  %u \n",runinfo->RUNHEADER_TIME);
458        printf(" RUN TRAILER absolute time is: %u \n",runinfo->RUNTRAILER_TIME);        printf(" RUN TRAILER absolute time is: %u \n",runinfo->RUNTRAILER_TIME);
459        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);
460      }//      }//
461      // Open Level0 file      // Open Level0 file
462      l0File = new TFile(fname.Data());      l0File = new TFile(fname.Data());
# Line 478  int OrbitalInfoCore(UInt_t run, TFile *f Line 498  int OrbitalInfoCore(UInt_t run, TFile *f
498        goto closeandexit;        goto closeandexit;
499      };      };
500      //      //
501    //     TTree *tp = (TTree*)l0File->Get("RunHeader");
502    //     tp->SetBranchAddress("Header", &eH);
503    //     tp->SetBranchAddress("RunHeader", &reh);
504    //     tp->GetEntry(0);
505    //     ph = eH->GetPscuHeader();
506    //     ULong_t TimeSync = reh->LAST_TIME_SYNC_INFO;
507    //     ULong_t ObtSync = reh->OBT_TIME_SYNC;    
508    //     if ( debug ) printf(" 1 TimeSync %lu ObtSync %lu DeltaOBT %lu\n",TimeSync,ObtSync,TimeSync-ObtSync);
509    //
510        ULong_t TimeSync = (ULong_t)dbtime->GetTimesync();
511        ULong_t ObtSync = (ULong_t)(dbtime->GetObt0()/1000);
512        ULong_t DeltaOBT = TimeSync - ObtSync;
513    
514        if ( debug ) printf(" 2 TimeSync %lu ObtSync %lu DeltaOBT %lu\n",(ULong_t)(dbtime->GetTimesync()/1000),(ULong_t)dbtime->GetObt0(),TimeSync-ObtSync);
515        
516        l0trm = (TTree*)l0File->Get("Mcmd");
517        neventsm = l0trm->GetEntries();
518        //    neventsm = 0;
519        //
520        if (neventsm == 0){
521          if ( debug ) printf("InclinationInfo - WARNING: No quaternions in this File");
522          //      l0File->Close();
523          code = 900;
524          //      goto closeandexit;
525        }
526        //
527        
528        l0trm->SetBranchAddress("Mcmd", &mcmdev);
529        //    l0trm->SetBranchAddress("Header", &eh);
530        //
531        //
532        //
533        UInt_t mctren = 0;    
534        UInt_t mcreen = 0;  
535        UInt_t numrec = 0;
536        //
537        Double_t upperqtime = 0;
538        Double_t lowerqtime = 0;
539        
540        Double_t incli = 0;
541        oi = 0;
542        UInt_t ooi = 0;
543        //
544        // init quaternions sync
545        //
546        Bool_t isf = true;
547        Int_t fgh = 0;
548        //
549      // run over all the events of the run      // run over all the events of the run
550      //      //
551      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");
552      //      //
553      for ( re = runinfo->EV_FROM; re < (runinfo->EV_FROM+runinfo->NEVENTS); re++){      for ( re = runinfo->EV_FROM; re < (runinfo->EV_FROM+runinfo->NEVENTS); re++){
554        
555        //        //
556        if ( procev%1000 == 0 && procev > 0 && verbose ) printf(" %iK \n",procev/1000);          if ( procev%1000 == 0 && procev > 0 && verbose ) printf(" %iK \n",procev/1000);  
557          if ( debug ) printf(" %i \n",procev);      
558        //        //
559        l0head->GetEntry(re);        l0head->GetEntry(re);
560        //        //
561        // absolute time of this event        // absolute time of this event
562        //        //
563        ph = eh->GetPscuHeader();        ph = eh->GetPscuHeader();
564        atime = dbtime->DBabsTime(ph->GetOrbitalTime());          atime = dbtime->DBabsTime(ph->GetOrbitalTime());
565        //        //
566        // paranoid check        // paranoid check
567        //        //
# Line 507  int OrbitalInfoCore(UInt_t run, TFile *f Line 577  int OrbitalInfoCore(UInt_t run, TFile *f
577        //        //
578        orbitalinfo->Clear();        orbitalinfo->Clear();
579        //        //
580        // CHANGE HERE!!!!        // Fill OBT, pkt_num and absTime
581        //        //      
582        orbitalinfo->absTime = atime;        //      ph = eh->GetPscuHeader();
       // EM: add OBT and plt_num infos from the header  
       ph = eh->GetPscuHeader();  
583        orbitalinfo->pkt_num = ph->GetCounter();        orbitalinfo->pkt_num = ph->GetCounter();
584        orbitalinfo->OBT = ph->GetOrbitalTime();        orbitalinfo->OBT = ph->GetOrbitalTime();
585          orbitalinfo->absTime = atime;
586        // If the absolute time of the event overpass the time of the        //
587        // tle, get a new tle.  GL_TLE::GetToTime() default to zero.        // Propagate the orbit from the tle time to atime, using SGP(D)4.
588          //
589          cCoordGeo coo;
590          float jyear=0;    
591        //        //
       // 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;  
   
592        if(atime >= gltle->GetToTime()) {        if(atime >= gltle->GetToTime()) {
593          if ( !gltle->Query(atime, dbc) ){          if ( !gltle->Query(atime, dbc) ){
594                        //      
595            // Compute the magnetic dipole moment.            // Compute the magnetic dipole moment.
596              //
597            UInt_t year, month, day, hour, min, sec;            UInt_t year, month, day, hour, min, sec;
598                        //
599            TTimeStamp t = TTimeStamp(atime, kTRUE);            TTimeStamp t = TTimeStamp(atime, kTRUE);
600            t.GetDate(kTRUE, 0, &year, &month, &day);            t.GetDate(kTRUE, 0, &year, &month, &day);
601            t.GetTime(kTRUE, 0, &hour, &min, &sec);            t.GetTime(kTRUE, 0, &hour, &min, &sec);
602            jyear = (float) year            jyear = (float) year
603              + (month*31.+ (float) day)/365.              + (month*31.+ (float) day)/365.
604              + (hour*3600.+min*60.+(float)sec)/(24*3600*365.);              + (hour*3600.+min*60.+(float)sec)/(24*3600*365.);
605                        //
606            feldcof_(&jyear, &dimo); // get dipole moment for year            feldcof_(&jyear, &dimo); // get dipole moment for year
607          } else {          } else {
608            code = -56;            code = -56;
609            goto closeandexit;            goto closeandexit;
610          };          };
611        }        }
612          coo = getCoo(atime, gltle->GetFromTime(), gltle->GetTle());
613        // Propagate the orbit from the tle time to atime, using SGP(D)4.        //
614        cCoordGeo coo = getCoo(atime, gltle->GetFromTime(), gltle->GetTle());        cOrbit orbits(*gltle->GetTle());
615          //
616          if ( debug ) printf(" I am Here \n");
617          //
618          // synchronize with quaternions data
619          //
620          if ( isf && neventsm>0 ){
621            if ( debug ) printf(" I am here \n");
622            //
623            // First event
624            //
625            isf = false;
626            upperqtime = atime;
627            lowerqtime = runinfo->RUNHEADER_TIME;
628            for ( ik = 0; ik < neventsm; ik++){
629              l0trm->GetEntry(ik);
630              tmpSize = mcmdev->Records->GetEntries();
631              numrec = tmpSize;
632              for (Int_t j3 = 0;j3<tmpSize;j3++){
633                mcmdrc = (pamela::McmdRecord*)mcmdev->Records->At(j3);
634                if ((int)mcmdrc->ID1 == 226){
635                  L_QQ_Q_l_upper->fill(mcmdrc->McmdData);
636                  for (UInt_t ui = 0; ui < 6; ui++){
637                    if (ui>0){
638                      if (L_QQ_Q_l_upper->time[ui]>L_QQ_Q_l_upper->time[0]){
639                        if (dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[ui]*1000-DeltaOBT*1000))<atime){
640                          upperqtime = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[ui]*1000-DeltaOBT*1000));
641                          orbits.getPosition((double) (upperqtime - gltle->GetFromTime())/60., &eCi);
642                          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]);
643                        }else {
644                          lowerqtime = upperqtime;
645                          upperqtime = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[ui]*1000-DeltaOBT*1000));
646                          orbits.getPosition((double) (upperqtime - gltle->GetFromTime())/60., &eCi);
647                          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]);
648                          mcreen = j3;
649                          mctren = ik;
650                          if(fgh==0){
651                            CopyQ(L_QQ_Q_l_lower,L_QQ_Q_l_upper);
652                            CopyAng(RYPang_lower,RYPang_upper);
653                          }
654                          oi=ui;
655                          goto closethisloop;
656                        }
657                        fgh++;
658                        CopyQ(L_QQ_Q_l_lower,L_QQ_Q_l_upper);
659                        CopyAng(RYPang_lower,RYPang_upper);
660                      }
661                    }else{
662                      if (dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000))<atime){
663                        upperqtime = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
664                        orbits.getPosition((double) (upperqtime - gltle->GetFromTime())/60., &eCi);
665                        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]);
666                      }
667                      else {
668                        lowerqtime = upperqtime;
669                        upperqtime = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
670                        orbits.getPosition((double) (upperqtime - gltle->GetFromTime())/60., &eCi);
671                        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]);
672                        mcreen = j3;
673                        mctren = ik;
674                        if(fgh==0){
675                          CopyQ(L_QQ_Q_l_lower,L_QQ_Q_l_upper);
676                          CopyAng(RYPang_lower,RYPang_upper);
677                          lowerqtime = atime-1;
678                        }
679                        oi=ui;
680                        goto closethisloop;
681                        //_0 = true;
682                      }
683                      fgh++;
684                      CopyQ(L_QQ_Q_l_lower,L_QQ_Q_l_upper);
685                      CopyAng(RYPang_lower,RYPang_upper);
686                      //_0 = true;
687                    };
688                    //cin>>grib;
689                  };
690                };
691              };
692            };
693          };
694        closethisloop:
695          //
696          if ( debug ) printf(" I am There \n");
697          //
698          if (((atime>(UInt_t)upperqtime)||(atime<(UInt_t)lowerqtime)) && neventsm>0 ){
699            if ( debug ) printf(" I am there \n");
700            //
701            lowerqtime = upperqtime;
702            UInt_t maxloop = 100000000;    
703            UInt_t mn = 0;
704            bool gh=false;
705            ooi=oi;
706            if ( verbose ) printf(" OrbitalInfoCore: sync with quaternions data upperqtime %u lowerqtime %u atime %u \n",(UInt_t)upperqtime,(UInt_t)lowerqtime,atime);
707            while (!gh){      
708              if ( mn > maxloop ){
709                if ( verbose ) printf(" OrbitalInfoCore: quaternions sync out of range! exiting\n");
710                gh = true;
711              };
712              mn++;
713              if (oi<5) oi++;
714              else oi=0;
715              if (oi==0){
716                mcreen++;
717                if (mcreen == numrec){
718                  mctren++;
719                  mcreen = 0;
720                  l0trm->GetEntry(mctren);
721                  numrec = mcmdev->Records->GetEntries();
722                }
723                CopyQ(L_QQ_Q_l_lower,L_QQ_Q_l_upper);
724                CopyAng(RYPang_lower,RYPang_upper);
725                mcmdrc = (pamela::McmdRecord*)mcmdev->Records->At(mcreen);
726                if ((int)mcmdrc->ID1 == 226){
727                  L_QQ_Q_l_upper->fill(mcmdrc->McmdData);
728                  upperqtime = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
729                  if (upperqtime<lowerqtime){
730                    upperqtime=runinfo->RUNTRAILER_TIME;
731                    CopyQ(L_QQ_Q_l_upper,L_QQ_Q_l_lower);
732                    CopyAng(RYPang_upper,RYPang_lower);
733                  }else{
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[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]);
736                  }
737                  //              re--;
738                  gh=true;
739                }
740              }else{
741                if ((Int_t)L_QQ_Q_l_upper->time[oi]>(Int_t)L_QQ_Q_l_upper->time[0]){
742                  upperqtime = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000));
743                  orbits.getPosition((double) (upperqtime - gltle->GetFromTime())/60., &eCi);
744                  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]);
745                  orbits.getPosition((double) (lowerqtime - gltle->GetFromTime())/60., &eCi);
746                  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]);
747                  //              re--;
748                  gh=true;
749                };
750              };
751            };
752            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);
753          };
754          //
755          if ( debug ) printf(" I am THIS \n");
756          //
757          // Fill in quaternions and angles
758          //
759          if ((atime<=(UInt_t)upperqtime)&&(atime>=(UInt_t)lowerqtime)&& neventsm>0){      
760            if ( debug ) printf(" I am this \n");
761            UInt_t tut = holeq(lowerqtime, upperqtime, L_QQ_Q_l_lower, L_QQ_Q_l_upper, oi);
762            if (oi == 0){
763              if ((tut!=5)||(tut!=6)){
764                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)));
765                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));
766                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)));
767                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));
768                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)));
769                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));
770                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)));
771                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));
772            
773                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)));
774                orbitalinfo->theta =  incli*atime+RYPang_upper->Tangazh-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
775                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)));
776                orbitalinfo->phi =  incli*atime+RYPang_upper->Ryskanie-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
777                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)));
778                orbitalinfo->etha =  incli*atime+RYPang_upper->Kren-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
779              }
780              if (tut==6){
781                if (fabs(RYPang_lower->Kren-RYPang_upper->Kren)<0.1){
782                  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)));
783                  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));
784                  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)));
785                  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));
786                  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)));
787                  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));
788                  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)));
789                  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));
790            
791                  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)));
792                  orbitalinfo->theta =  incli*atime+RYPang_upper->Tangazh-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
793                  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)));
794                  orbitalinfo->phi =  incli*atime+RYPang_upper->Ryskanie-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
795                  //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";
796                  //cin>>grib;
797                  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)));
798                  orbitalinfo->etha =  incli*atime+RYPang_upper->Kren-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
799                }
800              }
801            } else {
802              if((tut!=6)||(tut!=7)||(tut!=9)){
803                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)));
804                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));
805                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)));
806                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));
807                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)));
808                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));
809                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)));
810                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));
811            
812                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)));
813                orbitalinfo->theta =  incli*atime+RYPang_upper->Tangazh-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000));
814                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)));
815                orbitalinfo->phi =  incli*atime+RYPang_upper->Ryskanie-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000));
816                //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";
817                //cin>>grib;
818                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)));
819                orbitalinfo->etha =  incli*atime+RYPang_upper->Kren-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000));
820              }
821              if (tut==6){
822                if (fabs(RYPang_lower->Kren-RYPang_upper->Kren)<0.1){
823                  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)));
824                  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));
825                  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)));
826                  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));
827                  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)));
828                  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));
829                  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)));
830                  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));
831            
832                  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)));
833                  orbitalinfo->theta =  incli*atime+RYPang_upper->Tangazh-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000));
834                  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)));
835                  orbitalinfo->phi =  incli*atime+RYPang_upper->Ryskanie-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000));
836                  //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";
837                  //cin>>grib;
838                  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)));
839                  orbitalinfo->etha =  incli*atime+RYPang_upper->Kren-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000));
840                }
841              }            
842            }
843            //
844            orbitalinfo->mode = holeq(lowerqtime, upperqtime, L_QQ_Q_l_lower, L_QQ_Q_l_upper, oi);
845            //
846          } else {
847            if ( debug ) printf(" ops no incl! \n");
848            orbitalinfo->mode = -1;
849          };
850      
851          //
852          // fill orbital positions
853          //        
854        // Build coordinates in the right range.  We want to convert,        // Build coordinates in the right range.  We want to convert,
855        // longitude from (0, 2*pi) to (-180deg, 180deg).  Altitude is        // longitude from (0, 2*pi) to (-180deg, 180deg).  Altitude is
856        // in meters.        // in meters.
857        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);
858        lat = rad2deg(coo.m_Lat);        lat = rad2deg(coo.m_Lat);
859        alt = coo.m_Alt;        alt = coo.m_Alt;
860          //
   
       //      if((lon>180) || (lon<-180) || (lat>90) || (lat<-90) || (alt<0))  
       //        continue;  
861        if( lon<180 && lon>-180 && lat<90 && lat>-90 && alt>0 ){          if( lon<180 && lon>-180 && lat<90 && lat>-90 && alt>0 ){  
862            //      
863          orbitalinfo->lon = lon;          orbitalinfo->lon = lon;
864          orbitalinfo->lat = lat;          orbitalinfo->lat = lat;
865          orbitalinfo->alt = alt ;          orbitalinfo->alt = alt ;
866                    //
867          // compute mag field components and L shell.          // compute mag field components and L shell.
868            //
869          feldg_(&lat, &lon, &alt, &bnorth, &beast, &bdown, &babs);          feldg_(&lat, &lon, &alt, &bnorth, &beast, &bdown, &babs);
870          shellg_(&lat, &lon, &alt, &dimo, &xl, &icode, &bab1);          shellg_(&lat, &lon, &alt, &dimo, &xl, &icode, &bab1);
871          findb0_(&stps, &bdel, &value, &bequ, &rr0);          findb0_(&stps, &bdel, &value, &bequ, &rr0);
872                    //
873          orbitalinfo->Bnorth = bnorth;          orbitalinfo->Bnorth = bnorth;
874          orbitalinfo->Beast = beast;          orbitalinfo->Beast = beast;
875          orbitalinfo->Bdown = bdown;          orbitalinfo->Bdown = bdown;
876          orbitalinfo->Babs = babs;          orbitalinfo->Babs = babs;
877          orbitalinfo->BB0 = babs/bequ;          orbitalinfo->BB0 = babs/bequ;
878          orbitalinfo->L = xl;          orbitalinfo->L = xl;      
           
879          // Set Stormer vertical cutoff using L shell.          // Set Stormer vertical cutoff using L shell.
880          orbitalinfo->cutoff[0] = 14.9/(xl*xl);          orbitalinfo->cutoff[0] = 14.9/(xl*xl);
881        };          //
882                };      
883        // end EM        //
884          // Fill the class
885          //
886        OrbitalInfotr->Fill();        OrbitalInfotr->Fill();
887      }        //
888            }; // loop over the events in the run
889      //      //
890      // Here you may want to clear some variables before processing another run        // Here you may want to clear some variables before processing another run  
891      //      //
892      delete dbtime;      delete dbtime;
893        if ( L_QQ_Q_l_upper ) delete L_QQ_Q_l_upper;
894        if ( L_QQ_Q_l_lower ) delete L_QQ_Q_l_lower;
895        if ( RYPang_upper ) delete RYPang_upper;
896        if ( RYPang_lower ) delete RYPang_lower;
897    }; // process all the runs    }; // process all the runs
898    //    
899    if (verbose) printf("\n Finished processing data \n");    if (verbose) printf("\n Finished processing data \n");
900    //    //
901   closeandexit:   closeandexit:
# Line 663  cCoordGeo getCoo(UInt_t atime, UInt_t tl Line 970  cCoordGeo getCoo(UInt_t atime, UInt_t tl
970  {  {
971    cEci eci;    cEci eci;
972    cOrbit orbit(*tle);    cOrbit orbit(*tle);
     
973    orbit.getPosition((double) (atime - tletime)/60., &eci);    orbit.getPosition((double) (atime - tletime)/60., &eci);
974        
975    return eci.toGeo();    return eci.toGeo();
976  }  }
977    
978    // function of copyng of quatrnions classes
979    
980    void CopyQ(Quaternions *Q1, Quaternions *Q2){
981      for(UInt_t i = 0; i < 6; i++){
982        Q1->time[i]=Q2->time[i];
983        for (UInt_t j = 0; j < 4; j++)Q1->quat[i][j]=Q2->quat[i][j];
984      }
985      return;
986    }
987    
988    // functions of copyng InclinationInfo classes
989    
990    void CopyAng(InclinationInfo *A1, InclinationInfo *A2){
991      A1->Tangazh = A2->Tangazh;
992      A1->Ryskanie = A2->Ryskanie;
993      A1->Kren = A2->Kren;
994      return;
995    }
996    
997    UInt_t holeq(Double_t lower,Double_t upper,Quaternions *Qlower, Quaternions *Qupper, UInt_t f){
998      
999      UInt_t hole = 10;
1000      bool R10l = false;     // Sign of R10 mode in lower quaternions array
1001      bool R10u = false;     // Sign of R10 mode in upper quaternions array
1002      bool insm = false;     // Sign that we inside quaternions array
1003      bool mxtml = false;    // Sign of mixt mode in lower quaternions array
1004      bool mxtmu = false;    // Sign of mixt mode in upper quaternions array
1005      bool npasm = false;     // Sign of normall pass between R10 and non R10 or between non R10 and R10
1006      UInt_t NCQl = 6;       // Number of correct quaternions in lower array
1007      UInt_t NCQu = 6;       // Number of correct quaternions in upper array
1008      if (f>0){
1009        insm = true;
1010        if(Qupper->time[f]-Qupper->time[f-1]==30) R10u = false;
1011        if(Qupper->time[f]-Qupper->time[f-1]<1) R10u = true;
1012      }else{
1013        insm = false;
1014        if((Qlower->time[5]-Qlower->time[0]<2)&&(Qlower->time[1]-Qlower->time[0]<2)) R10l = true;
1015        if((Qupper->time[5]-Qupper->time[0]<2)&&(Qupper->time[1]-Qupper->time[0]<2)) R10u = true;
1016        if((Qlower->time[5]-Qlower->time[0]==150)&&(Qlower->time[1]-Qlower->time[0]==30)) R10l = false;
1017        if((Qupper->time[5]-Qupper->time[0]==150)&&(Qupper->time[1]-Qupper->time[0]==30)) R10u = false;
1018        if((Qlower->time[5]-Qlower->time[0]<2)&&(Qlower->time[1]-Qlower->time[0]==30)){
1019          mxtml = true;
1020          for(UInt_t i = 1; i < 6; i++){
1021            if(Qlower->time[i]-Qlower->time[0]==30*i) NCQl=i;
1022          }
1023        }
1024        if((Qupper->time[5]-Qupper->time[0]<2)&&(Qupper->time[1]-Qupper->time[0]==30)){
1025          mxtmu = true;
1026          for(UInt_t i = 1; i < 6; i++){
1027            if(Qupper->time[i]-Qupper->time[0]==30*i) NCQu=i;
1028          }
1029        }
1030      }
1031      
1032      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;
1033      
1034      
1035      if (R10u&&insm) hole=0; // best event R10
1036      if ((upper-lower<=5)&&(!insm)&&R10l&&R10u) hole = 1; // when first of 6 quaternions in array is correct
1037      if (((!R10u)&&insm)||((!insm)&&(!R10u)&&(!R10l)&&((upper-lower==210+(6-NCQl)*30)||(upper-lower==30)))) hole = 2; //non R10
1038      if (npasm&&(!insm)&&((R10l&&!R10u)||(R10u&&!R10l))) hole = 3; //normall pass from R10 to non R10 or from non R10 to R10
1039      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
1040      if ((upper-lower>=300)&&(!insm)&&((R10l&&!R10u)||(R10u&&!R10l))) hole = 5; //uneliminable hole between R10 and non R10 or between non R10 and R10
1041      if ((upper-lower>5)&&(upper-lower<=300)&&R10u&&R10l) hole = 6; // eliminable hole inside R10
1042      if ((upper-lower>300)&&R10u&&R10l) hole = 7; //uneliminable hole inside R10
1043      if ((upper-lower>210)&&(upper-lower<=1200)&&(!R10u)&&(!R10l)) hole = 8; //eliminable hole inside non R10
1044      if ((upper-lower>1200)&&!R10u&&!R10l) hole = 9; // uneliminable hole inside non R10
1045      return hole;
1046    }
1047    
Legend:
Removed from v.1.10  
changed lines
  Added in v.1.18
  ViewVC Help
Powered by ViewVC 1.1.23