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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.72 by mocchiut, Tue Apr 15 15:48:17 2014 UTC revision 1.73 by mocchiut, Fri Jun 6 14:18:20 2014 UTC
# Line 53  Line 53 
53  //  //
54  #include <ToFLevel2.h>  #include <ToFLevel2.h>
55  #include <TrkLevel2.h>  #include <TrkLevel2.h>
56    #include <ExtTrack.h> // new tracking code
57    
58  using namespace std;  using namespace std;
59    
# Line 279  int OrbitalInfoCore(UInt_t run, TFile *f Line 280  int OrbitalInfoCore(UInt_t run, TFile *f
280    vector<Float_t> recq2;    vector<Float_t> recq2;
281    vector<Float_t> recq3;    vector<Float_t> recq3;
282    Float_t Norm = 1;    Float_t Norm = 1;
283    
284      vector<UInt_t> RTtime1;
285      vector<UInt_t> RTtime2;
286      vector<Double_t> RTbank1;
287      vector<Double_t> RTbank2;
288      vector<Int_t> RTazim;
289      vector<Int_t> RTdir1;
290      vector<Int_t> RTdir2;
291      vector<Int_t> RTerrq;
292    
293      TClonesArray *tcNucleiTrk = NULL;
294      TClonesArray *tcExtNucleiTrk = NULL;
295      TClonesArray *tcExtTrk = NULL;
296      TClonesArray *tcNucleiTof = NULL;
297      TClonesArray *tcExtNucleiTof = NULL;
298      TClonesArray *tcExtTof = NULL;
299      TClonesArray *torbNucleiTrk = NULL;
300      TClonesArray *torbExtNucleiTrk = NULL;
301      TClonesArray *torbExtTrk = NULL;
302      Bool_t hasNucleiTrk = false;
303      Bool_t hasExtNucleiTrk = false;
304      Bool_t hasExtTrk = false;
305      Bool_t hasNucleiTof = false;
306      Bool_t hasExtNucleiTof = false;
307      Bool_t hasExtTof = false;
308    
309      ifstream in;
310      ifstream an;
311      //  ofstream mc;
312      //  TString gr;
313      Int_t parerror2=0;
314    
315    Int_t parerror=glparam0->Query_GL_PARAM(1,303,dbc); // parameters stored in DB in GL_PRAM table    Int_t parerror=glparam0->Query_GL_PARAM(1,303,dbc); // parameters stored in DB in GL_PRAM table
316    if ( verbose ) cout<<parerror<<"\t"<<(char*)(glparam0->PATH+glparam0->NAME).Data()<<endl;    if ( verbose ) cout<<parerror<<"\t"<<(char*)(glparam0->PATH+glparam0->NAME).Data()<<endl;
   ifstream in((char*)(glparam0->PATH+glparam0->NAME).Data(),ios::in);  
