/[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.69 by pam-mep, Fri Mar 28 20:47:15 2014 UTC revision 1.70 by mocchiut, Fri Apr 4 08:46:01 2014 UTC
# Line 368  int OrbitalInfoCore(UInt_t run, TFile *f Line 368  int OrbitalInfoCore(UInt_t run, TFile *f
368      //Geomagnetic coordinates calculations staff      //Geomagnetic coordinates calculations staff
369    
370    GMtype_CoordGeodetic location;    GMtype_CoordGeodetic location;
371    GMtype_CoordDipole GMlocation;    //  GMtype_CoordDipole GMlocation;
372    GMtype_Ellipsoid Ellip;    GMtype_Ellipsoid Ellip;
373    GMtype_Data G0, G1, H1;    GMtype_Data G0, G1, H1;
374                    
# Line 692  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      //      //
 <<<<<<< OrbitalInfoCore.cpp  
     if ( irun == 0 ){  
       if ( l0head->GetEntry(runinfo->EV_FROM) <= 0 ) throw -36;  
       //  
       // absolute time of first event of the run (it should not matter a lot)  
       //  
       ph = eh->GetPscuHeader();  
       atime = dbtime->DBabsTime(ph->GetOrbitalTime());  
         
       parerror=glparam->Query_GL_PARAM(atime-anni5,301,dbc); // parameters stored in DB in GL_PRAM table    
       if ( parerror<0 ) {  
         code = parerror;  
         goto closeandexit;  
     };  
       ltp1 = (Int_t)(glparam->PATH+glparam->NAME).Length();  
       if ( verbose ) printf(" Reading Earth's Magnetic Field parameter file: %s \n",(glparam->PATH+glparam->NAME).Data());  
       //  
       parerror=glparam2->Query_GL_PARAM(atime,301,dbc); // parameters stored in DB in GL_PRAM table    
       if ( parerror<0 ) {  
         code = parerror;  
         goto closeandexit;  
       };  
       ltp2 = (Int_t)(glparam2->PATH+glparam2->NAME).Length();  
       if ( verbose ) printf(" Reading Earth's Magnetic Field parameter file: %s \n",(glparam2->PATH+glparam2->NAME).Data());  
       //  
       parerror=glparam3->Query_GL_PARAM(atime,302,dbc); // parameters stored in DB in GL_PRAM table  
       if ( parerror<0 ) {  
         code = parerror;  
         goto closeandexit;  
       };  
       ltp3 = (Int_t)(glparam3->PATH+glparam3->NAME).Length();  
       if ( verbose ) printf(" Reading Earth's Magnetic Field parameter file: %s \n",(glparam3->PATH+glparam3->NAME).Data());  
       //  
       initize_((char *)(glparam->PATH+glparam->NAME).Data(),&ltp1,(char *)(glparam2->PATH+glparam2->NAME).Data(),&ltp2,(char *)(glparam3->PATH+glparam3->NAME).Data(),&ltp3);  
       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;  
       //  
 =======  
695      if ( !igrfloaded ){      if ( !igrfloaded ){
696    
697        if ( l0head->GetEntry(runinfo->EV_FROM) > 0 ){        if ( l0head->GetEntry(runinfo->EV_FROM) > 0 ){
# Line 765  int OrbitalInfoCore(UInt_t run, TFile *f Line 728  int OrbitalInfoCore(UInt_t run, TFile *f
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        }        }
 >>>>>>> 1.68  
733      }      }
734      //      //
735      // End IGRF stuff//      // End IGRF stuff//
# Line 1045  int OrbitalInfoCore(UInt_t run, TFile *f Line 1008  int OrbitalInfoCore(UInt_t run, TFile *f
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);            
           feldcof_(&jyear, &dimo); // get dipole moment for year  
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
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);            GM_TimeAdjustCoefs(year, jyear, G0, G1, H1, &Model);
# Line 1073  int OrbitalInfoCore(UInt_t run, TFile *f Line 1036  int OrbitalInfoCore(UInt_t run, TFile *f
1036          for ( ik = 0; ik < neventsm; ik++){  //number of macrocommad packets          for ( ik = 0; ik < neventsm; ik++){  //number of macrocommad packets
1037            if ( ch->GetEntry(ik) <= 0 ) throw -36;            if ( ch->GetEntry(ik) <= 0 ) throw -36;
1038            tmpSize = mcmdev->Records->GetEntries();            tmpSize = mcmdev->Records->GetEntries();
 <<<<<<< OrbitalInfoCore.cpp  
           numrec = tmpSize;  
           if ( debug ) cout << "packet number " << ik <<"\tnumber of subpackets is " << numrec << endl;  
 =======  
1039            //      numrec = tmpSize;            //      numrec = tmpSize;
1040  >>>>>>> 1.68            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
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]){
 <<<<<<< OrbitalInfoCore.cpp  
                         //if ( debug ) printf(" here1 %i \n",ui);  
 =======  
1050                        if ( debug ) printf(" here1 %i \n",ui);                        if ( debug ) printf(" here1 %i \n",ui);
 >>>>>>> 1.68  
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]){
 <<<<<<< OrbitalInfoCore.cpp  
1054                           // to avoid interpolation between bad quaternions arrays                           // to avoid interpolation between bad quaternions arrays
1055                           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){                           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){
1056                            Int_t sizeqmcmd = qtime.size();                            Int_t sizeqmcmd = qtime.size();
# Line 1114  int OrbitalInfoCore(UInt_t run, TFile *f Line 1068  int OrbitalInfoCore(UInt_t run, TFile *f
1068                            qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;                            qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;
1069                            qPitch[sizeqmcmd]=RYPang_upper->Tangazh;                            qPitch[sizeqmcmd]=RYPang_upper->Tangazh;
1070                           }                           }
 =======  
                         Int_t sizeqmcmd = qtime.size();  
                         inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);  
                         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;  
 >>>>>>> 1.68  
