/[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.84 by mocchiut, Thu Mar 26 14:55:38 2015 UTC revision 1.93 by mayorov, Thu Dec 27 13:12:28 2018 UTC
# Line 90  int OrbitalInfoCore(UInt_t run, TFile *f Line 90  int OrbitalInfoCore(UInt_t run, TFile *f
90      i = 0;      i = 0;
91      while ( i < OrbitalInfoargc ){      while ( i < OrbitalInfoargc ){
92        if ( !strcmp(OrbitalInfoargv[i],"-processFolder") ) {        if ( !strcmp(OrbitalInfoargv[i],"-processFolder") ) {
93          if ( OrbitalInfoargc < i+1 ){        if ( OrbitalInfoargc < i+1 ){
94            throw -3;      throw -3;
95          }    }
96          processFolder = (TString)OrbitalInfoargv[i+1];    processFolder = (TString)OrbitalInfoargv[i+1];
97          i++;    i++;
98        }        }
99        if ( (!strcmp(OrbitalInfoargv[i],"--debug")) || (!strcmp(OrbitalInfoargv[i],"-g")) ) {        if ( (!strcmp(OrbitalInfoargv[i],"--debug")) || (!strcmp(OrbitalInfoargv[i],"-g")) ) {
100          verbose = true;          verbose = true;
101          debug = true;          debug = true;
102        }        }
103        if ( (!strcmp(OrbitalInfoargv[i],"--verbose")) || (!strcmp(OrbitalInfoargv[i],"-v")) ) {        if ( (!strcmp(OrbitalInfoargv[i],"--verbose")) || (!strcmp(OrbitalInfoargv[i],"-v")) ) {
104          verbose = true;          verbose = true;
105        }        }
106        if ( (!strcmp(OrbitalInfoargv[i],"--standalone")) ) {        if ( (!strcmp(OrbitalInfoargv[i],"--standalone")) ) {
107          standalone = true;          standalone = true;
108        }        }
109        if ( (!strcmp(OrbitalInfoargv[i],"--calculate-pitch")) ) {        if ( (!strcmp(OrbitalInfoargv[i],"--calculate-pitch")) ) {
110          standalone = false;          standalone = false;
111        }        }
112        i++;        i++;
113      }      }
# Line 198  int OrbitalInfoCore(UInt_t run, TFile *f Line 198  int OrbitalInfoCore(UInt_t run, TFile *f
198    Float_t dimo = 0.0; // dipole moment (computed from dat files) // EM GCC 4.7    Float_t dimo = 0.0; // dipole moment (computed from dat files) // EM GCC 4.7
199    Float_t bnorth, beast, bdown, babs;    Float_t bnorth, beast, bdown, babs;
200    Float_t xl; // L value    Float_t xl; // L value
201    Float_t icode; // code value for L accuracy (see fortran code)    Int_t icode; // code value for L accuracy (see fortran code)
202    Float_t bab1; // What's  the difference with babs?    Float_t bab1; // What's  the difference with babs?
203    Float_t stps = 0.005; // step size for field line tracing    Float_t stps = 0.005; // step size for field line tracing
204    Float_t bdel = 0.01; // required accuracy    Float_t bdel = 0.01; // required accuracy
# Line 263  int OrbitalInfoCore(UInt_t run, TFile *f Line 263  int OrbitalInfoCore(UInt_t run, TFile *f
263    //  Double_t A2;    //  Double_t A2;
264    //  Double_t A3;    //  Double_t A3;
265    Double_t Px = 0;    Double_t Px = 0;
266    Double_t Py = 0;          Double_t Py = 0;  
267    Double_t Pz = 0;      Double_t Pz = 0;  
268    TTree *ttof = 0;    TTree *ttof = 0;
269    ToFLevel2 *tof = new ToFLevel2();    ToFLevel2 *tof = new ToFLevel2();
# Line 419  int OrbitalInfoCore(UInt_t run, TFile *f Line 419  int OrbitalInfoCore(UInt_t run, TFile *f
419    //  GMtype_CoordDipole GMlocation;    //  GMtype_CoordDipole GMlocation;
420    GMtype_Ellipsoid Ellip;    GMtype_Ellipsoid Ellip;
421    GMtype_Data G0, G1, H1;    GMtype_Data G0, G1, H1;
422              
423    //  { // 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    //  {  // 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
424    //    TString igpath="/data03/Malakhov/pam9Malakhov/installed10/calib/orb-param/";    //    TString igpath="/data03/Malakhov/pam9Malakhov/installed10/calib/orb-param/";
425    //  }    //  }
426    
# Line 706  int OrbitalInfoCore(UInt_t run, TFile *f Line 706  int OrbitalInfoCore(UInt_t run, TFile *f
706      //            //      
707      if ( nobefrun > 0 ){      if ( nobefrun > 0 ){
708        if (verbose){        if (verbose){
709          printf("\n Pre-processing: copying events from the old tree before the processed run\n");            printf("\n Pre-processing: copying events from the old tree before the processed run\n");  
710          printf(" Copying %u events in the file which are before the beginning of the run %u \n",nobefrun,run);          printf(" Copying %u events in the file which are before the beginning of the run %u \n",nobefrun,run);
711          printf(" Start copying at event number 0, end copying at event number %u \n",nobefrun);          printf(" Start copying at event number 0, end copying at event number %u \n",nobefrun);
712        }        }
713        for (UInt_t j = 0; j < nobefrun; j++){        for (UInt_t j = 0; j < nobefrun; j++){
714          //          //
715          if ( OrbitalInfotrclone->GetEntry(j) <= 0 ) throw -36;              if ( OrbitalInfotrclone->GetEntry(j) <= 0 ) throw -36;    
716          //          //
717          // copy orbitalinfoclone to mydec          // copy orbitalinfoclone to mydec
718          //          //
719          //      orbitalinfo->Clear();          //  orbitalinfo->Clear();
720          //          //
721          memcpy(&orbitalinfo,&orbitalinfoclone,sizeof(orbitalinfoclone));          memcpy(&orbitalinfo,&orbitalinfoclone,sizeof(orbitalinfoclone));
722          //          //
723          // Fill entry in the new tree          // Fill entry in the new tree
724          //          //
725          OrbitalInfotr->Fill();          OrbitalInfotr->Fill();
726          //          //
727        };        };
728        if (verbose) printf(" Finished successful copying!\n");        if (verbose) printf(" Finished successful copying!\n");
729      };      };
# Line 879  int OrbitalInfoCore(UInt_t run, TFile *f Line 879  int OrbitalInfoCore(UInt_t run, TFile *f
879        Row = pResult->Next();        Row = pResult->Next();
880        //        //
881        while ( Row ){        while ( Row ){
882          //          //
883          // store infos and exit          // store infos and exit
884          //          //
885          l0fid[i] = (UInt_t)atoll(Row->GetField(0));          l0fid[i] = (UInt_t)atoll(Row->GetField(0));
886          i--;          i--;
887          if (Row){  // memleak!          if (Row){  // memleak!
888            delete Row;            delete Row;
889            Row = 0;            Row = 0;
890          }          }
891          Row = pResult->Next();            Row = pResult->Next();  
892          //        //
893        }        }
894        if (Row) delete Row;        if (Row) delete Row;
895        pResult->Delete();        pResult->Delete();
# Line 906  int OrbitalInfoCore(UInt_t run, TFile *f Line 906  int OrbitalInfoCore(UInt_t run, TFile *f
906        Row = pResult->Next();        Row = pResult->Next();
907        //        //
908        while ( Row ){        while ( Row ){
909          //          //
910          // store infos and exit          // store infos and exit
911          //          //
912          l0fid[i] = (UInt_t)atoll(Row->GetField(0));          l0fid[i] = (UInt_t)atoll(Row->GetField(0));
913          i++;          i++;
914          if (Row){  // memleak!          if (Row){  // memleak!
915            delete Row;            delete Row;
916            Row = 0;            Row = 0;
917          }          }
918          Row = pResult->Next();            Row = pResult->Next();  
919          //          //
920        }        }
921        if (Row) delete Row;        if (Row) delete Row;
922        pResult->Delete();        pResult->Delete();
# Line 926  int OrbitalInfoCore(UInt_t run, TFile *f Line 926  int OrbitalInfoCore(UInt_t run, TFile *f
926      UInt_t previd = 0;      UInt_t previd = 0;
927      while ( i < 10 ){      while ( i < 10 ){
928        if ( l0fid[i] && previd != l0fid[i] ){        if ( l0fid[i] && previd != l0fid[i] ){
929          previd = l0fid[i];          previd = l0fid[i];
930          myquery.str("");          myquery.str("");
931          myquery << "select PATH,NAME from GL_ROOT where ID=" << l0fid[i] << " ;";          myquery << "select PATH,NAME from GL_ROOT where ID=" << l0fid[i] << " ;";
932          //          //
933          pResult = dbc->Query(myquery.str().c_str());          pResult = dbc->Query(myquery.str().c_str());
934          //          //
935          if( pResult ){          if( pResult ){
936            //            //
937            Row = pResult->Next();            Row = pResult->Next();
938            //            //
939            if ( debug ) printf(" Using inclination informations from file: %s \n",(((TString)gSystem->ExpandPathName(Row->GetField(0)))+"/"+(TString)Row->GetField(1)).Data());            if ( debug ) printf(" Using inclination informations from file: %s \n",(((TString)gSystem->ExpandPathName(Row->GetField(0)))+"/"+(TString)Row->GetField(1)).Data());
940            ch->Add(((TString)gSystem->ExpandPathName(Row->GetField(0)))+"/"+(TString)Row->GetField(1));            ch->Add(((TString)gSystem->ExpandPathName(Row->GetField(0)))+"/"+(TString)Row->GetField(1));
941            //            //
942            if (Row) delete Row;            if (Row) delete Row;
943            pResult->Delete();            pResult->Delete();
944          }          }
945        }        }
946        i++;        i++;
947      }      }
# Line 1001  int OrbitalInfoCore(UInt_t run, TFile *f Line 1001  int OrbitalInfoCore(UInt_t run, TFile *f
1001        //for ( re = runinfo->EV_FROM; re < (runinfo->EV_FROM+10); re++){        //for ( re = runinfo->EV_FROM; re < (runinfo->EV_FROM+10); re++){
1002    
1003        //        //
1004        if ( procev%1000 == 0 && procev > 0 && verbose ) printf(" %iK \n",procev/1000);          if ( procev%1000 == 0 && procev > 0 && verbose ) printf(" %iK \n",procev/1000);  
1005        if ( debug ) printf(" %i \n",procev);              if ( debug ) printf(" %i \n",procev);      
1006        //        //
1007        if ( l0head->GetEntry(re) <= 0 ) throw -36;        if ( l0head->GetEntry(re) <= 0 ) throw -36;
# Line 1015  int OrbitalInfoCore(UInt_t run, TFile *f Line 1015  int OrbitalInfoCore(UInt_t run, TFile *f
1015        // paranoid check        // paranoid check
1016        //        //
1017        if ( (atime > (runinfo->RUNTRAILER_TIME+1)) || (atime < (runinfo->RUNHEADER_TIME-1))  ) {        if ( (atime > (runinfo->RUNTRAILER_TIME+1)) || (atime < (runinfo->RUNHEADER_TIME-1))  ) {
1018          if (verbose) printf(" OrbitalInfo - WARNING: event at time outside the run time window, skipping it\n");          if (verbose) printf(" OrbitalInfo - WARNING: event at time outside the run time window, skipping it\n");
1019          jumped++;          jumped++;
1020          //      debug = true;          //  debug = true;
1021          continue;          continue;
1022        }        }
1023    
1024        // just for testing:        // just for testing:
# Line 1094  int OrbitalInfoCore(UInt_t run, TFile *f Line 1094  int OrbitalInfoCore(UInt_t run, TFile *f
1094        // retrieve tof informations        // retrieve tof informations
1095        //        //
1096        if ( !reprocall ){        if ( !reprocall ){
1097          itr = nobefrun + (re - evfrom - jumped);          itr = nobefrun + (re - evfrom - jumped);
1098          //itr = re-(46438+200241);          //itr = re-(46438+200241);
1099        } else {        } else {
1100          itr = runinfo->GetFirstEntry() + (re - evfrom - jumped);          itr = runinfo->GetFirstEntry() + (re - evfrom - jumped);
1101        };        };
1102        //        //
1103        if ( !standalone ){        if ( !standalone ){
1104          if ( itr > nevtofl2 ){            if ( itr > nevtofl2 ){  
1105            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);
1106            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);
1107            l0File->Close();            l0File->Close();
1108            code = -904;            code = -904;
1109            goto closeandexit;            goto closeandexit;
1110          };          };
1111          //          //
1112          tof->Clear();          tof->Clear();
1113          //          //
1114          // Clones array must be cleared before going on          // Clones array must be cleared before going on
1115          //          //
1116          if ( hasNucleiTof ){          if ( hasNucleiTof ){
# Line 1122  int OrbitalInfoCore(UInt_t run, TFile *f Line 1122  int OrbitalInfoCore(UInt_t run, TFile *f
1122          if ( hasExtTof ){          if ( hasExtTof ){
1123            tcExtTof->Delete();            tcExtTof->Delete();
1124          }          }
1125          //          //
1126          if ( verbose ) printf(" get tof tree entries... entry = %i in Level2 file\n",itr);          if ( verbose ) printf(" get tof tree entries... entry = %i in Level2 file\n",itr);
1127          if ( ttof->GetEntry(itr) <= 0 ){          if ( ttof->GetEntry(itr) <= 0 ){
1128            if ( verbose ) printf(" problems with tof tree entries... entry = %i in Level2 file\n",itr);            if ( verbose ) printf(" problems with tof tree entries... entry = %i in Level2 file\n",itr);
1129            if ( verbose ) printf(" nobefrun %u re %u evfrom %u jumped %u reprocall %i \n",nobefrun,re,evfrom,jumped,reprocall);            if ( verbose ) printf(" nobefrun %u re %u evfrom %u jumped %u reprocall %i \n",nobefrun,re,evfrom,jumped,reprocall);
1130            throw -36;            throw -36;
1131          }          }
1132          if ( verbose ) printf(" gat0\n");          if ( verbose ) printf(" gat0\n");
1133          //          //
1134        }        }
1135        //        //
1136        // retrieve tracker informations        // retrieve tracker informations
1137        //        //
1138        if ( !standalone ){        if ( !standalone ){
1139          if ( itr > nevtrkl2 ){            if ( itr > nevtrkl2 ){  
1140            if ( verbose ) printf(" OrbitalInfo - ERROR: no trk events with entry = %i in Level2 file\n",itr);            if ( verbose ) printf(" OrbitalInfo - ERROR: no trk events with entry = %i in Level2 file\n",itr);
1141            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);
1142            l0File->Close();            l0File->Close();
1143            code = -905;            code = -905;
1144            goto closeandexit;            goto closeandexit;
1145          }          }
1146          //          //
1147          if ( verbose ) printf(" gat1\n");          if ( verbose ) printf(" gat1\n");
1148          trke->Clear();          trke->Clear();
1149          //          //
1150          // Clones array must be cleared before going on          // Clones array must be cleared before going on
1151          //          //
1152          if ( hasNucleiTrk ){          if ( hasNucleiTrk ){
# Line 1167  int OrbitalInfoCore(UInt_t run, TFile *f Line 1167  int OrbitalInfoCore(UInt_t run, TFile *f
1167            if ( verbose ) printf(" gat7\n");            if ( verbose ) printf(" gat7\n");
1168            torbExtTrk->Delete();            torbExtTrk->Delete();
1169          }          }
1170          //          //
1171          if ( verbose ) printf(" get trk tree entries... entry = %i in Level2 file\n",itr);          if ( verbose ) printf(" get trk tree entries... entry = %i in Level2 file\n",itr);
1172          if ( ttrke->GetEntry(itr) <= 0 ) throw -36;          if ( ttrke->GetEntry(itr) <= 0 ) throw -36;
1173          //          //
1174        }        }
1175    
1176        //        //
# Line 1180  int OrbitalInfoCore(UInt_t run, TFile *f Line 1180  int OrbitalInfoCore(UInt_t run, TFile *f
1180        //        //
1181        if ( debug ) printf(" %i start processing \n",procev);              if ( debug ) printf(" %i start processing \n",procev);      
1182        orbitalinfo->Clear();        orbitalinfo->Clear();
1183    
1184        //        //
1185        OrbitalInfoTrkVar *t_orb = new OrbitalInfoTrkVar();        OrbitalInfoTrkVar *t_orb = new OrbitalInfoTrkVar();
1186        if( !(orbitalinfo->OrbitalInfoTrk) ) orbitalinfo->OrbitalInfoTrk = new TClonesArray("OrbitalInfoTrkVar",2);        if( !(orbitalinfo->OrbitalInfoTrk) ) orbitalinfo->OrbitalInfoTrk = new TClonesArray("OrbitalInfoTrkVar",2);
# Line 1203  int OrbitalInfoCore(UInt_t run, TFile *f Line 1204  int OrbitalInfoCore(UInt_t run, TFile *f
1204        //        //
1205        if ( debug ) printf(" %i sgp4 \n",procev);              if ( debug ) printf(" %i sgp4 \n",procev);      
1206        cCoordGeo coo;        cCoordGeo coo;
1207        Float_t jyear=0.;            Float_t jyear=0.;
1208        //        //
1209        if(atime >= gltle->GetToTime() || atime < gltle->GetFromTime() ) {  // AGH! bug when reprocessing??        if(atime >= gltle->GetToTime() || atime < gltle->GetFromTime() ) {  // AGH! bug when reprocessing??
1210          if ( !gltle->Query(atime, dbc) ){  
1211            //                if ( !gltle->Query(atime, dbc) ){
1212            // Compute the magnetic dipole moment.            //    
1213            //            // Compute the magnetic dipole moment.
1214            if ( debug ) printf(" %i compute magnetic dipole moment \n",procev);                  //
1215            UInt_t year, month, day, hour, min, sec;            if ( debug ) printf(" %i compute magnetic dipole moment \n",procev);      
1216            //            UInt_t year, month, day, hour, min, sec;
1217            TTimeStamp t = TTimeStamp(atime, kTRUE);            //
1218            t.GetDate(kTRUE, 0, &year, &month, &day);            TTimeStamp t = TTimeStamp(atime, kTRUE);
1219            t.GetTime(kTRUE, 0, &hour, &min, &sec);            t.GetDate(kTRUE, 0, &year, &month, &day);
1220            jyear = (float) year            t.GetTime(kTRUE, 0, &hour, &min, &sec);
1221              + (month*31.+ (float) day)/365.            jyear = (float) year
1222              + (hour*3600.+min*60.+(float)sec)/(24.*3600.*365.);              + (month*31.+ (float) day)/365.
1223            //              + (hour*3600.+min*60.+(float)sec)/(24.*3600.*365.);
1224            if ( debug ) printf(" %i compute magnetic dipole moment get dipole moment for year\n",procev);                        //
1225            if ( debug ) printf(" %i jyear %f dimo %f \n",procev,jyear,dimo);            if ( debug ) printf(" %i compute magnetic dipole moment get dipole moment for year\n",procev);            
1226            feldcof_(&jyear, &dimo); // get dipole moment for year            if ( debug ) printf(" %i jyear %f dimo %f \n",procev,jyear,dimo);
1227            if ( debug ) printf(" %i compute magnetic dipole moment end\n",procev);            feldcof_(&jyear, &dimo); // get dipole moment for year
1228              if ( debug ) printf(" %i compute magnetic dipole moment end\n",procev);
1229            //      GM_TimeAdjustCoefs(year, jyear, G0, G1, H1, &Model);  
1230            GM_TimeAdjustCoefs(GM_STARTYEAR, (jyear-currentYear+GM_STARTYEAR), G0, G1, H1, &Model);  // EM: input this way due to the new way of storing data into Gn,H1 and to avoid changing GM_Time...            //    GM_TimeAdjustCoefs(year, jyear, G0, G1, H1, &Model);
1231            GM_PoleLocation(Model, &Pole);            GM_TimeAdjustCoefs(GM_STARTYEAR, (jyear-currentYear+GM_STARTYEAR), G0, G1, H1, &Model);  // EM: input this way due to the new way of storing data into Gn,H1 and to avoid changing GM_Time...
1232                        GM_PoleLocation(Model, &Pole);
1233          } else {          } else {
1234            code = -56;            code = -56;
1235            goto closeandexit;            goto closeandexit;
1236          };          };
1237        }        }
1238        coo = getCoo(atime, gltle->GetFromTime(), gltle->GetTle());        coo = getCoo(atime, gltle->GetFromTime(), gltle->GetTle());
1239        //        //
# Line 1241  int OrbitalInfoCore(UInt_t run, TFile *f Line 1242  int OrbitalInfoCore(UInt_t run, TFile *f
1242        // synchronize with quaternions data        // synchronize with quaternions data
1243        //        //
1244        if ( isf && neventsm>0 ){        if ( isf && neventsm>0 ){
1245          //          //
1246          // First event          // First event
1247          //          //
1248          isf = false;          isf = false;
1249          //      upperqtime = atime;          //  upperqtime = atime;
1250          lowerqtime = runinfo->RUNHEADER_TIME;          lowerqtime = runinfo->RUNHEADER_TIME;
1251          for ( ik = 0; ik < neventsm; ik++){  //number of macrocommad packets          for ( ik = 0; ik < neventsm; ik++){  //number of macrocommad packets
1252            if ( ch->GetEntry(ik) <= 0 ) throw -36;            if ( ch->GetEntry(ik) <= 0 ) throw -36;
1253            tmpSize = mcmdev->Records->GetEntries();            tmpSize = mcmdev->Records->GetEntries();
1254            //      numrec = tmpSize;            //    numrec = tmpSize;
1255            if ( debug ) cout << "packet number " << ik <<"\tnumber of subpackets is " << tmpSize << endl;            if ( debug ) cout << "packet number " << ik <<"\tnumber of subpackets is " << tmpSize << endl;
1256            for (Int_t j3 = 0;j3<tmpSize;j3++){  //number of subpackets              for (Int_t j3 = 0;j3<tmpSize;j3++){  //number of subpackets
1257              mcmdrc = (pamela::McmdRecord*)mcmdev->Records->At(j3);                mcmdrc = (pamela::McmdRecord*)mcmdev->Records->At(j3);
1258              if ( mcmdrc ){ // missing inclination bug [8RED 090116]                if ( mcmdrc ){ // missing inclination bug [8RED 090116]
1259                if ( debug ) printf(" pluto \n");                  if ( debug ) printf(" pluto \n");
1260                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
1261                 L_QQ_Q_l_upper->fill(mcmdrc->McmdData);                    L_QQ_Q_l_upper->fill(mcmdrc->McmdData);
1262                  for (UInt_t ui = 0; ui < 6; ui++){                    for (UInt_t ui = 0; ui < 6; ui++){
1263                    if (ui>0){                      if (ui>0){
1264                      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]){
1265                        if ( debug ) printf(" here1 %i \n",ui);                          Double_t u_time = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[ui]*1000-DeltaOBT*1000));
1266                        Double_t u_time = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[ui]*1000-DeltaOBT*1000));                          Int_t recSize = recqtime.size();
1267                        Int_t recSize = recqtime.size();                          if(lowerqtime > recqtime[recSize-1]){
1268                        if(lowerqtime > recqtime[recSize-1]){                            // to avoid interpolation between bad quaternions arrays
1269                           // to avoid interpolation between bad quaternions arrays                            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){
1270                           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){                              Int_t sizeqmcmd = qtime.size();
1271                            Int_t sizeqmcmd = qtime.size();                              inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);
1272                            inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);                              qtime[sizeqmcmd]=u_time;
1273                            qtime[sizeqmcmd]=u_time;                              q0[sizeqmcmd]=L_QQ_Q_l_upper->quat[ui][0];
1274                            q0[sizeqmcmd]=L_QQ_Q_l_upper->quat[ui][0];                              q1[sizeqmcmd]=L_QQ_Q_l_upper->quat[ui][1];
1275                            q1[sizeqmcmd]=L_QQ_Q_l_upper->quat[ui][1];                              q2[sizeqmcmd]=L_QQ_Q_l_upper->quat[ui][2];
1276                            q2[sizeqmcmd]=L_QQ_Q_l_upper->quat[ui][2];                              q3[sizeqmcmd]=L_QQ_Q_l_upper->quat[ui][3];
1277                            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);
1278                            qmode[sizeqmcmd]=holeq(lowerqtime,qtime[sizeqmcmd],L_QQ_Q_l_lower,L_QQ_Q_l_upper,ui);                              lowerqtime = u_time;
1279                            lowerqtime = u_time;                              orbits.getPosition((double) (u_time - gltle->GetFromTime())/60., &eCi);
1280                            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]);
1281                            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;
1282                            qRoll[sizeqmcmd]=RYPang_upper->Kren;                              qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;
1283                            qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;                              qPitch[sizeqmcmd]=RYPang_upper->Tangazh;
1284                            qPitch[sizeqmcmd]=RYPang_upper->Tangazh;                            }
1285                           }                          }
1286                        }                          for(Int_t mu = nt;mu<recSize;mu++){
1287                        for(Int_t mu = nt;mu<recSize;mu++){                            if(recqtime[mu]>lowerqtime && recqtime[mu]<u_time){
1288                          if(recqtime[mu]>lowerqtime && recqtime[mu]<u_time){                              if(sqrt(pow(recq0[mu],2)+pow(recq1[mu],2)+pow(recq2[mu],2)+pow(recq3[mu],2))>0.99999){
1289                            if(sqrt(pow(recq0[mu],2)+pow(recq1[mu],2)+pow(recq2[mu],2)+pow(recq3[mu],2))>0.99999){                                nt=mu;
1290                              nt=mu;                                Int_t sizeqmcmd = qtime.size();
1291                              Int_t sizeqmcmd = qtime.size();                                inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);
1292                              inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);                                qtime[sizeqmcmd]=recqtime[mu];
1293                              qtime[sizeqmcmd]=recqtime[mu];                                q0[sizeqmcmd]=recq0[mu];
1294                              q0[sizeqmcmd]=recq0[mu];                                q1[sizeqmcmd]=recq1[mu];
1295                              q1[sizeqmcmd]=recq1[mu];                                q2[sizeqmcmd]=recq2[mu];
1296                              q2[sizeqmcmd]=recq2[mu];                                q3[sizeqmcmd]=recq3[mu];
1297                              q3[sizeqmcmd]=recq3[mu];                                qmode[sizeqmcmd]=-10;
1298                              qmode[sizeqmcmd]=-10;                                orbits.getPosition((double) (qtime[sizeqmcmd] - gltle->GetFromTime())/60., &eCi);
1299                              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]);
1300                              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;
1301                              qRoll[sizeqmcmd]=RYPang_upper->Kren;                                qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;
1302                              qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;                                qPitch[sizeqmcmd]=RYPang_upper->Tangazh;
1303                              qPitch[sizeqmcmd]=RYPang_upper->Tangazh;                              }
1304                            }                            }
1305                          }                            if(recqtime[mu]>=u_time){
1306                          if(recqtime[mu]>=u_time){                              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){
1307                            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){                                Int_t sizeqmcmd = qtime.size();
1308                              Int_t sizeqmcmd = qtime.size();                                inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);
1309                              inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);                                qtime[sizeqmcmd]=u_time;
1310                              qtime[sizeqmcmd]=u_time;                                q0[sizeqmcmd]=L_QQ_Q_l_upper->quat[ui][0];
1311                              q0[sizeqmcmd]=L_QQ_Q_l_upper->quat[ui][0];                                q1[sizeqmcmd]=L_QQ_Q_l_upper->quat[ui][1];
1312                              q1[sizeqmcmd]=L_QQ_Q_l_upper->quat[ui][1];                                q2[sizeqmcmd]=L_QQ_Q_l_upper->quat[ui][2];
1313                              q2[sizeqmcmd]=L_QQ_Q_l_upper->quat[ui][2];                                q3[sizeqmcmd]=L_QQ_Q_l_upper->quat[ui][3];
1314                              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);
1315                              qmode[sizeqmcmd]=holeq(lowerqtime,qtime[sizeqmcmd],L_QQ_Q_l_lower,L_QQ_Q_l_upper,ui);                                lowerqtime = u_time;
1316                              lowerqtime = u_time;                                orbits.getPosition((double) (u_time - gltle->GetFromTime())/60., &eCi);
1317                              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]);
1318                              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;
1319                              qRoll[sizeqmcmd]=RYPang_upper->Kren;                                qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;
1320                              qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;                                qPitch[sizeqmcmd]=RYPang_upper->Tangazh;
1321                              qPitch[sizeqmcmd]=RYPang_upper->Tangazh;                                break;
1322                              break;                              }
1323                            }                            }
1324                          }                          }
1325                        }                        }
1326                      }                      }else{
1327                    }else{                        //if ( debug ) printf(" here2 %i \n",ui);
1328                      if ( debug ) printf(" here2 %i \n",ui);                        Double_t u_time = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
1329                      Double_t u_time = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));                        if(lowerqtime>u_time)nt=0;
1330                      if(lowerqtime>u_time)nt=0;                        Int_t recSize = recqtime.size();
1331                      Int_t recSize = recqtime.size();                        if(lowerqtime > recqtime[recSize-1]){
1332                      if(lowerqtime > recqtime[recSize-1]){                          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){
1333                        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){                            Int_t sizeqmcmd = qtime.size();
1334                          Int_t sizeqmcmd = qtime.size();                            inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);
1335                          inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);                            qtime[sizeqmcmd]=u_time;
1336                          qtime[sizeqmcmd]=u_time;                            q0[sizeqmcmd]=L_QQ_Q_l_upper->quat[0][0];
1337                          q0[sizeqmcmd]=L_QQ_Q_l_upper->quat[0][0];                            q1[sizeqmcmd]=L_QQ_Q_l_upper->quat[0][1];
1338                          q1[sizeqmcmd]=L_QQ_Q_l_upper->quat[0][1];                            q2[sizeqmcmd]=L_QQ_Q_l_upper->quat[0][2];
1339                          q2[sizeqmcmd]=L_QQ_Q_l_upper->quat[0][2];                            q3[sizeqmcmd]=L_QQ_Q_l_upper->quat[0][3];
1340                          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);
1341                          qmode[sizeqmcmd]=holeq(lowerqtime,qtime[sizeqmcmd],L_QQ_Q_l_lower,L_QQ_Q_l_upper,ui);                            lowerqtime = u_time;
1342                          lowerqtime = u_time;                            orbits.getPosition((double) (u_time - gltle->GetFromTime())/60., &eCi);
1343                          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]);
1344                          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;
1345                          qRoll[sizeqmcmd]=RYPang_upper->Kren;                            qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;
1346                          qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;                            qPitch[sizeqmcmd]=RYPang_upper->Tangazh;
1347                          qPitch[sizeqmcmd]=RYPang_upper->Tangazh;                          }
1348                        }                        }
1349                      }                        for(Int_t mu = nt;mu<recSize;mu++){
1350                      for(Int_t mu = nt;mu<recSize;mu++){                          if(recqtime[mu]>lowerqtime && recqtime[mu]<u_time){
1351                        if(recqtime[mu]>lowerqtime && recqtime[mu]<u_time){                           if(sqrt(pow(recq0[mu],2)+pow(recq1[mu],2)+pow(recq2[mu],2)+pow(recq3[mu],2))>0.99999){
1352                           if(sqrt(pow(recq0[mu],2)+pow(recq1[mu],2)+pow(recq2[mu],2)+pow(recq3[mu],2))>0.99999){                             Int_t sizeqmcmd = qtime.size();
1353  //                         nt=mu;                             inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);
1354                             Int_t sizeqmcmd = qtime.size();                             qtime[sizeqmcmd]=recqtime[mu];
1355                             inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);                             q0[sizeqmcmd]=recq0[mu];
1356                             qtime[sizeqmcmd]=recqtime[mu];                             q1[sizeqmcmd]=recq1[mu];
1357                             q0[sizeqmcmd]=recq0[mu];                             q2[sizeqmcmd]=recq2[mu];
1358                             q1[sizeqmcmd]=recq1[mu];                             q3[sizeqmcmd]=recq3[mu];
1359                             q2[sizeqmcmd]=recq2[mu];                             qmode[sizeqmcmd]=-10;
1360                             q3[sizeqmcmd]=recq3[mu];                             orbits.getPosition((double) (qtime[sizeqmcmd] - gltle->GetFromTime())/60., &eCi);
1361                             qmode[sizeqmcmd]=-10;                             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]);
1362                             orbits.getPosition((double) (qtime[sizeqmcmd] - gltle->GetFromTime())/60., &eCi);                             qRoll[sizeqmcmd]=RYPang_upper->Kren;
1363                             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]);                             qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;
1364                             qRoll[sizeqmcmd]=RYPang_upper->Kren;                             qPitch[sizeqmcmd]=RYPang_upper->Tangazh;
1365                             qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;                           }
1366                             qPitch[sizeqmcmd]=RYPang_upper->Tangazh;                          }
1367  /* CHECK RECOVERED QUATERNIONS PROBLEM */                          if(recqtime[mu]>=u_time){
1368  if(recqtime[mu]>=1160987921.75 && recqtime[mu]<=1160987932.00){                            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){
1369    cout<<"We found it while checking all quaternions"<<"\t"<<recqtime[mu]<<endl;                              Int_t sizeqmcmd = qtime.size();
1370  }                              inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);
1371                           }                              qtime[sizeqmcmd]=u_time;
1372                        }                              q0[sizeqmcmd]=L_QQ_Q_l_upper->quat[0][0];
1373                        if(recqtime[mu]>=u_time){                              q1[sizeqmcmd]=L_QQ_Q_l_upper->quat[0][1];
1374                           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){                              q2[sizeqmcmd]=L_QQ_Q_l_upper->quat[0][2];
1375                             Int_t sizeqmcmd = qtime.size();                              q3[sizeqmcmd]=L_QQ_Q_l_upper->quat[0][3];
1376                             inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);                              qmode[sizeqmcmd]=holeq(lowerqtime,qtime[sizeqmcmd],L_QQ_Q_l_lower,L_QQ_Q_l_upper,ui);
1377                             qtime[sizeqmcmd]=u_time;                              lowerqtime = u_time;
1378                             q0[sizeqmcmd]=L_QQ_Q_l_upper->quat[0][0];                              orbits.getPosition((double) (u_time - gltle->GetFromTime())/60., &eCi);
1379                             q1[sizeqmcmd]=L_QQ_Q_l_upper->quat[0][1];                              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]);
1380                             q2[sizeqmcmd]=L_QQ_Q_l_upper->quat[0][2];                              qRoll[sizeqmcmd]=RYPang_upper->Kren;
1381                             q3[sizeqmcmd]=L_QQ_Q_l_upper->quat[0][3];                              qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;
1382                             qmode[sizeqmcmd]=holeq(lowerqtime,qtime[sizeqmcmd],L_QQ_Q_l_lower,L_QQ_Q_l_upper,ui);                              qPitch[sizeqmcmd]=RYPang_upper->Tangazh;
1383                             lowerqtime = u_time;                              CopyQ(L_QQ_Q_l_lower,L_QQ_Q_l_upper);
1384                             orbits.getPosition((double) (u_time - gltle->GetFromTime())/60., &eCi);                              break;
1385                             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]);                            }
1386                             qRoll[sizeqmcmd]=RYPang_upper->Kren;                          }
1387                             qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;                        }
1388                             qPitch[sizeqmcmd]=RYPang_upper->Tangazh;                      }
1389                             CopyQ(L_QQ_Q_l_lower,L_QQ_Q_l_upper);                    }
1390                             break;                  }
1391                           }                }
                       }  
                     }  
                   }  
                 }  
               }  
             }  
             //if ( debug ) cout << "subpacket " << j3 << "\t qtime = " << qtime[qtime.size()-1] << endl;  
           }  
         }  
           
         if(qtime.size()==0){                            // in case if no orientation information in data  
           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;  
1392              }              }
1393            }            }
1394          }    
1395                      if(qtime.size()==0){        // in case if no orientation information in data
1396                if ( debug ) cout << "qtime.size() = 0" << endl;
1397          if ( debug ) printf(" puffi \n");              for(UInt_t my=0;my<recqtime.size();my++){
1398          Double_t tmin = 9999999999.;                if(sqrt(pow(recq0[my],2)+pow(recq1[my],2)+pow(recq2[my],2)+pow(recq3[my],2))>0.99999){
1399          Double_t tmax = 0.;                  Int_t sizeqmcmd = qtime.size();
1400          for(UInt_t tre = 0;tre<qtime.size();tre++){                  inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);
1401            if(qtime[tre]>tmax)tmax = qtime[tre];                  qtime[sizeqmcmd]=recqtime[my];
1402            if(qtime[tre]<tmin)tmin = qtime[tre];                  q0[sizeqmcmd]=recq0[my];
1403          }                  q1[sizeqmcmd]=recq1[my];
1404          // sorting quaternions by time                  q2[sizeqmcmd]=recq2[my];
1405         Bool_t t = true;                  q3[sizeqmcmd]=recq3[my];
1406          while(t){                  qmode[sizeqmcmd]=-10;
1407           t=false;                  orbits.getPosition((double) (qtime[sizeqmcmd] - gltle->GetFromTime())/60., &eCi);
1408            for(UInt_t i=0;i<qtime.size()-1;i++){                  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]);
1409              if(qtime[i]>qtime[i+1]){                  qRoll[sizeqmcmd]=RYPang_upper->Kren;
1410                Double_t tmpr = qtime[i];                  qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;
1411                qtime[i]=qtime[i+1];                  qPitch[sizeqmcmd]=RYPang_upper->Tangazh;
1412                qtime[i+1] = tmpr;                }
1413                tmpr = q0[i];              }
1414                q0[i]=q0[i+1];            }
1415                q0[i+1] = tmpr;            //if ( debug ) printf(" puffi \n");
1416                tmpr = q1[i];            Double_t tmin = 9999999999.;
1417                q1[i]=q1[i+1];            Double_t tmax = 0.;
1418                q1[i+1] = tmpr;            for(UInt_t tre = 0;tre<qtime.size();tre++){
1419                tmpr = q2[i];              if(qtime[tre]>tmax)tmax = qtime[tre];
1420                q2[i]=q2[i+1];              if(qtime[tre]<tmin)tmin = qtime[tre];
1421                q2[i+1] = tmpr;            }
1422                tmpr = q3[i];            // sorting quaternions by time
1423                q3[i]=q3[i+1];            Bool_t t = true;
1424                q3[i+1] = tmpr;            while(t){
1425                tmpr = qRoll[i];              t=false;
1426                qRoll[i]=qRoll[i+1];              for(UInt_t i=0;i<qtime.size()-1;i++){
1427                qRoll[i+1] = tmpr;                if(qtime[i]>qtime[i+1]){
1428                tmpr = qYaw[i];                  Double_t tmpr = qtime[i];
1429                qYaw[i]=qYaw[i+1];                  qtime[i]=qtime[i+1];
1430                qYaw[i+1] = tmpr;                  qtime[i+1] = tmpr;
1431                tmpr = qPitch[i];                  tmpr = q0[i];
1432                qPitch[i]=qPitch[i+1];                  q0[i]=q0[i+1];
1433                qPitch[i+1] = tmpr;                  q0[i+1] = tmpr;
1434                  t=true;                  tmpr = q1[i];
1435              }                  q1[i]=q1[i+1];
1436            }                  q1[i+1] = tmpr;
1437          }                  tmpr = q2[i];
1438                    q2[i]=q2[i+1];
1439          if ( debug ){                  q2[i+1] = tmpr;
1440            cout << "we have loaded quaternions: size of quaternions set is "<< qtime.size() << endl;                  tmpr = q3[i];
1441           for(UInt_t i=0;i<qtime.size();i++) cout << qtime[i] << "\t";                  q3[i]=q3[i+1];
1442            cout << endl << endl;                  q3[i+1] = tmpr;
1443            Int_t lopu;                  tmpr = qRoll[i];
1444            cin >> lopu;                  qRoll[i]=qRoll[i+1];
1445         }                  qRoll[i+1] = tmpr;
1446                    tmpr = qYaw[i];
1447                    qYaw[i]=qYaw[i+1];
1448                    qYaw[i+1] = tmpr;
1449                    tmpr = qPitch[i];
1450                    qPitch[i]=qPitch[i+1];
1451                    qPitch[i+1] = tmpr;
1452                    t=true;
1453                  }
1454                }
1455              }
1456              if ( debug ){
1457                cout << "we have loaded quaternions: size of quaternions set is "<< qtime.size() << endl;
1458                for(UInt_t i=0;i<qtime.size();i++) cout << qtime[i] << "\t";
1459                cout << endl << endl;
1460                Int_t lopu;
1461                cin >> lopu;
1462              }
1463        } // if we processed first event        } // if we processed first event
   
