/[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.35 by mocchiut, Thu Dec 11 10:08:19 2008 UTC revision 1.42 by mocchiut, Wed Sep 30 12:22:31 2009 UTC
# Line 44  Line 44 
44  #include <OrbitalInfoCore.h>  #include <OrbitalInfoCore.h>
45  #include <InclinationInfo.h>  #include <InclinationInfo.h>
46    
47    
48  using namespace std;  using namespace std;
49    
50  //  //
# Line 146  int OrbitalInfoCore(UInt_t run, TFile *f Line 147  int OrbitalInfoCore(UInt_t run, TFile *f
147    //    //
148    TFile *l0File = 0;    TFile *l0File = 0;
149    TTree *l0tr = 0;    TTree *l0tr = 0;
150    TTree *l0trm = 0;    //  TTree *l0trm = 0;
151      TChain *ch = 0;
152    // EM: open also header branch    // EM: open also header branch
153    TBranch *l0head = 0;    TBranch *l0head = 0;
154    pamela::EventHeader *eh = 0;    pamela::EventHeader *eh = 0;
# Line 176  int OrbitalInfoCore(UInt_t run, TFile *f Line 178  int OrbitalInfoCore(UInt_t run, TFile *f
178    //    //
179    // IGRF stuff    // IGRF stuff
180    //    //
181    float dimo = 0.0; // dipole moment (computed from dat files)    Float_t dimo = 0.0; // dipole moment (computed from dat files)
182    float bnorth, beast, bdown, babs;    Float_t bnorth, beast, bdown, babs;
183    float xl; // L value    Float_t xl; // L value
184    float icode; // code value for L accuracy (see fortran code)    Float_t icode; // code value for L accuracy (see fortran code)
185    float bab1; // What's  the difference with babs?    Float_t bab1; // What's  the difference with babs?
186    float stps = 0.005; // step size for field line tracing    Float_t stps = 0.005; // step size for field line tracing
187    float bdel = 0.01; // required accuracy    Float_t bdel = 0.01; // required accuracy
188    float bequ;  // equatorial b value (also called b_0)    Float_t bequ;  // equatorial b value (also called b_0)
189    bool value = 0; // false if bequ is not the minimum b value    Bool_t value = 0; // false if bequ is not the minimum b value
190    float rr0; // equatorial radius normalized to earth radius    Float_t rr0; // equatorial radius normalized to earth radius
191    
192    //    //
193    // Working filename    // Working filename
# Line 209  int OrbitalInfoCore(UInt_t run, TFile *f Line 211  int OrbitalInfoCore(UInt_t run, TFile *f
211    OrbitalInfofolder << tempname.str().c_str();    OrbitalInfofolder << tempname.str().c_str();
212    tempname << "/OrbitalInfotree_run";    tempname << "/OrbitalInfotree_run";
213    tempname << run << ".root";      tempname << run << ".root";  
214      UInt_t totnorun = 0;
215    //    //
216    // DB classes    // DB classes
217    //    //
# Line 229  int OrbitalInfoCore(UInt_t run, TFile *f Line 232  int OrbitalInfoCore(UInt_t run, TFile *f
232    Int_t ltp2 = 0;    Int_t ltp2 = 0;
233    Int_t ltp3 = 0;    Int_t ltp3 = 0;
234    Int_t uno = 1;    Int_t uno = 1;
235    char *niente = " ";    const char *niente = " ";
236    GL_PARAM *glparam = new GL_PARAM();    GL_PARAM *glparam = new GL_PARAM();
237    GL_PARAM *glparam2 = new GL_PARAM();    GL_PARAM *glparam2 = new GL_PARAM();
238    Int_t parerror=glparam->Query_GL_PARAM(1,301,dbc); // parameters stored in DB in GL_PRAM table    Int_t parerror=glparam->Query_GL_PARAM(1,301,dbc); // parameters stored in DB in GL_PRAM table
# Line 342  int OrbitalInfoCore(UInt_t run, TFile *f Line 345  int OrbitalInfoCore(UInt_t run, TFile *f
345    // number of run to be processed    // number of run to be processed
346    //    //
347    numbofrun = runinfo->GetNoRun();    numbofrun = runinfo->GetNoRun();
348    UInt_t totnorun = runinfo->GetRunEntries();    totnorun = runinfo->GetRunEntries();
349    //    //
350    // Try to access the OrbitalInfo tree in the file, if it exists we are reprocessing data if not we are processing a new run    // Try to access the OrbitalInfo tree in the file, if it exists we are reprocessing data if not we are processing a new run
351    //    //
# Line 502  int OrbitalInfoCore(UInt_t run, TFile *f Line 505  int OrbitalInfoCore(UInt_t run, TFile *f
505      fname = ftmpname.str().c_str();      fname = ftmpname.str().c_str();
506      ftmpname.str("");      ftmpname.str("");
507      //      //
508      // print out informations      // print nout informations
509      //      //
510      totevent = runinfo->NEVENTS;      totevent = runinfo->NEVENTS;
511      evfrom = runinfo->EV_FROM;      evfrom = runinfo->EV_FROM;
# Line 569  int OrbitalInfoCore(UInt_t run, TFile *f Line 572  int OrbitalInfoCore(UInt_t run, TFile *f
572      ULong_t DeltaOBT = TimeSync - ObtSync;      ULong_t DeltaOBT = TimeSync - ObtSync;
573    
574      if ( debug ) printf(" 2 TimeSync %lu ObtSync %lu DeltaOBT %lu\n",(ULong_t)(dbtime->GetTimesync()/1000),(ULong_t)dbtime->GetObt0(),TimeSync-ObtSync);      if ( debug ) printf(" 2 TimeSync %lu ObtSync %lu DeltaOBT %lu\n",(ULong_t)(dbtime->GetTimesync()/1000),(ULong_t)dbtime->GetObt0(),TimeSync-ObtSync);
575            //
576      l0trm = (TTree*)l0File->Get("Mcmd");      // Read MCMDs from up to 11 files, 5 before and 5 after the present one in order to have some kind of inclination information
577      neventsm = l0trm->GetEntries();      //
578        ch = new TChain("Mcmd","Mcmd");
579        //
580        // look in the DB to find the closest files to this run
581        //
582        TSQLResult *pResult = 0;
583        TSQLRow *Row = 0;
584        stringstream myquery;
585        UInt_t l0fid[10];
586        Int_t i = 0;
587        memset(l0fid,0,10*sizeof(Int_t));
588        //
589        myquery.str("");
590        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;";
591        //
592        pResult = dbc->Query(myquery.str().c_str());
593        //
594        i = 9;
595        if( pResult ){
596          //
597          Row = pResult->Next();
598          //
599          while ( Row ){
600            //
601            // store infos and exit
602            //
603            l0fid[i] = (UInt_t)atoll(Row->GetField(0));
604            i--;
605            Row = pResult->Next();  
606            //
607          };
608          pResult->Delete();
609        };
610        //
611        myquery.str("");
612        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;";
613        //
614        pResult = dbc->Query(myquery.str().c_str());
615        //
616        i = 0;
617        if( pResult ){
618          //
619          Row = pResult->Next();
620          //
621          while ( Row ){
622            //
623            // store infos and exit
624            //
625            l0fid[i] = (UInt_t)atoll(Row->GetField(0));
626            i++;
627            Row = pResult->Next();  
628            //
629          };
630          pResult->Delete();
631        };
632        //
633        i = 0;
634        UInt_t previd = 0;
635        while ( i < 10 ){
636          if ( l0fid[i] && previd != l0fid[i] ){
637            previd = l0fid[i];
638            myquery.str("");
639            myquery << "select PATH,NAME from GL_ROOT where ID=" << l0fid[i] << " ;";
640            //
641            pResult = dbc->Query(myquery.str().c_str());
642            //
643            if( pResult ){
644              //
645              Row = pResult->Next();
646              //
647              if ( debug ) printf(" Using inclination informations from file: %s \n",(((TString)gSystem->ExpandPathName(Row->GetField(0)))+"/"+(TString)Row->GetField(1)).Data());
648              ch->Add(((TString)gSystem->ExpandPathName(Row->GetField(0)))+"/"+(TString)Row->GetField(1));
649              //
650              pResult->Delete();
651            };
652          };
653          i++;
654        };
655        //
656        //    l0trm = (TTree*)l0File->Get("Mcmd");
657        //    ch->ls();
658        ch->SetBranchAddress("Mcmd",&mcmdev);
659        //    printf(" entries %llu \n", ch->GetEntries());
660        //    l0trm = ch->GetTree();
661        //    neventsm = l0trm->GetEntries();
662        neventsm = ch->GetEntries();
663        if ( debug ) printf(" entries %u \n", neventsm);
664      //    neventsm = 0;      //    neventsm = 0;
665      //      //
666      if (neventsm == 0){      if (neventsm == 0){
# Line 582  int OrbitalInfoCore(UInt_t run, TFile *f Line 671  int OrbitalInfoCore(UInt_t run, TFile *f
671      }      }
672      //      //
673            
674      l0trm->SetBranchAddress("Mcmd", &mcmdev);      //    l0trm->SetBranchAddress("Mcmd", &mcmdev);
675      //    l0trm->SetBranchAddress("Header", &eh);      //    l0trm->SetBranchAddress("Header", &eh);
676      //      //
677      //      //
# Line 609  int OrbitalInfoCore(UInt_t run, TFile *f Line 698  int OrbitalInfoCore(UInt_t run, TFile *f
698      //      //
699       //       //
700      for ( re = runinfo->EV_FROM; re < (runinfo->EV_FROM+runinfo->NEVENTS); re++){      for ( re = runinfo->EV_FROM; re < (runinfo->EV_FROM+runinfo->NEVENTS); re++){
       
701        //        //
702        if ( procev%1000 == 0 && procev > 0 && verbose ) printf(" %iK \n",procev/1000);          if ( procev%1000 == 0 && procev > 0 && verbose ) printf(" %iK \n",procev/1000);  
703        if ( debug ) printf(" %i \n",procev);              if ( debug ) printf(" %i \n",procev);      
# Line 620  int OrbitalInfoCore(UInt_t run, TFile *f Line 708  int OrbitalInfoCore(UInt_t run, TFile *f
708        //        //
709        ph = eh->GetPscuHeader();        ph = eh->GetPscuHeader();
710        atime = dbtime->DBabsTime(ph->GetOrbitalTime());        atime = dbtime->DBabsTime(ph->GetOrbitalTime());
711          if ( debug ) printf(" %i absolute time \n",procev);      
712        //        //
713        // paranoid check        // paranoid check
714        //        //
# Line 659  int OrbitalInfoCore(UInt_t run, TFile *f Line 748  int OrbitalInfoCore(UInt_t run, TFile *f
748        //        //
749        // start processing        // start processing
750        //        //
751          if ( debug ) printf(" %i start processing \n",procev);      
752        orbitalinfo->Clear();        orbitalinfo->Clear();
753        //        //
754        OrbitalInfoTrkVar *t_orb = new OrbitalInfoTrkVar();        OrbitalInfoTrkVar *t_orb = new OrbitalInfoTrkVar();
# Line 670  int OrbitalInfoCore(UInt_t run, TFile *f Line 760  int OrbitalInfoCore(UInt_t run, TFile *f
760        orbitalinfo->pkt_num = ph->GetCounter();        orbitalinfo->pkt_num = ph->GetCounter();
761        orbitalinfo->OBT = ph->GetOrbitalTime();        orbitalinfo->OBT = ph->GetOrbitalTime();
762        orbitalinfo->absTime = atime;        orbitalinfo->absTime = atime;
763          if ( debug ) printf(" %i pktnum obt abstime \n",procev);      
764        //        //
765        // Propagate the orbit from the tle time to atime, using SGP(D)4.        // Propagate the orbit from the tle time to atime, using SGP(D)4.
766        //        //
767          if ( debug ) printf(" %i sgp4 \n",procev);      
768        cCoordGeo coo;        cCoordGeo coo;
769        float jyear=0;            Float_t jyear=0.;    
770        //        //
771        if(atime >= gltle->GetToTime()) {        if(atime >= gltle->GetToTime()) {
772          if ( !gltle->Query(atime, dbc) ){          if ( !gltle->Query(atime, dbc) ){
773            //                  //      
774            // Compute the magnetic dipole moment.            // Compute the magnetic dipole moment.
775            //            //
776              if ( debug ) printf(" %i compute magnetic dipole moment \n",procev);      
777            UInt_t year, month, day, hour, min, sec;            UInt_t year, month, day, hour, min, sec;
778            //            //
779            TTimeStamp t = TTimeStamp(atime, kTRUE);            TTimeStamp t = TTimeStamp(atime, kTRUE);
# Line 690  int OrbitalInfoCore(UInt_t run, TFile *f Line 783  int OrbitalInfoCore(UInt_t run, TFile *f
783              + (month*31.+ (float) day)/365.              + (month*31.+ (float) day)/365.
784              + (hour*3600.+min*60.+(float)sec)/(24*3600*365.);              + (hour*3600.+min*60.+(float)sec)/(24*3600*365.);
785            //            //
786              if ( debug ) printf(" %i compute magnetic dipole moment get dipole moment for year\n",procev);      
787            feldcof_(&jyear, &dimo); // get dipole moment for year            feldcof_(&jyear, &dimo); // get dipole moment for year
788              if ( debug ) printf(" %i compute magnetic dipole moment end\n",procev);      
789          } else {          } else {
790            code = -56;            code = -56;
791            goto closeandexit;            goto closeandexit;
# Line 713  int OrbitalInfoCore(UInt_t run, TFile *f Line 808  int OrbitalInfoCore(UInt_t run, TFile *f
808          upperqtime = atime;          upperqtime = atime;
809          lowerqtime = runinfo->RUNHEADER_TIME;          lowerqtime = runinfo->RUNHEADER_TIME;
810          for ( ik = 0; ik < neventsm; ik++){          for ( ik = 0; ik < neventsm; ik++){
811            l0trm->GetEntry(ik);            //      l0trm->GetEntry(ik);
812              ch->GetEntry(ik);
813            tmpSize = mcmdev->Records->GetEntries();            tmpSize = mcmdev->Records->GetEntries();
814            numrec = tmpSize;            numrec = tmpSize;
815            for (Int_t j3 = 0;j3<tmpSize;j3++){            for (Int_t j3 = 0;j3<tmpSize;j3++){
816              if ( debug ) printf(" eh eh eh \n");              if ( debug ) printf(" ik %i j3 %i eh eh eh \n",ik,j3);
817              mcmdrc = (pamela::McmdRecord*)mcmdev->Records->At(j3);              mcmdrc = (pamela::McmdRecord*)mcmdev->Records->At(j3);
818              if ((int)mcmdrc->ID1 == 226){              if ( mcmdrc ){ // missing inclination bug [8RED 090116]
819                L_QQ_Q_l_upper->fill(mcmdrc->McmdData);                if ((int)mcmdrc->ID1 == 226 && mcmdrc->Mcmd_Block_crc_ok == 1){
820                for (UInt_t ui = 0; ui < 6; ui++){                  L_QQ_Q_l_upper->fill(mcmdrc->McmdData);
821                  if (ui>0){                  for (UInt_t ui = 0; ui < 6; ui++){
822                    if (L_QQ_Q_l_upper->time[ui]>L_QQ_Q_l_upper->time[0]){                    if (ui>0){
823                      if (dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[ui]*1000-DeltaOBT*1000))<atime){                      if (L_QQ_Q_l_upper->time[ui]>L_QQ_Q_l_upper->time[0]){
824                        upperqtime = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[ui]*1000-DeltaOBT*1000));                        if (dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[ui]*1000-DeltaOBT*1000))<atime){
825                            upperqtime = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[ui]*1000-DeltaOBT*1000));
826                            orbits.getPosition((double) (upperqtime - gltle->GetFromTime())/60., &eCi);
827                            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]);
828                          }else {
829                            lowerqtime = upperqtime;
830                            upperqtime = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[ui]*1000-DeltaOBT*1000));
831                            orbits.getPosition((double) (upperqtime - gltle->GetFromTime())/60., &eCi);
832                            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]);
833                            mcreen = j3;
834                            mctren = ik;
835                            if(fgh==0){
836                              CopyQ(L_QQ_Q_l_lower,L_QQ_Q_l_upper);
837                              CopyAng(RYPang_lower,RYPang_upper);
838                            }
839                            oi=ui;
840                            goto closethisloop;
841                          }
842                          fgh++;
843                          CopyQ(L_QQ_Q_l_lower,L_QQ_Q_l_upper);
844                          CopyAng(RYPang_lower,RYPang_upper);
845                        }
846                      }else{
847                        if (dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000))<atime){
848                          upperqtime = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
849                        orbits.getPosition((double) (upperqtime - gltle->GetFromTime())/60., &eCi);                        orbits.getPosition((double) (upperqtime - gltle->GetFromTime())/60., &eCi);
850                        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]);                        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]);
851                      }else {                      }
852                        else {
853                        lowerqtime = upperqtime;                        lowerqtime = upperqtime;
854                        upperqtime = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[ui]*1000-DeltaOBT*1000));                        upperqtime = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
855                        orbits.getPosition((double) (upperqtime - gltle->GetFromTime())/60., &eCi);                        orbits.getPosition((double) (upperqtime - gltle->GetFromTime())/60., &eCi);
856                        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]);                        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]);
857                        mcreen = j3;                        mcreen = j3;
858                        mctren = ik;                        mctren = ik;
859                        if(fgh==0){                        if(fgh==0){
860                          CopyQ(L_QQ_Q_l_lower,L_QQ_Q_l_upper);                          CopyQ(L_QQ_Q_l_lower,L_QQ_Q_l_upper);
861                          CopyAng(RYPang_lower,RYPang_upper);                          CopyAng(RYPang_lower,RYPang_upper);
862                            lowerqtime = atime-1;
863                        }                        }
864                        oi=ui;                        oi=ui;
865                        goto closethisloop;                        goto closethisloop;
866                          //_0 = true;
867                      }                      }
868                      fgh++;                      fgh++;
869                      CopyQ(L_QQ_Q_l_lower,L_QQ_Q_l_upper);                      CopyQ(L_QQ_Q_l_lower,L_QQ_Q_l_upper);
870                      CopyAng(RYPang_lower,RYPang_upper);                      CopyAng(RYPang_lower,RYPang_upper);
                   }  
                 }else{  
                   if (dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000))<atime){  
                     upperqtime = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));  
                     orbits.getPosition((double) (upperqtime - gltle->GetFromTime())/60., &eCi);  
                     RYPang_upper->TransAngle(eCi.getPos().m_x,eCi.getPos().m_y,eCi.getPos().m_z,eCi.getVel().m_x,eCi.getVel().m_y,eCi.getVel().m_z,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]);  
                   }  
                   else {  
                     lowerqtime = upperqtime;  
                     upperqtime = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));  
                     orbits.getPosition((double) (upperqtime - gltle->GetFromTime())/60., &eCi);  
                     RYPang_upper->TransAngle(eCi.getPos().m_x,eCi.getPos().m_y,eCi.getPos().m_z,eCi.getVel().m_x,eCi.getVel().m_y,eCi.getVel().m_z,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]);  
                     mcreen = j3;  
                     mctren = ik;  
                     if(fgh==0){  
                       CopyQ(L_QQ_Q_l_lower,L_QQ_Q_l_upper);  
                       CopyAng(RYPang_lower,RYPang_upper);  
                       lowerqtime = atime-1;  
                     }  
                     oi=ui;  
                     goto closethisloop;  
