/[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.47 by mocchiut, Mon Feb 1 05:34:39 2010 UTC revision 1.74 by malakhov, Wed Jul 9 09:11:55 2014 UTC
# Line 9  Line 9 
9  // ROOT headers  // ROOT headers
10  //  //
11  //#include <TCanvas.h>  //#include <TCanvas.h>
12  //#include <TH2F.h> //for test only. Vitaly.  #include <TH2F.h> //for test only. Vitaly.
13    #include <TVector3.h>
14  //#include <TF1.h>  //#include <TF1.h>
15    
16  #include <TTree.h>  #include <TTree.h>
# Line 47  Line 48 
48  #include <OrbitalInfoCore.h>  #include <OrbitalInfoCore.h>
49  #include <InclinationInfo.h>  #include <InclinationInfo.h>
50    
51    //
52    // Tracker and ToF classes headers and definitions
53    //
54    #include <ToFLevel2.h>
55    #include <TrkLevel2.h>
56    #include <ExtTrack.h> // new tracking code
57    
58  using namespace std;  using namespace std;
59    
# Line 65  int OrbitalInfoCore(UInt_t run, TFile *f Line 72  int OrbitalInfoCore(UInt_t run, TFile *f
72    stringstream myquery;    stringstream myquery;
73    myquery.str("");    myquery.str("");
74    myquery << "SET time_zone='+0:00'";    myquery << "SET time_zone='+0:00'";
75    dbc->Query(myquery.str().c_str());    delete dbc->Query(myquery.str().c_str());
76    //    //
77    TString processFolder = Form("OrbitalInfoFolder_%u",run);    TString processFolder = Form("OrbitalInfoFolder_%u",run);
78    //    //
# Line 121  int OrbitalInfoCore(UInt_t run, TFile *f Line 128  int OrbitalInfoCore(UInt_t run, TFile *f
128    TTree *OrbitalInfotrclone = 0;    TTree *OrbitalInfotrclone = 0;
129    Bool_t reproc = false;    Bool_t reproc = false;
130    Bool_t reprocall = false;    Bool_t reprocall = false;
131      Bool_t igrfloaded = false;
132    UInt_t nobefrun = 0;    UInt_t nobefrun = 0;
133    UInt_t noaftrun = 0;    UInt_t noaftrun = 0;
134    UInt_t numbofrun = 0;    UInt_t numbofrun = 0;
# Line 128  int OrbitalInfoCore(UInt_t run, TFile *f Line 136  int OrbitalInfoCore(UInt_t run, TFile *f
136    TString fname;    TString fname;
137    UInt_t totfileentries = 0;    UInt_t totfileentries = 0;
138    UInt_t idRun = 0;    UInt_t idRun = 0;
139      UInt_t anni5 = 60 * 60 * 24 * 365 * 5 ;//1576800
140    //    //
141    // My variables. Vitaly.    // My variables. Vitaly.
142    //    //
# Line 180  int OrbitalInfoCore(UInt_t run, TFile *f Line 189  int OrbitalInfoCore(UInt_t run, TFile *f
189    //    //
190    // IGRF stuff    // IGRF stuff
191    //    //
192    Float_t dimo = 0.0; // dipole moment (computed from dat files)    Double_t dimo = 0.0; // dipole moment (computed from dat files) // EM GCC 4.7
193    Float_t bnorth, beast, bdown, babs;    Float_t bnorth, beast, bdown, babs;
194    Float_t xl; // L value    Float_t xl; // L value
195    Float_t icode; // code value for L accuracy (see fortran code)    Float_t icode; // code value for L accuracy (see fortran code)
# Line 223  int OrbitalInfoCore(UInt_t run, TFile *f Line 232  int OrbitalInfoCore(UInt_t run, TFile *f
232    //    //
233    //Quaternions classes    //Quaternions classes
234    //    //
235    Quaternions *L_QQ_Q_l_lower = new Quaternions();    Quaternions *L_QQ_Q_l_lower = 0;
236    InclinationInfo *RYPang_lower = new InclinationInfo();    InclinationInfo *RYPang_lower = 0;
237    Quaternions *L_QQ_Q_l_upper = new Quaternions();    Quaternions *L_QQ_Q_l_upper = 0;
238    InclinationInfo *RYPang_upper = new InclinationInfo();    InclinationInfo *RYPang_upper = 0;
239        
240    cEci eCi;    cEci eCi;
241        
242    // Initialize fortran routines!!!    // Initialize fortran routines!!!
243      Int_t ltp1 = 0;
244    Int_t ltp2 = 0;    Int_t ltp2 = 0;
245    Int_t ltp3 = 0;    Int_t ltp3 = 0;
246    Int_t uno = 1;    //  Int_t uno = 1;
247    const char *niente = " ";    //  const char *niente = " ";
248      GL_PARAM *glparam0 = new GL_PARAM();
249    GL_PARAM *glparam = new GL_PARAM();    GL_PARAM *glparam = new GL_PARAM();
250    GL_PARAM *glparam2 = new GL_PARAM();    GL_PARAM *glparam2 = new GL_PARAM();
251      GL_PARAM *glparam3 = new GL_PARAM();
252    
253    //    //
254    // Orientation variables. Vitaly    // Orientation variables. Vitaly
255    //    //
256    
257    UInt_t evfrom = 0;    UInt_t evfrom = 0;
258    UInt_t jumped = 0;    UInt_t jumped = 0;
259    Int_t itr = -1;        Int_t itr = -1;    
260    Double_t A1;    //  Double_t A1;
261    Double_t A2;    //  Double_t A2;
262    Double_t A3;    //  Double_t A3;
263    Double_t Px = 0;    Double_t Px = 0;
264    Double_t Py = 0;          Double_t Py = 0;      
265    Double_t Pz = 0;      Double_t Pz = 0;  
266    TTree *ttof = 0;    TTree *ttof = 0;
267    ToFLevel2 *tof = new ToFLevel2();    ToFLevel2 *tof = new ToFLevel2();
268      TTree *ttrke = 0;
269      TrkLevel2 *trke = new TrkLevel2();
270    OrientationInfo *PO = new OrientationInfo();    OrientationInfo *PO = new OrientationInfo();
271    Int_t nz = 6;    Int_t nz = 6;
272    Float_t zin[6];    Float_t zin[6];
273    Int_t nevtofl2 = 0;    Int_t nevtofl2 = 0;
274      Int_t nevtrkl2 = 0;
275    if ( verbose ) cout<<"Reading quaternions external file"<<endl;    if ( verbose ) cout<<"Reading quaternions external file"<<endl;
276    cout.setf(ios::fixed,ios::floatfield);      cout.setf(ios::fixed,ios::floatfield);  
277    /******Reading recovered quaternions...*********/    /******Reading recovered quaternions...*********/
# Line 265  int OrbitalInfoCore(UInt_t run, TFile *f Line 281  int OrbitalInfoCore(UInt_t run, TFile *f
281    vector<Float_t> recq2;    vector<Float_t> recq2;
282    vector<Float_t> recq3;    vector<Float_t> recq3;
283    Float_t Norm = 1;    Float_t Norm = 1;
284    Int_t parerror=glparam->Query_GL_PARAM(1,303,dbc); // parameters stored in DB in GL_PRAM table    
285    ifstream in((glparam->PATH+glparam->NAME).Data(),ios::in);    vector<UInt_t> RTtime1;
286      vector<UInt_t> RTtime2;
287      vector<Double_t> RTbank1;
288      vector<Double_t> RTbank2;
289      vector<Double_t> RTbpluto1;
290      vector<Double_t> RTbpluto2;
291      vector<Int_t> RTazim;
292      vector<Int_t> RTstart;
293      vector<Int_t> RTpluto2;
294      vector<Int_t> RTpluto1;
295      vector<Int_t> RTerrq;
296    
297      TClonesArray *tcNucleiTrk = NULL;
298      TClonesArray *tcExtNucleiTrk = NULL;
299      TClonesArray *tcExtTrk = NULL;
300      TClonesArray *tcNucleiTof = NULL;
301      TClonesArray *tcExtNucleiTof = NULL;
302      TClonesArray *tcExtTof = NULL;
303      TClonesArray *torbNucleiTrk = NULL;
304      TClonesArray *torbExtNucleiTrk = NULL;
305      TClonesArray *torbExtTrk = NULL;
306      Bool_t hasNucleiTrk = false;
307      Bool_t hasExtNucleiTrk = false;
308      Bool_t hasExtTrk = false;
309      Bool_t hasNucleiTof = false;
310      Bool_t hasExtNucleiTof = false;
311      Bool_t hasExtTof = false;
312    
313      ifstream in;
314      ifstream an;
315      //  ofstream mc;
316      //  TString gr;
317      Int_t parerror2=0;
318    
319      Int_t parerror=glparam0->Query_GL_PARAM(1,303,dbc); // parameters stored in DB in GL_PRAM table
320      if ( verbose ) cout<<parerror<<"\t"<<(char*)(glparam0->PATH+glparam0->NAME).Data()<<endl;
321    if ( parerror<0 ) {    if ( parerror<0 ) {
322      code = parerror;      code = parerror;
323      goto closeandexit;      goto closeandexit;
324    };    }
325      in.open((char*)(glparam0->PATH+glparam0->NAME).Data(),ios::in);
326    while(!in.eof()){    while(!in.eof()){
327      recqtime.resize(recqtime.size()+1);      recqtime.resize(recqtime.size()+1);
328      Int_t sizee = recqtime.size();      Int_t sizee = recqtime.size();
# Line 285  int OrbitalInfoCore(UInt_t run, TFile *f Line 337  int OrbitalInfoCore(UInt_t run, TFile *f
337      in>>recq3[sizee-1];      in>>recq3[sizee-1];
338      in>>Norm;      in>>Norm;
339    }    }
340      in.close();
341    if ( verbose ) cout<<"We have read recovered data"<<endl;    if ( verbose ) cout<<"We have read recovered data"<<endl;
342      if (debug) cout << "size of recovered quaterions data set is " << recqtime.size() << endl;
343    
344      if ( verbose ) cout<<"read Rotation Table"<<endl;
345      
346      parerror2=glparam0->Query_GL_PARAM(1,305,dbc);
347    
348    parerror=glparam->Query_GL_PARAM(1,301,dbc); // parameters stored in DB in GL_PRAM table      if ( verbose ) cout<<parerror2<<"\t"<<(char*)(glparam0->PATH+glparam0->NAME).Data()<<endl;
349    if ( parerror<0 ) {    if ( parerror2<0 ) {
350      code = parerror;      code = parerror;
351      goto closeandexit;      goto closeandexit;
352    };    }
353    ltp2 = (Int_t)(glparam->PATH+glparam->NAME).Length();    an.open((char*)(glparam0->PATH+glparam0->NAME).Data(),ios::in);
354    if ( verbose ) printf(" Reading Earth's Magnetic Field parameter file: %s \n",(glparam->PATH+glparam->NAME).Data());    while(!an.eof()){
355    //      RTtime1.resize(RTtime1.size()+1);
356    parerror=glparam2->Query_GL_PARAM(1,302,dbc); // parameters stored in DB in GL_PRAM table      Int_t sizee = RTtime1.size();
357    if ( parerror<0 ) {      RTbank1.resize(sizee+1);
358      code = parerror;      RTazim.resize(sizee+1);
359      goto closeandexit;      RTerrq.resize(sizee+1);
360    };      RTstart.resize(sizee+1);
361    ltp3 = (Int_t)(glparam2->PATH+glparam2->NAME).Length();      RTpluto1.resize(sizee+1);
362    if ( verbose ) printf(" Reading Earth's Magnetic Field parameter file: %s \n",(glparam2->PATH+glparam2->NAME).Data());      RTbpluto1.resize(sizee+1);
363    //      an>>RTtime1[sizee-1];
364    initize_((char *)niente,&uno,(char *)(glparam->PATH+glparam->NAME).Data(),&ltp2,(char *)(glparam2->PATH+glparam2->NAME).Data(),&ltp3);      an>>RTbank1[sizee-1];
365    //      an>>RTstart[sizee-1];
366    // End IGRF stuff//      an>>RTpluto1[sizee-1];
367    //      an>>RTbpluto1[sizee-1];
368        an>>RTazim[sizee-1];
369        an>>RTerrq[sizee-1];
370        if(sizee>1) {
371          RTtime2.resize(sizee+1);
372          RTbank2.resize(sizee+1);
373          RTpluto2.resize(sizee+1);
374          RTbpluto2.resize(sizee+1);
375          RTtime2[sizee-2]=RTtime1[sizee-1];
376          RTpluto2[sizee-2]=RTpluto1[sizee-1];
377          RTbank2[sizee-2]=RTbank1[sizee-1];
378          RTbpluto2[sizee-2]=RTbpluto1[sizee-1];
379        }
380      }
381      an.close();
382      //cout<<"put some number here"<<endl;
383      //Int_t yupi;
384      //cin>>yupi;
385      
386      if ( verbose ) cout<<"We have read Rotation Table"<<endl;
387        //Geomagnetic coordinates calculations staff
388    
389      GMtype_CoordGeodetic location;
390      //  GMtype_CoordDipole GMlocation;
391      GMtype_Ellipsoid Ellip;
392      GMtype_Data G0, G1, H1;
393            
394      //  { // 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
395      //    TString igpath="/data03/Malakhov/pam9Malakhov/installed10/calib/orb-param/";
396      //  }
397    
398      //  GM_ScanIGRF(glparam->PATH, &G0, &G1, &H1);
399      GM_ScanIGRF(dbc, &G0, &G1, &H1);
400    
401      //cout << G0.element[0] << "\t" << G1.element[0] << "\t" << H1.element[0] << endl;
402      //cout << G0.element[5] << "\t" << G1.element[5] << "\t" << H1.element[5] << endl;
403    
404      GM_SetEllipsoid(&Ellip);
405    
406      // IGRF stuff moved inside run loop!  
407    
408    for (Int_t ip=0;ip<nz;ip++){    for (Int_t ip=0;ip<nz;ip++){
409      zin[ip] = tof->GetZTOF(tof->GetToFPlaneID(ip));      zin[ip] = tof->GetZTOF(tof->GetToFPlaneID(ip));
410    };    };
411    //    //
412    if ( !standalone ){    if ( !standalone ){
413      //      //
414      // Does it contain the Tracker tree?      // Does it contain the Tracker and ToF trees?
415      //      //
416      ttof = (TTree*)file->Get("ToF");      ttof = (TTree*)file->Get("ToF");
417      if ( !ttof ) {      if ( !ttof ) {
418        if ( verbose ) printf(" OrbitalInfo - ERROR: no tof tree\n");        if ( verbose ) printf(" OrbitalInfo - ERROR: no tof tree\n");
419        code = -900;        code = -900;
420        goto closeandexit;        goto closeandexit;
421      };      }
422      ttof->SetBranchAddress("ToFLevel2",&tof);        ttof->SetBranchAddress("ToFLevel2",&tof);  
423      nevtofl2 = ttof->GetEntries();      nevtofl2 = ttof->GetEntries();
424    };  
425        //
426        // Look for extended tracking algorithm
427        //
428        if ( verbose ) printf("Look for extended and nuclei tracking algorithms in ToF\n");
429        // Nuclei tracking algorithm
430        Int_t checkAlgo = 0;
431        tcNucleiTof =  new TClonesArray("ToFTrkVar");
432        checkAlgo = ttof->SetBranchAddress("TrackNuclei",&tcNucleiTof);    
433        if ( !checkAlgo ){
434          if ( verbose ) printf(" Nuclei tracking algorithm ToF branch found! :D \n");
435          hasNucleiTof = true;
436        } else {
437          if ( verbose ) printf(" Nuclei tracking algorithm ToF branch not found :( !\n");
438          printf(" ok, this is not a problem (it depends on tracker settings) \n");
439          delete tcNucleiTof;
440        }
441        // Nuclei tracking algorithm using calorimeter points
442        tcExtNucleiTof =  new TClonesArray("ToFTrkVar");
443        checkAlgo = ttof->SetBranchAddress("RecoveredTrackNuclei",&tcExtNucleiTof);    
444        if ( !checkAlgo ){
445          if ( verbose ) printf(" Recovered nuclei tracking algorithm ToF branch found! :D \n");
446          hasExtNucleiTof = true;
447        } else {
448          if ( verbose ) printf(" Recovered nuclei tracking algorithm ToF branch not found :( !\n");
449          printf(" ok, this is not a problem (it depends on tracker settings) \n");
450          delete tcExtNucleiTof;
451        }
452        // Tracking algorithm using calorimeter points
453        tcExtTof =  new TClonesArray("ToFTrkVar");
454        checkAlgo = ttof->SetBranchAddress("RecoveredTrack",&tcExtTof);
455        if ( !checkAlgo ){
456          if ( verbose ) printf(" Recovered track algorithm ToF branch found! :D \n");
457          hasExtTof = true;
458        } else {
459          if ( verbose ) printf(" Recovered track algorithm ToF branch not found :( !\n");
460          printf(" ok, this is not a problem (it depends on tracker settings) \n");
461          delete tcExtTof;
462        }
463    
464        ttrke = (TTree*)file->Get("Tracker");
465        if ( !ttrke ) {
466          if ( verbose ) printf(" OrbitalInfo - ERROR: no trk tree\n");
467          code = -903;
468          goto closeandexit;
469        }
470        ttrke->SetBranchAddress("TrkLevel2",&trke);  
471        nevtrkl2 = ttrke->GetEntries();
472    
473        //
474        // Look for extended tracking algorithm
475        //
476        if ( verbose ) printf("Look for extended and nuclei tracking algorithms\n");
477        // Nuclei tracking algorithm
478        checkAlgo = 0;
479        tcNucleiTrk =  new TClonesArray("TrkTrack");
480        checkAlgo = ttrke->SetBranchAddress("TrackNuclei",&tcNucleiTrk);    
481        if ( !checkAlgo ){
482          if ( verbose ) printf(" Nuclei tracking algorithm branch found! :D \n");
483          hasNucleiTrk = true;
484        } else {
485          if ( verbose ) printf(" Nuclei tracking algorithm branch not found :( !\n");
486          printf(" ok, this is not a problem (it depends on tracker settings) \n");
487          delete tcNucleiTrk;
488        }
489        // Nuclei tracking algorithm using calorimeter points
490        tcExtNucleiTrk =  new TClonesArray("ExtTrack");
491        checkAlgo = ttrke->SetBranchAddress("RecoveredTrackNuclei",&tcExtNucleiTrk);    
492        if ( !checkAlgo ){
493          if ( verbose ) printf(" Recovered nuclei tracking algorithm branch found! :D \n");
494          hasExtNucleiTrk = true;
495        } else {
496          if ( verbose ) printf(" Recovered nuclei tracking algorithm branch not found :( !\n");
497          printf(" ok, this is not a problem (it depends on tracker settings) \n");
498          delete tcExtNucleiTrk;
499        }
500        // Tracking algorithm using calorimeter points
501        tcExtTrk =  new TClonesArray("ExtTrack");
502        checkAlgo = ttrke->SetBranchAddress("RecoveredTrack",&tcExtTrk);
503        if ( !checkAlgo ){
504          if ( verbose ) printf(" Recovered track algorithm branch found! :D \n");
505          hasExtTrk = true;
506        } else {
507          if ( verbose ) printf(" Recovered track algorithm branch not found :( !\n");
508          printf(" ok, this is not a problem (it depends on tracker settings) \n");
509          delete tcExtTrk;
510        }
511    
512        if ( (hasNucleiTrk && !hasNucleiTof) || (!hasNucleiTrk && hasNucleiTof) ||
513             (hasExtNucleiTrk && !hasExtNucleiTof) || (!hasExtNucleiTrk && hasExtNucleiTof) ||
514             (hasExtTrk && !hasExtTof) || (!hasExtTrk && hasExtTof)
515             ){
516          if ( verbose ) printf(" ERROR: Mismatch between tracker and tof tree branches concerning extended tracking algorithm(s)\n");
517          if ( debug ) printf("hasNucleiTrk %i hasExtNucleiTrk %i hasExtTrk %i \n",hasNucleiTrk,hasExtNucleiTrk,hasExtTrk);
518          if ( debug ) printf("hasNucleiTof %i hasExtNucleiTof %i hasExtTof %i \n",hasNucleiTof,hasExtNucleiTof,hasExtTof);
519          throw -901;
520        }
521    
522      }
523    //    //
524    // Let's start!    // Let's start!
525    //    //
# Line 413  int OrbitalInfoCore(UInt_t run, TFile *f Line 608  int OrbitalInfoCore(UInt_t run, TFile *f
608        //        //
609        reprocall = true;        reprocall = true;
610        //        //
611        if (verbose) printf("\n OrbitalInfo - WARNING: Reprocessing all runs\n");        if (verbose) printf("\n OrbitalInfo - WARNING: Reprocessing all runs\n Deleting old tree...\n");
612        //        //
613      } else {      } else {
614        //        //
# Line 431  int OrbitalInfoCore(UInt_t run, TFile *f Line 626  int OrbitalInfoCore(UInt_t run, TFile *f
626        tempOrbitalInfo = OrbitalInfotrclone->CloneTree(-1,"fast");        tempOrbitalInfo = OrbitalInfotrclone->CloneTree(-1,"fast");
627        tempOrbitalInfo->SetName("OrbitalInfo-old");        tempOrbitalInfo->SetName("OrbitalInfo-old");
628        tempfile->Write();        tempfile->Write();
629          tempOrbitalInfo->Delete();
630        tempfile->Close();          tempfile->Close();  
631      }      }
632      //      //
633      // Delete the old tree from old file and memory      // Delete the old tree from old file and memory
634      //      //
635        OrbitalInfotrclone->Clear();
636      OrbitalInfotrclone->Delete("all");      OrbitalInfotrclone->Delete("all");
637      //      //
638      if (verbose) printf(" ...done!\n");      if (verbose) printf(" ...done!\n");
# Line 450  int OrbitalInfoCore(UInt_t run, TFile *f Line 647  int OrbitalInfoCore(UInt_t run, TFile *f
647    orbitalinfo->Set();//ELENA **TEMPORANEO?**    orbitalinfo->Set();//ELENA **TEMPORANEO?**
648    OrbitalInfotr->Branch("OrbitalInfo","OrbitalInfo",&orbitalinfo);    OrbitalInfotr->Branch("OrbitalInfo","OrbitalInfo",&orbitalinfo);
649    //    //
650      // create new branches for new tracking algorithms
651      //
652      if ( hasNucleiTrk ){
653        torbNucleiTrk = new TClonesArray("OrbitalInfoTrkVar",1);
654        OrbitalInfotr->Branch("TrackNuclei",&torbNucleiTrk);
655      }
656      if ( hasExtNucleiTrk ){
657        torbExtNucleiTrk = new TClonesArray("OrbitalInfoTrkVar",1);
658        OrbitalInfotr->Branch("RecoveredTrackNuclei",&torbExtNucleiTrk);
659      }
660      if ( hasExtTrk ){
661        torbExtTrk = new TClonesArray("OrbitalInfoTrkVar",1);
662        OrbitalInfotr->Branch("RecoveredTrack",&torbExtTrk);
663      }
664    
665      //
666    if ( reproc && !reprocall ){    if ( reproc && !reprocall ){
667      //      //
668      //  open new file and retrieve also tree informations      //  open new file and retrieve also tree informations
# Line 481  int OrbitalInfoCore(UInt_t run, TFile *f Line 694  int OrbitalInfoCore(UInt_t run, TFile *f
694          //          //
695        };        };
696        if (verbose) printf(" Finished successful copying!\n");        if (verbose) printf(" Finished successful copying!\n");
697      };                };
698    };    };
699    //    //
700    //    //
# Line 492  int OrbitalInfoCore(UInt_t run, TFile *f Line 705  int OrbitalInfoCore(UInt_t run, TFile *f
705    // Loop over the run to be processed    // Loop over the run to be processed
706    //    //
707    for (UInt_t irun=0; irun < numbofrun; irun++){    for (UInt_t irun=0; irun < numbofrun; irun++){
708    
709        L_QQ_Q_l_lower = new Quaternions();
710        RYPang_lower = new InclinationInfo();
711        L_QQ_Q_l_upper = new Quaternions();
712        RYPang_upper = new InclinationInfo();
713    
714      //      //
715      // retrieve the first run ID to be processed using the RunInfo list      // retrieve the first run ID to be processed using the RunInfo list
716      //      //
# Line 554  int OrbitalInfoCore(UInt_t run, TFile *f Line 773  int OrbitalInfoCore(UInt_t run, TFile *f
773      //      //
774      //    if ( !totevent ) goto closeandexit;      //    if ( !totevent ) goto closeandexit;
775      // Open Level0 file      // Open Level0 file
776        if ( l0File ) l0File->Close();
777      l0File = new TFile(fname.Data());      l0File = new TFile(fname.Data());
778      if ( !l0File ) {      if ( !l0File ) {
779        if ( debug ) printf(" OrbitalInfo - ERROR: problems opening Level0 file\n");        if ( debug ) printf(" OrbitalInfo - ERROR: problems opening Level0 file\n");
# Line 592  int OrbitalInfoCore(UInt_t run, TFile *f Line 812  int OrbitalInfoCore(UInt_t run, TFile *f
812        code = -12;        code = -12;
813        goto closeandexit;        goto closeandexit;
814      };      };
815    
816        //
817        // open IGRF files and do it only once if we are processing a full level2 file
818        //
819        if ( !igrfloaded ){
820    
821          if ( l0head->GetEntry(runinfo->EV_FROM) > 0 ){
822            igrfloaded = true;
823            //
824            // absolute time of first event of the run (it should not matter a lot)
825            //
826            ph = eh->GetPscuHeader();
827            atime = dbtime->DBabsTime(ph->GetOrbitalTime());
828            
829            parerror=glparam->Query_GL_PARAM(atime-anni5,301,dbc); // parameters stored in DB in GL_PRAM table  
830            if ( parerror<0 ) {
831              code = parerror;
832              goto closeandexit;
833            }
834            ltp1 = (Int_t)(glparam->PATH+glparam->NAME).Length();
835            if ( verbose ) printf(" Reading Earth's Magnetic Field parameter file: %s \n",(glparam->PATH+glparam->NAME).Data());
836            //
837            parerror=glparam2->Query_GL_PARAM(atime,301,dbc); // parameters stored in DB in GL_PRAM table  
838            if ( parerror<0 ) {
839              code = parerror;
840              goto closeandexit;
841            }
842            ltp2 = (Int_t)(glparam2->PATH+glparam2->NAME).Length();
843            if ( verbose ) printf(" Reading Earth's Magnetic Field parameter file: %s \n",(glparam2->PATH+glparam2->NAME).Data());
844            //
845            parerror=glparam3->Query_GL_PARAM(atime,302,dbc); // parameters stored in DB in GL_PRAM table
846            if ( parerror<0 ) {
847              code = parerror;
848              goto closeandexit;
849            }
850            ltp3 = (Int_t)(glparam3->PATH+glparam3->NAME).Length();
851            if ( verbose ) printf(" Reading Earth's Magnetic Field parameter file: %s \n",(glparam3->PATH+glparam3->NAME).Data());
852            //
853            initize_((char *)(glparam->PATH+glparam->NAME).Data(),&ltp1,(char *)(glparam2->PATH+glparam2->NAME).Data(),&ltp2,(char *)(glparam3->PATH+glparam3->NAME).Data(),&ltp3);
854            //
855            if (debug) cout<<"initize: "<<(char *)(glparam->PATH+glparam->NAME).Data()<<"\t"<<(char *)(glparam2->PATH+glparam2->NAME).Data()<<"\t"<<(char *)(glparam3->PATH+glparam3->NAME).Data()<<endl;
856          }
857        }
858      //      //
859      //     TTree *tp = (TTree*)l0File->Get("RunHeader");      // End IGRF stuff//
     //     tp->SetBranchAddress("Header", &eH);  
     //     tp->SetBranchAddress("RunHeader", &reh);  
     //     tp->GetEntry(0);  
     //     ph = eH->GetPscuHeader();  
     //     ULong_t TimeSync = reh->LAST_TIME_SYNC_INFO;  
     //     ULong_t ObtSync = reh->OBT_TIME_SYNC;      
     //     if ( debug ) printf(" 1 TimeSync %lu ObtSync %lu DeltaOBT %lu\n",TimeSync,ObtSync,TimeSync-ObtSync);  
860      //      //
861    
862      ULong_t TimeSync = (ULong_t)dbtime->GetTimesync();      ULong_t TimeSync = (ULong_t)dbtime->GetTimesync();
863      ULong_t ObtSync = (ULong_t)(dbtime->GetObt0()/1000);      ULong_t ObtSync = (ULong_t)(dbtime->GetObt0()/1000);
864      ULong_t DeltaOBT = TimeSync - ObtSync;      ULong_t DeltaOBT = TimeSync - ObtSync;
# Line 688  int OrbitalInfoCore(UInt_t run, TFile *f Line 945  int OrbitalInfoCore(UInt_t run, TFile *f
945        i++;        i++;
946      };      };
947      //      //
     //    l0trm = (TTree*)l0File->Get("Mcmd");  
     //    ch->ls();  