317    if ( parerror<0 ) {    if ( parerror<0 ) {
318      code = parerror;      code = parerror;
319      //goto closeandexit;      goto closeandexit;
320    }    }
321      in.open((char*)(glparam0->PATH+glparam0->NAME).Data(),ios::in);
322    while(!in.eof()){    while(!in.eof()){
323      recqtime.resize(recqtime.size()+1);      recqtime.resize(recqtime.size()+1);
324      Int_t sizee = recqtime.size();      Int_t sizee = recqtime.size();
# Line 304  int OrbitalInfoCore(UInt_t run, TFile *f Line 337  int OrbitalInfoCore(UInt_t run, TFile *f
337    if ( verbose ) cout<<"We have read recovered data"<<endl;    if ( verbose ) cout<<"We have read recovered data"<<endl;
338    if (debug) cout << "size of recovered quaterions data set is " << recqtime.size() << endl;    if (debug) cout << "size of recovered quaterions data set is " << recqtime.size() << endl;
339    
340    vector<UInt_t> RTtime1;    // 10RED CHECK
341    vector<UInt_t> RTtime2;    
342    vector<Double_t> RTbank1;    //  TH2F* DIFFX = new TH2F("diffx","",100,0,100,90,0,90);
343    vector<Double_t> RTbank2;    //  TH2F* DIFFY = new TH2F("diffy","",100,0,100,90,0,90);
344    vector<Int_t> RTazim;    //  TH2F* DIFFZ = new TH2F("diffz","",100,0,100,90,0,90);
345    vector<Int_t> RTdir1;    
346    vector<Int_t> RTdir2;    //  gr = "methodscomparison_";
347    vector<Int_t> RTerrq;    //  gr+=run;
348      //  gr+=".txt";
349  // 10RED CHECK    //  mc.open(gr);
350      //  mc.setf(ios::fixed,ios::floatfield);
351  //  TH2F* DIFFX = new TH2F("diffx","",100,0,100,90,0,90);    // 10RED CHECK END
 //  TH2F* DIFFY = new TH2F("diffy","",100,0,100,90,0,90);  
 //  TH2F* DIFFZ = new TH2F("diffz","",100,0,100,90,0,90);  
   
   ofstream mc;  
   TString gr = "methodscomparison_";  
   gr+=run;  
   gr+=".txt";  
   mc.open(gr);  
   mc.setf(ios::fixed,ios::floatfield);  
 // 10RED CHECK END  
352    
353    if ( verbose ) cout<<"read Rotation Table"<<endl;    if ( verbose ) cout<<"read Rotation Table"<<endl;
354      
355    Int_t parerror2=glparam0->Query_GL_PARAM(1,305,dbc);    parerror2=glparam0->Query_GL_PARAM(1,305,dbc);
   ifstream an((char*)(glparam0->PATH+glparam0->NAME).Data(),ios::in);  
356    if ( verbose ) cout<<parerror2<<"\t"<<(char*)(glparam0->PATH+glparam0->NAME).Data()<<endl;    if ( verbose ) cout<<parerror2<<"\t"<<(char*)(glparam0->PATH+glparam0->NAME).Data()<<endl;
357  //  if ( parerror2<0 ) {    if ( parerror2<0 ) {
358  //    code = parerror;      code = parerror;
359      //goto closeandexit;      goto closeandexit;
360  //  }    }
361      an.open((char*)(glparam0->PATH+glparam0->NAME).Data(),ios::in);
362    
 //ifstream an("/data03/Malakhov/pam9Malakhov/installed10/calib/orb-param/RDBCC.txt",ios::in);  
363    while(!an.eof()){    while(!an.eof()){
364      RTtime1.resize(RTtime1.size()+1);      RTtime1.resize(RTtime1.size()+1);
365      Int_t sizee = RTtime1.size();      Int_t sizee = RTtime1.size();
# Line 392  int OrbitalInfoCore(UInt_t run, TFile *f Line 414  int OrbitalInfoCore(UInt_t run, TFile *f
414    //    //
415    if ( !standalone ){    if ( !standalone ){
416      //      //
417      // Does it contain the Tracker tree?      // Does it contain the Tracker and ToF trees?
418      //      //
419      ttof = (TTree*)file->Get("ToF");      ttof = (TTree*)file->Get("ToF");
420      if ( !ttof ) {      if ( !ttof ) {
# Line 403  int OrbitalInfoCore(UInt_t run, TFile *f Line 425  int OrbitalInfoCore(UInt_t run, TFile *f
425      ttof->SetBranchAddress("ToFLevel2",&tof);        ttof->SetBranchAddress("ToFLevel2",&tof);  
426      nevtofl2 = ttof->GetEntries();      nevtofl2 = ttof->GetEntries();
427    
428        //
429        // Look for extended tracking algorithm
430        //
431        if ( verbose ) printf("Look for extended and nuclei tracking algorithms in ToF\n");
432        // Nuclei tracking algorithm
433        Int_t checkAlgo = 0;
434        tcNucleiTof =  new TClonesArray("ToFTrkVar");
435        checkAlgo = ttof->SetBranchAddress("TrackNuclei",&tcNucleiTof);    
436        if ( !checkAlgo ){
437          if ( verbose ) printf(" Nuclei tracking algorithm ToF branch found! :D \n");
438          hasNucleiTof = true;
439        } else {
440          if ( verbose ) printf(" Nuclei tracking algorithm ToF branch not found :( !\n");
441          printf(" ok, this is not a problem (it depends on tracker settings) \n");
442          delete tcNucleiTof;
443        }
444        // Nuclei tracking algorithm using calorimeter points
445        tcExtNucleiTof =  new TClonesArray("ToFTrkVar");
446        checkAlgo = ttof->SetBranchAddress("RecoveredTrackNuclei",&tcExtNucleiTof);    
447        if ( !checkAlgo ){
448          if ( verbose ) printf(" Recovered nuclei tracking algorithm ToF branch found! :D \n");
449          hasExtNucleiTof = true;
450        } else {
451          if ( verbose ) printf(" Recovered nuclei tracking algorithm ToF branch not found :( !\n");
452          printf(" ok, this is not a problem (it depends on tracker settings) \n");
453          delete tcExtNucleiTof;
454        }
455        // Tracking algorithm using calorimeter points
456        tcExtTof =  new TClonesArray("ToFTrkVar");
457        checkAlgo = ttof->SetBranchAddress("RecoveredTrack",&tcExtTof);
458        if ( !checkAlgo ){
459          if ( verbose ) printf(" Recovered track algorithm ToF branch found! :D \n");
460          hasExtTof = true;
461        } else {
462          if ( verbose ) printf(" Recovered track algorithm ToF branch not found :( !\n");
463          printf(" ok, this is not a problem (it depends on tracker settings) \n");
464          delete tcExtTof;
465        }
466    
467      ttrke = (TTree*)file->Get("Tracker");      ttrke = (TTree*)file->Get("Tracker");
468      if ( !ttrke ) {      if ( !ttrke ) {
469        if ( verbose ) printf(" OrbitalInfo - ERROR: no trk tree\n");        if ( verbose ) printf(" OrbitalInfo - ERROR: no trk tree\n");
# Line 411  int OrbitalInfoCore(UInt_t run, TFile *f Line 472  int OrbitalInfoCore(UInt_t run, TFile *f
472      }      }
473      ttrke->SetBranchAddress("TrkLevel2",&trke);        ttrke->SetBranchAddress("TrkLevel2",&trke);  
474      nevtrkl2 = ttrke->GetEntries();      nevtrkl2 = ttrke->GetEntries();
475    
476        //
477        // Look for extended tracking algorithm
478        //
479        if ( verbose ) printf("Look for extended and nuclei tracking algorithms\n");
480        // Nuclei tracking algorithm
481        checkAlgo = 0;
482        tcNucleiTrk =  new TClonesArray("TrkTrack");
483        checkAlgo = ttrke->SetBranchAddress("TrackNuclei",&tcNucleiTrk);    
484        if ( !checkAlgo ){
485          if ( verbose ) printf(" Nuclei tracking algorithm branch found! :D \n");
486          hasNucleiTrk = true;
487        } else {
488          if ( verbose ) printf(" Nuclei tracking algorithm branch not found :( !\n");
489          printf(" ok, this is not a problem (it depends on tracker settings) \n");
490          delete tcNucleiTrk;
491        }
492        // Nuclei tracking algorithm using calorimeter points
493        tcExtNucleiTrk =  new TClonesArray("ExtTrack");
494        checkAlgo = ttrke->SetBranchAddress("RecoveredTrackNuclei",&tcExtNucleiTrk);    
495        if ( !checkAlgo ){
496          if ( verbose ) printf(" Recovered nuclei tracking algorithm branch found! :D \n");
497          hasExtNucleiTrk = true;
498        } else {
499          if ( verbose ) printf(" Recovered nuclei tracking algorithm branch not found :( !\n");
500          printf(" ok, this is not a problem (it depends on tracker settings) \n");
501          delete tcExtNucleiTrk;
502        }
503        // Tracking algorithm using calorimeter points
504        tcExtTrk =  new TClonesArray("ExtTrack");
505        checkAlgo = ttrke->SetBranchAddress("RecoveredTrack",&tcExtTrk);
506        if ( !checkAlgo ){
507          if ( verbose ) printf(" Recovered track algorithm branch found! :D \n");
508          hasExtTrk = true;
509        } else {
510          if ( verbose ) printf(" Recovered track algorithm branch not found :( !\n");
511          printf(" ok, this is not a problem (it depends on tracker settings) \n");
512          delete tcExtTrk;
513        }
514    
515        if ( (hasNucleiTrk && !hasNucleiTof) || (!hasNucleiTrk && hasNucleiTof) ||
516             (hasExtNucleiTrk && !hasExtNucleiTof) || (!hasExtNucleiTrk && hasExtNucleiTof) ||
517             (hasExtTrk && !hasExtTof) || (!hasExtTrk && hasExtTof)
518             ){
519          if ( verbose ) printf(" ERROR: Mismatch between tracker and tof tree branches concerning extended tracking algorithm(s)\n");
520          if ( debug ) printf("hasNucleiTrk %i hasExtNucleiTrk %i hasExtTrk %i \n",hasNucleiTrk,hasExtNucleiTrk,hasExtTrk);
521          if ( debug ) printf("hasNucleiTof %i hasExtNucleiTof %i hasExtTof %i \n",hasNucleiTof,hasExtNucleiTof,hasExtTof);
522          throw -901;
523        }
524    
525    }    }
526    //    //
527    // Let's start!    // Let's start!
# Line 539  int OrbitalInfoCore(UInt_t run, TFile *f Line 650  int OrbitalInfoCore(UInt_t run, TFile *f
650    orbitalinfo->Set();//ELENA **TEMPORANEO?**    orbitalinfo->Set();//ELENA **TEMPORANEO?**
651    OrbitalInfotr->Branch("OrbitalInfo","OrbitalInfo",&orbitalinfo);    OrbitalInfotr->Branch("OrbitalInfo","OrbitalInfo",&orbitalinfo);
652    //    //
653      // create new branches for new tracking algorithms
654      //
655      if ( hasNucleiTrk ){
656        torbNucleiTrk = new TClonesArray("OrbitalInfoTrkVar",1);
657        OrbitalInfotr->Branch("TrackNuclei",&torbNucleiTrk);
658      }
659      if ( hasExtNucleiTrk ){
660        torbExtNucleiTrk = new TClonesArray("OrbitalInfoTrkVar",1);
661        OrbitalInfotr->Branch("RecoveredTrackNuclei",&torbExtNucleiTrk);
662      }
663      if ( hasExtTrk ){
664        torbExtTrk = new TClonesArray("OrbitalInfoTrkVar",1);
665        OrbitalInfotr->Branch("RecoveredTrack",&torbExtTrk);
666      }
667    
668      //
669    if ( reproc && !reprocall ){    if ( reproc && !reprocall ){
670      //      //
671      //  open new file and retrieve also tree informations      //  open new file and retrieve also tree informations
# Line 735  int OrbitalInfoCore(UInt_t run, TFile *f Line 862  int OrbitalInfoCore(UInt_t run, TFile *f
862      // End IGRF stuff//      // End IGRF stuff//
863      //      //
864    
     //  
     //     TTree *tp = (TTree*)l0File->Get("RunHeader");  
     //     tp->SetBranchAddress("Header", &eH);  
     //     tp->SetBranchAddress("RunHeader", &reh);  
     //     tp->GetEntry(0);  
     //     ph = eH->GetPscuHeader();  
     //     ULong_t TimeSync = reh->LAST_TIME_SYNC_INFO;  
     //     ULong_t ObtSync = reh->OBT_TIME_SYNC;      
     //     if ( debug ) printf(" 1 TimeSync %lu ObtSync %lu DeltaOBT %lu\n",TimeSync,ObtSync,TimeSync-ObtSync);  
     //  
865      ULong_t TimeSync = (ULong_t)dbtime->GetTimesync();      ULong_t TimeSync = (ULong_t)dbtime->GetTimesync();
866      ULong_t ObtSync = (ULong_t)(dbtime->GetObt0()/1000);      ULong_t ObtSync = (ULong_t)(dbtime->GetObt0()/1000);
867      ULong_t DeltaOBT = TimeSync - ObtSync;      ULong_t DeltaOBT = TimeSync - ObtSync;
# Line 831  int OrbitalInfoCore(UInt_t run, TFile *f Line 948  int OrbitalInfoCore(UInt_t run, TFile *f
948        i++;        i++;
949      };      };
950      //      //
     //    l0trm = (TTree*)l0File->Get("Mcmd");  
     //    ch->ls();  
951      ch->SetBranchAddress("Mcmd",&mcmdev);      ch->SetBranchAddress("Mcmd",&mcmdev);
     //    printf(" entries %llu \n", ch->GetEntries());  
     //    l0trm = ch->GetTree();  
     //    neventsm = l0trm->GetEntries();  
952      neventsm = ch->GetEntries();      neventsm = ch->GetEntries();
953      if ( debug ) printf(" entries %u \n", neventsm);      if ( debug ) printf(" entries %u \n", neventsm);
     //    neventsm = 0;  
954      //      //
955      if (neventsm == 0){      if (neventsm == 0){
956        if ( debug ) printf("InclinationInfo - WARNING: No quaternions in this File");        if ( debug ) printf("InclinationInfo - WARNING: No quaternions in this File");
       //      l0File->Close();  
957        code = 900;        code = 900;
       //      goto closeandexit;  
958      }      }
959      //      //
960            Double_t lowerqtime = 0;    
     //    l0trm->SetBranchAddress("Mcmd", &mcmdev);  
     //    l0trm->SetBranchAddress("Header", &eh);  
     //  
     //  
     //  
   
 //    UInt_t mctren = 0;      
 //    UInt_t mcreen = 0;          
 //    UInt_t numrec = 0;  
     //  
     //    Double_t upperqtime = 0;  
     Double_t lowerqtime = 0;  
       
     //    Double_t incli = 0;  
     //    oi = 0;  
     //    UInt_t ooi = 0;  
961      //      //
962      // init quaternions information from mcmd-packets      // init quaternions information from mcmd-packets
963      //      //
964      Bool_t isf = true;      Bool_t isf = true;
     //    Int_t fgh = 0;  
965    
966      vector<Float_t> q0;      vector<Float_t> q0;
967      vector<Float_t> q1;      vector<Float_t> q1;
# Line 882  int OrbitalInfoCore(UInt_t run, TFile *f Line 974  int OrbitalInfoCore(UInt_t run, TFile *f
974      vector<Int_t> qmode;      vector<Int_t> qmode;
975    
976      Int_t nt = 0;      Int_t nt = 0;
       
977      UInt_t must = 0;      UInt_t must = 0;
978    
979      //      //
# Line 934  int OrbitalInfoCore(UInt_t run, TFile *f Line 1025  int OrbitalInfoCore(UInt_t run, TFile *f
1025          //          //
1026          tof->Clear();          tof->Clear();
1027          //          //
1028            // Clones array must be cleared before going on
1029            //
1030            if ( hasNucleiTof ){
1031              tcNucleiTof->Delete();
1032            }
1033            if ( hasExtNucleiTof ){
1034              tcExtNucleiTof->Delete();
1035            }          
1036            if ( hasExtTof ){
1037              tcExtTof->Delete();
1038            }
1039            //
1040            if ( verbose ) printf(" get tof tree entries... entry = %i in Level2 file\n",itr);
1041          if ( ttof->GetEntry(itr) <= 0 ){          if ( ttof->GetEntry(itr) <= 0 ){
1042           if ( verbose ) printf(" problems with tof tree entries... entry = %i in Level2 file\n",itr);            if ( verbose ) printf(" problems with tof tree entries... entry = %i in Level2 file\n",itr);
1043           if ( verbose ) printf(" nobefrun %u re %u evfrom %u jumped %u reprocall %i \n",nobefrun,re,evfrom,jumped,reprocall);            if ( verbose ) printf(" nobefrun %u re %u evfrom %u jumped %u reprocall %i \n",nobefrun,re,evfrom,jumped,reprocall);
1044           throw -36;            throw -36;
1045          }          }
1046            if ( verbose ) printf(" gat0\n");
1047          //          //
1048        }        }
1049        //        //
# Line 951  int OrbitalInfoCore(UInt_t run, TFile *f Line 1056  int OrbitalInfoCore(UInt_t run, TFile *f
1056            l0File->Close();            l0File->Close();
1057            code = -905;            code = -905;
1058            goto closeandexit;            goto closeandexit;
1059          };          }
1060          //          //
1061            if ( verbose ) printf(" gat1\n");
1062          trke->Clear();          trke->Clear();
1063          //          //
1064            // Clones array must be cleared before going on
1065            //
1066            if ( hasNucleiTrk ){
1067              if ( verbose ) printf(" gat2\n");
1068              tcNucleiTrk->Delete();
1069              if ( verbose ) printf(" gat3\n");
1070              torbNucleiTrk->Delete();
1071            }
1072            if ( hasExtNucleiTrk ){
1073              if ( verbose ) printf(" gat4\n");
1074              tcExtNucleiTrk->Delete();
1075              if ( verbose ) printf(" gat5\n");
1076              torbExtNucleiTrk->Delete();
1077            }          
1078            if ( hasExtTrk ){
1079              if ( verbose ) printf(" gat6\n");
1080              tcExtTrk->Delete();
1081              if ( verbose ) printf(" gat7\n");
1082              torbExtTrk->Delete();
1083            }
1084            //
1085            if ( verbose ) printf(" get trk tree entries... entry = %i in Level2 file\n",itr);
1086          if ( ttrke->GetEntry(itr) <= 0 ) throw -36;          if ( ttrke->GetEntry(itr) <= 0 ) throw -36;
1087          //          //
1088        }        }
1089    
   