871                      //_0 = true;                      //_0 = true;
872                    }                    };
873                    fgh++;                    //cin>>grib;
                   CopyQ(L_QQ_Q_l_lower,L_QQ_Q_l_upper);  
                   CopyAng(RYPang_lower,RYPang_upper);  
                   //_0 = true;  
874                  };                  };
                 //cin>>grib;  
875                };                };
876              };              };
877            };            };
# Line 789  int OrbitalInfoCore(UInt_t run, TFile *f Line 887  int OrbitalInfoCore(UInt_t run, TFile *f
887          lowerqtime = upperqtime;          lowerqtime = upperqtime;
888          Long64_t maxloop = 100000000LL;          Long64_t maxloop = 100000000LL;
889          Long64_t mn = 0;          Long64_t mn = 0;
890          bool gh=false;          Bool_t gh=false;
891          ooi=oi;          ooi=oi;
892          if ( verbose ) printf(" OrbitalInfoCore: sync with quaternions data upperqtime %u lowerqtime %u atime %u \n",(UInt_t)upperqtime,(UInt_t)lowerqtime,atime);          if ( verbose ) printf(" OrbitalInfoCore: sync with quaternions data upperqtime %u lowerqtime %u atime %u \n",(UInt_t)upperqtime,(UInt_t)lowerqtime,atime);
893          while (!gh){                while (!gh){      
# Line 807  int OrbitalInfoCore(UInt_t run, TFile *f Line 905  int OrbitalInfoCore(UInt_t run, TFile *f
905              if (mcreen == numrec){              if (mcreen == numrec){
906                mctren++;                mctren++;
907                mcreen = 0;                mcreen = 0;
908                l0trm->GetEntry(mctren);                //              l0trm->GetEntry(mctren);
909                  ch->GetEntry(mctren);
910                numrec = mcmdev->Records->GetEntries();                numrec = mcmdev->Records->GetEntries();
911              }              }
912              CopyQ(L_QQ_Q_l_lower,L_QQ_Q_l_upper);              CopyQ(L_QQ_Q_l_lower,L_QQ_Q_l_upper);
913              CopyAng(RYPang_lower,RYPang_upper);              CopyAng(RYPang_lower,RYPang_upper);
914              mcmdrc = (pamela::McmdRecord*)mcmdev->Records->At(mcreen);              mcmdrc = (pamela::McmdRecord*)mcmdev->Records->At(mcreen);
915              if ((int)mcmdrc->ID1 == 226){              if ( mcmdrc ){ // missing inclination bug [8RED 090116]
916                L_QQ_Q_l_upper->fill(mcmdrc->McmdData);                if ((int)mcmdrc->ID1 == 226 && mcmdrc->Mcmd_Block_crc_ok == 1){
917                upperqtime = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));                  L_QQ_Q_l_upper->fill(mcmdrc->McmdData);
918                if (upperqtime<lowerqtime){                  upperqtime = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[0]*1000-DeltaOBT*1000));
919                  upperqtime=runinfo->RUNTRAILER_TIME;                  if (upperqtime<lowerqtime){
920                  CopyQ(L_QQ_Q_l_upper,L_QQ_Q_l_lower);                    upperqtime=runinfo->RUNTRAILER_TIME;
921                  CopyAng(RYPang_upper,RYPang_lower);                    CopyQ(L_QQ_Q_l_upper,L_QQ_Q_l_lower);
922                }else{                    CopyAng(RYPang_upper,RYPang_lower);
923                  orbits.getPosition((double) (upperqtime - gltle->GetFromTime())/60., &eCi);                  }else{
924                  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]);                    orbits.getPosition((double) (upperqtime - gltle->GetFromTime())/60., &eCi);
925                      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]);
926                    }
927                    //            re--;
928                    gh=true;
929                }                }
930                //              re--;              };
               gh=true;  
             }  