948      ch->SetBranchAddress("Mcmd",&mcmdev);      ch->SetBranchAddress("Mcmd",&mcmdev);
     //    printf(" entries %llu \n", ch->GetEntries());  
     //    l0trm = ch->GetTree();  
     //    neventsm = l0trm->GetEntries();  
949      neventsm = ch->GetEntries();      neventsm = ch->GetEntries();
950      if ( debug ) printf(" entries %u \n", neventsm);      if ( debug ) printf(" entries %u \n", neventsm);
     //    neventsm = 0;  
951      //      //
952      if (neventsm == 0){      if (neventsm == 0){
953        if ( debug ) printf("InclinationInfo - WARNING: No quaternions in this File");        if ( debug ) printf("InclinationInfo - WARNING: No quaternions in this File");
       //      l0File->Close();  
954        code = 900;        code = 900;
       //      goto closeandexit;  
955      }      }
956      //      //
957            Double_t lowerqtime = 0;    
     //    l0trm->SetBranchAddress("Mcmd", &mcmdev);  
     //    l0trm->SetBranchAddress("Header", &eh);  
     //  
     //  
     //  
   
 //    UInt_t mctren = 0;      
 //    UInt_t mcreen = 0;          
     UInt_t numrec = 0;  
     //  
     Double_t upperqtime = 0;  
     Double_t lowerqtime = 0;  
       
 //    Double_t incli = 0;  
 //    oi = 0;  
 //    UInt_t ooi = 0;  
958      //      //
959      // init quaternions information from mcmd-packets      // init quaternions information from mcmd-packets
960      //      //
961      Bool_t isf = true;      Bool_t isf = true;
 //    Int_t fgh = 0;  
962    
963      vector<Float_t> q0;      vector<Float_t> q0;
964      vector<Float_t> q1;      vector<Float_t> q1;
# Line 739  int OrbitalInfoCore(UInt_t run, TFile *f Line 971  int OrbitalInfoCore(UInt_t run, TFile *f
971      vector<Int_t> qmode;      vector<Int_t> qmode;
972    
973      Int_t nt = 0;      Int_t nt = 0;
       
     //init sine-function interpolation  
       
     //cout<<"Sine coeficient initialisation..."<<endl;  
     vector<Sine> q0sine;  
     vector<Sine> q1sine;  
     vector<Sine> q2sine;  
     vector<Sine> q3sine;  
     vector<Sine> Yawsine;  
   
     /*TH2F* q0testing = new TH2F();  
       TH2F* q1testing = new TH2F();  
       TH2F* q2testing = new TH2F();  
       TH2F* q3testing = new TH2F();  
       TH2F* Pitchtesting = new TH2F();  
     */  
