/[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.66 by mocchiut, Mon Feb 3 16:45:14 2014 UTC revision 1.72 by mocchiut, Tue Apr 15 15:48:17 2014 UTC
# Line 48  Line 48 
48  #include <OrbitalInfoCore.h>  #include <OrbitalInfoCore.h>
49  #include <InclinationInfo.h>  #include <InclinationInfo.h>
50    
51    //
52    // Tracker and ToF classes headers and definitions
53    //
54    #include <ToFLevel2.h>
55    #include <TrkLevel2.h>
56    
57  using namespace std;  using namespace std;
58    
# Line 239  int OrbitalInfoCore(UInt_t run, TFile *f Line 244  int OrbitalInfoCore(UInt_t run, TFile *f
244    Int_t ltp3 = 0;    Int_t ltp3 = 0;
245    //  Int_t uno = 1;    //  Int_t uno = 1;
246    //  const char *niente = " ";    //  const char *niente = " ";
247      GL_PARAM *glparam0 = new GL_PARAM();
248    GL_PARAM *glparam = new GL_PARAM();    GL_PARAM *glparam = new GL_PARAM();
249    GL_PARAM *glparam2 = new GL_PARAM();    GL_PARAM *glparam2 = new GL_PARAM();
250    GL_PARAM *glparam3 = new GL_PARAM();    GL_PARAM *glparam3 = new GL_PARAM();
# Line 257  int OrbitalInfoCore(UInt_t run, TFile *f Line 263  int OrbitalInfoCore(UInt_t run, TFile *f
263    Double_t Pz = 0;      Double_t Pz = 0;  
264    TTree *ttof = 0;    TTree *ttof = 0;
265    ToFLevel2 *tof = new ToFLevel2();    ToFLevel2 *tof = new ToFLevel2();
266      TTree *ttrke = 0;
267      TrkLevel2 *trke = new TrkLevel2();
268    OrientationInfo *PO = new OrientationInfo();    OrientationInfo *PO = new OrientationInfo();
269    Int_t nz = 6;    Int_t nz = 6;
270    Float_t zin[6];    Float_t zin[6];
271    Int_t nevtofl2 = 0;    Int_t nevtofl2 = 0;
272      Int_t nevtrkl2 = 0;
273    if ( verbose ) cout<<"Reading quaternions external file"<<endl;    if ( verbose ) cout<<"Reading quaternions external file"<<endl;
274    cout.setf(ios::fixed,ios::floatfield);      cout.setf(ios::fixed,ios::floatfield);  
275    /******Reading recovered quaternions...*********/    /******Reading recovered quaternions...*********/
# Line 270  int OrbitalInfoCore(UInt_t run, TFile *f Line 279  int OrbitalInfoCore(UInt_t run, TFile *f
279    vector<Float_t> recq2;    vector<Float_t> recq2;
280    vector<Float_t> recq3;    vector<Float_t> recq3;
281    Float_t Norm = 1;    Float_t Norm = 1;
282    Int_t parerror=glparam->Query_GL_PARAM(1,303,dbc); // parameters stored in DB in GL_PRAM table      Int_t parerror=glparam0->Query_GL_PARAM(1,303,dbc); // parameters stored in DB in GL_PRAM table
283    ifstream in((glparam->PATH+glparam->NAME).Data(),ios::in);    if ( verbose ) cout<<parerror<<"\t"<<(char*)(glparam0->PATH+glparam0->NAME).Data()<<endl;
284      ifstream in((char*)(glparam0->PATH+glparam0->NAME).Data(),ios::in);
285    if ( parerror<0 ) {    if ( parerror<0 ) {
286      code = parerror;      code = parerror;
287      goto closeandexit;      //goto closeandexit;
288    }    }
289    while(!in.eof()){    while(!in.eof()){
290      recqtime.resize(recqtime.size()+1);      recqtime.resize(recqtime.size()+1);
# Line 290  int OrbitalInfoCore(UInt_t run, TFile *f Line 300  int OrbitalInfoCore(UInt_t run, TFile *f
300      in>>recq3[sizee-1];      in>>recq3[sizee-1];
301      in>>Norm;      in>>Norm;
302    }    }
303      in.close();
304    if ( verbose ) cout<<"We have read recovered data"<<endl;    if ( verbose ) cout<<"We have read recovered data"<<endl;
305      if (debug) cout << "size of recovered quaterions data set is " << recqtime.size() << endl;
306    
307      vector<UInt_t> RTtime1;
308      vector<UInt_t> RTtime2;
309      vector<Double_t> RTbank1;
310      vector<Double_t> RTbank2;
311      vector<Int_t> RTazim;
312      vector<Int_t> RTdir1;
313      vector<Int_t> RTdir2;
314      vector<Int_t> RTerrq;
315    
316    // 10RED CHECK
317    
318    //  TH2F* DIFFX = new TH2F("diffx","",100,0,100,90,0,90);
319    //  TH2F* DIFFY = new TH2F("diffy","",100,0,100,90,0,90);
320    //  TH2F* DIFFZ = new TH2F("diffz","",100,0,100,90,0,90);
321    
322      ofstream mc;
323      TString gr = "methodscomparison_";
324      gr+=run;
325      gr+=".txt";
326      mc.open(gr);
327      mc.setf(ios::fixed,ios::floatfield);
328    // 10RED CHECK END
329    
330      if ( verbose ) cout<<"read Rotation Table"<<endl;
331    
332      Int_t parerror2=glparam0->Query_GL_PARAM(1,305,dbc);
333      ifstream an((char*)(glparam0->PATH+glparam0->NAME).Data(),ios::in);
334      if ( verbose ) cout<<parerror2<<"\t"<<(char*)(glparam0->PATH+glparam0->NAME).Data()<<endl;
335    //  if ( parerror2<0 ) {
336    //    code = parerror;
337        //goto closeandexit;
338    //  }
339    
340    //ifstream an("/data03/Malakhov/pam9Malakhov/installed10/calib/orb-param/RDBCC.txt",ios::in);
341      while(!an.eof()){
342        RTtime1.resize(RTtime1.size()+1);
343        Int_t sizee = RTtime1.size();
344        RTbank1.resize(sizee+1);
345        RTazim.resize(sizee+1);
346        RTdir1.resize(sizee+1);
347        RTerrq.resize(sizee+1);
348        an>>RTtime1[sizee-1];
349        an>>RTbank1[sizee-1];
350        an>>RTazim[sizee-1];
351        an>>RTdir1[sizee-1];
352        an>>RTerrq[sizee-1];
353        if(sizee>1) {
354          RTtime2.resize(sizee+1);
355          RTbank2.resize(sizee+1);
356          RTdir2.resize(sizee+1);
357          RTtime2[sizee-2]=RTtime1[sizee-1];
358          RTbank2[sizee-2]=RTbank1[sizee-1];
359          RTdir2[sizee-2]=RTdir1[sizee-1];
360        }
361      }
362      an.close();
363      //cout<<"put some number here"<<endl;
364      //Int_t yupi;
365      //cin>>yupi;
366      
367      if ( verbose ) cout<<"We have read Rotation Table"<<endl;
368        //Geomagnetic coordinates calculations staff
369    
370      GMtype_CoordGeodetic location;
371      //  GMtype_CoordDipole GMlocation;
372      GMtype_Ellipsoid Ellip;
373      GMtype_Data G0, G1, H1;
374            
375      //  { // this braces is necessary to avoid jump to label 'closeandexit'  error   // but it is wrong since the variable "igpath" will not exist outside. To overcome the "jump to label 'closeandexit'  error" it is necessary to set the "igpath" before line 276
376      //    TString igpath="/data03/Malakhov/pam9Malakhov/installed10/calib/orb-param/";
377      //  }
378    
379      //  GM_ScanIGRF(glparam->PATH, &G0, &G1, &H1);
380      GM_ScanIGRF(dbc, &G0, &G1, &H1);
381    
382      //cout << G0.element[0] << "\t" << G1.element[0] << "\t" << H1.element[0] << endl;
383      //cout << G0.element[5] << "\t" << G1.element[5] << "\t" << H1.element[5] << endl;
384    
385      GM_SetEllipsoid(&Ellip);
386    
387    // IGRF stuff moved inside run loop!      // IGRF stuff moved inside run loop!  
388    
# Line 307  int OrbitalInfoCore(UInt_t run, TFile *f Line 399  int OrbitalInfoCore(UInt_t run, TFile *f
399        if ( verbose ) printf(" OrbitalInfo - ERROR: no tof tree\n");        if ( verbose ) printf(" OrbitalInfo - ERROR: no tof tree\n");
400        code = -900;        code = -900;
401        goto closeandexit;        goto closeandexit;
402      };      }
403      ttof->SetBranchAddress("ToFLevel2",&tof);        ttof->SetBranchAddress("ToFLevel2",&tof);  
404      nevtofl2 = ttof->GetEntries();      nevtofl2 = ttof->GetEntries();
405    };  
406        ttrke = (TTree*)file->Get("Tracker");
407        if ( !ttrke ) {
408          if ( verbose ) printf(" OrbitalInfo - ERROR: no trk tree\n");
409          code = -903;
410          goto closeandexit;
411        }
412        ttrke->SetBranchAddress("TrkLevel2",&trke);  
413        nevtrkl2 = ttrke->GetEntries();
414      }
415    //    //
416    // Let's start!    // Let's start!
417    //    //
# Line 591  int OrbitalInfoCore(UInt_t run, TFile *f Line 692  int OrbitalInfoCore(UInt_t run, TFile *f
692      //      //
693      // open IGRF files and do it only once if we are processing a full level2 file      // open IGRF files and do it only once if we are processing a full level2 file
694      //      //
695      if ( reprocall && !igrfloaded ){      if ( !igrfloaded ){
696    
697        if ( l0head->GetEntry(runinfo->EV_FROM) > 0 ){        if ( l0head->GetEntry(runinfo->EV_FROM) > 0 ){
698          igrfloaded = true;          igrfloaded = true;
699        //          //
700        // absolute time of first event of the run (it should not matter a lot)          // absolute time of first event of the run (it should not matter a lot)
701        //          //
702        ph = eh->GetPscuHeader();          ph = eh->GetPscuHeader();
703        atime = dbtime->DBabsTime(ph->GetOrbitalTime());          atime = dbtime->DBabsTime(ph->GetOrbitalTime());
704                  
705        parerror=glparam->Query_GL_PARAM(atime-anni5,301,dbc); // parameters stored in DB in GL_PRAM table            parerror=glparam->Query_GL_PARAM(atime-anni5,301,dbc); // parameters stored in DB in GL_PRAM table  
706        if ( parerror<0 ) {          if ( parerror<0 ) {
707          code = parerror;            code = parerror;
708          goto closeandexit;            goto closeandexit;
709        }          }
710        ltp1 = (Int_t)(glparam->PATH+glparam->NAME).Length();          ltp1 = (Int_t)(glparam->PATH+glparam->NAME).Length();
711        if ( verbose ) printf(" Reading Earth's Magnetic Field parameter file: %s \n",(glparam->PATH+glparam->NAME).Data());          if ( verbose ) printf(" Reading Earth's Magnetic Field parameter file: %s \n",(glparam->PATH+glparam->NAME).Data());
712        //          //
713        parerror=glparam2->Query_GL_PARAM(atime,301,dbc); // parameters stored in DB in GL_PRAM table            parerror=glparam2->Query_GL_PARAM(atime,301,dbc); // parameters stored in DB in GL_PRAM table  
714        if ( parerror<0 ) {          if ( parerror<0 ) {
715          code = parerror;            code = parerror;
716          goto closeandexit;            goto closeandexit;
717        }          }
718        ltp2 = (Int_t)(glparam2->PATH+glparam2->NAME).Length();          ltp2 = (Int_t)(glparam2->PATH+glparam2->NAME).Length();
719        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());
720        //          //
721        parerror=glparam3->Query_GL_PARAM(atime,302,dbc); // parameters stored in DB in GL_PRAM table          parerror=glparam3->Query_GL_PARAM(atime,302,dbc); // parameters stored in DB in GL_PRAM table
722        if ( parerror<0 ) {          if ( parerror<0 ) {
723          code = parerror;            code = parerror;
724          goto closeandexit;            goto closeandexit;
725        }          }
726        ltp3 = (Int_t)(glparam3->PATH+glparam3->NAME).Length();          ltp3 = (Int_t)(glparam3->PATH+glparam3->NAME).Length();
727        if ( verbose ) printf(" Reading Earth's Magnetic Field parameter file: %s \n",(glparam3->PATH+glparam3->NAME).Data());          if ( verbose ) printf(" Reading Earth's Magnetic Field parameter file: %s \n",(glparam3->PATH+glparam3->NAME).Data());
728        //          //
729        initize_((char *)(glparam->PATH+glparam->NAME).Data(),&ltp1,(char *)(glparam2->PATH+glparam2->NAME).Data(),&ltp2,(char *)(glparam3->PATH+glparam3->NAME).Data(),&ltp3);          initize_((char *)(glparam->PATH+glparam->NAME).Data(),&ltp1,(char *)(glparam2->PATH+glparam2->NAME).Data(),&ltp2,(char *)(glparam3->PATH+glparam3->NAME).Data(),&ltp3);
730        //          //
731            if (debug) cout<<"initize: "<<(char *)(glparam->PATH+glparam->NAME).Data()<<"\t"<<(char *)(glparam2->PATH+glparam2->NAME).Data()<<"\t"<<(char *)(glparam3->PATH+glparam3->NAME).Data()<<endl;
732        }        }
733      }      }
734      //      //
# Line 781  int OrbitalInfoCore(UInt_t run, TFile *f Line 883  int OrbitalInfoCore(UInt_t run, TFile *f
883    
884      Int_t nt = 0;      Int_t nt = 0;
885            
     //init sine-function interpolation  
       
     //cout<<"Sine coeficient initialisation..."<<endl;  
     vector<Sine> q0sine;  
     vector<Sine> q1sine;  
     vector<Sine> q2sine;  
     vector<Sine> q3sine;  
     vector<Sine> Yawsine;  
   
     /*TH2F* q0testing = new TH2F();  
       TH2F* q1testing = new TH2F();  
       TH2F* q2testing = new TH2F();  
       TH2F* q3testing = new TH2F();  
       TH2F* Pitchtesting = new TH2F();  
     */  
