/[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.21 by mocchiut, Mon May 7 04:33:52 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 470  int OrbitalInfoCore(UInt_t run, TFile *f Line 490  int OrbitalInfoCore(UInt_t run, TFile *f
490        code = -11;        code = -11;
491        goto closeandexit;        goto closeandexit;
492      };      };
493      //      //
494      if ( runinfo->EV_TO > nevents-1 ) {      if ( runinfo->EV_TO > nevents-1 ) {
495        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");
496        l0File->Close();        l0File->Close();
# 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        //        //
568        if ( (atime > runinfo->RUNTRAILER_TIME) || (atime < runinfo->RUNHEADER_TIME)  ) {        if ( (atime > runinfo->RUNTRAILER_TIME) || (atime < runinfo->RUNHEADER_TIME)  ) {
569          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");
570          debug = true;  //      debug = true;
571          continue;          continue;
572        }        }
573        //        //
# 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            Long64_t maxloop = 100000000LL;
703            Long64_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                neventsm = 0;
712              };
713              mn++;
714              if (oi<5) oi++;
715              else oi=0;
716              if (oi==0){
717                mcreen++;
718                if (mcreen == numrec){
719                  mctren++;
720                  mcreen = 0;
721                  l0trm->GetEntry(mctren);
722                  numrec = mcmdev->Records->GetEntries();
723                }
724                CopyQ(L_QQ_Q_l_lower,L_QQ_Q_l_upper);
725                CopyAng(RYPang_lower,RYPang_upper);
726                mcmdrc = (pamela::McmdRecord*)mcmdev->Records->At(mcreen);
727                if ((int)mcmdrc->ID1 == 226){
728                  L_QQ_Q_l_upper->fill(mcmdrc->McmdData);
729                  upperqtime = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
730                  if (upperqtime<lowerqtime){
731                    upperqtime=runinfo->RUNTRAILER_TIME;
732                    CopyQ(L_QQ_Q_l_upper,L_QQ_Q_l_lower);
733                    CopyAng(RYPang_upper,RYPang_lower);
734                  }else{
735                    orbits.getPosition((double) (upperqtime - gltle->GetFromTime())/60., &eCi);
736                    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]);
737                  }
738                  //              re--;
739                  gh=true;
740                }
741              }else{
742                if ((Int_t)L_QQ_Q_l_upper->time[oi]>(Int_t)L_QQ_Q_l_upper->time[0]){
743                  upperqtime = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000));
744                  orbits.getPosition((double) (upperqtime - gltle->GetFromTime())/60., &eCi);
745                  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]);
746                  orbits.getPosition((double) (lowerqtime - gltle->GetFromTime())/60., &eCi);
747                  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]);
748                  //              re--;
749                  gh=true;
750                };
751              };
752            };
753            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);
754          };
755          //
756          if ( debug ) printf(" I am THIS \n");
757          //
758          // Fill in quaternions and angles
759          //
760          if ((atime<=(UInt_t)upperqtime)&&(atime>=(UInt_t)lowerqtime)&& neventsm>0){      
761            if ( debug ) printf(" I am this \n");
762            UInt_t tut = holeq(lowerqtime, upperqtime, L_QQ_Q_l_lower, L_QQ_Q_l_upper, oi);
763            if (oi == 0){
764              if ((tut!=5)||(tut!=6)){
765                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)));
766                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));
767                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)));
768                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));
769                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)));
770                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));
771                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)));
772                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));
773            
774                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)));
775                orbitalinfo->theta =  incli*atime+RYPang_upper->Tangazh-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
776                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)));
777                orbitalinfo->phi =  incli*atime+RYPang_upper->Ryskanie-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
778                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)));
779                orbitalinfo->etha =  incli*atime+RYPang_upper->Kren-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
780              }
781              if (tut==6){
782                if (fabs(RYPang_lower->Kren-RYPang_upper->Kren)<0.1){
783                  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)));
784                  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));
785                  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)));
786                  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));
787                  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)));
788                  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));
789                  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)));
790                  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));
791            
792                  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)));
793                  orbitalinfo->theta =  incli*atime+RYPang_upper->Tangazh-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
794                  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)));
795                  orbitalinfo->phi =  incli*atime+RYPang_upper->Ryskanie-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
796                  //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";
797                  //cin>>grib;
798                  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)));
799                  orbitalinfo->etha =  incli*atime+RYPang_upper->Kren-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
800                }
801              }
802            } else {
803              if((tut!=6)||(tut!=7)||(tut!=9)){
804                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)));
805                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));
806                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)));
807                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));
808                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)));
809                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));
810                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)));
811                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));
812            
813                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)));
814                orbitalinfo->theta =  incli*atime+RYPang_upper->Tangazh-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000));
815                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)));
816                orbitalinfo->phi =  incli*atime+RYPang_upper->Ryskanie-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000));
817                //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";
818                //cin>>grib;
819                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)));
820                orbitalinfo->etha =  incli*atime+RYPang_upper->Kren-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000));
821              }
822              if (tut==6){
823                if (fabs(RYPang_lower->Kren-RYPang_upper->Kren)<0.1){
824                  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)));
825                  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));
826                  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)));
827                  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));
828                  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)));
829                  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));
830                  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)));
831                  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));
832            
833                  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)));
834                  orbitalinfo->theta =  incli*atime+RYPang_upper->Tangazh-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000));
835                  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)));
836                  orbitalinfo->phi =  incli*atime+RYPang_upper->Ryskanie-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000));
837                  //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";
838                  //cin>>grib;
839                  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)));
840                  orbitalinfo->etha =  incli*atime+RYPang_upper->Kren-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000));
841                }
842              }            
843            }
844            //
845            orbitalinfo->mode = holeq(lowerqtime, upperqtime, L_QQ_Q_l_lower, L_QQ_Q_l_upper, oi);
846            //
847          } else {
848            if ( debug ) printf(" ops no incl! \n");
849            orbitalinfo->mode = -1;
850          };
851      
852          //
853          // fill orbital positions
854          //        
855        // Build coordinates in the right range.  We want to convert,        // Build coordinates in the right range.  We want to convert,
856        // longitude from (0, 2*pi) to (-180deg, 180deg).  Altitude is        // longitude from (0, 2*pi) to (-180deg, 180deg).  Altitude is
857        // in meters.        // in meters.
858        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);
859        lat = rad2deg(coo.m_Lat);        lat = rad2deg(coo.m_Lat);
860        alt = coo.m_Alt;        alt = coo.m_Alt;
861          //
   
       //      if((lon>180) || (lon<-180) || (lat>90) || (lat<-90) || (alt<0))  
       //        continue;  