974      UInt_t must = 0;      UInt_t must = 0;
975    
976      //      //
# Line 800  int OrbitalInfoCore(UInt_t run, TFile *f Line 1016  int OrbitalInfoCore(UInt_t run, TFile *f
1016            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);
1017            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);
1018            l0File->Close();            l0File->Close();
1019            code = -901;            code = -904;
1020            goto closeandexit;            goto closeandexit;
1021          };          };
1022          //          //
1023          tof->Clear();          tof->Clear();
1024          //          //
1025          if ( ttof->GetEntry(itr) <= 0 ) throw -36;          // Clones array must be cleared before going on
1026            //
1027            if ( hasNucleiTof ){
1028              tcNucleiTof->Delete();
1029            }
1030            if ( hasExtNucleiTof ){
1031              tcExtNucleiTof->Delete();
1032            }          
1033            if ( hasExtTof ){
1034              tcExtTof->Delete();
1035            }
1036            //
1037            if ( verbose ) printf(" get tof tree entries... entry = %i in Level2 file\n",itr);
1038            if ( ttof->GetEntry(itr) <= 0 ){
1039              if ( verbose ) printf(" problems with tof tree entries... entry = %i in Level2 file\n",itr);
1040              if ( verbose ) printf(" nobefrun %u re %u evfrom %u jumped %u reprocall %i \n",nobefrun,re,evfrom,jumped,reprocall);
1041              throw -36;
1042            }
1043            if ( verbose ) printf(" gat0\n");
1044          //          //
1045        };        }
1046          //
1047          // retrieve tracker informations
1048          //
1049          if ( !standalone ){
1050            if ( itr > nevtrkl2 ){  
1051              if ( verbose ) printf(" OrbitalInfo - ERROR: no trk events with entry = %i in Level2 file\n",itr);
1052              if ( debug ) printf(" nobefrun %u re %u evfrom %u jumped %u reprocall %i \n",nobefrun,re,evfrom,jumped,reprocall);
1053              l0File->Close();
1054              code = -905;
1055              goto closeandexit;
1056            }
1057            //
1058            if ( verbose ) printf(" gat1\n");
1059            trke->Clear();
1060            //
1061            // Clones array must be cleared before going on
1062            //
1063            if ( hasNucleiTrk ){
1064              if ( verbose ) printf(" gat2\n");
1065              tcNucleiTrk->Delete();
1066              if ( verbose ) printf(" gat3\n");
1067              torbNucleiTrk->Delete();
1068            }
1069            if ( hasExtNucleiTrk ){
1070              if ( verbose ) printf(" gat4\n");
1071              tcExtNucleiTrk->Delete();
1072              if ( verbose ) printf(" gat5\n");
1073              torbExtNucleiTrk->Delete();
1074            }          
1075            if ( hasExtTrk ){
1076              if ( verbose ) printf(" gat6\n");
1077              tcExtTrk->Delete();
1078              if ( verbose ) printf(" gat7\n");
1079              torbExtTrk->Delete();
1080            }
1081            //
1082            if ( verbose ) printf(" get trk tree entries... entry = %i in Level2 file\n",itr);
1083            if ( ttrke->GetEntry(itr) <= 0 ) throw -36;
1084            //
1085          }
1086    
1087        //        //
1088        procev++;        procev++;
1089        //        //
# Line 820  int OrbitalInfoCore(UInt_t run, TFile *f Line 1095  int OrbitalInfoCore(UInt_t run, TFile *f
1095        OrbitalInfoTrkVar *t_orb = new OrbitalInfoTrkVar();        OrbitalInfoTrkVar *t_orb = new OrbitalInfoTrkVar();
1096        if( !(orbitalinfo->OrbitalInfoTrk) ) orbitalinfo->OrbitalInfoTrk = new TClonesArray("OrbitalInfoTrkVar",2);        if( !(orbitalinfo->OrbitalInfoTrk) ) orbitalinfo->OrbitalInfoTrk = new TClonesArray("OrbitalInfoTrkVar",2);
1097        TClonesArray &tor = *orbitalinfo->OrbitalInfoTrk;        TClonesArray &tor = *orbitalinfo->OrbitalInfoTrk;
1098    
1099          // Geomagnetic coordinates calculation variables
1100          GMtype_CoordSpherical CoordSpherical, DipoleSpherical;
1101          GMtype_CoordCartesian CoordCartesian, DipoleCartesian;
1102          GMtype_Model Model;
1103          GMtype_Pole Pole;
1104    
1105        //        //
1106        // Fill OBT, pkt_num and absTime        // Fill OBT, pkt_num and absTime
1107        //              //      
# Line 849  int OrbitalInfoCore(UInt_t run, TFile *f Line 1131  int OrbitalInfoCore(UInt_t run, TFile *f
1131              + (month*31.+ (float) day)/365.              + (month*31.+ (float) day)/365.
1132              + (hour*3600.+min*60.+(float)sec)/(24.*3600.*365.);              + (hour*3600.+min*60.+(float)sec)/(24.*3600.*365.);
1133            //            //
1134            if ( debug ) printf(" %i compute magnetic dipole moment get dipole moment for year\n",procev);                  if ( debug ) printf(" %i compute magnetic dipole moment get dipole moment for year\n",procev);            
1135              if ( debug ) printf(" %i jyear %f dimo %f \n",procev,jyear,dimo);
1136            feldcof_(&jyear, &dimo); // get dipole moment for year            feldcof_(&jyear, &dimo); // get dipole moment for year
1137            if ( debug ) printf(" %i compute magnetic dipole moment end\n",procev);                  if ( debug ) printf(" %i compute magnetic dipole moment end\n",procev);
1138    
1139              GM_TimeAdjustCoefs(year, jyear, G0, G1, H1, &Model);
1140              GM_PoleLocation(Model, &Pole);
1141              
1142          } else {          } else {
1143            code = -56;            code = -56;
1144            goto closeandexit;            goto closeandexit;
# Line 861  int OrbitalInfoCore(UInt_t run, TFile *f Line 1148  int OrbitalInfoCore(UInt_t run, TFile *f
1148        //        //
1149        cOrbit orbits(*gltle->GetTle());        cOrbit orbits(*gltle->GetTle());
1150        //        //
       if ( debug ) printf(" I am Here \n");  
       //  
1151        // synchronize with quaternions data        // synchronize with quaternions data
1152        //        //
1153        if ( isf && neventsm>0 ){        if ( isf && neventsm>0 ){
# Line 870  int OrbitalInfoCore(UInt_t run, TFile *f Line 1155  int OrbitalInfoCore(UInt_t run, TFile *f
1155          // First event          // First event
1156          //          //
1157          isf = false;          isf = false;
1158          upperqtime = atime;          //      upperqtime = atime;
1159          lowerqtime = runinfo->RUNHEADER_TIME;          lowerqtime = runinfo->RUNHEADER_TIME;
1160          for ( ik = 0; ik < neventsm; ik++){  //number of macrocommad packets          for ( ik = 0; ik < neventsm; ik++){  //number of macrocommad packets
1161            if ( ch->GetEntry(ik) <= 0 ) throw -36;            if ( ch->GetEntry(ik) <= 0 ) throw -36;
1162            tmpSize = mcmdev->Records->GetEntries();            tmpSize = mcmdev->Records->GetEntries();
1163            numrec = tmpSize;            //      numrec = tmpSize;
1164              if ( debug ) cout << "packet number " << ik <<"\tnumber of subpackets is " << tmpSize << endl;
1165            for (Int_t j3 = 0;j3<tmpSize;j3++){  //number of subpackets            for (Int_t j3 = 0;j3<tmpSize;j3++){  //number of subpackets
             if ( debug ) printf(" ik %i j3 %i eh eh eh \n",ik,j3);  
1166              mcmdrc = (pamela::McmdRecord*)mcmdev->Records->At(j3);              mcmdrc = (pamela::McmdRecord*)mcmdev->Records->At(j3);
1167              if ( mcmdrc ){ // missing inclination bug [8RED 090116]              if ( mcmdrc ){ // missing inclination bug [8RED 090116]
1168                if ( debug ) printf(" pluto \n");                if ( debug ) printf(" pluto \n");
1169                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
1170                  L_QQ_Q_l_upper->fill(mcmdrc->McmdData);                 L_QQ_Q_l_upper->fill(mcmdrc->McmdData);
1171                  for (UInt_t ui = 0; ui < 6; ui++){                  for (UInt_t ui = 0; ui < 6; ui++){
1172                    if (ui>0){                    if (ui>0){
1173                      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]){
1174                          if ( debug ) printf(" here1 %i \n",ui);                        if ( debug ) printf(" here1 %i \n",ui);
1175                        Double_t u_time = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[ui]*1000-DeltaOBT*1000));                        Double_t u_time = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[ui]*1000-DeltaOBT*1000));
1176                        Int_t recSize = recqtime.size();                        Int_t recSize = recqtime.size();
1177                        for(Int_t mu = nt;mu<recSize;mu++){                        if(lowerqtime > recqtime[recSize-1]){
1178                          if(recqtime[mu]>lowerqtime && recqtime[mu]<u_time){                           // to avoid interpolation between bad quaternions arrays
1179                            nt=mu;                           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();  
                           inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);  
                           qtime[sizeqmcmd]=recqtime[mu];  
                           q0[sizeqmcmd]=recq0[mu];  
                           q1[sizeqmcmd]=recq1[mu];  
                           q2[sizeqmcmd]=recq2[mu];  
                           q3[sizeqmcmd]=recq3[mu];  
                           qmode[sizeqmcmd]=-10;  
                           orbits.getPosition((double) (qtime[sizeqmcmd] - gltle->GetFromTime())/60., &eCi);  
                           RYPang_upper->TransAngle(eCi.getPos().m_x,eCi.getPos().m_y,eCi.getPos().m_z,eCi.getVel().m_x,eCi.getVel().m_y,eCi.getVel().m_z,recq0[mu],recq1[mu],recq2[mu],recq3[mu]);  
                           qRoll[sizeqmcmd]=RYPang_upper->Kren;  
                           qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;  
                           qPitch[sizeqmcmd]=RYPang_upper->Tangazh;  
                         }  
                         if(recqtime[mu]>=u_time){  
1180                            Int_t sizeqmcmd = qtime.size();                            Int_t sizeqmcmd = qtime.size();
1181                            inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);                            inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);
1182                            qtime[sizeqmcmd]=u_time;                            qtime[sizeqmcmd]=u_time;
# Line 921  int OrbitalInfoCore(UInt_t run, TFile *f Line 1191  int OrbitalInfoCore(UInt_t run, TFile *f
1191                            qRoll[sizeqmcmd]=RYPang_upper->Kren;                            qRoll[sizeqmcmd]=RYPang_upper->Kren;
1192                            qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;                            qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;
1193                            qPitch[sizeqmcmd]=RYPang_upper->Tangazh;                            qPitch[sizeqmcmd]=RYPang_upper->Tangazh;
1194                            break;                           }
1195                          }
1196                          for(Int_t mu = nt;mu<recSize;mu++){
1197                            if(recqtime[mu]>lowerqtime && recqtime[mu]<u_time){
1198                              if(sqrt(pow(recq0[mu],2)+pow(recq1[mu],2)+pow(recq2[mu],2)+pow(recq3[mu],2))>0.99999){
1199                                nt=mu;
1200                                Int_t sizeqmcmd = qtime.size();
1201                                inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);
1202                                qtime[sizeqmcmd]=recqtime[mu];
1203                                q0[sizeqmcmd]=recq0[mu];
1204                                q1[sizeqmcmd]=recq1[mu];
1205                                q2[sizeqmcmd]=recq2[mu];
1206                                q3[sizeqmcmd]=recq3[mu];
1207                                qmode[sizeqmcmd]=-10;
1208                                orbits.getPosition((double) (qtime[sizeqmcmd] - gltle->GetFromTime())/60., &eCi);
1209                                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]);
1210                                qRoll[sizeqmcmd]=RYPang_upper->Kren;
1211                                qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;
1212                                qPitch[sizeqmcmd]=RYPang_upper->Tangazh;
1213                              }
1214                            }
1215                            if(recqtime[mu]>=u_time){
1216                              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){
1217                                Int_t sizeqmcmd = qtime.size();
1218                                inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);
1219                                qtime[sizeqmcmd]=u_time;
1220                                q0[sizeqmcmd]=L_QQ_Q_l_upper->quat[ui][0];
1221                                q1[sizeqmcmd]=L_QQ_Q_l_upper->quat[ui][1];
1222                                q2[sizeqmcmd]=L_QQ_Q_l_upper->quat[ui][2];
1223                                q3[sizeqmcmd]=L_QQ_Q_l_upper->quat[ui][3];
1224                                qmode[sizeqmcmd]=holeq(lowerqtime,qtime[sizeqmcmd],L_QQ_Q_l_lower,L_QQ_Q_l_upper,ui);
1225                                lowerqtime = u_time;
1226                                orbits.getPosition((double) (u_time - gltle->GetFromTime())/60., &eCi);
1227                                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]);
1228                                qRoll[sizeqmcmd]=RYPang_upper->Kren;
1229                                qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;
1230                                qPitch[sizeqmcmd]=RYPang_upper->Tangazh;
1231                                break;
1232                              }
1233                          }                          }
1234                        }                        }
1235                      }                      }
1236                    }else{                    }else{
1237                          if ( debug ) printf(" here2 %i \n",ui);                      if ( debug ) printf(" here2 %i \n",ui);
1238                      Double_t u_time = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));                      Double_t u_time = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
1239                      if(lowerqtime>u_time)nt=0;                      if(lowerqtime>u_time)nt=0;
1240                      Int_t recSize = recqtime.size();                      Int_t recSize = recqtime.size();
1241                      for(Int_t mu = nt;mu<recSize;mu++){                      if(lowerqtime > recqtime[recSize-1]){
1242                        if(recqtime[mu]>lowerqtime && 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){
                         nt=mu;  
                         Int_t sizeqmcmd = qtime.size();  
                         inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);  
                         qtime[sizeqmcmd]=recqtime[mu];  
                         q0[sizeqmcmd]=recq0[mu];  
                         q1[sizeqmcmd]=recq1[mu];  
                         q2[sizeqmcmd]=recq2[mu];  
                         q3[sizeqmcmd]=recq3[mu];  
                         qmode[sizeqmcmd]=-10;  
                         orbits.getPosition((double) (qtime[sizeqmcmd] - gltle->GetFromTime())/60., &eCi);  
                         RYPang_upper->TransAngle(eCi.getPos().m_x,eCi.getPos().m_y,eCi.getPos().m_z,eCi.getVel().m_x,eCi.getVel().m_y,eCi.getVel().m_z,recq0[mu],recq1[mu],recq2[mu],recq3[mu]);  
                         qRoll[sizeqmcmd]=RYPang_upper->Kren;  
                         qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;  
                         qPitch[sizeqmcmd]=RYPang_upper->Tangazh;  
                       }  
                       if(recqtime[mu]>=u_time){  
1243                          Int_t sizeqmcmd = qtime.size();                          Int_t sizeqmcmd = qtime.size();
1244                          inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);                          inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);
1245                          qtime[sizeqmcmd]=u_time;                          qtime[sizeqmcmd]=u_time;
# Line 962  int OrbitalInfoCore(UInt_t run, TFile *f Line 1254  int OrbitalInfoCore(UInt_t run, TFile *f
1254                          qRoll[sizeqmcmd]=RYPang_upper->Kren;                          qRoll[sizeqmcmd]=RYPang_upper->Kren;
1255                          qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;                          qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;
1256                          qPitch[sizeqmcmd]=RYPang_upper->Tangazh;                          qPitch[sizeqmcmd]=RYPang_upper->Tangazh;
1257                          CopyQ(L_QQ_Q_l_lower,L_QQ_Q_l_upper);                        }
1258                          break;                      }
1259                        for(Int_t mu = nt;mu<recSize;mu++){
1260                          if(recqtime[mu]>lowerqtime && recqtime[mu]<u_time){
1261                             if(sqrt(pow(recq0[mu],2)+pow(recq1[mu],2)+pow(recq2[mu],2)+pow(recq3[mu],2))>0.99999){
1262                               nt=mu;
1263                               Int_t sizeqmcmd = qtime.size();
1264                               inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);
1265                               qtime[sizeqmcmd]=recqtime[mu];
1266                               q0[sizeqmcmd]=recq0[mu];
1267                               q1[sizeqmcmd]=recq1[mu];
1268                               q2[sizeqmcmd]=recq2[mu];
1269                               q3[sizeqmcmd]=recq3[mu];
1270                               qmode[sizeqmcmd]=-10;
1271                               orbits.getPosition((double) (qtime[sizeqmcmd] - gltle->GetFromTime())/60., &eCi);
1272                               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]);
1273                               qRoll[sizeqmcmd]=RYPang_upper->Kren;
1274                               qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;
1275                               qPitch[sizeqmcmd]=RYPang_upper->Tangazh;
1276                             }
1277                          }
1278                          if(recqtime[mu]>=u_time){
1279                             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){
1280                               Int_t sizeqmcmd = qtime.size();
1281                               inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);
1282                               qtime[sizeqmcmd]=u_time;
1283                               q0[sizeqmcmd]=L_QQ_Q_l_upper->quat[0][0];
1284                               q1[sizeqmcmd]=L_QQ_Q_l_upper->quat[0][1];
1285                               q2[sizeqmcmd]=L_QQ_Q_l_upper->quat[0][2];
1286                               q3[sizeqmcmd]=L_QQ_Q_l_upper->quat[0][3];
1287                               qmode[sizeqmcmd]=holeq(lowerqtime,qtime[sizeqmcmd],L_QQ_Q_l_lower,L_QQ_Q_l_upper,ui);
1288                               lowerqtime = u_time;
1289                               orbits.getPosition((double) (u_time - gltle->GetFromTime())/60., &eCi);
1290                               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]);
1291                               qRoll[sizeqmcmd]=RYPang_upper->Kren;
1292                               qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;
1293                               qPitch[sizeqmcmd]=RYPang_upper->Tangazh;
1294                               CopyQ(L_QQ_Q_l_lower,L_QQ_Q_l_upper);
1295                               break;
1296                             }
1297                        }                        }
1298                      }                      }
1299                    }                    }
1300                  }                  }
1301                }                }
1302              }              }
1303              if ( debug ) printf(" ciccio \n");              //if ( debug ) cout << "subpacket " << j3 << "\t qtime = " << qtime[qtime.size()-1] << endl;
1304            }            }
1305          }          }
1306                    
1307          if(qtime.size()==0){          if(qtime.size()==0){                            // in case if no orientation information in data
1308              for(Int_t my=0;my<recqtime.size();my++){            if ( debug ) cout << "qtime.size() = 0" << endl;
1309                  Int_t sizeqmcmd = qtime.size();            for(UInt_t my=0;my<recqtime.size();my++){
1310                  inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);              if(sqrt(pow(recq0[my],2)+pow(recq1[my],2)+pow(recq2[my],2)+pow(recq3[my],2))>0.99999){
1311                  qtime[sizeqmcmd]=recqtime[my];                Int_t sizeqmcmd = qtime.size();
1312                  q0[sizeqmcmd]=recq0[my];                inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);
1313                  q1[sizeqmcmd]=recq1[my];                qtime[sizeqmcmd]=recqtime[my];
1314                  q2[sizeqmcmd]=recq2[my];                q0[sizeqmcmd]=recq0[my];
1315                  q3[sizeqmcmd]=recq3[my];                q1[sizeqmcmd]=recq1[my];
1316                  qmode[sizeqmcmd]=-10;                q2[sizeqmcmd]=recq2[my];
1317                  orbits.getPosition((double) (qtime[sizeqmcmd] - gltle->GetFromTime())/60., &eCi);                q3[sizeqmcmd]=recq3[my];
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,recq0[my],recq1[my],recq2[my],recq3[my]);                qmode[sizeqmcmd]=-10;
1319                  qRoll[sizeqmcmd]=RYPang_upper->Kren;                orbits.getPosition((double) (qtime[sizeqmcmd] - gltle->GetFromTime())/60., &eCi);
1320                  qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;                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]);
1321                  qPitch[sizeqmcmd]=RYPang_upper->Tangazh;                qRoll[sizeqmcmd]=RYPang_upper->Kren;
1322              }                qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;
1323                  qPitch[sizeqmcmd]=RYPang_upper->Tangazh;
1324                }
1325              }
1326          }          }
1327                    
1328          if ( debug ) printf(" fuffi \n");  
         sineparam(q0sine,qtime,q0,qRoll,qPitch,0.60);  
         sineparam(q1sine,qtime,q1,qRoll,qPitch,0.82);  
         sineparam(q2sine,qtime,q2,qRoll,qPitch,0.82);  
         sineparam(q3sine,qtime,q3,qRoll,qPitch,0.60);  
         sineparam(Yawsine,qtime,qYaw,qRoll,qPitch,4);  