886      UInt_t must = 0;      UInt_t must = 0;
887    
888      //      //
# Line 841  int OrbitalInfoCore(UInt_t run, TFile *f Line 928  int OrbitalInfoCore(UInt_t run, TFile *f
928            if ( verbose ) printf(" OrbitalInfo - ERROR: no tof events with entry = %i in Level2 file\n",itr);            if ( verbose ) printf(" OrbitalInfo - ERROR: no tof events with entry = %i in Level2 file\n",itr);
929            if ( debug ) printf(" nobefrun %u re %u evfrom %u jumped %u reprocall %i \n",nobefrun,re,evfrom,jumped,reprocall);            if ( debug ) printf(" nobefrun %u re %u evfrom %u jumped %u reprocall %i \n",nobefrun,re,evfrom,jumped,reprocall);
930            l0File->Close();            l0File->Close();
931            code = -901;            code = -904;
932            goto closeandexit;            goto closeandexit;
933          };          };
934          //          //
935          tof->Clear();          tof->Clear();
936          //          //
937          if ( ttof->GetEntry(itr) <= 0 ) throw -36;          if ( ttof->GetEntry(itr) <= 0 ){
938             if ( verbose ) printf(" problems with tof tree entries... entry = %i in Level2 file\n",itr);
939             if ( verbose ) printf(" nobefrun %u re %u evfrom %u jumped %u reprocall %i \n",nobefrun,re,evfrom,jumped,reprocall);
940             throw -36;
941            }
942          //          //
943        };        }
944          //
945          // retrieve tracker informations
946          //
947          if ( !standalone ){
948            if ( itr > nevtrkl2 ){  
949              if ( verbose ) printf(" OrbitalInfo - ERROR: no trk events with entry = %i in Level2 file\n",itr);
950              if ( debug ) printf(" nobefrun %u re %u evfrom %u jumped %u reprocall %i \n",nobefrun,re,evfrom,jumped,reprocall);
951              l0File->Close();
952              code = -905;
953              goto closeandexit;
954            };
955            //
956            trke->Clear();
957            //
958            if ( ttrke->GetEntry(itr) <= 0 ) throw -36;
959            //
960          }
961    
962    
963        //        //
964        procev++;        procev++;
965        //        //
# Line 861  int OrbitalInfoCore(UInt_t run, TFile *f Line 971  int OrbitalInfoCore(UInt_t run, TFile *f
971        OrbitalInfoTrkVar *t_orb = new OrbitalInfoTrkVar();        OrbitalInfoTrkVar *t_orb = new OrbitalInfoTrkVar();
972        if( !(orbitalinfo->OrbitalInfoTrk) ) orbitalinfo->OrbitalInfoTrk = new TClonesArray("OrbitalInfoTrkVar",2);        if( !(orbitalinfo->OrbitalInfoTrk) ) orbitalinfo->OrbitalInfoTrk = new TClonesArray("OrbitalInfoTrkVar",2);
973        TClonesArray &tor = *orbitalinfo->OrbitalInfoTrk;        TClonesArray &tor = *orbitalinfo->OrbitalInfoTrk;
974    
975          // Geomagnetic coordinates calculation variables
976          GMtype_CoordSpherical CoordSpherical, DipoleSpherical;
977          GMtype_CoordCartesian CoordCartesian, DipoleCartesian;
978          GMtype_Model Model;
979          GMtype_Pole Pole;
980    
981        //        //
982        // Fill OBT, pkt_num and absTime        // Fill OBT, pkt_num and absTime
983        //              //      
# Line 890  int OrbitalInfoCore(UInt_t run, TFile *f Line 1007  int OrbitalInfoCore(UInt_t run, TFile *f
1007              + (month*31.+ (float) day)/365.              + (month*31.+ (float) day)/365.
1008              + (hour*3600.+min*60.+(float)sec)/(24.*3600.*365.);              + (hour*3600.+min*60.+(float)sec)/(24.*3600.*365.);
1009            //            //
1010            if ( debug ) printf(" %i compute magnetic dipole moment get dipole moment for year\n",procev);                  if ( debug ) printf(" %i compute magnetic dipole moment get dipole moment for year\n",procev);            
1011            if ( debug ) printf(" %i jyear %f dimo %f \n",procev,jyear,dimo);                  if ( debug ) printf(" %i jyear %f dimo %f \n",procev,jyear,dimo);
1012            feldcof_(&jyear, &dimo); // get dipole moment for year            feldcof_(&jyear, &dimo); // get dipole moment for year
1013            if ( debug ) printf(" %i compute magnetic dipole moment end\n",procev);                  if ( debug ) printf(" %i compute magnetic dipole moment end\n",procev);
1014    
1015              GM_TimeAdjustCoefs(year, jyear, G0, G1, H1, &Model);
1016              GM_PoleLocation(Model, &Pole);
1017              
1018          } else {          } else {
1019            code = -56;            code = -56;
1020            goto closeandexit;            goto closeandexit;
# Line 903  int OrbitalInfoCore(UInt_t run, TFile *f Line 1024  int OrbitalInfoCore(UInt_t run, TFile *f
1024        //        //
1025        cOrbit orbits(*gltle->GetTle());        cOrbit orbits(*gltle->GetTle());
1026        //        //
       if ( debug ) printf(" I am Here \n");  
       //  
1027        // synchronize with quaternions data        // synchronize with quaternions data
1028        //        //
1029        if ( isf && neventsm>0 ){        if ( isf && neventsm>0 ){
# Line 918  int OrbitalInfoCore(UInt_t run, TFile *f Line 1037  int OrbitalInfoCore(UInt_t run, TFile *f
1037            if ( ch->GetEntry(ik) <= 0 ) throw -36;            if ( ch->GetEntry(ik) <= 0 ) throw -36;
1038            tmpSize = mcmdev->Records->GetEntries();            tmpSize = mcmdev->Records->GetEntries();
1039            //      numrec = tmpSize;            //      numrec = tmpSize;
1040              if ( debug ) cout << "packet number " << ik <<"\tnumber of subpackets is " << tmpSize << endl;
1041            for (Int_t j3 = 0;j3<tmpSize;j3++){  //number of subpackets            for (Int_t j3 = 0;j3<tmpSize;j3++){  //number of subpackets
             if ( debug ) printf(" ik %i j3 %i eh eh eh \n",ik,j3);  
1042              mcmdrc = (pamela::McmdRecord*)mcmdev->Records->At(j3);              mcmdrc = (pamela::McmdRecord*)mcmdev->Records->At(j3);
1043              if ( mcmdrc ){ // missing inclination bug [8RED 090116]              if ( mcmdrc ){ // missing inclination bug [8RED 090116]
1044                if ( debug ) printf(" pluto \n");                if ( debug ) printf(" pluto \n");
1045                if ((int)mcmdrc->ID1 == 226 && mcmdrc->Mcmd_Block_crc_ok == 1){ //Check that it is Inclination Packet                if ((int)mcmdrc->ID1 == 226 && mcmdrc->Mcmd_Block_crc_ok == 1){ //Check that it is Inclination Packet
1046                  L_QQ_Q_l_upper->fill(mcmdrc->McmdData);                 L_QQ_Q_l_upper->fill(mcmdrc->McmdData);
1047                  for (UInt_t ui = 0; ui < 6; ui++){                  for (UInt_t ui = 0; ui < 6; ui++){
1048                    if (ui>0){                    if (ui>0){
1049                      if (L_QQ_Q_l_upper->time[ui]>L_QQ_Q_l_upper->time[0]){                      if (L_QQ_Q_l_upper->time[ui]>L_QQ_Q_l_upper->time[0]){
# Line 932  int OrbitalInfoCore(UInt_t run, TFile *f Line 1051  int OrbitalInfoCore(UInt_t run, TFile *f
1051                        Double_t u_time = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[ui]*1000-DeltaOBT*1000));                        Double_t u_time = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[ui]*1000-DeltaOBT*1000));
1052                        Int_t recSize = recqtime.size();                        Int_t recSize = recqtime.size();
1053                        if(lowerqtime > recqtime[recSize-1]){                        if(lowerqtime > recqtime[recSize-1]){
1054                          Int_t sizeqmcmd = qtime.size();                           // to avoid interpolation between bad quaternions arrays
1055                          inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);                           if(sqrt(pow(L_QQ_Q_l_upper->quat[ui][0],2)+pow(L_QQ_Q_l_upper->quat[ui][1],2)+pow(L_QQ_Q_l_upper->quat[ui][2],2)+pow(L_QQ_Q_l_upper->quat[ui][3],2))>0.99999){
                         qtime[sizeqmcmd]=u_time;  
                         q0[sizeqmcmd]=L_QQ_Q_l_upper->quat[ui][0];  
                         q1[sizeqmcmd]=L_QQ_Q_l_upper->quat[ui][1];  
                         q2[sizeqmcmd]=L_QQ_Q_l_upper->quat[ui][2];  
                         q3[sizeqmcmd]=L_QQ_Q_l_upper->quat[ui][3];  
                         qmode[sizeqmcmd]=holeq(lowerqtime,qtime[sizeqmcmd],L_QQ_Q_l_lower,L_QQ_Q_l_upper,ui);  
                         lowerqtime = u_time;  
                         orbits.getPosition((double) (u_time - gltle->GetFromTime())/60., &eCi);  
                         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]);  
                         qRoll[sizeqmcmd]=RYPang_upper->Kren;  
                         qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;  
                         qPitch[sizeqmcmd]=RYPang_upper->Tangazh;  
                       }  
                       for(Int_t mu = nt;mu<recSize;mu++){  
                         if(recqtime[mu]>lowerqtime && recqtime[mu]<u_time){  
                           nt=mu;  
                           Int_t sizeqmcmd = qtime.size();  
                           inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);  
                           qtime[sizeqmcmd]=recqtime[mu];  
                           q0[sizeqmcmd]=recq0[mu];  
                           q1[sizeqmcmd]=recq1[mu];  
                           q2[sizeqmcmd]=recq2[mu];  
                           q3[sizeqmcmd]=recq3[mu];  
                           qmode[sizeqmcmd]=-10;  
                           orbits.getPosition((double) (qtime[sizeqmcmd] - gltle->GetFromTime())/60., &eCi);  
                           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,recq0[mu],recq1[mu],recq2[mu],recq3[mu]);  
                           qRoll[sizeqmcmd]=RYPang_upper->Kren;  
                           qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;  
                           qPitch[sizeqmcmd]=RYPang_upper->Tangazh;  
                         }  
                         if(recqtime[mu]>=u_time){  
1056                            Int_t sizeqmcmd = qtime.size();                            Int_t sizeqmcmd = qtime.size();
1057                            inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);                            inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);
1058                            qtime[sizeqmcmd]=u_time;                            qtime[sizeqmcmd]=u_time;
# Line 979  int OrbitalInfoCore(UInt_t run, TFile *f Line 1067  int OrbitalInfoCore(UInt_t run, TFile *f
1067                            qRoll[sizeqmcmd]=RYPang_upper->Kren;                            qRoll[sizeqmcmd]=RYPang_upper->Kren;
1068                            qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;                            qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;
1069                            qPitch[sizeqmcmd]=RYPang_upper->Tangazh;                            qPitch[sizeqmcmd]=RYPang_upper->Tangazh;
1070                            break;                           }
1071                          }
1072                          for(Int_t mu = nt;mu<recSize;mu++){
1073                            if(recqtime[mu]>lowerqtime && recqtime[mu]<u_time){
1074                              if(sqrt(pow(recq0[mu],2)+pow(recq1[mu],2)+pow(recq2[mu],2)+pow(recq3[mu],2))>0.99999){
1075                                nt=mu;
1076                                Int_t sizeqmcmd = qtime.size();
1077                                inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);
1078                                qtime[sizeqmcmd]=recqtime[mu];
1079                                q0[sizeqmcmd]=recq0[mu];
1080                                q1[sizeqmcmd]=recq1[mu];
1081                                q2[sizeqmcmd]=recq2[mu];
1082                                q3[sizeqmcmd]=recq3[mu];
1083                                qmode[sizeqmcmd]=-10;
1084                                orbits.getPosition((double) (qtime[sizeqmcmd] - gltle->GetFromTime())/60., &eCi);
1085                                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,recq0[mu],recq1[mu],recq2[mu],recq3[mu]);
1086                                qRoll[sizeqmcmd]=RYPang_upper->Kren;
1087                                qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;
1088                                qPitch[sizeqmcmd]=RYPang_upper->Tangazh;
1089                              }
1090                            }
1091                            if(recqtime[mu]>=u_time){
1092                              if(sqrt(pow(L_QQ_Q_l_upper->quat[ui][0],2)+pow(L_QQ_Q_l_upper->quat[ui][1],2)+pow(L_QQ_Q_l_upper->quat[ui][2],2)+pow(L_QQ_Q_l_upper->quat[ui][3],2))>0.99999){
1093                                Int_t sizeqmcmd = qtime.size();
1094                                inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);
1095                                qtime[sizeqmcmd]=u_time;
1096                                q0[sizeqmcmd]=L_QQ_Q_l_upper->quat[ui][0];
1097                                q1[sizeqmcmd]=L_QQ_Q_l_upper->quat[ui][1];
1098                                q2[sizeqmcmd]=L_QQ_Q_l_upper->quat[ui][2];
1099                                q3[sizeqmcmd]=L_QQ_Q_l_upper->quat[ui][3];
1100                                qmode[sizeqmcmd]=holeq(lowerqtime,qtime[sizeqmcmd],L_QQ_Q_l_lower,L_QQ_Q_l_upper,ui);
1101                                lowerqtime = u_time;
1102                                orbits.getPosition((double) (u_time - gltle->GetFromTime())/60., &eCi);
1103                                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]);
1104                                qRoll[sizeqmcmd]=RYPang_upper->Kren;
1105                                qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;
1106                                qPitch[sizeqmcmd]=RYPang_upper->Tangazh;
1107                                break;
1108                              }
1109                          }                          }
1110                        }                        }
1111                      }                      }
# Line 989  int OrbitalInfoCore(UInt_t run, TFile *f Line 1115  int OrbitalInfoCore(UInt_t run, TFile *f
1115                      if(lowerqtime>u_time)nt=0;                      if(lowerqtime>u_time)nt=0;
1116                      Int_t recSize = recqtime.size();                      Int_t recSize = recqtime.size();
1117                      if(lowerqtime > recqtime[recSize-1]){                      if(lowerqtime > recqtime[recSize-1]){
1118                        Int_t sizeqmcmd = qtime.size();                        if(sqrt(pow(L_QQ_Q_l_upper->quat[ui][0],2)+pow(L_QQ_Q_l_upper->quat[ui][1],2)+pow(L_QQ_Q_l_upper->quat[ui][2],2)+pow(L_QQ_Q_l_upper->quat[ui][3],2))>0.99999){
                       inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);  
                       qtime[sizeqmcmd]=u_time;  
                       q0[sizeqmcmd]=L_QQ_Q_l_upper->quat[0][0];  
                       q1[sizeqmcmd]=L_QQ_Q_l_upper->quat[0][1];  
                       q2[sizeqmcmd]=L_QQ_Q_l_upper->quat[0][2];  
                       q3[sizeqmcmd]=L_QQ_Q_l_upper->quat[0][3];  
                       qmode[sizeqmcmd]=holeq(lowerqtime,qtime[sizeqmcmd],L_QQ_Q_l_lower,L_QQ_Q_l_upper,ui);  
                       lowerqtime = u_time;  
                       orbits.getPosition((double) (u_time - gltle->GetFromTime())/60., &eCi);  
                       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]);  
                       qRoll[sizeqmcmd]=RYPang_upper->Kren;  
                       qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;  
                       qPitch[sizeqmcmd]=RYPang_upper->Tangazh;  
                     }  
                     for(Int_t mu = nt;mu<recSize;mu++){  
                       if(recqtime[mu]>lowerqtime && recqtime[mu]<u_time){  
                         nt=mu;  
                         Int_t sizeqmcmd = qtime.size();  
                         inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);  
                         qtime[sizeqmcmd]=recqtime[mu];  
                         q0[sizeqmcmd]=recq0[mu];  
                         q1[sizeqmcmd]=recq1[mu];  
                         q2[sizeqmcmd]=recq2[mu];  
                         q3[sizeqmcmd]=recq3[mu];  
                         qmode[sizeqmcmd]=-10;  
                         orbits.getPosition((double) (qtime[sizeqmcmd] - gltle->GetFromTime())/60., &eCi);  
                         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,recq0[mu],recq1[mu],recq2[mu],recq3[mu]);  
                         qRoll[sizeqmcmd]=RYPang_upper->Kren;  
                         qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;  
                         qPitch[sizeqmcmd]=RYPang_upper->Tangazh;  
                       }  
                       if(recqtime[mu]>=u_time){  
1119                          Int_t sizeqmcmd = qtime.size();                          Int_t sizeqmcmd = qtime.size();
1120                          inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);                          inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);
1121                          qtime[sizeqmcmd]=u_time;                          qtime[sizeqmcmd]=u_time;
# Line 1036  int OrbitalInfoCore(UInt_t run, TFile *f Line 1130  int OrbitalInfoCore(UInt_t run, TFile *f
1130                          qRoll[sizeqmcmd]=RYPang_upper->Kren;                          qRoll[sizeqmcmd]=RYPang_upper->Kren;
1131                          qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;                          qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;
1132                          qPitch[sizeqmcmd]=RYPang_upper->Tangazh;                          qPitch[sizeqmcmd]=RYPang_upper->Tangazh;
1133                          CopyQ(L_QQ_Q_l_lower,L_QQ_Q_l_upper);                        }
1134                          break;                      }
1135                        for(Int_t mu = nt;mu<recSize;mu++){
1136                          if(recqtime[mu]>lowerqtime && recqtime[mu]<u_time){
1137                             if(sqrt(pow(recq0[mu],2)+pow(recq1[mu],2)+pow(recq2[mu],2)+pow(recq3[mu],2))>0.99999){
1138                               nt=mu;
1139                               Int_t sizeqmcmd = qtime.size();
1140                               inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);
1141                               qtime[sizeqmcmd]=recqtime[mu];
1142                               q0[sizeqmcmd]=recq0[mu];
1143                               q1[sizeqmcmd]=recq1[mu];
1144                               q2[sizeqmcmd]=recq2[mu];
1145                               q3[sizeqmcmd]=recq3[mu];
1146                               qmode[sizeqmcmd]=-10;
1147                               orbits.getPosition((double) (qtime[sizeqmcmd] - gltle->GetFromTime())/60., &eCi);
1148                               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,recq0[mu],recq1[mu],recq2[mu],recq3[mu]);
1149                               qRoll[sizeqmcmd]=RYPang_upper->Kren;
1150                               qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;
1151                               qPitch[sizeqmcmd]=RYPang_upper->Tangazh;
1152                             }
1153                          }
1154                          if(recqtime[mu]>=u_time){
1155                             if(sqrt(pow(L_QQ_Q_l_upper->quat[0][0],2)+pow(L_QQ_Q_l_upper->quat[0][1],2)+pow(L_QQ_Q_l_upper->quat[0][2],2)+pow(L_QQ_Q_l_upper->quat[0][3],2))>0.99999){
1156                               Int_t sizeqmcmd = qtime.size();
1157                               inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);
1158                               qtime[sizeqmcmd]=u_time;
1159                               q0[sizeqmcmd]=L_QQ_Q_l_upper->quat[0][0];
1160                               q1[sizeqmcmd]=L_QQ_Q_l_upper->quat[0][1];
1161                               q2[sizeqmcmd]=L_QQ_Q_l_upper->quat[0][2];
1162                               q3[sizeqmcmd]=L_QQ_Q_l_upper->quat[0][3];
1163                               qmode[sizeqmcmd]=holeq(lowerqtime,qtime[sizeqmcmd],L_QQ_Q_l_lower,L_QQ_Q_l_upper,ui);
1164                               lowerqtime = u_time;
1165                               orbits.getPosition((double) (u_time - gltle->GetFromTime())/60., &eCi);
1166                               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]);
1167                               qRoll[sizeqmcmd]=RYPang_upper->Kren;
1168                               qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;
1169                               qPitch[sizeqmcmd]=RYPang_upper->Tangazh;
1170                               CopyQ(L_QQ_Q_l_lower,L_QQ_Q_l_upper);
1171                               break;
1172                             }
1173                        }                        }
1174                      }                      }
1175                    }                    }
1176                  }                  }
1177                }                }
1178              }              }
1179              if ( debug ) printf(" ciccio \n");              //if ( debug ) cout << "subpacket " << j3 << "\t qtime = " << qtime[qtime.size()-1] << endl;
1180            }            }
1181          }          }
1182                    
1183          if(qtime.size()==0){          if(qtime.size()==0){                            // in case if no orientation information in data
1184              if ( debug ) cout << "qtime.size() = 0" << endl;
1185            for(UInt_t my=0;my<recqtime.size();my++){            for(UInt_t my=0;my<recqtime.size();my++){
1186              Int_t sizeqmcmd = qtime.size();              if(sqrt(pow(recq0[my],2)+pow(recq1[my],2)+pow(recq2[my],2)+pow(recq3[my],2))>0.99999){
1187              inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);                Int_t sizeqmcmd = qtime.size();
1188              qtime[sizeqmcmd]=recqtime[my];                inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);
1189              q0[sizeqmcmd]=recq0[my];                qtime[sizeqmcmd]=recqtime[my];
1190              q1[sizeqmcmd]=recq1[my];                q0[sizeqmcmd]=recq0[my];
1191              q2[sizeqmcmd]=recq2[my];                q1[sizeqmcmd]=recq1[my];
1192              q3[sizeqmcmd]=recq3[my];                q2[sizeqmcmd]=recq2[my];
1193              qmode[sizeqmcmd]=-10;                q3[sizeqmcmd]=recq3[my];
1194              orbits.getPosition((double) (qtime[sizeqmcmd] - gltle->GetFromTime())/60., &eCi);                qmode[sizeqmcmd]=-10;
1195              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,recq0[my],recq1[my],recq2[my],recq3[my]);                orbits.getPosition((double) (qtime[sizeqmcmd] - gltle->GetFromTime())/60., &eCi);
1196              qRoll[sizeqmcmd]=RYPang_upper->Kren;                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,recq0[my],recq1[my],recq2[my],recq3[my]);
1197              qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;                qRoll[sizeqmcmd]=RYPang_upper->Kren;
1198              qPitch[sizeqmcmd]=RYPang_upper->Tangazh;                qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;
1199                  qPitch[sizeqmcmd]=RYPang_upper->Tangazh;
1200                }
1201            }            }
1202          }          }
1203                    
         if ( debug ) printf(" fuffi \n");  