1090        //        //
1091        procev++;        procev++;
1092        //        //
# Line 1455  int OrbitalInfoCore(UInt_t run, TFile *f Line 1582  int OrbitalInfoCore(UInt_t run, TFile *f
1582              if(atime<=RTtime2[mu] && atime>=RTtime1[mu]){              if(atime<=RTtime2[mu] && atime>=RTtime1[mu]){
1583                  // search for angle betwean velosity and direction to north in tangential to Earth surfase plane in satellite position                  // search for angle betwean velosity and direction to north in tangential to Earth surfase plane in satellite position
1584                Double_t tlat=orbitalinfo->lat;                Double_t tlat=orbitalinfo->lat;
1585  /*            Double_t phint=(163.7-0.0002387*tlat-0.005802*tlat*tlat-0.005802e-7*tlat*tlat*tlat-1.776e-6*tlat*tlat*tlat*tlat+1.395e-10*tlat*tlat*tlat*tlat*tlat);                Double_t kar=(RTbank2[mu]-RTbank1[mu])/(RTtime2[mu]-RTtime1[mu]);
               Double_t phin=TMath::Abs(90.0*(1+diro)-phint);  
               Double_t phi=TMath::Abs(90.0*(1-diro)-TMath::RadToDeg()*atan(TMath::Abs(tan(TMath::DegToRad()*phin))/sqrt(1+pow(tan(TMath::DegToRad()*tlat),2))));  
                     
                 //Get vectors of Satellite reference frame axis in GEO in satndard case (No rotations, all Euler angles equals to 0)  
               TVector3 XDij(0,sin(TMath::DegToRad()*phi),cos(TMath::DegToRad()*phi));  
               TVector3 YDij(1,0,0);  
               TVector3 ZDij(0,sin(TMath::DegToRad()*(phi+90)),cos(TMath::DegToRad()*(phi+90.0)));  
   
                 //Get Vectors to rotate about  
               TVector3 B1 = V;  
               B1.RotateZ(-lon*TMath::DegToRad());  
               B1.RotateY(lat*TMath::DegToRad());  
               Float_t elipangle=TMath::ACos((pow(B1.Y(),2)+pow(B1.Z(),2))/B1.Mag()/sqrt(pow(B1.Y(),2)+pow(B1.Z(),2)));  
               TVector3 Tre(0,B1.Y(),B1.Z());  
               if(B1.X()<0) elipangle=-elipangle;  
               TVector3 Vperp=B1;                // axis to rotate around initial Dij on ellip and spitch angles  
               Vperp.RotateX(TMath::Pi()/2.);  
               Vperp.SetX(0);  
 */            Double_t kar=(RTbank2[mu]-RTbank1[mu])/(RTtime2[mu]-RTtime1[mu]);  
1586                Double_t bak=RTbank1[mu]-kar*RTtime1[mu];                Double_t bak=RTbank1[mu]-kar*RTtime1[mu];
1587                Double_t bank=kar*atime+bak;                Double_t bank=kar*atime+bak;
1588                Float_t spitch = 0.00001;  // temprary not zero to avoid problem with tranzition from Euler angles to orientation matrix                Float_t spitch = 0.00001;  // temprary not zero to avoid problem with tranzition from Euler angles to orientation matrix
# Line 1487  int OrbitalInfoCore(UInt_t run, TFile *f Line 1595  int OrbitalInfoCore(UInt_t run, TFile *f
1595                   Float_t bva=spitch1-kva*RTtime1[mu];                   Float_t bva=spitch1-kva*RTtime1[mu];
1596                   spitch=kva*atime+bva;                   spitch=kva*atime+bva;
1597                }                }
1598  /*            //spitch=0.0;  
               //Rotations future Dij matrix on ellip and spitch angles  
               XDij.Rotate(-elipangle-spitch,Vperp);  
               YDij.Rotate(-elipangle-spitch,Vperp);  
               ZDij.Rotate(-elipangle-spitch,Vperp);  
                   
                 //Rotation on bank angle;  
               if(TMath::Abs(bank)>0.5){  
                  XDij.Rotate(TMath::DegToRad()*bank,B1);  
                  YDij.Rotate(TMath::DegToRad()*bank,B1);  
                  ZDij.Rotate(TMath::DegToRad()*bank,B1);  
               }  
               Dij(0,0)=XDij.X(); Dij(1,0)=XDij.Y(); Dij(2,0)=XDij.Z();  
               Dij(0,1)=YDij.X(); Dij(1,1)=YDij.Y(); Dij(2,1)=YDij.Z();  
               Dij(0,2)=ZDij.X(); Dij(1,2)=ZDij.Y(); Dij(2,2)=ZDij.Z();  
 */  