1329          if ( debug ) printf(" puffi \n");          if ( debug ) printf(" puffi \n");
1330          Double_t tmin = 9999999999.;          Double_t tmin = 9999999999.;
1331          Double_t tmax = 0.;          Double_t tmax = 0.;
# Line 1005  int OrbitalInfoCore(UInt_t run, TFile *f Line 1333  int OrbitalInfoCore(UInt_t run, TFile *f
1333            if(qtime[tre]>tmax)tmax = qtime[tre];            if(qtime[tre]>tmax)tmax = qtime[tre];
1334            if(qtime[tre]<tmin)tmin = qtime[tre];            if(qtime[tre]<tmin)tmin = qtime[tre];
1335          }          }
1336          if ( debug ) printf(" gnfuffi \n");          // sorting quaternions by time
1337           Bool_t t = true;
1338            while(t){
1339             t=false;
1340              for(UInt_t i=0;i<qtime.size()-1;i++){
1341                if(qtime[i]>qtime[i+1]){
1342                  Double_t tmpr = qtime[i];
1343                  qtime[i]=qtime[i+1];
1344                  qtime[i+1] = tmpr;
1345                  tmpr = q0[i];
1346                  q0[i]=q0[i+1];
1347                  q0[i+1] = tmpr;
1348                  tmpr = q1[i];
1349                  q1[i]=q1[i+1];
1350                  q1[i+1] = tmpr;
1351                  tmpr = q2[i];
1352                  q2[i]=q2[i+1];
1353                  q2[i+1] = tmpr;
1354                  tmpr = q3[i];
1355                  q3[i]=q3[i+1];
1356                  q3[i+1] = tmpr;
1357                  tmpr = qRoll[i];
1358                  qRoll[i]=qRoll[i+1];
1359                  qRoll[i+1] = tmpr;
1360                  tmpr = qYaw[i];
1361                  qYaw[i]=qYaw[i+1];
1362                  qYaw[i+1] = tmpr;
1363                  tmpr = qPitch[i];
1364                  qPitch[i]=qPitch[i+1];
1365                  qPitch[i+1] = tmpr;
1366                    t=true;
1367                }
1368              }
1369            }
1370    
1371            if ( debug ){
1372              cout << "we have loaded quaternions: size of quaternions set is "<< qtime.size() << endl;
1373             for(UInt_t i=0;i<qtime.size();i++) cout << qtime[i] << "\t";
1374              cout << endl << endl;
1375              Int_t lopu;
1376              cin >> lopu;
1377           }
1378    
         //q0testing->SetName("q0testing");  
         //q1testing->SetName("q1testing");  
         //q2testing->SetName("q2testing");  
         //q3testing->SetName("q3testing");  
           
 //      Int_t ss=10.*(tmax-tmin);  
         //q0testing->SetBins(ss,tmin,tmax,1000,-1.,1.);  
         //Pitchtesting->SetBins(ss,tmin,tmax,1000,-40.,40.);  
   
 //      for(Int_t tre = 0;tre<qtime.size();tre++){  
           //cout<<"q0["<<tre<<" = "<<q0[tre]<<endl;  
           //q0testing->Fill(qtime[tre],q0[tre]);  
           //q1testing->Fill(qtime[tre],q1[tre]);  
           //Pitchtesting->Fill(qtime[tre],qPitch[tre],100);  
           //if(qmode[tre] == -10)Pitchtesting->Fill(qtime[tre],10,100);  
           //q2testing->Fill(qtime[tre],q2[tre],100);  
           //q3testing->Fill(qtime[tre],q3[tre],100);  
 //      }  
           
         //for(Int_t tre=0;tre<q0sine.size();tre++)cout<<q1sine[tre].A<<"*sin("<<q1sine[tre].b<<"x+"<<q1sine[tre].c<<")\t time start: "<<q1sine[tre].startPoint<<"\ttime end: "<<q1sine[tre].finishPoint<<endl;  
         //for(Int_t tre=0;tre<q0sine.size();tre++)cout<<q1sine[tre].A<<"*sin("<<q1sine[tre].b<<"x+"<<q1sine[tre].c<<")\t time start: "<<q0sine[tre].startPoint<<"\ttime end: "<<q0sine[tre].finishPoint<<endl;  
