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