1071                        }                        }
1072                        for(Int_t mu = nt;mu<recSize;mu++){                        for(Int_t mu = nt;mu<recSize;mu++){
1073                          if(recqtime[mu]>lowerqtime && recqtime[mu]<u_time){                          if(recqtime[mu]>lowerqtime && recqtime[mu]<u_time){
# Line 1172  int OrbitalInfoCore(UInt_t run, TFile *f Line 1110  int OrbitalInfoCore(UInt_t run, TFile *f
1110                        }                        }
1111                      }                      }
1112                    }else{                    }else{
 <<<<<<< OrbitalInfoCore.cpp  
                         //if ( debug ) printf(" here2 %i \n",ui);  
 =======  
1113                      if ( debug ) printf(" here2 %i \n",ui);                      if ( debug ) printf(" here2 %i \n",ui);
 >>>>>>> 1.68  
1114                      Double_t u_time = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));                      Double_t u_time = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
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]){
 <<<<<<< OrbitalInfoCore.cpp  
1118                        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){                        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){
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);
# Line 1198  int OrbitalInfoCore(UInt_t run, TFile *f Line 1131  int OrbitalInfoCore(UInt_t run, TFile *f
1131                          qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;                          qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;
1132                          qPitch[sizeqmcmd]=RYPang_upper->Tangazh;                          qPitch[sizeqmcmd]=RYPang_upper->Tangazh;
1133                        }                        }
 =======  
                       Int_t sizeqmcmd = qtime.size();  
                       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;  
 >>>>>>> 1.68  
1134                      }                      }
1135                      for(Int_t mu = nt;mu<recSize;mu++){                      for(Int_t mu = nt;mu<recSize;mu++){
1136                        if(recqtime[mu]>lowerqtime && recqtime[mu]<u_time){                        if(recqtime[mu]>lowerqtime && recqtime[mu]<u_time){
# Line 1263  int OrbitalInfoCore(UInt_t run, TFile *f Line 1180  int OrbitalInfoCore(UInt_t run, TFile *f
1180            }            }
1181          }          }
1182                    
 <<<<<<< OrbitalInfoCore.cpp  