1379        } // if we processed first event        } // if we processed first event
1380    
1381                
1382        //Filling Inclination information        //Filling Inclination information
1383        Double_t incli = 0;        Double_t incli = 0;
1384        if ( qtime.size() > 1 ){        if ( qtime.size() > 1 ){
1385        for(UInt_t mu = must;mu<qtime.size()-1;mu++){          if ( debug ) cout << "ok quaternions is exist and mu = " << must << endl;
1386          if ( debug ) printf(" ??grfuffi %i sixe %i must %i \n",mu,qtime.size()-1,must);          if ( debug ) cout << "qtimes[ " << qtime[0] << " , " << qtime[qtime.size()-1] << " ]\tatime = "<<atime<<endl;
1387          if(qtime[mu+1]>qtime[mu]){          for(UInt_t mu = must;mu<qtime.size()-1;mu++){
1388            if ( debug ) printf(" grfuffi2 %i \n",mu);            if ( debug ) printf(" ??grfuffi %i sixe %i must %i \n",mu,qtime.size()-1,must);
1389            if(atime<=qtime[mu+1] && atime>=qtime[mu]){            if(qtime[mu+1]>qtime[mu]){
1390              must = mu;              if ( debug ) cout << "qtime[" << mu << "] = " << qtime[mu] << "\tqtime[" << mu+1 << "] = " << qtime[mu+1] << "\tatime = " << atime << endl;
1391              incli = (qPitch[mu+1]-qPitch[mu])/(qtime[mu+1]-qtime[mu]);              if(atime<=qtime[mu+1] && atime>=qtime[mu]){
1392              orbitalinfo->theta =  incli*atime+qPitch[mu+1]-incli*qtime[mu+1];                if ( debug ) cout << "here we have found proper quaternions for interpolation: mu = "<<mu<<endl;
1393              incli = (qRoll[mu+1]-qRoll[mu])/(qtime[mu+1]-qtime[mu]);                must = mu;
1394              orbitalinfo->etha =  incli*atime+qRoll[mu+1]-incli*qtime[mu+1];                incli = (qPitch[mu+1]-qPitch[mu])/(qtime[mu+1]-qtime[mu]);
1395              incli = (qYaw[mu+1]-qYaw[mu])/(qtime[mu+1]-qtime[mu]);                orbitalinfo->theta =  incli*atime+qPitch[mu+1]-incli*qtime[mu+1];
1396              orbitalinfo->phi =  incli*atime+qYaw[mu+1]-incli*qtime[mu+1];                incli = (qRoll[mu+1]-qRoll[mu])/(qtime[mu+1]-qtime[mu]);
1397                              orbitalinfo->etha =  incli*atime+qRoll[mu+1]-incli*qtime[mu+1];
1398              incli = (q0[mu+1]-q0[mu])/(qtime[mu+1]-qtime[mu]);                incli = (qYaw[mu+1]-qYaw[mu])/(qtime[mu+1]-qtime[mu]);
1399              orbitalinfo->q0t =  incli*atime+q0[mu+1]-incli*qtime[mu+1];                orbitalinfo->phi =  incli*atime+qYaw[mu+1]-incli*qtime[mu+1];
1400              incli = (q1[mu+1]-q1[mu])/(qtime[mu+1]-qtime[mu]);                
1401              orbitalinfo->q1t =  incli*atime+q1[mu+1]-incli*qtime[mu+1];                incli = (q0[mu+1]-q0[mu])/(qtime[mu+1]-qtime[mu]);
1402              incli = (q2[mu+1]-q2[mu])/(qtime[mu+1]-qtime[mu]);                orbitalinfo->q0 =  incli*atime+q0[mu+1]-incli*qtime[mu+1];
1403              orbitalinfo->q2t =  incli*atime+q2[mu+1]-incli*qtime[mu+1];                incli = (q1[mu+1]-q1[mu])/(qtime[mu+1]-qtime[mu]);
1404              incli = (q3[mu+1]-q3[mu])/(qtime[mu+1]-qtime[mu]);                orbitalinfo->q1 =  incli*atime+q1[mu+1]-incli*qtime[mu+1];
1405              orbitalinfo->q3t =  incli*atime+q3[mu+1]-incli*qtime[mu+1];                incli = (q2[mu+1]-q2[mu])/(qtime[mu+1]-qtime[mu]);
1406                              orbitalinfo->q2 =  incli*atime+q2[mu+1]-incli*qtime[mu+1];
1407              orbitalinfo->TimeGap = qtime[mu+1]-qtime[mu];                incli = (q3[mu+1]-q3[mu])/(qtime[mu+1]-qtime[mu]);
1408              orbitalinfo->mode = qmode[mu+1];                orbitalinfo->q3 =  incli*atime+q3[mu+1]-incli*qtime[mu+1];
1409              if(qmode[mu+1]==-10) orbitalinfo->R10r = true;else orbitalinfo->R10r = false;                Float_t tg = (qtime[mu+1]-qtime[mu])/1000.0;
1410              if(qmode[mu+1]==-10 || qmode[mu+1]==0 || qmode[mu+1]==1 || qmode[mu+1]==3 || qmode[mu+1]==4 || qmode[mu+1]==6){                if(tg>=1) tg=0.00;
1411                //linear interpolation                orbitalinfo->TimeGap = TMath::Min(TMath::Abs(qtime[mu+1])-atime,TMath::Abs(atime-qtime[mu]))+tg;//qtime[mu+1]-qtime[mu];
1412                incli = (q0[mu+1]-q0[mu])/(qtime[mu+1]-qtime[mu]);                orbitalinfo->mode = qmode[mu+1];
1413                orbitalinfo->q0 =  incli*atime+q0[mu+1]-incli*qtime[mu+1];                //if(atime==qtime[mu] || atime==qtime[mu+1]) orbitalinfo->qkind = 0; else orbitalinfo->qkind=1;
1414                incli = (q1[mu+1]-q1[mu])/(qtime[mu+1]-qtime[mu]);                //if(qmode[mu+1]==-10) orbitalinfo->R10r = true;else orbitalinfo->R10r = false;
1415                orbitalinfo->q1 =  incli*atime+q1[mu+1]-incli*qtime[mu+1];                if ( debug ) printf(" grfuffi4 %i \n",mu);
1416                incli = (q2[mu+1]-q2[mu])/(qtime[mu+1]-qtime[mu]);                break;
1417                orbitalinfo->q2 =  incli*atime+q2[mu+1]-incli*qtime[mu+1];              }
1418                incli = (q3[mu+1]-q3[mu])/(qtime[mu+1]-qtime[mu]);            }
1419                orbitalinfo->q3 =  incli*atime+q3[mu+1]-incli*qtime[mu+1];          }
             }else{  
               //sine interpolation  
               for(UInt_t mt=0;mt<q0sine.size();mt++){  
                 if(atime<=q0sine[mt].finishPoint && atime>=q0sine[mt].startPoint){  
                   if(!q0sine[mt].NeedFit)orbitalinfo->q0=q0sine[mt].A*sin(q0sine[mt].b*atime+q0sine[mt].c);else{  
                     incli = (q0[mu+1]-q0[mu])/(qtime[mu+1]-qtime[mu]);  
                     orbitalinfo->q0 =  incli*atime+q0[mu+1]-incli*qtime[mu+1];  
                   }  
                 }  
                 if(atime<=q1sine[mt].finishPoint && atime>=q1sine[mt].startPoint){  
                   if(!q1sine[mt].NeedFit)orbitalinfo->q1=q1sine[mt].A*sin(q1sine[mt].b*atime+q1sine[mt].c);else{  
                     incli = (q1[mu+1]-q1[mu])/(qtime[mu+1]-qtime[mu]);  
                     orbitalinfo->q1 =  incli*atime+q1[mu+1]-incli*qtime[mu+1];  
                   }  
                 }  
                 if(atime<=q2sine[mt].finishPoint && atime>=q2sine[mt].startPoint){  
                   if(!q2sine[mt].NeedFit)orbitalinfo->q2=q0sine[mt].A*sin(q2sine[mt].b*atime+q2sine[mt].c);else{  
                     incli = (q2[mu+1]-q2[mu])/(qtime[mu+1]-qtime[mu]);  
                     orbitalinfo->q2 =  incli*atime+q2[mu+1]-incli*qtime[mu+1];  
                   }  
                 }  
                 if(atime<=q3sine[mt].finishPoint && atime>=q3sine[mt].startPoint){  
                   if(!q3sine[mt].NeedFit)orbitalinfo->q3=q0sine[mt].A*sin(q3sine[mt].b*atime+q3sine[mt].c);else{  
                     incli = (q3[mu+1]-q3[mu])/(qtime[mu+1]-qtime[mu]);  
                     orbitalinfo->q3 =  incli*atime+q3[mu+1]-incli*qtime[mu+1];  
                   }  
                 }  
                 if(atime<=Yawsine[mt].finishPoint && atime>=Yawsine[mt].startPoint){  
                   if(!Yawsine[mt].NeedFit)orbitalinfo->phi=Yawsine[mt].A*sin(Yawsine[mt].b*atime+Yawsine[mt].c);else{  
                     incli = (qYaw[mu+1]-qYaw[mu])/(qtime[mu+1]-qtime[mu]);  
                     orbitalinfo->phi =  incli*atime+qYaw[mu+1]-incli*qtime[mu+1];  
                   }  
                 }  
               }  
             }  
             //q0testing->Fill(atime,orbitalinfo->q0,100);  
             //q1testing->Fill(atime,orbitalinfo->q1,100);  
             //Pitchtesting->Fill(atime,orbitalinfo->etha);  
             //q2testing->Fill(atime,orbitalinfo->q2);  
             //q3testing->Fill(atime,orbitalinfo->q3);  
             break;  
           }  
         }  
       }  
1420        }        }
1421          if ( debug ) printf(" grfuffi5  \n");
1422        //        //
1423        // ops no inclination information        // ops no inclination information
1424        //        //
1425          
1426        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 ){
1427            if ( debug ) cout << "ops no iclination information" << endl;
1428          orbitalinfo->mode = 10;          orbitalinfo->mode = 10;
1429          orbitalinfo->q0 = -1000.;          orbitalinfo->q0 = -1000.;
1430          orbitalinfo->q1 = -1000.;          orbitalinfo->q1 = -1000.;
# Line 1126  int OrbitalInfoCore(UInt_t run, TFile *f Line 1433  int OrbitalInfoCore(UInt_t run, TFile *f
1433          orbitalinfo->etha = -1000.;          orbitalinfo->etha = -1000.;
1434          orbitalinfo->phi = -1000.;          orbitalinfo->phi = -1000.;
1435          orbitalinfo->theta = -1000.;          orbitalinfo->theta = -1000.;
1436        };          orbitalinfo->TimeGap = -1000.;
1437            //orbitalinfo->qkind = -1000;
1438            
1439            //      if ( debug ){
1440            //        Int_t lopu;
1441            //         cin >> lopu;
1442            //      }
1443            if ( debug ) printf(" grfuffi6 \n");
1444          }
1445        //        //
1446          if ( debug ) printf(" filling \n");
1447        // #########################################################################################################################          // #########################################################################################################################  
1448        //        //
1449        // fill orbital positions        // fill orbital positions
# Line 1138  int OrbitalInfoCore(UInt_t run, TFile *f Line 1454  int OrbitalInfoCore(UInt_t run, TFile *f
1454        lon = (coo.m_Lon > M_PI) ? rad2deg(coo.m_Lon - 2*M_PI) : rad2deg(coo.m_Lon);        lon = (coo.m_Lon > M_PI) ? rad2deg(coo.m_Lon - 2*M_PI) : rad2deg(coo.m_Lon);
1455        lat = rad2deg(coo.m_Lat);        lat = rad2deg(coo.m_Lat);
1456        alt = coo.m_Alt;        alt = coo.m_Alt;
1457        //  
1458          cOrbit orbits2(*gltle->GetTle());
1459          orbits2.getPosition((double) (atime - gltle->GetFromTime())/60., &eCi);
1460          //      Float_t x=eCi.getPos().m_x;
1461          //      Float_t y=eCi.getPos().m_y;
1462          //      Float_t z=eCi.getPos().m_z;
1463          
1464          TVector3 V(eCi.getVel().m_x,eCi.getVel().m_y,eCi.getVel().m_z);
1465          TVector3 Pos(eCi.getPos().m_x,eCi.getPos().m_y,eCi.getPos().m_z);
1466          
1467          Float_t dlon=Pos.Phi()*TMath::RadToDeg()-lon;
1468          
1469          Pos.RotateZ(-dlon*TMath::DegToRad());
1470          V.RotateZ(-dlon*TMath::DegToRad());
1471          Float_t diro;
1472          if(V.Z()>0) diro=1; else diro=-1;
1473          
1474          // 10REDNEW
1475          Int_t errq=0;
1476          Int_t azim=0;
1477          Int_t MU=0;
1478          for(UInt_t mu = must;mu<RTtime2.size()-1;mu++){
1479            if(atime<=RTstart[mu+1] && atime>=RTstart[mu]){
1480              errq=RTerrq[mu];
1481              azim=RTazim[mu];
1482              MU=mu;
1483              break;
1484            }
1485          }
1486          orbitalinfo->errq = errq;
1487          orbitalinfo->azim = azim;
1488          orbitalinfo->qkind = 0;
1489          
1490          if ( debug ) printf(" coord done \n");
1491        if( lon<180 && lon>-180 && lat<90 && lat>-90 && alt>0 ){          if( lon<180 && lon>-180 && lat<90 && lat>-90 && alt>0 ){  
1492          //                //      
1493          orbitalinfo->lon = lon;          orbitalinfo->lon = lon;
1494          orbitalinfo->lat = lat;          orbitalinfo->lat = lat;
1495          orbitalinfo->alt = alt ;          orbitalinfo->alt = alt;
1496            orbitalinfo->V = V;
1497    
1498            //      GMtype_CoordGeodetic  location;
1499            location.lambda = lon;
1500            location.phi = lat;
1501            location.HeightAboveEllipsoid = alt;
1502    
1503            GM_GeodeticToSpherical(Ellip, location, &CoordSpherical);
1504            GM_SphericalToCartesian(CoordSpherical,  &CoordCartesian);
1505            GM_EarthCartToDipoleCartCD(Pole, CoordCartesian, &DipoleCartesian);
1506            GM_CartesianToSpherical(DipoleCartesian, &DipoleSpherical);
1507            orbitalinfo->londip = DipoleSpherical.lambda;
1508            orbitalinfo->latdip = DipoleSpherical.phig;
1509    
1510            if(debug)cout<<"geodetic:\t"<<lon<<"\t"<<lat<<"\tgeomagnetic:\t"<<orbitalinfo->londip<<"\t"<<orbitalinfo->latdip<<endl;
1511    
1512          //          //
1513          // compute mag field components and L shell.          // compute mag field components and L shell.
1514          //          //
1515            if ( debug ) printf(" call igrf feldg \n");
1516          feldg_(&lat, &lon, &alt, &bnorth, &beast, &bdown, &babs);          feldg_(&lat, &lon, &alt, &bnorth, &beast, &bdown, &babs);
1517            if ( debug ) printf(" call igrf shellg \n");
1518          shellg_(&lat, &lon, &alt, &dimo, &xl, &icode, &bab1);          shellg_(&lat, &lon, &alt, &dimo, &xl, &icode, &bab1);
1519            if ( debug ) printf(" call igrf findb \n");
1520          findb0_(&stps, &bdel, &value, &bequ, &rr0);          findb0_(&stps, &bdel, &value, &bequ, &rr0);
1521          //          //
1522            if ( debug ) printf(" done igrf \n");
1523          orbitalinfo->Bnorth = bnorth;          orbitalinfo->Bnorth = bnorth;
1524          orbitalinfo->Beast = beast;          orbitalinfo->Beast = beast;
1525          orbitalinfo->Bdown = bdown;          orbitalinfo->Bdown = bdown;
1526          orbitalinfo->Babs = babs;          orbitalinfo->Babs = babs;
1527            orbitalinfo->M = dimo;
1528          orbitalinfo->BB0 = babs/bequ;          orbitalinfo->BB0 = babs/bequ;
1529          orbitalinfo->L = xl;                orbitalinfo->L = xl;      
1530          // Set Stormer vertical cutoff using L shell.          // Set Stormer vertical cutoff using L shell.
1531          orbitalinfo->cutoffsvl = 14.9/(xl*xl);          orbitalinfo->cutoffsvl = 14.295 / (xl*xl); //
1532            if(debug)cout << "L = " << xl << "\tM = " << dimo << "\tvertical cutoff:  "<< orbitalinfo->cutoffsvl << endl;
1533    
1534            /*
1535              ---------- Forwarded message ----------
1536              Date: Wed, 09 May 2012 12:16:47 +0200
1537              From: Alessandro Bruno <alessandro.bruno@ba.infn.it>
1538              To: Mirko Boezio <mirko.boezio@ts.infn.it>
1539              Cc: Francesco S. Cafagna <Francesco.Cafagna@ba.infn.it>
1540              Subject: Störmer vertical cutoff
1541    
1542              Ciao Mirko,
1543              volevo segnalarti che il valore dello Störmer vertical cutoff nel Level2 è
1544              sovrastimato di circa il 4%.
1545              Dopo un'approfondita analisi con l'IGRF-05 abbiamo ricavano un valore pari
1546              a: 14.295 / L^2 anzichè 14.9 / L^2, valore obsoleto in quanto riferito agli
1547              anni '50.
1548            */
1549            //14.9/(xl*xl);
1550            orbitalinfo->igrf_icode = icode;
1551          //          //
1552        };              }      
1553        //        //
1554        if ( debug ) printf(" pitch angle \n");        if ( debug ) printf(" pitch angle \n");
1555        //        //
1556        // pitch angles        // pitch angles
1557        //        //
1558        if ( orbitalinfo->mode != 10 && orbitalinfo->mode != 5 && orbitalinfo->mode !=7 && orbitalinfo->mode != 9 ){        if( orbitalinfo->TimeGap>0){
         //  
         Float_t Bx = -orbitalinfo->Bdown;                       //don't need for PamExp ExpOnly for all geography areas  
         Float_t By = orbitalinfo->Beast;                        //don't need for PamExp ExpOnly for all geography areas  
         Float_t Bz = orbitalinfo->Bnorth;                       //don't need for PamExp ExpOnly for all geography areas  
1559          //          //
1560          TMatrixD Fij = PO->ECItoGreenwich(PO->QuatoECI(orbitalinfo->q0,orbitalinfo->q1,orbitalinfo->q2,orbitalinfo->q3),orbitalinfo->absTime);          if ( debug ) printf(" timegap %f \n",orbitalinfo->TimeGap);
1561          TMatrixD Dij = PO->GreenwichtoGEO(orbitalinfo->lat,orbitalinfo->lon,Fij);          Float_t Bx = -orbitalinfo->Bdown;
1562            Float_t By = orbitalinfo->Beast;
1563            Float_t Bz = orbitalinfo->Bnorth;
1564    
1565            TMatrixD Qiji(3,3);
1566            TMatrixD Qij = PO->QuatoECI(orbitalinfo->q0,orbitalinfo->q1,orbitalinfo->q2,orbitalinfo->q3);
1567            TMatrixD Dij = PO->ECItoGEO(Qij,orbitalinfo->absTime,orbitalinfo->lat,orbitalinfo->lon);
1568    
1569    //10REDNEW
1570            /* If initial orientation data have reason to be inaccurate */
1571            Float_t tg = 0.00;
1572           Float_t tmptg;
1573           if(MU!=0){
1574    //      if(orbitalinfo->TimeGap>0 && errq==0 && azim==0){               // 10RED CHECK  (comparison between three metod of recovering orientation)
1575           if((atime>=RTstart[MU] && atime<RTstart[MU+1] && RTbank1[MU]==0 && RTbank2[MU]==0 && TMath::Abs(orbitalinfo->etha)>0.01) || ((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)))){
1576            /*  found in Rotation Table this data for this time interval*/
1577            if(atime<RTtime1[0])
1578              orbitalinfo->azim = 5;        //means that RotationTable no started yet
1579           else{
1580                    // search for angle betwean velosity and direction to north in tangential to Earth surfase plane in satellite position
1581                  Double_t bank=RTstart[MU];
1582                  Double_t tlat=orbitalinfo->lat;
1583    
1584                  tg=modf(orbitalinfo->TimeGap,&tg)*1000;
1585    
1586                  if(atime>=RTpluto1[MU] && atime<=RTpluto2[MU]){
1587                    Double_t kar=(RTbank2[MU]-RTbank1[MU])/(RTtime2[MU]-RTtime1[MU]);
1588                    Double_t bak=RTbank1[MU]-kar*RTtime1[MU];
1589                    bank=kar*atime+bak;
1590                  }
1591                  if(atime>=RTstart[MU] && atime<RTpluto1[MU]){
1592                     Double_t s_dBdt2=(RTbpluto1[MU]-RTbank1[MU])/(RTpluto1[MU]-RTstart[MU]);
1593                     Double_t s_t2=((Double_t)RTpluto1[MU]+(Double_t)RTstart[MU])/2. - RTstart[MU];
1594                     Double_t s_t1=RTstart[MU]-RTstart[MU];
1595                     Double_t s_k=s_dBdt2/(s_t2-s_t1);
1596                     Double_t s_b=-s_k*s_t1;
1597                     Double_t s_t3=RTpluto1[MU]-RTstart[MU];
1598                     Double_t s_b3=RTbpluto1[MU];
1599                     Double_t s_c=s_b3-0.5*s_k*s_t3*s_t3 -s_b*s_t3;
1600                     bank=0.5*s_k*(atime-RTstart[MU])*(atime-RTstart[MU]) + s_b*(atime-RTstart[MU]) + s_c;
1601                 }
1602                  if(atime>RTpluto2[MU] && atime<=RTstart[MU+1]){
1603                     Double_t s_dBdt2=(RTbpluto2[MU] - RTbank2[MU])/(RTpluto2[MU]-RTstart[MU+1]);
1604                     Double_t s_t2=((Double_t)RTpluto2[MU]+(Double_t)RTstart[MU+1])/2. - RTstart[MU];
1605                     Double_t s_t1=RTstart[MU+1]-RTstart[MU];
1606                     Double_t s_k=s_dBdt2/(s_t2-s_t1);
1607                     Double_t s_b=-s_k*s_t1;
1608                     Double_t s_t3=RTpluto2[MU]-RTstart[MU];
1609                     Double_t s_b3=RTbpluto2[MU];
1610                     Double_t s_c=s_b3-0.5*s_k*s_t3*s_t3 -s_b*s_t3;
1611                   bank=0.5*s_k*(atime-RTstart[MU])*(atime-RTstart[MU]) + s_b*(atime-RTstart[MU]) + s_c;
1612                 }
1613                  orbitalinfo->etha=bank;
1614                  Double_t spitch = 0.00001;  // temprary not zero to avoid problem with tranzition from Euler angles to orientation matrix
1615    
1616                    //Estimations of pitch angle of satellite
1617                  if(TMath::Abs(bank)>0.7){
1618                     Float_t spitch1=TMath::DegToRad()*0.7*diro;//RTdir1[MU];
1619                     Float_t spitch2=TMath::DegToRad()*0.7*diro;//RTdir2[MU];
1620                     Float_t kva=(spitch2-spitch1)/(RTtime2[MU]-RTtime1[MU]);
1621                     Float_t bva=spitch1-kva*RTtime1[MU];
1622                     spitch=kva*atime+bva;
1623                  }
1624    
1625                  //Calculate Yaw angle accordingly with fit, see picture FitYaw.jpg
1626                  Double_t yaw=0.00001;  // temprary not zero to avoid problem with tranzition from Euler angles to orientation matrix
1627                  if(TMath::Abs(tlat)<70)
1628                    yaw = -3.7e-8*tlat*tlat*tlat*tlat + 1.4e-7*tlat*tlat*tlat - 0.0005*tlat*tlat - 0.00025*tlat + 3.6;
1629                  yaw = diro*yaw;   //because should be different sign for ascending and descending orbits!
1630                  orbitalinfo->phi=yaw;
1631    
1632                  if(TMath::Abs(bank)>0.5 && TMath::Abs(yaw-orbitalinfo->phi)<3.0) yaw=orbitalinfo->phi;
1633    
1634    //            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
1635                  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
1636                  orbitalinfo->qkind = 1;
1637    
1638              //Qij = PO->GEOtoECI(Dij,orbitalinfo->absTime,orbitalinfo->lat,orbitalinfo->lon); // to convert from Dij to Qij
1639    
1640            } // end of if(atime<RTtime1[0]
1641            } // end of f(((orbitalinfo->TimeGap>60.0 && TMath...
1642           } // end of MU~=0
1643    
1644            TMatrixD qij = PO->ColPermutation(Qij);
1645            TMatrixD Fij = PO->ECItoGreenwich(Qij,orbitalinfo->absTime);
1646            TMatrixD Gij = PO->ColPermutation(Fij);
1647            Dij = PO->ECItoGEO(Qij,orbitalinfo->absTime,orbitalinfo->lat,orbitalinfo->lon);
1648          TMatrixD Iij = PO->ColPermutation(Dij);          TMatrixD Iij = PO->ColPermutation(Dij);
1649          //          TVector3 SP = PO->GetSunPosition(orbitalinfo->absTime);
1650            // go to Pamela reference frame from Resurs reference frame
1651            Float_t tmpy = SP.Y();
1652            SP.SetY(SP.Z());
1653            SP.SetZ(-tmpy);
1654            TVector3 SunZenith;
1655            SunZenith.SetMagThetaPhi(1,23.439281*TMath::DegToRad(),TMath::Pi()/2.);
1656            TVector3 SunMag = -SP;
1657            SunMag.Rotate(-45*TMath::DegToRad(),SunZenith);
1658            tmpy=SunMag.Y();
1659            SunMag.SetY(SunMag.Z());
1660            SunMag.SetZ(-tmpy);
1661    
1662          orbitalinfo->Iij.ResizeTo(Iij);          orbitalinfo->Iij.ResizeTo(Iij);
1663          orbitalinfo->Iij = Iij;          orbitalinfo->Iij = Iij;
1664          //          //
1665          A1 = Iij(0,2);          //      A1 = Iij(0,2);
1666          A2 = Iij(1,2);          //      A2 = Iij(1,2);
1667          A3 = Iij(2,2);          //      A3 = Iij(2,2);
1668          //                //
1669          //      orbitalinfo->pamzenitangle = (Float_t)PO->GetPitchAngle(1,0,0,A1,A2,A3);                        // Angle between zenit and Pamela's main axiz          //      orbitalinfo->pamzenitangle = (Float_t)PO->GetPitchAngle(1,0,0,A1,A2,A3);                        // Angle between zenit and Pamela's main axiz
1670          //      orbitalinfo->pamBangle = (Float_t)PO->GetPitchAngle(A1,A2,A3,Bx,By,Bz);                 // Angle between Pamela's main axiz and B          //      orbitalinfo->pamBangle = (Float_t)PO->GetPitchAngle(A1,A2,A3,Bx,By,Bz);                 // Angle between Pamela's main axiz and B
1671          //          //
1672          if ( !standalone && tof->ntrk() > 0 ){          if ( debug ) printf(" matrixes done  \n");
1673            if ( !standalone ){
1674              if ( debug ) printf(" !standalone \n");
1675            //            //
1676              // Standard tracking algorithm
1677              //
1678            Int_t nn = 0;            Int_t nn = 0;
1679              if ( verbose ) printf(" standard tracking \n");
1680            for(Int_t nt=0; nt < tof->ntrk(); nt++){              for(Int_t nt=0; nt < tof->ntrk(); nt++){  
1681              //              //
1682              ToFTrkVar *ptt = tof->GetToFTrkVar(nt);              ToFTrkVar *ptt = tof->GetToFTrkVar(nt);
1683                if (debug) cout<<"tof->ntrk() = "<<tof->ntrk()<<"\tptt->trkseqno = "<<ptt->trkseqno<<"\ttrke->ntrk() = "<<trke->ntrk()<<endl;
1684              Double_t E11x = ptt->xtr_tof[0]; // tr->x[0];              Double_t E11x = ptt->xtr_tof[0]; // tr->x[0];
1685              Double_t E11y = ptt->ytr_tof[0]; //tr->y[0];              Double_t E11y = ptt->ytr_tof[0]; //tr->y[0];
1686              Double_t E11z = zin[0];              Double_t E11z = zin[0];
# Line 1199  int OrbitalInfoCore(UInt_t run, TFile *f Line 1688  int OrbitalInfoCore(UInt_t run, TFile *f
1688              Double_t E22y = ptt->ytr_tof[3];//tr->y[3];              Double_t E22y = ptt->ytr_tof[3];//tr->y[3];
1689              Double_t E22z = zin[3];              Double_t E22z = zin[3];
1690              if ( (E11x < 100. && E11y < 100. && E22x < 100. && E22y < 100.) || ptt->trkseqno != -1  ){              if ( (E11x < 100. && E11y < 100. && E22x < 100. && E22y < 100.) || ptt->trkseqno != -1  ){
1691                  TrkTrack *mytrack = trke->GetStoredTrack(ptt->trkseqno);
1692                  Float_t rig=1/mytrack->GetDeflection();
1693                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));
1694                //              Double_t MyAzim = TMath::RadToDeg()*atan(TMath::Abs(E22y-E11y)/TMath::Abs(E22x-E11x));                //
               //              if(E22x-E11x>=0 && E22y-E11y <0) MyAzim =  360. - MyAzim;  
               //              if(E22x-E11x>=0 && E22y-E11y >=0) MyAzim = MyAzim;  
               //              if(E22x-E11x<0 && E22y-E11y >0) MyAzim = 180. - MyAzim;  
               //              if(E22x-E11x<0 && E22y-E11y <0) MyAzim = 180. + MyAzim;  
1695                Px = (E22x-E11x)/norm;                Px = (E22x-E11x)/norm;
1696                Py = (E22y-E11y)/norm;                Py = (E22y-E11y)/norm;
1697                Pz = (E22z-E11z)/norm;                Pz = (E22z-E11z)/norm;
# Line 1215  int OrbitalInfoCore(UInt_t run, TFile *f Line 1702  int OrbitalInfoCore(UInt_t run, TFile *f
1702                t_orb->Eij.ResizeTo(Eij);                t_orb->Eij.ResizeTo(Eij);
1703                t_orb->Eij = Eij;                t_orb->Eij = Eij;
1704                //                //
1705                TMatrixD Sij = PO->PamelatoGEO(Fij,Px,Py,Pz);                TMatrixD Sij = PO->PamelatoGEO(Gij,Px,Py,Pz);
1706                t_orb->Sij.ResizeTo(Sij);                t_orb->Sij.ResizeTo(Sij);
1707                t_orb->Sij = Sij;                t_orb->Sij = Sij;
1708                //                            //            
1709                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);
1710                //                //
1711                  // SunPosition in instrumental reference frame
1712                  TMatrixD Kij = PO->PamelatoGEO(qij,Px,Py,Pz);
1713                  TMatrixD Lij = PO->PamelatoGEO(qij,0,0,1);
1714                  t_orb->sunangle=(Float_t)PO->GetPitchAngle(Kij(0,0),Kij(1,0),Kij(2,0),-SP.X(),-SP.Y(),-SP.Z());
1715                  t_orb->sunmagangle=(Float_t)PO->GetPitchAngle(Kij(0,0),Kij(1,0),Kij(2,0),SunMag.X(),SunMag.Y(),SunMag.Z());
1716                //                //
               Double_t omega = PO->GetPitchAngle(Eij(0,0),Eij(1,0),Eij(2,0),cos(orbitalinfo->lon+TMath::Pi()/2)-sin(orbitalinfo->lon+TMath::Pi()/2),cos(orbitalinfo->lon+TMath::Pi()/2)+sin(orbitalinfo->lon+TMath::Pi()/2),1);  
1717                //                //
1718                t_orb->cutoff = 59.3/(pow(orbitalinfo->L,2)*pow((1+sqrt(1-pow(orbitalinfo->L,-3/2)*cos(omega))),2));                Double_t omega = PO->GetPitchAngle(-Eij(0,0),-Eij(1,0),-Eij(2,0),1,0,0) * TMath::DegToRad();
1719                  TVector3 Bxy(0,By,Bz);
1720                  TVector3 Exy(0,-Eij(1,0),-Eij(2,0));
1721                  Double_t dzeta=Bxy.Angle(Exy);
1722                  if (-Eij(1,0) < 0) dzeta=2.0*TMath::Pi() - dzeta;
1723                  
1724                  if(debug) cout << "omega = "<<omega*TMath::RadToDeg()<<"\tdzeta = "<<dzeta*TMath::RadToDeg()<<endl;
1725    
1726                  // Formula from D.F. Smart *, M.A. Shea [2005]; A review of geomagnetic cutoff rigidities for earth-orbiting spacecraft
1727                  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));
1728                  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));
1729                  if (debug) cout << "R = " << rig << "\tcutoff = " << t_orb->cutoff << endl;
1730    
1731                //                //
1732                if ( t_orb->pitch != t_orb->pitch ) t_orb->pitch = -1000.;                if ( t_orb->pitch != t_orb->pitch ) t_orb->pitch = -1000.;
1733                if ( t_orb->cutoff != t_orb->cutoff ) t_orb->cutoff = -1000.;                if ( t_orb->cutoff != t_orb->cutoff ) t_orb->cutoff = -1000.;
1734                  if ( t_orb->sunangle != t_orb->sunangle ) t_orb->sunangle = -1000.;
1735                  if ( t_orb->sunmagangle != t_orb->sunmagangle ) t_orb->sunmagangle = -1000.;
1736                //                //
1737                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);
1738                //                //
# Line 1236  int OrbitalInfoCore(UInt_t run, TFile *f Line 1741  int OrbitalInfoCore(UInt_t run, TFile *f
1741                //                //
1742                t_orb->Clear();                t_orb->Clear();
1743                //                //
1744              };              }
1745              //              //
1746            };            } // end standard tracking algorithm
1747    
1748              //
1749              // Code for extended tracking algorithm:
1750              //
1751              if ( hasNucleiTrk ){
1752                Int_t ttentry = 0;
1753                if ( verbose ) printf(" hasNucleiTrk \n");
1754                for(Int_t nt=0; nt < tcNucleiTof->GetEntries() ; nt++){  
1755                  //
1756                  if ( verbose ) printf(" got1\n");
1757                  ToFTrkVar *ptt = (ToFTrkVar*)(tcNucleiTof->At(nt));
1758                  if (verbose) cout<<" tcNucleiTof->GetEntries() = "<<tcNucleiTof->GetEntries()<<"\tptt->trkseqno = "<<ptt->trkseqno<<endl;
1759                  Double_t E11x = ptt->xtr_tof[0]; // tr->x[0];
1760                  Double_t E11y = ptt->ytr_tof[0]; //tr->y[0];
1761                  Double_t E11z = zin[0];
1762                  Double_t E22x = ptt->xtr_tof[3];//tr->x[3];
1763                  Double_t E22y = ptt->ytr_tof[3];//tr->y[3];
1764                  Double_t E22z = zin[3];
1765                  if ( (E11x < 100. && E11y < 100. && E22x < 100. && E22y < 100.) || ptt->trkseqno != -1  ){
1766                    TrkTrack *mytrack = (TrkTrack*)(tcNucleiTrk->At(ptt->trkseqno));
1767                    if ( verbose ) printf(" got tcNucleiTrk \n");
1768                    Float_t rig=1/mytrack->GetDeflection();
1769                    Double_t norm = sqrt(pow(E22x-E11x,2)+pow(E22y-E11y,2)+pow(E22z-E11z,2));
1770                    //
1771                    Px = (E22x-E11x)/norm;
1772                    Py = (E22y-E11y)/norm;
1773                    Pz = (E22z-E11z)/norm;
1774                    //
1775                    t_orb->trkseqno = ptt->trkseqno;
1776                    //
1777                    TMatrixD Eij = PO->PamelatoGEO(Iij,Px,Py,Pz);
1778                    t_orb->Eij.ResizeTo(Eij);      
1779                    t_orb->Eij = Eij;      
1780                    //
1781                    TMatrixD Sij = PO->PamelatoGEO(Gij,Px,Py,Pz);
1782                    t_orb->Sij.ResizeTo(Sij);
1783                    t_orb->Sij = Sij;
1784                    //          
1785                    t_orb->pitch = (Float_t)PO->GetPitchAngle(Eij(0,0),Eij(1,0),Eij(2,0),Bx,By,Bz);
1786                    //
1787                    // SunPosition in instrumental reference frame
1788                    TMatrixD Kij = PO->PamelatoGEO(qij,Px,Py,Pz);
1789                    TMatrixD Lij = PO->PamelatoGEO(qij,0,0,1);
1790                    t_orb->sunangle=(Float_t)PO->GetPitchAngle(Kij(0,0),Kij(1,0),Kij(2,0),-SP.X(),-SP.Y(),-SP.Z());
1791                    t_orb->sunmagangle=(Float_t)PO->GetPitchAngle(Kij(0,0),Kij(1,0),Kij(2,0),SunMag.X(),SunMag.Y(),SunMag.Z());
1792                    //
1793                    //
1794                    Double_t omega = PO->GetPitchAngle(-Eij(0,0),-Eij(1,0),-Eij(2,0),1,0,0) * TMath::DegToRad();
1795                    TVector3 Bxy(0,By,Bz);
1796                    TVector3 Exy(0,-Eij(1,0),-Eij(2,0));
1797                    Double_t dzeta=Bxy.Angle(Exy);
1798                    if (-Eij(1,0) < 0) dzeta=2.0*TMath::Pi() - dzeta;
1799                    
1800                    if(debug) cout << "omega = "<<omega*TMath::RadToDeg()<<"\tdzeta = "<<dzeta*TMath::RadToDeg()<<endl;
1801                    
1802                    // Formula from D.F. Smart *, M.A. Shea [2005]; A review of geomagnetic cutoff rigidities for earth-orbiting spacecraft
1803                    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));
1804                    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));
1805                    if (debug) cout << "R = " << rig << "\tcutoff = " << t_orb->cutoff << endl;
1806                    
1807                    //
1808                    if ( t_orb->pitch != t_orb->pitch ) t_orb->pitch = -1000.;
1809                    if ( t_orb->cutoff != t_orb->cutoff ) t_orb->cutoff = -1000.;
1810                    if ( t_orb->sunangle != t_orb->sunangle ) t_orb->sunangle = -1000.;
1811                    if ( t_orb->sunmagangle != t_orb->sunmagangle ) t_orb->sunmagangle = -1000.;
1812                    //
1813                    if ( debug ) printf(" orbitalinfo->cutoffsvl %f vitaly %f \n",orbitalinfo->cutoffsvl,t_orb->cutoff);
1814                    //
1815                    TClonesArray &tt1 = *torbNucleiTrk;
1816                    new(tt1[ttentry]) OrbitalInfoTrkVar(*t_orb);
1817                    ttentry++;
1818                    //
1819                    t_orb->Clear();
1820                    //
1821                  }
1822                  //
1823                }
1824              } // end standard tracking algorithm: nuclei
1825              if ( hasExtNucleiTrk ){
1826                Int_t ttentry = 0;
1827                if ( verbose ) printf(" hasExtNucleiTrk \n");
1828                for(Int_t nt=0; nt < tcExtNucleiTof->GetEntries() ; nt++){  
1829                  //
1830                  if ( verbose ) printf(" got2\n");
1831                  ToFTrkVar *ptt = (ToFTrkVar*)(tcExtNucleiTof->At(nt));
1832                  if (verbose) cout<<" tcExtNucleiTof->GetEntries() = "<<tcExtNucleiTof->GetEntries()<<"\tptt->trkseqno = "<<ptt->trkseqno<<endl;
1833                  Double_t E11x = ptt->xtr_tof[0]; // tr->x[0];
1834                  Double_t E11y = ptt->ytr_tof[0]; //tr->y[0];
1835                  Double_t E11z = zin[0];
1836                  Double_t E22x = ptt->xtr_tof[3];//tr->x[3];
1837                  Double_t E22y = ptt->ytr_tof[3];//tr->y[3];
1838                  Double_t E22z = zin[3];
1839                  if ( (E11x < 100. && E11y < 100. && E22x < 100. && E22y < 100.) || ptt->trkseqno != -1  ){
1840                    ExtTrack *mytrack = (ExtTrack*)(tcExtNucleiTrk->At(ptt->trkseqno));
1841                    if ( verbose ) printf(" got tcExtNucleiTrk \n");
1842                    Float_t rig=1/mytrack->GetDeflection();
1843                    Double_t norm = sqrt(pow(E22x-E11x,2)+pow(E22y-E11y,2)+pow(E22z-E11z,2));
1844                    //
1845                    Px = (E22x-E11x)/norm;
1846                    Py = (E22y-E11y)/norm;
1847                    Pz = (E22z-E11z)/norm;
1848                    //
1849                    t_orb->trkseqno = ptt->trkseqno;
1850                    //
1851                    TMatrixD Eij = PO->PamelatoGEO(Iij,Px,Py,Pz);
1852                    t_orb->Eij.ResizeTo(Eij);      
1853                    t_orb->Eij = Eij;      
1854                    //
1855                    TMatrixD Sij = PO->PamelatoGEO(Gij,Px,Py,Pz);
1856                    t_orb->Sij.ResizeTo(Sij);
1857                    t_orb->Sij = Sij;
1858                    //          
1859                    t_orb->pitch = (Float_t)PO->GetPitchAngle(Eij(0,0),Eij(1,0),Eij(2,0),Bx,By,Bz);
1860                    //
1861                    // SunPosition in instrumental reference frame
1862                    TMatrixD Kij = PO->PamelatoGEO(qij,Px,Py,Pz);
1863                    TMatrixD Lij = PO->PamelatoGEO(qij,0,0,1);
1864                    t_orb->sunangle=(Float_t)PO->GetPitchAngle(Kij(0,0),Kij(1,0),Kij(2,0),-SP.X(),-SP.Y(),-SP.Z());
1865                    t_orb->sunmagangle=(Float_t)PO->GetPitchAngle(Kij(0,0),Kij(1,0),Kij(2,0),SunMag.X(),SunMag.Y(),SunMag.Z());
1866                    //
1867                    //
1868                    Double_t omega = PO->GetPitchAngle(-Eij(0,0),-Eij(1,0),-Eij(2,0),1,0,0) * TMath::DegToRad();
1869                    TVector3 Bxy(0,By,Bz);
1870                    TVector3 Exy(0,-Eij(1,0),-Eij(2,0));
1871                    Double_t dzeta=Bxy.Angle(Exy);
1872                    if (-Eij(1,0) < 0) dzeta=2.0*TMath::Pi() - dzeta;
1873                    
1874                    if(debug) cout << "omega = "<<omega*TMath::RadToDeg()<<"\tdzeta = "<<dzeta*TMath::RadToDeg()<<endl;
1875                    
1876                    // Formula from D.F. Smart *, M.A. Shea [2005]; A review of geomagnetic cutoff rigidities for earth-orbiting spacecraft
1877                    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));
1878                    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));
1879                    if (debug) cout << "R = " << rig << "\tcutoff = " << t_orb->cutoff << endl;
1880                    
1881                    //
1882                    if ( t_orb->pitch != t_orb->pitch ) t_orb->pitch = -1000.;
1883                    if ( t_orb->cutoff != t_orb->cutoff ) t_orb->cutoff = -1000.;
1884                    if ( t_orb->sunangle != t_orb->sunangle ) t_orb->sunangle = -1000.;
1885                    if ( t_orb->sunmagangle != t_orb->sunmagangle ) t_orb->sunmagangle = -1000.;
1886                    //
1887                    if ( debug ) printf(" orbitalinfo->cutoffsvl %f vitaly %f \n",orbitalinfo->cutoffsvl,t_orb->cutoff);
1888                    //
1889                    TClonesArray &tt2 = *torbExtNucleiTrk;
1890                    new(tt2[ttentry]) OrbitalInfoTrkVar(*t_orb);
1891                    ttentry++;
1892                    //
1893                    t_orb->Clear();
1894                    //
1895                  }
1896                  //
1897                }
1898              } // end standard tracking algorithm: nuclei extended
1899             if ( hasExtTrk ){
1900                Int_t ttentry = 0;
1901                if ( verbose ) printf(" hasExtTrk \n");
1902                for(Int_t nt=0; nt < tcExtTof->GetEntries() ; nt++){  
1903                  //
1904                  if ( verbose ) printf(" got3\n");
1905                  ToFTrkVar *ptt = (ToFTrkVar*)(tcExtTof->At(nt));
1906                  if (verbose) cout<<" tcExtTof->GetEntries() = "<<tcExtTof->GetEntries()<<"\tptt->trkseqno = "<<ptt->trkseqno<<endl;
1907                  Double_t E11x = ptt->xtr_tof[0]; // tr->x[0];
1908                  Double_t E11y = ptt->ytr_tof[0]; //tr->y[0];
1909                  Double_t E11z = zin[0];
1910                  Double_t E22x = ptt->xtr_tof[3];//tr->x[3];
1911                  Double_t E22y = ptt->ytr_tof[3];//tr->y[3];
1912                  Double_t E22z = zin[3];
1913                  if ( (E11x < 100. && E11y < 100. && E22x < 100. && E22y < 100.) || ptt->trkseqno != -1  ){
1914                    ExtTrack *mytrack = (ExtTrack*)(tcExtTrk->At(ptt->trkseqno));
1915                    if ( verbose ) printf(" got tcExtTrk \n");
1916                    Float_t rig=1/mytrack->GetDeflection();
1917                    Double_t norm = sqrt(pow(E22x-E11x,2)+pow(E22y-E11y,2)+pow(E22z-E11z,2));
1918                    //
1919                    Px = (E22x-E11x)/norm;
1920                    Py = (E22y-E11y)/norm;
1921                    Pz = (E22z-E11z)/norm;
1922                    //
1923                    t_orb->trkseqno = ptt->trkseqno;
1924                    //
1925                    TMatrixD Eij = PO->PamelatoGEO(Iij,Px,Py,Pz);
1926                    t_orb->Eij.ResizeTo(Eij);      
1927                    t_orb->Eij = Eij;      
1928                    //
1929                    TMatrixD Sij = PO->PamelatoGEO(Gij,Px,Py,Pz);
1930                    t_orb->Sij.ResizeTo(Sij);
1931                    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);
1934                    //
1935                    // SunPosition in instrumental reference frame
1936                    TMatrixD Kij = PO->PamelatoGEO(qij,Px,Py,Pz);
1937                    TMatrixD Lij = PO->PamelatoGEO(qij,0,0,1);
1938                    t_orb->sunangle=(Float_t)PO->GetPitchAngle(Kij(0,0),Kij(1,0),Kij(2,0),-SP.X(),-SP.Y(),-SP.Z());
1939                    t_orb->sunmagangle=(Float_t)PO->GetPitchAngle(Kij(0,0),Kij(1,0),Kij(2,0),SunMag.X(),SunMag.Y(),SunMag.Z());
1940                    //
1941                    //
1942                    Double_t omega = PO->GetPitchAngle(-Eij(0,0),-Eij(1,0),-Eij(2,0),1,0,0) * TMath::DegToRad();
1943                    TVector3 Bxy(0,By,Bz);
1944                    TVector3 Exy(0,-Eij(1,0),-Eij(2,0));
1945                    Double_t dzeta=Bxy.Angle(Exy);
1946                    if (-Eij(1,0) < 0) dzeta=2.0*TMath::Pi() - dzeta;
1947                    
1948                    if(debug) cout << "omega = "<<omega*TMath::RadToDeg()<<"\tdzeta = "<<dzeta*TMath::RadToDeg()<<endl;
1949                    
1950                    // Formula from D.F. Smart *, M.A. Shea [2005]; A review of geomagnetic cutoff rigidities for earth-orbiting spacecraft
1951                    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));
1952                    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));
1953                    if (debug) cout << "R = " << rig << "\tcutoff = " << t_orb->cutoff << endl;
1954                    
1955                    //
1956                    if ( t_orb->pitch != t_orb->pitch ) t_orb->pitch = -1000.;
1957                    if ( t_orb->cutoff != t_orb->cutoff ) t_orb->cutoff = -1000.;
1958                    if ( t_orb->sunangle != t_orb->sunangle ) t_orb->sunangle = -1000.;
1959                    if ( t_orb->sunmagangle != t_orb->sunmagangle ) t_orb->sunmagangle = -1000.;
1960                    //
1961                    if ( debug ) printf(" orbitalinfo->cutoffsvl %f vitaly %f \n",orbitalinfo->cutoffsvl,t_orb->cutoff);
1962                    //
1963                    TClonesArray &tt3 = *torbExtTrk;
1964                    new(tt3[ttentry]) OrbitalInfoTrkVar(*t_orb);
1965                    ttentry++;
1966                    //
1967                    t_orb->Clear();
1968                    //
1969                  }
1970                  //
1971                }
1972              } // end standard tracking algorithm: extended
1973    
1974          } else {          } else {
1975            if ( debug ) printf(" mmm... mode %u standalone %i ntrk %i \n",orbitalinfo->mode,standalone,tof->ntrk());            if ( debug ) printf(" mmm... mode %u standalone  \n",orbitalinfo->mode);
1976          };          }
1977          //          //
1978        } else {        } else { // HERE IT MISS ALL CODE FOR EXTENDED TRACKING!
1979          if ( !standalone && tof->ntrk() > 0 ){          if ( !standalone ){
1980            //            //
1981              if ( verbose ) printf(" no orb info for tracks \n");
1982            Int_t nn = 0;            Int_t nn = 0;
1983            for(Int_t nt=0; nt < tof->ntrk(); nt++){              for(Int_t nt=0; nt < tof->ntrk(); nt++){  
1984              //              //
# Line 1260  int OrbitalInfoCore(UInt_t run, TFile *f Line 1993  int OrbitalInfoCore(UInt_t run, TFile *f
1993                //                            //            
1994                t_orb->pitch = -1000.;                t_orb->pitch = -1000.;
1995                //                //
1996                  t_orb->sunangle = -1000.;
1997                  //
1998                  t_orb->sunmagangle = -1000;
1999                  //
2000                t_orb->cutoff = -1000.;                t_orb->cutoff = -1000.;
2001                //                //
2002                new(tor[nn]) OrbitalInfoTrkVar(*t_orb);                new(tor[nn]) OrbitalInfoTrkVar(*t_orb);
# Line 1267  int OrbitalInfoCore(UInt_t run, TFile *f Line 2004  int OrbitalInfoCore(UInt_t run, TFile *f
2004                //                //
2005                t_orb->Clear();                t_orb->Clear();
2006                //                //
2007              };              }
2008              //              //
2009            };                }
2010          };            //
2011        };            // Code for extended tracking algorithm:
2012              //
2013              if ( hasNucleiTrk ){
2014                Int_t ttentry = 0;
2015                if ( verbose ) printf(" hasNucleiTrk \n");
2016                for(Int_t nt=0; nt < tcNucleiTof->GetEntries() ; nt++){  
2017                  //  
2018                  ToFTrkVar *ptt = (ToFTrkVar*)(tcNucleiTof->At(nt));
2019                  if ( ptt->trkseqno != -1  ){
2020                    //
2021                    t_orb->trkseqno = ptt->trkseqno;
2022                    //
2023                    t_orb->Eij = 0;
2024                    //
2025                    t_orb->Sij = 0;
2026                    //          
2027                    t_orb->pitch = -1000.;
2028                    //
2029                    t_orb->sunangle = -1000.;
2030                    //
2031                    t_orb->sunmagangle = -1000;
2032                    //
2033                    t_orb->cutoff = -1000.;
2034                    //
2035                    TClonesArray &tz1 = *torbNucleiTrk;
2036                    new(tz1[ttentry]) OrbitalInfoTrkVar(*t_orb);
2037                    ttentry++;
2038                    //
2039                    t_orb->Clear();
2040                    //
2041                  }
2042                  //
2043                }
2044              }
2045              if ( hasExtNucleiTrk ){
2046                Int_t ttentry = 0;
2047                if ( verbose ) printf(" hasExtNucleiTrk \n");
2048                for(Int_t nt=0; nt < tcExtNucleiTof->GetEntries() ; nt++){  
2049                  //
2050                  if ( verbose ) printf(" got2\n");
2051                  ToFTrkVar *ptt = (ToFTrkVar*)(tcExtNucleiTof->At(nt));
2052                  if ( ptt->trkseqno != -1  ){
2053                    //
2054                    t_orb->trkseqno = ptt->trkseqno;
2055                    //
2056                    t_orb->Eij = 0;
2057                    //
2058                    t_orb->Sij = 0;
2059                    //          
2060                    t_orb->pitch = -1000.;
2061                    //
2062                    t_orb->sunangle = -1000.;
2063                    //
2064                    t_orb->sunmagangle = -1000;
2065                    //
2066                    t_orb->cutoff = -1000.;
2067                    //
2068                    TClonesArray &tz2 = *torbExtNucleiTrk;
2069                    new(tz2[ttentry]) OrbitalInfoTrkVar(*t_orb);
2070                    ttentry++;
2071                    //
2072                    t_orb->Clear();
2073                    //
2074                  }
2075                  //
2076                }
2077              }
2078              if ( hasExtTrk ){
2079                Int_t ttentry = 0;
2080                if ( verbose ) printf(" hasExtTrk \n");
2081                for(Int_t nt=0; nt < tcExtTof->GetEntries() ; nt++){  
2082                  //
2083                  if ( verbose ) printf(" got3\n");
2084                  ToFTrkVar *ptt = (ToFTrkVar*)(tcExtTof->At(nt));
2085                  if ( ptt->trkseqno != -1  ){
2086                    //
2087                    t_orb->trkseqno = ptt->trkseqno;
2088                    //
2089                    t_orb->Eij = 0;
2090                    //
2091                    t_orb->Sij = 0;
2092                    //          
2093                    t_orb->pitch = -1000.;
2094                    //
2095                    t_orb->sunangle = -1000.;
2096                    //
2097                    t_orb->sunmagangle = -1000;
2098                    //
2099                    t_orb->cutoff = -1000.;
2100                    //
2101                    TClonesArray &tz3 = *torbExtNucleiTrk;
2102                    new(tz3[ttentry]) OrbitalInfoTrkVar(*t_orb);
2103                    ttentry++;
2104                    //
2105                    t_orb->Clear();
2106                    //
2107                  }
2108                  //
2109                }
2110              }
2111            }
2112          } // if( orbitalinfo->TimeGap>0){
2113        //        //
2114        // Fill the class        // Fill the class
2115        //        //
# Line 1284  int OrbitalInfoCore(UInt_t run, TFile *f Line 2122  int OrbitalInfoCore(UInt_t run, TFile *f
2122      // Here you may want to clear some variables before processing another run        // Here you may want to clear some variables before processing another run  
2123      //      //
2124    
2125      //gStyle->SetOptStat(000000);      if ( verbose ) printf(" Clear before new run \n");
     //gStyle->SetPalette(1);  
       
     /*TCanvas* c1 = new TCanvas("c1","",1200,800);  
     //c1->Divide(1,4);  
     c1->cd(1);  
     //q0testing->Draw("colz");  
     //c1->cd(2);  
     //q1testing->Draw("colz");  
     //c1->cd(3);  
     Pitchtesting->Draw("colz");  
     //c1->cd(4);  
     //q3testing->Draw("colz");  
     c1->SaveAs("9.Rollhyst.png");  
     delete c1;*/  
   