1464                
1465        //Filling Inclination information        //Filling Inclination information
1466        Double_t incli = 0;        Double_t incli = 0;
1467        if ( qtime.size() > 1 ){        if ( qtime.size() > 1 ){
1468          if ( debug ) cout << "ok quaternions is exist and mu = " << must << endl;          if(atime<qtime[0]){
1469          if ( debug ) cout << "qtimes[ " << qtime[0] << " , " << qtime[qtime.size()-1] << " ]\tatime = "<<atime<<endl;            for(UInt_t mu = 1;mu<qtime.size()-1;mu++){
1470                if(qtime[mu]>qtime[0]){
1471                  incli = (qPitch[mu]-qPitch[0])/(qtime[mu]-qtime[0]);
1472                  orbitalinfo->theta =  incli*atime+qPitch[mu]-incli*qtime[mu];
1473                  incli = (qRoll[mu]-qRoll[0])/(qtime[mu]-qtime[0]);
1474                  orbitalinfo->etha =  incli*atime+qRoll[mu]-incli*qtime[mu];
1475                  incli = (qYaw[mu]-qYaw[0])/(qtime[mu]-qtime[0]);
1476                  orbitalinfo->phi =  incli*atime+qYaw[mu]-incli*qtime[mu];
1477                  
1478                  incli = (q0[mu]-q0[0])/(qtime[mu]-qtime[0]);
1479                  orbitalinfo->q0 =  incli*atime+q0[mu]-incli*qtime[mu];
1480                  incli = (q1[mu]-q1[0])/(qtime[mu]-qtime[0]);
1481                  orbitalinfo->q1 =  incli*atime+q1[mu]-incli*qtime[mu];
1482                  incli = (q2[mu]-q2[0])/(qtime[mu]-qtime[0]);
1483                  orbitalinfo->q2 =  incli*atime+q2[mu]-incli*qtime[mu];
1484                  incli = (q3[mu]-q3[0])/(qtime[mu]-qtime[0]);
1485                  orbitalinfo->q3 =  incli*atime+q3[mu]-incli*qtime[mu];
1486                  orbitalinfo->TimeGap=qtime[0]-atime;
1487                  break;
1488                }
1489              }
1490            }
1491            Float_t eend=qtime.size()-1;
1492            if(atime>qtime[eend]){
1493              for(UInt_t mu=eend-1;mu>=0;mu--){
1494                if(qtime[mu]<qtime[eend]){
1495                  incli = (qPitch[eend]-qPitch[mu])/(qtime[eend]-qtime[mu]);
1496                  orbitalinfo->theta =  incli*atime+qPitch[eend]-incli*qtime[eend];
1497                  incli = (qRoll[eend]-qRoll[mu])/(qtime[eend]-qtime[mu]);
1498                  orbitalinfo->etha =  incli*atime+qRoll[eend]-incli*qtime[eend];
1499                  incli = (qYaw[eend]-qYaw[mu])/(qtime[eend]-qtime[mu]);
1500                  orbitalinfo->phi =  incli*atime+qYaw[eend]-incli*qtime[eend];
1501              
1502                  incli = (q0[eend]-q0[mu])/(qtime[eend]-qtime[mu]);
1503                  orbitalinfo->q0 =  incli*atime+q0[eend]-incli*qtime[eend];
1504                  incli = (q1[eend]-q1[mu])/(qtime[eend]-qtime[mu]);
1505                  orbitalinfo->q1 =  incli*atime+q1[eend]-incli*qtime[eend];
1506                  incli = (q2[eend]-q2[mu])/(qtime[eend]-qtime[mu]);
1507                  orbitalinfo->q2 =  incli*atime+q2[eend]-incli*qtime[eend];
1508                  incli = (q3[eend]-q3[mu])/(qtime[eend]-qtime[mu]);
1509                  orbitalinfo->q3 =  incli*atime+q3[eend]-incli*qtime[eend];
1510                  break;
1511                }
1512              }
1513            }
1514          for(UInt_t mu = must;mu<qtime.size()-1;mu++){          for(UInt_t mu = must;mu<qtime.size()-1;mu++){
1515            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);
1516            if(qtime[mu+1]>qtime[mu]){            if(qtime[mu+1]>qtime[mu]){
# Line 1502  if(recqtime[mu]>=1160987921.75 && recqti Line 1535  if(recqtime[mu]>=1160987921.75 && recqti
1535                orbitalinfo->q3 =  incli*atime+q3[mu+1]-incli*qtime[mu+1];                orbitalinfo->q3 =  incli*atime+q3[mu+1]-incli*qtime[mu+1];
1536                Float_t tg = (qtime[mu+1]-qtime[mu])/1000.0;                Float_t tg = (qtime[mu+1]-qtime[mu])/1000.0;
1537                if(tg>=1) tg=0.00;                if(tg>=1) tg=0.00;
1538                orbitalinfo->TimeGap = TMath::Min(TMath::Abs(qtime[mu+1])-atime,TMath::Abs(atime-qtime[mu]))+tg;//qtime[mu+1]-qtime[mu];                orbitalinfo->TimeGap = TMath::Min(TMath::Abs(qtime[mu+1]-atime),TMath::Abs(atime-qtime[mu]))+tg;//qtime[mu+1]-qtime[mu];
1539                orbitalinfo->mode = qmode[mu+1];                orbitalinfo->mode = qmode[mu+1];
1540                //if(atime==qtime[mu] || atime==qtime[mu+1]) orbitalinfo->qkind = 0; else orbitalinfo->qkind=1;                //if(atime==qtime[mu] || atime==qtime[mu+1]) orbitalinfo->qkind = 0; else orbitalinfo->qkind=1;
1541                //if(qmode[mu+1]==-10) orbitalinfo->R10r = true;else orbitalinfo->R10r = false;                //if(qmode[mu+1]==-10) orbitalinfo->R10r = true;else orbitalinfo->R10r = false;
# Line 1518  if(recqtime[mu]>=1160987921.75 && recqti Line 1551  if(recqtime[mu]>=1160987921.75 && recqti
1551        //        //
1552                
1553        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 ){
1554          if ( debug ) cout << "ops no iclination information" << endl;          if (debug) cout << "Warning: no iclination information "<< endl;
1555          orbitalinfo->mode = 10;            orbitalinfo->mode = 10;
1556          orbitalinfo->q0 = -1000.;            orbitalinfo->q0 = -1000.;
1557          orbitalinfo->q1 = -1000.;            orbitalinfo->q1 = -1000.;
1558          orbitalinfo->q2 = -1000.;            orbitalinfo->q2 = -1000.;
1559          orbitalinfo->q3 = -1000.;            orbitalinfo->q3 = -1000.;
1560          orbitalinfo->etha = -1000.;            orbitalinfo->etha = -1000.;
1561          orbitalinfo->phi = -1000.;            orbitalinfo->phi = -1000.;
1562          orbitalinfo->theta = -1000.;            orbitalinfo->theta = -1000.;
1563          orbitalinfo->TimeGap = -1000.;            orbitalinfo->TimeGap = -1000.;
1564          //orbitalinfo->qkind = -1000;            TMatrixD Iij(3,3);
1565              Iij(0,0)=0; Iij(0,1)=0; Iij(0,2)=0;
1566              Iij(1,0)=0; Iij(1,1)=0; Iij(1,2)=0;
1567              Iij(2,0)=0; Iij(2,1)=0; Iij(2,2)=0;
1568              Iij.Zero();
1569              orbitalinfo->Iij.ResizeTo(Iij);
1570              orbitalinfo->Iij = Iij;
1571              //orbitalinfo->qkind = -1000;
1572                    
1573          //      if ( debug ){            //  if ( debug ){
1574          //        Int_t lopu;            //    Int_t lopu;
1575          //         cin >> lopu;            //         cin >> lopu;
1576          //      }            //  }
1577          if ( debug ) printf(" grfuffi6 \n");            if ( debug ) printf(" grfuffi6 \n");
       }  
       //  
       if ( debug ) printf(" filling \n");  
       // #########################################################################################################################    
       //  
       // fill orbital positions  
       //          
       // Build coordinates in the right range.  We want to convert,  
       // longitude from (0, 2*pi) to (-180deg, 180deg).  Altitude is  
       // in meters.  
       lon = (coo.m_Lon > M_PI) ? rad2deg(coo.m_Lon - 2*M_PI) : rad2deg(coo.m_Lon);  
       lat = rad2deg(coo.m_Lat);  
       alt = coo.m_Alt;  
   
       cOrbit orbits2(*gltle->GetTle());  
       orbits2.getPosition((double) (atime - gltle->GetFromTime())/60., &eCi);  
       //      Float_t x=eCi.getPos().m_x;  
       //      Float_t y=eCi.getPos().m_y;  
       //      Float_t z=eCi.getPos().m_z;  
         
       TVector3 V(eCi.getVel().m_x,eCi.getVel().m_y,eCi.getVel().m_z);  
       TVector3 Pos(eCi.getPos().m_x,eCi.getPos().m_y,eCi.getPos().m_z);  
         
       Float_t dlon=Pos.Phi()*TMath::RadToDeg()-lon;  
         
       Pos.RotateZ(-dlon*TMath::DegToRad());  
       V.RotateZ(-dlon*TMath::DegToRad());  
       Float_t diro;  
       if(V.Z()>0) diro=1; else diro=-1;  
         
       // 10REDNEW  
       Int_t errq=0;  
       Int_t azim=0;  
       Int_t qual=0;  
       Int_t MU=0;  
       for(UInt_t mu = 0;mu<RTtime2.size()-1;mu++){  
         if(atime<RTstart[mu+1] && atime>=RTstart[mu]){  
           errq=RTerrq[mu];  
           azim=RTazim[mu];  
           qual=RTqual[mu];  
           MU=mu;  
           break;  
1578          }          }
1579        }          //
1580        orbitalinfo->errq = errq;          if ( debug ) printf(" filling \n");
1581        orbitalinfo->azim = azim;          // #########################################################################################################################  
1582        orbitalinfo->rtqual=qual;          //
1583        orbitalinfo->qkind = 0;          // fill orbital positions
1584                  //        
1585        if ( debug ) printf(" coord done \n");          // Build coordinates in the right range.  We want to convert,
1586        if( lon<180 && lon>-180 && lat<90 && lat>-90 && alt>0 ){            // longitude from (0, 2*pi) to (-180deg, 180deg).  Altitude is
1587          //                // in meters.
1588          orbitalinfo->lon = lon;          lon = (coo.m_Lon > M_PI) ? rad2deg(coo.m_Lon - 2*M_PI) : rad2deg(coo.m_Lon);
1589          orbitalinfo->lat = lat;          lat = rad2deg(coo.m_Lat);
1590          orbitalinfo->alt = alt;          alt = coo.m_Alt;
1591          orbitalinfo->V = V;  
1592            cOrbit orbits2(*gltle->GetTle());
1593          //      GMtype_CoordGeodetic  location;          orbits2.getPosition((double) (atime - gltle->GetFromTime())/60., &eCi);
1594          location.lambda = lon;          //      Float_t x=eCi.getPos().m_x;
1595          location.phi = lat;          //      Float_t y=eCi.getPos().m_y;
1596          location.HeightAboveEllipsoid = alt;          //      Float_t z=eCi.getPos().m_z;
1597            
1598          GM_GeodeticToSpherical(Ellip, location, &CoordSpherical);          TVector3 V(eCi.getVel().m_x,eCi.getVel().m_y,eCi.getVel().m_z);
1599          GM_SphericalToCartesian(CoordSpherical,  &CoordCartesian);          TVector3 Pos(eCi.getPos().m_x,eCi.getPos().m_y,eCi.getPos().m_z);
1600          GM_EarthCartToDipoleCartCD(Pole, CoordCartesian, &DipoleCartesian);          
1601          GM_CartesianToSpherical(DipoleCartesian, &DipoleSpherical);          Float_t dlon=Pos.Phi()*TMath::RadToDeg()-lon;
1602          orbitalinfo->londip = DipoleSpherical.lambda;          
1603          orbitalinfo->latdip = DipoleSpherical.phig;          Pos.RotateZ(-dlon*TMath::DegToRad());
1604            V.RotateZ(-dlon*TMath::DegToRad());
1605          if(debug)cout<<"geodetic:\t"<<lon<<"\t"<<lat<<"\tgeomagnetic:\t"<<orbitalinfo->londip<<"\t"<<orbitalinfo->latdip<<endl;          Float_t diro;
1606            if(V.Z()>0) diro=1; else diro=-1;
1607          //          
1608          // compute mag field components and L shell.          // 10REDNEW
1609          //          Int_t errq=0;
1610          if ( debug ) printf(" call igrf feldg \n");          Int_t azim=0;
1611          feldg_(&lat, &lon, &alt, &bnorth, &beast, &bdown, &babs);          Int_t qual=0;
1612          if ( debug ) printf(" call igrf shellg \n");          Int_t MU=0;
1613          shellg_(&lat, &lon, &alt, &dimo, &xl, &icode, &bab1);          for(UInt_t mu = 0;mu<RTtime2.size()-1;mu++){
1614          if ( debug ) printf(" call igrf findb \n");            if(atime<RTstart[mu+1] && atime>=RTstart[mu]){
1615          findb0_(&stps, &bdel, &value, &bequ, &rr0);              errq=RTerrq[mu];
1616          //              azim=RTazim[mu];
1617          if ( debug ) printf(" done igrf \n");              qual=RTqual[mu];
1618          orbitalinfo->Bnorth = bnorth;              MU=mu;
1619          orbitalinfo->Beast = beast;              break;
1620          orbitalinfo->Bdown = bdown;            }
1621          orbitalinfo->Babs = babs;          }
1622          orbitalinfo->M = dimo;          orbitalinfo->errq = errq;
1623          orbitalinfo->BB0 = babs/bequ;          orbitalinfo->azim = azim;
1624          orbitalinfo->L = xl;                orbitalinfo->rtqual=qual;
1625          // Set Stormer vertical cutoff using L shell.          orbitalinfo->qkind = 0;
1626          orbitalinfo->cutoffsvl = 14.295 / (xl*xl); //          
1627          if(debug)cout << "L = " << xl << "\tM = " << dimo << "\tvertical cutoff:  "<< orbitalinfo->cutoffsvl << endl;          if ( debug ) printf(" coord done \n");
1628            if( lon<180 && lon>-180 && lat<90 && lat>-90 && alt>0 ){  
1629                      orbitalinfo->lon = lon;
1630  //           ---------- Forwarded message ----------            orbitalinfo->lat = lat;
1631  //           Date: Wed, 09 May 2012 12:16:47 +0200            orbitalinfo->alt = alt;
1632  //           From: Alessandro Bruno <alessandro.bruno@ba.infn.it>            orbitalinfo->V = V;
1633  //           To: Mirko Boezio <mirko.boezio@ts.infn.it>  
1634  //           Cc: Francesco S. Cafagna <Francesco.Cafagna@ba.infn.it>            //  GMtype_CoordGeodetic  location;
1635  //           Subject: Störmer vertical cutoff            location.lambda = lon;
1636              location.phi = lat;
1637  //           Ciao Mirko,            location.HeightAboveEllipsoid = alt;
1638  //           volevo segnalarti che il valore dello Störmer vertical cutoff nel Level2 è  
1639  //           sovrastimato di circa il 4%.            GM_GeodeticToSpherical(Ellip, location, &CoordSpherical);
1640  //           Dopo un'approfondita analisi con l'IGRF-05 abbiamo ricavano un valore pari            GM_SphericalToCartesian(CoordSpherical,  &CoordCartesian);
1641  //           a: 14.295 / L^2 anzichè 14.9 / L^2, valore obsoleto in quanto riferito agli            GM_EarthCartToDipoleCartCD(Pole, CoordCartesian, &DipoleCartesian);
1642  //           anni '50.            GM_CartesianToSpherical(DipoleCartesian, &DipoleSpherical);
1643  //                    orbitalinfo->londip = DipoleSpherical.lambda;
1644          //14.9/(xl*xl);            orbitalinfo->latdip = DipoleSpherical.phig;
1645          orbitalinfo->igrf_icode = icode;  
1646          //            if(debug)cout<<"geodetic:\t"<<lon<<"\t"<<lat<<"\tgeomagnetic:\t"<<orbitalinfo->londip<<"\t"<<orbitalinfo->latdip<<endl;
1647        }        
1648        //            //
1649        if ( debug ) printf(" pitch angle \n");            // compute mag field components and L shell.
1650        //            //
1651        // pitch angles            if ( debug ) printf(" call igrf feldg \n");
1652        //            feldg_(&lat, &lon, &alt, &bnorth, &beast, &bdown, &babs);
1653        if( orbitalinfo->TimeGap>0){            if ( debug ) printf(" call igrf shellg \n");
1654          //            shellg_(&lat, &lon, &alt, &dimo, &xl, &icode, &bab1);
1655          if ( debug ) printf(" timegap %f \n",orbitalinfo->TimeGap);            if ( debug ) printf(" call igrf findb \n");
1656          Float_t Bx = -orbitalinfo->Bdown;            findb0_(&stps, &bdel, &value, &bequ, &rr0);
1657          Float_t By = orbitalinfo->Beast;            //
1658          Float_t Bz = orbitalinfo->Bnorth;            if ( debug ) printf(" done igrf \n");
1659                orbitalinfo->Bnorth = bnorth;
1660  //      TMatrixD Qiji(3,3);              orbitalinfo->Beast = beast;
1661          TMatrixD Qij = PO->QuatoECI(orbitalinfo->q0,orbitalinfo->q1,orbitalinfo->q2,orbitalinfo->q3);              orbitalinfo->Bdown = bdown;
1662          TMatrixD Dij = PO->ECItoGEO(Qij,orbitalinfo->absTime,orbitalinfo->lat,orbitalinfo->lon);              orbitalinfo->Babs = babs;
1663                orbitalinfo->M = dimo;
1664  //10REDNEW              orbitalinfo->BB0 = babs/bequ;
1665          // If initial orientation data have reason to be inaccurate              orbitalinfo->L = xl;      
1666         Float_t tg = 0.00;              // Set Stormer vertical cutoff using L shell.
1667         Float_t tmptg;              orbitalinfo->cutoffsvl = 14.295 / (xl*xl); //
1668         if(MU!=0){              if(debug)cout << "L = " << xl << "\tM = " << dimo << "\tvertical cutoff:  "<< orbitalinfo->cutoffsvl << endl;
1669  //      if(orbitalinfo->TimeGap>0 && errq==0 && azim==0){               // 10RED CHECK  (comparison between three metod of recovering orientation)                      
1670         if((atime>=RTstart[MU] && atime<RTstart[MU+1] && RTbank1[MU]==0 && RTbank2[MU]==0 && TMath::Abs(orbitalinfo->etha)>0.1) || ((RTbank1[MU]!=0 || RTbank2[MU]!=0) && atime>=RTstart[MU] && atime<RTstart[MU+1] && azim==0 && (errq!=0 || orbitalinfo->TimeGap>10.0 || ((modf(orbitalinfo->TimeGap,&tmptg)*1000>10 || modf(orbitalinfo->TimeGap,&tmptg)*1000==0.0) && orbitalinfo->TimeGap>2.0)))){              //           ---------- Forwarded message ----------
1671          //found in Rotation Table this data for this time interval              //           Date: Wed, 09 May 2012 12:16:47 +0200
1672          if(atime<RTtime1[0])              //           From: Alessandro Bruno <alessandro.bruno@ba.infn.it>
1673            orbitalinfo->azim = 5;        //means that RotationTable no started yet              //           To: Mirko Boezio <mirko.boezio@ts.infn.it>
1674         else{              //           Cc: Francesco S. Cafagna <Francesco.Cafagna@ba.infn.it>
1675                  // search for angle betwean velosity and direction to north in tangential to Earth surfase plane in satellite position              //           Subject: Störmer vertical cutoff
1676                Double_t bank=RTstart[MU];  
1677                Double_t tlat=orbitalinfo->lat;              //           Ciao Mirko,
1678                //           volevo segnalarti che il valore dello Störmer vertical cutoff nel Level2 è
1679                tg=modf(orbitalinfo->TimeGap,&tg)*1000;              //           sovrastimato di circa il 4%.
1680                //           Dopo un'approfondita analisi con l'IGRF-05 abbiamo ricavano un valore pari
1681                if(atime>=RTpluto1[MU] && atime<=RTpluto2[MU]){              //           a: 14.295 / L^2 anzichè 14.9 / L^2, valore obsoleto in quanto riferito agli
1682                  Double_t kar=(RTbank2[MU]-RTbank1[MU])/(RTtime2[MU]-RTtime1[MU]);              //           anni '50.
1683                  Double_t bak=RTbank1[MU]-kar*RTtime1[MU];              //        
1684                  bank=kar*atime+bak;              //           14.9/(xl*xl);
1685                }              orbitalinfo->igrf_icode = (Float_t)icode;
1686                if(atime>=RTstart[MU] && atime<RTpluto1[MU]){              //
1687                   Double_t s_dBdt2=(RTbpluto1[MU]-RTbank1[MU])/(Int_t)(RTpluto1[MU]-RTstart[MU]);            }  //check lon lat alt      
1688                   Double_t s_t2=((Double_t)RTpluto1[MU]+(Double_t)RTstart[MU])/2. - RTstart[MU];            //
1689                   Double_t s_t1=RTstart[MU]-RTstart[MU];            if ( debug ) printf(" pitch angle \n");
1690                   Double_t s_k=s_dBdt2/(s_t2-s_t1);            //
1691                   Double_t s_b=-s_k*s_t1;            // pitch angles
1692                   Double_t s_t3=RTpluto1[MU]-RTstart[MU];            //
1693                   Double_t s_b3=RTbpluto1[MU];            if( orbitalinfo->TimeGap>=0){
1694                   Double_t s_c=s_b3-0.5*s_k*s_t3*s_t3 -s_b*s_t3;              //
1695                   bank=0.5*s_k*(atime-RTstart[MU])*(atime-RTstart[MU]) + s_b*(atime-RTstart[MU]) + s_c;              if ( debug ) printf(" timegap %f \n",orbitalinfo->TimeGap);
1696               }              Float_t Bx = -orbitalinfo->Bdown;
1697                if(atime>RTpluto2[MU] && atime<=RTstart[MU+1]){              Float_t By = orbitalinfo->Beast;
1698                   Double_t s_dBdt2=(RTbpluto2[MU] - RTbank2[MU])/(Int_t)(RTpluto2[MU]-RTstart[MU+1]);              Float_t Bz = orbitalinfo->Bnorth;
1699                   Double_t s_t2=((Double_t)RTpluto2[MU]+(Double_t)RTstart[MU+1])/2. - RTstart[MU];  
1700                   Double_t s_t1=RTstart[MU+1]-RTstart[MU];            //  TMatrixD Qiji(3,3);
1701                   Double_t s_k=s_dBdt2/(s_t2-s_t1);              TMatrixD Qij = PO->QuatoECI(orbitalinfo->q0,orbitalinfo->q1,orbitalinfo->q2,orbitalinfo->q3);
1702                   Double_t s_b=-s_k*s_t1;              TMatrixD Dij = PO->ECItoGEO(Qij,orbitalinfo->absTime,orbitalinfo->lat,orbitalinfo->lon);
1703                   Double_t s_t3=RTpluto2[MU]-RTstart[MU];  
1704                   Double_t s_b3=RTbpluto2[MU];              //10REDNEW
1705                   Double_t s_c=s_b3-0.5*s_k*s_t3*s_t3 -s_b*s_t3;              // If initial orientation data have reason to be inaccurate
1706                 bank=0.5*s_k*(atime-RTstart[MU])*(atime-RTstart[MU]) + s_b*(atime-RTstart[MU]) + s_c;              Float_t tg = 0.00;
1707               }              Float_t tmptg;
1708                if(TMath::Abs(orbitalinfo->etha-bank)>0.1){              Bool_t tgpar=false;
1709                  orbitalinfo->etha=bank;              Bool_t tgpar0=false;
1710                  Double_t spitch = 0.00001;  // temprary not zero to avoid problem with tranzition from Euler angles to orientation matrix              if (orbitalinfo->TimeGap>10.0 || ((modf(orbitalinfo->TimeGap,&tmptg)*1000>10 || modf(orbitalinfo->TimeGap,&tmptg)*1000==0.0) && orbitalinfo->TimeGap>2.0)) tgpar=true;
1711                if (orbitalinfo->TimeGap>180.0) tgpar0=true;
1712                  //Estimations of pitch angle of satellite              if(MU!=0){
1713                  if(TMath::Abs(bank)>0.7){              //  if(orbitalinfo->TimeGap>0 && errq==0 && azim==0){    // 10RED CHECK  (comparison between three metod of recovering orientation)
1714                     Float_t spitch1=TMath::DegToRad()*0.7*diro;//RTdir1[MU];                if((atime>=RTstart[MU] && atime<RTstart[MU+1] && RTbank1[MU]==0 && RTbank2[MU]==0 && (TMath::Abs(orbitalinfo->etha)>0.1 || tgpar0)) || ((RTbank1[MU]!=0 || RTbank2[MU]!=0) && atime>=RTstart[MU] && atime<RTstart[MU+1] && azim==0 && (errq!=0 || tgpar))){
1715                     Float_t spitch2=TMath::DegToRad()*0.7*diro;//RTdir2[MU];                  //found in Rotation Table this data for this time interval
1716                     Float_t kva=(spitch2-spitch1)/(RTtime2[MU]-RTtime1[MU]);                  if(atime<RTtime1[0])
1717                     Float_t bva=spitch1-kva*RTtime1[MU];                    orbitalinfo->azim = 5;  //means that RotationTable no started yet
1718                     spitch=kva*atime+bva;                  else{
1719                  }                    // search for angle betwean velosity and direction to north in tangential to Earth surfase plane in satellite position
1720                      Double_t bank=RTstart[MU];
1721                  //Calculate Yaw angle accordingly with fit, see picture FitYaw.jpg                    Double_t tlat=orbitalinfo->lat;
1722                  Double_t yaw=0.00001;  // temprary not zero to avoid problem with tranzition from Euler angles to orientation matrix  
1723                  if(TMath::Abs(tlat)<70)                    tg=modf(orbitalinfo->TimeGap,&tg)*1000;
1724                    yaw = -3.7e-8*tlat*tlat*tlat*tlat + 1.4e-7*tlat*tlat*tlat - 0.0005*tlat*tlat - 0.00025*tlat + 3.6;  
1725                  yaw = diro*yaw; //because should be different sign for ascending and descending orbits!                    if(atime>=RTpluto1[MU] && atime<=RTpluto2[MU]){
1726                  orbitalinfo->phi=yaw;                      Double_t kar=(RTbank2[MU]-RTbank1[MU])/(RTtime2[MU]-RTtime1[MU]);
1727                        Double_t bak=RTbank1[MU]-kar*RTtime1[MU];
1728                  if(TMath::Abs(bank)>0.5 && TMath::Abs(yaw-orbitalinfo->phi)<3.0) yaw=orbitalinfo->phi;                      bank=kar*atime+bak;
1729                      }
1730  //              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                    if(atime>=RTstart[MU] && atime<RTpluto1[MU]){
1731                  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                      Double_t s_dBdt2=(RTbpluto1[MU]-RTbank1[MU])/(Int_t)(RTpluto1[MU]-RTstart[MU]);
1732                  orbitalinfo->qkind = 1;                      Double_t s_t2=((Double_t)RTpluto1[MU]+(Double_t)RTstart[MU])/2. - RTstart[MU];
1733               }                      Double_t s_t1=RTstart[MU]-RTstart[MU];
1734                        Double_t s_k=s_dBdt2/(s_t2-s_t1);
1735            //Qij = PO->GEOtoECI(Dij,orbitalinfo->absTime,orbitalinfo->lat,orbitalinfo->lon); // to convert from Dij to Qij                      Double_t s_b=-s_k*s_t1;
1736                        Double_t s_t3=RTpluto1[MU]-RTstart[MU];
1737          } // end of if(atime<RTtime1[0]                      Double_t s_b3=RTbpluto1[MU];
1738          } // end of (((orbitalinfo->TimeGap>60.0 && TMath...                      Double_t s_c=s_b3-0.5*s_k*s_t3*s_t3 -s_b*s_t3;
1739         } // end of MU~=0                      bank=0.5*s_k*(atime-RTstart[MU])*(atime-RTstart[MU]) + s_b*(atime-RTstart[MU]) + s_c;
1740                      }
1741          TMatrixD qij = PO->ColPermutation(Qij);                    if(atime>RTpluto2[MU] && atime<=RTstart[MU+1]){
1742          TMatrixD Fij = PO->ECItoGreenwich(Qij,orbitalinfo->absTime);                      Double_t s_dBdt2=(RTbpluto2[MU] - RTbank2[MU])/(Int_t)(RTpluto2[MU]-RTstart[MU+1]);
1743          TMatrixD Gij = PO->ColPermutation(Fij);                      Double_t s_t2=((Double_t)RTpluto2[MU]+(Double_t)RTstart[MU+1])/2. - RTstart[MU];
1744          Dij = PO->ECItoGEO(Qij,orbitalinfo->absTime,orbitalinfo->lat,orbitalinfo->lon);                      Double_t s_t1=RTstart[MU+1]-RTstart[MU];
1745          TMatrixD Iij = PO->ColPermutation(Dij);                      Double_t s_k=s_dBdt2/(s_t2-s_t1);
1746          TVector3 SP = PO->GetSunPosition(orbitalinfo->absTime);                      Double_t s_b=-s_k*s_t1;
1747          // go to Pamela reference frame from Resurs reference frame                      Double_t s_t3=RTpluto2[MU]-RTstart[MU];
1748          Float_t tmpy = SP.Y();                      Double_t s_b3=RTbpluto2[MU];
1749          SP.SetY(SP.Z());                      Double_t s_c=s_b3-0.5*s_k*s_t3*s_t3 -s_b*s_t3;
1750          SP.SetZ(-tmpy);                      bank=0.5*s_k*(atime-RTstart[MU])*(atime-RTstart[MU]) + s_b*(atime-RTstart[MU]) + s_c;
1751          TVector3 SunZenith;                    }
1752          SunZenith.SetMagThetaPhi(1,23.439281*TMath::DegToRad(),TMath::Pi()/2.);                    if(TMath::Abs(orbitalinfo->etha-bank)>0.1){
1753          TVector3 SunMag = -SP;                      orbitalinfo->etha=bank;
1754          SunMag.Rotate(-45*TMath::DegToRad(),SunZenith);                      Double_t spitch = 0.00001;  // temprary not zero to avoid problem with tranzition from Euler angles to orientation matrix
1755          tmpy=SunMag.Y();  
1756          SunMag.SetY(SunMag.Z());                      //Estimations of pitch angle of satellite
1757          SunMag.SetZ(-tmpy);                      if(TMath::Abs(bank)>0.7){
1758                          Float_t spitch1=TMath::DegToRad()*0.7*diro;//RTdir1[MU];
1759          orbitalinfo->Iij.ResizeTo(Iij);                        Float_t spitch2=TMath::DegToRad()*0.7*diro;//RTdir2[MU];
1760          orbitalinfo->Iij = Iij;                        Float_t kva=(spitch2-spitch1)/(RTtime2[MU]-RTtime1[MU]);
1761          //                        Float_t bva=spitch1-kva*RTtime1[MU];
1762          //      A1 = Iij(0,2);                        spitch=kva*atime+bva;
1763          //      A2 = Iij(1,2);                      }
1764          //      A3 = Iij(2,2);  
1765          //                      //Calculate Yaw angle accordingly with fit, see picture FitYaw.jpg
1766          //      orbitalinfo->pamzenitangle = (Float_t)PO->GetPitchAngle(1,0,0,A1,A2,A3);                        // Angle between zenit and Pamela's main axiz                      Double_t yaw=0.00001;  // temprary not zero to avoid problem with tranzition from Euler angles to orientation matrix
1767          //      orbitalinfo->pamBangle = (Float_t)PO->GetPitchAngle(A1,A2,A3,Bx,By,Bz);                 // Angle between Pamela's main axiz and B                      if(TMath::Abs(tlat)<70)
1768          //                        yaw = -3.7e-8*tlat*tlat*tlat*tlat + 1.4e-7*tlat*tlat*tlat - 0.0005*tlat*tlat - 0.00025*tlat + 3.6;
1769                        yaw = diro*yaw;  //because should be different sign for ascending and descending orbits!
1770                        orbitalinfo->phi=yaw;
1771    
1772                        if(TMath::Abs(bank)>0.5 && TMath::Abs(yaw-orbitalinfo->phi)<3.0) yaw=orbitalinfo->phi;
1773    
1774                        //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
1775                        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
1776                        orbitalinfo->qkind = 1;
1777                      }
1778                      //Qij = PO->GEOtoECI(Dij,orbitalinfo->absTime,orbitalinfo->lat,orbitalinfo->lon); // to convert from Dij to Qij
1779                    } // end of if(atime<RTtime1[0]
1780                  } // end of (((orbitalinfo->TimeGap>60.0 && TMath...
1781                } // end of MU~=0
1782    
1783                TMatrixD qij = PO->ColPermutation(Qij);
1784                TMatrixD Fij = PO->ECItoGreenwich(Qij,orbitalinfo->absTime);
1785                TMatrixD Gij = PO->ColPermutation(Fij);
1786                Dij = PO->ECItoGEO(Qij,orbitalinfo->absTime,orbitalinfo->lat,orbitalinfo->lon);
1787                TMatrixD Iij = PO->ColPermutation(Dij);
1788    
1789                TVector3 SP = PO->GetSunPosition(orbitalinfo->absTime);
1790                // go to Pamela reference frame from Resurs reference frame
1791                Float_t tmpy = SP.Y();
1792                SP.SetY(SP.Z());
1793                SP.SetZ(-tmpy);
1794                TVector3 SunZenith;
1795                SunZenith.SetMagThetaPhi(1,23.439281*TMath::DegToRad(),TMath::Pi()/2.);
1796                TVector3 SunMag = -SP;
1797                SunMag.Rotate(-45*TMath::DegToRad(),SunZenith);
1798                tmpy=SunMag.Y();
1799                SunMag.SetY(SunMag.Z());
1800                SunMag.SetZ(-tmpy);
1801    
1802                orbitalinfo->Iij.ResizeTo(Iij);
1803                orbitalinfo->Iij = Iij;
1804    
1805                Bool_t saso=true;
1806                if (orbitalinfo->qkind==1) saso=true;
1807                if (orbitalinfo->qkind==0 && orbitalinfo->azim>0 && orbitalinfo->azim!=5 && tgpar) saso=false;
1808                if (orbitalinfo->qkind==0 && orbitalinfo->azim==5 && TMath::Abs(orbitalinfo->etha)>0.1 && tgpar) saso=false;
1809                if (orbitalinfo->qkind==0 && orbitalinfo->azim==5 && TMath::Abs(orbitalinfo->etha)<=0.1 && tgpar0) saso=false;
1810                if (saso) orbitalinfo->mode=orbitalinfo->rtqual; else orbitalinfo->mode=2;
1811    
1812      //
1813            //  A1 = Iij(0,2);
1814            //  A2 = Iij(1,2);
1815            //  A3 = Iij(2,2);
1816      //
1817      //  orbitalinfo->pamzenitangle = (Float_t)PO->GetPitchAngle(1,0,0,A1,A2,A3);            // Angle between zenit and Pamela's main axiz
1818      //  orbitalinfo->pamBangle = (Float_t)PO->GetPitchAngle(A1,A2,A3,Bx,By,Bz);      // Angle between Pamela's main axiz and B
1819      //
1820          if ( debug ) printf(" matrixes done  \n");          if ( debug ) printf(" matrixes done  \n");
1821          if ( !standalone ){    if ( !standalone ){
1822            if ( debug ) printf(" !standalone \n");            if ( debug ) printf(" !standalone \n");
1823            //      //
1824            // Standard tracking algorithm            // Standard tracking algorithm
1825            //            //
1826            Int_t nn = 0;      Int_t nn = 0;
1827            if ( verbose ) printf(" standard tracking \n");            if ( verbose ) printf(" standard tracking \n");
1828            for(Int_t nt=0; nt < tof->ntrk(); nt++){        for(Int_t nt=0; nt < tof->ntrk(); nt++){  
1829              //        //
1830              ToFTrkVar *ptt = tof->GetToFTrkVar(nt);        ToFTrkVar *ptt = tof->GetToFTrkVar(nt);
1831              if (debug) cout<<"tof->ntrk() = "<<tof->ntrk()<<"\tptt->trkseqno = "<<ptt->trkseqno<<"\ttrke->ntrk() = "<<trke->ntrk()<<endl;        if (debug) cout<<"tof->ntrk() = "<<tof->ntrk()<<"\tptt->trkseqno = "<<ptt->trkseqno<<"\ttrke->ntrk() = "<<trke->ntrk()<<endl;
1832              Double_t E11x = ptt->xtr_tof[0]; // tr->x[0];        Double_t E11x = ptt->xtr_tof[0]; // tr->x[0];
1833              Double_t E11y = ptt->ytr_tof[0]; //tr->y[0];        Double_t E11y = ptt->ytr_tof[0]; //tr->y[0];
1834              Double_t E11z = zin[0];        Double_t E11z = zin[0];
1835              Double_t E22x = ptt->xtr_tof[3];//tr->x[3];        Double_t E22x = ptt->xtr_tof[3];//tr->x[3];
1836              Double_t E22y = ptt->ytr_tof[3];//tr->y[3];        Double_t E22y = ptt->ytr_tof[3];//tr->y[3];
1837              Double_t E22z = zin[3];        Double_t E22z = zin[3];
1838              if ( (E11x < 100. && E11y < 100. && E22x < 100. && E22y < 100.) || ptt->trkseqno != -1  ){        if ( (E11x < 100. && E11y < 100. && E22x < 100. && E22y < 100.) || ptt->trkseqno != -1  ){
1839                TrkTrack *mytrack = trke->GetStoredTrack(ptt->trkseqno);          TrkTrack *mytrack = trke->GetStoredTrack(ptt->trkseqno);
1840                Float_t rig=1/mytrack->GetDeflection();          Float_t rig=1/mytrack->GetDeflection();
1841                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));
1842                //                //
1843                Px = (E22x-E11x)/norm;          Px = (E22x-E11x)/norm;
1844                Py = (E22y-E11y)/norm;          Py = (E22y-E11y)/norm;
1845                Pz = (E22z-E11z)/norm;          Pz = (E22z-E11z)/norm;
1846                //          //
1847                t_orb->trkseqno = ptt->trkseqno;          t_orb->trkseqno = ptt->trkseqno;
1848                //          //
1849                TMatrixD Eij = PO->PamelatoGEO(Iij,Px,Py,Pz);          TMatrixD Eij = PO->PamelatoGEO(Iij,Px,Py,Pz);
1850                t_orb->Eij.ResizeTo(Eij);          t_orb->Eij.ResizeTo(Eij);  
1851                t_orb->Eij = Eij;          t_orb->Eij = Eij;  
1852                //          //
1853                TMatrixD Sij = PO->PamelatoGEO(Gij,Px,Py,Pz);          TMatrixD Sij = PO->PamelatoGEO(Gij,Px,Py,Pz);
1854                t_orb->Sij.ResizeTo(Sij);          t_orb->Sij.ResizeTo(Sij);
1855                t_orb->Sij = Sij;          t_orb->Sij = Sij;
1856                //                      //      
1857                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);
1858                //          //
1859                // SunPosition in instrumental reference frame                // SunPosition in instrumental reference frame
1860                TMatrixD Kij = PO->PamelatoGEO(qij,Px,Py,Pz);          TMatrixD Kij = PO->PamelatoGEO(qij,Px,Py,Pz);
1861                TMatrixD Lij = PO->PamelatoGEO(qij,0,0,1);          TMatrixD Lij = PO->PamelatoGEO(qij,0,0,1);
1862                t_orb->sunangle=(Float_t)PO->GetPitchAngle(Kij(0,0),Kij(1,0),Kij(2,0),-SP.X(),-SP.Y(),-SP.Z());          t_orb->sunangle=(Float_t)PO->GetPitchAngle(Kij(0,0),Kij(1,0),Kij(2,0),-SP.X(),-SP.Y(),-SP.Z());
1863                t_orb->sunmagangle=(Float_t)PO->GetPitchAngle(Kij(0,0),Kij(1,0),Kij(2,0),SunMag.X(),SunMag.Y(),SunMag.Z());          t_orb->sunmagangle=(Float_t)PO->GetPitchAngle(Kij(0,0),Kij(1,0),Kij(2,0),SunMag.X(),SunMag.Y(),SunMag.Z());
1864                //          //
1865                //          //
1866                Double_t omega = PO->GetPitchAngle(-Eij(0,0),-Eij(1,0),-Eij(2,0),1,0,0) * TMath::DegToRad();          Double_t omega = PO->GetPitchAngle(-Eij(0,0),-Eij(1,0),-Eij(2,0),1,0,0) * TMath::DegToRad();
1867                TVector3 Bxy(0,By,Bz);          TVector3 Bxy(0,By,Bz);
1868                TVector3 Exy(0,-Eij(1,0),-Eij(2,0));          TVector3 Exy(0,-Eij(1,0),-Eij(2,0));
1869                Double_t dzeta=Bxy.Angle(Exy);          Double_t dzeta=Bxy.Angle(Exy);
1870                if (-Eij(1,0) < 0) dzeta=2.0*TMath::Pi() - dzeta;          if (-Eij(1,0) < 0) dzeta=2.0*TMath::Pi() - dzeta;
1871                                
1872                if(debug) cout << "omega = "<<omega*TMath::RadToDeg()<<"\tdzeta = "<<dzeta*TMath::RadToDeg()<<endl;                if(debug) cout << "omega = "<<omega*TMath::RadToDeg()<<"\tdzeta = "<<dzeta*TMath::RadToDeg()<<endl;
1873    
1874                // Formula from D.F. Smart *, M.A. Shea [2005]; A review of geomagnetic cutoff rigidities for earth-orbiting spacecraft                // Formula from D.F. Smart *, M.A. Shea [2005]; A review of geomagnetic cutoff rigidities for earth-orbiting spacecraft
1875                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));          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));
1876                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));          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));
1877                if (debug) cout << "R = " << rig << "\tcutoff = " << t_orb->cutoff << endl;          if (debug) cout << "R = " << rig << "\tcutoff = " << t_orb->cutoff << endl;
1878    
1879                //          //
1880                if ( t_orb->pitch != t_orb->pitch ) t_orb->pitch = -1000.;          if ( t_orb->pitch != t_orb->pitch ) t_orb->pitch = -1000.;
1881                if ( t_orb->cutoff != t_orb->cutoff ) t_orb->cutoff = -1000.;          if ( t_orb->cutoff != t_orb->cutoff ) t_orb->cutoff = -1000.;
1882                if ( t_orb->sunangle != t_orb->sunangle ) t_orb->sunangle = -1000.;          if ( t_orb->sunangle != t_orb->sunangle ) t_orb->sunangle = -1000.;
1883                if ( t_orb->sunmagangle != t_orb->sunmagangle ) t_orb->sunmagangle = -1000.;          if ( t_orb->sunmagangle != t_orb->sunmagangle ) t_orb->sunmagangle = -1000.;
1884                //          //
1885                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);
1886                //           //
1887                new(tor[nn]) OrbitalInfoTrkVar(*t_orb);           new(tor[nn]) OrbitalInfoTrkVar(*t_orb);
1888                nn++;          nn++;
1889                //           //
1890                t_orb->Clear();           t_orb->Clear();
1891                //          //
1892              }        }
1893              //        //
1894            } // end standard tracking algorithm      } // end standard tracking algorithm
1895    
1896            //            //
1897            // Code for extended tracking algorithm:            // Code for extended tracking algorithm:
# Line 1874  if(recqtime[mu]>=1160987921.75 && recqti Line 1923  if(recqtime[mu]>=1160987921.75 && recqti
1923                  t_orb->trkseqno = ptt->trkseqno;                  t_orb->trkseqno = ptt->trkseqno;
1924                  //                  //
1925                  TMatrixD Eij = PO->PamelatoGEO(Iij,Px,Py,Pz);                  TMatrixD Eij = PO->PamelatoGEO(Iij,Px,Py,Pz);
1926                  t_orb->Eij.ResizeTo(Eij);                        t_orb->Eij.ResizeTo(Eij);  
1927                  t_orb->Eij = Eij;                        t_orb->Eij = Eij;  
1928                  //                  //
1929                  TMatrixD Sij = PO->PamelatoGEO(Gij,Px,Py,Pz);                  TMatrixD Sij = PO->PamelatoGEO(Gij,Px,Py,Pz);
1930                  t_orb->Sij.ResizeTo(Sij);                  t_orb->Sij.ResizeTo(Sij);
1931                  t_orb->Sij = Sij;                  t_orb->Sij = Sij;
1932                  //                            //      
1933                  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);
1934                  //                  //
1935                  // SunPosition in instrumental reference frame                  // SunPosition in instrumental reference frame
# Line 1948  if(recqtime[mu]>=1160987921.75 && recqti Line 1997  if(recqtime[mu]>=1160987921.75 && recqti
1997                  t_orb->trkseqno = ptt->trkseqno;                  t_orb->trkseqno = ptt->trkseqno;
1998                  //                  //
1999                  TMatrixD Eij = PO->PamelatoGEO(Iij,Px,Py,Pz);                  TMatrixD Eij = PO->PamelatoGEO(Iij,Px,Py,Pz);
2000                  t_orb->Eij.ResizeTo(Eij);                        t_orb->Eij.ResizeTo(Eij);  
2001                  t_orb->Eij = Eij;                        t_orb->Eij = Eij;  
2002                  //                  //
2003                  TMatrixD Sij = PO->PamelatoGEO(Gij,Px,Py,Pz);                  TMatrixD Sij = PO->PamelatoGEO(Gij,Px,Py,Pz);
2004                  t_orb->Sij.ResizeTo(Sij);                  t_orb->Sij.ResizeTo(Sij);
2005                  t_orb->Sij = Sij;                  t_orb->Sij = Sij;
2006                  //                            //      
2007                  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);
2008                  //                  //
2009                  // SunPosition in instrumental reference frame                  // SunPosition in instrumental reference frame
# Line 2022  if(recqtime[mu]>=1160987921.75 && recqti Line 2071  if(recqtime[mu]>=1160987921.75 && recqti
2071                  t_orb->trkseqno = ptt->trkseqno;                  t_orb->trkseqno = ptt->trkseqno;
2072                  //                  //
2073                  TMatrixD Eij = PO->PamelatoGEO(Iij,Px,Py,Pz);                  TMatrixD Eij = PO->PamelatoGEO(Iij,Px,Py,Pz);
2074                  t_orb->Eij.ResizeTo(Eij);                        t_orb->Eij.ResizeTo(Eij);  
2075                  t_orb->Eij = Eij;                        t_orb->Eij = Eij;  
2076                  //                  //
2077                  TMatrixD Sij = PO->PamelatoGEO(Gij,Px,Py,Pz);                  TMatrixD Sij = PO->PamelatoGEO(Gij,Px,Py,Pz);
2078                  t_orb->Sij.ResizeTo(Sij);                  t_orb->Sij.ResizeTo(Sij);
2079                  t_orb->Sij = Sij;                  t_orb->Sij = Sij;
2080                  //                            //      
2081                  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);
2082                  //                  //
2083                  // SunPosition in instrumental reference frame                  // SunPosition in instrumental reference frame
# Line 2070  if(recqtime[mu]>=1160987921.75 && recqti Line 2119  if(recqtime[mu]>=1160987921.75 && recqti
2119              }              }
2120            } // end standard tracking algorithm: extended            } // end standard tracking algorithm: extended
2121    
2122          } else {    } else {
2123            if ( debug ) printf(" mmm... mode %u standalone  \n",orbitalinfo->mode);      if ( debug ) printf(" mmm... mode %u standalone  \n",orbitalinfo->mode);
2124          }    }
2125          //    //
2126        } else { // HERE IT MISS ALL CODE FOR EXTENDED TRACKING!        } else { // HERE IT MISS ALL CODE FOR EXTENDED TRACKING!
2127          if ( !standalone ){    if ( !standalone ){
2128            //      //
2129            if ( verbose ) printf(" no orb info for tracks \n");            if ( verbose ) printf(" no orb info for tracks \n");
2130            Int_t nn = 0;      Int_t nn = 0;
2131            for(Int_t nt=0; nt < tof->ntrk(); nt++){        for(Int_t nt=0; nt < tof->ntrk(); nt++){  
2132              //        //
2133              ToFTrkVar *ptt = tof->GetToFTrkVar(nt);        ToFTrkVar *ptt = tof->GetToFTrkVar(nt);
2134              if ( ptt->trkseqno != -1  ){        if ( ptt->trkseqno != -1  ){
2135                //          //
2136                t_orb->trkseqno = ptt->trkseqno;          t_orb->trkseqno = ptt->trkseqno;
2137                //          //
2138                t_orb->Eij = 0;            TMatrixD Iij(3,1);
2139                //          Iij(0,0)=0.; Iij(1,0)=0.; Iij(2,0)=1.;
2140                t_orb->Sij = 0;          //Iij(1,0)=0; Iij(1,1)=0; Iij(1,2)=0;
2141                //                      //Iij(2,0)=0; Iij(2,1)=0; Iij(2,2)=0;
2142                t_orb->pitch = -1000.;          //Iij.Zero();
2143                //          t_orb->Eij.ResizeTo(Iij);
2144                t_orb->sunangle = -1000.;          t_orb->Sij.ResizeTo(Iij);
2145                //          t_orb->Eij = Iij;
2146                t_orb->sunmagangle = -1000;          //
2147                //          t_orb->Sij = Iij;
2148                t_orb->cutoff = -1000.;          //      
2149                //          t_orb->pitch = -1000.;
2150                new(tor[nn]) OrbitalInfoTrkVar(*t_orb);          //
2151                nn++;          t_orb->sunangle = -1000.;
2152                //          //
2153                t_orb->Clear();          t_orb->sunmagangle = -1000;
2154                //          //
2155              }          t_orb->cutoff = -1000.;
2156              //          //
2157            }           new(tor[nn]) OrbitalInfoTrkVar(*t_orb);
2158            nn++;
2159             //
2160             t_orb->Clear();
2161            //
2162          }
2163          //
2164        }
2165            //            //
2166            // Code for extended tracking algorithm:            // Code for extended tracking algorithm:
2167            //            //
# Line 2119  if(recqtime[mu]>=1160987921.75 && recqti Line 2175  if(recqtime[mu]>=1160987921.75 && recqti
2175                  //                  //
2176                  t_orb->trkseqno = ptt->trkseqno;                  t_orb->trkseqno = ptt->trkseqno;
2177                  //                  //
2178                  t_orb->Eij = 0;                  TMatrixD Iij(3,1);
2179                    Iij(0,0)=0.; Iij(1,0)=0.; Iij(2,0)=1.;
2180                    //Iij(1,0)=0; Iij(1,1)=0; Iij(1,2)=0;
2181                    //Iij(2,0)=0; Iij(2,1)=0; Iij(2,2)=0;
2182                    //Iij.Zero();
2183                    t_orb->Eij.ResizeTo(Iij);
2184                    t_orb->Sij.ResizeTo(Iij);
2185                    t_orb->Eij = Iij;  
2186                  //                  //
2187                  t_orb->Sij = 0;                  t_orb->Sij = Iij;
2188                  //                            //      
2189                  t_orb->pitch = -1000.;                  t_orb->pitch = -1000.;
2190                  //                  //
2191                  t_orb->sunangle = -1000.;                  t_orb->sunangle = -1000.;
# Line 2152  if(recqtime[mu]>=1160987921.75 && recqti Line 2215  if(recqtime[mu]>=1160987921.75 && recqti
2215                  //                  //
2216                  t_orb->trkseqno = ptt->trkseqno;                  t_orb->trkseqno = ptt->trkseqno;
2217                  //                  //
2218                  t_orb->Eij = 0;                  TMatrixD Iij(3,1);
2219                    Iij(0,0)=0.; Iij(1,0)=0.; Iij(2,0)=1.;
2220                    //Iij(1,0)=0; Iij(1,1)=0; Iij(1,2)=0;
2221                    //Iij(2,0)=0; Iij(2,1)=0; Iij(2,2)=0;
2222                    //Iij.Zero();
2223                    t_orb->Eij.ResizeTo(Iij);
2224                    t_orb->Sij.ResizeTo(Iij);
2225                    t_orb->Eij = Iij;  
2226                  //                  //
2227                  t_orb->Sij = 0;                  t_orb->Sij = Iij;
2228                  //                            //      
2229                  t_orb->pitch = -1000.;                  t_orb->pitch = -1000.;
2230                  //                  //
2231                  t_orb->sunangle = -1000.;                  t_orb->sunangle = -1000.;
# Line 2175  if(recqtime[mu]>=1160987921.75 && recqti Line 2245  if(recqtime[mu]>=1160987921.75 && recqti
2245              }              }
2246            }            }
2247            if ( hasExtTrk ){            if ( hasExtTrk ){
2248              Int_t ttentry = 0;              Int_t ttentry = 0;
2249              if ( verbose ) printf(" hasExtTrk \n");              if ( verbose ) printf(" hasExtTrk \n");
2250              for(Int_t nt=0; nt < tcExtTof->GetEntries() ; nt++){                for(Int_t nt=0; nt < tcExtTof->GetEntries() ; nt++){  
2251                //                //
# Line 2185  if(recqtime[mu]>=1160987921.75 && recqti Line 2255  if(recqtime[mu]>=1160987921.75 && recqti
2255                  //                  //
2256                  t_orb->trkseqno = ptt->trkseqno;                  t_orb->trkseqno = ptt->trkseqno;
2257                  //                  //
2258                  t_orb->Eij = 0;                  TMatrixD Iij(3,1);
2259                    Iij(0,0)=0.; Iij(1,0)=0.; Iij(2,0)=1.;
2260                    //Iij(1,0)=0; Iij(1,1)=0; Iij(1,2)=0;
2261                    //Iij(2,0)=0; Iij(2,1)=0; Iij(2,2)=0;
2262                    //Iij.Zero();
2263                    t_orb->Eij.ResizeTo(Iij);
2264                    t_orb->Sij.ResizeTo(Iij);
2265                    t_orb->Eij = Iij;  
2266                  //                  //
2267                  t_orb->Sij = 0;                  t_orb->Sij = Iij;
2268                  //                            //      
2269                  t_orb->pitch = -1000.;                  t_orb->pitch = -1000.;
2270                  //                  //
2271                  t_orb->sunangle = -1000.;                  t_orb->sunangle = -1000.;
# Line 2197  if(recqtime[mu]>=1160987921.75 && recqti Line 2274  if(recqtime[mu]>=1160987921.75 && recqti
2274                  //                  //
2275                  t_orb->cutoff = -1000.;                  t_orb->cutoff = -1000.;
2276                  //                  //
2277                  TClonesArray &tz3 = *torbExtNucleiTrk;                  TClonesArray &tz3 = *torbExtTrk;
2278                  new(tz3[ttentry]) OrbitalInfoTrkVar(*t_orb);                  new(tz3[ttentry]) OrbitalInfoTrkVar(*t_orb);
2279                  ttentry++;                  ttentry++;
2280                  //                  //
# Line 2207  if(recqtime[mu]>=1160987921.75 && recqti Line 2284  if(recqtime[mu]>=1160987921.75 && recqti
2284                //                //
2285              }              }
2286            }            }
2287          }          }
2288        } // if( orbitalinfo->TimeGap>0){        }  // if( orbitalinfo->TimeGap>0){
2289        //        //
2290        // Fill the class        // Fill the class
2291        //        //
# Line 2296  if(recqtime[mu]>=1160987921.75 && recqti Line 2373  if(recqtime[mu]>=1160987921.75 && recqti
2373    if ( !reprocall && reproc && code >= 0 ){    if ( !reprocall && reproc && code >= 0 ){
2374      if ( totfileentries > noaftrun ){      if ( totfileentries > noaftrun ){
2375        if (verbose){        if (verbose){
2376          printf("\n Post-processing: copying events from the old tree after the processed run\n");      printf("\n Post-processing: copying events from the old tree after the processed run\n");  
2377          printf(" Copying %i events in the file which are after the end of the run %i \n",(int)(totfileentries-noaftrun),(int)run);    printf(" Copying %i events in the file which are after the end of the run %i \n",(int)(totfileentries-noaftrun),(int)run);
2378          printf(" Start copying at event number %i end copying at event number %i \n",(int)noaftrun,(int)totfileentries);    printf(" Start copying at event number %i end copying at event number %i \n",(int)noaftrun,(int)totfileentries);
2379        }        }
2380        for (UInt_t j = noaftrun; j < totfileentries; j++ ){        for (UInt_t j = noaftrun; j < totfileentries; j++ ){
2381          //    //
2382          // Get entry from old tree    // Get entry from old tree
2383          //    //
2384          if ( OrbitalInfotrclone->GetEntry(j) <= 0 ) throw -36;        if ( OrbitalInfotrclone->GetEntry(j) <= 0 ) throw -36;    
2385          //    //
2386          // copy orbitalinfoclone to OrbitalInfo    // copy orbitalinfoclone to OrbitalInfo
2387          //    //
2388          //      orbitalinfo->Clear();          //  orbitalinfo->Clear();
2389          //    //
2390          memcpy(&orbitalinfo,&orbitalinfoclone,sizeof(orbitalinfoclone));    memcpy(&orbitalinfo,&orbitalinfoclone,sizeof(orbitalinfoclone));
2391          //    //
2392          // Fill entry in the new tree    // Fill entry in the new tree
2393          //    //
2394          OrbitalInfotr->Fill();    OrbitalInfotr->Fill();
2395        };        };
2396        if (verbose) printf(" Finished successful copying!\n");        if (verbose) printf(" Finished successful copying!\n");
2397      };      };
# Line 2473  UInt_t holeq(Double_t lower,Double_t upp Line 2550  UInt_t holeq(Double_t lower,Double_t upp
2550      if((Qlower->time[5]-Qlower->time[0]<2)&&(Qlower->time[1]-Qlower->time[0]==30)){      if((Qlower->time[5]-Qlower->time[0]<2)&&(Qlower->time[1]-Qlower->time[0]==30)){
2551        //      mxtml = true;        //      mxtml = true;
2552        for(UInt_t i = 1; i < 6; i++){        for(UInt_t i = 1; i < 6; i++){
2553          if(Qlower->time[i]-Qlower->time[0]==30*i) NCQl=i;    if(Qlower->time[i]-Qlower->time[0]==30*i) NCQl=i;
2554        }        }
2555      }      }
2556      //    if((Qupper->time[5]-Qupper->time[0]<2)&&(Qupper->time[1]-Qupper->time[0]==30)){      //    if((Qupper->time[5]-Qupper->time[0]<2)&&(Qupper->time[1]-Qupper->time[0]==30)){
2557        //      mxtmu = true;        //      mxtmu = true;
2558        //      for(UInt_t i = 1; i < 6; i++){        //      for(UInt_t i = 1; i < 6; i++){
2559        //        if(Qupper->time[i]-Qupper->time[0]==30*i) NCQu=i;        //  if(Qupper->time[i]-Qupper->time[0]==30*i) NCQu=i;
2560        //      }        //      }
2561      //    }      //    }
2562    }    }
# Line 2522  void GM_ScanIGRF(TSQLServer *dbc, GMtype Line 2599  void GM_ScanIGRF(TSQLServer *dbc, GMtype
2599    if ( parerror<0 ) {    if ( parerror<0 ) {
2600      throw -902;      throw -902;
2601    }    }
2602          /*This function scans inputs G0, G1, and H1 of the IGRF table into 3 data arrays*/    /*This function scans inputs G0, G1, and H1 of the IGRF table into 3 data arrays*/
2603    //    TString SATH="/data03/Malakhov/pam9Malakhov/installed10/calib/orb-param/";    //  TString SATH="/data03/Malakhov/pam9Malakhov/installed10/calib/orb-param/";
2604          int i;    int i;
2605          double temp;    double temp;
2606          char buffer[200];    char buffer[200];
2607          FILE *IGRF;    FILE *IGRF;
2608          IGRF = fopen((glp->PATH+glp->NAME).Data(), "r");    IGRF = fopen((glp->PATH+glp->NAME).Data(), "r");
2609          //      IGRF = fopen(PATH+"IGRF.tab", "r");          //  IGRF = fopen(PATH+"IGRF.tab", "r");
2610          G0->size = 25;    G0->size = 25;
2611          G1->size = 25;    G1->size = 25;
2612          H1->size = 25;    H1->size = 25;
2613          for( i = 0; i < 4; i++)    for( i = 0; i < 4; i++)
2614          {    {
2615                  fgets(buffer, 200, IGRF);      fgets(buffer, 200, IGRF);
2616          }    }
2617          fscanf(IGRF, "g 1 0 %lf ", &G0->element[0]);    fscanf(IGRF, "g 1 0 %lf ", &G0->element[0]);
2618          for(i = 1; i <= 22; i++)    for(i = 1; i <= 22; i++)
2619          {    {
2620                  fscanf(IGRF ,"%lf ", &G0->element[i]);      fscanf(IGRF ,"%lf ", &G0->element[i]);
2621          }    }
2622          fscanf(IGRF ,"%lf\n", &temp);    fscanf(IGRF ,"%lf\n", &temp);
2623          G0->element[23] = temp * 5 + G0->element[22];    G0->element[23] = temp * 5 + G0->element[22];
2624          G0->element[24] = G0->element[23] + 5 * temp;    G0->element[24] = G0->element[23] + 5 * temp;
2625          fscanf(IGRF, "g 1 1 %lf ", &G1->element[0]);    fscanf(IGRF, "g 1 1 %lf ", &G1->element[0]);
2626          for(i = 1; i <= 22; i++)    for(i = 1; i <= 22; i++)
2627          {    {
2628                  fscanf( IGRF, "%lf ", &G1->element[i]);      fscanf( IGRF, "%lf ", &G1->element[i]);
2629          }    }
2630          fscanf(IGRF, "%lf\n", &temp);    fscanf(IGRF, "%lf\n", &temp);
2631          G1->element[23] = temp * 5 + G1->element[22];    G1->element[23] = temp * 5 + G1->element[22];
2632          G1->element[24] = temp * 5 + G1->element[23];    G1->element[24] = temp * 5 + G1->element[23];
2633          fscanf(IGRF, "h 1 1 %lf ", &H1->element[0]);    fscanf(IGRF, "h 1 1 %lf ", &H1->element[0]);
2634          for(i = 1; i <= 22; i++)    for(i = 1; i <= 22; i++)
2635          {    {
2636                  fscanf( IGRF, "%lf ", &H1->element[i]);      fscanf( IGRF, "%lf ", &H1->element[i]);
2637          }    }
2638          fscanf(IGRF, "%lf\n", &temp);    fscanf(IGRF, "%lf\n", &temp);
2639          H1->element[23] = temp * 5 + H1->element[22];    H1->element[23] = temp * 5 + H1->element[22];
2640          H1->element[24] = temp * 5 + H1->element[23];    H1->element[24] = temp * 5 + H1->element[23];
2641    if ( glp ) delete glp;    if ( glp ) delete glp;
2642    /*    /*
2643    printf("############################## SCAN IGRF ######################################\n");    printf("############################## SCAN IGRF ######################################\n");
# Line 2630  void GM_SetIGRF(Int_t isSecular, TString Line 2707  void GM_SetIGRF(Int_t isSecular, TString
2707    
2708  void GM_SetEllipsoid(GMtype_Ellipsoid *Ellip)  void GM_SetEllipsoid(GMtype_Ellipsoid *Ellip)
2709  {  {
2710          /*This function sets the WGS84 reference ellipsoid to its default values*/    /*This function sets the WGS84 reference ellipsoid to its default values*/
2711          Ellip->a        =                       6378.137; /*semi-major axis of the ellipsoid in */    Ellip->a  =      6378.137; /*semi-major axis of the ellipsoid in */
2712          Ellip->b        =                       6356.7523142;/*semi-minor axis of the ellipsoid in */    Ellip->b  =      6356.7523142;/*semi-minor axis of the ellipsoid in */
2713          Ellip->fla      =                       1/298.257223563;/* flattening */    Ellip->fla  =      1/298.257223563;/* flattening */
2714          Ellip->eps      =                       sqrt(1- ( Ellip->b *    Ellip->b) / (Ellip->a * Ellip->a ));  /*first eccentricity */    Ellip->eps  =      sqrt(1- ( Ellip->b *  Ellip->b) / (Ellip->a * Ellip->a ));  /*first eccentricity */
2715          Ellip->epssq    =                       (Ellip->eps * Ellip->eps);   /*first eccentricity squared */    Ellip->epssq  =      (Ellip->eps * Ellip->eps);   /*first eccentricity squared */
2716          Ellip->re       =                       6371.2;/* Earth's radius */    Ellip->re  =      6371.2;/* Earth's radius */
2717  } /*GM_SetEllipsoid*/  } /*GM_SetEllipsoid*/
2718    
2719    
2720  void GM_EarthCartToDipoleCartCD(GMtype_Pole Pole, GMtype_CoordCartesian EarthCoord, GMtype_CoordCartesian *DipoleCoords)  void GM_EarthCartToDipoleCartCD(GMtype_Pole Pole, GMtype_CoordCartesian EarthCoord, GMtype_CoordCartesian *DipoleCoords)
2721  {  {
2722          /*This function converts from Earth centered cartesian coordinates to dipole centered cartesian coordinates*/    /*This function converts from Earth centered cartesian coordinates to dipole centered cartesian coordinates*/
2723          double X, Y, Z, CosPhi, SinPhi, CosLambda, SinLambda;    double X, Y, Z, CosPhi, SinPhi, CosLambda, SinLambda;
2724          CosPhi = cos(TMath::DegToRad()*Pole.phi);    CosPhi = cos(TMath::DegToRad()*Pole.phi);
2725          SinPhi = sin(TMath::DegToRad()*Pole.phi);    SinPhi = sin(TMath::DegToRad()*Pole.phi);
2726          CosLambda = cos(TMath::DegToRad()*Pole.lambda);    CosLambda = cos(TMath::DegToRad()*Pole.lambda);
2727          SinLambda = sin(TMath::DegToRad()*Pole.lambda);    SinLambda = sin(TMath::DegToRad()*Pole.lambda);
2728          X = EarthCoord.x;    X = EarthCoord.x;
2729          Y = EarthCoord.y;    Y = EarthCoord.y;
2730          Z = EarthCoord.z;    Z = EarthCoord.z;
2731              
2732          /*These equations are taken from a document by Wallace H. Campbell*/    /*These equations are taken from a document by Wallace H. Campbell*/
2733          DipoleCoords->x = X * CosPhi * CosLambda + Y * CosPhi * SinLambda - Z * SinPhi;    DipoleCoords->x = X * CosPhi * CosLambda + Y * CosPhi * SinLambda - Z * SinPhi;
2734          DipoleCoords->y = -X * SinLambda + Y * CosLambda;    DipoleCoords->y = -X * SinLambda + Y * CosLambda;
2735          DipoleCoords->z = X * SinPhi * CosLambda + Y * SinPhi * SinLambda + Z * CosPhi;    DipoleCoords->z = X * SinPhi * CosLambda + Y * SinPhi * SinLambda + Z * CosPhi;
2736  } /*GM_EarthCartToDipoleCartCD*/  } /*GM_EarthCartToDipoleCartCD*/
2737    
2738  void GM_GeodeticToSpherical(GMtype_Ellipsoid Ellip, GMtype_CoordGeodetic CoordGeodetic, GMtype_CoordSpherical *CoordSpherical)  void GM_GeodeticToSpherical(GMtype_Ellipsoid Ellip, GMtype_CoordGeodetic CoordGeodetic, GMtype_CoordSpherical *CoordSpherical)
2739  {  {
2740          double CosLat, SinLat, rc, xp, zp; /*all local variables */    double CosLat, SinLat, rc, xp, zp; /*all local variables */
2741          /*    /*
2742          ** Convert geodetic coordinates, (defined by the WGS-84    ** Convert geodetic coordinates, (defined by the WGS-84
2743          ** reference ellipsoid), to Earth Centered Earth Fixed Cartesian    ** reference ellipsoid), to Earth Centered Earth Fixed Cartesian
2744          ** coordinates, and then to spherical coordinates.    ** coordinates, and then to spherical coordinates.
2745          */    */
2746    
2747          CosLat = cos(TMath::DegToRad()*CoordGeodetic.phi);    CosLat = cos(TMath::DegToRad()*CoordGeodetic.phi);
2748          SinLat = sin(TMath::DegToRad()*CoordGeodetic.phi);    SinLat = sin(TMath::DegToRad()*CoordGeodetic.phi);
2749    
2750          /* compute the local radius of curvature on the WGS-84 reference ellipsoid */    /* compute the local radius of curvature on the WGS-84 reference ellipsoid */
2751    
2752          rc = Ellip.a / sqrt(1.0 - Ellip.epssq * SinLat * SinLat);    rc = Ellip.a / sqrt(1.0 - Ellip.epssq * SinLat * SinLat);
2753    
2754          /* compute ECEF Cartesian coordinates of specified point (for longitude=0) */    /* compute ECEF Cartesian coordinates of specified point (for longitude=0) */
2755    
2756          xp = (rc + CoordGeodetic.HeightAboveEllipsoid) * CosLat;    xp = (rc + CoordGeodetic.HeightAboveEllipsoid) * CosLat;
2757          zp = (rc*(1.0 - Ellip.epssq) + CoordGeodetic.HeightAboveEllipsoid) * SinLat;    zp = (rc*(1.0 - Ellip.epssq) + CoordGeodetic.HeightAboveEllipsoid) * SinLat;
2758    
2759          /* compute spherical radius and angle lambda and phi of specified point */    /* compute spherical radius and angle lambda and phi of specified point */
2760    
2761          CoordSpherical->r = sqrt(xp * xp + zp * zp);    CoordSpherical->r = sqrt(xp * xp + zp * zp);
2762          CoordSpherical->phig = TMath::RadToDeg()*asin(zp / CoordSpherical->r);     /* geocentric latitude */    CoordSpherical->phig = TMath::RadToDeg()*asin(zp / CoordSpherical->r);     /* geocentric latitude */
2763          CoordSpherical->lambda = CoordGeodetic.lambda;                   /* longitude */    CoordSpherical->lambda = CoordGeodetic.lambda;                   /* longitude */
2764  } /*GM_GeodeticToSpherical*/  } /*GM_GeodeticToSpherical*/
2765    
2766  void GM_PoleLocation(GMtype_Model Model, GMtype_Pole *Pole)  void GM_PoleLocation(GMtype_Model Model, GMtype_Pole *Pole)
2767  {  {
2768          /*This function finds the location of the north magnetic pole in spherical coordinates.  The equations are    /*This function finds the location of the north magnetic pole in spherical coordinates.  The equations are
2769          **from Wallace H. Campbell's Introduction to Geomagnetic Fields*/    **from Wallace H. Campbell's Introduction to Geomagnetic Fields*/
2770    
2771          Pole->phi = TMath::RadToDeg()*-atan(sqrt(Model.h1 * Model.h1 + Model.g1 * Model.g1)/Model.g0);    Pole->phi = TMath::RadToDeg()*-atan(sqrt(Model.h1 * Model.h1 + Model.g1 * Model.g1)/Model.g0);
2772          Pole->lambda = TMath::RadToDeg()*atan(Model.h1/Model.g1);    Pole->lambda = TMath::RadToDeg()*atan(Model.h1/Model.g1);
2773  } /*GM_PoleLocation*/  } /*GM_PoleLocation*/
2774    
2775  void GM_SphericalToCartesian(GMtype_CoordSpherical CoordSpherical, GMtype_CoordCartesian *CoordCartesian)  void GM_SphericalToCartesian(GMtype_CoordSpherical CoordSpherical, GMtype_CoordCartesian *CoordCartesian)
2776  {  {
2777          /*This function converts spherical coordinates into Cartesian coordinates*/    /*This function converts spherical coordinates into Cartesian coordinates*/
2778          double CosPhi = cos(TMath::DegToRad()*CoordSpherical.phig);    double CosPhi = cos(TMath::DegToRad()*CoordSpherical.phig);
2779          double SinPhi = sin(TMath::DegToRad()*CoordSpherical.phig);    double SinPhi = sin(TMath::DegToRad()*CoordSpherical.phig);
2780          double CosLambda = cos(TMath::DegToRad()*CoordSpherical.lambda);    double CosLambda = cos(TMath::DegToRad()*CoordSpherical.lambda);
2781          double SinLambda = sin(TMath::DegToRad()*CoordSpherical.lambda);    double SinLambda = sin(TMath::DegToRad()*CoordSpherical.lambda);
2782              
2783          CoordCartesian->x = CoordSpherical.r * CosPhi * CosLambda;    CoordCartesian->x = CoordSpherical.r * CosPhi * CosLambda;
2784          CoordCartesian->y = CoordSpherical.r * CosPhi * SinLambda;    CoordCartesian->y = CoordSpherical.r * CosPhi * SinLambda;
2785          CoordCartesian->z = CoordSpherical.r * SinPhi;    CoordCartesian->z = CoordSpherical.r * SinPhi;
2786  } /*GM_SphericalToCartesian*/  } /*GM_SphericalToCartesian*/
2787    
2788  void GM_TimeAdjustCoefs(Float_t year, Float_t jyear, GMtype_Data g0d, GMtype_Data g1d, GMtype_Data h1d, GMtype_Model *Model)  void GM_TimeAdjustCoefs(Float_t year, Float_t jyear, GMtype_Data g0d, GMtype_Data g1d, GMtype_Data h1d, GMtype_Model *Model)
2789  {  {
2790          /*This function calls GM_LinearInterpolation for the coefficients to estimate the value of the    /*This function calls GM_LinearInterpolation for the coefficients to estimate the value of the
2791          **coefficient for the given date*/    **coefficient for the given date*/
2792          int index;    int index;
2793          double x;    double x;
2794          index = (year - GM_STARTYEAR) / 5;    index = (year - GM_STARTYEAR) / 5;
2795          x = (jyear - GM_STARTYEAR) / 5.;    x = (jyear - GM_STARTYEAR) / 5.;
2796          Model->g0 = GM_LinearInterpolation(index, index+1, g0d.element[index], g0d.element[index+1], x);    Model->g0 = GM_LinearInterpolation(index, index+1, g0d.element[index], g0d.element[index+1], x);
2797          Model->g1 = GM_LinearInterpolation(index, index+1, g1d.element[index], g1d.element[index+1], x);    Model->g1 = GM_LinearInterpolation(index, index+1, g1d.element[index], g1d.element[index+1], x);
2798          Model->h1 = GM_LinearInterpolation(index, index+1, h1d.element[index], h1d.element[index+1], x);    Model->h1 = GM_LinearInterpolation(index, index+1, h1d.element[index], h1d.element[index+1], x);
2799  } /*GM_TimeAdjustCoefs*/  } /*GM_TimeAdjustCoefs*/
2800    
2801  double GM_LinearInterpolation(double x1, double x2, double y1, double y2, double x)  double GM_LinearInterpolation(double x1, double x2, double y1, double y2, double x)
2802  {  {
2803          /*This function takes a linear interpolation between two given points for x*/    /*This function takes a linear interpolation between two given points for x*/
2804          double weight, y;    double weight, y;
2805          weight  = (x - x1) / (x2 - x1);    weight  = (x - x1) / (x2 - x1);
2806          y = y1 * (1. - weight) + y2 * weight;    y = y1 * (1. - weight) + y2 * weight;
2807          //        printf(" x1 %f x2 %f y1 %f y2 %f x %f ==> y %f \n",x1,x2,y1,y2,x,y);          //        printf(" x1 %f x2 %f y1 %f y2 %f x %f ==> y %f \n",x1,x2,y1,y2,x,y);
2808          return y;    return y;
2809  }/*GM_LinearInterpolation*/  }/*GM_LinearInterpolation*/
2810    
2811  void GM_CartesianToSpherical(GMtype_CoordCartesian CoordCartesian, GMtype_CoordSpherical *CoordSpherical)  void GM_CartesianToSpherical(GMtype_CoordCartesian CoordCartesian, GMtype_CoordSpherical *CoordSpherical)
2812  {  {
2813          /*This function converts a point from Cartesian coordinates into spherical coordinates*/    /*This function converts a point from Cartesian coordinates into spherical coordinates*/
2814          double X, Y, Z;    double X, Y, Z;
2815              
2816          X = CoordCartesian.x;    X = CoordCartesian.x;
2817          Y = CoordCartesian.y;    Y = CoordCartesian.y;
2818          Z = CoordCartesian.z;    Z = CoordCartesian.z;
2819    
2820          CoordSpherical->r = sqrt(X * X + Y * Y + Z * Z);    CoordSpherical->r = sqrt(X * X + Y * Y + Z * Z);
2821          CoordSpherical->phig = TMath::RadToDeg()*asin(Z / (CoordSpherical->r));    CoordSpherical->phig = TMath::RadToDeg()*asin(Z / (CoordSpherical->r));
2822          CoordSpherical->lambda = TMath::RadToDeg()*atan2(Y, X);    CoordSpherical->lambda = TMath::RadToDeg()*atan2(Y, X);
2823  } /*GM_CartesianToSpherical*/  } /*GM_CartesianToSpherical*/
Legend:
Removed from v.1.84  
changed lines
  Added in v.1.93
  ViewVC Help
Powered by ViewVC 1.1.23