931            }else{            }else{
932              if ((Int_t)L_QQ_Q_l_upper->time[oi]>(Int_t)L_QQ_Q_l_upper->time[0]){              if ((Int_t)L_QQ_Q_l_upper->time[oi]>(Int_t)L_QQ_Q_l_upper->time[0]){
933                upperqtime = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000));                upperqtime = dbtime->DBabsTime((UInt_t)(L_QQ_Q_l_upper->time[oi]*1000-DeltaOBT*1000));
# Line 985  int OrbitalInfoCore(UInt_t run, TFile *f Line 1086  int OrbitalInfoCore(UInt_t run, TFile *f
1086          //          //
1087        };              };      
1088        //        //
1089          if ( debug ) printf(" pitch angle \n");
1090          //
1091        // pitch angles        // pitch angles
1092        //        //
1093        if ( orbitalinfo->mode != 10 && orbitalinfo->mode != 5 && orbitalinfo->mode !=7 && orbitalinfo->mode != 9 ){        if ( orbitalinfo->mode != 10 && orbitalinfo->mode != 5 && orbitalinfo->mode !=7 && orbitalinfo->mode != 9 ){
# Line 1021  int OrbitalInfoCore(UInt_t run, TFile *f Line 1124  int OrbitalInfoCore(UInt_t run, TFile *f
1124              Double_t E22z = zin[3];              Double_t E22z = zin[3];
1125              if ( (E11x < 100. && E11y < 100. && E22x < 100. && E22y < 100.) || ptt->trkseqno != -1  ){              if ( (E11x < 100. && E11y < 100. && E22x < 100. && E22y < 100.) || ptt->trkseqno != -1  ){
1126                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));
1127                Double_t MyAzim = TMath::RadToDeg()*atan(TMath::Abs(E22y-E11y)/TMath::Abs(E22x-E11x));  //            Double_t MyAzim = TMath::RadToDeg()*atan(TMath::Abs(E22y-E11y)/TMath::Abs(E22x-E11x));
1128                if(E22x-E11x>=0 && E22y-E11y <0) MyAzim =  360. - MyAzim;  //            if(E22x-E11x>=0 && E22y-E11y <0) MyAzim =  360. - MyAzim;
1129                if(E22x-E11x>=0 && E22y-E11y >=0) MyAzim = MyAzim;  //            if(E22x-E11x>=0 && E22y-E11y >=0) MyAzim = MyAzim;
1130                if(E22x-E11x<0 && E22y-E11y >0) MyAzim = 180. - MyAzim;  //            if(E22x-E11x<0 && E22y-E11y >0) MyAzim = 180. - MyAzim;
1131                if(E22x-E11x<0 && E22y-E11y <0) MyAzim = 180. + MyAzim;  //            if(E22x-E11x<0 && E22y-E11y <0) MyAzim = 180. + MyAzim;
1132                Px = (E22x-E11x)/norm;                Px = (E22x-E11x)/norm;
1133                Py = (E22y-E11y)/norm;                Py = (E22y-E11y)/norm;
1134                Pz = (E22z-E11z)/norm;                Pz = (E22z-E11z)/norm;
# Line 1047  int OrbitalInfoCore(UInt_t run, TFile *f Line 1150  int OrbitalInfoCore(UInt_t run, TFile *f
1150                //                //
1151                t_orb->cutoff = 59.3/(pow(orbitalinfo->L,2)*pow((1+sqrt(1-pow(orbitalinfo->L,-3/2)*cos(omega))),2));                t_orb->cutoff = 59.3/(pow(orbitalinfo->L,2)*pow((1+sqrt(1-pow(orbitalinfo->L,-3/2)*cos(omega))),2));
1152                //                //
1153                  if ( t_orb->pitch != t_orb->pitch ) t_orb->pitch = -1000.;
1154                  if ( t_orb->cutoff != t_orb->cutoff ) t_orb->cutoff = -1000.;
1155                  //
1156                if ( debug ) printf(" orbitalinfo->cutoffsvl %f vitaly %f \n",orbitalinfo->cutoffsvl,t_orb->cutoff);                if ( debug ) printf(" orbitalinfo->cutoffsvl %f vitaly %f \n",orbitalinfo->cutoffsvl,t_orb->cutoff);
1157                //                //
1158                new(tor[nn]) OrbitalInfoTrkVar(*t_orb);                new(tor[nn]) OrbitalInfoTrkVar(*t_orb);
# Line 1217  void CopyAng(InclinationInfo *A1, Inclin Line 1323  void CopyAng(InclinationInfo *A1, Inclin
1323  UInt_t holeq(Double_t lower,Double_t upper,Quaternions *Qlower, Quaternions *Qupper, UInt_t f){  UInt_t holeq(Double_t lower,Double_t upper,Quaternions *Qlower, Quaternions *Qupper, UInt_t f){
1324        
1325    UInt_t hole = 10;    UInt_t hole = 10;
1326    bool R10l = false;     // Sign of R10 mode in lower quaternions array    Bool_t R10l = false;     // Sign of R10 mode in lower quaternions array
1327    bool R10u = false;     // Sign of R10 mode in upper quaternions array    Bool_t R10u = false;     // Sign of R10 mode in upper quaternions array
1328    bool insm = false;     // Sign that we inside quaternions array    Bool_t insm = false;     // Sign that we inside quaternions array
1329    bool mxtml = false;    // Sign of mixt mode in lower quaternions array    Bool_t mxtml = false;    // Sign of mixt mode in lower quaternions array
1330    bool mxtmu = false;    // Sign of mixt mode in upper quaternions array    Bool_t mxtmu = false;    // Sign of mixt mode in upper quaternions array
1331    bool npasm = false;     // Sign of normall pass between R10 and non R10 or between non R10 and R10    Bool_t npasm = false;     // Sign of normall pass between R10 and non R10 or between non R10 and R10
1332    UInt_t NCQl = 6;       // Number of correct quaternions in lower array    UInt_t NCQl = 6;       // Number of correct quaternions in lower array
1333    UInt_t NCQu = 6;       // Number of correct quaternions in upper array    UInt_t NCQu = 6;       // Number of correct quaternions in upper array
1334    if (f>0){    if (f>0){
Legend:
Removed from v.1.35  
changed lines
  Added in v.1.42
  ViewVC Help
Powered by ViewVC 1.1.23