/[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.3 by mocchiut, Fri Jun 30 09:21:58 2006 UTC revision 1.35 by mocchiut, Thu Dec 11 10:08:19 2008 UTC
# Line 13  Line 13 
13  #include <TClassEdit.h>  #include <TClassEdit.h>
14  #include <TObject.h>  #include <TObject.h>
15  #include <TList.h>  #include <TList.h>
16  #include <TArrayL.h>  #include <TArrayI.h>
17  #include <TSystem.h>  #include <TSystem.h>
18  #include <TSystemDirectory.h>  #include <TSystemDirectory.h>
19  #include <TString.h>  #include <TString.h>
# Line 23  Line 23 
23  #include <TSQLRow.h>  #include <TSQLRow.h>
24  #include <TSQLResult.h>  #include <TSQLResult.h>
25  //  //
26    // RunInfo header
27    //
28    #include <RunInfo.h>
29    #include <GLTables.h>
30    //
31  // YODA headers  // YODA headers
32  //  //
33  #include <PamelaRun.h>  #include <PamelaRun.h>
 #include <RegistryEvent.h>  
34  #include <PscuHeader.h>  #include <PscuHeader.h>
35  #include <PscuEvent.h>  #include <PscuEvent.h>
36  #include <EventHeader.h>  #include <EventHeader.h>
37  #include <RegistryEvent.h>  #include <mcmd/McmdEvent.h>
38  //  #include <mcmd/McmdRecord.h>
 // RunInfo header  
 //  
 #include <RunInfo.h>  
 #include <GLTables.h>  
39  //  //
40  // This program headers  // This program headers
41  //  //
42  #include <OrbitalInfo.h>  #include <OrbitalInfo.h>
 #include <OrbitalInfoCore.h>  
43  #include <OrbitalInfoVerl2.h>  #include <OrbitalInfoVerl2.h>
44    #include <OrbitalInfoCore.h>
45    #include <InclinationInfo.h>
46    
47  using namespace std;  using namespace std;
48    
# Line 49  using namespace std; Line 50  using namespace std;
50  // CORE ROUTINE  // CORE ROUTINE
51  //  //
52  //  //
53    int OrbitalInfoCore(UInt_t run, TFile *file, GL_TABLES *glt, Int_t OrbitalInfoargc, char *OrbitalInfoargv[]){
54  int OrbitalInfoCore(ULong64_t run, TFile *file, TSQLServer *dbc, Int_t OrbitalInfoargc, char *OrbitalInfoargv[]){    //
55    Int_t i = 0;    Int_t i = 0;
56      TString host = glt->CGetHost();
57      TString user = glt->CGetUser();
58      TString psw = glt->CGetPsw();
59      TSQLServer *dbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data());
60      //
61      stringstream myquery;
62      myquery.str("");
63      myquery << "SET time_zone='+0:00'";
64      dbc->Query(myquery.str().c_str());
65    //    //
66    TString processFolder = "OrbitalInfoFolder";    TString processFolder = Form("OrbitalInfoFolder_%u",run);
67    //    //
68    // Set these to true to have a very verbose output.    // Set these to true to have a very verbose output.
69    //    //
70    Bool_t debug = false;    Bool_t debug = false;
71    //    //
72    Bool_t verbose = false;    Bool_t verbose = false;
73      //
74      Bool_t standalone = false;
75      //
76    if ( OrbitalInfoargc > 0 ){    if ( OrbitalInfoargc > 0 ){
77      i = 0;      i = 0;
78      while ( i < OrbitalInfoargc ){      while ( i < OrbitalInfoargc ){
# Line 73  int OrbitalInfoCore(ULong64_t run, TFile Line 85  int OrbitalInfoCore(ULong64_t run, TFile
85        };        };
86        if ( (!strcmp(OrbitalInfoargv[i],"--debug")) || (!strcmp(OrbitalInfoargv[i],"-g")) ) {        if ( (!strcmp(OrbitalInfoargv[i],"--debug")) || (!strcmp(OrbitalInfoargv[i],"-g")) ) {
87          verbose = true;          verbose = true;
88            debug = true;
89        };        };
90        if ( (!strcmp(OrbitalInfoargv[i],"--verbose")) || (!strcmp(OrbitalInfoargv[i],"-v")) ) {        if ( (!strcmp(OrbitalInfoargv[i],"--verbose")) || (!strcmp(OrbitalInfoargv[i],"-v")) ) {
91          verbose = true;          verbose = true;
92        };        };
93          if ( (!strcmp(OrbitalInfoargv[i],"--standalone")) ) {
94            standalone = true;
95          };
96          if ( (!strcmp(OrbitalInfoargv[i],"--calculate-pitch")) ) {
97            standalone = false;
98          };
99        i++;        i++;
100      };      };
101    };    };
# Line 84  int OrbitalInfoCore(ULong64_t run, TFile Line 103  int OrbitalInfoCore(ULong64_t run, TFile
103    const char* outDir = gSystem->DirName(gSystem->DirName(file->GetPath()));    const char* outDir = gSystem->DirName(gSystem->DirName(file->GetPath()));
104    //    //
105    TTree *OrbitalInfotr = 0;    TTree *OrbitalInfotr = 0;
106    Long64_t nevents = 0LL;    UInt_t nevents = 0;
107      UInt_t neventsm = 0;
108    //    //
109    // variables needed to reprocess data    // variables needed to reprocess data
110    //    //
111      Long64_t maxsize = 10000000000LL;  
112      TTree::SetMaxTreeSize(maxsize);
113      //
114    TString OrbitalInfoversion;    TString OrbitalInfoversion;
115    ItoRunInfo *runinfo = 0;    ItoRunInfo *runinfo = 0;
116    TArrayL *runlist = 0;    TArrayI *runlist = 0;
117    TTree *OrbitalInfotrclone = 0;    TTree *OrbitalInfotrclone = 0;
118    Bool_t reproc = false;    Bool_t reproc = false;
119    Bool_t reprocall = false;    Bool_t reprocall = false;
# Line 99  int OrbitalInfoCore(ULong64_t run, TFile Line 122  int OrbitalInfoCore(ULong64_t run, TFile
122    UInt_t numbofrun = 0;    UInt_t numbofrun = 0;
123    stringstream ftmpname;    stringstream ftmpname;
124    TString fname;    TString fname;
125    Long64_t totfileentries = 0ULL;    UInt_t totfileentries = 0;
126    Long64_t idRun = 0LL;    UInt_t idRun = 0;
127      //
128      // My variables. Vitaly.
129      //
130      //  UInt_t iev = 0;
131      //  UInt_t j3 = 0;
132      UInt_t oi = 0;
133      Int_t tmpSize = 0;
134    //    //
135    // variables needed to handle error signals    // variables needed to handle error signals
136    //    //
# Line 116  int OrbitalInfoCore(ULong64_t run, TFile Line 146  int OrbitalInfoCore(ULong64_t run, TFile
146    //    //
147    TFile *l0File = 0;    TFile *l0File = 0;
148    TTree *l0tr = 0;    TTree *l0tr = 0;
149    TBranch *l0registry = 0;    TTree *l0trm = 0;
   pamela::RegistryEvent *l0reg=0;  
150    // EM: open also header branch    // EM: open also header branch
151    TBranch *l0head = 0;    TBranch *l0head = 0;
152    pamela::EventHeader *eh = 0;    pamela::EventHeader *eh = 0;
153    pamela::PscuHeader *ph = 0;    pamela::PscuHeader *ph = 0;
154      pamela::McmdEvent *mcmdev = 0;
155      pamela::McmdRecord *mcmdrc = 0;
156    // end EM    // end EM
157      
158      //  pamela::RunHeaderEvent *reh = new pamela::RunHeaderEvent;
159      //  pamela::EventHeader    *eH  = new pamela::EventHeader;
160      
161    //    //
162    // Define other basic variables    // Define other basic variables
163    //    //
# Line 131  int OrbitalInfoCore(ULong64_t run, TFile Line 166  int OrbitalInfoCore(ULong64_t run, TFile
166    stringstream file3;    stringstream file3;
167    stringstream qy;    stringstream qy;
168    Int_t totevent = 0;    Int_t totevent = 0;
169    ULong64_t atime = 0ULL;    UInt_t atime = 0;
170    Int_t ei = 0;    UInt_t re = 0;
171    Int_t re = 0;    UInt_t ik = 0;
172    
173      // Position
174      Float_t lon, lat, alt;
175    
176      //
177      // IGRF stuff
178      //
179      float dimo = 0.0; // dipole moment (computed from dat files)
180      float bnorth, beast, bdown, babs;
181      float xl; // L value
182      float icode; // code value for L accuracy (see fortran code)
183      float bab1; // What's  the difference with babs?
184      float stps = 0.005; // step size for field line tracing
185      float bdel = 0.01; // required accuracy
186      float bequ;  // equatorial b value (also called b_0)
187      bool value = 0; // false if bequ is not the minimum b value
188      float rr0; // equatorial radius normalized to earth radius
189    
190    //    //
191    // Working filename    // Working filename
192    //    //
# Line 148  int OrbitalInfoCore(ULong64_t run, TFile Line 201  int OrbitalInfoCore(ULong64_t run, TFile
201    TTree *tempOrbitalInfo = 0;    TTree *tempOrbitalInfo = 0;
202    stringstream tempname;    stringstream tempname;
203    stringstream OrbitalInfofolder;    stringstream OrbitalInfofolder;
204      Bool_t myfold = false;
205    tempname.str("");    tempname.str("");
206    tempname << outDir;    tempname << outDir;
207    tempname << "/" << processFolder.Data();    tempname << "/" << processFolder.Data();
208    OrbitalInfofolder.str("");    OrbitalInfofolder.str("");
209    OrbitalInfofolder << tempname.str().c_str();    OrbitalInfofolder << tempname.str().c_str();
   gSystem->MakeDirectory(OrbitalInfofolder.str().c_str());  