1599                //Calculate Yaw angle accordingly with fit, see picture FitYaw.jpg                //Calculate Yaw angle accordingly with fit, see picture FitYaw.jpg
1600                Double_t yaw=0.00001;  // temprary not zero to avoid problem with tranzition from Euler angles to orientation matrix                Double_t yaw=0.00001;  // temprary not zero to avoid problem with tranzition from Euler angles to orientation matrix
1601                if(TMath::Abs(tlat)<70)                if(TMath::Abs(tlat)<70)
# Line 1553  int OrbitalInfoCore(UInt_t run, TFile *f Line 1646  int OrbitalInfoCore(UInt_t run, TFile *f
1646          //      orbitalinfo->pamBangle = (Float_t)PO->GetPitchAngle(A1,A2,A3,Bx,By,Bz);                 // Angle between Pamela's main axiz and B          //      orbitalinfo->pamBangle = (Float_t)PO->GetPitchAngle(A1,A2,A3,Bx,By,Bz);                 // Angle between Pamela's main axiz and B
1647          //          //
1648          if ( debug ) printf(" matrixes done  \n");          if ( debug ) printf(" matrixes done  \n");
1649          if ( !standalone && tof->ntrk() > 0 ){          if ( !standalone ){
1650            if ( debug ) printf(" !standalone \n");            if ( debug ) printf(" !standalone \n");
1651            //            //
1652              // Standard tracking algorithm
1653              //
1654            Int_t nn = 0;            Int_t nn = 0;
1655              if ( verbose ) printf(" standard tracking \n");
1656            for(Int_t nt=0; nt < tof->ntrk(); nt++){              for(Int_t nt=0; nt < tof->ntrk(); nt++){  
1657              //              //
1658              ToFTrkVar *ptt = tof->GetToFTrkVar(nt);              ToFTrkVar *ptt = tof->GetToFTrkVar(nt);
# Line 1571  int OrbitalInfoCore(UInt_t run, TFile *f Line 1667  int OrbitalInfoCore(UInt_t run, TFile *f
1667                TrkTrack *mytrack = trke->GetStoredTrack(ptt->trkseqno);                TrkTrack *mytrack = trke->GetStoredTrack(ptt->trkseqno);
1668                Float_t rig=1/mytrack->GetDeflection();                Float_t rig=1/mytrack->GetDeflection();
1669                Double_t norm = sqrt(pow(E22x-E11x,2)+pow(E22y-E11y,2)+pow(E22z-E11z,2));                Double_t norm = sqrt(pow(E22x-E11x,2)+pow(E22y-E11y,2)+pow(E22z-E11z,2));
1670                //              Double_t MyAzim = TMath::RadToDeg()*atan(TMath::Abs(E22y-E11y)/TMath::Abs(E22x-E11x));                //
               //              if(E22x-E11x>=0 && E22y-E11y <0) MyAzim =  360. - MyAzim;  
               //              if(E22x-E11x>=0 && E22y-E11y >=0) MyAzim = MyAzim;  
               //              if(E22x-E11x<0 && E22y-E11y >0) MyAzim = 180. - MyAzim;  
               //              if(E22x-E11x<0 && E22y-E11y <0) MyAzim = 180. + MyAzim;  
1671                Px = (E22x-E11x)/norm;                Px = (E22x-E11x)/norm;
1672                Py = (E22y-E11y)/norm;                Py = (E22y-E11y)/norm;
1673                Pz = (E22z-E11z)/norm;                Pz = (E22z-E11z)/norm;
# Line 1592  int OrbitalInfoCore(UInt_t run, TFile *f Line 1684  int OrbitalInfoCore(UInt_t run, TFile *f
1684                //                            //            
1685                t_orb->pitch = (Float_t)PO->GetPitchAngle(Eij(0,0),Eij(1,0),Eij(2,0),Bx,By,Bz);                t_orb->pitch = (Float_t)PO->GetPitchAngle(Eij(0,0),Eij(1,0),Eij(2,0),Bx,By,Bz);
1686                //                //
1687                  // SunPosition in instrumental reference frame                // SunPosition in instrumental reference frame
1688                TMatrixD Kij = PO->PamelatoGEO(qij,Px,Py,Pz);                TMatrixD Kij = PO->PamelatoGEO(qij,Px,Py,Pz);
1689                TMatrixD Lij = PO->PamelatoGEO(qij,0,0,1);                TMatrixD Lij = PO->PamelatoGEO(qij,0,0,1);
1690                t_orb->sunangle=(Float_t)PO->GetPitchAngle(Kij(0,0),Kij(1,0),Kij(2,0),-SP.X(),-SP.Y(),-SP.Z());                t_orb->sunangle=(Float_t)PO->GetPitchAngle(Kij(0,0),Kij(1,0),Kij(2,0),-SP.X(),-SP.Y(),-SP.Z());
1691                t_orb->sunmagangle=(Float_t)PO->GetPitchAngle(Kij(0,0),Kij(1,0),Kij(2,0),SunMag.X(),SunMag.Y(),SunMag.Z());                t_orb->sunmagangle=(Float_t)PO->GetPitchAngle(Kij(0,0),Kij(1,0),Kij(2,0),SunMag.X(),SunMag.Y(),SunMag.Z());
   
1692                //                //
1693                //                //
               //Double_t omega = PO->GetPitchAngle(Eij(0,0),Eij(1,0),Eij(2,0),cos(orbitalinfo->lon+TMath::Pi()/2)-sin(orbitalinfo->lon+TMath::Pi()/2),cos(orbitalinfo->lon+TMath::Pi()/2)+sin(orbitalinfo->lon+TMath::Pi()/2),1);  
1694                Double_t omega = PO->GetPitchAngle(-Eij(0,0),-Eij(1,0),-Eij(2,0),1,0,0) * TMath::DegToRad();                Double_t omega = PO->GetPitchAngle(-Eij(0,0),-Eij(1,0),-Eij(2,0),1,0,0) * TMath::DegToRad();
1695                TVector3 Bxy(0,By,Bz);                TVector3 Bxy(0,By,Bz);
1696                TVector3 Exy(0,-Eij(1,0),-Eij(2,0));                TVector3 Exy(0,-Eij(1,0),-Eij(2,0));
1697                Double_t dzeta=Bxy.Angle(Exy);                Double_t dzeta=Bxy.Angle(Exy);
1698                if (-Eij(1,0) < 0) dzeta=2.0*TMath::Pi() - dzeta;                if (-Eij(1,0) < 0) dzeta=2.0*TMath::Pi() - dzeta;
1699                  
1700                  if(debug) cout << "omega = "<<omega*TMath::RadToDeg()<<"\tdzeta = "<<dzeta*TMath::RadToDeg()<<endl;
1701    
1702                  if(debug) cout << "omega = "<<omega*TMath::RadToDeg()<<"\tdzeta = "<<dzeta*TMath::RadToDeg()<<endl;                // Formula from D.F. Smart *, M.A. Shea [2005]; A review of geomagnetic cutoff rigidities for earth-orbiting spacecraft
   
                 // Formula from D.F. Smart *, M.A. Shea [2005]; A review of geomagnetic cutoff rigidities for earth-orbiting spacecraft  
1703                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));                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));
1704                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));                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));
1705                if (debug) cout << "R = " << rig << "\tcutoff = " << t_orb->cutoff << endl;                if (debug) cout << "R = " << rig << "\tcutoff = " << t_orb->cutoff << endl;
1706    
               //t_orb->cutoff = 59.3/(pow(orbitalinfo->L,2)*pow((1+sqrt(1-pow(orbitalinfo->L,-3/2)*cos(omega))),2));  
   
