/[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.39 by mocchiut, Wed Mar 18 08:54:41 2009 UTC revision 1.40 by mocchiut, Tue Aug 4 14:01:26 2009 UTC
# Line 44  Line 44 
44  #include <OrbitalInfoCore.h>  #include <OrbitalInfoCore.h>
45  #include <InclinationInfo.h>  #include <InclinationInfo.h>
46    
47    
48  using namespace std;  using namespace std;
49    
50  //  //
# Line 177  int OrbitalInfoCore(UInt_t run, TFile *f Line 178  int OrbitalInfoCore(UInt_t run, TFile *f
178    //    //
179    // IGRF stuff    // IGRF stuff
180    //    //
181    float dimo = 0.0; // dipole moment (computed from dat files)    Float_t dimo = 0.0; // dipole moment (computed from dat files)
182    float bnorth, beast, bdown, babs;    Float_t bnorth, beast, bdown, babs;
183    float xl; // L value    Float_t xl; // L value
184    float icode; // code value for L accuracy (see fortran code)    Float_t icode; // code value for L accuracy (see fortran code)
185    float bab1; // What's  the difference with babs?    Float_t bab1; // What's  the difference with babs?
186    float stps = 0.005; // step size for field line tracing    Float_t stps = 0.005; // step size for field line tracing
187    float bdel = 0.01; // required accuracy    Float_t bdel = 0.01; // required accuracy
188    float bequ;  // equatorial b value (also called b_0)    Float_t bequ;  // equatorial b value (also called b_0)
189    bool value = 0; // false if bequ is not the minimum b value    Bool_t value = 0; // false if bequ is not the minimum b value
190    float rr0; // equatorial radius normalized to earth radius    Float_t rr0; // equatorial radius normalized to earth radius
191    
192    //    //
193    // Working filename    // Working filename
# Line 210  int OrbitalInfoCore(UInt_t run, TFile *f Line 211  int OrbitalInfoCore(UInt_t run, TFile *f
211    OrbitalInfofolder << tempname.str().c_str();    OrbitalInfofolder << tempname.str().c_str();
212    tempname << "/OrbitalInfotree_run";    tempname << "/OrbitalInfotree_run";
213    tempname << run << ".root";      tempname << run << ".root";  
214      UInt_t totnorun = 0;
215    //    //
216    // DB classes    // DB classes
217    //    //
# Line 343  int OrbitalInfoCore(UInt_t run, TFile *f Line 345  int OrbitalInfoCore(UInt_t run, TFile *f
345    // number of run to be processed    // number of run to be processed
346    //    //
347    numbofrun = runinfo->GetNoRun();    numbofrun = runinfo->GetNoRun();
348    UInt_t totnorun = runinfo->GetRunEntries();    totnorun = runinfo->GetRunEntries();
349    //    //
350    // Try to access the OrbitalInfo tree in the file, if it exists we are reprocessing data if not we are processing a new run    // Try to access the OrbitalInfo tree in the file, if it exists we are reprocessing data if not we are processing a new run
351    //    //
# Line 503  int OrbitalInfoCore(UInt_t run, TFile *f Line 505  int OrbitalInfoCore(UInt_t run, TFile *f
505      fname = ftmpname.str().c_str();      fname = ftmpname.str().c_str();
506      ftmpname.str("");      ftmpname.str("");
507      //      //
508      // print out informations      // print nout informations
509      //      //
510      totevent = runinfo->NEVENTS;      totevent = runinfo->NEVENTS;
511      evfrom = runinfo->EV_FROM;      evfrom = runinfo->EV_FROM;
# Line 706  int OrbitalInfoCore(UInt_t run, TFile *f Line 708  int OrbitalInfoCore(UInt_t run, TFile *f
708        //        //
709        ph = eh->GetPscuHeader();        ph = eh->GetPscuHeader();
710        atime = dbtime->DBabsTime(ph->GetOrbitalTime());        atime = dbtime->DBabsTime(ph->GetOrbitalTime());
711          if ( debug ) printf(" %i absolute time \n",procev);      
712        //        //
713        // paranoid check        // paranoid check
714        //        //
# Line 745  int OrbitalInfoCore(UInt_t run, TFile *f Line 748  int OrbitalInfoCore(UInt_t run, TFile *f
748        //        //
749        // start processing        // start processing
750        //        //
751          if ( debug ) printf(" %i start processing \n",procev);      
752        orbitalinfo->Clear();        orbitalinfo->Clear();
753        //        //
754        OrbitalInfoTrkVar *t_orb = new OrbitalInfoTrkVar();        OrbitalInfoTrkVar *t_orb = new OrbitalInfoTrkVar();
# Line 756  int OrbitalInfoCore(UInt_t run, TFile *f Line 760  int OrbitalInfoCore(UInt_t run, TFile *f
760        orbitalinfo->pkt_num = ph->GetCounter();        orbitalinfo->pkt_num = ph->GetCounter();
761        orbitalinfo->OBT = ph->GetOrbitalTime();        orbitalinfo->OBT = ph->GetOrbitalTime();
762        orbitalinfo->absTime = atime;        orbitalinfo->absTime = atime;
763          if ( debug ) printf(" %i pktnum obt abstime \n",procev);      
764        //        //
765        // Propagate the orbit from the tle time to atime, using SGP(D)4.        // Propagate the orbit from the tle time to atime, using SGP(D)4.
766        //        //
767          if ( debug ) printf(" %i sgp4 \n",procev);      
768        cCoordGeo coo;        cCoordGeo coo;
769        float jyear=0;            Float_t jyear=0.;    
770        //        //
771        if(atime >= gltle->GetToTime()) {        if(atime >= gltle->GetToTime()) {
772          if ( !gltle->Query(atime, dbc) ){          if ( !gltle->Query(atime, dbc) ){
773            //                  //      
774            // Compute the magnetic dipole moment.            // Compute the magnetic dipole moment.
775            //            //
776              if ( debug ) printf(" %i compute magnetic dipole moment \n",procev);      
777            UInt_t year, month, day, hour, min, sec;            UInt_t year, month, day, hour, min, sec;
778            //            //
779            TTimeStamp t = TTimeStamp(atime, kTRUE);            TTimeStamp t = TTimeStamp(atime, kTRUE);
# Line 776  int OrbitalInfoCore(UInt_t run, TFile *f Line 783  int OrbitalInfoCore(UInt_t run, TFile *f
783              + (month*31.+ (float) day)/365.              + (month*31.+ (float) day)/365.
784              + (hour*3600.+min*60.+(float)sec)/(24*3600*365.);              + (hour*3600.+min*60.+(float)sec)/(24*3600*365.);
785            //            //
786              if ( debug ) printf(" %i compute magnetic dipole moment get dipole moment for year\n",procev);      
787            feldcof_(&jyear, &dimo); // get dipole moment for year            feldcof_(&jyear, &dimo); // get dipole moment for year
788              if ( debug ) printf(" %i compute magnetic dipole moment end\n",procev);      
789          } else {          } else {
790            code = -56;            code = -56;
791            goto closeandexit;            goto closeandexit;
# Line 878  int OrbitalInfoCore(UInt_t run, TFile *f Line 887  int OrbitalInfoCore(UInt_t run, TFile *f
887          lowerqtime = upperqtime;          lowerqtime = upperqtime;
888          Long64_t maxloop = 100000000LL;          Long64_t maxloop = 100000000LL;
889          Long64_t mn = 0;          Long64_t mn = 0;
890          bool gh=false;          Bool_t gh=false;
891          ooi=oi;          ooi=oi;
892          if ( verbose ) printf(" OrbitalInfoCore: sync with quaternions data upperqtime %u lowerqtime %u atime %u \n",(UInt_t)upperqtime,(UInt_t)lowerqtime,atime);          if ( verbose ) printf(" OrbitalInfoCore: sync with quaternions data upperqtime %u lowerqtime %u atime %u \n",(UInt_t)upperqtime,(UInt_t)lowerqtime,atime);
893          while (!gh){                while (!gh){      
# Line 1077  int OrbitalInfoCore(UInt_t run, TFile *f Line 1086  int OrbitalInfoCore(UInt_t run, TFile *f
1086          //          //
1087        };              };      
1088        //        //
1089          if ( debug ) printf(" pitch angle \n");
1090          //
1091        // pitch angles        // pitch angles
1092        //        //
1093        if ( orbitalinfo->mode != 10 && orbitalinfo->mode != 5 && orbitalinfo->mode !=7 && orbitalinfo->mode != 9 ){        if ( orbitalinfo->mode != 10 && orbitalinfo->mode != 5 && orbitalinfo->mode !=7 && orbitalinfo->mode != 9 ){
# Line 1312  void CopyAng(InclinationInfo *A1, Inclin Line 1323  void CopyAng(InclinationInfo *A1, Inclin
1323  UInt_t holeq(Double_t lower,Double_t upper,Quaternions *Qlower, Quaternions *Qupper, UInt_t f){  UInt_t holeq(Double_t lower,Double_t upper,Quaternions *Qlower, Quaternions *Qupper, UInt_t f){
1324        
1325    UInt_t hole = 10;    UInt_t hole = 10;
1326    bool R10l = false;     // Sign of R10 mode in lower quaternions array    Bool_t R10l = false;     // Sign of R10 mode in lower quaternions array
1327    bool R10u = false;     // Sign of R10 mode in upper quaternions array    Bool_t R10u = false;     // Sign of R10 mode in upper quaternions array
1328    bool insm = false;     // Sign that we inside quaternions array    Bool_t insm = false;     // Sign that we inside quaternions array
1329    bool mxtml = false;    // Sign of mixt mode in lower quaternions array    Bool_t mxtml = false;    // Sign of mixt mode in lower quaternions array
1330    bool mxtmu = false;    // Sign of mixt mode in upper quaternions array    Bool_t mxtmu = false;    // Sign of mixt mode in upper quaternions array
1331    bool npasm = false;     // Sign of normall pass between R10 and non R10 or between non R10 and R10    Bool_t npasm = false;     // Sign of normall pass between R10 and non R10 or between non R10 and R10
1332    UInt_t NCQl = 6;       // Number of correct quaternions in lower array    UInt_t NCQl = 6;       // Number of correct quaternions in lower array
1333    UInt_t NCQu = 6;       // Number of correct quaternions in upper array    UInt_t NCQu = 6;       // Number of correct quaternions in upper array
1334    if (f>0){    if (f>0){
Legend:
Removed from v.1.39  
changed lines
  Added in v.1.40
  ViewVC Help
Powered by ViewVC 1.1.23