1183          if(qtime.size()==0){                            // in case if no orientation information in data          if(qtime.size()==0){                            // in case if no orientation information in data
1184             //if ( debug ) cout << "qtime.size() = 0" << endl;            if ( debug ) cout << "qtime.size() = 0" << endl;
             for(UInt_t my=0;my<recqtime.size();my++){  
               if(sqrt(pow(recq0[my],2)+pow(recq1[my],2)+pow(recq2[my],2)+pow(recq3[my],2))>0.99999){  
                 Int_t sizeqmcmd = qtime.size();  
                 inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);  
                 qtime[sizeqmcmd]=recqtime[my];  
                 q0[sizeqmcmd]=recq0[my];  
                 q1[sizeqmcmd]=recq1[my];  
                 q2[sizeqmcmd]=recq2[my];  
                 q3[sizeqmcmd]=recq3[my];  
                 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[my],recq1[my],recq2[my],recq3[my]);  
                 qRoll[sizeqmcmd]=RYPang_upper->Kren;  
                 qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;  
                 qPitch[sizeqmcmd]=RYPang_upper->Tangazh;  
               }  
             }  
 =======  
         if(qtime.size()==0){  
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            }            }
 >>>>>>> 1.68  
1202          }          }
1203                    
1204    
# Line 1360  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 ){
 <<<<<<< OrbitalInfoCore.cpp  
1261          if ( debug ) cout << "ok quaternions is exist and mu = " << must << endl;          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;          if ( debug ) cout << "qtimes[ " << qtime[0] << " , " << qtime[qtime.size()-1] << " ]\tatime = "<<atime<<endl;
       for(UInt_t mu = must;mu<qtime.size()-1;mu++){  
         //if ( debug ) printf(" ??grfuffi %i sixe %i must %i \n",mu,qtime.size()-1,must);  
         if(qtime[mu+1]>qtime[mu]){  
           if ( debug ) cout << "qtime[" << mu << "] = " << qtime[mu] << "\tqtime[" << mu+1 << "] = " << qtime[mu+1] << "\tatime = " << atime << endl;  
           if(atime<=qtime[mu+1] && atime>=qtime[mu]){  
             if ( debug ) cout << "here we have found proper quaternions for interpolation: mu = "<<mu<<endl;  
             must = mu;  
             incli = (qPitch[mu+1]-qPitch[mu])/(qtime[mu+1]-qtime[mu]);  
             orbitalinfo->theta =  incli*atime+qPitch[mu+1]-incli*qtime[mu+1];  
             incli = (qRoll[mu+1]-qRoll[mu])/(qtime[mu+1]-qtime[mu]);  
             orbitalinfo->etha =  incli*atime+qRoll[mu+1]-incli*qtime[mu+1];  
             incli = (qYaw[mu+1]-qYaw[mu])/(qtime[mu+1]-qtime[mu]);  
             orbitalinfo->phi =  incli*atime+qYaw[mu+1]-incli*qtime[mu+1];  
               
             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];  
             Float_t tg = (qtime[mu+1]-qtime[mu])/1000.;  
             if(tg>=1) tg=0.00;  
             orbitalinfo->TimeGap = TMath::Min(TMath::Abs(qtime[mu+1])-atime,TMath::Abs(atime-qtime[mu]))+tg;//qtime[mu+1]-qtime[mu];  
             orbitalinfo->mode = qmode[mu+1];  
   
             //if(atime==qtime[mu] || atime==qtime[mu+1]) orbitalinfo->qkind = 0; else orbitalinfo->qkind=1;  
             //if(qmode[mu+1]==-10) orbitalinfo->R10r = true;else orbitalinfo->R10r = false;  
   
             break;  
           }  
         }  
       }  
 =======  
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 1420  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            }            }
1297          }          }
 >>>>>>> 1.68  
1298        }        }
1299        if ( debug ) printf(" grfuffi5  \n");        if ( debug ) printf(" grfuffi5  \n");
1300        //        //
# Line 1497  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.;
 <<<<<<< OrbitalInfoCore.cpp  