210    tempname << "/OrbitalInfotree_run";    tempname << "/OrbitalInfotree_run";
211    tempname << run << ".root";      tempname << run << ".root";  
212    //    //
213    // DB classes    // DB classes
214    //    //
215    GL_ROOT *glroot = new GL_ROOT();    GL_ROOT *glroot = new GL_ROOT();
216      GL_TIMESYNC *dbtime = 0;
217      GL_TLE *gltle = new GL_TLE();
218      //
219      //Quaternions classes
220      //
221      Quaternions *L_QQ_Q_l_lower = new Quaternions();
222      InclinationInfo *RYPang_lower = new InclinationInfo();
223      Quaternions *L_QQ_Q_l_upper = new Quaternions();
224      InclinationInfo *RYPang_upper = new InclinationInfo();
225      
226      cEci eCi;
227      
228      // Initialize fortran routines!!!
229      Int_t ltp2 = 0;
230      Int_t ltp3 = 0;
231      Int_t uno = 1;
232      char *niente = " ";
233      GL_PARAM *glparam = new GL_PARAM();
234      GL_PARAM *glparam2 = new GL_PARAM();
235      Int_t parerror=glparam->Query_GL_PARAM(1,301,dbc); // parameters stored in DB in GL_PRAM table
236      //
237      // Orientation variables
238      //
239      UInt_t evfrom = 0;
240      UInt_t jumped = 0;
241      Int_t itr = -1;    
242      Double_t A1;
243      Double_t A2;
244      Double_t A3;
245      Double_t Px = 0;
246      Double_t Py = 0;      
247      Double_t Pz = 0;  
248      TTree *ttof = 0;
249      ToFLevel2 *tof = new ToFLevel2();
250      OrientationInfo *PO = new OrientationInfo();
251      Int_t nz = 6;
252      Float_t zin[6];
253      Int_t nevtofl2 = 0;
254      //  
255      if ( parerror<0 ) {
256        code = parerror;
257        goto closeandexit;
258      };
259      ltp2 = (Int_t)(glparam->PATH+glparam->NAME).Length();
260      if ( verbose ) printf(" Reading Earth's Magnetic Field parameter file: %s \n",(glparam->PATH+glparam->NAME).Data());
261      //
262      parerror=glparam2->Query_GL_PARAM(1,302,dbc); // parameters stored in DB in GL_PRAM table
263      if ( parerror<0 ) {
264        code = parerror;
265        goto closeandexit;
266      };
267      ltp3 = (Int_t)(glparam2->PATH+glparam2->NAME).Length();
268      if ( verbose ) printf(" Reading Earth's Magnetic Field parameter file: %s \n",(glparam2->PATH+glparam2->NAME).Data());
269      //
270      initize_((char *)niente,&uno,(char *)(glparam->PATH+glparam->NAME).Data(),&ltp2,(char *)(glparam2->PATH+glparam2->NAME).Data(),&ltp3);
271      //
272      // End IGRF stuff//
273      //
274      for (Int_t ip=0;ip<nz;ip++){
275        zin[ip] = tof->GetZTOF(tof->GetToFPlaneID(ip));
276      };
277      //
278      if ( !standalone ){
279        //
280        // Does it contain the Tracker tree?
281        //
282        ttof = (TTree*)file->Get("ToF");
283        if ( !ttof ) {
284          if ( verbose ) printf(" OrbitalInfo - ERROR: no tof tree\n");
285          code = -900;
286          goto closeandexit;
287        };
288        ttof->SetBranchAddress("ToFLevel2",&tof);  
289        nevtofl2 = ttof->GetEntries();
290      };
291    //    //
292    // Let's start!    // Let's start!
293    //    //
# Line 167  int OrbitalInfoCore(ULong64_t run, TFile Line 295  int OrbitalInfoCore(ULong64_t run, TFile
295    // if run != 0 we must process only that run but first we have to check if the tree MyDetector2 already exist in the file    // if run != 0 we must process only that run but first we have to check if the tree MyDetector2 already exist in the file
296    // if it exists we are reprocessing data and we must delete that entries, if not we must create it.    // if it exists we are reprocessing data and we must delete that entries, if not we must create it.
297    //    //
298    if ( run == 0ULL )  reproc = true;    if ( run == 0 )  reproc = true;
299    //    //
300    //    //
301    // Output file is "outputfile"    // Output file is "outputfile"
# Line 214  int OrbitalInfoCore(ULong64_t run, TFile Line 342  int OrbitalInfoCore(ULong64_t run, TFile
342    // number of run to be processed    // number of run to be processed
343    //    //
344    numbofrun = runinfo->GetNoRun();    numbofrun = runinfo->GetNoRun();
345      UInt_t totnorun = runinfo->GetRunEntries();
346    //    //
347    // Try to access the OrbitalInfo tree in the file, if it exists we are reprocessing data if not we are processing a new run    // Try to access the OrbitalInfo tree in the file, if it exists we are reprocessing data if not we are processing a new run
348    //    //
# Line 224  int OrbitalInfoCore(ULong64_t run, TFile Line 353  int OrbitalInfoCore(ULong64_t run, TFile
353      // tree does not exist, we are not reprocessing      // tree does not exist, we are not reprocessing
354      //      //
355      reproc = false;      reproc = false;
356      if ( run == 0ULL ){      if ( run == 0 ){
357        if (verbose) printf(" OrbitalInfo - WARNING: you are reprocessing data but OrbitalInfo tree does not exist!\n");        if (verbose) printf(" OrbitalInfo - WARNING: you are reprocessing data but OrbitalInfo tree does not exist!\n");
358      }      }
359      if ( runinfo->IsReprocessing() && run != 0ULL ) {      if ( runinfo->IsReprocessing() && run != 0 ) {
360        if (verbose) printf(" OrbitalInfo - WARNING: it seems you are not reprocessing data but OrbitalInfo\n versioning information already exists in RunInfo.\n");        if (verbose) printf(" OrbitalInfo - WARNING: it seems you are not reprocessing data but OrbitalInfo\n versioning information already exists in RunInfo.\n");
361      }      }
362    } else {    } else {
363      //      //
364      // tree exists, we are reprocessing data. Are we reprocessing a single run or all the file?      // tree exists, we are reprocessing data. Are we reprocessing a single run or all the file?
365      //      //
366        OrbitalInfotrclone->SetAutoSave(900000000000000LL);
367      reproc = true;      reproc = true;
368      //      //
369      //      //
370      if (verbose) printf("\n Preparing the pre-processing...\n");      if (verbose) printf("\n Preparing the pre-processing...\n");
371      //      //
372      if ( run == 0ULL ){      if ( run == 0 || totnorun == 1 ){
373        //        //
374        // we are reprocessing all the file        // we are reprocessing all the file
375        // if we are reprocessing everything we don't need to copy any old event and we can just work with the new tree and delete the old one immediately        // if we are reprocessing everything we don't need to copy any old event and we can just work with the new tree and delete the old one immediately
# Line 254  int OrbitalInfoCore(ULong64_t run, TFile Line 384  int OrbitalInfoCore(ULong64_t run, TFile
384        //        //
385        reprocall = false;        reprocall = false;
386        //        //
387        if (verbose) printf("\n OrbitalInfo - WARNING: Reprocessing run number %llu \n",run);        if (verbose) printf("\n OrbitalInfo - WARNING: Reprocessing run number %u \n",run);
388        //        //
389        // copying old tree to a new file        // copying old tree to a new file
390        //        //
391          gSystem->MakeDirectory(OrbitalInfofolder.str().c_str());
392          myfold = true;
393        tempfile = new TFile(tempname.str().c_str(),"RECREATE");        tempfile = new TFile(tempname.str().c_str(),"RECREATE");
394        tempOrbitalInfo = OrbitalInfotrclone->CloneTree(-1,"fast");        tempOrbitalInfo = OrbitalInfotrclone->CloneTree(-1,"fast");
395        tempOrbitalInfo->SetName("OrbitalInfo-old");        tempOrbitalInfo->SetName("OrbitalInfo-old");
# Line 277  int OrbitalInfoCore(ULong64_t run, TFile Line 409  int OrbitalInfoCore(ULong64_t run, TFile
409    //    //
410    file->cd();    file->cd();
411    OrbitalInfotr = new TTree("OrbitalInfo-new","PAMELA OrbitalInfo data");    OrbitalInfotr = new TTree("OrbitalInfo-new","PAMELA OrbitalInfo data");
412      OrbitalInfotr->SetAutoSave(900000000000000LL);
413      orbitalinfo->Set();//ELENA **TEMPORANEO?**
414    OrbitalInfotr->Branch("OrbitalInfo","OrbitalInfo",&orbitalinfo);    OrbitalInfotr->Branch("OrbitalInfo","OrbitalInfo",&orbitalinfo);
415    //    //
416    if ( reproc && !reprocall ){    if ( reproc && !reprocall ){
# Line 285  int OrbitalInfoCore(ULong64_t run, TFile Line 419  int OrbitalInfoCore(ULong64_t run, TFile
419      //      //
420      tempfile = new TFile(tempname.str().c_str(),"READ");      tempfile = new TFile(tempname.str().c_str(),"READ");
421      OrbitalInfotrclone = (TTree*)tempfile->Get("OrbitalInfo-old");      OrbitalInfotrclone = (TTree*)tempfile->Get("OrbitalInfo-old");
422        OrbitalInfotrclone->SetAutoSave(900000000000000LL);
423      OrbitalInfotrclone->SetBranchAddress("OrbitalInfo",&orbitalinfoclone);      OrbitalInfotrclone->SetBranchAddress("OrbitalInfo",&orbitalinfoclone);
424      //            //      
425      if ( nobefrun > 0 ){      if ( nobefrun > 0 ){
426        if (verbose){        if (verbose){
427        printf("\n Pre-processing: copying events from the old tree before the processed run\n");            printf("\n Pre-processing: copying events from the old tree before the processed run\n");  
428        printf(" Copying %u events in the file which are before the beginning of the run %llu \n",nobefrun,run);          printf(" Copying %u events in the file which are before the beginning of the run %u \n",nobefrun,run);
429        printf(" Start copying at event number 0, end copying at event number %u \n",nobefrun);          printf(" Start copying at event number 0, end copying at event number %u \n",nobefrun);
430        }        }
431        for (UInt_t j = 0; j < nobefrun; j++){        for (UInt_t j = 0; j < nobefrun; j++){
432          //          //
# Line 299  int OrbitalInfoCore(ULong64_t run, TFile Line 434  int OrbitalInfoCore(ULong64_t run, TFile
434          //          //
435          // copy orbitalinfoclone to mydec          // copy orbitalinfoclone to mydec
436          //          //
         //      orbitalinfo = new OrbitalInfo();  
437          orbitalinfo->Clear();          orbitalinfo->Clear();
438            //
439          memcpy(&orbitalinfo,&orbitalinfoclone,sizeof(orbitalinfoclone));          memcpy(&orbitalinfo,&orbitalinfoclone,sizeof(orbitalinfoclone));
440          //          //
441          // Fill entry in the new tree          // Fill entry in the new tree
# Line 322  int OrbitalInfoCore(ULong64_t run, TFile Line 457  int OrbitalInfoCore(ULong64_t run, TFile
457      //      //
458      // retrieve the first run ID to be processed using the RunInfo list      // retrieve the first run ID to be processed using the RunInfo list
459      //      //
460        
461      idRun = runlist->At(irun);      idRun = runlist->At(irun);
462      if (verbose){      if (verbose){
463        printf("\n\n\n ####################################################################### \n");        printf("\n\n\n ####################################################################### \n");
# Line 329  int OrbitalInfoCore(ULong64_t run, TFile Line 465  int OrbitalInfoCore(ULong64_t run, TFile
465        printf(" ####################################################################### \n\n\n");        printf(" ####################################################################### \n\n\n");
466      }      }
467      //      //
468      runinfo->ID_REG_RUN = 0ULL;      runinfo->ID_ROOT_L0 = 0;
469      //      //
470      // store in the runinfo class the GL_RUN variables for our run      // store in the runinfo class the GL_RUN variables for our run
471      //      //
472      sgnl = 0;      sgnl = 0;
473      sgnl = runinfo->GetRunInfo(idRun);      sgnl = runinfo->GetRunInfo(idRun);
474      if ( sgnl ){      if ( sgnl ){
475        //printf("\n OrbitalInfo - ERROR: RunInfo exited with non-zero status\n");        if ( debug ) printf("\n OrbitalInfo - ERROR: RunInfo exited with non-zero status\n");
476        code = sgnl;        code = sgnl;
477        goto closeandexit;        goto closeandexit;
478      } else {      } else {
# Line 345  int OrbitalInfoCore(ULong64_t run, TFile Line 481  int OrbitalInfoCore(ULong64_t run, TFile
481      //      //
482      // now you can access that variables using the RunInfo class this way runinfo->ID_REG_RUN      // now you can access that variables using the RunInfo class this way runinfo->ID_REG_RUN
483      //      //
484      if ( runinfo->ID_REG_RUN == 0 ){      if ( runinfo->ID_ROOT_L0 == 0 ){
485        //printf("\n OrbitalInfo - ERROR: no run with ID_RUN = %i \n\n Exiting... \n\n",(int)idRun);        if ( debug ) printf("\n OrbitalInfo - ERROR: no run with ID_RUN = %u \n\n Exiting... \n\n",idRun);
486        code = -5;        code = -5;
487        goto closeandexit;            goto closeandexit;    
488      };      };
489      //      //
490        // prepare the timesync for the db
491        //
492        dbtime = new GL_TIMESYNC(runinfo->ID_ROOT_L0,"ID",dbc);
493      
494        //
495      // Search in the DB the path and name of the LEVEL0 file to be processed.      // Search in the DB the path and name of the LEVEL0 file to be processed.
496      //      //
497      glroot->Query_GL_ROOT(runinfo->ID_REG_RUN,dbc);      glroot->Query_GL_ROOT(runinfo->ID_ROOT_L0,dbc);
498      //      //
499      ftmpname.str("");      ftmpname.str("");
500      ftmpname << glroot->PATH.Data() << "/";      ftmpname << glroot->PATH.Data() << "/";
501      ftmpname << glroot->NAME.Data();      ftmpname << glroot->NAME.Data();
502      fname = ftmpname.str().c_str();      fname = ftmpname.str().c_str();
503        ftmpname.str("");
504      //      //
505      // print out informations      // print out informations
506      //      //
507      totevent = runinfo->EV_REG_PHYS_TO - runinfo->EV_REG_PHYS_FROM + 1;      totevent = runinfo->NEVENTS;
508        evfrom = runinfo->EV_FROM;
509        //cout<<"totevents = "<<totevent<<"\n";
510      if (verbose){      if (verbose){
511        printf("\n LEVEL0 data file: %s \n",fname.Data());        printf("\n LEVEL0 data file: %s \n",fname.Data());
512        printf(" RUN HEADER absolute time is:  %llu \n",runinfo->RUNHEADER_TIME);        printf(" RUN HEADER absolute time is:  %u \n",runinfo->RUNHEADER_TIME);
513        printf(" RUN TRAILER absolute time is: %llu \n",runinfo->RUNTRAILER_TIME);        printf(" RUN TRAILER absolute time is: %u \n",runinfo->RUNTRAILER_TIME);
514        printf(" %i events to be processed for run %llu: from %i to %i (reg entries)\n\n",totevent,idRun,runinfo->EV_REG_PHYS_FROM,runinfo->EV_REG_PHYS_TO);        printf(" %i events to be processed for run %u: from %i to %i \n\n",totevent,idRun,runinfo->EV_FROM+1,runinfo->EV_FROM+totevent);
515      }//      }//
516        //
517        //    if ( !totevent ) goto closeandexit;
518      // Open Level0 file      // Open Level0 file
519      l0File = new TFile(fname.Data());      l0File = new TFile(fname.Data());
520      if ( !l0File ) {      if ( !l0File ) {
521        //printf(" OrbitalInfo - ERROR: problems opening Level0 file\n");        if ( debug ) printf(" OrbitalInfo - ERROR: problems opening Level0 file\n");
522        code = -6;        code = -6;
523        goto closeandexit;        goto closeandexit;
524      };      };
525      l0tr = (TTree*)l0File->Get("Physics");      l0tr = (TTree*)l0File->Get("Physics");
526      if ( !l0tr ) {      if ( !l0tr ) {
527        //printf(" OrbitalInfo - ERROR: no Physics tree in Level0 file\n");        if ( debug ) printf(" OrbitalInfo - ERROR: no Physics tree in Level0 file\n");
528        l0File->Close();        l0File->Close();
529        code = -7;        code = -7;
530        goto closeandexit;        goto closeandexit;
# Line 386  int OrbitalInfoCore(ULong64_t run, TFile Line 532  int OrbitalInfoCore(ULong64_t run, TFile
532      // EM: open header branch as well      // EM: open header branch as well
533      l0head = l0tr->GetBranch("Header");      l0head = l0tr->GetBranch("Header");
534      if ( !l0head ) {      if ( !l0head ) {
535        //if ( verbose ) printf(" CALORIMETER - ERROR: no Header branch in Level0 tree\n");        if ( debug ) printf(" OrbitalInfo - ERROR: no Header branch in Level0 tree\n");
536        l0File->Close();        l0File->Close();
537        code = -8;        code = -8;
538        goto closeandexit;            goto closeandexit;    
539      };      };
540      l0tr->SetBranchAddress("Header", &eh);      l0tr->SetBranchAddress("Header", &eh);
541      // end EM      // end EM
542      l0registry = l0tr->GetBranch("Registry");      nevents = l0head->GetEntries();
     if ( !l0registry ) {  
       //printf(" OrbitalInfo - ERROR: no Registry branch in Level0 tree\n");  
       l0File->Close();  
       code = -9;  
       goto closeandexit;      
     };  
     //  
     l0tr->SetBranchAddress("Registry", &l0reg);  
     //  
     nevents = l0registry->GetEntries();  
543      //      //
544      if ( nevents < 1 ) {      if ( nevents < 1 && totevent ) {
545        //printf(" OrbitalInfo - ERROR: Level0 file is empty\n\n");        if ( debug ) printf(" OrbitalInfo - ERROR: Level0 file is empty\n\n");
546        l0File->Close();        l0File->Close();
547        code = -11;        code = -11;
548        goto closeandexit;        goto closeandexit;
549      };      };
550      //      //
551      if ( runinfo->EV_REG_PHYS_TO > nevents-1 ) {      if ( runinfo->EV_TO > nevents-1 && totevent ) {
552        //printf(" OrbitalInfo - ERROR: too few entries in the registry tree\n");        if ( debug ) printf(" OrbitalInfo - ERROR: too few entries in the registry tree\n");
553        l0File->Close();        l0File->Close();
554        code = -12;        code = -12;
555        goto closeandexit;        goto closeandexit;
556      };      };
557      //      //
558    //     TTree *tp = (TTree*)l0File->Get("RunHeader");
559    //     tp->SetBranchAddress("Header", &eH);
560    //     tp->SetBranchAddress("RunHeader", &reh);
561    //     tp->GetEntry(0);
562    //     ph = eH->GetPscuHeader();
563    //     ULong_t TimeSync = reh->LAST_TIME_SYNC_INFO;
564    //     ULong_t ObtSync = reh->OBT_TIME_SYNC;    
565    //     if ( debug ) printf(" 1 TimeSync %lu ObtSync %lu DeltaOBT %lu\n",TimeSync,ObtSync,TimeSync-ObtSync);
566    //
567        ULong_t TimeSync = (ULong_t)dbtime->GetTimesync();
568        ULong_t ObtSync = (ULong_t)(dbtime->GetObt0()/1000);
569        ULong_t DeltaOBT = TimeSync - ObtSync;
570    
571        if ( debug ) printf(" 2 TimeSync %lu ObtSync %lu DeltaOBT %lu\n",(ULong_t)(dbtime->GetTimesync()/1000),(ULong_t)dbtime->GetObt0(),TimeSync-ObtSync);
572        
573        l0trm = (TTree*)l0File->Get("Mcmd");
574        neventsm = l0trm->GetEntries();
575        //    neventsm = 0;
576        //
577        if (neventsm == 0){
578          if ( debug ) printf("InclinationInfo - WARNING: No quaternions in this File");
579          //      l0File->Close();
580          code = 900;
581          //      goto closeandexit;
582        }
583        //
584        
585        l0trm->SetBranchAddress("Mcmd", &mcmdev);
586        //    l0trm->SetBranchAddress("Header", &eh);
587        //
588        //
589        //
590        UInt_t mctren = 0;    
591        UInt_t mcreen = 0;  
592        UInt_t numrec = 0;
593        //
594        Double_t upperqtime = 0;
595        Double_t lowerqtime = 0;
596        
597        Double_t incli = 0;
598        oi = 0;
599        UInt_t ooi = 0;
600        //
601        // init quaternions sync
602        //
603        Bool_t isf = true;
604        Int_t fgh = 0;
605        //
606      // run over all the events of the run      // run over all the events of the run
607      //      //
608      if (verbose) printf("\n Ready to start! \n\n Processed events: \n\n");      if (verbose) printf("\n Ready to start! \n\n Processed events: \n\n");
609      //      //
610      for ( re = runinfo->EV_REG_PHYS_FROM; re <= runinfo->EV_REG_PHYS_TO; re++){       //
611        for ( re = runinfo->EV_FROM; re < (runinfo->EV_FROM+runinfo->NEVENTS); re++){
612        
613        //        //
614        if ( procev%1000 == 0 && procev > 0 && verbose) printf(" %iK \n",procev/1000);            if ( procev%1000 == 0 && procev > 0 && verbose ) printf(" %iK \n",procev/1000);  
615          if ( debug ) printf(" %i \n",procev);      
616        //        //
617        l0registry->GetEntry(re);        l0head->GetEntry(re);
618        //        //
619        // absolute time of this event        // absolute time of this event
620        //        //
621        atime = l0reg->absTime;        ph = eh->GetPscuHeader();
622        //        atime = dbtime->DBabsTime(ph->GetOrbitalTime());
       // physics events is at entry number ei where  
       //  
       ei = l0reg->event;  
623        //        //
624        // paranoid check        // paranoid check
625        //        //
626        if ( (atime > runinfo->RUNTRAILER_TIME) || (atime < runinfo->RUNHEADER_TIME)  ) {        if ( (atime > (runinfo->RUNTRAILER_TIME+1)) || (atime < (runinfo->RUNHEADER_TIME-1))  ) {
627          if (verbose) printf(" OrbitalInfo - WARNING: event at time outside the run time window, skipping it\n");          if (verbose) printf(" OrbitalInfo - WARNING: event at time outside the run time window, skipping it\n");
628          goto jumpev;          jumped++;
629    //      debug = true;
630            continue;
631          }
632    
633          //
634          // retrieve tof informations
635          //
636          if ( !reprocall ){
637            itr = nobefrun + (re - evfrom - jumped);
638            //itr = re-(46438+200241);
639          } else {
640            itr = runinfo->GetFirstEntry() + (re - evfrom - jumped);
641          };
642          //
643          if ( !standalone ){
644            if ( itr > nevtofl2 ){  
645              if ( verbose ) printf(" OrbitalInfo - ERROR: no tof events with entry = %i in Level2 file\n",itr);
646              if ( debug ) printf(" nobefrun %u re %u evfrom %u jumped %u reprocall %i \n",nobefrun,re,evfrom,jumped,reprocall);
647              l0File->Close();
648              code = -901;
649              goto closeandexit;
650            };
651            //
652            tof->Clear();
653            //
654            ttof->GetEntry(itr);
655            //
656        };        };
657        //        //
658        procev++;        procev++;
# Line 449  int OrbitalInfoCore(ULong64_t run, TFile Line 660  int OrbitalInfoCore(ULong64_t run, TFile
660        // start processing        // start processing
661        //        //
662        orbitalinfo->Clear();        orbitalinfo->Clear();
663        //orbitalinfo = new OrbitalInfo();        //
664        orbitalinfo->absTime = l0reg->absTime;        OrbitalInfoTrkVar *t_orb = new OrbitalInfoTrkVar();
665        // EM: add OBT and plt_num infos from the header        if( !(orbitalinfo->OrbitalInfoTrk) ) orbitalinfo->OrbitalInfoTrk = new TClonesArray("OrbitalInfoTrkVar",2);
666        l0head->GetEntry(ei);        TClonesArray &tor = *orbitalinfo->OrbitalInfoTrk;
667        ph = eh->GetPscuHeader();        //
668          // Fill OBT, pkt_num and absTime
669          //      
670        orbitalinfo->pkt_num = ph->GetCounter();        orbitalinfo->pkt_num = ph->GetCounter();
671        orbitalinfo->OBT = ph->GetOrbitalTime();        orbitalinfo->OBT = ph->GetOrbitalTime();
672        // end EM        orbitalinfo->absTime = atime;
673          //
674          // Propagate the orbit from the tle time to atime, using SGP(D)4.
675          //
676          cCoordGeo coo;
677          float jyear=0;    
678          //
679          if(atime >= gltle->GetToTime()) {
680            if ( !gltle->Query(atime, dbc) ){
681              //      
682              // Compute the magnetic dipole moment.
683              //
684              UInt_t year, month, day, hour, min, sec;
685              //
686              TTimeStamp t = TTimeStamp(atime, kTRUE);
687              t.GetDate(kTRUE, 0, &year, &month, &day);
688              t.GetTime(kTRUE, 0, &hour, &min, &sec);
689              jyear = (float) year
690                + (month*31.+ (float) day)/365.
691                + (hour*3600.+min*60.+(float)sec)/(24*3600*365.);
692              //
693              feldcof_(&jyear, &dimo); // get dipole moment for year
694            } else {
695              code = -56;
696              goto closeandexit;
697            };
698          }
699          coo = getCoo(atime, gltle->GetFromTime(), gltle->GetTle());
700          //
701          cOrbit orbits(*gltle->GetTle());
702          //
703          if ( debug ) printf(" I am Here \n");
704          //
705          // synchronize with quaternions data
706          //
707          if ( isf && neventsm>0 ){
708            if ( debug ) printf(" I am here \n");
709            //
710            // First event
711            //
712            isf = false;
713            upperqtime = atime;
714            lowerqtime = runinfo->RUNHEADER_TIME;
715            for ( ik = 0; ik < neventsm; ik++){
716              l0trm->GetEntry(ik);
717              tmpSize = mcmdev->Records->GetEntries();
718              numrec = tmpSize;
719              for (Int_t j3 = 0;j3<tmpSize;j3++){
720                if ( debug ) printf(" eh eh eh \n");
721                mcmdrc = (pamela::McmdRecord*)mcmdev->Records->At(j3);
722                if ((int)mcmdrc->ID1 == 226){
723                  L_QQ_Q_l_upper->fill(mcmdrc->McmdData);
724                  for (UInt_t ui = 0; ui < 6; ui++){
725                    if (ui>0){
726                      if (L_QQ_Q_l_upper->time[ui]>L_QQ_Q_l_upper->time[0]){
727                        if (dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[ui]*1000-DeltaOBT*1000))<atime){
728                          upperqtime = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[ui]*1000-DeltaOBT*1000));
729                          orbits.getPosition((double) (upperqtime - gltle->GetFromTime())/60., &eCi);
730                          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]);
731                        }else {
732                          lowerqtime = upperqtime;
733                          upperqtime = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[ui]*1000-DeltaOBT*1000));
734                          orbits.getPosition((double) (upperqtime - gltle->GetFromTime())/60., &eCi);
735                          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]);
736                          mcreen = j3;
737                          mctren = ik;
738                          if(fgh==0){
739                            CopyQ(L_QQ_Q_l_lower,L_QQ_Q_l_upper);
740                            CopyAng(RYPang_lower,RYPang_upper);
741                          }
742                          oi=ui;
743                          goto closethisloop;
744                        }
745                        fgh++;
746                        CopyQ(L_QQ_Q_l_lower,L_QQ_Q_l_upper);
747                        CopyAng(RYPang_lower,RYPang_upper);
748                      }
749                    }else{
750                      if (dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000))<atime){
751                        upperqtime = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
752                        orbits.getPosition((double) (upperqtime - gltle->GetFromTime())/60., &eCi);
753                        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]);
754                      }
755                      else {
756                        lowerqtime = upperqtime;
757                        upperqtime = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
758                        orbits.getPosition((double) (upperqtime - gltle->GetFromTime())/60., &eCi);
759                        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]);
760                        mcreen = j3;
761                        mctren = ik;
762                        if(fgh==0){
763                          CopyQ(L_QQ_Q_l_lower,L_QQ_Q_l_upper);
764                          CopyAng(RYPang_lower,RYPang_upper);
765                          lowerqtime = atime-1;
766                        }
767                        oi=ui;
768                        goto closethisloop;
769                        //_0 = true;
770                      }
771                      fgh++;
772                      CopyQ(L_QQ_Q_l_lower,L_QQ_Q_l_upper);
773                      CopyAng(RYPang_lower,RYPang_upper);
774                      //_0 = true;
775                    };
776                    //cin>>grib;
777                  };
778                };
779              };
780            };
781          };
782        closethisloop:
783          //
784          if ( debug ) printf(" I am There \n");
785          //
786          if (((atime>(UInt_t)upperqtime)||(atime<(UInt_t)lowerqtime)) && neventsm>0 ){
787            if ( debug ) printf(" I am there \n");
788            //
789            lowerqtime = upperqtime;
790            Long64_t maxloop = 100000000LL;
791            Long64_t mn = 0;
792            bool gh=false;
793            ooi=oi;
794            if ( verbose ) printf(" OrbitalInfoCore: sync with quaternions data upperqtime %u lowerqtime %u atime %u \n",(UInt_t)upperqtime,(UInt_t)lowerqtime,atime);
795            while (!gh){      
796              if ( mn > maxloop ){
797                if ( verbose ) printf(" OrbitalInfoCore: quaternions sync out of range! exiting\n");
798                gh = true;
799                neventsm = 0;
800              };
801              mn++;
802              if (oi<5) oi++;
803              else oi=0;
804              if (oi==0 && numrec > 0){
805                if ( debug ) printf(" mumble \n");
806                mcreen++;
807                if (mcreen == numrec){
808                  mctren++;
809                  mcreen = 0;
810                  l0trm->GetEntry(mctren);
811                  numrec = mcmdev->Records->GetEntries();
812                }
813                CopyQ(L_QQ_Q_l_lower,L_QQ_Q_l_upper);
814                CopyAng(RYPang_lower,RYPang_upper);
815                mcmdrc = (pamela::McmdRecord*)mcmdev->Records->At(mcreen);
816                if ((int)mcmdrc->ID1 == 226){
817                  L_QQ_Q_l_upper->fill(mcmdrc->McmdData);
818                  upperqtime = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
819                  if (upperqtime<lowerqtime){
820                    upperqtime=runinfo->RUNTRAILER_TIME;
821                    CopyQ(L_QQ_Q_l_upper,L_QQ_Q_l_lower);
822                    CopyAng(RYPang_upper,RYPang_lower);
823                  }else{
824                    orbits.getPosition((double) (upperqtime - gltle->GetFromTime())/60., &eCi);
825                    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]);
826                  }
827                  //              re--;
828                  gh=true;
829                }
830              }else{
831                if ((Int_t)L_QQ_Q_l_upper->time[oi]>(Int_t)L_QQ_Q_l_upper->time[0]){
832                  upperqtime = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000));
833                  orbits.getPosition((double) (upperqtime - gltle->GetFromTime())/60., &eCi);
834                  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[oi][0],L_QQ_Q_l_upper->quat[oi][1],L_QQ_Q_l_upper->quat[oi][2],L_QQ_Q_l_upper->quat[oi][3]);
835                  orbits.getPosition((double) (lowerqtime - gltle->GetFromTime())/60., &eCi);
836                  RYPang_lower->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[oi-1][0],L_QQ_Q_l_upper->quat[oi-1][1],L_QQ_Q_l_upper->quat[oi-1][2],L_QQ_Q_l_upper->quat[oi-1][3]);
837                  //              re--;
838                  gh=true;
839                };
840              };
841            };
842            if ( verbose ) printf(" OrbitalInfoCore: sync with quaternions data now we have upperqtime %u lowerqtime %u atime %u \n",(UInt_t)upperqtime,(UInt_t)lowerqtime,atime);
843          };
844          //
845          if ( debug ) printf(" I am THIS \n");
846          //
847          // Fill in quaternions and angles
848          //
849          if ((atime<=(UInt_t)upperqtime)&&(atime>=(UInt_t)lowerqtime)&& neventsm>0){      
850            if ( debug ) printf(" I am this \n");
851            UInt_t tut = holeq(lowerqtime, upperqtime, L_QQ_Q_l_lower, L_QQ_Q_l_upper, oi);
852            if (oi == 0){
853              if ((tut!=5)||(tut!=6)){
854                incli = (L_QQ_Q_l_upper->quat[0][0]-L_QQ_Q_l_lower->quat[ooi][0])/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_lower->time[ooi]*1000-DeltaOBT*1000)));
855                orbitalinfo->q0 =  incli*atime+L_QQ_Q_l_upper->quat[0][0]-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
856                incli =     (L_QQ_Q_l_upper->quat[0][1]-L_QQ_Q_l_lower->quat[ooi][1])/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_lower->time[ooi]*1000-DeltaOBT*1000)));
857                orbitalinfo->q1 =  incli*atime+L_QQ_Q_l_upper->quat[0][1]-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
858                incli = (L_QQ_Q_l_upper->quat[0][2]-L_QQ_Q_l_lower->quat[ooi][2])/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_lower->time[ooi]*1000-DeltaOBT*1000)));
859                orbitalinfo->q2 =  incli*atime+L_QQ_Q_l_upper->quat[0][2]-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
860                incli = (L_QQ_Q_l_upper->quat[0][3]-L_QQ_Q_l_lower->quat[ooi][3])/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_lower->time[ooi]*1000-DeltaOBT*1000)));
861                orbitalinfo->q3 =  incli*atime+L_QQ_Q_l_upper->quat[0][3]-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
862            
863                incli = (RYPang_upper->Tangazh-RYPang_lower->Tangazh)/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_lower->time[ooi]*1000-DeltaOBT*1000)));
864                orbitalinfo->theta =  incli*atime+RYPang_upper->Tangazh-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
865                incli = (RYPang_upper->Ryskanie-RYPang_lower->Ryskanie)/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_lower->time[ooi]*1000-DeltaOBT*1000)));
866                orbitalinfo->phi =  incli*atime+RYPang_upper->Ryskanie-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
867                incli = (RYPang_upper->Kren-RYPang_lower->Kren)/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_lower->time[ooi]*1000)));
868                orbitalinfo->etha =  incli*atime+RYPang_upper->Kren-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
869              }
870              if (tut==6){
871                if (fabs(RYPang_lower->Kren-RYPang_upper->Kren)<0.1){
872                  incli = (L_QQ_Q_l_upper->quat[0][0]-L_QQ_Q_l_lower->quat[ooi][0])/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_lower->time[ooi]*1000-DeltaOBT*1000)));
873                  orbitalinfo->q0 =  incli*atime+L_QQ_Q_l_upper->quat[0][0]-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
874                  incli =           (L_QQ_Q_l_upper->quat[0][1]-L_QQ_Q_l_lower->quat[ooi][1])/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_lower->time[ooi]*1000-DeltaOBT*1000)));
875                  orbitalinfo->q1 =  incli*atime+L_QQ_Q_l_upper->quat[0][1]-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
876                  incli = (L_QQ_Q_l_upper->quat[0][2]-L_QQ_Q_l_lower->quat[ooi][2])/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_lower->time[ooi]*1000-DeltaOBT*1000)));
877                  orbitalinfo->q2 =  incli*atime+L_QQ_Q_l_upper->quat[0][2]-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
878                  incli = (L_QQ_Q_l_upper->quat[0][3]-L_QQ_Q_l_lower->quat[ooi][3])/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_lower->time[ooi]*1000-DeltaOBT*1000)));
879                  orbitalinfo->q3 =  incli*atime+L_QQ_Q_l_upper->quat[0][3]-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
880            
881                  incli = (RYPang_upper->Tangazh-RYPang_lower->Tangazh)/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_lower->time[ooi]*1000-DeltaOBT*1000)));
882                  orbitalinfo->theta =  incli*atime+RYPang_upper->Tangazh-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
883                  incli = (RYPang_upper->Ryskanie-RYPang_lower->Ryskanie)/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_lower->time[ooi]*1000-DeltaOBT*1000)));
884                  orbitalinfo->phi =  incli*atime+RYPang_upper->Ryskanie-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
885                  //cout<<"upper = "<<RYPang_upper->Ryskanie<<" lower = "<<RYPang_lower->Ryskanie<<" timeupper[0] = "<<L_QQ_Q_l_upper->time[0]-5500000<<" timelower["<<ooi<<"] = "<<L_QQ_Q_l_lower->time[ooi]-5500000<<" Ryscanie = "<<orbitalinfo->phi<<" incli = "<<incli<<" upper-lower = "<<RYPang_upper->Ryskanie-RYPang_lower->Ryskanie<<" Dtime = "<<dbtime->DBabsTime((UInt_t)L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000)<<"-"<<dbtime->DBabsTime((UInt_t)L_QQ_Q_l_lower->time[ooi]*1000-DeltaOBT*1000)<<" atime = "<<atime<<"\n";
886                  //cin>>grib;
887                  incli = (RYPang_upper->Kren-RYPang_lower->Kren)/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_lower->time[ooi]*1000-DeltaOBT*1000)));
888                  orbitalinfo->etha =  incli*atime+RYPang_upper->Kren-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
889                }
890              }
891            } else {
892              if((tut!=6)||(tut!=7)||(tut!=9)){
893                incli = (L_QQ_Q_l_upper->quat[oi][0]-L_QQ_Q_l_upper->quat[oi-1][0])/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi-1]*1000-DeltaOBT*1000)));
894                orbitalinfo->q0 =  incli*atime+L_QQ_Q_l_upper->quat[oi][0]-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000));
895                incli = (L_QQ_Q_l_upper->quat[oi][1]-L_QQ_Q_l_upper->quat[oi-1][1])/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi-1]*1000-DeltaOBT*1000)));
896                orbitalinfo->q1 =  incli*atime+L_QQ_Q_l_upper->quat[oi][1]-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000));
897                incli = (L_QQ_Q_l_upper->quat[oi][2]-L_QQ_Q_l_upper->quat[oi-1][2])/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi-1]*1000-DeltaOBT*1000)));
898                orbitalinfo->q2 =  incli*atime+L_QQ_Q_l_upper->quat[oi][2]-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000));
899                incli = (L_QQ_Q_l_upper->quat[oi][3]-L_QQ_Q_l_upper->quat[oi-1][3])/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi-1]*1000-DeltaOBT*1000)));
900                orbitalinfo->q3 =  incli*atime+L_QQ_Q_l_upper->quat[oi][3]-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000));
901            
902                incli = (RYPang_upper->Tangazh-RYPang_lower->Tangazh)/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi-1]*1000-DeltaOBT*1000)));
903                orbitalinfo->theta =  incli*atime+RYPang_upper->Tangazh-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000));
904                incli = (RYPang_upper->Ryskanie-RYPang_lower->Ryskanie)/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi-1]*1000-DeltaOBT*1000)));
905                orbitalinfo->phi =  incli*atime+RYPang_upper->Ryskanie-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000));
906                //cout<<"upper = "<<RYPang_upper->Ryskanie<<" lower = "<<RYPang_lower->Ryskanie<<" timeupper["<<oi<<"] = "<<L_QQ_Q_l_upper->time[oi]-5500000<<" timelower["<<oi-1<<"] = "<<L_QQ_Q_l_lower->time[oi-1]-5500000<<" Ryscanie = "<<orbitalinfo->phi<<" incli = "<<incli<<" upper-lower = "<<RYPang_upper->Ryskanie-RYPang_lower->Ryskanie<<" Dtime = "<<dbtime->DBabsTime((UInt_t)L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000)<<"-"<<dbtime->DBabsTime((UInt_t)L_QQ_Q_l_lower->time[oi-1]*1000-DeltaOBT*1000)<<" atime = "<<atime<<"\n";
907                //cin>>grib;
908                incli = (RYPang_upper->Kren-RYPang_lower->Kren)/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi-1]*1000-DeltaOBT*1000)));
909                orbitalinfo->etha =  incli*atime+RYPang_upper->Kren-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000));
910              }
911              if (tut==6){
912                if (fabs(RYPang_lower->Kren-RYPang_upper->Kren)<0.1){
913                  incli = (L_QQ_Q_l_upper->quat[oi][0]-L_QQ_Q_l_upper->quat[oi-1][0])/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi-1]*1000-DeltaOBT*1000)));
914                  orbitalinfo->q0 =  incli*atime+L_QQ_Q_l_upper->quat[oi][0]-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000));
915                  incli = (L_QQ_Q_l_upper->quat[oi][1]-L_QQ_Q_l_upper->quat[oi-1][1])/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi-1]*1000-DeltaOBT*1000)));
916                  orbitalinfo->q1 =  incli*atime+L_QQ_Q_l_upper->quat[oi][1]-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000));
917                  incli = (L_QQ_Q_l_upper->quat[oi][2]-L_QQ_Q_l_upper->quat[oi-1][2])/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi-1]*1000-DeltaOBT*1000)));
918                  orbitalinfo->q2 =  incli*atime+L_QQ_Q_l_upper->quat[oi][2]-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000));
919                  incli = (L_QQ_Q_l_upper->quat[oi][3]-L_QQ_Q_l_upper->quat[oi-1][3])/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi-1]*1000-DeltaOBT*1000)));
920                  orbitalinfo->q3 =  incli*atime+L_QQ_Q_l_upper->quat[oi][3]-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000));
921            
922                  incli = (RYPang_upper->Tangazh-RYPang_lower->Tangazh)/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi-1]*1000-DeltaOBT*1000)));
923                  orbitalinfo->theta =  incli*atime+RYPang_upper->Tangazh-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000));
924                  incli = (RYPang_upper->Ryskanie-RYPang_lower->Ryskanie)/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi-1]*1000-DeltaOBT*1000)));
925                  orbitalinfo->phi =  incli*atime+RYPang_upper->Ryskanie-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000));
926                  //cout<<"upper = "<<RYPang_upper->Ryskanie<<" lower = "<<RYPang_lower->Ryskanie<<" timeupper["<<oi<<"] = "<<L_QQ_Q_l_upper->time[oi]-5500000<<" timelower["<<oi-1<<"] = "<<L_QQ_Q_l_lower->time[oi-1]-5500000<<" Ryscanie = "<<orbitalinfo->phi<<" incli = "<<incli<<" upper-lower = "<<RYPang_upper->Ryskanie-RYPang_lower->Ryskanie<<" Dtime = "<<dbtime->DBabsTime((UInt_t)L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000)<<"-"<<dbtime->DBabsTime((UInt_t)L_QQ_Q_l_lower->time[oi-1]*1000-DeltaOBT*1000)<<" atime = "<<atime<<"\n";
927                  //cin>>grib;
928                  incli = (RYPang_upper->Kren-RYPang_lower->Kren)/(dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000))-dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi-1]*1000-DeltaOBT*1000)));
929                  orbitalinfo->etha =  incli*atime+RYPang_upper->Kren-incli*dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000));
930                }
931              }            
932            }
933            //
934            orbitalinfo->mode = holeq(lowerqtime, upperqtime, L_QQ_Q_l_lower, L_QQ_Q_l_upper, oi);
935            //
936          } else {
937            if ( debug ) printf(" ops no incl! \n");
938            orbitalinfo->mode = 10;
939          };
940          //
941          // ops no inclination information
942          //
943          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 ){
944            orbitalinfo->mode = 10;
945            orbitalinfo->q0 = -1000.;
946            orbitalinfo->q1 = -1000.;
947            orbitalinfo->q2 = -1000.;
948            orbitalinfo->q3 = -1000.;
949            orbitalinfo->etha = -1000.;
950            orbitalinfo->phi = -1000.;
951            orbitalinfo->theta = -1000.;
952          };
953          //
954          // #########################################################################################################################  
955          //
956          // fill orbital positions
957          //        
958          // Build coordinates in the right range.  We want to convert,
959          // longitude from (0, 2*pi) to (-180deg, 180deg).  Altitude is
960          // in meters.
961          lon = (coo.m_Lon > M_PI) ? rad2deg(coo.m_Lon - 2*M_PI) : rad2deg(coo.m_Lon);
962          lat = rad2deg(coo.m_Lat);
963          alt = coo.m_Alt;
964          //
965          if( lon<180 && lon>-180 && lat<90 && lat>-90 && alt>0 ){  
966            //      
967            orbitalinfo->lon = lon;
968            orbitalinfo->lat = lat;
969            orbitalinfo->alt = alt ;
970            //
971            // compute mag field components and L shell.
972            //
973            feldg_(&lat, &lon, &alt, &bnorth, &beast, &bdown, &babs);
974            shellg_(&lat, &lon, &alt, &dimo, &xl, &icode, &bab1);
975            findb0_(&stps, &bdel, &value, &bequ, &rr0);
976            //
977            orbitalinfo->Bnorth = bnorth;
978            orbitalinfo->Beast = beast;
979            orbitalinfo->Bdown = bdown;
980            orbitalinfo->Babs = babs;
981            orbitalinfo->BB0 = babs/bequ;
982            orbitalinfo->L = xl;      
983            // Set Stormer vertical cutoff using L shell.
984            orbitalinfo->cutoffsvl = 14.9/(xl*xl);
985            //
986          };      
987          //
988          // pitch angles
989          //
990          if ( orbitalinfo->mode != 10 && orbitalinfo->mode != 5 && orbitalinfo->mode !=7 && orbitalinfo->mode != 9 ){
991            //
992            Float_t Bx = -orbitalinfo->Bdown;                       //don't need for PamExp ExpOnly for all geography areas
993            Float_t By = orbitalinfo->Beast;                        //don't need for PamExp ExpOnly for all geography areas
994            Float_t Bz = orbitalinfo->Bnorth;                       //don't need for PamExp ExpOnly for all geography areas
995            //
996            TMatrixD Fij = PO->ECItoGreenwich(PO->QuatoECI(orbitalinfo->q0,orbitalinfo->q1,orbitalinfo->q2,orbitalinfo->q3),orbitalinfo->absTime);
997            TMatrixD Dij = PO->GreenwichtoGEO(orbitalinfo->lat,orbitalinfo->lon,Fij);
998            TMatrixD Iij = PO->ColPermutation(Dij);
999            //
1000            orbitalinfo->Iij.ResizeTo(Iij);
1001            orbitalinfo->Iij = Iij;
1002            //
1003            A1 = Iij(0,2);
1004            A2 = Iij(1,2);
1005            A3 = Iij(2,2);
1006            //      
1007            //      orbitalinfo->pamzenitangle = (Float_t)PO->GetPitchAngle(1,0,0,A1,A2,A3);                        // Angle between zenit and Pamela's main axiz
1008            //      orbitalinfo->pamBangle = (Float_t)PO->GetPitchAngle(A1,A2,A3,Bx,By,Bz);                 // Angle between Pamela's main axiz and B
1009            //
1010            if ( !standalone && tof->ntrk() > 0 ){
1011              //
1012              Int_t nn = 0;
1013              for(Int_t nt=0; nt < tof->ntrk(); nt++){  
1014                //
1015                ToFTrkVar *ptt = tof->GetToFTrkVar(nt);
1016                Double_t E11x = ptt->xtr_tof[0]; // tr->x[0];
1017                Double_t E11y = ptt->ytr_tof[0]; //tr->y[0];
1018                Double_t E11z = zin[0];
1019                Double_t E22x = ptt->xtr_tof[3];//tr->x[3];
1020                Double_t E22y = ptt->ytr_tof[3];//tr->y[3];
1021                Double_t E22z = zin[3];
1022                if ( (E11x < 100. && E11y < 100. && E22x < 100. && E22y < 100.) || ptt->trkseqno != -1  ){
1023                  Double_t norm = sqrt(pow(E22x-E11x,2)+pow(E22y-E11y,2)+pow(E22z-E11z,2));
1024                  Double_t MyAzim = TMath::RadToDeg()*atan(TMath::Abs(E22y-E11y)/TMath::Abs(E22x-E11x));
1025                  if(E22x-E11x>=0 && E22y-E11y <0) MyAzim =  360. - MyAzim;
1026                  if(E22x-E11x>=0 && E22y-E11y >=0) MyAzim = MyAzim;
1027                  if(E22x-E11x<0 && E22y-E11y >0) MyAzim = 180. - MyAzim;
1028                  if(E22x-E11x<0 && E22y-E11y <0) MyAzim = 180. + MyAzim;
1029                  Px = (E22x-E11x)/norm;
1030                  Py = (E22y-E11y)/norm;
1031                  Pz = (E22z-E11z)/norm;
1032                  //
1033                  t_orb->trkseqno = ptt->trkseqno;
1034                  //
1035                  TMatrixD Eij = PO->PamelatoGEO(Iij,Px,Py,Pz);
1036                  t_orb->Eij.ResizeTo(Eij);
1037                  t_orb->Eij = Eij;
1038                  //
1039                  TMatrixD Sij = PO->PamelatoGEO(Fij,Px,Py,Pz);
1040                  t_orb->Sij.ResizeTo(Sij);
1041                  t_orb->Sij = Sij;
1042                  //            
1043                  t_orb->pitch = (Float_t)PO->GetPitchAngle(Eij(0,0),Eij(1,0),Eij(2,0),Bx,By,Bz);
1044                  //
1045                  //
1046                  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);
1047                  //
1048                  t_orb->cutoff = 59.3/(pow(orbitalinfo->L,2)*pow((1+sqrt(1-pow(orbitalinfo->L,-3/2)*cos(omega))),2));
1049                  //
1050                  if ( debug ) printf(" orbitalinfo->cutoffsvl %f vitaly %f \n",orbitalinfo->cutoffsvl,t_orb->cutoff);
1051                  //
1052                  new(tor[nn]) OrbitalInfoTrkVar(*t_orb);
1053                  nn++;
1054                  //
1055                  t_orb->Clear();
1056                  //
1057                };
1058                //
1059              };
1060            } else {
1061              if ( debug ) printf(" mmm... mode %u standalone %i ntrk %i \n",orbitalinfo->mode,standalone,tof->ntrk());
1062            };
1063            //
1064          } else {
1065            if ( !standalone && tof->ntrk() > 0 ){
1066              //
1067              Int_t nn = 0;
1068              for(Int_t nt=0; nt < tof->ntrk(); nt++){  
1069                //
1070                ToFTrkVar *ptt = tof->GetToFTrkVar(nt);
1071                if ( ptt->trkseqno != -1  ){
1072                  //
1073                  t_orb->trkseqno = ptt->trkseqno;
1074                  //
1075                  t_orb->Eij = 0;  
1076                  //
1077                  t_orb->Sij = 0;
1078                  //            
1079                  t_orb->pitch = -1000.;
1080                  //
1081                  t_orb->cutoff = -1000.;
1082                  //
1083                  new(tor[nn]) OrbitalInfoTrkVar(*t_orb);
1084                  nn++;
1085                  //
1086                  t_orb->Clear();
1087                  //
1088                };
1089                //
1090              };    
1091            };
1092          };
1093          //
1094          // Fill the class
1095          //
1096        OrbitalInfotr->Fill();        OrbitalInfotr->Fill();
       //      