2126      delete dbtime;      delete dbtime;
2127      if ( L_QQ_Q_l_upper ) delete L_QQ_Q_l_upper;  
2128        if ( mcmdrc ) mcmdrc->Clear();
2129        mcmdrc = 0;
2130        
2131        if ( verbose ) printf(" Clear before new run1 \n");
2132      if ( L_QQ_Q_l_lower ) delete L_QQ_Q_l_lower;      if ( L_QQ_Q_l_lower ) delete L_QQ_Q_l_lower;
2133        if ( verbose ) printf(" Clear before new run2 \n");
2134        if ( L_QQ_Q_l_upper ) delete L_QQ_Q_l_upper;
2135        if ( verbose ) printf(" Clear before new run3 \n");
2136      if ( RYPang_upper ) delete RYPang_upper;      if ( RYPang_upper ) delete RYPang_upper;
2137        if ( verbose ) printf(" Clear before new run4 \n");
2138      if ( RYPang_lower ) delete RYPang_lower;      if ( RYPang_lower ) delete RYPang_lower;
2139    
2140        if ( l0tr ) l0tr->Delete();
2141        
2142        if ( verbose ) printf(" End run \n");
2143    
2144    }; // process all the runs    }; // process all the runs
2145        
2146    if (verbose) printf("\n Finished processing data \n");    if (verbose) printf("\n Finished processing data \n");
# Line 1338  int OrbitalInfoCore(UInt_t run, TFile *f Line 2174  int OrbitalInfoCore(UInt_t run, TFile *f
2174        };        };
2175        if (verbose) printf(" Finished successful copying!\n");        if (verbose) printf(" Finished successful copying!\n");
2176      };      };
2177        //if ( OrbitalInfotrclone )    OrbitalInfotrclone->Clear();        
2178        //if ( OrbitalInfotrclone )    OrbitalInfotrclone->Delete();        
2179    };    };
2180    //    //
2181    // Close files, delete old tree(s), write and close level2 file    // Close files, delete old tree(s), write and close level2 file
2182    //    //
2183    
2184    if ( l0File ) l0File->Close();    if ( l0File ) l0File->Close();
   if ( tempfile ) tempfile->Close();              