1314          orbitalinfo->TimeGap = -1000.;          orbitalinfo->TimeGap = -1000.;
1315          //orbitalinfo->qkind = -1000;          //orbitalinfo->qkind = -1000;
1316        };          
1317          if ( debug ){          //      if ( debug ){
1318            Int_t lopu;          //        Int_t lopu;
1319           cin >> lopu;          //         cin >> lopu;
1320          }          //      }
 =======  
1321          if ( debug ) printf(" grfuffi6 \n");          if ( debug ) printf(" grfuffi6 \n");
1322        }        }
 >>>>>>> 1.68  
1323        //        //
1324        if ( debug ) printf(" filling \n");        if ( debug ) printf(" filling \n");
1325        // #########################################################################################################################          // #########################################################################################################################  
# Line 1521  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;
 <<<<<<< OrbitalInfoCore.cpp  
1335    
1336        cOrbit orbits2(*gltle->GetTle());        cOrbit orbits2(*gltle->GetTle());
1337        orbits2.getPosition((double) (atime - gltle->GetFromTime())/60., &eCi);        orbits2.getPosition((double) (atime - gltle->GetFromTime())/60., &eCi);
1338        Float_t x=eCi.getPos().m_x;        //      Float_t x=eCi.getPos().m_x;
1339        Float_t y=eCi.getPos().m_y;        //      Float_t y=eCi.getPos().m_y;
1340        Float_t z=eCi.getPos().m_z;        //      Float_t z=eCi.getPos().m_z;
1341          
1342        TVector3 V(eCi.getVel().m_x,eCi.getVel().m_y,eCi.getVel().m_z);        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);        TVector3 Pos(eCi.getPos().m_x,eCi.getPos().m_y,eCi.getPos().m_z);
1344          
1345        Float_t dlon=Pos.Phi()*TMath::RadToDeg()-lon;        Float_t dlon=Pos.Phi()*TMath::RadToDeg()-lon;
1346          
1347        Pos.RotateZ(-dlon*TMath::DegToRad());        Pos.RotateZ(-dlon*TMath::DegToRad());
1348        V.RotateZ(-dlon*TMath::DegToRad());        V.RotateZ(-dlon*TMath::DegToRad());
1349        Float_t diro;        Float_t diro;
1350        if(V.Z()>0) diro=1; else diro=-1;        if(V.Z()>0) diro=1; else diro=-1;
1351          
1352  // 10REDNEW        // 10REDNEW
1353        Int_t errq=0;        Int_t errq=0;
1354        Int_t azim=0;;        Int_t azim=0;;
1355       for(UInt_t mu = must;mu<RTtime2.size()-1;mu++){        for(UInt_t mu = must;mu<RTtime2.size()-1;mu++){
1356         if(atime<=RTtime2[mu] && atime>=RTtime1[mu]){          if(atime<=RTtime2[mu] && atime>=RTtime1[mu]){
1357           errq=RTerrq[mu];            errq=RTerrq[mu];
1358           azim=RTazim[mu];            azim=RTazim[mu];
1359         }          }
1360       }        }
1361       orbitalinfo->errq = errq;        orbitalinfo->errq = errq;
1362       orbitalinfo->azim = azim;        orbitalinfo->azim = azim;
1363       orbitalinfo->qkind = 0;        orbitalinfo->qkind = 0;
1364          
      if( lon<180 && lon>-180 && lat<90 && lat>-90 && alt>0 ){    
 =======  
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 ){  
 >>>>>>> 1.68  
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;          orbitalinfo->V = V;
1372    
1373          GMtype_CoordGeodetic  location;          //      GMtype_CoordGeodetic  location;
1374          location.lambda = lon;          location.lambda = lon;
1375          location.phi = lat;          location.phi = lat;
1376          location.HeightAboveEllipsoid = alt;          location.HeightAboveEllipsoid = alt;
Legend:
Removed from v.1.69  
changed lines
  Added in v.1.70
  ViewVC Help
Powered by ViewVC 1.1.23