1097        //        //
1098      jumpev:        delete t_orb;
       debug = false;  
1099        //        //
1100      };      }; // loop over the events in the run
1101      //      //
1102      // Here you may want to clear some variables before processing another run        // Here you may want to clear some variables before processing another run  
1103      //      //
1104      ei = 0;      delete dbtime;
1105        if ( L_QQ_Q_l_upper ) delete L_QQ_Q_l_upper;
1106        if ( L_QQ_Q_l_lower ) delete L_QQ_Q_l_lower;
1107        if ( RYPang_upper ) delete RYPang_upper;
1108        if ( RYPang_lower ) delete RYPang_lower;
1109    }; // process all the runs    }; // process all the runs
1110    //    
1111    if (verbose) printf("\n Finished processing data \n");    if (verbose) printf("\n Finished processing data \n");
1112    //    //
1113   closeandexit:   closeandexit:
# Line 492  int OrbitalInfoCore(ULong64_t run, TFile Line 1130  int OrbitalInfoCore(ULong64_t run, TFile
1130          // copy orbitalinfoclone to OrbitalInfo          // copy orbitalinfoclone to OrbitalInfo
1131          //          //
1132          orbitalinfo->Clear();          orbitalinfo->Clear();
1133          //      orbitalinfo = new OrbitalInfo();          //
1134          memcpy(&orbitalinfo,&orbitalinfoclone,sizeof(orbitalinfoclone));          memcpy(&orbitalinfo,&orbitalinfoclone,sizeof(orbitalinfoclone));
1135          //          //
1136          // Fill entry in the new tree          // Fill entry in the new tree
# Line 507  int OrbitalInfoCore(ULong64_t run, TFile Line 1145  int OrbitalInfoCore(ULong64_t run, TFile
1145    //    //
1146    if ( l0File ) l0File->Close();    if ( l0File ) l0File->Close();
1147    if ( tempfile ) tempfile->Close();                if ( tempfile ) tempfile->Close();            
1148    gSystem->Unlink(tempname.str().c_str());    if ( myfold ) gSystem->Unlink(tempname.str().c_str());
1149      //
   //if ( code < 0 ) printf("\n OrbitalInfo - ERROR: an error occurred, try to save anyway...\n");  
   //printf("\n Writing and closing rootple\n");  