1707                //                //
1708                if ( t_orb->pitch != t_orb->pitch ) t_orb->pitch = -1000.;                if ( t_orb->pitch != t_orb->pitch ) t_orb->pitch = -1000.;
1709                if ( t_orb->cutoff != t_orb->cutoff ) t_orb->cutoff = -1000.;                if ( t_orb->cutoff != t_orb->cutoff ) t_orb->cutoff = -1000.;
# Line 1629  int OrbitalInfoCore(UInt_t run, TFile *f Line 1717  int OrbitalInfoCore(UInt_t run, TFile *f
1717                //                //
1718                t_orb->Clear();                t_orb->Clear();
1719                //                //
1720              };              }
1721              //              //
1722            };            } // end standard tracking algorithm
1723    
1724              //
1725              // Code for extended tracking algorithm:
1726              //
1727              if ( hasNucleiTrk ){
1728                Int_t ttentry = 0;
1729                if ( verbose ) printf(" hasNucleiTrk \n");
1730                for(Int_t nt=0; nt < tcNucleiTof->GetEntries() ; nt++){  
1731                  //
1732                  if ( verbose ) printf(" got1\n");
1733                  ToFTrkVar *ptt = (ToFTrkVar*)(tcNucleiTof->At(nt));
1734                  if (verbose) cout<<" tcNucleiTof->GetEntries() = "<<tcNucleiTof->GetEntries()<<"\tptt->trkseqno = "<<ptt->trkseqno<<endl;
1735                  Double_t E11x = ptt->xtr_tof[0]; // tr->x[0];
1736                  Double_t E11y = ptt->ytr_tof[0]; //tr->y[0];
1737                  Double_t E11z = zin[0];
1738                  Double_t E22x = ptt->xtr_tof[3];//tr->x[3];
1739                  Double_t E22y = ptt->ytr_tof[3];//tr->y[3];
1740                  Double_t E22z = zin[3];
1741                  if ( (E11x < 100. && E11y < 100. && E22x < 100. && E22y < 100.) || ptt->trkseqno != -1  ){
1742                    TrkTrack *mytrack = (TrkTrack*)(tcNucleiTrk->At(ptt->trkseqno));
1743                    if ( verbose ) printf(" got tcNucleiTrk \n");
1744                    Float_t rig=1/mytrack->GetDeflection();
1745                    Double_t norm = sqrt(pow(E22x-E11x,2)+pow(E22y-E11y,2)+pow(E22z-E11z,2));
1746                    //
1747                    Px = (E22x-E11x)/norm;
1748                    Py = (E22y-E11y)/norm;
1749                    Pz = (E22z-E11z)/norm;
1750                    //
1751                    t_orb->trkseqno = ptt->trkseqno;
1752                    //
1753                    TMatrixD Eij = PO->PamelatoGEO(Iij,Px,Py,Pz);
1754                    t_orb->Eij.ResizeTo(Eij);      
1755                    t_orb->Eij = Eij;      
1756                    //
1757                    TMatrixD Sij = PO->PamelatoGEO(Gij,Px,Py,Pz);
1758                    t_orb->Sij.ResizeTo(Sij);
1759                    t_orb->Sij = Sij;
1760                    //          
1761                    t_orb->pitch = (Float_t)PO->GetPitchAngle(Eij(0,0),Eij(1,0),Eij(2,0),Bx,By,Bz);
1762                    //
1763                    // SunPosition in instrumental reference frame
1764                    TMatrixD Kij = PO->PamelatoGEO(qij,Px,Py,Pz);
1765                    TMatrixD Lij = PO->PamelatoGEO(qij,0,0,1);
1766                    t_orb->sunangle=(Float_t)PO->GetPitchAngle(Kij(0,0),Kij(1,0),Kij(2,0),-SP.X(),-SP.Y(),-SP.Z());
1767                    t_orb->sunmagangle=(Float_t)PO->GetPitchAngle(Kij(0,0),Kij(1,0),Kij(2,0),SunMag.X(),SunMag.Y(),SunMag.Z());
1768                    //
1769                    //
1770                    Double_t omega = PO->GetPitchAngle(-Eij(0,0),-Eij(1,0),-Eij(2,0),1,0,0) * TMath::DegToRad();
1771                    TVector3 Bxy(0,By,Bz);
1772                    TVector3 Exy(0,-Eij(1,0),-Eij(2,0));
1773                    Double_t dzeta=Bxy.Angle(Exy);
1774                    if (-Eij(1,0) < 0) dzeta=2.0*TMath::Pi() - dzeta;
1775                    
1776                    if(debug) cout << "omega = "<<omega*TMath::RadToDeg()<<"\tdzeta = "<<dzeta*TMath::RadToDeg()<<endl;
1777                    
1778                    // Formula from D.F. Smart *, M.A. Shea [2005]; A review of geomagnetic cutoff rigidities for earth-orbiting spacecraft
1779                    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));
1780                    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));
1781                    if (debug) cout << "R = " << rig << "\tcutoff = " << t_orb->cutoff << endl;
1782                    
1783                    //
1784                    if ( t_orb->pitch != t_orb->pitch ) t_orb->pitch = -1000.;
1785                    if ( t_orb->cutoff != t_orb->cutoff ) t_orb->cutoff = -1000.;
1786                    if ( t_orb->sunangle != t_orb->sunangle ) t_orb->sunangle = -1000.;
1787                    if ( t_orb->sunmagangle != t_orb->sunmagangle ) t_orb->sunmagangle = -1000.;
1788                    //
1789                    if ( debug ) printf(" orbitalinfo->cutoffsvl %f vitaly %f \n",orbitalinfo->cutoffsvl,t_orb->cutoff);
1790                    //
1791                    TClonesArray &tt1 = *torbNucleiTrk;
1792                    new(tt1[ttentry]) OrbitalInfoTrkVar(*t_orb);
1793                    ttentry++;
1794                    //
1795                    t_orb->Clear();
1796                    //
1797                  }
1798                  //
1799                }
1800              } // end standard tracking algorithm: nuclei
1801              if ( hasExtNucleiTrk ){
1802                Int_t ttentry = 0;
1803                if ( verbose ) printf(" hasExtNucleiTrk \n");
1804                for(Int_t nt=0; nt < tcExtNucleiTof->GetEntries() ; nt++){  
1805                  //
1806                  if ( verbose ) printf(" got2\n");
1807                  ToFTrkVar *ptt = (ToFTrkVar*)(tcExtNucleiTof->At(nt));
1808                  if (verbose) cout<<" tcExtNucleiTof->GetEntries() = "<<tcExtNucleiTof->GetEntries()<<"\tptt->trkseqno = "<<ptt->trkseqno<<endl;
1809                  Double_t E11x = ptt->xtr_tof[0]; // tr->x[0];
1810                  Double_t E11y = ptt->ytr_tof[0]; //tr->y[0];
1811                  Double_t E11z = zin[0];
1812                  Double_t E22x = ptt->xtr_tof[3];//tr->x[3];
1813                  Double_t E22y = ptt->ytr_tof[3];//tr->y[3];
1814                  Double_t E22z = zin[3];
1815                  if ( (E11x < 100. && E11y < 100. && E22x < 100. && E22y < 100.) || ptt->trkseqno != -1  ){
1816                    ExtTrack *mytrack = (ExtTrack*)(tcExtNucleiTrk->At(ptt->trkseqno));
1817                    if ( verbose ) printf(" got tcExtNucleiTrk \n");
1818                    Float_t rig=1/mytrack->GetDeflection();
1819                    Double_t norm = sqrt(pow(E22x-E11x,2)+pow(E22y-E11y,2)+pow(E22z-E11z,2));
1820                    //
1821                    Px = (E22x-E11x)/norm;
1822                    Py = (E22y-E11y)/norm;
1823                    Pz = (E22z-E11z)/norm;
1824                    //
1825                    t_orb->trkseqno = ptt->trkseqno;
1826                    //
1827                    TMatrixD Eij = PO->PamelatoGEO(Iij,Px,Py,Pz);
1828                    t_orb->Eij.ResizeTo(Eij);      
1829                    t_orb->Eij = Eij;      
1830                    //
1831                    TMatrixD Sij = PO->PamelatoGEO(Gij,Px,Py,Pz);
1832                    t_orb->Sij.ResizeTo(Sij);
1833                    t_orb->Sij = Sij;
1834                    //          
1835                    t_orb->pitch = (Float_t)PO->GetPitchAngle(Eij(0,0),Eij(1,0),Eij(2,0),Bx,By,Bz);
1836                    //
1837                    // SunPosition in instrumental reference frame
1838                    TMatrixD Kij = PO->PamelatoGEO(qij,Px,Py,Pz);
1839                    TMatrixD Lij = PO->PamelatoGEO(qij,0,0,1);
1840                    t_orb->sunangle=(Float_t)PO->GetPitchAngle(Kij(0,0),Kij(1,0),Kij(2,0),-SP.X(),-SP.Y(),-SP.Z());
1841                    t_orb->sunmagangle=(Float_t)PO->GetPitchAngle(Kij(0,0),Kij(1,0),Kij(2,0),SunMag.X(),SunMag.Y(),SunMag.Z());
1842                    //
1843                    //
1844                    Double_t omega = PO->GetPitchAngle(-Eij(0,0),-Eij(1,0),-Eij(2,0),1,0,0) * TMath::DegToRad();
1845                    TVector3 Bxy(0,By,Bz);
1846                    TVector3 Exy(0,-Eij(1,0),-Eij(2,0));
1847                    Double_t dzeta=Bxy.Angle(Exy);
1848                    if (-Eij(1,0) < 0) dzeta=2.0*TMath::Pi() - dzeta;
1849                    
1850                    if(debug) cout << "omega = "<<omega*TMath::RadToDeg()<<"\tdzeta = "<<dzeta*TMath::RadToDeg()<<endl;
1851                    
1852                    // Formula from D.F. Smart *, M.A. Shea [2005]; A review of geomagnetic cutoff rigidities for earth-orbiting spacecraft
1853                    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));
1854                    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));
1855                    if (debug) cout << "R = " << rig << "\tcutoff = " << t_orb->cutoff << endl;
1856                    
1857                    //
1858                    if ( t_orb->pitch != t_orb->pitch ) t_orb->pitch = -1000.;
1859                    if ( t_orb->cutoff != t_orb->cutoff ) t_orb->cutoff = -1000.;
1860                    if ( t_orb->sunangle != t_orb->sunangle ) t_orb->sunangle = -1000.;
1861                    if ( t_orb->sunmagangle != t_orb->sunmagangle ) t_orb->sunmagangle = -1000.;
1862                    //
1863                    if ( debug ) printf(" orbitalinfo->cutoffsvl %f vitaly %f \n",orbitalinfo->cutoffsvl,t_orb->cutoff);
1864                    //
1865                    TClonesArray &tt2 = *torbExtNucleiTrk;
1866                    new(tt2[ttentry]) OrbitalInfoTrkVar(*t_orb);
1867                    ttentry++;
1868                    //
1869                    t_orb->Clear();
1870                    //
1871                  }
1872                  //
1873                }
1874              } // end standard tracking algorithm: nuclei extended
1875             if ( hasExtTrk ){
1876                Int_t ttentry = 0;
1877                if ( verbose ) printf(" hasExtTrk \n");
1878                for(Int_t nt=0; nt < tcExtTof->GetEntries() ; nt++){  
1879                  //
1880                  if ( verbose ) printf(" got3\n");
1881                  ToFTrkVar *ptt = (ToFTrkVar*)(tcExtTof->At(nt));
1882                  if (verbose) cout<<" tcExtTof->GetEntries() = "<<tcExtTof->GetEntries()<<"\tptt->trkseqno = "<<ptt->trkseqno<<endl;
1883                  Double_t E11x = ptt->xtr_tof[0]; // tr->x[0];
1884                  Double_t E11y = ptt->ytr_tof[0]; //tr->y[0];
1885                  Double_t E11z = zin[0];
1886                  Double_t E22x = ptt->xtr_tof[3];//tr->x[3];
1887                  Double_t E22y = ptt->ytr_tof[3];//tr->y[3];
1888                  Double_t E22z = zin[3];
1889                  if ( (E11x < 100. && E11y < 100. && E22x < 100. && E22y < 100.) || ptt->trkseqno != -1  ){
1890                    ExtTrack *mytrack = (ExtTrack*)(tcExtTrk->At(ptt->trkseqno));
1891                    if ( verbose ) printf(" got tcExtTrk \n");
1892                    Float_t rig=1/mytrack->GetDeflection();
1893                    Double_t norm = sqrt(pow(E22x-E11x,2)+pow(E22y-E11y,2)+pow(E22z-E11z,2));
1894                    //
1895                    Px = (E22x-E11x)/norm;
1896                    Py = (E22y-E11y)/norm;
1897                    Pz = (E22z-E11z)/norm;
1898                    //
1899                    t_orb->trkseqno = ptt->trkseqno;
1900                    //
1901                    TMatrixD Eij = PO->PamelatoGEO(Iij,Px,Py,Pz);
1902                    t_orb->Eij.ResizeTo(Eij);      
1903                    t_orb->Eij = Eij;      
1904                    //
1905                    TMatrixD Sij = PO->PamelatoGEO(Gij,Px,Py,Pz);
1906                    t_orb->Sij.ResizeTo(Sij);
1907                    t_orb->Sij = Sij;
1908                    //          
1909                    t_orb->pitch = (Float_t)PO->GetPitchAngle(Eij(0,0),Eij(1,0),Eij(2,0),Bx,By,Bz);
1910                    //
1911                    // SunPosition in instrumental reference frame
1912                    TMatrixD Kij = PO->PamelatoGEO(qij,Px,Py,Pz);
1913                    TMatrixD Lij = PO->PamelatoGEO(qij,0,0,1);
1914                    t_orb->sunangle=(Float_t)PO->GetPitchAngle(Kij(0,0),Kij(1,0),Kij(2,0),-SP.X(),-SP.Y(),-SP.Z());
1915                    t_orb->sunmagangle=(Float_t)PO->GetPitchAngle(Kij(0,0),Kij(1,0),Kij(2,0),SunMag.X(),SunMag.Y(),SunMag.Z());
1916                    //
1917                    //
1918                    Double_t omega = PO->GetPitchAngle(-Eij(0,0),-Eij(1,0),-Eij(2,0),1,0,0) * TMath::DegToRad();
1919                    TVector3 Bxy(0,By,Bz);
1920                    TVector3 Exy(0,-Eij(1,0),-Eij(2,0));
1921                    Double_t dzeta=Bxy.Angle(Exy);
1922                    if (-Eij(1,0) < 0) dzeta=2.0*TMath::Pi() - dzeta;
1923                    
1924                    if(debug) cout << "omega = "<<omega*TMath::RadToDeg()<<"\tdzeta = "<<dzeta*TMath::RadToDeg()<<endl;
1925                    
1926                    // Formula from D.F. Smart *, M.A. Shea [2005]; A review of geomagnetic cutoff rigidities for earth-orbiting spacecraft
1927                    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));
1928                    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));
1929                    if (debug) cout << "R = " << rig << "\tcutoff = " << t_orb->cutoff << endl;
1930                    
1931                    //
1932                    if ( t_orb->pitch != t_orb->pitch ) t_orb->pitch = -1000.;
1933                    if ( t_orb->cutoff != t_orb->cutoff ) t_orb->cutoff = -1000.;
1934                    if ( t_orb->sunangle != t_orb->sunangle ) t_orb->sunangle = -1000.;
1935                    if ( t_orb->sunmagangle != t_orb->sunmagangle ) t_orb->sunmagangle = -1000.;
1936                    //
1937                    if ( debug ) printf(" orbitalinfo->cutoffsvl %f vitaly %f \n",orbitalinfo->cutoffsvl,t_orb->cutoff);
1938                    //
1939                    TClonesArray &tt3 = *torbExtTrk;
1940                    new(tt3[ttentry]) OrbitalInfoTrkVar(*t_orb);
1941                    ttentry++;
1942                    //
1943                    t_orb->Clear();
1944                    //
1945                  }
1946                  //
1947                }
1948              } // end standard tracking algorithm: extended
1949    
1950          } else {          } else {
1951            if ( debug ) printf(" mmm... mode %u standalone  \n",orbitalinfo->mode);            if ( debug ) printf(" mmm... mode %u standalone  \n",orbitalinfo->mode);
1952          }          }
1953          //          //
1954        } else {        } else { // HERE IT MISS ALL CODE FOR EXTENDED TRACKING!
1955          if ( !standalone && tof->ntrk() > 0 ){          if ( !standalone ){
1956            //            //
1957              if ( verbose ) printf(" no orb info for tracks \n");
1958            Int_t nn = 0;            Int_t nn = 0;
1959            for(Int_t nt=0; nt < tof->ntrk(); nt++){              for(Int_t nt=0; nt < tof->ntrk(); nt++){  
1960              //              //
# Line 1664  int OrbitalInfoCore(UInt_t run, TFile *f Line 1980  int OrbitalInfoCore(UInt_t run, TFile *f
1980                //                //
1981                t_orb->Clear();                t_orb->Clear();
1982                //                //
1983              };              }
1984              //              //
1985            };                }
1986          };            //
1987        };        // if( orbitalinfo->TimeGap>0){            // Code for extended tracking algorithm:
1988              //
1989              if ( hasNucleiTrk ){
1990                Int_t ttentry = 0;
1991                if ( verbose ) printf(" hasNucleiTrk \n");
1992                for(Int_t nt=0; nt < tcNucleiTof->GetEntries() ; nt++){  
1993                  //  
1994                  ToFTrkVar *ptt = (ToFTrkVar*)(tcNucleiTof->At(nt));
1995                  if ( ptt->trkseqno != -1  ){
1996                    //
1997                    t_orb->trkseqno = ptt->trkseqno;
1998                    //
1999                    t_orb->Eij = 0;
2000                    //
2001                    t_orb->Sij = 0;
2002                    //          
2003                    t_orb->pitch = -1000.;
2004                    //
2005                    t_orb->sunangle = -1000.;
2006                    //
2007                    t_orb->sunmagangle = -1000;
2008                    //
2009                    t_orb->cutoff = -1000.;
2010                    //
2011                    TClonesArray &tz1 = *torbNucleiTrk;
2012                    new(tz1[ttentry]) OrbitalInfoTrkVar(*t_orb);
2013                    ttentry++;
2014                    //
2015                    t_orb->Clear();
2016                    //
2017                  }
2018                  //
2019                }
2020              }
2021              if ( hasExtNucleiTrk ){
2022                Int_t ttentry = 0;
2023                if ( verbose ) printf(" hasExtNucleiTrk \n");
2024                for(Int_t nt=0; nt < tcExtNucleiTof->GetEntries() ; nt++){  
2025                  //
2026                  if ( verbose ) printf(" got2\n");
2027                  ToFTrkVar *ptt = (ToFTrkVar*)(tcExtNucleiTof->At(nt));
2028                  if ( ptt->trkseqno != -1  ){
2029                    //
2030                    t_orb->trkseqno = ptt->trkseqno;
2031                    //
2032                    t_orb->Eij = 0;
2033                    //
2034                    t_orb->Sij = 0;
2035                    //          
2036                    t_orb->pitch = -1000.;
2037                    //
2038                    t_orb->sunangle = -1000.;
2039                    //
2040                    t_orb->sunmagangle = -1000;
2041                    //
2042                    t_orb->cutoff = -1000.;
2043                    //
2044                    TClonesArray &tz2 = *torbExtNucleiTrk;
2045                    new(tz2[ttentry]) OrbitalInfoTrkVar(*t_orb);
2046                    ttentry++;
2047                    //
2048                    t_orb->Clear();
2049                    //
2050                  }
2051                  //
2052                }
2053              }
2054              if ( hasExtTrk ){
2055                Int_t ttentry = 0;
2056                if ( verbose ) printf(" hasExtTrk \n");
2057                for(Int_t nt=0; nt < tcExtTof->GetEntries() ; nt++){  
2058                  //
2059                  if ( verbose ) printf(" got3\n");
2060                  ToFTrkVar *ptt = (ToFTrkVar*)(tcExtTof->At(nt));
2061                  if ( ptt->trkseqno != -1  ){
2062                    //
2063                    t_orb->trkseqno = ptt->trkseqno;
2064                    //
2065                    t_orb->Eij = 0;
2066                    //
2067                    t_orb->Sij = 0;
2068                    //          
2069                    t_orb->pitch = -1000.;
2070                    //
2071                    t_orb->sunangle = -1000.;
2072                    //
2073                    t_orb->sunmagangle = -1000;
2074                    //
2075                    t_orb->cutoff = -1000.;
2076                    //
2077                    TClonesArray &tz3 = *torbExtNucleiTrk;
2078                    new(tz3[ttentry]) OrbitalInfoTrkVar(*t_orb);
2079                    ttentry++;
2080                    //
2081                    t_orb->Clear();
2082                    //
2083                  }
2084                  //
2085                }
2086              }
2087            }
2088          } // if( orbitalinfo->TimeGap>0){
2089        //        //
2090        // Fill the class        // Fill the class
2091        //        //
# Line 1775  int OrbitalInfoCore(UInt_t run, TFile *f Line 2192  int OrbitalInfoCore(UInt_t run, TFile *f
2192    if ( runinfo ) runinfo->Close();        if ( runinfo ) runinfo->Close();    
2193    if ( runinfo ) delete runinfo;    if ( runinfo ) delete runinfo;
2194    
2195      if ( tcNucleiTrk ){
2196        tcNucleiTrk->Delete();
2197        delete tcNucleiTrk;
2198        tcNucleiTrk = NULL;
2199      }
2200      if ( tcExtNucleiTrk ){
2201        tcExtNucleiTrk->Delete();
2202        delete tcExtNucleiTrk;
2203        tcExtNucleiTrk = NULL;
2204      }
2205      if ( tcExtTrk ){
2206        tcExtTrk->Delete();
2207        delete tcExtTrk;
2208        tcExtTrk = NULL;
2209      }
2210    
2211      if ( tcNucleiTof ){
2212        tcNucleiTof->Delete();
2213        delete tcNucleiTof;
2214        tcNucleiTof = NULL;
2215      }
2216      if ( tcExtNucleiTof ){
2217        tcExtNucleiTof->Delete();
2218        delete tcExtNucleiTof;
2219        tcExtNucleiTof = NULL;
2220      }
2221      if ( tcExtTof ){
2222        tcExtTof->Delete();
2223        delete tcExtTof;
2224        tcExtTrk = NULL;
2225      }
2226    
2227    
2228    if ( tof ) delete tof;    if ( tof ) delete tof;
2229    if ( trke ) delete trke;    if ( trke ) delete trke;
2230    
Legend:
Removed from v.1.72  
changed lines
  Added in v.1.73
  ViewVC Help
Powered by ViewVC 1.1.23