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