2185    if ( myfold ) gSystem->Unlink(tempname.str().c_str());    if ( myfold ) gSystem->Unlink(tempname.str().c_str());
2186    //    //
   if ( runinfo ) runinfo->Close();      
2187    if ( OrbitalInfotr ) OrbitalInfotr->SetName("OrbitalInfo");        if ( OrbitalInfotr ) OrbitalInfotr->SetName("OrbitalInfo");    
   if ( tof ) tof->Delete();  
   if ( ttof ) ttof->Delete();  
2188    //    //
2189    if ( file ){    if ( file ){
2190      file->cd();      file->cd();
2191      file->Write();      if ( OrbitalInfotr ) OrbitalInfotr->Write("OrbitalInfo", TObject::kOverwrite); // 10 RED bug fixed
2192    };    };
2193    //    //
2194      if (verbose) printf("\n Exiting...\n");
2195    
2196    if ( myfold ) gSystem->Unlink(OrbitalInfofolder.str().c_str());    if ( myfold ) gSystem->Unlink(OrbitalInfofolder.str().c_str());
2197    //    //
2198    // the end    // the end
# Line 1364  int OrbitalInfoCore(UInt_t run, TFile *f Line 2201  int OrbitalInfoCore(UInt_t run, TFile *f
2201      dbc->Close();      dbc->Close();
2202      delete dbc;      delete dbc;
2203    };    };
   if (verbose) printf("\n Exiting...\n");  
   if(OrbitalInfotr)OrbitalInfotr->Delete();  