1204    
1205          if ( debug ) printf(" puffi \n");          if ( debug ) printf(" puffi \n");
1206          Double_t tmin = 9999999999.;          Double_t tmin = 9999999999.;
# Line 1075  int OrbitalInfoCore(UInt_t run, TFile *f Line 1209  int OrbitalInfoCore(UInt_t run, TFile *f
1209            if(qtime[tre]>tmax)tmax = qtime[tre];            if(qtime[tre]>tmax)tmax = qtime[tre];
1210            if(qtime[tre]<tmin)tmin = qtime[tre];            if(qtime[tre]<tmin)tmin = qtime[tre];
1211          }          }
1212          if ( debug ) printf(" gnfuffi \n");          // sorting quaternions by time
1213           Bool_t t = true;
1214            while(t){
1215             t=false;
1216              for(UInt_t i=0;i<qtime.size()-1;i++){
1217                if(qtime[i]>qtime[i+1]){
1218                  Double_t tmpr = qtime[i];
1219                  qtime[i]=qtime[i+1];
1220                  qtime[i+1] = tmpr;
1221                  tmpr = q0[i];
1222                  q0[i]=q0[i+1];
1223                  q0[i+1] = tmpr;
1224                  tmpr = q1[i];
1225                  q1[i]=q1[i+1];
1226                  q1[i+1] = tmpr;
1227                  tmpr = q2[i];
1228                  q2[i]=q2[i+1];
1229                  q2[i+1] = tmpr;
1230                  tmpr = q3[i];
1231                  q3[i]=q3[i+1];
1232                  q3[i+1] = tmpr;
1233                  tmpr = qRoll[i];
1234                  qRoll[i]=qRoll[i+1];
1235                  qRoll[i+1] = tmpr;
1236                  tmpr = qYaw[i];
1237                  qYaw[i]=qYaw[i+1];
1238                  qYaw[i+1] = tmpr;
1239                  tmpr = qPitch[i];
1240                  qPitch[i]=qPitch[i+1];
1241                  qPitch[i+1] = tmpr;
1242                    t=true;
1243                }
1244              }
1245            }
1246    
1247            if ( debug ){
1248              cout << "we have loaded quaternions: size of quaternions set is "<< qtime.size() << endl;
1249             for(UInt_t i=0;i<qtime.size();i++) cout << qtime[i] << "\t";
1250              cout << endl << endl;
1251              Int_t lopu;
1252              cin >> lopu;
1253           }
1254    
1255        } // if we processed first event        } // if we processed first event
1256    
# Line 1083  int OrbitalInfoCore(UInt_t run, TFile *f Line 1258  int OrbitalInfoCore(UInt_t run, TFile *f
1258        //Filling Inclination information        //Filling Inclination information
1259        Double_t incli = 0;        Double_t incli = 0;
1260        if ( qtime.size() > 1 ){        if ( qtime.size() > 1 ){
1261            if ( debug ) cout << "ok quaternions is exist and mu = " << must << endl;
1262            if ( debug ) cout << "qtimes[ " << qtime[0] << " , " << qtime[qtime.size()-1] << " ]\tatime = "<<atime<<endl;
1263          for(UInt_t mu = must;mu<qtime.size()-1;mu++){          for(UInt_t mu = must;mu<qtime.size()-1;mu++){
1264            if ( debug ) printf(" ??grfuffi %i sixe %i must %i \n",mu,qtime.size()-1,must);            if ( debug ) printf(" ??grfuffi %i sixe %i must %i \n",mu,qtime.size()-1,must);
1265            if(qtime[mu+1]>qtime[mu]){            if(qtime[mu+1]>qtime[mu]){
1266              if ( debug ) printf(" grfuffi2 %i \n",mu);              if ( debug ) cout << "qtime[" << mu << "] = " << qtime[mu] << "\tqtime[" << mu+1 << "] = " << qtime[mu+1] << "\tatime = " << atime << endl;
1267              if(atime<=qtime[mu+1] && atime>=qtime[mu]){              if(atime<=qtime[mu+1] && atime>=qtime[mu]){
1268                  if ( debug ) cout << "here we have found proper quaternions for interpolation: mu = "<<mu<<endl;
1269                must = mu;                must = mu;
               if ( debug ) printf(" grfuffi3 %i \n",mu);  
1270                incli = (qPitch[mu+1]-qPitch[mu])/(qtime[mu+1]-qtime[mu]);                incli = (qPitch[mu+1]-qPitch[mu])/(qtime[mu+1]-qtime[mu]);
1271                orbitalinfo->theta =  incli*atime+qPitch[mu+1]-incli*qtime[mu+1];                orbitalinfo->theta =  incli*atime+qPitch[mu+1]-incli*qtime[mu+1];
1272                incli = (qRoll[mu+1]-qRoll[mu])/(qtime[mu+1]-qtime[mu]);                incli = (qRoll[mu+1]-qRoll[mu])/(qtime[mu+1]-qtime[mu]);
# Line 1105  int OrbitalInfoCore(UInt_t run, TFile *f Line 1282  int OrbitalInfoCore(UInt_t run, TFile *f
1282                orbitalinfo->q2 =  incli*atime+q2[mu+1]-incli*qtime[mu+1];                orbitalinfo->q2 =  incli*atime+q2[mu+1]-incli*qtime[mu+1];
1283                incli = (q3[mu+1]-q3[mu])/(qtime[mu+1]-qtime[mu]);                incli = (q3[mu+1]-q3[mu])/(qtime[mu+1]-qtime[mu]);
1284                orbitalinfo->q3 =  incli*atime+q3[mu+1]-incli*qtime[mu+1];                orbitalinfo->q3 =  incli*atime+q3[mu+1]-incli*qtime[mu+1];
1285                                Float_t tg = (qtime[mu+1]-qtime[mu])/1000.;
1286                orbitalinfo->TimeGap = qtime[mu+1]-qtime[mu];                if(tg>=1) tg=0.00;
1287                  orbitalinfo->TimeGap = TMath::Min(TMath::Abs(qtime[mu+1])-atime,TMath::Abs(atime-qtime[mu]))+tg;//qtime[mu+1]-qtime[mu];
1288                orbitalinfo->mode = qmode[mu+1];                orbitalinfo->mode = qmode[mu+1];
1289                  
1290                  //if(atime==qtime[mu] || atime==qtime[mu+1]) orbitalinfo->qkind = 0; else orbitalinfo->qkind=1;
1291                //if(qmode[mu+1]==-10) orbitalinfo->R10r = true;else orbitalinfo->R10r = false;                //if(qmode[mu+1]==-10) orbitalinfo->R10r = true;else orbitalinfo->R10r = false;
               //reserved for next versions Vitaly.  
               /*if(qmode[mu+1]==-10 || qmode[mu+1]==0 || qmode[mu+1]==1 || qmode[mu+1]==3 || qmode[mu+1]==4 || qmode[mu+1]==6){  
               //linear interpolation  
               incli = (q0[mu+1]-q0[mu])/(qtime[mu+1]-qtime[mu]);  
               orbitalinfo->q0 =  incli*atime+q0[mu+1]-incli*qtime[mu+1];  
               incli = (q1[mu+1]-q1[mu])/(qtime[mu+1]-qtime[mu]);  
               orbitalinfo->q1 =  incli*atime+q1[mu+1]-incli*qtime[mu+1];  
               incli = (q2[mu+1]-q2[mu])/(qtime[mu+1]-qtime[mu]);  
               orbitalinfo->q2 =  incli*atime+q2[mu+1]-incli*qtime[mu+1];  
               incli = (q3[mu+1]-q3[mu])/(qtime[mu+1]-qtime[mu]);  
               orbitalinfo->q3 =  incli*atime+q3[mu+1]-incli*qtime[mu+1];  
               }else{  
               //sine interpolation  
               for(UInt_t mt=0;mt<q0sine.size();mt++){  
               if(atime<=q0sine[mt].finishPoint && atime>=q0sine[mt].startPoint){  
               if(!q0sine[mt].NeedFit)orbitalinfo->q0=q0sine[mt].A*sin(q0sine[mt].b*atime+q0sine[mt].c);else{  
               incli = (q0[mu+1]-q0[mu])/(qtime[mu+1]-qtime[mu]);  
               orbitalinfo->q0 =  incli*atime+q0[mu+1]-incli*qtime[mu+1];  
               }  
               }  
               if(atime<=q1sine[mt].finishPoint && atime>=q1sine[mt].startPoint){  
               if(!q1sine[mt].NeedFit)orbitalinfo->q1=q1sine[mt].A*sin(q1sine[mt].b*atime+q1sine[mt].c);else{  
               incli = (q1[mu+1]-q1[mu])/(qtime[mu+1]-qtime[mu]);  
               orbitalinfo->q1 =  incli*atime+q1[mu+1]-incli*qtime[mu+1];  
               }  
               }  
               if(atime<=q2sine[mt].finishPoint && atime>=q2sine[mt].startPoint){  
               if(!q2sine[mt].NeedFit)orbitalinfo->q2=q0sine[mt].A*sin(q2sine[mt].b*atime+q2sine[mt].c);else{  
               incli = (q2[mu+1]-q2[mu])/(qtime[mu+1]-qtime[mu]);  
               orbitalinfo->q2 =  incli*atime+q2[mu+1]-incli*qtime[mu+1];  
               }  
               }  
               if(atime<=q3sine[mt].finishPoint && atime>=q3sine[mt].startPoint){  
               if(!q3sine[mt].NeedFit)orbitalinfo->q3=q0sine[mt].A*sin(q3sine[mt].b*atime+q3sine[mt].c);else{  
               incli = (q3[mu+1]-q3[mu])/(qtime[mu+1]-qtime[mu]);  
               orbitalinfo->q3 =  incli*atime+q3[mu+1]-incli*qtime[mu+1];  
               }  
               }  
               if(atime<=Yawsine[mt].finishPoint && atime>=Yawsine[mt].startPoint){  
               if(!Yawsine[mt].NeedFit)orbitalinfo->phi=Yawsine[mt].A*sin(Yawsine[mt].b*atime+Yawsine[mt].c);else{  
               incli = (qYaw[mu+1]-qYaw[mu])/(qtime[mu+1]-qtime[mu]);  
               orbitalinfo->phi =  incli*atime+qYaw[mu+1]-incli*qtime[mu+1];  
               }  
               }  
               }  
               }*/  
               //q0testing->Fill(atime,orbitalinfo->q0,100);  
               //q1testing->Fill(atime,orbitalinfo->q1,100);  
               //Pitchtesting->Fill(atime,orbitalinfo->etha);  
               //q2testing->Fill(atime,orbitalinfo->q2);  
               //q3testing->Fill(atime,orbitalinfo->q3);  
1292                if ( debug ) printf(" grfuffi4 %i \n",mu);                if ( debug ) printf(" grfuffi4 %i \n",mu);
1293                  
1294                break;                break;
1295              }              }
1296            }            }
# Line 1172  int OrbitalInfoCore(UInt_t run, TFile *f Line 1302  int OrbitalInfoCore(UInt_t run, TFile *f
1302        //        //
1303                
1304        if ( orbitalinfo->q0< -999 || orbitalinfo->q1 < -999 || orbitalinfo->q2 < -999 || orbitalinfo->q3 < -999 || orbitalinfo->q0 != orbitalinfo->q0 || orbitalinfo->q1 != orbitalinfo->q1 || orbitalinfo->q2 != orbitalinfo->q2 || orbitalinfo->q3 != orbitalinfo->q3 ){        if ( orbitalinfo->q0< -999 || orbitalinfo->q1 < -999 || orbitalinfo->q2 < -999 || orbitalinfo->q3 < -999 || orbitalinfo->q0 != orbitalinfo->q0 || orbitalinfo->q1 != orbitalinfo->q1 || orbitalinfo->q2 != orbitalinfo->q2 || orbitalinfo->q3 != orbitalinfo->q3 ){
1305            if ( debug ) cout << "ops no iclination information" << endl;
1306          orbitalinfo->mode = 10;          orbitalinfo->mode = 10;
1307          orbitalinfo->q0 = -1000.;          orbitalinfo->q0 = -1000.;
1308          orbitalinfo->q1 = -1000.;          orbitalinfo->q1 = -1000.;
# Line 1180  int OrbitalInfoCore(UInt_t run, TFile *f Line 1311  int OrbitalInfoCore(UInt_t run, TFile *f
1311          orbitalinfo->etha = -1000.;          orbitalinfo->etha = -1000.;
1312          orbitalinfo->phi = -1000.;          orbitalinfo->phi = -1000.;
1313          orbitalinfo->theta = -1000.;          orbitalinfo->theta = -1000.;
1314            orbitalinfo->TimeGap = -1000.;
1315            //orbitalinfo->qkind = -1000;
1316            
1317            //      if ( debug ){
1318            //        Int_t lopu;
1319            //         cin >> lopu;
1320            //      }
1321          if ( debug ) printf(" grfuffi6 \n");          if ( debug ) printf(" grfuffi6 \n");
1322        }        }
1323        //        //
# Line 1194  int OrbitalInfoCore(UInt_t run, TFile *f Line 1332  int OrbitalInfoCore(UInt_t run, TFile *f
1332        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);
1333        lat = rad2deg(coo.m_Lat);        lat = rad2deg(coo.m_Lat);
1334        alt = coo.m_Alt;        alt = coo.m_Alt;
1335    
1336          cOrbit orbits2(*gltle->GetTle());
1337          orbits2.getPosition((double) (atime - gltle->GetFromTime())/60., &eCi);
1338          //      Float_t x=eCi.getPos().m_x;
1339          //      Float_t y=eCi.getPos().m_y;
1340          //      Float_t z=eCi.getPos().m_z;
1341          
1342          TVector3 V(eCi.getVel().m_x,eCi.getVel().m_y,eCi.getVel().m_z);
1343          TVector3 Pos(eCi.getPos().m_x,eCi.getPos().m_y,eCi.getPos().m_z);
1344          
1345          Float_t dlon=Pos.Phi()*TMath::RadToDeg()-lon;
1346          
1347          Pos.RotateZ(-dlon*TMath::DegToRad());
1348          V.RotateZ(-dlon*TMath::DegToRad());
1349          Float_t diro;
1350          if(V.Z()>0) diro=1; else diro=-1;
1351          
1352          // 10REDNEW
1353          Int_t errq=0;
1354          Int_t azim=0;;
1355          for(UInt_t mu = must;mu<RTtime2.size()-1;mu++){
1356            if(atime<=RTtime2[mu] && atime>=RTtime1[mu]){
1357              errq=RTerrq[mu];
1358              azim=RTazim[mu];
1359            }
1360          }
1361          orbitalinfo->errq = errq;
1362          orbitalinfo->azim = azim;
1363          orbitalinfo->qkind = 0;
1364          
1365        if ( debug ) printf(" coord done \n");        if ( debug ) printf(" coord done \n");
       //  
1366        if( lon<180 && lon>-180 && lat<90 && lat>-90 && alt>0 ){          if( lon<180 && lon>-180 && lat<90 && lat>-90 && alt>0 ){  
1367          //                //      
1368          orbitalinfo->lon = lon;          orbitalinfo->lon = lon;
1369          orbitalinfo->lat = lat;          orbitalinfo->lat = lat;
1370          orbitalinfo->alt = alt ;          orbitalinfo->alt = alt;
1371            orbitalinfo->V = V;
1372    
1373            //      GMtype_CoordGeodetic  location;
1374            location.lambda = lon;
1375            location.phi = lat;
1376            location.HeightAboveEllipsoid = alt;
1377    
1378            GM_GeodeticToSpherical(Ellip, location, &CoordSpherical);
1379            GM_SphericalToCartesian(CoordSpherical,  &CoordCartesian);
1380            GM_EarthCartToDipoleCartCD(Pole, CoordCartesian, &DipoleCartesian);
1381            GM_CartesianToSpherical(DipoleCartesian, &DipoleSpherical);
1382            orbitalinfo->londip = DipoleSpherical.lambda;
1383            orbitalinfo->latdip = DipoleSpherical.phig;
1384    
1385            if(debug)cout<<"geodetic:\t"<<lon<<"\t"<<lat<<"\tgeomagnetic:\t"<<orbitalinfo->londip<<"\t"<<orbitalinfo->latdip<<endl;
1386    
1387          //          //
1388          // compute mag field components and L shell.          // compute mag field components and L shell.
1389          //          //
# Line 1221  int OrbitalInfoCore(UInt_t run, TFile *f Line 1404  int OrbitalInfoCore(UInt_t run, TFile *f
1404          orbitalinfo->L = xl;                orbitalinfo->L = xl;      
1405          // Set Stormer vertical cutoff using L shell.          // Set Stormer vertical cutoff using L shell.
1406          orbitalinfo->cutoffsvl = 14.295 / (xl*xl); //          orbitalinfo->cutoffsvl = 14.295 / (xl*xl); //
1407            if(debug)cout << "L = " << xl << "\tM = " << dimo << "\tvertical cutoff:  "<< orbitalinfo->cutoffsvl << endl;
1408    
1409          /*          /*
1410            ---------- Forwarded message ----------            ---------- Forwarded message ----------
1411            Date: Wed, 09 May 2012 12:16:47 +0200            Date: Wed, 09 May 2012 12:16:47 +0200
# Line 1245  int OrbitalInfoCore(UInt_t run, TFile *f Line 1430  int OrbitalInfoCore(UInt_t run, TFile *f
1430        //        //
1431        // pitch angles        // pitch angles
1432        //        //
1433        //if ( orbitalinfo->mode != 10 && orbitalinfo->mode != 5 && orbitalinfo->mode !=7 && orbitalinfo->mode != 9 ){        if( orbitalinfo->TimeGap>0){
       if( orbitalinfo->TimeGap>0 && orbitalinfo->TimeGap<2000000){  
1434          //          //
1435          if ( debug ) printf(" timegap %f \n",orbitalinfo->TimeGap);          if ( debug ) printf(" timegap %f \n",orbitalinfo->TimeGap);
1436          Float_t Bx = -orbitalinfo->Bdown;          Float_t Bx = -orbitalinfo->Bdown;
1437          Float_t By = orbitalinfo->Beast;          Float_t By = orbitalinfo->Beast;
1438          Float_t Bz = orbitalinfo->Bnorth;          Float_t Bz = orbitalinfo->Bnorth;
1439          //  
1440          TMatrixD Fij = PO->ECItoGreenwich(PO->QuatoECI(orbitalinfo->q0,orbitalinfo->q1,orbitalinfo->q2,orbitalinfo->q3),orbitalinfo->absTime);          TMatrixD Qiji(3,3);
1441            TMatrixD Qij = PO->QuatoECI(orbitalinfo->q0,orbitalinfo->q1,orbitalinfo->q2,orbitalinfo->q3);
1442            TMatrixD Dij = PO->ECItoGEO(Qij,orbitalinfo->absTime,orbitalinfo->lat,orbitalinfo->lon);
1443    
1444    //10REDNEW
1445            /* If initial orientation data have reason to be inaccurate */
1446            Double_t tg = 0;
1447            if ( debug ) cout<<modf(orbitalinfo->TimeGap,&tg)<<endl;
1448    //      if(orbitalinfo->TimeGap>0 && errq==0 && azim==0){               // 10RED CHECK  (comparison between three metod of recovering orientation)
1449            if(((orbitalinfo->TimeGap>60.0 && TMath::Abs(orbitalinfo->etha)>0.5) || errq!=0 || modf(orbitalinfo->TimeGap,&tg)*1000>700 || modf(orbitalinfo->TimeGap,&tg)*1000==0.0 ) && azim==0){           //Standard condition to use this;               One of these two cases should be commented
1450            /*  found in Rotation Table this data for this time interval*/
1451            if(atime<RTtime1[0])
1452              orbitalinfo->azim = 5;        //means that RotationTable no started yet
1453           else{
1454              for(UInt_t mu = must;mu<RTtime2.size()-1;mu++){
1455                if(atime<=RTtime2[mu] && atime>=RTtime1[mu]){
1456                    // search for angle betwean velosity and direction to north in tangential to Earth surfase plane in satellite position
1457                  Double_t tlat=orbitalinfo->lat;
1458    /*            Double_t phint=(163.7-0.0002387*tlat-0.005802*tlat*tlat-0.005802e-7*tlat*tlat*tlat-1.776e-6*tlat*tlat*tlat*tlat+1.395e-10*tlat*tlat*tlat*tlat*tlat);
1459                  Double_t phin=TMath::Abs(90.0*(1+diro)-phint);
1460                  Double_t phi=TMath::Abs(90.0*(1-diro)-TMath::RadToDeg()*atan(TMath::Abs(tan(TMath::DegToRad()*phin))/sqrt(1+pow(tan(TMath::DegToRad()*tlat),2))));
1461                      
1462                    //Get vectors of Satellite reference frame axis in GEO in satndard case (No rotations, all Euler angles equals to 0)
1463                  TVector3 XDij(0,sin(TMath::DegToRad()*phi),cos(TMath::DegToRad()*phi));
1464                  TVector3 YDij(1,0,0);
1465                  TVector3 ZDij(0,sin(TMath::DegToRad()*(phi+90)),cos(TMath::DegToRad()*(phi+90.0)));
1466    
1467                    //Get Vectors to rotate about
1468                  TVector3 B1 = V;
1469                  B1.RotateZ(-lon*TMath::DegToRad());
1470                  B1.RotateY(lat*TMath::DegToRad());
1471                  Float_t elipangle=TMath::ACos((pow(B1.Y(),2)+pow(B1.Z(),2))/B1.Mag()/sqrt(pow(B1.Y(),2)+pow(B1.Z(),2)));
1472                  TVector3 Tre(0,B1.Y(),B1.Z());
1473                  if(B1.X()<0) elipangle=-elipangle;
1474                  TVector3 Vperp=B1;                // axis to rotate around initial Dij on ellip and spitch angles
1475                  Vperp.RotateX(TMath::Pi()/2.);
1476                  Vperp.SetX(0);
1477    */            Double_t kar=(RTbank2[mu]-RTbank1[mu])/(RTtime2[mu]-RTtime1[mu]);
1478                  Double_t bak=RTbank1[mu]-kar*RTtime1[mu];
1479                  Double_t bank=kar*atime+bak;
1480                  Float_t spitch = 0.00001;  // temprary not zero to avoid problem with tranzition from Euler angles to orientation matrix
1481    
1482                    //Estimations of pitch angle of satellite
1483                  if(TMath::Abs(bank)>0.7){
1484                     Float_t spitch1=TMath::DegToRad()*0.7*RTdir1[mu];
1485                     Float_t spitch2=TMath::DegToRad()*0.7*RTdir2[mu];
1486                     Float_t kva=(spitch2-spitch1)/(RTtime2[mu]-RTtime1[mu]);
1487                     Float_t bva=spitch1-kva*RTtime1[mu];
1488                     spitch=kva*atime+bva;
1489                  }
1490    /*            //spitch=0.0;
1491                  //Rotations future Dij matrix on ellip and spitch angles
1492                  XDij.Rotate(-elipangle-spitch,Vperp);
1493                  YDij.Rotate(-elipangle-spitch,Vperp);
1494                  ZDij.Rotate(-elipangle-spitch,Vperp);
1495                    
1496                    //Rotation on bank angle;
1497                  if(TMath::Abs(bank)>0.5){
1498                     XDij.Rotate(TMath::DegToRad()*bank,B1);
1499                     YDij.Rotate(TMath::DegToRad()*bank,B1);
1500                     ZDij.Rotate(TMath::DegToRad()*bank,B1);
1501                  }
1502                  Dij(0,0)=XDij.X(); Dij(1,0)=XDij.Y(); Dij(2,0)=XDij.Z();
1503                  Dij(0,1)=YDij.X(); Dij(1,1)=YDij.Y(); Dij(2,1)=YDij.Z();
1504                  Dij(0,2)=ZDij.X(); Dij(1,2)=ZDij.Y(); Dij(2,2)=ZDij.Z();
1505    */
1506                  //Calculate Yaw angle accordingly with fit, see picture FitYaw.jpg
1507                  Double_t yaw=0.00001;  // temprary not zero to avoid problem with tranzition from Euler angles to orientation matrix
1508                  if(TMath::Abs(tlat)<70)
1509                        yaw = -3.7e-8*tlat*tlat*tlat*tlat + 1.4e-7*tlat*tlat*tlat - 0.0005*tlat*tlat - 0.00025*tlat + 3.6;
1510                  yaw = diro*yaw;   //because should be different sign for ascending and descending orbits!
1511    
1512                  if(TMath::Abs(bank)>0.5 && TMath::Abs(yaw-orbitalinfo->phi)<3.0) yaw=orbitalinfo->phi;
1513    
1514    //            Qiji = PO->EulertoEci(eCi.getPos().m_x,eCi.getPos().m_y,eCi.getPos().m_z,eCi.getVel().m_x,eCi.getVel().m_y,eCi.getVel().m_z,bank,yaw,spitch);             // 10RED CHECK
1515                  Qij = PO->EulertoEci(eCi.getPos().m_x,eCi.getPos().m_y,eCi.getPos().m_z,eCi.getVel().m_x,eCi.getVel().m_y,eCi.getVel().m_z,bank,yaw,spitch);              // STANDARD
1516                  orbitalinfo->qkind = 1;
1517    
1518                  break;
1519                }
1520              }     // enf of  loop for(UInt_t mu = must;mu<RTtime2.size()-1;mu++){
1521    
1522              //Qij = PO->GEOtoECI(Dij,orbitalinfo->absTime,orbitalinfo->lat,orbitalinfo->lon); // to convert from Dij to Qij
1523    
1524            } // end of if(atime<RTtime1[0]
1525            } // end of f(((orbitalinfo->TimeGap>60.0 && TMath...
1526    
1527            TMatrixD qij = PO->ColPermutation(Qij);
1528            TMatrixD Fij = PO->ECItoGreenwich(Qij,orbitalinfo->absTime);
1529          TMatrixD Gij = PO->ColPermutation(Fij);          TMatrixD Gij = PO->ColPermutation(Fij);
1530          TMatrixD Dij = PO->GreenwichtoGEO(orbitalinfo->lat,orbitalinfo->lon,Fij);          Dij = PO->ECItoGEO(Qij,orbitalinfo->absTime,orbitalinfo->lat,orbitalinfo->lon);
1531          TMatrixD Iij = PO->ColPermutation(Dij);          TMatrixD Iij = PO->ColPermutation(Dij);
1532          //          TVector3 SP = PO->GetSunPosition(orbitalinfo->absTime);
1533            // go to Pamela reference frame from Resurs reference frame
1534            Float_t tmpy = SP.Y();
1535            SP.SetY(SP.Z());
1536            SP.SetZ(-tmpy);
1537            TVector3 SunZenith;
1538            SunZenith.SetMagThetaPhi(1,23.439281*TMath::DegToRad(),TMath::Pi()/2.);
1539            TVector3 SunMag = -SP;
1540            SunMag.Rotate(-45*TMath::DegToRad(),SunZenith);
1541            tmpy=SunMag.Y();
1542            SunMag.SetY(SunMag.Z());
1543            SunMag.SetZ(-tmpy);
1544    
1545          orbitalinfo->Iij.ResizeTo(Iij);          orbitalinfo->Iij.ResizeTo(Iij);
1546          orbitalinfo->Iij = Iij;          orbitalinfo->Iij = Iij;
1547          //          //
# Line 1276  int OrbitalInfoCore(UInt_t run, TFile *f Line 1560  int OrbitalInfoCore(UInt_t run, TFile *f
1560            for(Int_t nt=0; nt < tof->ntrk(); nt++){              for(Int_t nt=0; nt < tof->ntrk(); nt++){  
1561              //              //
1562              ToFTrkVar *ptt = tof->GetToFTrkVar(nt);              ToFTrkVar *ptt = tof->GetToFTrkVar(nt);
1563                if (debug) cout<<"tof->ntrk() = "<<tof->ntrk()<<"\tptt->trkseqno = "<<ptt->trkseqno<<"\ttrke->ntrk() = "<<trke->ntrk()<<endl;
1564              Double_t E11x = ptt->xtr_tof[0]; // tr->x[0];              Double_t E11x = ptt->xtr_tof[0]; // tr->x[0];
1565              Double_t E11y = ptt->ytr_tof[0]; //tr->y[0];              Double_t E11y = ptt->ytr_tof[0]; //tr->y[0];
1566              Double_t E11z = zin[0];              Double_t E11z = zin[0];
# Line 1283  int OrbitalInfoCore(UInt_t run, TFile *f Line 1568  int OrbitalInfoCore(UInt_t run, TFile *f
1568              Double_t E22y = ptt->ytr_tof[3];//tr->y[3];              Double_t E22y = ptt->ytr_tof[3];//tr->y[3];
1569              Double_t E22z = zin[3];              Double_t E22z = zin[3];
1570              if ( (E11x < 100. && E11y < 100. && E22x < 100. && E22y < 100.) || ptt->trkseqno != -1  ){              if ( (E11x < 100. && E11y < 100. && E22x < 100. && E22y < 100.) || ptt->trkseqno != -1  ){
1571                  TrkTrack *mytrack = trke->GetStoredTrack(ptt->trkseqno);
1572                  Float_t rig=1/mytrack->GetDeflection();
1573                Double_t norm = sqrt(pow(E22x-E11x,2)+pow(E22y-E11y,2)+pow(E22z-E11z,2));                Double_t norm = sqrt(pow(E22x-E11x,2)+pow(E22y-E11y,2)+pow(E22z-E11z,2));
1574                //              Double_t MyAzim = TMath::RadToDeg()*atan(TMath::Abs(E22y-E11y)/TMath::Abs(E22x-E11x));                //              Double_t MyAzim = TMath::RadToDeg()*atan(TMath::Abs(E22y-E11y)/TMath::Abs(E22x-E11x));
1575                //              if(E22x-E11x>=0 && E22y-E11y <0) MyAzim =  360. - MyAzim;                //              if(E22x-E11x>=0 && E22y-E11y <0) MyAzim =  360. - MyAzim;
# Line 1305  int OrbitalInfoCore(UInt_t run, TFile *f Line 1592  int OrbitalInfoCore(UInt_t run, TFile *f
1592                //                            //            
1593                t_orb->pitch = (Float_t)PO->GetPitchAngle(Eij(0,0),Eij(1,0),Eij(2,0),Bx,By,Bz);                t_orb->pitch = (Float_t)PO->GetPitchAngle(Eij(0,0),Eij(1,0),Eij(2,0),Bx,By,Bz);
1594                //                //
1595                    // SunPosition in instrumental reference frame
1596                  TMatrixD Kij = PO->PamelatoGEO(qij,Px,Py,Pz);
1597                  TMatrixD Lij = PO->PamelatoGEO(qij,0,0,1);
1598                  t_orb->sunangle=(Float_t)PO->GetPitchAngle(Kij(0,0),Kij(1,0),Kij(2,0),-SP.X(),-SP.Y(),-SP.Z());
1599                  t_orb->sunmagangle=(Float_t)PO->GetPitchAngle(Kij(0,0),Kij(1,0),Kij(2,0),SunMag.X(),SunMag.Y(),SunMag.Z());
1600    
1601                //                //
               Double_t omega = PO->GetPitchAngle(Eij(0,0),Eij(1,0),Eij(2,0),cos(orbitalinfo->lon+TMath::Pi()/2)-sin(orbitalinfo->lon+TMath::Pi()/2),cos(orbitalinfo->lon+TMath::Pi()/2)+sin(orbitalinfo->lon+TMath::Pi()/2),1);  
1602                //                //
1603                t_orb->cutoff = 59.3/(pow(orbitalinfo->L,2)*pow((1+sqrt(1-pow(orbitalinfo->L,-3/2)*cos(omega))),2));                //Double_t omega = PO->GetPitchAngle(Eij(0,0),Eij(1,0),Eij(2,0),cos(orbitalinfo->lon+TMath::Pi()/2)-sin(orbitalinfo->lon+TMath::Pi()/2),cos(orbitalinfo->lon+TMath::Pi()/2)+sin(orbitalinfo->lon+TMath::Pi()/2),1);
1604                  Double_t omega = PO->GetPitchAngle(-Eij(0,0),-Eij(1,0),-Eij(2,0),1,0,0) * TMath::DegToRad();
1605                  TVector3 Bxy(0,By,Bz);
1606                  TVector3 Exy(0,-Eij(1,0),-Eij(2,0));
1607                  Double_t dzeta=Bxy.Angle(Exy);
1608                  if (-Eij(1,0) < 0) dzeta=2.0*TMath::Pi() - dzeta;
1609    
1610                    if(debug) cout << "omega = "<<omega*TMath::RadToDeg()<<"\tdzeta = "<<dzeta*TMath::RadToDeg()<<endl;
1611    
1612                    // Formula from D.F. Smart *, M.A. Shea [2005]; A review of geomagnetic cutoff rigidities for earth-orbiting spacecraft
1613                  if(rig>=0) t_orb->cutoff = 59.3/(pow(orbitalinfo->L,2)*pow(1+sqrt(1-sin(omega)*sin(dzeta)*pow(cos(orbitalinfo->lat*TMath::DegToRad()),3)),2));
1614                  else t_orb->cutoff = 59.3/(pow(orbitalinfo->L,2)*pow(1+sqrt(1-sin(omega)*sin(TMath::Pi()+dzeta)*pow(cos(orbitalinfo->lat*TMath::DegToRad()),3)),2));
1615                  if (debug) cout << "R = " << rig << "\tcutoff = " << t_orb->cutoff << endl;
1616    
1617                  //t_orb->cutoff = 59.3/(pow(orbitalinfo->L,2)*pow((1+sqrt(1-pow(orbitalinfo->L,-3/2)*cos(omega))),2));
1618    
1619                //                //
1620                if ( t_orb->pitch != t_orb->pitch ) t_orb->pitch = -1000.;                if ( t_orb->pitch != t_orb->pitch ) t_orb->pitch = -1000.;
1621                if ( t_orb->cutoff != t_orb->cutoff ) t_orb->cutoff = -1000.;                if ( t_orb->cutoff != t_orb->cutoff ) t_orb->cutoff = -1000.;
1622                  if ( t_orb->sunangle != t_orb->sunangle ) t_orb->sunangle = -1000.;
1623                  if ( t_orb->sunmagangle != t_orb->sunmagangle ) t_orb->sunmagangle = -1000.;
1624                //                //
1625                if ( debug ) printf(" orbitalinfo->cutoffsvl %f vitaly %f \n",orbitalinfo->cutoffsvl,t_orb->cutoff);                if ( debug ) printf(" orbitalinfo->cutoffsvl %f vitaly %f \n",orbitalinfo->cutoffsvl,t_orb->cutoff);
1626                //                //
# Line 1344  int OrbitalInfoCore(UInt_t run, TFile *f Line 1653  int OrbitalInfoCore(UInt_t run, TFile *f
1653                //                            //            
1654                t_orb->pitch = -1000.;                t_orb->pitch = -1000.;
1655                //                //
1656                  t_orb->sunangle = -1000.;
1657                  //
1658                  t_orb->sunmagangle = -1000;
1659                  //
1660                t_orb->cutoff = -1000.;                t_orb->cutoff = -1000.;
1661                //                //
1662                new(tor[nn]) OrbitalInfoTrkVar(*t_orb);                new(tor[nn]) OrbitalInfoTrkVar(*t_orb);
# Line 1355  int OrbitalInfoCore(UInt_t run, TFile *f Line 1668  int OrbitalInfoCore(UInt_t run, TFile *f
1668              //              //
1669            };                };    
1670          };          };
1671        };        };        // if( orbitalinfo->TimeGap>0){
1672        //        //
1673        // Fill the class        // Fill the class
1674        //        //
# Line 1368  int OrbitalInfoCore(UInt_t run, TFile *f Line 1681  int OrbitalInfoCore(UInt_t run, TFile *f
1681      // Here you may want to clear some variables before processing another run        // Here you may want to clear some variables before processing another run  
1682      //      //
1683    
     //gStyle->SetOptStat(000000);  
     //gStyle->SetPalette(1);  
       
     /*TCanvas* c1 = new TCanvas("c1","",1200,800);  
     //c1->Divide(1,4);  
     c1->cd(1);  
     //q0testing->Draw("colz");  
     //c1->cd(2);  
     //q1testing->Draw("colz");  
     //c1->cd(3);  
     Pitchtesting->Draw("colz");  
     //c1->cd(4);  
     //q3testing->Draw("colz");  
     c1->SaveAs("9.Rollhyst.png");  
     delete c1;*/  
   
1684      if ( verbose ) printf(" Clear before new run \n");      if ( verbose ) printf(" Clear before new run \n");
1685      delete dbtime;      delete dbtime;
1686    
# Line 1442  int OrbitalInfoCore(UInt_t run, TFile *f Line 1739  int OrbitalInfoCore(UInt_t run, TFile *f
1739    //    //
1740    // Close files, delete old tree(s), write and close level2 file    // Close files, delete old tree(s), write and close level2 file
1741    //    //
1742    
1743    if ( l0File ) l0File->Close();    if ( l0File ) l0File->Close();
1744    if ( myfold ) gSystem->Unlink(tempname.str().c_str());    if ( myfold ) gSystem->Unlink(tempname.str().c_str());
1745    //    //
# Line 1477  int OrbitalInfoCore(UInt_t run, TFile *f Line 1775  int OrbitalInfoCore(UInt_t run, TFile *f
1775    if ( runinfo ) runinfo->Close();        if ( runinfo ) runinfo->Close();    
1776    if ( runinfo ) delete runinfo;    if ( runinfo ) delete runinfo;
1777    
1778      if ( tof ) delete tof;
1779      if ( trke ) delete trke;
1780    
1781    if ( debug ){      if ( debug ){  
1782    cout << "1   0x" << OrbitalInfotr << endl;    cout << "1   0x" << OrbitalInfotr << endl;
1783    cout << "2   0x" << OrbitalInfotrclone << endl;    cout << "2   0x" << OrbitalInfotrclone << endl;
# Line 1591  void inclresize(vector<Double_t>& t,vect Line 1892  void inclresize(vector<Double_t>& t,vect
1892    Yaw.resize(sizee);    Yaw.resize(sizee);
1893  }  }
1894    
1895  //Find fitting sine functions for q0,q1,q2,q3 and Yaw-angle;  // geomagnetic calculation staff
1896  void sineparam(vector<Sine>& qsine, vector<Double_t>& qtime, vector<Float_t>& q, vector<Float_t>& Roll, vector<Float_t>& Pitch, Float_t limsin){  
1897    UInt_t mulast = 0;  //void GM_ScanIGRF(TString PATH, GMtype_Data *G0, GMtype_Data *G1, GMtype_Data *H1)
1898    UInt_t munow = 0;  void GM_ScanIGRF(TSQLServer *dbc, GMtype_Data *G0, GMtype_Data *G1, GMtype_Data *H1)
1899    UInt_t munext = 0;  {
1900    Bool_t increase = false;    GL_PARAM *glp = new GL_PARAM();
1901    Bool_t decrease = false;    Int_t parerror=glp->Query_GL_PARAM(1,304,dbc); // parameters stored in DB in GL_PRAM table  
1902    Bool_t Max_is_defined = false;    if ( parerror<0 ) {
1903    Bool_t Start_point_is_defined = false;      throw -902;
1904    Bool_t Period_is_defined = false;    }
1905    Bool_t Large_gap = false;          /*This function scans inputs G0, G1, and H1 of the IGRF table into 3 data arrays*/
1906    Bool_t normal_way = true;    //    TString SATH="/data03/Malakhov/pam9Malakhov/installed10/calib/orb-param/";
1907    Bool_t small_gap_on_ridge = false;          int i;
1908    Double_t t1 = 0;          double temp;
1909    Double_t t1A = 0;          char buffer[200];
1910    Int_t sinesize = 0;          FILE *IGRF;
1911    Int_t nfi = 0;          IGRF = fopen((glp->PATH+glp->NAME).Data(), "r");
1912    for(UInt_t mu = 0;mu<qtime.size();mu++){          //      IGRF = fopen(PATH+"IGRF.tab", "r");
1913      //cout<<"Roll["<<mu<<"] = "<<Roll[mu]<<endl;          G0->size = 25;
1914      if(TMath::Abs(Roll[mu])<1. && TMath::Abs(Pitch[mu])<1. && TMath::Abs(q[mu])<limsin){          G1->size = 25;
1915      //cout<<"q["<<mu<<endl<<"] = "<<q[mu]<<endl;          H1->size = 25;
1916      if(mulast!=0 && munow!=0 && munext!=0){mulast=munow;munow=munext;munext=mu;}          for( i = 0; i < 4; i++)
1917      if(munext==0 && munow!=0)munext=mu;          {
1918      if(munow==0 && mulast!=0)munow=mu;                  fgets(buffer, 200, IGRF);
     if(mulast==0)mulast=mu;  
       
     //cout<<"mulast = "<<mulast<<"\tmunow = "<<munow<<"\tmunext = "<<munext<<endl;  
     //Int_t ref;  
     //cin>>ref;  
     if(TMath::Abs(q[munow])>TMath::Abs(q[mulast]) && TMath::Abs(q[munow])>TMath::Abs(q[munext]) && q[mulast]*q[munext]>0 && qtime[munext]-qtime[mulast]>400)small_gap_on_ridge = true;  
     if(munext>mulast && (qtime[munext]-qtime[mulast]>=2000 || qtime[munext]-qtime[mulast]<0)){if(Large_gap){normal_way = false;Large_gap = false;}else{Large_gap = true;normal_way = false;}}else normal_way = true;  
     //if(normal_way)cout<<"Normal_Way"<<endl;  
     if(Large_gap || small_gap_on_ridge){  
       //cout<<"Large gap..."<<endl;  
       //if(small_gap_on_ridge)cout<<"small gap..."<<endl;  
       //cout<<"q["<<mulast<<"] = "<<q[mulast]<<"\tq["<<munow<<"] = "<<q[munow]<<"\tq["<<munext<<"] = "<<q[munext]<<endl;  
       //cout<<"qtime["<<mulast<<"] = "<<qtime[mulast]<<"\tqtime["<<munow<<"] = "<<qtime[munow]<<"\tqtime["<<munext<<"] = "<<qtime[munext]<<endl;  
       increase = false;  
       decrease = false;  
       if(nfi>0){  
         qsine.resize(qsine.size()-1);  
         sinesize = qsine.size();  
         //cout<<"nfi was larger then zero"<<endl;  
       }else{  
         //cout<<"nfi was not larger then zero :( nfi = "<<nfi<<endl;  
         //cout<<"qsine.size = "<<qsine.size()<<endl;  
         if(!Period_is_defined){  
           //cout<<"Period was defined"<<endl;  
           if(qsine.size()>1){  
             qsine[sinesize-1].b = qsine[sinesize-2].b;  
             qsine[sinesize-1].c = qsine[sinesize-2].c;  
           }else{  
             qsine[sinesize-1].b = TMath::Pi()/1591.54;  
             qsine[sinesize-1].c = qsine[sinesize-1].startPoint;  
           }  
         }  
         if(!Max_is_defined){  
           //cout<<"Max was already defined"<<endl;  
           if(qsine.size()>1)qsine[sinesize-1].A = qsine[sinesize-2].A;else qsine[sinesize-1].A = limsin;  
         }  
         qsine[sinesize-1].NeedFit = true;  
       }  
       qsine[sinesize-1].finishPoint = qtime[munow];  
       //cout<<"finish point before large gap = "<<qtime[munow]<<endl;  
       nfi = 0;  
       Max_is_defined = false;  
       Start_point_is_defined = false;  
       Period_is_defined = false;  
       small_gap_on_ridge = false;  
     }  
     //cout<<"Slope "<<increase<<"\t"<<decrease<<endl;  
     //cout<<"mulast = "<<mulast<<"\tmunow = "<<munow<<"\tmunext = "<<munext<<endl;  
     if((munext>munow) && (munow>mulast) && normal_way){  
       if(!increase && !decrease){  
         //cout<<"Normal way have started"<<endl;  
         qsine.resize(qsine.size()+1);  
         sinesize = qsine.size();  
         qsine[sinesize-1].startPoint=qtime[mulast];  
         if(q[munext]>q[munow] && q[munow]>q[mulast]) increase = true;  
         if(q[munext]<q[munow] && q[munow]<q[mulast]) decrease = true;  
       }  
       //if(TMath::Abs(q[munow])>TMath::Abs(q[mulast]) && TMath::Abs(q[munow])>TMath::Abs(q[munext]) && TMath::Abs(q[munow])>limsin && qtime[munow]-qtime[mulast]>=400 || qtime[munext]-qtime[munow]>=400){small_gap_on_ridge = true;mu--;continue;}  
       if(TMath::Abs(q[munow])>TMath::Abs(q[mulast]) && TMath::Abs(q[munow])>TMath::Abs(q[munext]) && TMath::Abs(q[munow])>0.9*limsin && qtime[munow]-qtime[mulast]<400 && qtime[munext]-qtime[munow]<400){  
         //cout<<"Max point is qtime = "<<qtime[munow]<<"\tq = "<<q[munow]<<endl;  
         if(q[munow]>q[mulast]){  
           increase = false;  
           decrease = true;  
1919          }          }
1920          if(q[munow]<q[mulast]){          fscanf(IGRF, "g 1 0 %lf ", &G0->element[0]);
1921            increase = true;          for(i = 1; i <= 22; i++)
1922            decrease = false;          {
1923                    fscanf(IGRF ,"%lf ", &G0->element[i]);
1924          }          }
1925          if(Max_is_defined && !Start_point_is_defined){          fscanf(IGRF ,"%lf\n", &temp);
1926            Double_t qPer = qtime[munow]-t1A;          G0->element[23] = temp * 5 + G0->element[22];
1927            if(qPer>1000){          G0->element[24] = G0->element[23] + 5 * temp;
1928              //cout<<"qsine["<<sinesize-1<<"] = "<<qPer<<" = "<<qtime[munow]<<" - "<<t1A<<"\tlim = "<<limsin<<endl;          fscanf(IGRF, "g 1 1 %lf ", &G1->element[0]);
1929              qsine[sinesize-1].b=TMath::Pi()/qPer;          for(i = 1; i <= 22; i++)
1930              if(decrease)qsine[sinesize-1].c=-qsine[sinesize-1].b*t1A;          {
1931              if(increase)qsine[sinesize-1].c=-qsine[sinesize-1].b*(t1A-qPer);                  fscanf( IGRF, "%lf ", &G1->element[i]);
             Period_is_defined = true;  
           }  
1932          }          }
1933          Max_is_defined = true;          fscanf(IGRF, "%lf\n", &temp);
1934          qsine[sinesize-1].A = TMath::Abs(q[munow]);          G1->element[23] = temp * 5 + G1->element[22];
1935          if(Start_point_is_defined && Period_is_defined){          G1->element[24] = temp * 5 + G1->element[23];
1936            qsine[sinesize-1].finishPoint = qtime[munow];          fscanf(IGRF, "h 1 1 %lf ", &H1->element[0]);
1937            nfi++;          for(i = 1; i <= 22; i++)
1938            qsine[sinesize-1].NeedFit = false;          {
1939            Max_is_defined = false;                  fscanf( IGRF, "%lf ", &H1->element[i]);
           Start_point_is_defined = false;  
           Period_is_defined = false;  
           qsine.resize(qsine.size()+1);  
           sinesize = qsine.size();  
           qsine[sinesize-1].startPoint = qtime[munow];  
1940          }          }
1941          if(!Start_point_is_defined) t1A=qtime[munow];          fscanf(IGRF, "%lf\n", &temp);
1942        }          H1->element[23] = temp * 5 + H1->element[22];
1943        //if((q[munow]>=0 && q[mulast]<=0) || (q[munow]<=0 && q[mulast]>=0))cout<<"cross zero point diference = "<<qtime[munext] - qtime[mulast]<<"\tqlast = "<<qtime[mulast]<<"\tqnow = "<<qtime[munow]<<"\tqnext = "<<qtime[munext]<<endl;          H1->element[24] = temp * 5 + H1->element[23];
1944        if(((q[munow]>=0 && q[mulast]<=0) || (q[munow]<=0 && q[mulast]>=0)) && qtime[munow]-qtime[mulast]<2000 && qtime[munext]-qtime[munow]<2000){    if ( glp ) delete glp;
1945          Double_t tcrosszero = 0;  } /*GM_ScanIGRF*/
         //cout<<"cross zero point...qtime = "<<qtime[munow]<<endl;  
         if(q[munow]==0.) tcrosszero = qtime[munow];else  
           if(q[mulast]==0.)tcrosszero = qtime[mulast];else{  
             Double_t k_ = (q[munow]-q[mulast])/(qtime[munow]-qtime[mulast]);  
             Double_t b_ = q[munow]-k_*qtime[munow];  
             tcrosszero = -b_/k_;  
           }  
         if(Start_point_is_defined){  
           //cout<<"Start Point allready defined"<<endl;  
           Double_t qPer = tcrosszero - t1;  
           qsine[sinesize-1].b = TMath::Pi()/qPer;  
           //cout<<"qsine["<<sinesize-1<<"].b = "<<TMath::Pi()/qPer<<endl;  
           Period_is_defined = true;  
           Float_t x0 = 0;  
           if(decrease)x0 = t1;  
           if(increase)x0 = tcrosszero;  
           qsine[sinesize-1].c = -qsine[sinesize-1].b*x0;  
           if(Max_is_defined){  
             //cout<<"Max was previous defined"<<endl;  
             qsine[sinesize-1].finishPoint = qtime[munow];  
             nfi++;  
             qsine[sinesize-1].NeedFit = false;  
             Max_is_defined = false;  
             t1 = tcrosszero;  
             Start_point_is_defined = true;  
             Period_is_defined = false;  
             qsine.resize(qsine.size()+1);  
             sinesize = qsine.size();  
             qsine[sinesize-1].startPoint = qtime[munow];  
           }  
         }else{  
           t1 = tcrosszero;  
           Start_point_is_defined = true;  
         }  
       }  
     }  
     }  
   }  
1946    
1947    //cout<<"FINISH SINE INTERPOLATION FUNCTION..."<<endl<<endl;  void GM_SetEllipsoid(GMtype_Ellipsoid *Ellip)
1948  }  {
1949            /*This function sets the WGS84 reference ellipsoid to its default values*/
1950            Ellip->a        =                       6378.137; /*semi-major axis of the ellipsoid in */
1951            Ellip->b        =                       6356.7523142;/*semi-minor axis of the ellipsoid in */
1952            Ellip->fla      =                       1/298.257223563;/* flattening */
1953            Ellip->eps      =                       sqrt(1- ( Ellip->b *    Ellip->b) / (Ellip->a * Ellip->a ));  /*first eccentricity */
1954            Ellip->epssq    =                       (Ellip->eps * Ellip->eps);   /*first eccentricity squared */
1955            Ellip->re       =                       6371.2;/* Earth's radius */
1956    } /*GM_SetEllipsoid*/
1957    
1958    
1959    void GM_EarthCartToDipoleCartCD(GMtype_Pole Pole, GMtype_CoordCartesian EarthCoord, GMtype_CoordCartesian *DipoleCoords)
1960    {
1961            /*This function converts from Earth centered cartesian coordinates to dipole centered cartesian coordinates*/
1962            double X, Y, Z, CosPhi, SinPhi, CosLambda, SinLambda;
1963            CosPhi = cos(TMath::DegToRad()*Pole.phi);
1964            SinPhi = sin(TMath::DegToRad()*Pole.phi);
1965            CosLambda = cos(TMath::DegToRad()*Pole.lambda);
1966            SinLambda = sin(TMath::DegToRad()*Pole.lambda);
1967            X = EarthCoord.x;
1968            Y = EarthCoord.y;
1969            Z = EarthCoord.z;
1970            
1971            /*These equations are taken from a document by Wallace H. Campbell*/
1972            DipoleCoords->x = X * CosPhi * CosLambda + Y * CosPhi * SinLambda - Z * SinPhi;
1973            DipoleCoords->y = -X * SinLambda + Y * CosLambda;
1974            DipoleCoords->z = X * SinPhi * CosLambda + Y * SinPhi * SinLambda + Z * CosPhi;
1975    } /*GM_EarthCartToDipoleCartCD*/
1976    
1977    void GM_GeodeticToSpherical(GMtype_Ellipsoid Ellip, GMtype_CoordGeodetic CoordGeodetic, GMtype_CoordSpherical *CoordSpherical)
1978    {
1979            double CosLat, SinLat, rc, xp, zp; /*all local variables */
1980            /*
1981            ** Convert geodetic coordinates, (defined by the WGS-84
1982            ** reference ellipsoid), to Earth Centered Earth Fixed Cartesian
1983            ** coordinates, and then to spherical coordinates.
1984            */
1985    
1986            CosLat = cos(TMath::DegToRad()*CoordGeodetic.phi);
1987            SinLat = sin(TMath::DegToRad()*CoordGeodetic.phi);
1988    
1989            /* compute the local radius of curvature on the WGS-84 reference ellipsoid */
1990    
1991            rc = Ellip.a / sqrt(1.0 - Ellip.epssq * SinLat * SinLat);
1992    
1993            /* compute ECEF Cartesian coordinates of specified point (for longitude=0) */
1994    
1995            xp = (rc + CoordGeodetic.HeightAboveEllipsoid) * CosLat;
1996            zp = (rc*(1.0 - Ellip.epssq) + CoordGeodetic.HeightAboveEllipsoid) * SinLat;
1997    
1998            /* compute spherical radius and angle lambda and phi of specified point */
1999    
2000            CoordSpherical->r = sqrt(xp * xp + zp * zp);
2001            CoordSpherical->phig = TMath::RadToDeg()*asin(zp / CoordSpherical->r);     /* geocentric latitude */
2002            CoordSpherical->lambda = CoordGeodetic.lambda;                   /* longitude */
2003    } /*GM_GeodeticToSpherical*/
2004    
2005    void GM_PoleLocation(GMtype_Model Model, GMtype_Pole *Pole)
2006    {
2007            /*This function finds the location of the north magnetic pole in spherical coordinates.  The equations are
2008            **from Wallace H. Campbell's Introduction to Geomagnetic Fields*/
2009    
2010            Pole->phi = TMath::RadToDeg()*-atan(sqrt(Model.h1 * Model.h1 + Model.g1 * Model.g1)/Model.g0);
2011            Pole->lambda = TMath::RadToDeg()*atan(Model.h1/Model.g1);
2012    } /*GM_PoleLocation*/
2013    
2014    void GM_SphericalToCartesian(GMtype_CoordSpherical CoordSpherical, GMtype_CoordCartesian *CoordCartesian)
2015    {
2016            /*This function converts spherical coordinates into Cartesian coordinates*/
2017            double CosPhi = cos(TMath::DegToRad()*CoordSpherical.phig);
2018            double SinPhi = sin(TMath::DegToRad()*CoordSpherical.phig);
2019            double CosLambda = cos(TMath::DegToRad()*CoordSpherical.lambda);
2020            double SinLambda = sin(TMath::DegToRad()*CoordSpherical.lambda);
2021            
2022            CoordCartesian->x = CoordSpherical.r * CosPhi * CosLambda;
2023            CoordCartesian->y = CoordSpherical.r * CosPhi * SinLambda;
2024            CoordCartesian->z = CoordSpherical.r * SinPhi;
2025    } /*GM_SphericalToCartesian*/
2026    
2027    void GM_TimeAdjustCoefs(Float_t year, Float_t jyear, GMtype_Data g0d, GMtype_Data g1d, GMtype_Data h1d, GMtype_Model *Model)
2028    {
2029            /*This function calls GM_LinearInterpolation for the coefficients to estimate the value of the
2030            **coefficient for the given date*/
2031            int index;
2032            double x;
2033            index = (year - GM_STARTYEAR) / 5;
2034            x = (jyear - GM_STARTYEAR) / 5;
2035            Model->g0 = GM_LinearInterpolation(index, index+1, g0d.element[index], g0d.element[index+1], x);
2036            Model->g1 = GM_LinearInterpolation(index, index+1, g1d.element[index], g1d.element[index+1], x);
2037            Model->h1 = GM_LinearInterpolation(index, index+1, h1d.element[index], h1d.element[index+1], x);
2038    } /*GM_TimeAdjustCoefs*/
2039    
2040    double GM_LinearInterpolation(double x1, double x2, double y1, double y2, double x)
2041    {
2042            /*This function takes a linear interpolation between two given points for x*/
2043            double weight, y;
2044            weight  = (x - x1) / (x2 - x1);
2045            y = y1 * (1 - weight) + y2 * weight;
2046            return y;
2047    }/*GM_LinearInterpolation*/
2048    
2049    void GM_CartesianToSpherical(GMtype_CoordCartesian CoordCartesian, GMtype_CoordSpherical *CoordSpherical)
2050    {
2051            /*This function converts a point from Cartesian coordinates into spherical coordinates*/
2052            double X, Y, Z;
2053            
2054            X = CoordCartesian.x;
2055            Y = CoordCartesian.y;
2056            Z = CoordCartesian.z;
2057    
2058            CoordSpherical->r = sqrt(X * X + Y * Y + Z * Z);
2059            CoordSpherical->phig = TMath::RadToDeg()*asin(Z / (CoordSpherical->r));
2060            CoordSpherical->lambda = TMath::RadToDeg()*atan2(Y, X);
2061    } /*GM_CartesianToSpherical*/
Legend:
Removed from v.1.66  
changed lines
  Added in v.1.72
  ViewVC Help
Powered by ViewVC 1.1.23