862        if( lon<180 && lon>-180 && lat<90 && lat>-90 && alt>0 ){          if( lon<180 && lon>-180 && lat<90 && lat>-90 && alt>0 ){  
863            //      
864          orbitalinfo->lon = lon;          orbitalinfo->lon = lon;
865          orbitalinfo->lat = lat;          orbitalinfo->lat = lat;
866          orbitalinfo->alt = alt ;          orbitalinfo->alt = alt ;
867                    //
868          // compute mag field components and L shell.          // compute mag field components and L shell.
869            //
870          feldg_(&lat, &lon, &alt, &bnorth, &beast, &bdown, &babs);          feldg_(&lat, &lon, &alt, &bnorth, &beast, &bdown, &babs);
871          shellg_(&lat, &lon, &alt, &dimo, &xl, &icode, &bab1);          shellg_(&lat, &lon, &alt, &dimo, &xl, &icode, &bab1);
872          findb0_(&stps, &bdel, &value, &bequ, &rr0);          findb0_(&stps, &bdel, &value, &bequ, &rr0);
873                    //
874          orbitalinfo->Bnorth = bnorth;          orbitalinfo->Bnorth = bnorth;
875          orbitalinfo->Beast = beast;          orbitalinfo->Beast = beast;
876          orbitalinfo->Bdown = bdown;          orbitalinfo->Bdown = bdown;
877          orbitalinfo->Babs = babs;          orbitalinfo->Babs = babs;
878          orbitalinfo->BB0 = babs/bequ;          orbitalinfo->BB0 = babs/bequ;
879          orbitalinfo->L = xl;          orbitalinfo->L = xl;      
           
880          // Set Stormer vertical cutoff using L shell.          // Set Stormer vertical cutoff using L shell.
881          orbitalinfo->cutoff[0] = 14.9/(xl*xl);          orbitalinfo->cutoff[0] = 14.9/(xl*xl);
882        };          //
883                };      
884        // end EM        //
885          // Fill the class
886          //
887        OrbitalInfotr->Fill();        OrbitalInfotr->Fill();
888      }        //
889            }; // loop over the events in the run
890      //      //
891      // Here you may want to clear some variables before processing another run        // Here you may want to clear some variables before processing another run  
892      //      //
893      delete dbtime;      delete dbtime;
894        if ( L_QQ_Q_l_upper ) delete L_QQ_Q_l_upper;
895        if ( L_QQ_Q_l_lower ) delete L_QQ_Q_l_lower;
896        if ( RYPang_upper ) delete RYPang_upper;
897        if ( RYPang_lower ) delete RYPang_lower;
898    }; // process all the runs    }; // process all the runs
899    //    
900    if (verbose) printf("\n Finished processing data \n");    if (verbose) printf("\n Finished processing data \n");
901    //    //
902   closeandexit:   closeandexit:
# Line 663  cCoordGeo getCoo(UInt_t atime, UInt_t tl Line 971  cCoordGeo getCoo(UInt_t atime, UInt_t tl
971  {  {
972    cEci eci;    cEci eci;
973    cOrbit orbit(*tle);    cOrbit orbit(*tle);
     
974    orbit.getPosition((double) (atime - tletime)/60., &eci);    orbit.getPosition((double) (atime - tletime)/60., &eci);
975        
976    return eci.toGeo();    return eci.toGeo();
977  }  }
978    
979    // function of copyng of quatrnions classes
980    
981    void CopyQ(Quaternions *Q1, Quaternions *Q2){
982      for(UInt_t i = 0; i < 6; i++){
983        Q1->time[i]=Q2->time[i];
984        for (UInt_t j = 0; j < 4; j++)Q1->quat[i][j]=Q2->quat[i][j];
985      }
986      return;
987    }
988    
989    // functions of copyng InclinationInfo classes
990    
991    void CopyAng(InclinationInfo *A1, InclinationInfo *A2){
992      A1->Tangazh = A2->Tangazh;
993      A1->Ryskanie = A2->Ryskanie;
994      A1->Kren = A2->Kren;
995      return;
996    }
997    
998    UInt_t holeq(Double_t lower,Double_t upper,Quaternions *Qlower, Quaternions *Qupper, UInt_t f){
999      
1000      UInt_t hole = 10;
1001      bool R10l = false;     // Sign of R10 mode in lower quaternions array
1002      bool R10u = false;     // Sign of R10 mode in upper quaternions array
1003      bool insm = false;     // Sign that we inside quaternions array
1004      bool mxtml = false;    // Sign of mixt mode in lower quaternions array
1005      bool mxtmu = false;    // Sign of mixt mode in upper quaternions array
1006      bool npasm = false;     // Sign of normall pass between R10 and non R10 or between non R10 and R10
1007      UInt_t NCQl = 6;       // Number of correct quaternions in lower array
1008      UInt_t NCQu = 6;       // Number of correct quaternions in upper array
1009      if (f>0){
1010        insm = true;
1011        if(Qupper->time[f]-Qupper->time[f-1]==30) R10u = false;
1012        if(Qupper->time[f]-Qupper->time[f-1]<1) R10u = true;
1013      }else{
1014        insm = false;
1015        if((Qlower->time[5]-Qlower->time[0]<2)&&(Qlower->time[1]-Qlower->time[0]<2)) R10l = true;
1016        if((Qupper->time[5]-Qupper->time[0]<2)&&(Qupper->time[1]-Qupper->time[0]<2)) R10u = true;
1017        if((Qlower->time[5]-Qlower->time[0]==150)&&(Qlower->time[1]-Qlower->time[0]==30)) R10l = false;
1018        if((Qupper->time[5]-Qupper->time[0]==150)&&(Qupper->time[1]-Qupper->time[0]==30)) R10u = false;
1019        if((Qlower->time[5]-Qlower->time[0]<2)&&(Qlower->time[1]-Qlower->time[0]==30)){
1020          mxtml = true;
1021          for(UInt_t i = 1; i < 6; i++){
1022            if(Qlower->time[i]-Qlower->time[0]==30*i) NCQl=i;
1023          }
1024        }
1025        if((Qupper->time[5]-Qupper->time[0]<2)&&(Qupper->time[1]-Qupper->time[0]==30)){
1026          mxtmu = true;
1027          for(UInt_t i = 1; i < 6; i++){
1028            if(Qupper->time[i]-Qupper->time[0]==30*i) NCQu=i;
1029          }
1030        }
1031      }
1032      
1033      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;
1034      
1035      
1036      if (R10u&&insm) hole=0; // best event R10
1037      if ((upper-lower<=5)&&(!insm)&&R10l&&R10u) hole = 1; // when first of 6 quaternions in array is correct
1038      if (((!R10u)&&insm)||((!insm)&&(!R10u)&&(!R10l)&&((upper-lower==210+(6-NCQl)*30)||(upper-lower==30)))) hole = 2; //non R10
1039      if (npasm&&(!insm)&&((R10l&&!R10u)||(R10u&&!R10l))) hole = 3; //normall pass from R10 to non R10 or from non R10 to R10
1040      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
1041      if ((upper-lower>=300)&&(!insm)&&((R10l&&!R10u)||(R10u&&!R10l))) hole = 5; //uneliminable hole between R10 and non R10 or between non R10 and R10
1042      if ((upper-lower>5)&&(upper-lower<=300)&&R10u&&R10l) hole = 6; // eliminable hole inside R10
1043      if ((upper-lower>300)&&R10u&&R10l) hole = 7; //uneliminable hole inside R10
1044      if ((upper-lower>210)&&(upper-lower<=1200)&&(!R10u)&&(!R10l)) hole = 8; //eliminable hole inside non R10
1045      if ((upper-lower>1200)&&!R10u&&!R10l) hole = 9; // uneliminable hole inside non R10
1046      return hole;
1047    }
1048    
Legend:
Removed from v.1.10  
changed lines
  Added in v.1.21
  ViewVC Help
Powered by ViewVC 1.1.23