1150    if ( runinfo ) runinfo->Close();        if ( runinfo ) runinfo->Close();    
1151    if ( OrbitalInfotr ) OrbitalInfotr->SetName("OrbitalInfo");        if ( OrbitalInfotr ) OrbitalInfotr->SetName("OrbitalInfo");    
1152      if ( tof ) tof->Delete();
1153      if ( ttof ) ttof->Delete();
1154      //
1155    if ( file ){    if ( file ){
1156      file->cd();      file->cd();
1157      file->Write();      file->Write();
1158    };    };
1159    //    //
1160    gSystem->Unlink(OrbitalInfofolder.str().c_str());    if ( myfold ) gSystem->Unlink(OrbitalInfofolder.str().c_str());
1161    //    //
1162    // the end    // the end
1163    //    //
1164      if ( dbc ){
1165        dbc->Close();
1166        delete dbc;
1167      };
1168    if (verbose) printf("\n Exiting...\n");    if (verbose) printf("\n Exiting...\n");
1169    if(OrbitalInfotr)OrbitalInfotr->Delete();    if(OrbitalInfotr)OrbitalInfotr->Delete();
1170      //
1171      if ( PO ) delete PO;
1172      if ( orbitalinfo ) delete orbitalinfo;
1173      if ( orbitalinfoclone ) delete orbitalinfoclone;
1174      if ( glroot ) delete glroot;
1175      if ( runinfo ) delete runinfo;
1176      //
1177    if(code < 0)  throw code;    if(code < 0)  throw code;
1178    return(code);    return(code);
1179  }  }
1180    
1181    
1182    //
1183    // Returns the cCoordGeo structure holding the geographical
1184    // coordinates for the event (see sgp4.h).
1185    //
1186    // atime is the abstime of the event in UTC unix time.
1187    // tletime is the time of the tle in UTC unix time.
1188    // tle is the previous and nearest tle (compared to atime).
1189    cCoordGeo getCoo(UInt_t atime, UInt_t tletime, cTle *tle)
1190    {
1191      cEci eci;
1192      cOrbit orbit(*tle);
1193      orbit.getPosition((double) (atime - tletime)/60., &eci);
1194      
1195      return eci.toGeo();
1196    }
1197    
1198    // function of copyng of quatrnions classes
1199    
1200    void CopyQ(Quaternions *Q1, Quaternions *Q2){
1201      for(UInt_t i = 0; i < 6; i++){
1202        Q1->time[i]=Q2->time[i];
1203        for (UInt_t j = 0; j < 4; j++)Q1->quat[i][j]=Q2->quat[i][j];
1204      }
1205      return;
1206    }
1207    
1208    // functions of copyng InclinationInfo classes
1209    
1210    void CopyAng(InclinationInfo *A1, InclinationInfo *A2){
1211      A1->Tangazh = A2->Tangazh;
1212      A1->Ryskanie = A2->Ryskanie;
1213      A1->Kren = A2->Kren;
1214      return;
1215    }
1216    
1217    UInt_t holeq(Double_t lower,Double_t upper,Quaternions *Qlower, Quaternions *Qupper, UInt_t f){
1218      
1219      UInt_t hole = 10;
1220      bool R10l = false;     // Sign of R10 mode in lower quaternions array
1221      bool R10u = false;     // Sign of R10 mode in upper quaternions array
1222      bool insm = false;     // Sign that we inside quaternions array
1223      bool mxtml = false;    // Sign of mixt mode in lower quaternions array
1224      bool mxtmu = false;    // Sign of mixt mode in upper quaternions array
1225      bool npasm = false;     // Sign of normall pass between R10 and non R10 or between non R10 and R10
1226      UInt_t NCQl = 6;       // Number of correct quaternions in lower array
1227      UInt_t NCQu = 6;       // Number of correct quaternions in upper array
1228      if (f>0){
1229        insm = true;
1230        if(Qupper->time[f]-Qupper->time[f-1]==30) R10u = false;
1231        if(Qupper->time[f]-Qupper->time[f-1]<1) R10u = true;
1232      }else{
1233        insm = false;
1234        if((Qlower->time[5]-Qlower->time[0]<2)&&(Qlower->time[1]-Qlower->time[0]<2)) R10l = true;
1235        if((Qupper->time[5]-Qupper->time[0]<2)&&(Qupper->time[1]-Qupper->time[0]<2)) R10u = true;
1236        if((Qlower->time[5]-Qlower->time[0]==150)&&(Qlower->time[1]-Qlower->time[0]==30)) R10l = false;
1237        if((Qupper->time[5]-Qupper->time[0]==150)&&(Qupper->time[1]-Qupper->time[0]==30)) R10u = false;
1238        if((Qlower->time[5]-Qlower->time[0]<2)&&(Qlower->time[1]-Qlower->time[0]==30)){
1239          mxtml = true;
1240          for(UInt_t i = 1; i < 6; i++){
1241            if(Qlower->time[i]-Qlower->time[0]==30*i) NCQl=i;
1242          }
1243        }
1244        if((Qupper->time[5]-Qupper->time[0]<2)&&(Qupper->time[1]-Qupper->time[0]==30)){
1245          mxtmu = true;
1246          for(UInt_t i = 1; i < 6; i++){
1247            if(Qupper->time[i]-Qupper->time[0]==30*i) NCQu=i;
1248          }
1249        }
1250      }
1251      
1252      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;
1253      
1254      
1255      if (R10u&&insm) hole=0; // best event R10
1256      if ((upper-lower<=5)&&(!insm)&&R10l&&R10u) hole = 1; // when first of 6 quaternions in array is correct
1257      if (((!R10u)&&insm)||((!insm)&&(!R10u)&&(!R10l)&&((upper-lower==210+(6-NCQl)*30)||(upper-lower==30)))) hole = 2; //non R10
1258      if (npasm&&(!insm)&&((R10l&&!R10u)||(R10u&&!R10l))) hole = 3; //normall pass from R10 to non R10 or from non R10 to R10
1259      if ((!npasm)&&(upper-lower<=300)&&(!insm)&&((R10l&&!R10u)||(R10u&&!R10l))) hole = 4; // eliminable hole between R10 and non R10 or between non R10 and R10
1260      if ((upper-lower>=300)&&(!insm)&&((R10l&&!R10u)||(R10u&&!R10l))) hole = 5; //uneliminable hole between R10 and non R10 or between non R10 and R10
1261      if ((upper-lower>5)&&(upper-lower<=300)&&R10u&&R10l) hole = 6; // eliminable hole inside R10
1262      if ((upper-lower>300)&&R10u&&R10l) hole = 7; //uneliminable hole inside R10
1263      if ((upper-lower>210)&&(upper-lower<=1200)&&(!R10u)&&(!R10l)) hole = 8; //eliminable hole inside non R10
1264      if ((upper-lower>1200)&&!R10u&&!R10l) hole = 9; // uneliminable hole inside non R10
1265      return hole;
1266    }
1267    
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.35
  ViewVC Help
Powered by ViewVC 1.1.23