2204    //    //
2205      if (verbose) printf("\n Exiting...\n");
2206      if ( tempfile ) tempfile->Close();            
2207      
2208    if ( PO ) delete PO;    if ( PO ) delete PO;
2209    if ( orbitalinfo ) delete orbitalinfo;    if ( gltle ) delete gltle;
2210    if ( orbitalinfoclone ) delete orbitalinfoclone;    if ( glparam ) delete glparam;
2211      if ( glparam2 ) delete glparam2;
2212      if ( glparam3 ) delete glparam3;
2213      if (verbose) printf("\n Exiting3...\n");
2214    if ( glroot ) delete glroot;    if ( glroot ) delete glroot;
2215      if (verbose) printf("\n Exiting4...\n");
2216      if ( runinfo ) runinfo->Close();    
2217    if ( runinfo ) delete runinfo;    if ( runinfo ) delete runinfo;
2218    
2219      if ( tcNucleiTrk ){
2220        tcNucleiTrk->Delete();
2221        delete tcNucleiTrk;
2222        tcNucleiTrk = NULL;
2223      }
2224      if ( tcExtNucleiTrk ){
2225        tcExtNucleiTrk->Delete();
2226        delete tcExtNucleiTrk;
2227        tcExtNucleiTrk = NULL;
2228      }
2229      if ( tcExtTrk ){
2230        tcExtTrk->Delete();
2231        delete tcExtTrk;
2232        tcExtTrk = NULL;
2233      }
2234    
2235      if ( tcNucleiTof ){
2236        tcNucleiTof->Delete();
2237        delete tcNucleiTof;
2238        tcNucleiTof = NULL;
2239      }
2240      if ( tcExtNucleiTof ){
2241        tcExtNucleiTof->Delete();
2242        delete tcExtNucleiTof;
2243        tcExtNucleiTof = NULL;
2244      }
2245      if ( tcExtTof ){
2246        tcExtTof->Delete();
2247        delete tcExtTof;
2248        tcExtTrk = NULL;
2249      }
2250    
2251    
2252      if ( tof ) delete tof;
2253      if ( trke ) delete trke;
2254    
2255      if ( debug ){  
2256      cout << "1   0x" << OrbitalInfotr << endl;
2257      cout << "2   0x" << OrbitalInfotrclone << endl;
2258      cout << "3   0x" << l0tr << endl;
2259      cout << "4   0x" << tempOrbitalInfo << endl;
2260      cout << "5   0x" << ttof << endl;
2261      }
2262      //
2263      if ( debug )  file->ls();
2264    //    //
2265    if(code < 0)  throw code;    if(code < 0)  throw code;
2266    return(code);    return(code);
# Line 1419  UInt_t holeq(Double_t lower,Double_t upp Line 2308  UInt_t holeq(Double_t lower,Double_t upp
2308    Bool_t R10l = false;     // Sign of R10 mode in lower quaternions array    Bool_t R10l = false;     // Sign of R10 mode in lower quaternions array
2309    Bool_t R10u = false;     // Sign of R10 mode in upper quaternions array    Bool_t R10u = false;     // Sign of R10 mode in upper quaternions array
2310    Bool_t insm = false;     // Sign that we inside quaternions array    Bool_t insm = false;     // Sign that we inside quaternions array
2311    Bool_t mxtml = false;    // Sign of mixt mode in lower quaternions array    //  Bool_t mxtml = false;    // Sign of mixt mode in lower quaternions array
2312    Bool_t mxtmu = false;    // Sign of mixt mode in upper quaternions array    //  Bool_t mxtmu = false;    // Sign of mixt mode in upper quaternions array
2313    Bool_t 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
2314    UInt_t NCQl = 6;       // Number of correct quaternions in lower array    UInt_t NCQl = 6;       // Number of correct quaternions in lower array
2315    UInt_t NCQu = 6;       // Number of correct quaternions in upper array    //  UInt_t NCQu = 6;       // Number of correct quaternions in upper array
2316    if (f>0){    if (f>0){
2317      insm = true;      insm = true;
2318      if(Qupper->time[f]-Qupper->time[f-1]==30) R10u = false;      if(Qupper->time[f]-Qupper->time[f-1]==30) R10u = false;
# Line 1435  UInt_t holeq(Double_t lower,Double_t upp Line 2324  UInt_t holeq(Double_t lower,Double_t upp
2324      if((Qlower->time[5]-Qlower->time[0]==150)&&(Qlower->time[1]-Qlower->time[0]==30)) R10l = false;      if((Qlower->time[5]-Qlower->time[0]==150)&&(Qlower->time[1]-Qlower->time[0]==30)) R10l = false;
2325      if((Qupper->time[5]-Qupper->time[0]==150)&&(Qupper->time[1]-Qupper->time[0]==30)) R10u = false;      if((Qupper->time[5]-Qupper->time[0]==150)&&(Qupper->time[1]-Qupper->time[0]==30)) R10u = false;
2326      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)){
2327        mxtml = true;        //      mxtml = true;
2328        for(UInt_t i = 1; i < 6; i++){        for(UInt_t i = 1; i < 6; i++){
2329          if(Qlower->time[i]-Qlower->time[0]==30*i) NCQl=i;          if(Qlower->time[i]-Qlower->time[0]==30*i) NCQl=i;
2330        }        }
2331      }      }
2332      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)){
2333        mxtmu = true;        //      mxtmu = true;
2334        for(UInt_t i = 1; i < 6; i++){        //      for(UInt_t i = 1; i < 6; i++){
2335          if(Qupper->time[i]-Qupper->time[0]==30*i) NCQu=i;        //        if(Qupper->time[i]-Qupper->time[0]==30*i) NCQu=i;
2336        }        //      }
2337      }      //    }
2338    }    }
2339        
2340    if(((upper-lower==1.5)||(upper-lower==3.)||(upper-lower==30.)||(upper-lower==31.5)||(upper-lower==33.)||(upper-lower==181.5)||(upper-lower==210.)||(upper-lower==211.5))&&!insm) npasm = true;    if(((upper-lower==1.5)||(upper-lower==3.)||(upper-lower==30.)||(upper-lower==31.5)||(upper-lower==33.)||(upper-lower==181.5)||(upper-lower==210.)||(upper-lower==211.5))&&!insm) npasm = true;
# Line 1477  void inclresize(vector<Double_t>& t,vect Line 2366  void inclresize(vector<Double_t>& t,vect
2366    Yaw.resize(sizee);    Yaw.resize(sizee);
2367  }  }
2368    
2369  //Find fitting sine functions for q0,q1,q2,q3 and Yaw-angle;  // geomagnetic calculation staff
2370  void sineparam(vector<Sine>& qsine, vector<Double_t>& qtime, vector<Float_t>& q, vector<Float_t>& Roll, vector<Float_t>& Pitch, Float_t limsin){  
2371    UInt_t mulast = 0;  //void GM_ScanIGRF(TString PATH, GMtype_Data *G0, GMtype_Data *G1, GMtype_Data *H1)
2372    UInt_t munow = 0;  void GM_ScanIGRF(TSQLServer *dbc, GMtype_Data *G0, GMtype_Data *G1, GMtype_Data *H1)
2373    UInt_t munext = 0;  {
2374    Bool_t increase = false;    GL_PARAM *glp = new GL_PARAM();
2375    Bool_t decrease = false;    Int_t parerror=glp->Query_GL_PARAM(1,304,dbc); // parameters stored in DB in GL_PRAM table  
2376    Bool_t Max_is_defined = false;    if ( parerror<0 ) {
2377    Bool_t Start_point_is_defined = false;      throw -902;
2378    Bool_t Period_is_defined = false;    }
2379    Bool_t Large_gap = false;          /*This function scans inputs G0, G1, and H1 of the IGRF table into 3 data arrays*/
2380    Bool_t normal_way = true;    //    TString SATH="/data03/Malakhov/pam9Malakhov/installed10/calib/orb-param/";
2381    Bool_t small_gap_on_ridge = false;          int i;
2382    Double_t t1 = 0;          double temp;
2383    Double_t t1A = 0;          char buffer[200];
2384    Int_t sinesize = 0;          FILE *IGRF;
2385    Int_t nfi = 0;          IGRF = fopen((glp->PATH+glp->NAME).Data(), "r");
2386    for(UInt_t mu = 0;mu<qtime.size();mu++){          //      IGRF = fopen(PATH+"IGRF.tab", "r");
2387      //cout<<"Roll["<<mu<<"] = "<<Roll[mu]<<endl;          G0->size = 25;
2388      if(TMath::Abs(Roll[mu])<1. && TMath::Abs(Pitch[mu])<1. && TMath::Abs(q[mu])<limsin){          G1->size = 25;
2389      //cout<<"q["<<mu<<endl<<"] = "<<q[mu]<<endl;          H1->size = 25;
2390      if(mulast!=0 && munow!=0 && munext!=0){mulast=munow;munow=munext;munext=mu;}          for( i = 0; i < 4; i++)
2391      if(munext==0 && munow!=0)munext=mu;          {
2392      if(munow==0 && mulast!=0)munow=mu;                  fgets(buffer, 200, IGRF);
     if(mulast==0)mulast=mu;  
       
     //cout<<"mulast = "<<mulast<<"\tmunow = "<<munow<<"\tmunext = "<<munext<<endl;  
     //Int_t ref;  
     //cin>>ref;  
     if(TMath::Abs(q[munow])>TMath::Abs(q[mulast]) && TMath::Abs(q[munow])>TMath::Abs(q[munext]) && q[mulast]*q[munext]>0 && qtime[munext]-qtime[mulast]>400)small_gap_on_ridge = true;  
     if(munext>mulast && (qtime[munext]-qtime[mulast]>=2000 || qtime[munext]-qtime[mulast]<0)){if(Large_gap){normal_way = false;Large_gap = false;}else{Large_gap = true;normal_way = false;}}else normal_way = true;  
     //if(normal_way)cout<<"Normal_Way"<<endl;  
     if(Large_gap || small_gap_on_ridge){  
       //cout<<"Large gap..."<<endl;  
       //if(small_gap_on_ridge)cout<<"small gap..."<<endl;  
       //cout<<"q["<<mulast<<"] = "<<q[mulast]<<"\tq["<<munow<<"] = "<<q[munow]<<"\tq["<<munext<<"] = "<<q[munext]<<endl;  
       //cout<<"qtime["<<mulast<<"] = "<<qtime[mulast]<<"\tqtime["<<munow<<"] = "<<qtime[munow]<<"\tqtime["<<munext<<"] = "<<qtime[munext]<<endl;  
       increase = false;  
       decrease = false;  
       if(nfi>0){  
         qsine.resize(qsine.size()-1);  
         sinesize = qsine.size();  
         //cout<<"nfi was larger then zero"<<endl;  
       }else{  
         //cout<<"nfi was not larger then zero :( nfi = "<<nfi<<endl;  
         //cout<<"qsine.size = "<<qsine.size()<<endl;  
         if(!Period_is_defined){  
           //cout<<"Period was defined"<<endl;  
           if(qsine.size()>1){  
             qsine[sinesize-1].b = qsine[sinesize-2].b;  
             qsine[sinesize-1].c = qsine[sinesize-2].c;  
           }else{  
             qsine[sinesize-1].b = TMath::Pi()/1591.54;  
             qsine[sinesize-1].c = qsine[sinesize-1].startPoint;  
           }  
         }  
         if(!Max_is_defined){  
           //cout<<"Max was already defined"<<endl;  
           if(qsine.size()>1)qsine[sinesize-1].A = qsine[sinesize-2].A;else qsine[sinesize-1].A = limsin;  
         }  
         qsine[sinesize-1].NeedFit = true;  
       }  
       qsine[sinesize-1].finishPoint = qtime[munow];  
       //cout<<"finish point before large gap = "<<qtime[munow]<<endl;  
       nfi = 0;  
       Max_is_defined = false;  
       Start_point_is_defined = false;  
       Period_is_defined = false;  
       small_gap_on_ridge = false;  
     }  
     //cout<<"Slope "<<increase<<"\t"<<decrease<<endl;  
     //cout<<"mulast = "<<mulast<<"\tmunow = "<<munow<<"\tmunext = "<<munext<<endl;  
     if((munext>munow) && (munow>mulast) && normal_way){  
       if(!increase && !decrease){  
         //cout<<"Normal way have started"<<endl;  
         qsine.resize(qsine.size()+1);  
         sinesize = qsine.size();  
         qsine[sinesize-1].startPoint=qtime[mulast];  
         if(q[munext]>q[munow] && q[munow]>q[mulast]) increase = true;  
         if(q[munext]<q[munow] && q[munow]<q[mulast]) decrease = true;  
       }  
       //if(TMath::Abs(q[munow])>TMath::Abs(q[mulast]) && TMath::Abs(q[munow])>TMath::Abs(q[munext]) && TMath::Abs(q[munow])>limsin && qtime[munow]-qtime[mulast]>=400 || qtime[munext]-qtime[munow]>=400){small_gap_on_ridge = true;mu--;continue;}  
       if(TMath::Abs(q[munow])>TMath::Abs(q[mulast]) && TMath::Abs(q[munow])>TMath::Abs(q[munext]) && TMath::Abs(q[munow])>0.9*limsin && qtime[munow]-qtime[mulast]<400 && qtime[munext]-qtime[munow]<400){  
         //cout<<"Max point is qtime = "<<qtime[munow]<<"\tq = "<<q[munow]<<endl;  
         if(q[munow]>q[mulast]){  
           increase = false;  
           decrease = true;  
         }  
         if(q[munow]<q[mulast]){  
           increase = true;  
           decrease = false;  
2393          }          }
2394          if(Max_is_defined && !Start_point_is_defined){          fscanf(IGRF, "g 1 0 %lf ", &G0->element[0]);
2395            Double_t qPer = qtime[munow]-t1A;          for(i = 1; i <= 22; i++)
2396            if(qPer>1000){          {
2397              //cout<<"qsine["<<sinesize-1<<"] = "<<qPer<<" = "<<qtime[munow]<<" - "<<t1A<<"\tlim = "<<limsin<<endl;                  fscanf(IGRF ,"%lf ", &G0->element[i]);
             qsine[sinesize-1].b=TMath::Pi()/qPer;  
             if(decrease)qsine[sinesize-1].c=-qsine[sinesize-1].b*t1A;  
             if(increase)qsine[sinesize-1].c=-qsine[sinesize-1].b*(t1A-qPer);  
             Period_is_defined = true;  
           }  
2398          }          }
2399          Max_is_defined = true;          fscanf(IGRF ,"%lf\n", &temp);
2400          qsine[sinesize-1].A = TMath::Abs(q[munow]);          G0->element[23] = temp * 5 + G0->element[22];
2401          if(Start_point_is_defined && Period_is_defined){          G0->element[24] = G0->element[23] + 5 * temp;
2402            qsine[sinesize-1].finishPoint = qtime[munow];          fscanf(IGRF, "g 1 1 %lf ", &G1->element[0]);
2403            nfi++;          for(i = 1; i <= 22; i++)
2404            qsine[sinesize-1].NeedFit = false;          {
2405            Max_is_defined = false;                  fscanf( IGRF, "%lf ", &G1->element[i]);
           Start_point_is_defined = false;  
           Period_is_defined = false;  
           qsine.resize(qsine.size()+1);  
           sinesize = qsine.size();  
           qsine[sinesize-1].startPoint = qtime[munow];  
2406          }          }
2407          if(!Start_point_is_defined) t1A=qtime[munow];          fscanf(IGRF, "%lf\n", &temp);
2408        }          G1->element[23] = temp * 5 + G1->element[22];
2409        //if((q[munow]>=0 && q[mulast]<=0) || (q[munow]<=0 && q[mulast]>=0))cout<<"cross zero point diference = "<<qtime[munext] - qtime[mulast]<<"\tqlast = "<<qtime[mulast]<<"\tqnow = "<<qtime[munow]<<"\tqnext = "<<qtime[munext]<<endl;          G1->element[24] = temp * 5 + G1->element[23];
2410        if(((q[munow]>=0 && q[mulast]<=0) || (q[munow]<=0 && q[mulast]>=0)) && qtime[munow]-qtime[mulast]<2000 && qtime[munext]-qtime[munow]<2000){          fscanf(IGRF, "h 1 1 %lf ", &H1->element[0]);
2411          Double_t tcrosszero = 0;          for(i = 1; i <= 22; i++)
2412          //cout<<"cross zero point...qtime = "<<qtime[munow]<<endl;          {
2413          if(q[munow]==0.) tcrosszero = qtime[munow];else                  fscanf( IGRF, "%lf ", &H1->element[i]);
           if(q[mulast]==0.)tcrosszero = qtime[mulast];else{  
             Double_t k_ = (q[munow]-q[mulast])/(qtime[munow]-qtime[mulast]);  
             Double_t b_ = q[munow]-k_*qtime[munow];  
             tcrosszero = -b_/k_;  
           }  
         if(Start_point_is_defined){  
           //cout<<"Start Point allready defined"<<endl;  
           Double_t qPer = tcrosszero - t1;  
           qsine[sinesize-1].b = TMath::Pi()/qPer;  
           //cout<<"qsine["<<sinesize-1<<"].b = "<<TMath::Pi()/qPer<<endl;  
           Period_is_defined = true;  
           Float_t x0 = 0;  
           if(decrease)x0 = t1;  
           if(increase)x0 = tcrosszero;  
           qsine[sinesize-1].c = -qsine[sinesize-1].b*x0;  
           if(Max_is_defined){  
             //cout<<"Max was previous defined"<<endl;  
             qsine[sinesize-1].finishPoint = qtime[munow];  
             nfi++;  
             qsine[sinesize-1].NeedFit = false;  
             Max_is_defined = false;  
             t1 = tcrosszero;  
             Start_point_is_defined = true;  
             Period_is_defined = false;  
             qsine.resize(qsine.size()+1);  
             sinesize = qsine.size();  
             qsine[sinesize-1].startPoint = qtime[munow];  
           }  
         }else{  
           t1 = tcrosszero;  
           Start_point_is_defined = true;  
2414          }          }
2415        }          fscanf(IGRF, "%lf\n", &temp);
2416      }          H1->element[23] = temp * 5 + H1->element[22];
2417      }          H1->element[24] = temp * 5 + H1->element[23];
2418    }    if ( glp ) delete glp;
2419    } /*GM_ScanIGRF*/
2420    
2421    //cout<<"FINISH SINE INTERPOLATION FUNCTION..."<<endl<<endl;  void GM_SetEllipsoid(GMtype_Ellipsoid *Ellip)
2422  }  {
2423            /*This function sets the WGS84 reference ellipsoid to its default values*/
2424            Ellip->a        =                       6378.137; /*semi-major axis of the ellipsoid in */
2425            Ellip->b        =                       6356.7523142;/*semi-minor axis of the ellipsoid in */
2426            Ellip->fla      =                       1/298.257223563;/* flattening */
2427            Ellip->eps      =                       sqrt(1- ( Ellip->b *    Ellip->b) / (Ellip->a * Ellip->a ));  /*first eccentricity */
2428            Ellip->epssq    =                       (Ellip->eps * Ellip->eps);   /*first eccentricity squared */
2429            Ellip->re       =                       6371.2;/* Earth's radius */
2430    } /*GM_SetEllipsoid*/
2431    
2432    
2433    void GM_EarthCartToDipoleCartCD(GMtype_Pole Pole, GMtype_CoordCartesian EarthCoord, GMtype_CoordCartesian *DipoleCoords)
2434    {
2435            /*This function converts from Earth centered cartesian coordinates to dipole centered cartesian coordinates*/
2436            double X, Y, Z, CosPhi, SinPhi, CosLambda, SinLambda;
2437            CosPhi = cos(TMath::DegToRad()*Pole.phi);
2438            SinPhi = sin(TMath::DegToRad()*Pole.phi);
2439            CosLambda = cos(TMath::DegToRad()*Pole.lambda);
2440            SinLambda = sin(TMath::DegToRad()*Pole.lambda);
2441            X = EarthCoord.x;
2442            Y = EarthCoord.y;
2443            Z = EarthCoord.z;
2444            
2445            /*These equations are taken from a document by Wallace H. Campbell*/
2446            DipoleCoords->x = X * CosPhi * CosLambda + Y * CosPhi * SinLambda - Z * SinPhi;
2447            DipoleCoords->y = -X * SinLambda + Y * CosLambda;
2448            DipoleCoords->z = X * SinPhi * CosLambda + Y * SinPhi * SinLambda + Z * CosPhi;
2449    } /*GM_EarthCartToDipoleCartCD*/
2450    
2451    void GM_GeodeticToSpherical(GMtype_Ellipsoid Ellip, GMtype_CoordGeodetic CoordGeodetic, GMtype_CoordSpherical *CoordSpherical)
2452    {
2453            double CosLat, SinLat, rc, xp, zp; /*all local variables */
2454            /*
2455            ** Convert geodetic coordinates, (defined by the WGS-84
2456            ** reference ellipsoid), to Earth Centered Earth Fixed Cartesian
2457            ** coordinates, and then to spherical coordinates.
2458            */
2459    
2460            CosLat = cos(TMath::DegToRad()*CoordGeodetic.phi);
2461            SinLat = sin(TMath::DegToRad()*CoordGeodetic.phi);
2462    
2463            /* compute the local radius of curvature on the WGS-84 reference ellipsoid */
2464    
2465            rc = Ellip.a / sqrt(1.0 - Ellip.epssq * SinLat * SinLat);
2466    
2467            /* compute ECEF Cartesian coordinates of specified point (for longitude=0) */
2468    
2469            xp = (rc + CoordGeodetic.HeightAboveEllipsoid) * CosLat;
2470            zp = (rc*(1.0 - Ellip.epssq) + CoordGeodetic.HeightAboveEllipsoid) * SinLat;
2471    
2472            /* compute spherical radius and angle lambda and phi of specified point */
2473    
2474            CoordSpherical->r = sqrt(xp * xp + zp * zp);
2475            CoordSpherical->phig = TMath::RadToDeg()*asin(zp / CoordSpherical->r);     /* geocentric latitude */
2476            CoordSpherical->lambda = CoordGeodetic.lambda;                   /* longitude */
2477    } /*GM_GeodeticToSpherical*/
2478    
2479    void GM_PoleLocation(GMtype_Model Model, GMtype_Pole *Pole)
2480    {
2481            /*This function finds the location of the north magnetic pole in spherical coordinates.  The equations are
2482            **from Wallace H. Campbell's Introduction to Geomagnetic Fields*/
2483    
2484            Pole->phi = TMath::RadToDeg()*-atan(sqrt(Model.h1 * Model.h1 + Model.g1 * Model.g1)/Model.g0);
2485            Pole->lambda = TMath::RadToDeg()*atan(Model.h1/Model.g1);
2486    } /*GM_PoleLocation*/
2487    
2488    void GM_SphericalToCartesian(GMtype_CoordSpherical CoordSpherical, GMtype_CoordCartesian *CoordCartesian)
2489    {
2490            /*This function converts spherical coordinates into Cartesian coordinates*/
2491            double CosPhi = cos(TMath::DegToRad()*CoordSpherical.phig);
2492            double SinPhi = sin(TMath::DegToRad()*CoordSpherical.phig);
2493            double CosLambda = cos(TMath::DegToRad()*CoordSpherical.lambda);
2494            double SinLambda = sin(TMath::DegToRad()*CoordSpherical.lambda);
2495            
2496            CoordCartesian->x = CoordSpherical.r * CosPhi * CosLambda;
2497            CoordCartesian->y = CoordSpherical.r * CosPhi * SinLambda;
2498            CoordCartesian->z = CoordSpherical.r * SinPhi;
2499    } /*GM_SphericalToCartesian*/
2500    
2501    void GM_TimeAdjustCoefs(Float_t year, Float_t jyear, GMtype_Data g0d, GMtype_Data g1d, GMtype_Data h1d, GMtype_Model *Model)
2502    {
2503            /*This function calls GM_LinearInterpolation for the coefficients to estimate the value of the
2504            **coefficient for the given date*/
2505            int index;
2506            double x;
2507            index = (year - GM_STARTYEAR) / 5;
2508            x = (jyear - GM_STARTYEAR) / 5;
2509            Model->g0 = GM_LinearInterpolation(index, index+1, g0d.element[index], g0d.element[index+1], x);
2510            Model->g1 = GM_LinearInterpolation(index, index+1, g1d.element[index], g1d.element[index+1], x);
2511            Model->h1 = GM_LinearInterpolation(index, index+1, h1d.element[index], h1d.element[index+1], x);
2512    } /*GM_TimeAdjustCoefs*/
2513    
2514    double GM_LinearInterpolation(double x1, double x2, double y1, double y2, double x)
2515    {
2516            /*This function takes a linear interpolation between two given points for x*/
2517            double weight, y;
2518            weight  = (x - x1) / (x2 - x1);
2519            y = y1 * (1 - weight) + y2 * weight;
2520            return y;
2521    }/*GM_LinearInterpolation*/
2522    
2523    void GM_CartesianToSpherical(GMtype_CoordCartesian CoordCartesian, GMtype_CoordSpherical *CoordSpherical)
2524    {
2525            /*This function converts a point from Cartesian coordinates into spherical coordinates*/
2526            double X, Y, Z;
2527            
2528            X = CoordCartesian.x;
2529            Y = CoordCartesian.y;
2530            Z = CoordCartesian.z;
2531    
2532            CoordSpherical->r = sqrt(X * X + Y * Y + Z * Z);
2533            CoordSpherical->phig = TMath::RadToDeg()*asin(Z / (CoordSpherical->r));
2534            CoordSpherical->lambda = TMath::RadToDeg()*atan2(Y, X);
2535    } /*GM_CartesianToSpherical*/
Legend:
Removed from v.1.47  
changed lines
  Added in v.1.74
  ViewVC Help
Powered by ViewVC 1.1.23