/[PAMELA software]/DarthVader/OrbitalInfo/src/OrbitalInfoCore.cpp
ViewVC logotype

Annotation of /DarthVader/OrbitalInfo/src/OrbitalInfoCore.cpp

Parent Directory Parent Directory | Revision Log Revision Log

Revision 1.93 - (hide annotations) (download)
Thu Dec 27 13:12:28 2018 UTC (6 years, 2 months ago) by mayorov
Branch: MAIN
CVS Tags: HEAD
Changes since 1.92: +879 -864 lines 
Fixed bug with EXT track variables initialization
1 mocchiut 1.1 // C/C++ headers
2     //
3     #include <fstream>
4     #include <string.h>
5     #include <iostream>
6     #include <cstring>
7     #include <stdio.h>
8     //
9     // ROOT headers
10     //
11 mocchiut 1.44 //#include <TCanvas.h>
12 pam-mep 1.60 #include <TH2F.h> //for test only. Vitaly.
13     #include <TVector3.h>
14 mocchiut 1.44 //#include <TF1.h>
15    
16 mocchiut 1.1 #include <TTree.h>
17     #include <TClassEdit.h>
18     #include <TObject.h>
19     #include <TList.h>
20 mocchiut 1.5 #include <TArrayI.h>
21 mocchiut 1.1 #include <TSystem.h>
22     #include <TSystemDirectory.h>
23     #include <TString.h>
24     #include <TFile.h>
25     #include <TClass.h>
26     #include <TSQLServer.h>
27     #include <TSQLRow.h>
28     #include <TSQLResult.h>
29 mocchiut 1.81 #include <TObjectTable.h>
30 mocchiut 1.1 //
31 mocchiut 1.8 // RunInfo header
32     //
33     #include <RunInfo.h>
34     #include <GLTables.h>
35     //
36 mocchiut 1.1 // YODA headers
37     //
38     #include <PamelaRun.h>
39     #include <PscuHeader.h>
40     #include <PscuEvent.h>
41     #include <EventHeader.h>
42 mocchiut 1.15 #include <mcmd/McmdEvent.h>
43     #include <mcmd/McmdRecord.h>
44 mocchiut 1.1 //
45     // This program headers
46     //
47     #include <OrbitalInfo.h>
48 mocchiut 1.7 #include <OrbitalInfoVerl2.h>
49 mocchiut 1.1 #include <OrbitalInfoCore.h>
50 mocchiut 1.15 #include <InclinationInfo.h>
51 mocchiut 1.1
52 pam-mep 1.69 //
53     // Tracker and ToF classes headers and definitions
54     //
55     #include <ToFLevel2.h>
56     #include <TrkLevel2.h>
57 mocchiut 1.73 #include <ExtTrack.h> // new tracking code
58 mocchiut 1.40
59 mocchiut 1.1 using namespace std;
60    
61     //
62     // CORE ROUTINE
63     //
64     //
65 mocchiut 1.81 int OrbitalInfoCore(UInt_t run, TFile *file, GL_TABLES *glt, Int_t OrbitalInfoargc, char *OrbitalInfoargv[]){
66 mocchiut 1.15 //
67 mocchiut 1.1 Int_t i = 0;
68 mocchiut 1.25 TString host = glt->CGetHost();
69     TString user = glt->CGetUser();
70     TString psw = glt->CGetPsw();
71     TSQLServer *dbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data());
72 mocchiut 1.1 //
73 mocchiut 1.26 stringstream myquery;
74     myquery.str("");
75 mocchiut 1.79 myquery << "SET time_zone='+0:00';";
76 mocchiut 1.59 delete dbc->Query(myquery.str().c_str());
77 mocchiut 1.79 delete dbc->Query("SET sql_mode = 'NO_UNSIGNED_SUBTRACTION';");
78 mocchiut 1.26 //
79 mocchiut 1.16 TString processFolder = Form("OrbitalInfoFolder_%u",run);
80 mocchiut 1.1 //
81     // Set these to true to have a very verbose output.
82     //
83     Bool_t debug = false;
84     //
85     Bool_t verbose = false;
86 mocchiut 1.30 //
87 mocchiut 1.31 Bool_t standalone = false;
88 mocchiut 1.30 //
89 mocchiut 1.1 if ( OrbitalInfoargc > 0 ){
90     i = 0;
91     while ( i < OrbitalInfoargc ){
92     if ( !strcmp(OrbitalInfoargv[i],"-processFolder") ) {
93 mayorov 1.93 if ( OrbitalInfoargc < i+1 ){
94     throw -3;
95     }
96     processFolder = (TString)OrbitalInfoargv[i+1];
97     i++;
98 mocchiut 1.81 }
99 mocchiut 1.1 if ( (!strcmp(OrbitalInfoargv[i],"--debug")) || (!strcmp(OrbitalInfoargv[i],"-g")) ) {
100 mayorov 1.93 verbose = true;
101     debug = true;
102 mocchiut 1.81 }
103 mocchiut 1.1 if ( (!strcmp(OrbitalInfoargv[i],"--verbose")) || (!strcmp(OrbitalInfoargv[i],"-v")) ) {
104 mayorov 1.93 verbose = true;
105 mocchiut 1.81 }
106 mocchiut 1.30 if ( (!strcmp(OrbitalInfoargv[i],"--standalone")) ) {
107 mayorov 1.93 standalone = true;
108 mocchiut 1.81 }
109 mocchiut 1.30 if ( (!strcmp(OrbitalInfoargv[i],"--calculate-pitch")) ) {
110 mayorov 1.93 standalone = false;
111 mocchiut 1.81 }
112 mocchiut 1.1 i++;
113 mocchiut 1.81 }
114     }
115     if ( debug ){
116     printf("START\n");
117     gObjectTable->Print();
118     }
119 mocchiut 1.1 //
120     const char* outDir = gSystem->DirName(gSystem->DirName(file->GetPath()));
121     //
122     TTree *OrbitalInfotr = 0;
123 mocchiut 1.5 UInt_t nevents = 0;
124 mocchiut 1.15 UInt_t neventsm = 0;
125 mocchiut 1.1 //
126     // variables needed to reprocess data
127     //
128 mocchiut 1.6 Long64_t maxsize = 10000000000LL;
129     TTree::SetMaxTreeSize(maxsize);
130     //
131 mocchiut 1.1 TString OrbitalInfoversion;
132     ItoRunInfo *runinfo = 0;
133 mocchiut 1.5 TArrayI *runlist = 0;
134 mocchiut 1.1 TTree *OrbitalInfotrclone = 0;
135     Bool_t reproc = false;
136     Bool_t reprocall = false;
137 mocchiut 1.66 Bool_t igrfloaded = false;
138 mocchiut 1.1 UInt_t nobefrun = 0;
139     UInt_t noaftrun = 0;
140     UInt_t numbofrun = 0;
141     stringstream ftmpname;
142     TString fname;
143 mocchiut 1.5 UInt_t totfileentries = 0;
144 mocchiut 1.15 UInt_t idRun = 0;
145 pam-fi 1.53 UInt_t anni5 = 60 * 60 * 24 * 365 * 5 ;//1576800
146 mocchiut 1.15 //
147     // My variables. Vitaly.
148     //
149 mocchiut 1.44 // UInt_t oi = 0;
150 mocchiut 1.15 Int_t tmpSize = 0;
151 mocchiut 1.1 //
152     // variables needed to handle error signals
153     //
154     Int_t code = 0;
155     Int_t sgnl;
156     //
157     // OrbitalInfo classes
158     //
159     OrbitalInfo *orbitalinfo = new OrbitalInfo();
160     OrbitalInfo *orbitalinfoclone = new OrbitalInfo();
161 mocchiut 1.44
162 mocchiut 1.1 //
163     // define variables for opening and reading level0 file
164     //
165     TFile *l0File = 0;
166     TTree *l0tr = 0;
167 mocchiut 1.37 // TTree *l0trm = 0;
168     TChain *ch = 0;
169 mocchiut 1.2 // EM: open also header branch
170     TBranch *l0head = 0;
171     pamela::EventHeader *eh = 0;
172     pamela::PscuHeader *ph = 0;
173 mocchiut 1.15 pamela::McmdEvent *mcmdev = 0;
174     pamela::McmdRecord *mcmdrc = 0;
175 mocchiut 1.2 // end EM
176 mocchiut 1.15
177     // pamela::RunHeaderEvent *reh = new pamela::RunHeaderEvent;
178     // pamela::EventHeader *eH = new pamela::EventHeader;
179    
180 mocchiut 1.1 //
181     // Define other basic variables
182     //
183     UInt_t procev = 0;
184     stringstream file2;
185     stringstream file3;
186     stringstream qy;
187     Int_t totevent = 0;
188 mocchiut 1.5 UInt_t atime = 0;
189     UInt_t re = 0;
190 mocchiut 1.15 UInt_t ik = 0;
191 mocchiut 1.7
192     // Position
193     Float_t lon, lat, alt;
194    
195     //
196     // IGRF stuff
197     //
198 mocchiut 1.81 Float_t dimo = 0.0; // dipole moment (computed from dat files) // EM GCC 4.7
199 mocchiut 1.40 Float_t bnorth, beast, bdown, babs;
200     Float_t xl; // L value
201 mocchiut 1.85 Int_t icode; // code value for L accuracy (see fortran code)
202 mocchiut 1.40 Float_t bab1; // What's the difference with babs?
203     Float_t stps = 0.005; // step size for field line tracing
204     Float_t bdel = 0.01; // required accuracy
205     Float_t bequ; // equatorial b value (also called b_0)
206     Bool_t value = 0; // false if bequ is not the minimum b value
207     Float_t rr0; // equatorial radius normalized to earth radius
208 mocchiut 1.7
209 mocchiut 1.1 //
210     // Working filename
211     //
212     TString outputfile;
213     stringstream name;
214     name.str("");
215     name << outDir << "/";
216     //
217     // temporary file and folder
218     //
219     TFile *tempfile = 0;
220     TTree *tempOrbitalInfo = 0;
221     stringstream tempname;
222     stringstream OrbitalInfofolder;
223 mocchiut 1.23 Bool_t myfold = false;
224 mocchiut 1.1 tempname.str("");
225     tempname << outDir;
226     tempname << "/" << processFolder.Data();
227     OrbitalInfofolder.str("");
228     OrbitalInfofolder << tempname.str().c_str();
229     tempname << "/OrbitalInfotree_run";
230     tempname << run << ".root";
231 mocchiut 1.40 UInt_t totnorun = 0;
232 mocchiut 1.1 //
233     // DB classes
234     //
235     GL_ROOT *glroot = new GL_ROOT();
236 mocchiut 1.5 GL_TIMESYNC *dbtime = 0;
237 mocchiut 1.7 GL_TLE *gltle = new GL_TLE();
238 mocchiut 1.15 //
239     //Quaternions classes
240     //
241 mocchiut 1.57 Quaternions *L_QQ_Q_l_lower = 0;
242     InclinationInfo *RYPang_lower = 0;
243     Quaternions *L_QQ_Q_l_upper = 0;
244     InclinationInfo *RYPang_upper = 0;
245 mocchiut 1.15
246     cEci eCi;
247    
248 mocchiut 1.7 // Initialize fortran routines!!!
249 mocchiut 1.52 Int_t ltp1 = 0;
250 mocchiut 1.7 Int_t ltp2 = 0;
251 pam-mep 1.69 GL_PARAM *glparam0 = new GL_PARAM();
252 mocchiut 1.7 GL_PARAM *glparam = new GL_PARAM();
253 mocchiut 1.10 GL_PARAM *glparam2 = new GL_PARAM();
254 mocchiut 1.44
255 mocchiut 1.30 //
256 mocchiut 1.44 // Orientation variables. Vitaly
257 mocchiut 1.30 //
258 malakhov 1.74
259 mocchiut 1.30 UInt_t evfrom = 0;
260     UInt_t jumped = 0;
261     Int_t itr = -1;
262 mocchiut 1.64 // Double_t A1;
263     // Double_t A2;
264     // Double_t A3;
265 mocchiut 1.30 Double_t Px = 0;
266 mayorov 1.93 Double_t Py = 0;
267 mocchiut 1.30 Double_t Pz = 0;
268     TTree *ttof = 0;
269     ToFLevel2 *tof = new ToFLevel2();
270 pam-mep 1.69 TTree *ttrke = 0;
271     TrkLevel2 *trke = new TrkLevel2();
272 mocchiut 1.30 OrientationInfo *PO = new OrientationInfo();
273     Int_t nz = 6;
274     Float_t zin[6];
275     Int_t nevtofl2 = 0;
276 pam-mep 1.69 Int_t nevtrkl2 = 0;
277 mocchiut 1.44 if ( verbose ) cout<<"Reading quaternions external file"<<endl;
278     cout.setf(ios::fixed,ios::floatfield);
279     /******Reading recovered quaternions...*********/
280     vector<Double_t> recqtime;
281     vector<Float_t> recq0;
282     vector<Float_t> recq1;
283     vector<Float_t> recq2;
284     vector<Float_t> recq3;
285     Float_t Norm = 1;
286 mocchiut 1.81 recqtime.reserve(1500000);
287     recq0.reserve(1500000);
288     recq1.reserve(1500000);
289     recq2.reserve(1500000);
290     recq3.reserve(1500000);
291 mocchiut 1.73
292     vector<UInt_t> RTtime1;
293     vector<UInt_t> RTtime2;
294     vector<Double_t> RTbank1;
295     vector<Double_t> RTbank2;
296 malakhov 1.74 vector<Double_t> RTbpluto1;
297     vector<Double_t> RTbpluto2;
298 mocchiut 1.73 vector<Int_t> RTazim;
299 mocchiut 1.75 vector<UInt_t> RTstart;
300     vector<UInt_t> RTpluto2;
301     vector<UInt_t> RTpluto1;
302 mocchiut 1.73 vector<Int_t> RTerrq;
303 malakhov 1.82 vector<Int_t> RTqual;
304 mocchiut 1.81 RTtime1.reserve(200000);
305     RTtime2.reserve(200000);
306     RTbank1.reserve(200000);
307     RTbank2.reserve(200000);
308     RTbpluto1.reserve(200000);
309     RTbpluto2.reserve(200000);
310     RTazim.reserve(200000);
311     RTstart.reserve(200000);
312     RTpluto1.reserve(200000);
313     RTpluto2.reserve(200000);
314     RTerrq.reserve(200000);
315 malakhov 1.82 RTqual.reserve(200000);
316 mocchiut 1.73
317     TClonesArray *tcNucleiTrk = NULL;
318     TClonesArray *tcExtNucleiTrk = NULL;
319     TClonesArray *tcExtTrk = NULL;
320     TClonesArray *tcNucleiTof = NULL;
321     TClonesArray *tcExtNucleiTof = NULL;
322     TClonesArray *tcExtTof = NULL;
323     TClonesArray *torbNucleiTrk = NULL;
324     TClonesArray *torbExtNucleiTrk = NULL;
325     TClonesArray *torbExtTrk = NULL;
326     Bool_t hasNucleiTrk = false;
327     Bool_t hasExtNucleiTrk = false;
328     Bool_t hasExtTrk = false;
329     Bool_t hasNucleiTof = false;
330     Bool_t hasExtNucleiTof = false;
331     Bool_t hasExtTof = false;
332    
333     ifstream in;
334     ifstream an;
335     // ofstream mc;
336     // TString gr;
337     Int_t parerror2=0;
338    
339 pam-mep 1.69 Int_t parerror=glparam0->Query_GL_PARAM(1,303,dbc); // parameters stored in DB in GL_PRAM table
340 mocchiut 1.72 if ( verbose ) cout<<parerror<<"\t"<<(char*)(glparam0->PATH+glparam0->NAME).Data()<<endl;
341 mocchiut 1.44 if ( parerror<0 ) {
342     code = parerror;
343 mocchiut 1.73 goto closeandexit;
344 pam-mep 1.54 }
345 mocchiut 1.73 in.open((char*)(glparam0->PATH+glparam0->NAME).Data(),ios::in);
346 mocchiut 1.44 while(!in.eof()){
347     recqtime.resize(recqtime.size()+1);
348     Int_t sizee = recqtime.size();
349     recq0.resize(sizee);
350     recq1.resize(sizee);
351     recq2.resize(sizee);
352     recq3.resize(sizee);
353     in>>recqtime[sizee-1];
354     in>>recq0[sizee-1];
355     in>>recq1[sizee-1];
356     in>>recq2[sizee-1];
357     in>>recq3[sizee-1];
358     in>>Norm;
359 malakhov 1.82 /* CHECK RECOVERED QUATERNIONS PROBLEM
360     if(recqtime[sizee-1]>=1160987921.75 && recqtime[sizee-1]<=1160987932.00){
361     cout<<"We found it at start"<<"\t"<<recqtime[sizee-1]<<endl;
362     } */
363 mocchiut 1.44 }
364 pam-mep 1.69 in.close();
365 mocchiut 1.44 if ( verbose ) cout<<"We have read recovered data"<<endl;
366 pam-mep 1.69 if (debug) cout << "size of recovered quaterions data set is " << recqtime.size() << endl;
367 mocchiut 1.81 if ( debug ) printf(" RQ size %i RQ capacity %i \n",(int)recqtime.size(),(int)recqtime.capacity());
368    
369 pam-mep 1.69 if ( verbose ) cout<<"read Rotation Table"<<endl;
370 mocchiut 1.73
371     parerror2=glparam0->Query_GL_PARAM(1,305,dbc);
372 malakhov 1.74
373 mocchiut 1.71 if ( verbose ) cout<<parerror2<<"\t"<<(char*)(glparam0->PATH+glparam0->NAME).Data()<<endl;
374 mocchiut 1.73 if ( parerror2<0 ) {
375     code = parerror;
376     goto closeandexit;
377     }
378     an.open((char*)(glparam0->PATH+glparam0->NAME).Data(),ios::in);
379 pam-mep 1.69 while(!an.eof()){
380     RTtime1.resize(RTtime1.size()+1);
381     Int_t sizee = RTtime1.size();
382     RTbank1.resize(sizee+1);
383     RTazim.resize(sizee+1);
384     RTerrq.resize(sizee+1);
385 malakhov 1.74 RTstart.resize(sizee+1);
386     RTpluto1.resize(sizee+1);
387     RTbpluto1.resize(sizee+1);
388 malakhov 1.82 RTqual.resize(sizee+1);
389 pam-mep 1.69 an>>RTtime1[sizee-1];
390     an>>RTbank1[sizee-1];
391 malakhov 1.74 an>>RTstart[sizee-1];
392     an>>RTpluto1[sizee-1];
393     an>>RTbpluto1[sizee-1];
394 pam-mep 1.69 an>>RTazim[sizee-1];
395     an>>RTerrq[sizee-1];
396 malakhov 1.82 an>>RTqual[sizee-1];
397 pam-mep 1.69 if(sizee>1) {
398     RTtime2.resize(sizee+1);
399     RTbank2.resize(sizee+1);
400 malakhov 1.74 RTpluto2.resize(sizee+1);
401     RTbpluto2.resize(sizee+1);
402 pam-mep 1.69 RTtime2[sizee-2]=RTtime1[sizee-1];
403 malakhov 1.74 RTpluto2[sizee-2]=RTpluto1[sizee-1];
404 pam-mep 1.69 RTbank2[sizee-2]=RTbank1[sizee-1];
405 malakhov 1.74 RTbpluto2[sizee-2]=RTbpluto1[sizee-1];
406 pam-mep 1.69 }
407     }
408     an.close();
409     //cout<<"put some number here"<<endl;
410     //Int_t yupi;
411     //cin>>yupi;
412    
413     if ( verbose ) cout<<"We have read Rotation Table"<<endl;
414     //Geomagnetic coordinates calculations staff
415 mocchiut 1.81
416     if ( debug ) printf(" RT size %i RT capacity %i \n",(int)RTtime2.size(),(int)RTtime2.capacity());
417 pam-mep 1.69
418     GMtype_CoordGeodetic location;
419 mocchiut 1.70 // GMtype_CoordDipole GMlocation;
420 pam-mep 1.69 GMtype_Ellipsoid Ellip;
421     GMtype_Data G0, G1, H1;
422 mayorov 1.93
423     // { // this braces is necessary to avoid jump to label 'closeandexit' error // but it is wrong since the variable "igpath" will not exist outside. To overcome the "jump to label 'closeandexit' error" it is necessary to set the "igpath" before line 276
424 pam-mep 1.69 // TString igpath="/data03/Malakhov/pam9Malakhov/installed10/calib/orb-param/";
425     // }
426    
427     //cout << G0.element[0] << "\t" << G1.element[0] << "\t" << H1.element[0] << endl;
428     //cout << G0.element[5] << "\t" << G1.element[5] << "\t" << H1.element[5] << endl;
429    
430     GM_SetEllipsoid(&Ellip);
431 mocchiut 1.44
432 mocchiut 1.55 // IGRF stuff moved inside run loop!
433    
434 mocchiut 1.30 for (Int_t ip=0;ip<nz;ip++){
435     zin[ip] = tof->GetZTOF(tof->GetToFPlaneID(ip));
436     };
437     //
438     if ( !standalone ){
439     //
440 mocchiut 1.73 // Does it contain the Tracker and ToF trees?
441 mocchiut 1.30 //
442     ttof = (TTree*)file->Get("ToF");
443     if ( !ttof ) {
444     if ( verbose ) printf(" OrbitalInfo - ERROR: no tof tree\n");
445     code = -900;
446     goto closeandexit;
447 pam-mep 1.69 }
448 mocchiut 1.30 ttof->SetBranchAddress("ToFLevel2",&tof);
449     nevtofl2 = ttof->GetEntries();
450 pam-mep 1.69
451 mocchiut 1.73 //
452     // Look for extended tracking algorithm
453     //
454     if ( verbose ) printf("Look for extended and nuclei tracking algorithms in ToF\n");
455     // Nuclei tracking algorithm
456     Int_t checkAlgo = 0;
457     tcNucleiTof = new TClonesArray("ToFTrkVar");
458     checkAlgo = ttof->SetBranchAddress("TrackNuclei",&tcNucleiTof);
459     if ( !checkAlgo ){
460     if ( verbose ) printf(" Nuclei tracking algorithm ToF branch found! :D \n");
461     hasNucleiTof = true;
462     } else {
463     if ( verbose ) printf(" Nuclei tracking algorithm ToF branch not found :( !\n");
464     printf(" ok, this is not a problem (it depends on tracker settings) \n");
465     delete tcNucleiTof;
466 mocchiutti 1.77 tcNucleiTof=NULL; // 10RED reprocessing bug
467 mocchiut 1.73 }
468     // Nuclei tracking algorithm using calorimeter points
469     tcExtNucleiTof = new TClonesArray("ToFTrkVar");
470     checkAlgo = ttof->SetBranchAddress("RecoveredTrackNuclei",&tcExtNucleiTof);
471     if ( !checkAlgo ){
472     if ( verbose ) printf(" Recovered nuclei tracking algorithm ToF branch found! :D \n");
473     hasExtNucleiTof = true;
474     } else {
475     if ( verbose ) printf(" Recovered nuclei tracking algorithm ToF branch not found :( !\n");
476     printf(" ok, this is not a problem (it depends on tracker settings) \n");
477     delete tcExtNucleiTof;
478 mocchiutti 1.77 tcExtNucleiTof=NULL; // 10RED reprocessing bug
479 mocchiut 1.73 }
480     // Tracking algorithm using calorimeter points
481     tcExtTof = new TClonesArray("ToFTrkVar");
482     checkAlgo = ttof->SetBranchAddress("RecoveredTrack",&tcExtTof);
483     if ( !checkAlgo ){
484     if ( verbose ) printf(" Recovered track algorithm ToF branch found! :D \n");
485     hasExtTof = true;
486     } else {
487     if ( verbose ) printf(" Recovered track algorithm ToF branch not found :( !\n");
488     printf(" ok, this is not a problem (it depends on tracker settings) \n");
489     delete tcExtTof;
490 mocchiutti 1.77 tcExtTof=NULL; // 10RED reprocessing bug
491 mocchiut 1.73 }
492    
493 pam-mep 1.69 ttrke = (TTree*)file->Get("Tracker");
494     if ( !ttrke ) {
495     if ( verbose ) printf(" OrbitalInfo - ERROR: no trk tree\n");
496     code = -903;
497     goto closeandexit;
498     }
499     ttrke->SetBranchAddress("TrkLevel2",&trke);
500     nevtrkl2 = ttrke->GetEntries();
501 mocchiut 1.73
502     //
503     // Look for extended tracking algorithm
504     //
505     if ( verbose ) printf("Look for extended and nuclei tracking algorithms\n");
506     // Nuclei tracking algorithm
507     checkAlgo = 0;
508     tcNucleiTrk = new TClonesArray("TrkTrack");
509     checkAlgo = ttrke->SetBranchAddress("TrackNuclei",&tcNucleiTrk);
510     if ( !checkAlgo ){
511     if ( verbose ) printf(" Nuclei tracking algorithm branch found! :D \n");
512     hasNucleiTrk = true;
513     } else {
514     if ( verbose ) printf(" Nuclei tracking algorithm branch not found :( !\n");
515     printf(" ok, this is not a problem (it depends on tracker settings) \n");
516     delete tcNucleiTrk;
517 mocchiutti 1.77 tcNucleiTrk=NULL; // 10RED reprocessing bug
518 mocchiut 1.73 }
519     // Nuclei tracking algorithm using calorimeter points
520     tcExtNucleiTrk = new TClonesArray("ExtTrack");
521     checkAlgo = ttrke->SetBranchAddress("RecoveredTrackNuclei",&tcExtNucleiTrk);
522     if ( !checkAlgo ){
523     if ( verbose ) printf(" Recovered nuclei tracking algorithm branch found! :D \n");
524     hasExtNucleiTrk = true;
525     } else {
526     if ( verbose ) printf(" Recovered nuclei tracking algorithm branch not found :( !\n");
527     printf(" ok, this is not a problem (it depends on tracker settings) \n");
528     delete tcExtNucleiTrk;
529 mocchiutti 1.77 tcExtNucleiTrk=NULL; // 10RED reprocessing bug
530 mocchiut 1.73 }
531     // Tracking algorithm using calorimeter points
532     tcExtTrk = new TClonesArray("ExtTrack");
533     checkAlgo = ttrke->SetBranchAddress("RecoveredTrack",&tcExtTrk);
534     if ( !checkAlgo ){
535     if ( verbose ) printf(" Recovered track algorithm branch found! :D \n");
536     hasExtTrk = true;
537     } else {
538     if ( verbose ) printf(" Recovered track algorithm branch not found :( !\n");
539     printf(" ok, this is not a problem (it depends on tracker settings) \n");
540     delete tcExtTrk;
541 mocchiutti 1.77 tcExtTrk=NULL; // 10RED reprocessing bug
542 mocchiut 1.73 }
543    
544     if ( (hasNucleiTrk && !hasNucleiTof) || (!hasNucleiTrk && hasNucleiTof) ||
545     (hasExtNucleiTrk && !hasExtNucleiTof) || (!hasExtNucleiTrk && hasExtNucleiTof) ||
546     (hasExtTrk && !hasExtTof) || (!hasExtTrk && hasExtTof)
547     ){
548     if ( verbose ) printf(" ERROR: Mismatch between tracker and tof tree branches concerning extended tracking algorithm(s)\n");
549     if ( debug ) printf("hasNucleiTrk %i hasExtNucleiTrk %i hasExtTrk %i \n",hasNucleiTrk,hasExtNucleiTrk,hasExtTrk);
550     if ( debug ) printf("hasNucleiTof %i hasExtNucleiTof %i hasExtTof %i \n",hasNucleiTof,hasExtNucleiTof,hasExtTof);
551     throw -901;
552     }
553    
554 pam-mep 1.69 }
555 mocchiut 1.1 //
556     // Let's start!
557     //
558     // As a first thing we must check what we have to do: if run = 0 we must process all events in the file has been passed
559     // if run != 0 we must process only that run but first we have to check if the tree MyDetector2 already exist in the file
560     // if it exists we are reprocessing data and we must delete that entries, if not we must create it.
561     //
562 mocchiut 1.5 if ( run == 0 ) reproc = true;
563 mocchiut 1.1 //
564     //
565     // Output file is "outputfile"
566     //
567     if ( !file->IsOpen() ){
568     //printf(" OrbitalInfo - ERROR: cannot open file for writing\n");
569     throw -901;
570     };
571     //
572     // Retrieve GL_RUN variables from the level2 file
573     //
574     OrbitalInfoversion = OrbitalInfoInfo(false); // we should decide how to handle versioning system
575     //
576     // create an interface to RunInfo called "runinfo"
577     //
578     runinfo = new ItoRunInfo(file);
579     //
580     // open "Run" tree in level2 file, if not existing return an error (sngl != 0)
581     //
582     sgnl = 0;
583     sgnl = runinfo->Update(run, "ORB", OrbitalInfoversion);
584     //sgnl = runinfo->Read(run);
585    
586     if ( sgnl ){
587     //printf("OrbitalInfo - ERROR: RunInfo exited with non-zero status\n");
588     code = sgnl;
589     goto closeandexit;
590     } else {
591     sgnl = 0;
592     };
593     //
594     // number of events in the file BEFORE the first event of our run
595     //
596     nobefrun = runinfo->GetFirstEntry();
597     //
598     // total number of events in the file
599     //
600     totfileentries = runinfo->GetFileEntries();
601     //
602     // first file entry AFTER the last event of our run
603     //
604     noaftrun = runinfo->GetLastEntry() + 1;
605     //
606     // number of run to be processed
607     //
608     numbofrun = runinfo->GetNoRun();
609 mocchiut 1.40 totnorun = runinfo->GetRunEntries();
610 mocchiut 1.1 //
611     // Try to access the OrbitalInfo tree in the file, if it exists we are reprocessing data if not we are processing a new run
612     //
613     OrbitalInfotrclone = (TTree*)file->Get("OrbitalInfo");
614     //
615     if ( !OrbitalInfotrclone ){
616     //
617     // tree does not exist, we are not reprocessing
618     //
619     reproc = false;
620 mocchiut 1.5 if ( run == 0 ){
621 mocchiut 1.1 if (verbose) printf(" OrbitalInfo - WARNING: you are reprocessing data but OrbitalInfo tree does not exist!\n");
622     }
623 mocchiut 1.5 if ( runinfo->IsReprocessing() && run != 0 ) {
624 mocchiut 1.1 if (verbose) printf(" OrbitalInfo - WARNING: it seems you are not reprocessing data but OrbitalInfo\n versioning information already exists in RunInfo.\n");
625     }
626     } else {
627     //
628     // tree exists, we are reprocessing data. Are we reprocessing a single run or all the file?
629     //
630 mocchiut 1.6 OrbitalInfotrclone->SetAutoSave(900000000000000LL);
631 mocchiut 1.1 reproc = true;
632     //
633     //
634     if (verbose) printf("\n Preparing the pre-processing...\n");
635     //
636 mocchiut 1.22 if ( run == 0 || totnorun == 1 ){
637 mocchiut 1.1 //
638     // we are reprocessing all the file
639     // 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
640     //
641     reprocall = true;
642     //
643 mocchiut 1.57 if (verbose) printf("\n OrbitalInfo - WARNING: Reprocessing all runs\n Deleting old tree...\n");
644 mocchiut 1.1 //
645     } else {
646     //
647     // we are reprocessing a single run, we must copy to the new tree the events in the file which preceed the first event of the run
648     //
649     reprocall = false;
650     //
651 mocchiut 1.5 if (verbose) printf("\n OrbitalInfo - WARNING: Reprocessing run number %u \n",run);
652 mocchiut 1.1 //
653     // copying old tree to a new file
654     //
655 mocchiut 1.23 gSystem->MakeDirectory(OrbitalInfofolder.str().c_str());
656     myfold = true;
657 mocchiut 1.1 tempfile = new TFile(tempname.str().c_str(),"RECREATE");
658     tempOrbitalInfo = OrbitalInfotrclone->CloneTree(-1,"fast");
659     tempOrbitalInfo->SetName("OrbitalInfo-old");
660     tempfile->Write();
661 mocchiut 1.57 tempOrbitalInfo->Delete();
662 mocchiut 1.1 tempfile->Close();
663     }
664     //
665     // Delete the old tree from old file and memory
666     //
667 mocchiut 1.57 OrbitalInfotrclone->Clear();
668 mocchiut 1.1 OrbitalInfotrclone->Delete("all");
669     //
670     if (verbose) printf(" ...done!\n");
671     //
672     };
673     //
674     // create mydetector tree mydect
675     //
676     file->cd();
677     OrbitalInfotr = new TTree("OrbitalInfo-new","PAMELA OrbitalInfo data");
678 mocchiut 1.6 OrbitalInfotr->SetAutoSave(900000000000000LL);
679 mocchiut 1.30 orbitalinfo->Set();//ELENA **TEMPORANEO?**
680 mocchiut 1.1 OrbitalInfotr->Branch("OrbitalInfo","OrbitalInfo",&orbitalinfo);
681     //
682 mocchiut 1.73 // create new branches for new tracking algorithms
683     //
684     if ( hasNucleiTrk ){
685     torbNucleiTrk = new TClonesArray("OrbitalInfoTrkVar",1);
686     OrbitalInfotr->Branch("TrackNuclei",&torbNucleiTrk);
687     }
688     if ( hasExtNucleiTrk ){
689     torbExtNucleiTrk = new TClonesArray("OrbitalInfoTrkVar",1);
690     OrbitalInfotr->Branch("RecoveredTrackNuclei",&torbExtNucleiTrk);
691     }
692     if ( hasExtTrk ){
693     torbExtTrk = new TClonesArray("OrbitalInfoTrkVar",1);
694     OrbitalInfotr->Branch("RecoveredTrack",&torbExtTrk);
695     }
696    
697     //
698 mocchiut 1.1 if ( reproc && !reprocall ){
699     //
700     // open new file and retrieve also tree informations
701     //
702     tempfile = new TFile(tempname.str().c_str(),"READ");
703     OrbitalInfotrclone = (TTree*)tempfile->Get("OrbitalInfo-old");
704 mocchiut 1.6 OrbitalInfotrclone->SetAutoSave(900000000000000LL);
705 mocchiut 1.1 OrbitalInfotrclone->SetBranchAddress("OrbitalInfo",&orbitalinfoclone);
706     //
707     if ( nobefrun > 0 ){
708     if (verbose){
709 mayorov 1.93 printf("\n Pre-processing: copying events from the old tree before the processed run\n");
710     printf(" Copying %u events in the file which are before the beginning of the run %u \n",nobefrun,run);
711     printf(" Start copying at event number 0, end copying at event number %u \n",nobefrun);
712 mocchiut 1.1 }
713     for (UInt_t j = 0; j < nobefrun; j++){
714 mayorov 1.93 //
715     if ( OrbitalInfotrclone->GetEntry(j) <= 0 ) throw -36;
716     //
717     // copy orbitalinfoclone to mydec
718     //
719     // orbitalinfo->Clear();
720     //
721     memcpy(&orbitalinfo,&orbitalinfoclone,sizeof(orbitalinfoclone));
722     //
723     // Fill entry in the new tree
724     //
725     OrbitalInfotr->Fill();
726     //
727 mocchiut 1.1 };
728     if (verbose) printf(" Finished successful copying!\n");
729 mocchiut 1.57 };
730 mocchiut 1.1 };
731     //
732 mocchiut 1.44 //
733 mocchiut 1.1 // Get the list of run to be processed, if only one run has to be processed the list will contain one entry only.
734     //
735     runlist = runinfo->GetRunList();
736 mocchiut 1.81 if ( debug ){
737     printf("BEFORE LOOP ON RUN\n");
738     gObjectTable->Print();
739     }
740 mocchiut 1.1 //
741     // Loop over the run to be processed
742     //
743 mocchiut 1.81 for (UInt_t irun=0; irun < numbofrun; irun++){ //===>
744 mocchiut 1.57
745     L_QQ_Q_l_lower = new Quaternions();
746     RYPang_lower = new InclinationInfo();
747     L_QQ_Q_l_upper = new Quaternions();
748     RYPang_upper = new InclinationInfo();
749    
750 mocchiut 1.1 //
751     // retrieve the first run ID to be processed using the RunInfo list
752     //
753 mocchiut 1.15
754 mocchiut 1.1 idRun = runlist->At(irun);
755     if (verbose){
756     printf("\n\n\n ####################################################################### \n");
757     printf(" PROCESSING RUN NUMBER %i \n",(int)idRun);
758     printf(" ####################################################################### \n\n\n");
759     }
760     //
761 mocchiut 1.5 runinfo->ID_ROOT_L0 = 0;
762 mocchiut 1.1 //
763     // store in the runinfo class the GL_RUN variables for our run
764     //
765     sgnl = 0;
766     sgnl = runinfo->GetRunInfo(idRun);
767     if ( sgnl ){
768 mocchiut 1.5 if ( debug ) printf("\n OrbitalInfo - ERROR: RunInfo exited with non-zero status\n");
769 mocchiut 1.1 code = sgnl;
770     goto closeandexit;
771     } else {
772     sgnl = 0;
773     };
774     //
775     // now you can access that variables using the RunInfo class this way runinfo->ID_REG_RUN
776     //
777 mocchiut 1.5 if ( runinfo->ID_ROOT_L0 == 0 ){
778     if ( debug ) printf("\n OrbitalInfo - ERROR: no run with ID_RUN = %u \n\n Exiting... \n\n",idRun);
779 mocchiut 1.1 code = -5;
780     goto closeandexit;
781     };
782     //
783 mocchiut 1.5 // prepare the timesync for the db
784     //
785     dbtime = new GL_TIMESYNC(runinfo->ID_ROOT_L0,"ID",dbc);
786 mocchiut 1.15
787 mocchiut 1.5 //
788 mocchiut 1.1 // Search in the DB the path and name of the LEVEL0 file to be processed.
789     //
790 mocchiut 1.5 glroot->Query_GL_ROOT(runinfo->ID_ROOT_L0,dbc);
791 mocchiut 1.1 //
792     ftmpname.str("");
793     ftmpname << glroot->PATH.Data() << "/";
794     ftmpname << glroot->NAME.Data();
795     fname = ftmpname.str().c_str();
796 mocchiut 1.15 ftmpname.str("");
797 mocchiut 1.1 //
798 mocchiut 1.40 // print nout informations
799 mocchiut 1.1 //
800 mocchiut 1.5 totevent = runinfo->NEVENTS;
801 mocchiut 1.30 evfrom = runinfo->EV_FROM;
802 mocchiut 1.15 //cout<<"totevents = "<<totevent<<"\n";
803 mocchiut 1.1 if (verbose){
804     printf("\n LEVEL0 data file: %s \n",fname.Data());
805 mocchiut 1.5 printf(" RUN HEADER absolute time is: %u \n",runinfo->RUNHEADER_TIME);
806     printf(" RUN TRAILER absolute time is: %u \n",runinfo->RUNTRAILER_TIME);
807 mocchiut 1.15 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);
808 mocchiut 1.1 }//
809 mocchiut 1.27 //
810 mocchiut 1.28 // if ( !totevent ) goto closeandexit;
811 mocchiut 1.1 // Open Level0 file
812 mocchiut 1.59 if ( l0File ) l0File->Close();
813 mocchiut 1.1 l0File = new TFile(fname.Data());
814     if ( !l0File ) {
815 mocchiut 1.5 if ( debug ) printf(" OrbitalInfo - ERROR: problems opening Level0 file\n");
816 mocchiut 1.1 code = -6;
817     goto closeandexit;
818     };
819     l0tr = (TTree*)l0File->Get("Physics");
820     if ( !l0tr ) {
821 mocchiut 1.5 if ( debug ) printf(" OrbitalInfo - ERROR: no Physics tree in Level0 file\n");
822 mocchiut 1.1 l0File->Close();
823     code = -7;
824     goto closeandexit;
825     };
826 mocchiut 1.2 // EM: open header branch as well
827     l0head = l0tr->GetBranch("Header");
828     if ( !l0head ) {
829 mocchiut 1.5 if ( debug ) printf(" OrbitalInfo - ERROR: no Header branch in Level0 tree\n");
830 mocchiut 1.2 l0File->Close();
831     code = -8;
832     goto closeandexit;
833     };
834     l0tr->SetBranchAddress("Header", &eh);
835     // end EM
836 mocchiut 1.5 nevents = l0head->GetEntries();
837 mocchiut 1.1 //
838 mocchiut 1.28 if ( nevents < 1 && totevent ) {
839 mocchiut 1.5 if ( debug ) printf(" OrbitalInfo - ERROR: Level0 file is empty\n\n");
840 mocchiut 1.1 l0File->Close();
841     code = -11;
842     goto closeandexit;
843     };
844 mocchiut 1.20 //
845 mocchiut 1.28 if ( runinfo->EV_TO > nevents-1 && totevent ) {
846 mocchiut 1.5 if ( debug ) printf(" OrbitalInfo - ERROR: too few entries in the registry tree\n");
847 mocchiut 1.1 l0File->Close();
848     code = -12;
849     goto closeandexit;
850     };
851 mocchiut 1.55
852 mocchiut 1.15 ULong_t TimeSync = (ULong_t)dbtime->GetTimesync();
853     ULong_t ObtSync = (ULong_t)(dbtime->GetObt0()/1000);
854     ULong_t DeltaOBT = TimeSync - ObtSync;
855    
856     if ( debug ) printf(" 2 TimeSync %lu ObtSync %lu DeltaOBT %lu\n",(ULong_t)(dbtime->GetTimesync()/1000),(ULong_t)dbtime->GetObt0(),TimeSync-ObtSync);
857 mocchiut 1.37 //
858     // Read MCMDs from up to 11 files, 5 before and 5 after the present one in order to have some kind of inclination information
859     //
860     ch = new TChain("Mcmd","Mcmd");
861     //
862     // look in the DB to find the closest files to this run
863     //
864     TSQLResult *pResult = 0;
865     TSQLRow *Row = 0;
866     stringstream myquery;
867     UInt_t l0fid[10];
868     Int_t i = 0;
869     memset(l0fid,0,10*sizeof(Int_t));
870     //
871     myquery.str("");
872     myquery << "select ID_ROOT_L0 from GL_RUN where RUNHEADER_TIME<=" << runinfo->RUNHEADER_TIME << " group by ID_ROOT_L0 order by RUNHEADER_TIME desc limit 5;";
873     //
874     pResult = dbc->Query(myquery.str().c_str());
875     //
876     i = 9;
877     if( pResult ){
878     //
879     Row = pResult->Next();
880     //
881     while ( Row ){
882 mayorov 1.93 //
883     // store infos and exit
884     //
885     l0fid[i] = (UInt_t)atoll(Row->GetField(0));
886     i--;
887 mocchiut 1.81 if (Row){ // memleak!
888     delete Row;
889     Row = 0;
890     }
891 mayorov 1.93 Row = pResult->Next();
892     //
893 mocchiut 1.81 }
894     if (Row) delete Row;
895 mocchiut 1.37 pResult->Delete();
896 mocchiut 1.81 }
897 mocchiut 1.37 //
898     myquery.str("");
899     myquery << "select ID_ROOT_L0 from GL_RUN where RUNHEADER_TIME>" << runinfo->RUNHEADER_TIME << " group by ID_ROOT_L0 order by RUNHEADER_TIME asc limit 5;";
900     //
901     pResult = dbc->Query(myquery.str().c_str());
902     //
903     i = 0;
904     if( pResult ){
905     //
906     Row = pResult->Next();
907     //
908     while ( Row ){
909 mayorov 1.93 //
910     // store infos and exit
911     //
912     l0fid[i] = (UInt_t)atoll(Row->GetField(0));
913     i++;
914 mocchiut 1.81 if (Row){ // memleak!
915     delete Row;
916     Row = 0;
917     }
918 mayorov 1.93 Row = pResult->Next();
919     //
920 mocchiut 1.81 }
921     if (Row) delete Row;
922 mocchiut 1.37 pResult->Delete();
923 mocchiut 1.81 }
924 mocchiut 1.37 //
925     i = 0;
926     UInt_t previd = 0;
927     while ( i < 10 ){
928     if ( l0fid[i] && previd != l0fid[i] ){
929 mayorov 1.93 previd = l0fid[i];
930     myquery.str("");
931     myquery << "select PATH,NAME from GL_ROOT where ID=" << l0fid[i] << " ;";
932     //
933     pResult = dbc->Query(myquery.str().c_str());
934     //
935     if( pResult ){
936     //
937     Row = pResult->Next();
938     //
939     if ( debug ) printf(" Using inclination informations from file: %s \n",(((TString)gSystem->ExpandPathName(Row->GetField(0)))+"/"+(TString)Row->GetField(1)).Data());
940     ch->Add(((TString)gSystem->ExpandPathName(Row->GetField(0)))+"/"+(TString)Row->GetField(1));
941     //
942 mocchiut 1.81 if (Row) delete Row;
943 mayorov 1.93 pResult->Delete();
944     }
945 mocchiut 1.81 }
946 mocchiut 1.37 i++;
947 mocchiut 1.81 }
948 mocchiut 1.37 //
949     ch->SetBranchAddress("Mcmd",&mcmdev);
950     neventsm = ch->GetEntries();
951     if ( debug ) printf(" entries %u \n", neventsm);
952 mocchiut 1.15 //
953     if (neventsm == 0){
954 mocchiut 1.18 if ( debug ) printf("InclinationInfo - WARNING: No quaternions in this File");
955     code = 900;
956 mocchiut 1.15 }
957     //
958 mocchiut 1.73 Double_t lowerqtime = 0;
959 mocchiut 1.15 //
960 mocchiut 1.44 // init quaternions information from mcmd-packets
961 mocchiut 1.15 //
962     Bool_t isf = true;
963 mocchiut 1.44
964     vector<Float_t> q0;
965     vector<Float_t> q1;
966     vector<Float_t> q2;
967     vector<Float_t> q3;
968     vector<Double_t> qtime;
969     vector<Float_t> qPitch;
970     vector<Float_t> qRoll;
971     vector<Float_t> qYaw;
972     vector<Int_t> qmode;
973    
974 mocchiut 1.81 q0.reserve(4096);
975     q1.reserve(4096);
976     q2.reserve(4096);
977     q3.reserve(4096);
978     qtime.reserve(4096);
979     qPitch.reserve(4096);
980     qRoll.reserve(4096);
981     qYaw.reserve(4096);
982     qmode.reserve(4096);
983     if ( debug ) printf(" q0 capa %i \n",(int)q0.capacity());
984 mocchiut 1.44 Int_t nt = 0;
985     UInt_t must = 0;
986    
987 mocchiut 1.81 Int_t currentYear = 0;
988     Int_t previousYear = 0;
989    
990 mocchiut 1.15 //
991 mocchiut 1.1 // run over all the events of the run
992     //
993     if (verbose) printf("\n Ready to start! \n\n Processed events: \n\n");
994 mocchiut 1.81 if ( debug ){
995     printf("BEFORE LOOP ON EVENTS\n");
996     gObjectTable->Print();
997     }
998 mocchiut 1.1 //
999 mocchiut 1.44 //
1000 mocchiut 1.5 for ( re = runinfo->EV_FROM; re < (runinfo->EV_FROM+runinfo->NEVENTS); re++){
1001 mocchiut 1.81 //for ( re = runinfo->EV_FROM; re < (runinfo->EV_FROM+10); re++){
1002    
1003 mocchiut 1.1 //
1004 mayorov 1.93 if ( procev%1000 == 0 && procev > 0 && verbose ) printf(" %iK \n",procev/1000);
1005 mocchiut 1.15 if ( debug ) printf(" %i \n",procev);
1006 mocchiut 1.1 //
1007 mocchiut 1.43 if ( l0head->GetEntry(re) <= 0 ) throw -36;
1008 mocchiut 1.1 //
1009     // absolute time of this event
1010     //
1011 mocchiut 1.5 ph = eh->GetPscuHeader();
1012 mocchiut 1.15 atime = dbtime->DBabsTime(ph->GetOrbitalTime());
1013 mocchiut 1.40 if ( debug ) printf(" %i absolute time \n",procev);
1014 mocchiut 1.1 //
1015     // paranoid check
1016     //
1017 mocchiut 1.34 if ( (atime > (runinfo->RUNTRAILER_TIME+1)) || (atime < (runinfo->RUNHEADER_TIME-1)) ) {
1018 mayorov 1.93 if (verbose) printf(" OrbitalInfo - WARNING: event at time outside the run time window, skipping it\n");
1019     jumped++;
1020     // debug = true;
1021     continue;
1022 mocchiut 1.7 }
1023 mocchiut 1.30
1024 mocchiut 1.81 // just for testing:
1025     // if (re >= 5+runinfo->EV_FROM) atime += anni5;
1026     // if (re >= 7+runinfo->EV_FROM) atime += anni5;
1027     // if (re >= 9+runinfo->EV_FROM) atime += anni5;
1028    
1029     //
1030     // open IGRF files and do it only once if we are processing a full level2 file
1031     //
1032     Float_t kkyear;
1033     UInt_t kyear, kmonth, kday, khour, kmin, ksec;
1034     //
1035     TTimeStamp kt = TTimeStamp(atime, kTRUE);
1036     kt.GetDate(kTRUE, 0, &kyear, &kmonth, &kday);
1037     kt.GetTime(kTRUE, 0, &khour, &kmin, &ksec);
1038     kkyear = (float) kyear
1039     + (kmonth*31.+ (float) kday)/365.
1040     + (khour*3600.+kmin*60.+(float)ksec)/(24.*3600.*365.);
1041     currentYear = int(kkyear/5.) * 5;
1042     if ( debug ) printf(" prevy %i curry %i igrfloaded %i \n",previousYear,currentYear,igrfloaded);
1043     if ( currentYear != previousYear ) igrfloaded = false;
1044     previousYear = currentYear;
1045     if ( debug ) printf(" prevy %i curry %i igrfloaded %i \n",previousYear,currentYear,igrfloaded);
1046     //
1047     if ( !igrfloaded ){
1048    
1049     igrfloaded = true;
1050    
1051     parerror=glparam->Query_GL_PARAM(atime,302,dbc); // parameters stored in DB in GL_PRAM table
1052     if ( parerror<0 ) {
1053     code = parerror;
1054     goto closeandexit;
1055     }
1056     ltp1 = (Int_t)(glparam->PATH+glparam->NAME).Length();
1057     if ( verbose ) printf(" Reading Earth's Magnetic Field parameter file: %s \n",(glparam->PATH+glparam->NAME).Data());
1058     //
1059     if ( glparam->NAME.EndsWith("s.txt") || glparam->NAME.EndsWith("s.dat") ){
1060     if ( verbose ) printf("ERROR: Current date is beyond the latest secular variation file time span. Please update IGRF files to process data\n");
1061     throw -906;
1062     }
1063     //
1064     int isSecular = false;
1065     //
1066     parerror=glparam2->Query_GL_PARAM(atime+anni5,302,dbc); // parameters stored in DB in GL_PRAM table
1067     if ( parerror<0 ) {
1068     code = parerror;
1069     goto closeandexit;
1070     }
1071     ltp2 = (Int_t)(glparam2->PATH+glparam2->NAME).Length();
1072     if ( verbose ) printf(" Reading Earth's Magnetic Field parameter file: %s \n",(glparam2->PATH+glparam2->NAME).Data());
1073     if ( glparam2->NAME.EndsWith("s.txt") || glparam2->NAME.EndsWith("s.dat") ){
1074     isSecular = true;
1075     if ( verbose ) printf(" Using secular variation file and hence fortran subroutine 'extrapolation'\n");
1076     } else {
1077     if ( verbose ) printf(" Using two field measurement files and hence fortran subroutine 'interpolation'\n");
1078     }
1079     //
1080     initize_(&isSecular,(char *)(glparam->PATH+glparam->NAME).Data(),&ltp1,(char *)(glparam2->PATH+glparam2->NAME).Data(),&ltp2);
1081     //
1082     if (debug) cout<<"initize: "<<(char *)(glparam->PATH+glparam->NAME).Data()<<"\t"<<(char *)(glparam2->PATH+glparam2->NAME).Data()<<"\t isSecular? "<<isSecular<<endl;
1083    
1084     // GM_ScanIGRF(dbc, &G0, &G1, &H1);
1085     TString igrfFile1 = glparam->PATH+glparam->NAME;
1086     TString igrfFile2 = glparam2->PATH+glparam2->NAME;
1087     GM_SetIGRF(isSecular,igrfFile1,igrfFile2, &G0, &G1, &H1);
1088     }
1089     //
1090     // End IGRF stuff//
1091     //
1092    
1093 mocchiut 1.30 //
1094     // retrieve tof informations
1095     //
1096     if ( !reprocall ){
1097 mayorov 1.93 itr = nobefrun + (re - evfrom - jumped);
1098     //itr = re-(46438+200241);
1099 mocchiut 1.30 } else {
1100 mayorov 1.93 itr = runinfo->GetFirstEntry() + (re - evfrom - jumped);
1101 mocchiut 1.30 };
1102     //
1103     if ( !standalone ){
1104 mayorov 1.93 if ( itr > nevtofl2 ){
1105     if ( verbose ) printf(" OrbitalInfo - ERROR: no tof events with entry = %i in Level2 file\n",itr);
1106     if ( debug ) printf(" nobefrun %u re %u evfrom %u jumped %u reprocall %i \n",nobefrun,re,evfrom,jumped,reprocall);
1107     l0File->Close();
1108     code = -904;
1109     goto closeandexit;
1110     };
1111     //
1112     tof->Clear();
1113     //
1114 mocchiut 1.73 // Clones array must be cleared before going on
1115     //
1116     if ( hasNucleiTof ){
1117     tcNucleiTof->Delete();
1118     }
1119     if ( hasExtNucleiTof ){
1120     tcExtNucleiTof->Delete();
1121     }
1122     if ( hasExtTof ){
1123     tcExtTof->Delete();
1124     }
1125 mayorov 1.93 //
1126 mocchiut 1.73 if ( verbose ) printf(" get tof tree entries... entry = %i in Level2 file\n",itr);
1127 mayorov 1.93 if ( ttof->GetEntry(itr) <= 0 ){
1128 mocchiut 1.73 if ( verbose ) printf(" problems with tof tree entries... entry = %i in Level2 file\n",itr);
1129     if ( verbose ) printf(" nobefrun %u re %u evfrom %u jumped %u reprocall %i \n",nobefrun,re,evfrom,jumped,reprocall);
1130     throw -36;
1131 mayorov 1.93 }
1132 mocchiut 1.73 if ( verbose ) printf(" gat0\n");
1133 mayorov 1.93 //
1134 pam-mep 1.69 }
1135     //
1136     // retrieve tracker informations
1137     //
1138     if ( !standalone ){
1139 mayorov 1.93 if ( itr > nevtrkl2 ){
1140     if ( verbose ) printf(" OrbitalInfo - ERROR: no trk events with entry = %i in Level2 file\n",itr);
1141     if ( debug ) printf(" nobefrun %u re %u evfrom %u jumped %u reprocall %i \n",nobefrun,re,evfrom,jumped,reprocall);
1142     l0File->Close();
1143     code = -905;
1144     goto closeandexit;
1145     }
1146     //
1147 mocchiut 1.73 if ( verbose ) printf(" gat1\n");
1148 mayorov 1.93 trke->Clear();
1149     //
1150 mocchiut 1.73 // Clones array must be cleared before going on
1151     //
1152     if ( hasNucleiTrk ){
1153     if ( verbose ) printf(" gat2\n");
1154     tcNucleiTrk->Delete();
1155     if ( verbose ) printf(" gat3\n");
1156     torbNucleiTrk->Delete();
1157     }
1158     if ( hasExtNucleiTrk ){
1159     if ( verbose ) printf(" gat4\n");
1160     tcExtNucleiTrk->Delete();
1161     if ( verbose ) printf(" gat5\n");
1162     torbExtNucleiTrk->Delete();
1163     }
1164     if ( hasExtTrk ){
1165     if ( verbose ) printf(" gat6\n");
1166     tcExtTrk->Delete();
1167     if ( verbose ) printf(" gat7\n");
1168     torbExtTrk->Delete();
1169     }
1170 mayorov 1.93 //
1171 mocchiut 1.73 if ( verbose ) printf(" get trk tree entries... entry = %i in Level2 file\n",itr);
1172 mayorov 1.93 if ( ttrke->GetEntry(itr) <= 0 ) throw -36;
1173     //
1174 pam-mep 1.69 }
1175    
1176 mocchiut 1.1 //
1177     procev++;
1178     //
1179     // start processing
1180     //
1181 mocchiut 1.40 if ( debug ) printf(" %i start processing \n",procev);
1182 mocchiut 1.3 orbitalinfo->Clear();
1183 pamela 1.92
1184 mocchiut 1.5 //
1185 mocchiut 1.30 OrbitalInfoTrkVar *t_orb = new OrbitalInfoTrkVar();
1186     if( !(orbitalinfo->OrbitalInfoTrk) ) orbitalinfo->OrbitalInfoTrk = new TClonesArray("OrbitalInfoTrkVar",2);
1187     TClonesArray &tor = *orbitalinfo->OrbitalInfoTrk;
1188 pam-mep 1.69
1189     // Geomagnetic coordinates calculation variables
1190     GMtype_CoordSpherical CoordSpherical, DipoleSpherical;
1191     GMtype_CoordCartesian CoordCartesian, DipoleCartesian;
1192     GMtype_Model Model;
1193     GMtype_Pole Pole;
1194    
1195 mocchiut 1.30 //
1196 mocchiut 1.15 // Fill OBT, pkt_num and absTime
1197     //
1198 mocchiut 1.2 orbitalinfo->pkt_num = ph->GetCounter();
1199     orbitalinfo->OBT = ph->GetOrbitalTime();
1200 mocchiut 1.15 orbitalinfo->absTime = atime;
1201 mocchiut 1.40 if ( debug ) printf(" %i pktnum obt abstime \n",procev);
1202 mocchiut 1.15 //
1203     // Propagate the orbit from the tle time to atime, using SGP(D)4.
1204     //
1205 mocchiut 1.40 if ( debug ) printf(" %i sgp4 \n",procev);
1206 mocchiut 1.15 cCoordGeo coo;
1207 pamela 1.92 Float_t jyear=0.;
1208 mocchiut 1.7 //
1209 mocchiut 1.84 if(atime >= gltle->GetToTime() || atime < gltle->GetFromTime() ) { // AGH! bug when reprocessing??
1210 pamela 1.92
1211 mayorov 1.93 if ( !gltle->Query(atime, dbc) ){
1212     //
1213     // Compute the magnetic dipole moment.
1214     //
1215     if ( debug ) printf(" %i compute magnetic dipole moment \n",procev);
1216     UInt_t year, month, day, hour, min, sec;
1217     //
1218     TTimeStamp t = TTimeStamp(atime, kTRUE);
1219     t.GetDate(kTRUE, 0, &year, &month, &day);
1220     t.GetTime(kTRUE, 0, &hour, &min, &sec);
1221     jyear = (float) year
1222     + (month*31.+ (float) day)/365.
1223     + (hour*3600.+min*60.+(float)sec)/(24.*3600.*365.);
1224     //
1225     if ( debug ) printf(" %i compute magnetic dipole moment get dipole moment for year\n",procev);
1226     if ( debug ) printf(" %i jyear %f dimo %f \n",procev,jyear,dimo);
1227     feldcof_(&jyear, &dimo); // get dipole moment for year
1228     if ( debug ) printf(" %i compute magnetic dipole moment end\n",procev);
1229    
1230     // GM_TimeAdjustCoefs(year, jyear, G0, G1, H1, &Model);
1231     GM_TimeAdjustCoefs(GM_STARTYEAR, (jyear-currentYear+GM_STARTYEAR), G0, G1, H1, &Model); // EM: input this way due to the new way of storing data into Gn,H1 and to avoid changing GM_Time...
1232     GM_PoleLocation(Model, &Pole);
1233     } else {
1234     code = -56;
1235     goto closeandexit;
1236     };
1237 mocchiut 1.7 }
1238 mocchiut 1.15 coo = getCoo(atime, gltle->GetFromTime(), gltle->GetTle());
1239     //
1240     cOrbit orbits(*gltle->GetTle());
1241     //
1242     // synchronize with quaternions data
1243     //
1244 mocchiut 1.18 if ( isf && neventsm>0 ){
1245 mayorov 1.93 //
1246     // First event
1247     //
1248     isf = false;
1249     // upperqtime = atime;
1250     lowerqtime = runinfo->RUNHEADER_TIME;
1251     for ( ik = 0; ik < neventsm; ik++){ //number of macrocommad packets
1252     if ( ch->GetEntry(ik) <= 0 ) throw -36;
1253     tmpSize = mcmdev->Records->GetEntries();
1254     // numrec = tmpSize;
1255 mocchiut 1.70 if ( debug ) cout << "packet number " << ik <<"\tnumber of subpackets is " << tmpSize << endl;
1256 mayorov 1.93 for (Int_t j3 = 0;j3<tmpSize;j3++){ //number of subpackets
1257     mcmdrc = (pamela::McmdRecord*)mcmdev->Records->At(j3);
1258     if ( mcmdrc ){ // missing inclination bug [8RED 090116]
1259     if ( debug ) printf(" pluto \n");
1260     if ((int)mcmdrc->ID1 == 226 && mcmdrc->Mcmd_Block_crc_ok == 1){ //Check that it is Inclination Packet
1261     L_QQ_Q_l_upper->fill(mcmdrc->McmdData);
1262     for (UInt_t ui = 0; ui < 6; ui++){
1263     if (ui>0){
1264     if (L_QQ_Q_l_upper->time[ui]>L_QQ_Q_l_upper->time[0]){
1265     Double_t u_time = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[ui]*1000-DeltaOBT*1000));
1266     Int_t recSize = recqtime.size();
1267     if(lowerqtime > recqtime[recSize-1]){
1268     // to avoid interpolation between bad quaternions arrays
1269     if(sqrt(pow(L_QQ_Q_l_upper->quat[ui][0],2)+pow(L_QQ_Q_l_upper->quat[ui][1],2)+pow(L_QQ_Q_l_upper->quat[ui][2],2)+pow(L_QQ_Q_l_upper->quat[ui][3],2))>0.99999){
1270     Int_t sizeqmcmd = qtime.size();
1271     inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);
1272     qtime[sizeqmcmd]=u_time;
1273     q0[sizeqmcmd]=L_QQ_Q_l_upper->quat[ui][0];
1274     q1[sizeqmcmd]=L_QQ_Q_l_upper->quat[ui][1];
1275     q2[sizeqmcmd]=L_QQ_Q_l_upper->quat[ui][2];
1276     q3[sizeqmcmd]=L_QQ_Q_l_upper->quat[ui][3];
1277     qmode[sizeqmcmd]=holeq(lowerqtime,qtime[sizeqmcmd],L_QQ_Q_l_lower,L_QQ_Q_l_upper,ui);
1278     lowerqtime = u_time;
1279     orbits.getPosition((double) (u_time - gltle->GetFromTime())/60., &eCi);
1280     RYPang_upper->TransAngle(eCi.getPos().m_x,eCi.getPos().m_y,eCi.getPos().m_z,eCi.getVel().m_x,eCi.getVel().m_y,eCi.getVel().m_z,L_QQ_Q_l_upper->quat[ui][0],L_QQ_Q_l_upper->quat[ui][1],L_QQ_Q_l_upper->quat[ui][2],L_QQ_Q_l_upper->quat[ui][3]);
1281     qRoll[sizeqmcmd]=RYPang_upper->Kren;
1282     qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;
1283     qPitch[sizeqmcmd]=RYPang_upper->Tangazh;
1284     }
1285     }
1286     for(Int_t mu = nt;mu<recSize;mu++){
1287     if(recqtime[mu]>lowerqtime && recqtime[mu]<u_time){
1288     if(sqrt(pow(recq0[mu],2)+pow(recq1[mu],2)+pow(recq2[mu],2)+pow(recq3[mu],2))>0.99999){
1289     nt=mu;
1290     Int_t sizeqmcmd = qtime.size();
1291     inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);
1292     qtime[sizeqmcmd]=recqtime[mu];
1293     q0[sizeqmcmd]=recq0[mu];
1294     q1[sizeqmcmd]=recq1[mu];
1295     q2[sizeqmcmd]=recq2[mu];
1296     q3[sizeqmcmd]=recq3[mu];
1297     qmode[sizeqmcmd]=-10;
1298     orbits.getPosition((double) (qtime[sizeqmcmd] - gltle->GetFromTime())/60., &eCi);
1299     RYPang_upper->TransAngle(eCi.getPos().m_x,eCi.getPos().m_y,eCi.getPos().m_z,eCi.getVel().m_x,eCi.getVel().m_y,eCi.getVel().m_z,recq0[mu],recq1[mu],recq2[mu],recq3[mu]);
1300     qRoll[sizeqmcmd]=RYPang_upper->Kren;
1301     qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;
1302     qPitch[sizeqmcmd]=RYPang_upper->Tangazh;
1303     }
1304     }
1305     if(recqtime[mu]>=u_time){
1306     if(sqrt(pow(L_QQ_Q_l_upper->quat[ui][0],2)+pow(L_QQ_Q_l_upper->quat[ui][1],2)+pow(L_QQ_Q_l_upper->quat[ui][2],2)+pow(L_QQ_Q_l_upper->quat[ui][3],2))>0.99999){
1307     Int_t sizeqmcmd = qtime.size();
1308     inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);
1309     qtime[sizeqmcmd]=u_time;
1310     q0[sizeqmcmd]=L_QQ_Q_l_upper->quat[ui][0];
1311     q1[sizeqmcmd]=L_QQ_Q_l_upper->quat[ui][1];
1312     q2[sizeqmcmd]=L_QQ_Q_l_upper->quat[ui][2];
1313     q3[sizeqmcmd]=L_QQ_Q_l_upper->quat[ui][3];
1314     qmode[sizeqmcmd]=holeq(lowerqtime,qtime[sizeqmcmd],L_QQ_Q_l_lower,L_QQ_Q_l_upper,ui);
1315     lowerqtime = u_time;
1316     orbits.getPosition((double) (u_time - gltle->GetFromTime())/60., &eCi);
1317     RYPang_upper->TransAngle(eCi.getPos().m_x,eCi.getPos().m_y,eCi.getPos().m_z,eCi.getVel().m_x,eCi.getVel().m_y,eCi.getVel().m_z,L_QQ_Q_l_upper->quat[ui][0],L_QQ_Q_l_upper->quat[ui][1],L_QQ_Q_l_upper->quat[ui][2],L_QQ_Q_l_upper->quat[ui][3]);
1318     qRoll[sizeqmcmd]=RYPang_upper->Kren;
1319     qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;
1320     qPitch[sizeqmcmd]=RYPang_upper->Tangazh;
1321     break;
1322     }
1323     }
1324     }
1325     }
1326     }else{
1327     //if ( debug ) printf(" here2 %i \n",ui);
1328     Double_t u_time = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
1329     if(lowerqtime>u_time)nt=0;
1330     Int_t recSize = recqtime.size();
1331     if(lowerqtime > recqtime[recSize-1]){
1332     if(sqrt(pow(L_QQ_Q_l_upper->quat[ui][0],2)+pow(L_QQ_Q_l_upper->quat[ui][1],2)+pow(L_QQ_Q_l_upper->quat[ui][2],2)+pow(L_QQ_Q_l_upper->quat[ui][3],2))>0.99999){
1333     Int_t sizeqmcmd = qtime.size();
1334     inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);
1335     qtime[sizeqmcmd]=u_time;
1336     q0[sizeqmcmd]=L_QQ_Q_l_upper->quat[0][0];
1337     q1[sizeqmcmd]=L_QQ_Q_l_upper->quat[0][1];
1338     q2[sizeqmcmd]=L_QQ_Q_l_upper->quat[0][2];
1339     q3[sizeqmcmd]=L_QQ_Q_l_upper->quat[0][3];
1340     qmode[sizeqmcmd]=holeq(lowerqtime,qtime[sizeqmcmd],L_QQ_Q_l_lower,L_QQ_Q_l_upper,ui);
1341     lowerqtime = u_time;
1342     orbits.getPosition((double) (u_time - gltle->GetFromTime())/60., &eCi);
1343     RYPang_upper->TransAngle(eCi.getPos().m_x,eCi.getPos().m_y,eCi.getPos().m_z,eCi.getVel().m_x,eCi.getVel().m_y,eCi.getVel().m_z,L_QQ_Q_l_upper->quat[0][0],L_QQ_Q_l_upper->quat[0][1],L_QQ_Q_l_upper->quat[0][2],L_QQ_Q_l_upper->quat[0][3]);
1344     qRoll[sizeqmcmd]=RYPang_upper->Kren;
1345     qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;
1346     qPitch[sizeqmcmd]=RYPang_upper->Tangazh;
1347     }
1348     }
1349     for(Int_t mu = nt;mu<recSize;mu++){
1350     if(recqtime[mu]>lowerqtime && recqtime[mu]<u_time){
1351     if(sqrt(pow(recq0[mu],2)+pow(recq1[mu],2)+pow(recq2[mu],2)+pow(recq3[mu],2))>0.99999){
1352     Int_t sizeqmcmd = qtime.size();
1353     inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);
1354     qtime[sizeqmcmd]=recqtime[mu];
1355     q0[sizeqmcmd]=recq0[mu];
1356     q1[sizeqmcmd]=recq1[mu];
1357     q2[sizeqmcmd]=recq2[mu];
1358     q3[sizeqmcmd]=recq3[mu];
1359     qmode[sizeqmcmd]=-10;
1360     orbits.getPosition((double) (qtime[sizeqmcmd] - gltle->GetFromTime())/60., &eCi);
1361     RYPang_upper->TransAngle(eCi.getPos().m_x,eCi.getPos().m_y,eCi.getPos().m_z,eCi.getVel().m_x,eCi.getVel().m_y,eCi.getVel().m_z,recq0[mu],recq1[mu],recq2[mu],recq3[mu]);
1362     qRoll[sizeqmcmd]=RYPang_upper->Kren;
1363     qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;
1364     qPitch[sizeqmcmd]=RYPang_upper->Tangazh;
1365     }
1366     }
1367     if(recqtime[mu]>=u_time){
1368     if(sqrt(pow(L_QQ_Q_l_upper->quat[0][0],2)+pow(L_QQ_Q_l_upper->quat[0][1],2)+pow(L_QQ_Q_l_upper->quat[0][2],2)+pow(L_QQ_Q_l_upper->quat[0][3],2))>0.99999){
1369     Int_t sizeqmcmd = qtime.size();
1370     inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);
1371     qtime[sizeqmcmd]=u_time;
1372     q0[sizeqmcmd]=L_QQ_Q_l_upper->quat[0][0];
1373     q1[sizeqmcmd]=L_QQ_Q_l_upper->quat[0][1];
1374     q2[sizeqmcmd]=L_QQ_Q_l_upper->quat[0][2];
1375     q3[sizeqmcmd]=L_QQ_Q_l_upper->quat[0][3];
1376     qmode[sizeqmcmd]=holeq(lowerqtime,qtime[sizeqmcmd],L_QQ_Q_l_lower,L_QQ_Q_l_upper,ui);
1377     lowerqtime = u_time;
1378     orbits.getPosition((double) (u_time - gltle->GetFromTime())/60., &eCi);
1379     RYPang_upper->TransAngle(eCi.getPos().m_x,eCi.getPos().m_y,eCi.getPos().m_z,eCi.getVel().m_x,eCi.getVel().m_y,eCi.getVel().m_z,L_QQ_Q_l_upper->quat[0][0],L_QQ_Q_l_upper->quat[0][1],L_QQ_Q_l_upper->quat[0][2],L_QQ_Q_l_upper->quat[0][3]);
1380     qRoll[sizeqmcmd]=RYPang_upper->Kren;
1381     qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;
1382     qPitch[sizeqmcmd]=RYPang_upper->Tangazh;
1383     CopyQ(L_QQ_Q_l_lower,L_QQ_Q_l_upper);
1384     break;
1385     }
1386     }
1387     }
1388     }
1389     }
1390     }
1391     }
1392     }
1393     }
1394    
1395     if(qtime.size()==0){ // in case if no orientation information in data
1396     if ( debug ) cout << "qtime.size() = 0" << endl;
1397     for(UInt_t my=0;my<recqtime.size();my++){
1398     if(sqrt(pow(recq0[my],2)+pow(recq1[my],2)+pow(recq2[my],2)+pow(recq3[my],2))>0.99999){
1399     Int_t sizeqmcmd = qtime.size();
1400     inclresize(qtime,q0,q1,q2,q3,qmode,qRoll,qPitch,qYaw);
1401     qtime[sizeqmcmd]=recqtime[my];
1402     q0[sizeqmcmd]=recq0[my];
1403     q1[sizeqmcmd]=recq1[my];
1404     q2[sizeqmcmd]=recq2[my];
1405     q3[sizeqmcmd]=recq3[my];
1406     qmode[sizeqmcmd]=-10;
1407     orbits.getPosition((double) (qtime[sizeqmcmd] - gltle->GetFromTime())/60., &eCi);
1408     RYPang_upper->TransAngle(eCi.getPos().m_x,eCi.getPos().m_y,eCi.getPos().m_z,eCi.getVel().m_x,eCi.getVel().m_y,eCi.getVel().m_z,recq0[my],recq1[my],recq2[my],recq3[my]);
1409     qRoll[sizeqmcmd]=RYPang_upper->Kren;
1410     qYaw[sizeqmcmd]=RYPang_upper->Ryskanie;
1411     qPitch[sizeqmcmd]=RYPang_upper->Tangazh;
1412     }
1413     }
1414     }
1415     //if ( debug ) printf(" puffi \n");
1416     Double_t tmin = 9999999999.;
1417     Double_t tmax = 0.;
1418     for(UInt_t tre = 0;tre<qtime.size();tre++){
1419     if(qtime[tre]>tmax)tmax = qtime[tre];
1420     if(qtime[tre]<tmin)tmin = qtime[tre];
1421     }
1422     // sorting quaternions by time
1423     Bool_t t = true;
1424     while(t){
1425     t=false;
1426     for(UInt_t i=0;i<qtime.size()-1;i++){
1427     if(qtime[i]>qtime[i+1]){
1428     Double_t tmpr = qtime[i];
1429     qtime[i]=qtime[i+1];
1430     qtime[i+1] = tmpr;
1431     tmpr = q0[i];
1432     q0[i]=q0[i+1];
1433     q0[i+1] = tmpr;
1434     tmpr = q1[i];
1435     q1[i]=q1[i+1];
1436     q1[i+1] = tmpr;
1437     tmpr = q2[i];
1438     q2[i]=q2[i+1];
1439     q2[i+1] = tmpr;
1440     tmpr = q3[i];
1441     q3[i]=q3[i+1];
1442     q3[i+1] = tmpr;
1443     tmpr = qRoll[i];
1444     qRoll[i]=qRoll[i+1];
1445     qRoll[i+1] = tmpr;
1446     tmpr = qYaw[i];
1447     qYaw[i]=qYaw[i+1];
1448     qYaw[i+1] = tmpr;
1449     tmpr = qPitch[i];
1450     qPitch[i]=qPitch[i+1];
1451     qPitch[i+1] = tmpr;
1452     t=true;
1453     }
1454 mocchiut 1.70 }
1455 mocchiut 1.65 }
1456 mayorov 1.93 if ( debug ){
1457     cout << "we have loaded quaternions: size of quaternions set is "<< qtime.size() << endl;
1458     for(UInt_t i=0;i<qtime.size();i++) cout << qtime[i] << "\t";
1459     cout << endl << endl;
1460     Int_t lopu;
1461     cin >> lopu;
1462     }
1463 mocchiut 1.44 } // if we processed first event
1464    
1465     //Filling Inclination information
1466     Double_t incli = 0;
1467 mocchiut 1.46 if ( qtime.size() > 1 ){
1468 pamela 1.92 if(atime<qtime[0]){
1469     for(UInt_t mu = 1;mu<qtime.size()-1;mu++){
1470     if(qtime[mu]>qtime[0]){
1471     incli = (qPitch[mu]-qPitch[0])/(qtime[mu]-qtime[0]);
1472     orbitalinfo->theta = incli*atime+qPitch[mu]-incli*qtime[mu];
1473     incli = (qRoll[mu]-qRoll[0])/(qtime[mu]-qtime[0]);
1474     orbitalinfo->etha = incli*atime+qRoll[mu]-incli*qtime[mu];
1475     incli = (qYaw[mu]-qYaw[0])/(qtime[mu]-qtime[0]);
1476     orbitalinfo->phi = incli*atime+qYaw[mu]-incli*qtime[mu];
1477    
1478     incli = (q0[mu]-q0[0])/(qtime[mu]-qtime[0]);
1479     orbitalinfo->q0 = incli*atime+q0[mu]-incli*qtime[mu];
1480     incli = (q1[mu]-q1[0])/(qtime[mu]-qtime[0]);
1481     orbitalinfo->q1 = incli*atime+q1[mu]-incli*qtime[mu];
1482     incli = (q2[mu]-q2[0])/(qtime[mu]-qtime[0]);
1483     orbitalinfo->q2 = incli*atime+q2[mu]-incli*qtime[mu];
1484     incli = (q3[mu]-q3[0])/(qtime[mu]-qtime[0]);
1485     orbitalinfo->q3 = incli*atime+q3[mu]-incli*qtime[mu];
1486     orbitalinfo->TimeGap=qtime[0]-atime;
1487     break;
1488     }
1489     }
1490     }
1491     Float_t eend=qtime.size()-1;
1492     if(atime>qtime[eend]){
1493     for(UInt_t mu=eend-1;mu>=0;mu--){
1494     if(qtime[mu]<qtime[eend]){
1495     incli = (qPitch[eend]-qPitch[mu])/(qtime[eend]-qtime[mu]);
1496     orbitalinfo->theta = incli*atime+qPitch[eend]-incli*qtime[eend];
1497     incli = (qRoll[eend]-qRoll[mu])/(qtime[eend]-qtime[mu]);
1498     orbitalinfo->etha = incli*atime+qRoll[eend]-incli*qtime[eend];
1499     incli = (qYaw[eend]-qYaw[mu])/(qtime[eend]-qtime[mu]);
1500     orbitalinfo->phi = incli*atime+qYaw[eend]-incli*qtime[eend];
1501    
1502     incli = (q0[eend]-q0[mu])/(qtime[eend]-qtime[mu]);
1503     orbitalinfo->q0 = incli*atime+q0[eend]-incli*qtime[eend];
1504     incli = (q1[eend]-q1[mu])/(qtime[eend]-qtime[mu]);
1505     orbitalinfo->q1 = incli*atime+q1[eend]-incli*qtime[eend];
1506     incli = (q2[eend]-q2[mu])/(qtime[eend]-qtime[mu]);
1507     orbitalinfo->q2 = incli*atime+q2[eend]-incli*qtime[eend];
1508     incli = (q3[eend]-q3[mu])/(qtime[eend]-qtime[mu]);
1509     orbitalinfo->q3 = incli*atime+q3[eend]-incli*qtime[eend];
1510     break;
1511     }
1512     }
1513     }
1514 mocchiut 1.65 for(UInt_t mu = must;mu<qtime.size()-1;mu++){
1515     if ( debug ) printf(" ??grfuffi %i sixe %i must %i \n",mu,qtime.size()-1,must);
1516     if(qtime[mu+1]>qtime[mu]){
1517 mocchiut 1.70 if ( debug ) cout << "qtime[" << mu << "] = " << qtime[mu] << "\tqtime[" << mu+1 << "] = " << qtime[mu+1] << "\tatime = " << atime << endl;
1518 mocchiut 1.65 if(atime<=qtime[mu+1] && atime>=qtime[mu]){
1519 mocchiut 1.70 if ( debug ) cout << "here we have found proper quaternions for interpolation: mu = "<<mu<<endl;
1520 mocchiut 1.65 must = mu;
1521     incli = (qPitch[mu+1]-qPitch[mu])/(qtime[mu+1]-qtime[mu]);
1522     orbitalinfo->theta = incli*atime+qPitch[mu+1]-incli*qtime[mu+1];
1523     incli = (qRoll[mu+1]-qRoll[mu])/(qtime[mu+1]-qtime[mu]);
1524     orbitalinfo->etha = incli*atime+qRoll[mu+1]-incli*qtime[mu+1];
1525     incli = (qYaw[mu+1]-qYaw[mu])/(qtime[mu+1]-qtime[mu]);
1526     orbitalinfo->phi = incli*atime+qYaw[mu+1]-incli*qtime[mu+1];
1527    
1528     incli = (q0[mu+1]-q0[mu])/(qtime[mu+1]-qtime[mu]);
1529     orbitalinfo->q0 = incli*atime+q0[mu+1]-incli*qtime[mu+1];
1530     incli = (q1[mu+1]-q1[mu])/(qtime[mu+1]-qtime[mu]);
1531     orbitalinfo->q1 = incli*atime+q1[mu+1]-incli*qtime[mu+1];
1532     incli = (q2[mu+1]-q2[mu])/(qtime[mu+1]-qtime[mu]);
1533     orbitalinfo->q2 = incli*atime+q2[mu+1]-incli*qtime[mu+1];
1534     incli = (q3[mu+1]-q3[mu])/(qtime[mu+1]-qtime[mu]);
1535     orbitalinfo->q3 = incli*atime+q3[mu+1]-incli*qtime[mu+1];
1536 malakhov 1.74 Float_t tg = (qtime[mu+1]-qtime[mu])/1000.0;
1537 mocchiut 1.70 if(tg>=1) tg=0.00;
1538 pamela 1.92 orbitalinfo->TimeGap = TMath::Min(TMath::Abs(qtime[mu+1]-atime),TMath::Abs(atime-qtime[mu]))+tg;//qtime[mu+1]-qtime[mu];
1539 mocchiut 1.70 orbitalinfo->mode = qmode[mu+1];
1540     //if(atime==qtime[mu] || atime==qtime[mu+1]) orbitalinfo->qkind = 0; else orbitalinfo->qkind=1;
1541 mocchiut 1.65 //if(qmode[mu+1]==-10) orbitalinfo->R10r = true;else orbitalinfo->R10r = false;
1542     if ( debug ) printf(" grfuffi4 %i \n",mu);
1543     break;
1544     }
1545     }
1546     }
1547 mocchiut 1.46 }
1548 mocchiut 1.65 if ( debug ) printf(" grfuffi5 \n");
1549 mocchiut 1.30 //
1550     // ops no inclination information
1551     //
1552 mocchiut 1.65
1553 mocchiut 1.30 if ( orbitalinfo->q0< -999 || orbitalinfo->q1 < -999 || orbitalinfo->q2 < -999 || orbitalinfo->q3 < -999 || orbitalinfo->q0 != orbitalinfo->q0 || orbitalinfo->q1 != orbitalinfo->q1 || orbitalinfo->q2 != orbitalinfo->q2 || orbitalinfo->q3 != orbitalinfo->q3 ){
1554 mayorov 1.93 if (debug) cout << "Warning: no iclination information "<< endl;
1555     orbitalinfo->mode = 10;
1556     orbitalinfo->q0 = -1000.;
1557     orbitalinfo->q1 = -1000.;
1558     orbitalinfo->q2 = -1000.;
1559     orbitalinfo->q3 = -1000.;
1560     orbitalinfo->etha = -1000.;
1561     orbitalinfo->phi = -1000.;
1562     orbitalinfo->theta = -1000.;
1563     orbitalinfo->TimeGap = -1000.;
1564     TMatrixD Iij(3,3);
1565     Iij(0,0)=0; Iij(0,1)=0; Iij(0,2)=0;
1566     Iij(1,0)=0; Iij(1,1)=0; Iij(1,2)=0;
1567     Iij(2,0)=0; Iij(2,1)=0; Iij(2,2)=0;
1568     Iij.Zero();
1569     orbitalinfo->Iij.ResizeTo(Iij);
1570     orbitalinfo->Iij = Iij;
1571     //orbitalinfo->qkind = -1000;
1572    
1573     // if ( debug ){
1574     // Int_t lopu;
1575     // cin >> lopu;
1576     // }
1577     if ( debug ) printf(" grfuffi6 \n");
1578     }
1579     //
1580     if ( debug ) printf(" filling \n");
1581     // #########################################################################################################################
1582     //
1583     // fill orbital positions
1584     //
1585     // Build coordinates in the right range. We want to convert,
1586     // longitude from (0, 2*pi) to (-180deg, 180deg). Altitude is
1587     // in meters.
1588     lon = (coo.m_Lon > M_PI) ? rad2deg(coo.m_Lon - 2*M_PI) : rad2deg(coo.m_Lon);
1589     lat = rad2deg(coo.m_Lat);
1590     alt = coo.m_Alt;
1591    
1592     cOrbit orbits2(*gltle->GetTle());
1593     orbits2.getPosition((double) (atime - gltle->GetFromTime())/60., &eCi);
1594     // Float_t x=eCi.getPos().m_x;
1595     // Float_t y=eCi.getPos().m_y;
1596     // Float_t z=eCi.getPos().m_z;
1597    
1598     TVector3 V(eCi.getVel().m_x,eCi.getVel().m_y,eCi.getVel().m_z);
1599     TVector3 Pos(eCi.getPos().m_x,eCi.getPos().m_y,eCi.getPos().m_z);
1600    
1601     Float_t dlon=Pos.Phi()*TMath::RadToDeg()-lon;
1602 mocchiut 1.70
1603 mayorov 1.93 Pos.RotateZ(-dlon*TMath::DegToRad());
1604     V.RotateZ(-dlon*TMath::DegToRad());
1605     Float_t diro;
1606     if(V.Z()>0) diro=1; else diro=-1;
1607    
1608     // 10REDNEW
1609     Int_t errq=0;
1610     Int_t azim=0;
1611     Int_t qual=0;
1612     Int_t MU=0;
1613     for(UInt_t mu = 0;mu<RTtime2.size()-1;mu++){
1614     if(atime<RTstart[mu+1] && atime>=RTstart[mu]){
1615     errq=RTerrq[mu];
1616     azim=RTazim[mu];
1617     qual=RTqual[mu];
1618     MU=mu;
1619     break;
1620     }
1621 mocchiut 1.70 }
1622 mayorov 1.93 orbitalinfo->errq = errq;
1623     orbitalinfo->azim = azim;
1624     orbitalinfo->rtqual=qual;
1625     orbitalinfo->qkind = 0;
1626    
1627     if ( debug ) printf(" coord done \n");
1628     if( lon<180 && lon>-180 && lat<90 && lat>-90 && alt>0 ){
1629     orbitalinfo->lon = lon;
1630     orbitalinfo->lat = lat;
1631     orbitalinfo->alt = alt;
1632     orbitalinfo->V = V;
1633    
1634     // GMtype_CoordGeodetic location;
1635     location.lambda = lon;
1636     location.phi = lat;
1637     location.HeightAboveEllipsoid = alt;
1638    
1639     GM_GeodeticToSpherical(Ellip, location, &CoordSpherical);
1640     GM_SphericalToCartesian(CoordSpherical, &CoordCartesian);
1641     GM_EarthCartToDipoleCartCD(Pole, CoordCartesian, &DipoleCartesian);
1642     GM_CartesianToSpherical(DipoleCartesian, &DipoleSpherical);
1643     orbitalinfo->londip = DipoleSpherical.lambda;
1644     orbitalinfo->latdip = DipoleSpherical.phig;
1645    
1646     if(debug)cout<<"geodetic:\t"<<lon<<"\t"<<lat<<"\tgeomagnetic:\t"<<orbitalinfo->londip<<"\t"<<orbitalinfo->latdip<<endl;
1647    
1648     //
1649     // compute mag field components and L shell.
1650     //
1651     if ( debug ) printf(" call igrf feldg \n");
1652     feldg_(&lat, &lon, &alt, &bnorth, &beast, &bdown, &babs);
1653     if ( debug ) printf(" call igrf shellg \n");
1654     shellg_(&lat, &lon, &alt, &dimo, &xl, &icode, &bab1);
1655     if ( debug ) printf(" call igrf findb \n");
1656     findb0_(&stps, &bdel, &value, &bequ, &rr0);
1657     //
1658     if ( debug ) printf(" done igrf \n");
1659     orbitalinfo->Bnorth = bnorth;
1660     orbitalinfo->Beast = beast;
1661     orbitalinfo->Bdown = bdown;
1662     orbitalinfo->Babs = babs;
1663     orbitalinfo->M = dimo;
1664     orbitalinfo->BB0 = babs/bequ;
1665     orbitalinfo->L = xl;
1666     // Set Stormer vertical cutoff using L shell.
1667     orbitalinfo->cutoffsvl = 14.295 / (xl*xl); //
1668     if(debug)cout << "L = " << xl << "\tM = " << dimo << "\tvertical cutoff: "<< orbitalinfo->cutoffsvl << endl;
1669    
1670     // ---------- Forwarded message ----------
1671     // Date: Wed, 09 May 2012 12:16:47 +0200
1672     // From: Alessandro Bruno <alessandro.bruno@ba.infn.it>
1673     // To: Mirko Boezio <mirko.boezio@ts.infn.it>
1674     // Cc: Francesco S. Cafagna <Francesco.Cafagna@ba.infn.it>
1675     // Subject: Störmer vertical cutoff
1676    
1677     // Ciao Mirko,
1678     // volevo segnalarti che il valore dello Störmer vertical cutoff nel Level2 è
1679     // sovrastimato di circa il 4%.
1680     // Dopo un'approfondita analisi con l'IGRF-05 abbiamo ricavano un valore pari
1681     // a: 14.295 / L^2 anzichè 14.9 / L^2, valore obsoleto in quanto riferito agli
1682     // anni '50.
1683     //
1684     // 14.9/(xl*xl);
1685     orbitalinfo->igrf_icode = (Float_t)icode;
1686     //
1687     } //check lon lat alt
1688     //
1689     if ( debug ) printf(" pitch angle \n");
1690     //
1691     // pitch angles
1692     //
1693     if( orbitalinfo->TimeGap>=0){
1694     //
1695     if ( debug ) printf(" timegap %f \n",orbitalinfo->TimeGap);
1696     Float_t Bx = -orbitalinfo->Bdown;
1697     Float_t By = orbitalinfo->Beast;
1698     Float_t Bz = orbitalinfo->Bnorth;
1699    
1700     // TMatrixD Qiji(3,3);
1701     TMatrixD Qij = PO->QuatoECI(orbitalinfo->q0,orbitalinfo->q1,orbitalinfo->q2,orbitalinfo->q3);
1702     TMatrixD Dij = PO->ECItoGEO(Qij,orbitalinfo->absTime,orbitalinfo->lat,orbitalinfo->lon);
1703    
1704     //10REDNEW
1705     // If initial orientation data have reason to be inaccurate
1706     Float_t tg = 0.00;
1707     Float_t tmptg;
1708     Bool_t tgpar=false;
1709     Bool_t tgpar0=false;
1710     if (orbitalinfo->TimeGap>10.0 || ((modf(orbitalinfo->TimeGap,&tmptg)*1000>10 || modf(orbitalinfo->TimeGap,&tmptg)*1000==0.0) && orbitalinfo->TimeGap>2.0)) tgpar=true;
1711     if (orbitalinfo->TimeGap>180.0) tgpar0=true;
1712     if(MU!=0){
1713     // if(orbitalinfo->TimeGap>0 && errq==0 && azim==0){ // 10RED CHECK (comparison between three metod of recovering orientation)
1714     if((atime>=RTstart[MU] && atime<RTstart[MU+1] && RTbank1[MU]==0 && RTbank2[MU]==0 && (TMath::Abs(orbitalinfo->etha)>0.1 || tgpar0)) || ((RTbank1[MU]!=0 || RTbank2[MU]!=0) && atime>=RTstart[MU] && atime<RTstart[MU+1] && azim==0 && (errq!=0 || tgpar))){
1715     //found in Rotation Table this data for this time interval
1716     if(atime<RTtime1[0])
1717     orbitalinfo->azim = 5; //means that RotationTable no started yet
1718     else{
1719     // search for angle betwean velosity and direction to north in tangential to Earth surfase plane in satellite position
1720     Double_t bank=RTstart[MU];
1721     Double_t tlat=orbitalinfo->lat;
1722    
1723     tg=modf(orbitalinfo->TimeGap,&tg)*1000;
1724    
1725     if(atime>=RTpluto1[MU] && atime<=RTpluto2[MU]){
1726     Double_t kar=(RTbank2[MU]-RTbank1[MU])/(RTtime2[MU]-RTtime1[MU]);
1727     Double_t bak=RTbank1[MU]-kar*RTtime1[MU];
1728     bank=kar*atime+bak;
1729     }
1730     if(atime>=RTstart[MU] && atime<RTpluto1[MU]){
1731     Double_t s_dBdt2=(RTbpluto1[MU]-RTbank1[MU])/(Int_t)(RTpluto1[MU]-RTstart[MU]);
1732     Double_t s_t2=((Double_t)RTpluto1[MU]+(Double_t)RTstart[MU])/2. - RTstart[MU];
1733     Double_t s_t1=RTstart[MU]-RTstart[MU];
1734     Double_t s_k=s_dBdt2/(s_t2-s_t1);
1735     Double_t s_b=-s_k*s_t1;
1736     Double_t s_t3=RTpluto1[MU]-RTstart[MU];
1737     Double_t s_b3=RTbpluto1[MU];
1738     Double_t s_c=s_b3-0.5*s_k*s_t3*s_t3 -s_b*s_t3;
1739     bank=0.5*s_k*(atime-RTstart[MU])*(atime-RTstart[MU]) + s_b*(atime-RTstart[MU]) + s_c;
1740     }
1741     if(atime>RTpluto2[MU] && atime<=RTstart[MU+1]){
1742     Double_t s_dBdt2=(RTbpluto2[MU] - RTbank2[MU])/(Int_t)(RTpluto2[MU]-RTstart[MU+1]);
1743     Double_t s_t2=((Double_t)RTpluto2[MU]+(Double_t)RTstart[MU+1])/2. - RTstart[MU];
1744     Double_t s_t1=RTstart[MU+1]-RTstart[MU];
1745     Double_t s_k=s_dBdt2/(s_t2-s_t1);
1746     Double_t s_b=-s_k*s_t1;
1747     Double_t s_t3=RTpluto2[MU]-RTstart[MU];
1748     Double_t s_b3=RTbpluto2[MU];
1749     Double_t s_c=s_b3-0.5*s_k*s_t3*s_t3 -s_b*s_t3;
1750     bank=0.5*s_k*(atime-RTstart[MU])*(atime-RTstart[MU]) + s_b*(atime-RTstart[MU]) + s_c;
1751     }
1752     if(TMath::Abs(orbitalinfo->etha-bank)>0.1){
1753     orbitalinfo->etha=bank;
1754     Double_t spitch = 0.00001; // temprary not zero to avoid problem with tranzition from Euler angles to orientation matrix
1755    
1756     //Estimations of pitch angle of satellite
1757     if(TMath::Abs(bank)>0.7){
1758     Float_t spitch1=TMath::DegToRad()*0.7*diro;//RTdir1[MU];
1759     Float_t spitch2=TMath::DegToRad()*0.7*diro;//RTdir2[MU];
1760     Float_t kva=(spitch2-spitch1)/(RTtime2[MU]-RTtime1[MU]);
1761     Float_t bva=spitch1-kva*RTtime1[MU];
1762     spitch=kva*atime+bva;
1763     }
1764    
1765     //Calculate Yaw angle accordingly with fit, see picture FitYaw.jpg
1766     Double_t yaw=0.00001; // temprary not zero to avoid problem with tranzition from Euler angles to orientation matrix
1767     if(TMath::Abs(tlat)<70)
1768     yaw = -3.7e-8*tlat*tlat*tlat*tlat + 1.4e-7*tlat*tlat*tlat - 0.0005*tlat*tlat - 0.00025*tlat + 3.6;
1769     yaw = diro*yaw; //because should be different sign for ascending and descending orbits!
1770     orbitalinfo->phi=yaw;
1771    
1772     if(TMath::Abs(bank)>0.5 && TMath::Abs(yaw-orbitalinfo->phi)<3.0) yaw=orbitalinfo->phi;
1773    
1774     //Qiji = PO->EulertoEci(eCi.getPos().m_x,eCi.getPos().m_y,eCi.getPos().m_z,eCi.getVel().m_x,eCi.getVel().m_y,eCi.getVel().m_z,bank,yaw,spitch); // 10RED CHECK
1775     Qij = PO->EulertoEci(eCi.getPos().m_x,eCi.getPos().m_y,eCi.getPos().m_z,eCi.getVel().m_x,eCi.getVel().m_y,eCi.getVel().m_z,bank,yaw,spitch); // STANDARD
1776     orbitalinfo->qkind = 1;
1777     }
1778     //Qij = PO->GEOtoECI(Dij,orbitalinfo->absTime,orbitalinfo->lat,orbitalinfo->lon); // to convert from Dij to Qij
1779     } // end of if(atime<RTtime1[0]
1780     } // end of (((orbitalinfo->TimeGap>60.0 && TMath...
1781     } // end of MU~=0
1782    
1783     TMatrixD qij = PO->ColPermutation(Qij);
1784     TMatrixD Fij = PO->ECItoGreenwich(Qij,orbitalinfo->absTime);
1785     TMatrixD Gij = PO->ColPermutation(Fij);
1786     Dij = PO->ECItoGEO(Qij,orbitalinfo->absTime,orbitalinfo->lat,orbitalinfo->lon);
1787     TMatrixD Iij = PO->ColPermutation(Dij);
1788    
1789     TVector3 SP = PO->GetSunPosition(orbitalinfo->absTime);
1790     // go to Pamela reference frame from Resurs reference frame
1791     Float_t tmpy = SP.Y();
1792     SP.SetY(SP.Z());
1793     SP.SetZ(-tmpy);
1794     TVector3 SunZenith;
1795     SunZenith.SetMagThetaPhi(1,23.439281*TMath::DegToRad(),TMath::Pi()/2.);
1796     TVector3 SunMag = -SP;
1797     SunMag.Rotate(-45*TMath::DegToRad(),SunZenith);
1798     tmpy=SunMag.Y();
1799     SunMag.SetY(SunMag.Z());
1800     SunMag.SetZ(-tmpy);
1801    
1802     orbitalinfo->Iij.ResizeTo(Iij);
1803     orbitalinfo->Iij = Iij;
1804    
1805     Bool_t saso=true;
1806     if (orbitalinfo->qkind==1) saso=true;
1807     if (orbitalinfo->qkind==0 && orbitalinfo->azim>0 && orbitalinfo->azim!=5 && tgpar) saso=false;
1808     if (orbitalinfo->qkind==0 && orbitalinfo->azim==5 && TMath::Abs(orbitalinfo->etha)>0.1 && tgpar) saso=false;
1809     if (orbitalinfo->qkind==0 && orbitalinfo->azim==5 && TMath::Abs(orbitalinfo->etha)<=0.1 && tgpar0) saso=false;
1810     if (saso) orbitalinfo->mode=orbitalinfo->rtqual; else orbitalinfo->mode=2;
1811    
1812     //
1813     // A1 = Iij(0,2);
1814     // A2 = Iij(1,2);
1815     // A3 = Iij(2,2);
1816     //
1817     // orbitalinfo->pamzenitangle = (Float_t)PO->GetPitchAngle(1,0,0,A1,A2,A3); // Angle between zenit and Pamela's main axiz
1818     // orbitalinfo->pamBangle = (Float_t)PO->GetPitchAngle(A1,A2,A3,Bx,By,Bz); // Angle between Pamela's main axiz and B
1819     //
1820 mocchiut 1.65 if ( debug ) printf(" matrixes done \n");
1821 mayorov 1.93 if ( !standalone ){
1822 mocchiut 1.65 if ( debug ) printf(" !standalone \n");
1823 mayorov 1.93 //
1824 mocchiut 1.73 // Standard tracking algorithm
1825     //
1826 mayorov 1.93 Int_t nn = 0;
1827 mocchiut 1.73 if ( verbose ) printf(" standard tracking \n");
1828 mayorov 1.93 for(Int_t nt=0; nt < tof->ntrk(); nt++){
1829     //
1830     ToFTrkVar *ptt = tof->GetToFTrkVar(nt);
1831     if (debug) cout<<"tof->ntrk() = "<<tof->ntrk()<<"\tptt->trkseqno = "<<ptt->trkseqno<<"\ttrke->ntrk() = "<<trke->ntrk()<<endl;
1832     Double_t E11x = ptt->xtr_tof[0]; // tr->x[0];
1833     Double_t E11y = ptt->ytr_tof[0]; //tr->y[0];
1834     Double_t E11z = zin[0];
1835     Double_t E22x = ptt->xtr_tof[3];//tr->x[3];
1836     Double_t E22y = ptt->ytr_tof[3];//tr->y[3];
1837     Double_t E22z = zin[3];
1838     if ( (E11x < 100. && E11y < 100. && E22x < 100. && E22y < 100.) || ptt->trkseqno != -1 ){
1839     TrkTrack *mytrack = trke->GetStoredTrack(ptt->trkseqno);
1840     Float_t rig=1/mytrack->GetDeflection();
1841     Double_t norm = sqrt(pow(E22x-E11x,2)+pow(E22y-E11y,2)+pow(E22z-E11z,2));
1842 mocchiut 1.73 //
1843 mayorov 1.93 Px = (E22x-E11x)/norm;
1844     Py = (E22y-E11y)/norm;
1845     Pz = (E22z-E11z)/norm;
1846     //
1847     t_orb->trkseqno = ptt->trkseqno;
1848     //
1849     TMatrixD Eij = PO->PamelatoGEO(Iij,Px,Py,Pz);
1850     t_orb->Eij.ResizeTo(Eij);
1851     t_orb->Eij = Eij;
1852     //
1853     TMatrixD Sij = PO->PamelatoGEO(Gij,Px,Py,Pz);
1854     t_orb->Sij.ResizeTo(Sij);
1855     t_orb->Sij = Sij;
1856     //
1857     t_orb->pitch = (Float_t)PO->GetPitchAngle(Eij(0,0),Eij(1,0),Eij(2,0),Bx,By,Bz);
1858     //
1859 mocchiut 1.73 // SunPosition in instrumental reference frame
1860 mayorov 1.93 TMatrixD Kij = PO->PamelatoGEO(qij,Px,Py,Pz);
1861     TMatrixD Lij = PO->PamelatoGEO(qij,0,0,1);
1862     t_orb->sunangle=(Float_t)PO->GetPitchAngle(Kij(0,0),Kij(1,0),Kij(2,0),-SP.X(),-SP.Y(),-SP.Z());
1863     t_orb->sunmagangle=(Float_t)PO->GetPitchAngle(Kij(0,0),Kij(1,0),Kij(2,0),SunMag.X(),SunMag.Y(),SunMag.Z());
1864     //
1865     //
1866     Double_t omega = PO->GetPitchAngle(-Eij(0,0),-Eij(1,0),-Eij(2,0),1,0,0) * TMath::DegToRad();
1867     TVector3 Bxy(0,By,Bz);
1868     TVector3 Exy(0,-Eij(1,0),-Eij(2,0));
1869     Double_t dzeta=Bxy.Angle(Exy);
1870     if (-Eij(1,0) < 0) dzeta=2.0*TMath::Pi() - dzeta;
1871 mocchiut 1.73
1872     if(debug) cout << "omega = "<<omega*TMath::RadToDeg()<<"\tdzeta = "<<dzeta*TMath::RadToDeg()<<endl;
1873 pam-mep 1.69
1874 mocchiut 1.73 // Formula from D.F. Smart *, M.A. Shea [2005]; A review of geomagnetic cutoff rigidities for earth-orbiting spacecraft
1875 mayorov 1.93 if(rig>=0) t_orb->cutoff = 59.3/(pow(orbitalinfo->L,2)*pow(1+sqrt(1-sin(omega)*sin(dzeta)*pow(cos(orbitalinfo->lat*TMath::DegToRad()),3)),2));
1876     else t_orb->cutoff = 59.3/(pow(orbitalinfo->L,2)*pow(1+sqrt(1-sin(omega)*sin(TMath::Pi()+dzeta)*pow(cos(orbitalinfo->lat*TMath::DegToRad()),3)),2));
1877     if (debug) cout << "R = " << rig << "\tcutoff = " << t_orb->cutoff << endl;
1878    
1879     //
1880     if ( t_orb->pitch != t_orb->pitch ) t_orb->pitch = -1000.;
1881     if ( t_orb->cutoff != t_orb->cutoff ) t_orb->cutoff = -1000.;
1882     if ( t_orb->sunangle != t_orb->sunangle ) t_orb->sunangle = -1000.;
1883     if ( t_orb->sunmagangle != t_orb->sunmagangle ) t_orb->sunmagangle = -1000.;
1884     //
1885     if ( debug ) printf(" orbitalinfo->cutoffsvl %f vitaly %f \n",orbitalinfo->cutoffsvl,t_orb->cutoff);
1886     //
1887     new(tor[nn]) OrbitalInfoTrkVar(*t_orb);
1888     nn++;
1889     //
1890     t_orb->Clear();
1891     //
1892     }
1893     //
1894     } // end standard tracking algorithm
1895 mocchiut 1.73
1896     //
1897     // Code for extended tracking algorithm:
1898     //
1899     if ( hasNucleiTrk ){
1900     Int_t ttentry = 0;
1901     if ( verbose ) printf(" hasNucleiTrk \n");
1902     for(Int_t nt=0; nt < tcNucleiTof->GetEntries() ; nt++){
1903     //
1904     if ( verbose ) printf(" got1\n");
1905     ToFTrkVar *ptt = (ToFTrkVar*)(tcNucleiTof->At(nt));
1906     if (verbose) cout<<" tcNucleiTof->GetEntries() = "<<tcNucleiTof->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     TrkTrack *mytrack = (TrkTrack*)(tcNucleiTrk->At(ptt->trkseqno));
1915     if ( verbose ) printf(" got tcNucleiTrk \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 mayorov 1.93 t_orb->Eij.ResizeTo(Eij);
1927     t_orb->Eij = Eij;
1928 mocchiut 1.73 //
1929     TMatrixD Sij = PO->PamelatoGEO(Gij,Px,Py,Pz);
1930     t_orb->Sij.ResizeTo(Sij);
1931     t_orb->Sij = Sij;
1932 mayorov 1.93 //
1933 mocchiut 1.73 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 &tt1 = *torbNucleiTrk;
1964     new(tt1[ttentry]) OrbitalInfoTrkVar(*t_orb);
1965     ttentry++;
1966     //
1967     t_orb->Clear();
1968     //
1969     }
1970     //
1971     }
1972     } // end standard tracking algorithm: nuclei
1973     if ( hasExtNucleiTrk ){
1974     Int_t ttentry = 0;
1975     if ( verbose ) printf(" hasExtNucleiTrk \n");
1976     for(Int_t nt=0; nt < tcExtNucleiTof->GetEntries() ; nt++){
1977     //
1978     if ( verbose ) printf(" got2\n");
1979     ToFTrkVar *ptt = (ToFTrkVar*)(tcExtNucleiTof->At(nt));
1980     if (verbose) cout<<" tcExtNucleiTof->GetEntries() = "<<tcExtNucleiTof->GetEntries()<<"\tptt->trkseqno = "<<ptt->trkseqno<<endl;
1981     Double_t E11x = ptt->xtr_tof[0]; // tr->x[0];
1982     Double_t E11y = ptt->ytr_tof[0]; //tr->y[0];
1983     Double_t E11z = zin[0];
1984     Double_t E22x = ptt->xtr_tof[3];//tr->x[3];
1985     Double_t E22y = ptt->ytr_tof[3];//tr->y[3];
1986     Double_t E22z = zin[3];
1987     if ( (E11x < 100. && E11y < 100. && E22x < 100. && E22y < 100.) || ptt->trkseqno != -1 ){
1988     ExtTrack *mytrack = (ExtTrack*)(tcExtNucleiTrk->At(ptt->trkseqno));
1989     if ( verbose ) printf(" got tcExtNucleiTrk \n");
1990     Float_t rig=1/mytrack->GetDeflection();
1991     Double_t norm = sqrt(pow(E22x-E11x,2)+pow(E22y-E11y,2)+pow(E22z-E11z,2));
1992     //
1993     Px = (E22x-E11x)/norm;
1994     Py = (E22y-E11y)/norm;
1995     Pz = (E22z-E11z)/norm;
1996     //
1997     t_orb->trkseqno = ptt->trkseqno;
1998     //
1999     TMatrixD Eij = PO->PamelatoGEO(Iij,Px,Py,Pz);
2000 mayorov 1.93 t_orb->Eij.ResizeTo(Eij);
2001     t_orb->Eij = Eij;
2002 mocchiut 1.73 //
2003     TMatrixD Sij = PO->PamelatoGEO(Gij,Px,Py,Pz);
2004     t_orb->Sij.ResizeTo(Sij);
2005     t_orb->Sij = Sij;
2006 mayorov 1.93 //
2007 mocchiut 1.73 t_orb->pitch = (Float_t)PO->GetPitchAngle(Eij(0,0),Eij(1,0),Eij(2,0),Bx,By,Bz);
2008     //
2009     // SunPosition in instrumental reference frame
2010     TMatrixD Kij = PO->PamelatoGEO(qij,Px,Py,Pz);
2011     TMatrixD Lij = PO->PamelatoGEO(qij,0,0,1);
2012     t_orb->sunangle=(Float_t)PO->GetPitchAngle(Kij(0,0),Kij(1,0),Kij(2,0),-SP.X(),-SP.Y(),-SP.Z());
2013     t_orb->sunmagangle=(Float_t)PO->GetPitchAngle(Kij(0,0),Kij(1,0),Kij(2,0),SunMag.X(),SunMag.Y(),SunMag.Z());
2014     //
2015     //
2016     Double_t omega = PO->GetPitchAngle(-Eij(0,0),-Eij(1,0),-Eij(2,0),1,0,0) * TMath::DegToRad();
2017     TVector3 Bxy(0,By,Bz);
2018     TVector3 Exy(0,-Eij(1,0),-Eij(2,0));
2019     Double_t dzeta=Bxy.Angle(Exy);
2020     if (-Eij(1,0) < 0) dzeta=2.0*TMath::Pi() - dzeta;
2021    
2022     if(debug) cout << "omega = "<<omega*TMath::RadToDeg()<<"\tdzeta = "<<dzeta*TMath::RadToDeg()<<endl;
2023    
2024     // Formula from D.F. Smart *, M.A. Shea [2005]; A review of geomagnetic cutoff rigidities for earth-orbiting spacecraft
2025     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));
2026     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));
2027     if (debug) cout << "R = " << rig << "\tcutoff = " << t_orb->cutoff << endl;
2028    
2029     //
2030     if ( t_orb->pitch != t_orb->pitch ) t_orb->pitch = -1000.;
2031     if ( t_orb->cutoff != t_orb->cutoff ) t_orb->cutoff = -1000.;
2032     if ( t_orb->sunangle != t_orb->sunangle ) t_orb->sunangle = -1000.;
2033     if ( t_orb->sunmagangle != t_orb->sunmagangle ) t_orb->sunmagangle = -1000.;
2034     //
2035     if ( debug ) printf(" orbitalinfo->cutoffsvl %f vitaly %f \n",orbitalinfo->cutoffsvl,t_orb->cutoff);
2036     //
2037     TClonesArray &tt2 = *torbExtNucleiTrk;
2038     new(tt2[ttentry]) OrbitalInfoTrkVar(*t_orb);
2039     ttentry++;
2040     //
2041     t_orb->Clear();
2042     //
2043     }
2044     //
2045     }
2046     } // end standard tracking algorithm: nuclei extended
2047     if ( hasExtTrk ){
2048     Int_t ttentry = 0;
2049     if ( verbose ) printf(" hasExtTrk \n");
2050     for(Int_t nt=0; nt < tcExtTof->GetEntries() ; nt++){
2051     //
2052     if ( verbose ) printf(" got3\n");
2053     ToFTrkVar *ptt = (ToFTrkVar*)(tcExtTof->At(nt));
2054     if (verbose) cout<<" tcExtTof->GetEntries() = "<<tcExtTof->GetEntries()<<"\tptt->trkseqno = "<<ptt->trkseqno<<endl;
2055     Double_t E11x = ptt->xtr_tof[0]; // tr->x[0];
2056     Double_t E11y = ptt->ytr_tof[0]; //tr->y[0];
2057     Double_t E11z = zin[0];
2058     Double_t E22x = ptt->xtr_tof[3];//tr->x[3];
2059     Double_t E22y = ptt->ytr_tof[3];//tr->y[3];
2060     Double_t E22z = zin[3];
2061     if ( (E11x < 100. && E11y < 100. && E22x < 100. && E22y < 100.) || ptt->trkseqno != -1 ){
2062     ExtTrack *mytrack = (ExtTrack*)(tcExtTrk->At(ptt->trkseqno));
2063     if ( verbose ) printf(" got tcExtTrk \n");
2064     Float_t rig=1/mytrack->GetDeflection();
2065     Double_t norm = sqrt(pow(E22x-E11x,2)+pow(E22y-E11y,2)+pow(E22z-E11z,2));
2066     //
2067     Px = (E22x-E11x)/norm;
2068     Py = (E22y-E11y)/norm;
2069     Pz = (E22z-E11z)/norm;
2070     //
2071     t_orb->trkseqno = ptt->trkseqno;
2072     //
2073     TMatrixD Eij = PO->PamelatoGEO(Iij,Px,Py,Pz);
2074 mayorov 1.93 t_orb->Eij.ResizeTo(Eij);
2075     t_orb->Eij = Eij;
2076 mocchiut 1.73 //
2077     TMatrixD Sij = PO->PamelatoGEO(Gij,Px,Py,Pz);
2078     t_orb->Sij.ResizeTo(Sij);
2079     t_orb->Sij = Sij;
2080 mayorov 1.93 //
2081 mocchiut 1.73 t_orb->pitch = (Float_t)PO->GetPitchAngle(Eij(0,0),Eij(1,0),Eij(2,0),Bx,By,Bz);
2082     //
2083     // SunPosition in instrumental reference frame
2084     TMatrixD Kij = PO->PamelatoGEO(qij,Px,Py,Pz);
2085     TMatrixD Lij = PO->PamelatoGEO(qij,0,0,1);
2086     t_orb->sunangle=(Float_t)PO->GetPitchAngle(Kij(0,0),Kij(1,0),Kij(2,0),-SP.X(),-SP.Y(),-SP.Z());
2087     t_orb->sunmagangle=(Float_t)PO->GetPitchAngle(Kij(0,0),Kij(1,0),Kij(2,0),SunMag.X(),SunMag.Y(),SunMag.Z());
2088     //
2089     //
2090     Double_t omega = PO->GetPitchAngle(-Eij(0,0),-Eij(1,0),-Eij(2,0),1,0,0) * TMath::DegToRad();
2091     TVector3 Bxy(0,By,Bz);
2092     TVector3 Exy(0,-Eij(1,0),-Eij(2,0));
2093     Double_t dzeta=Bxy.Angle(Exy);
2094     if (-Eij(1,0) < 0) dzeta=2.0*TMath::Pi() - dzeta;
2095    
2096     if(debug) cout << "omega = "<<omega*TMath::RadToDeg()<<"\tdzeta = "<<dzeta*TMath::RadToDeg()<<endl;
2097    
2098     // Formula from D.F. Smart *, M.A. Shea [2005]; A review of geomagnetic cutoff rigidities for earth-orbiting spacecraft
2099     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));
2100     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));
2101     if (debug) cout << "R = " << rig << "\tcutoff = " << t_orb->cutoff << endl;
2102    
2103     //
2104     if ( t_orb->pitch != t_orb->pitch ) t_orb->pitch = -1000.;
2105     if ( t_orb->cutoff != t_orb->cutoff ) t_orb->cutoff = -1000.;
2106     if ( t_orb->sunangle != t_orb->sunangle ) t_orb->sunangle = -1000.;
2107     if ( t_orb->sunmagangle != t_orb->sunmagangle ) t_orb->sunmagangle = -1000.;
2108     //
2109     if ( debug ) printf(" orbitalinfo->cutoffsvl %f vitaly %f \n",orbitalinfo->cutoffsvl,t_orb->cutoff);
2110     //
2111     TClonesArray &tt3 = *torbExtTrk;
2112     new(tt3[ttentry]) OrbitalInfoTrkVar(*t_orb);
2113     ttentry++;
2114     //
2115     t_orb->Clear();
2116     //
2117     }
2118     //
2119     }
2120     } // end standard tracking algorithm: extended
2121    
2122 mayorov 1.93 } else {
2123     if ( debug ) printf(" mmm... mode %u standalone \n",orbitalinfo->mode);
2124     }
2125     //
2126 mocchiut 1.73 } else { // HERE IT MISS ALL CODE FOR EXTENDED TRACKING!
2127 mayorov 1.93 if ( !standalone ){
2128     //
2129 mocchiut 1.73 if ( verbose ) printf(" no orb info for tracks \n");
2130 mayorov 1.93 Int_t nn = 0;
2131     for(Int_t nt=0; nt < tof->ntrk(); nt++){
2132     //
2133     ToFTrkVar *ptt = tof->GetToFTrkVar(nt);
2134     if ( ptt->trkseqno != -1 ){
2135     //
2136     t_orb->trkseqno = ptt->trkseqno;
2137     //
2138     TMatrixD Iij(3,1);
2139     Iij(0,0)=0.; Iij(1,0)=0.; Iij(2,0)=1.;
2140     //Iij(1,0)=0; Iij(1,1)=0; Iij(1,2)=0;
2141     //Iij(2,0)=0; Iij(2,1)=0; Iij(2,2)=0;
2142     //Iij.Zero();
2143     t_orb->Eij.ResizeTo(Iij);
2144     t_orb->Sij.ResizeTo(Iij);
2145     t_orb->Eij = Iij;
2146     //
2147     t_orb->Sij = Iij;
2148     //
2149     t_orb->pitch = -1000.;
2150     //
2151     t_orb->sunangle = -1000.;
2152     //
2153     t_orb->sunmagangle = -1000;
2154     //
2155     t_orb->cutoff = -1000.;
2156     //
2157     new(tor[nn]) OrbitalInfoTrkVar(*t_orb);
2158     nn++;
2159     //
2160     t_orb->Clear();
2161     //
2162     }
2163     //
2164     }
2165 mocchiut 1.73 //
2166     // Code for extended tracking algorithm:
2167     //
2168     if ( hasNucleiTrk ){
2169     Int_t ttentry = 0;
2170     if ( verbose ) printf(" hasNucleiTrk \n");
2171     for(Int_t nt=0; nt < tcNucleiTof->GetEntries() ; nt++){
2172     //
2173     ToFTrkVar *ptt = (ToFTrkVar*)(tcNucleiTof->At(nt));
2174     if ( ptt->trkseqno != -1 ){
2175     //
2176     t_orb->trkseqno = ptt->trkseqno;
2177     //
2178 mayorov 1.93 TMatrixD Iij(3,1);
2179     Iij(0,0)=0.; Iij(1,0)=0.; Iij(2,0)=1.;
2180     //Iij(1,0)=0; Iij(1,1)=0; Iij(1,2)=0;
2181     //Iij(2,0)=0; Iij(2,1)=0; Iij(2,2)=0;
2182     //Iij.Zero();
2183     t_orb->Eij.ResizeTo(Iij);
2184     t_orb->Sij.ResizeTo(Iij);
2185     t_orb->Eij = Iij;
2186 mocchiut 1.73 //
2187 mayorov 1.93 t_orb->Sij = Iij;
2188     //
2189 mocchiut 1.73 t_orb->pitch = -1000.;
2190     //
2191     t_orb->sunangle = -1000.;
2192     //
2193     t_orb->sunmagangle = -1000;
2194     //
2195     t_orb->cutoff = -1000.;
2196     //
2197     TClonesArray &tz1 = *torbNucleiTrk;
2198     new(tz1[ttentry]) OrbitalInfoTrkVar(*t_orb);
2199     ttentry++;
2200     //
2201     t_orb->Clear();
2202     //
2203     }
2204     //
2205     }
2206     }
2207     if ( hasExtNucleiTrk ){
2208     Int_t ttentry = 0;
2209     if ( verbose ) printf(" hasExtNucleiTrk \n");
2210     for(Int_t nt=0; nt < tcExtNucleiTof->GetEntries() ; nt++){
2211     //
2212     if ( verbose ) printf(" got2\n");
2213     ToFTrkVar *ptt = (ToFTrkVar*)(tcExtNucleiTof->At(nt));
2214     if ( ptt->trkseqno != -1 ){
2215     //
2216     t_orb->trkseqno = ptt->trkseqno;
2217     //
2218 mayorov 1.93 TMatrixD Iij(3,1);
2219     Iij(0,0)=0.; Iij(1,0)=0.; Iij(2,0)=1.;
2220     //Iij(1,0)=0; Iij(1,1)=0; Iij(1,2)=0;
2221     //Iij(2,0)=0; Iij(2,1)=0; Iij(2,2)=0;
2222     //Iij.Zero();
2223     t_orb->Eij.ResizeTo(Iij);
2224     t_orb->Sij.ResizeTo(Iij);
2225     t_orb->Eij = Iij;
2226 mocchiut 1.73 //
2227 mayorov 1.93 t_orb->Sij = Iij;
2228     //
2229 mocchiut 1.73 t_orb->pitch = -1000.;
2230     //
2231     t_orb->sunangle = -1000.;
2232     //
2233     t_orb->sunmagangle = -1000;
2234     //
2235     t_orb->cutoff = -1000.;
2236     //
2237     TClonesArray &tz2 = *torbExtNucleiTrk;
2238     new(tz2[ttentry]) OrbitalInfoTrkVar(*t_orb);
2239     ttentry++;
2240     //
2241     t_orb->Clear();
2242     //
2243     }
2244     //
2245     }
2246     }
2247     if ( hasExtTrk ){
2248 mayorov 1.93 Int_t ttentry = 0;
2249 mocchiut 1.73 if ( verbose ) printf(" hasExtTrk \n");
2250     for(Int_t nt=0; nt < tcExtTof->GetEntries() ; nt++){
2251     //
2252     if ( verbose ) printf(" got3\n");
2253     ToFTrkVar *ptt = (ToFTrkVar*)(tcExtTof->At(nt));
2254     if ( ptt->trkseqno != -1 ){
2255     //
2256     t_orb->trkseqno = ptt->trkseqno;
2257     //
2258 mayorov 1.93 TMatrixD Iij(3,1);
2259     Iij(0,0)=0.; Iij(1,0)=0.; Iij(2,0)=1.;
2260     //Iij(1,0)=0; Iij(1,1)=0; Iij(1,2)=0;
2261     //Iij(2,0)=0; Iij(2,1)=0; Iij(2,2)=0;
2262     //Iij.Zero();
2263     t_orb->Eij.ResizeTo(Iij);
2264     t_orb->Sij.ResizeTo(Iij);
2265     t_orb->Eij = Iij;
2266 mocchiut 1.73 //
2267 mayorov 1.93 t_orb->Sij = Iij;
2268     //
2269 mocchiut 1.73 t_orb->pitch = -1000.;
2270     //
2271     t_orb->sunangle = -1000.;
2272     //
2273     t_orb->sunmagangle = -1000;
2274     //
2275     t_orb->cutoff = -1000.;
2276     //
2277 mayorov 1.93 TClonesArray &tz3 = *torbExtTrk;
2278 mocchiut 1.73 new(tz3[ttentry]) OrbitalInfoTrkVar(*t_orb);
2279     ttentry++;
2280     //
2281     t_orb->Clear();
2282     //
2283     }
2284     //
2285     }
2286     }
2287 mayorov 1.93 }
2288     } // if( orbitalinfo->TimeGap>0){
2289 mocchiut 1.30 //
2290 mocchiut 1.15 // Fill the class
2291 pamelaprod 1.11 //
2292 mocchiut 1.1 OrbitalInfotr->Fill();
2293 mocchiut 1.15 //
2294 mocchiut 1.81 // tor.Clear("C"); // memory leak?
2295     tor.Delete(); // memory leak?
2296 mocchiut 1.30 delete t_orb;
2297     //
2298 mocchiut 1.81 // printf(" q0 size %i q0 capacity %i \n",(int)q0.size(),(int)q0.capacity());
2299     } // loop over the events in the run
2300    
2301    
2302 mocchiut 1.1 //
2303     // Here you may want to clear some variables before processing another run
2304     //
2305 mocchiut 1.44
2306 mocchiut 1.81 // OrbitalInfotr->FlushBaskets();
2307    
2308 mocchiut 1.57 if ( verbose ) printf(" Clear before new run \n");
2309 mocchiut 1.5 delete dbtime;
2310 mocchiut 1.57
2311 mocchiut 1.62 if ( mcmdrc ) mcmdrc->Clear();
2312 mocchiut 1.57 mcmdrc = 0;
2313    
2314     if ( verbose ) printf(" Clear before new run1 \n");
2315     if ( L_QQ_Q_l_lower ) delete L_QQ_Q_l_lower;
2316     if ( verbose ) printf(" Clear before new run2 \n");
2317 mocchiut 1.18 if ( L_QQ_Q_l_upper ) delete L_QQ_Q_l_upper;
2318 mocchiut 1.57 if ( verbose ) printf(" Clear before new run3 \n");
2319 mocchiut 1.18 if ( RYPang_upper ) delete RYPang_upper;
2320 mocchiut 1.57 if ( verbose ) printf(" Clear before new run4 \n");
2321 mocchiut 1.18 if ( RYPang_lower ) delete RYPang_lower;
2322 mocchiut 1.57
2323 mocchiut 1.81
2324     if ( l0tr ){
2325     if ( verbose ) printf(" delete l0tr\n");
2326     l0tr->Delete();
2327     l0tr = 0;
2328     }
2329     // if ( l0head ){
2330     // printf(" delete l0head\n");
2331     // l0head->Reset();
2332     // delete l0head;
2333     // printf(" delete l0head done\n");
2334     // l0head = 0;
2335     // }
2336     if (eh) {
2337     if ( verbose ) printf(" delete eh\n");
2338     delete eh;
2339     eh = 0;
2340     }
2341    
2342     if ( verbose ) printf(" close file \n");
2343     if ( l0File ) l0File->Close("R");
2344 mocchiut 1.57 if ( verbose ) printf(" End run \n");
2345    
2346 mocchiut 1.81 q0.clear();
2347     q1.clear();
2348     q2.clear();
2349     q3.clear();
2350     qtime.clear();
2351     qPitch.clear();
2352     qRoll.clear();
2353     qYaw.clear();
2354     qmode.clear();
2355    
2356     if (ch){
2357     if ( verbose ) printf(" delete ch\n");
2358     ch->Delete();
2359     ch = 0;
2360     }
2361     } // process all the runs <===
2362     if ( debug ){
2363     printf("AFTER LOOP ON RUNs\n");
2364     gObjectTable->Print();
2365     }
2366     //
2367 mocchiut 1.1 if (verbose) printf("\n Finished processing data \n");
2368     //
2369     closeandexit:
2370     //
2371     // we have finished processing the run(s). If we processed a single run now we must copy all the events after our run from the old tree to the new one and delete the old tree.
2372     //
2373     if ( !reprocall && reproc && code >= 0 ){
2374     if ( totfileentries > noaftrun ){
2375     if (verbose){
2376 mayorov 1.93 printf("\n Post-processing: copying events from the old tree after the processed run\n");
2377     printf(" Copying %i events in the file which are after the end of the run %i \n",(int)(totfileentries-noaftrun),(int)run);
2378     printf(" Start copying at event number %i end copying at event number %i \n",(int)noaftrun,(int)totfileentries);
2379 mocchiut 1.1 }
2380     for (UInt_t j = noaftrun; j < totfileentries; j++ ){
2381 mayorov 1.93 //
2382     // Get entry from old tree
2383     //
2384     if ( OrbitalInfotrclone->GetEntry(j) <= 0 ) throw -36;
2385     //
2386     // copy orbitalinfoclone to OrbitalInfo
2387     //
2388     // orbitalinfo->Clear();
2389     //
2390     memcpy(&orbitalinfo,&orbitalinfoclone,sizeof(orbitalinfoclone));
2391     //
2392     // Fill entry in the new tree
2393     //
2394     OrbitalInfotr->Fill();
2395 mocchiut 1.1 };
2396     if (verbose) printf(" Finished successful copying!\n");
2397     };
2398 mocchiut 1.57 //if ( OrbitalInfotrclone ) OrbitalInfotrclone->Clear();
2399     //if ( OrbitalInfotrclone ) OrbitalInfotrclone->Delete();
2400 mocchiut 1.1 };
2401     //
2402     // Close files, delete old tree(s), write and close level2 file
2403     //
2404 pam-mep 1.69
2405 mocchiut 1.1 if ( l0File ) l0File->Close();
2406 mocchiut 1.23 if ( myfold ) gSystem->Unlink(tempname.str().c_str());
2407 mocchiut 1.5 //
2408 mocchiut 1.1 if ( OrbitalInfotr ) OrbitalInfotr->SetName("OrbitalInfo");
2409 mocchiut 1.30 //
2410 mocchiut 1.1 if ( file ){
2411     file->cd();
2412 mocchiut 1.63 if ( OrbitalInfotr ) OrbitalInfotr->Write("OrbitalInfo", TObject::kOverwrite); // 10 RED bug fixed
2413 mocchiut 1.1 };
2414     //
2415 mocchiut 1.57 if (verbose) printf("\n Exiting...\n");
2416    
2417 mocchiut 1.23 if ( myfold ) gSystem->Unlink(OrbitalInfofolder.str().c_str());
2418 mocchiut 1.1 //
2419     // the end
2420     //
2421 mocchiut 1.25 if ( dbc ){
2422     dbc->Close();
2423     delete dbc;
2424     };
2425 mocchiut 1.57 //
2426 mocchiut 1.1 if (verbose) printf("\n Exiting...\n");
2427 mocchiut 1.57 if ( tempfile ) tempfile->Close();
2428    
2429 mocchiut 1.30 if ( PO ) delete PO;
2430 mocchiut 1.57 if ( gltle ) delete gltle;
2431     if ( glparam ) delete glparam;
2432     if ( glparam2 ) delete glparam2;
2433     if (verbose) printf("\n Exiting3...\n");
2434 mocchiut 1.4 if ( glroot ) delete glroot;
2435 mocchiut 1.57 if (verbose) printf("\n Exiting4...\n");
2436     if ( runinfo ) runinfo->Close();
2437 mocchiut 1.4 if ( runinfo ) delete runinfo;
2438 mocchiut 1.58
2439 mocchiut 1.73 if ( tcNucleiTrk ){
2440     tcNucleiTrk->Delete();
2441     delete tcNucleiTrk;
2442     tcNucleiTrk = NULL;
2443     }
2444     if ( tcExtNucleiTrk ){
2445     tcExtNucleiTrk->Delete();
2446     delete tcExtNucleiTrk;
2447     tcExtNucleiTrk = NULL;
2448     }
2449     if ( tcExtTrk ){
2450     tcExtTrk->Delete();
2451     delete tcExtTrk;
2452     tcExtTrk = NULL;
2453     }
2454    
2455     if ( tcNucleiTof ){
2456     tcNucleiTof->Delete();
2457     delete tcNucleiTof;
2458     tcNucleiTof = NULL;
2459     }
2460     if ( tcExtNucleiTof ){
2461     tcExtNucleiTof->Delete();
2462     delete tcExtNucleiTof;
2463     tcExtNucleiTof = NULL;
2464     }
2465     if ( tcExtTof ){
2466     tcExtTof->Delete();
2467     delete tcExtTof;
2468     tcExtTrk = NULL;
2469     }
2470    
2471    
2472 pam-mep 1.69 if ( tof ) delete tof;
2473     if ( trke ) delete trke;
2474    
2475 mocchiut 1.58 if ( debug ){
2476 mocchiut 1.57 cout << "1 0x" << OrbitalInfotr << endl;
2477     cout << "2 0x" << OrbitalInfotrclone << endl;
2478     cout << "3 0x" << l0tr << endl;
2479     cout << "4 0x" << tempOrbitalInfo << endl;
2480     cout << "5 0x" << ttof << endl;
2481 mocchiut 1.58 }
2482 mocchiut 1.57 //
2483 mocchiut 1.58 if ( debug ) file->ls();
2484 mocchiut 1.4 //
2485 mocchiut 1.81 if ( debug ){
2486     printf("BEFORE EXITING\n");
2487     gObjectTable->Print();
2488     }
2489 mocchiut 1.1 if(code < 0) throw code;
2490     return(code);
2491     }
2492    
2493 mocchiut 1.7
2494     //
2495     // Returns the cCoordGeo structure holding the geographical
2496     // coordinates for the event (see sgp4.h).
2497     //
2498     // atime is the abstime of the event in UTC unix time.
2499     // tletime is the time of the tle in UTC unix time.
2500     // tle is the previous and nearest tle (compared to atime).
2501     cCoordGeo getCoo(UInt_t atime, UInt_t tletime, cTle *tle)
2502     {
2503     cEci eci;
2504     cOrbit orbit(*tle);
2505     orbit.getPosition((double) (atime - tletime)/60., &eci);
2506    
2507     return eci.toGeo();
2508     }
2509 mocchiut 1.15
2510     // function of copyng of quatrnions classes
2511    
2512     void CopyQ(Quaternions *Q1, Quaternions *Q2){
2513     for(UInt_t i = 0; i < 6; i++){
2514     Q1->time[i]=Q2->time[i];
2515     for (UInt_t j = 0; j < 4; j++)Q1->quat[i][j]=Q2->quat[i][j];
2516     }
2517     return;
2518     }
2519    
2520     // functions of copyng InclinationInfo classes
2521    
2522     void CopyAng(InclinationInfo *A1, InclinationInfo *A2){
2523     A1->Tangazh = A2->Tangazh;
2524     A1->Ryskanie = A2->Ryskanie;
2525     A1->Kren = A2->Kren;
2526     return;
2527     }
2528    
2529     UInt_t holeq(Double_t lower,Double_t upper,Quaternions *Qlower, Quaternions *Qupper, UInt_t f){
2530    
2531     UInt_t hole = 10;
2532 mocchiut 1.40 Bool_t R10l = false; // Sign of R10 mode in lower quaternions array
2533     Bool_t R10u = false; // Sign of R10 mode in upper quaternions array
2534     Bool_t insm = false; // Sign that we inside quaternions array
2535 mocchiut 1.64 // Bool_t mxtml = false; // Sign of mixt mode in lower quaternions array
2536     // Bool_t mxtmu = false; // Sign of mixt mode in upper quaternions array
2537 mocchiut 1.40 Bool_t npasm = false; // Sign of normall pass between R10 and non R10 or between non R10 and R10
2538 mocchiut 1.15 UInt_t NCQl = 6; // Number of correct quaternions in lower array
2539 mocchiut 1.64 // UInt_t NCQu = 6; // Number of correct quaternions in upper array
2540 mocchiut 1.15 if (f>0){
2541     insm = true;
2542     if(Qupper->time[f]-Qupper->time[f-1]==30) R10u = false;
2543     if(Qupper->time[f]-Qupper->time[f-1]<1) R10u = true;
2544     }else{
2545     insm = false;
2546     if((Qlower->time[5]-Qlower->time[0]<2)&&(Qlower->time[1]-Qlower->time[0]<2)) R10l = true;
2547     if((Qupper->time[5]-Qupper->time[0]<2)&&(Qupper->time[1]-Qupper->time[0]<2)) R10u = true;
2548     if((Qlower->time[5]-Qlower->time[0]==150)&&(Qlower->time[1]-Qlower->time[0]==30)) R10l = false;
2549     if((Qupper->time[5]-Qupper->time[0]==150)&&(Qupper->time[1]-Qupper->time[0]==30)) R10u = false;
2550     if((Qlower->time[5]-Qlower->time[0]<2)&&(Qlower->time[1]-Qlower->time[0]==30)){
2551 mocchiut 1.64 // mxtml = true;
2552 mocchiut 1.15 for(UInt_t i = 1; i < 6; i++){
2553 mayorov 1.93 if(Qlower->time[i]-Qlower->time[0]==30*i) NCQl=i;
2554 mocchiut 1.15 }
2555     }
2556 mocchiut 1.64 // if((Qupper->time[5]-Qupper->time[0]<2)&&(Qupper->time[1]-Qupper->time[0]==30)){
2557     // mxtmu = true;
2558     // for(UInt_t i = 1; i < 6; i++){
2559 mayorov 1.93 // if(Qupper->time[i]-Qupper->time[0]==30*i) NCQu=i;
2560 mocchiut 1.64 // }
2561     // }
2562 mocchiut 1.15 }
2563    
2564     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;
2565    
2566    
2567     if (R10u&&insm) hole=0; // best event R10
2568     if ((upper-lower<=5)&&(!insm)&&R10l&&R10u) hole = 1; // when first of 6 quaternions in array is correct
2569     if (((!R10u)&&insm)||((!insm)&&(!R10u)&&(!R10l)&&((upper-lower==210+(6-NCQl)*30)||(upper-lower==30)))) hole = 2; //non R10
2570     if (npasm&&(!insm)&&((R10l&&!R10u)||(R10u&&!R10l))) hole = 3; //normall pass from R10 to non R10 or from non R10 to R10
2571     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
2572     if ((upper-lower>=300)&&(!insm)&&((R10l&&!R10u)||(R10u&&!R10l))) hole = 5; //uneliminable hole between R10 and non R10 or between non R10 and R10
2573     if ((upper-lower>5)&&(upper-lower<=300)&&R10u&&R10l) hole = 6; // eliminable hole inside R10
2574     if ((upper-lower>300)&&R10u&&R10l) hole = 7; //uneliminable hole inside R10
2575     if ((upper-lower>210)&&(upper-lower<=1200)&&(!R10u)&&(!R10l)) hole = 8; //eliminable hole inside non R10
2576     if ((upper-lower>1200)&&!R10u&&!R10l) hole = 9; // uneliminable hole inside non R10
2577     return hole;
2578     }
2579    
2580 mocchiut 1.44 void inclresize(vector<Double_t>& t,vector<Float_t>& q0,vector<Float_t>& q1,vector<Float_t>& q2,vector<Float_t>& q3,vector<Int_t>& mode,vector<Float_t>& Roll,vector<Float_t>& Pitch,vector<Float_t>& Yaw){
2581     Int_t sizee = t.size()+1;
2582     t.resize(sizee);
2583     q0.resize(sizee);
2584     q1.resize(sizee);
2585     q2.resize(sizee);
2586     q3.resize(sizee);
2587     mode.resize(sizee);
2588     Roll.resize(sizee);
2589     Pitch.resize(sizee);
2590     Yaw.resize(sizee);
2591     }
2592    
2593 pam-mep 1.69 // geomagnetic calculation staff
2594    
2595     void GM_ScanIGRF(TSQLServer *dbc, GMtype_Data *G0, GMtype_Data *G1, GMtype_Data *H1)
2596     {
2597     GL_PARAM *glp = new GL_PARAM();
2598     Int_t parerror=glp->Query_GL_PARAM(1,304,dbc); // parameters stored in DB in GL_PRAM table
2599     if ( parerror<0 ) {
2600     throw -902;
2601     }
2602 mayorov 1.93 /*This function scans inputs G0, G1, and H1 of the IGRF table into 3 data arrays*/
2603     // TString SATH="/data03/Malakhov/pam9Malakhov/installed10/calib/orb-param/";
2604     int i;
2605     double temp;
2606     char buffer[200];
2607     FILE *IGRF;
2608     IGRF = fopen((glp->PATH+glp->NAME).Data(), "r");
2609     // IGRF = fopen(PATH+"IGRF.tab", "r");
2610     G0->size = 25;
2611     G1->size = 25;
2612     H1->size = 25;
2613     for( i = 0; i < 4; i++)
2614     {
2615     fgets(buffer, 200, IGRF);
2616     }
2617     fscanf(IGRF, "g 1 0 %lf ", &G0->element[0]);
2618     for(i = 1; i <= 22; i++)
2619     {
2620     fscanf(IGRF ,"%lf ", &G0->element[i]);
2621     }
2622     fscanf(IGRF ,"%lf\n", &temp);
2623     G0->element[23] = temp * 5 + G0->element[22];
2624     G0->element[24] = G0->element[23] + 5 * temp;
2625     fscanf(IGRF, "g 1 1 %lf ", &G1->element[0]);
2626     for(i = 1; i <= 22; i++)
2627     {
2628     fscanf( IGRF, "%lf ", &G1->element[i]);
2629     }
2630     fscanf(IGRF, "%lf\n", &temp);
2631     G1->element[23] = temp * 5 + G1->element[22];
2632     G1->element[24] = temp * 5 + G1->element[23];
2633     fscanf(IGRF, "h 1 1 %lf ", &H1->element[0]);
2634     for(i = 1; i <= 22; i++)
2635     {
2636     fscanf( IGRF, "%lf ", &H1->element[i]);
2637     }
2638     fscanf(IGRF, "%lf\n", &temp);
2639     H1->element[23] = temp * 5 + H1->element[22];
2640     H1->element[24] = temp * 5 + H1->element[23];
2641 pam-mep 1.69 if ( glp ) delete glp;
2642 mocchiut 1.81 /*
2643     printf("############################## SCAN IGRF ######################################\n");
2644     printf(" G0 G1 H1\n");
2645     printf(" size %10i %10i %10i \n",G0->size,G1->size,H1->size);
2646     for ( i = 0; i < 30; i++){
2647     printf("%5i %10.2f %10.2f %10.2f \n",i,G0->element[i],G1->element[i],H1->element[i]);
2648     }
2649     printf("###############################################################################\n");
2650     */
2651     } /*GM_ScanIGRF*/
2652    
2653    
2654    
2655    
2656     void GM_SetIGRF(Int_t isSecular, TString ifile1, TString ifile2, GMtype_Data *G0, GMtype_Data *G1, GMtype_Data *H1)
2657     {
2658     /*This function scans inputs G0, G1, and H1 of the IGRF table into 3 data arrays*/
2659     int i;
2660     double temp,temp2;
2661     int it1,it2;
2662     char buffer[200];
2663     FILE *IGRF;
2664     G0->size = 2;
2665     G1->size = 2;
2666     H1->size = 2;
2667    
2668     for( i = 0; i < 30; i++){
2669     G0->element[i] = 0.;
2670     G1->element[i] = 0.;
2671     H1->element[i] = 0.;
2672     }
2673    
2674     IGRF = fopen(ifile1.Data(), "r");
2675     for( i = 0; i < 2; i++){
2676     fgets(buffer, 200, IGRF);
2677     }
2678     fscanf(IGRF, "%3i%3i%12lf%11lf",&it1,&it2, &G0->element[0],&temp);
2679     fscanf(IGRF, "%3i%3i%12lf%11lf",&it1,&it2, &G1->element[0],&H1->element[0]);
2680     fclose(IGRF);
2681    
2682     IGRF = fopen(ifile2.Data(), "r");
2683     for( i = 0; i < 2; i++){
2684     fgets(buffer, 200, IGRF);
2685     }
2686     if ( isSecular ){
2687     fscanf(IGRF, "%3i%3i%12lf%11lf",&it1,&it2,&temp,&temp2);
2688     G0->element[1] = temp * 5. + G0->element[0];
2689     fscanf(IGRF, "%3i%3i%12lf%11lf",&it1,&it2,&temp,&temp2);
2690     G1->element[1] = temp * 5. + G1->element[0];
2691     H1->element[1] = temp2 * 5. + H1->element[0];
2692     } else {
2693     fscanf(IGRF, "%3i%3i%12lf%11lf",&it1,&it2, &G0->element[1],&temp);
2694     fscanf(IGRF, "%3i%3i%12lf%11lf",&it1,&it2, &G1->element[1],&H1->element[1]);
2695     }
2696     fclose(IGRF);
2697     /*
2698     printf("############################## SCAN IGRF ######################################\n");
2699     printf(" G0 G1 H1\n");
2700     printf(" size %10i %10i %10i \n",G0->size,G1->size,H1->size);
2701     for ( i = 0; i < 30; i++){
2702     printf("%5i %10.2f %10.2f %10.2f \n",i,G0->element[i],G1->element[i],H1->element[i]);
2703     }
2704     printf("###############################################################################\n");
2705     */
2706 pam-mep 1.69 } /*GM_ScanIGRF*/
2707    
2708     void GM_SetEllipsoid(GMtype_Ellipsoid *Ellip)
2709     {
2710 mayorov 1.93 /*This function sets the WGS84 reference ellipsoid to its default values*/
2711     Ellip->a = 6378.137; /*semi-major axis of the ellipsoid in */
2712     Ellip->b = 6356.7523142;/*semi-minor axis of the ellipsoid in */
2713     Ellip->fla = 1/298.257223563;/* flattening */
2714     Ellip->eps = sqrt(1- ( Ellip->b * Ellip->b) / (Ellip->a * Ellip->a )); /*first eccentricity */
2715     Ellip->epssq = (Ellip->eps * Ellip->eps); /*first eccentricity squared */
2716     Ellip->re = 6371.2;/* Earth's radius */
2717 pam-mep 1.69 } /*GM_SetEllipsoid*/
2718    
2719    
2720     void GM_EarthCartToDipoleCartCD(GMtype_Pole Pole, GMtype_CoordCartesian EarthCoord, GMtype_CoordCartesian *DipoleCoords)
2721     {
2722 mayorov 1.93 /*This function converts from Earth centered cartesian coordinates to dipole centered cartesian coordinates*/
2723     double X, Y, Z, CosPhi, SinPhi, CosLambda, SinLambda;
2724     CosPhi = cos(TMath::DegToRad()*Pole.phi);
2725     SinPhi = sin(TMath::DegToRad()*Pole.phi);
2726     CosLambda = cos(TMath::DegToRad()*Pole.lambda);
2727     SinLambda = sin(TMath::DegToRad()*Pole.lambda);
2728     X = EarthCoord.x;
2729     Y = EarthCoord.y;
2730     Z = EarthCoord.z;
2731    
2732     /*These equations are taken from a document by Wallace H. Campbell*/
2733     DipoleCoords->x = X * CosPhi * CosLambda + Y * CosPhi * SinLambda - Z * SinPhi;
2734     DipoleCoords->y = -X * SinLambda + Y * CosLambda;
2735     DipoleCoords->z = X * SinPhi * CosLambda + Y * SinPhi * SinLambda + Z * CosPhi;
2736 pam-mep 1.69 } /*GM_EarthCartToDipoleCartCD*/
2737    
2738     void GM_GeodeticToSpherical(GMtype_Ellipsoid Ellip, GMtype_CoordGeodetic CoordGeodetic, GMtype_CoordSpherical *CoordSpherical)
2739     {
2740 mayorov 1.93 double CosLat, SinLat, rc, xp, zp; /*all local variables */
2741     /*
2742     ** Convert geodetic coordinates, (defined by the WGS-84
2743     ** reference ellipsoid), to Earth Centered Earth Fixed Cartesian
2744     ** coordinates, and then to spherical coordinates.
2745     */
2746 pam-mep 1.69
2747 mayorov 1.93 CosLat = cos(TMath::DegToRad()*CoordGeodetic.phi);
2748     SinLat = sin(TMath::DegToRad()*CoordGeodetic.phi);
2749 pam-mep 1.69
2750 mayorov 1.93 /* compute the local radius of curvature on the WGS-84 reference ellipsoid */
2751 pam-mep 1.69
2752 mayorov 1.93 rc = Ellip.a / sqrt(1.0 - Ellip.epssq * SinLat * SinLat);
2753 pam-mep 1.69
2754 mayorov 1.93 /* compute ECEF Cartesian coordinates of specified point (for longitude=0) */
2755 pam-mep 1.69
2756 mayorov 1.93 xp = (rc + CoordGeodetic.HeightAboveEllipsoid) * CosLat;
2757     zp = (rc*(1.0 - Ellip.epssq) + CoordGeodetic.HeightAboveEllipsoid) * SinLat;
2758 pam-mep 1.69
2759 mayorov 1.93 /* compute spherical radius and angle lambda and phi of specified point */
2760 pam-mep 1.69
2761 mayorov 1.93 CoordSpherical->r = sqrt(xp * xp + zp * zp);
2762     CoordSpherical->phig = TMath::RadToDeg()*asin(zp / CoordSpherical->r); /* geocentric latitude */
2763     CoordSpherical->lambda = CoordGeodetic.lambda; /* longitude */
2764 pam-mep 1.69 } /*GM_GeodeticToSpherical*/
2765    
2766     void GM_PoleLocation(GMtype_Model Model, GMtype_Pole *Pole)
2767     {
2768 mayorov 1.93 /*This function finds the location of the north magnetic pole in spherical coordinates. The equations are
2769     **from Wallace H. Campbell's Introduction to Geomagnetic Fields*/
2770 pam-mep 1.69
2771 mayorov 1.93 Pole->phi = TMath::RadToDeg()*-atan(sqrt(Model.h1 * Model.h1 + Model.g1 * Model.g1)/Model.g0);
2772     Pole->lambda = TMath::RadToDeg()*atan(Model.h1/Model.g1);
2773 pam-mep 1.69 } /*GM_PoleLocation*/
2774    
2775     void GM_SphericalToCartesian(GMtype_CoordSpherical CoordSpherical, GMtype_CoordCartesian *CoordCartesian)
2776     {
2777 mayorov 1.93 /*This function converts spherical coordinates into Cartesian coordinates*/
2778     double CosPhi = cos(TMath::DegToRad()*CoordSpherical.phig);
2779     double SinPhi = sin(TMath::DegToRad()*CoordSpherical.phig);
2780     double CosLambda = cos(TMath::DegToRad()*CoordSpherical.lambda);
2781     double SinLambda = sin(TMath::DegToRad()*CoordSpherical.lambda);
2782    
2783     CoordCartesian->x = CoordSpherical.r * CosPhi * CosLambda;
2784     CoordCartesian->y = CoordSpherical.r * CosPhi * SinLambda;
2785     CoordCartesian->z = CoordSpherical.r * SinPhi;
2786 pam-mep 1.69 } /*GM_SphericalToCartesian*/
2787    
2788     void GM_TimeAdjustCoefs(Float_t year, Float_t jyear, GMtype_Data g0d, GMtype_Data g1d, GMtype_Data h1d, GMtype_Model *Model)
2789     {
2790 mayorov 1.93 /*This function calls GM_LinearInterpolation for the coefficients to estimate the value of the
2791     **coefficient for the given date*/
2792     int index;
2793     double x;
2794     index = (year - GM_STARTYEAR) / 5;
2795     x = (jyear - GM_STARTYEAR) / 5.;
2796     Model->g0 = GM_LinearInterpolation(index, index+1, g0d.element[index], g0d.element[index+1], x);
2797     Model->g1 = GM_LinearInterpolation(index, index+1, g1d.element[index], g1d.element[index+1], x);
2798     Model->h1 = GM_LinearInterpolation(index, index+1, h1d.element[index], h1d.element[index+1], x);
2799 pam-mep 1.69 } /*GM_TimeAdjustCoefs*/
2800    
2801     double GM_LinearInterpolation(double x1, double x2, double y1, double y2, double x)
2802     {
2803 mayorov 1.93 /*This function takes a linear interpolation between two given points for x*/
2804     double weight, y;
2805     weight = (x - x1) / (x2 - x1);
2806     y = y1 * (1. - weight) + y2 * weight;
2807 mocchiut 1.81 // printf(" x1 %f x2 %f y1 %f y2 %f x %f ==> y %f \n",x1,x2,y1,y2,x,y);
2808 mayorov 1.93 return y;
2809 pam-mep 1.69 }/*GM_LinearInterpolation*/
2810 mocchiut 1.44
2811 pam-mep 1.69 void GM_CartesianToSpherical(GMtype_CoordCartesian CoordCartesian, GMtype_CoordSpherical *CoordSpherical)
2812     {
2813 mayorov 1.93 /*This function converts a point from Cartesian coordinates into spherical coordinates*/
2814     double X, Y, Z;
2815    
2816     X = CoordCartesian.x;
2817     Y = CoordCartesian.y;
2818     Z = CoordCartesian.z;
2819    
2820     CoordSpherical->r = sqrt(X * X + Y * Y + Z * Z);
2821     CoordSpherical->phig = TMath::RadToDeg()*asin(Z / (CoordSpherical->r));
2822     CoordSpherical->lambda = TMath::RadToDeg()*atan2(Y, X);
2823 pam-mep 1.69 } /*GM_CartesianToSpherical*/
  ViewVC Help
Powered by ViewVC 1.1.23