/[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.17 by mocchiut, Fri Apr 27 10:35:35 2007 UTC revision 1.35 by mocchiut, Thu Dec 11 10:08:19 2008 UTC
# Line 50  using namespace std; Line 50  using namespace std;
50  // CORE ROUTINE  // CORE ROUTINE
51  //  //
52  //  //
53  int OrbitalInfoCore(UInt_t run, TFile *file, TSQLServer *dbc, Int_t OrbitalInfoargc, char *OrbitalInfoargv[]){  int OrbitalInfoCore(UInt_t run, TFile *file, GL_TABLES *glt, Int_t OrbitalInfoargc, char *OrbitalInfoargv[]){
54    //    //
55    Int_t i = 0;    Int_t i = 0;
56      TString host = glt->CGetHost();
57      TString user = glt->CGetUser();
58      TString psw = glt->CGetPsw();
59      TSQLServer *dbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data());
60      //
61      stringstream myquery;
62      myquery.str("");
63      myquery << "SET time_zone='+0:00'";
64      dbc->Query(myquery.str().c_str());
65    //    //
66    TString processFolder = Form("OrbitalInfoFolder_%u",run);    TString processFolder = Form("OrbitalInfoFolder_%u",run);
67    //    //
# Line 61  int OrbitalInfoCore(UInt_t run, TFile *f Line 70  int OrbitalInfoCore(UInt_t run, TFile *f
70    Bool_t debug = false;    Bool_t debug = false;
71    //    //
72    Bool_t verbose = false;    Bool_t verbose = false;
73      //
74      Bool_t standalone = false;
75      //
76    if ( OrbitalInfoargc > 0 ){    if ( OrbitalInfoargc > 0 ){
77      i = 0;      i = 0;
78      while ( i < OrbitalInfoargc ){      while ( i < OrbitalInfoargc ){
# Line 79  int OrbitalInfoCore(UInt_t run, TFile *f Line 90  int OrbitalInfoCore(UInt_t run, TFile *f
90        if ( (!strcmp(OrbitalInfoargv[i],"--verbose")) || (!strcmp(OrbitalInfoargv[i],"-v")) ) {        if ( (!strcmp(OrbitalInfoargv[i],"--verbose")) || (!strcmp(OrbitalInfoargv[i],"-v")) ) {
91          verbose = true;          verbose = true;
92        };        };
93          if ( (!strcmp(OrbitalInfoargv[i],"--standalone")) ) {
94            standalone = true;
95          };
96          if ( (!strcmp(OrbitalInfoargv[i],"--calculate-pitch")) ) {
97            standalone = false;
98          };
99        i++;        i++;
100      };      };
101    };    };
# Line 184  int OrbitalInfoCore(UInt_t run, TFile *f Line 201  int OrbitalInfoCore(UInt_t run, TFile *f
201    TTree *tempOrbitalInfo = 0;    TTree *tempOrbitalInfo = 0;
202    stringstream tempname;    stringstream tempname;
203    stringstream OrbitalInfofolder;    stringstream OrbitalInfofolder;
204      Bool_t myfold = false;
205    tempname.str("");    tempname.str("");
206    tempname << outDir;    tempname << outDir;
207    tempname << "/" << processFolder.Data();    tempname << "/" << processFolder.Data();
208    OrbitalInfofolder.str("");    OrbitalInfofolder.str("");
209    OrbitalInfofolder << tempname.str().c_str();    OrbitalInfofolder << tempname.str().c_str();
   gSystem->MakeDirectory(OrbitalInfofolder.str().c_str());  
210    tempname << "/OrbitalInfotree_run";    tempname << "/OrbitalInfotree_run";
211    tempname << run << ".root";      tempname << run << ".root";  
212    //    //
# Line 216  int OrbitalInfoCore(UInt_t run, TFile *f Line 233  int OrbitalInfoCore(UInt_t run, TFile *f
233    GL_PARAM *glparam = new GL_PARAM();    GL_PARAM *glparam = new GL_PARAM();
234    GL_PARAM *glparam2 = new GL_PARAM();    GL_PARAM *glparam2 = new GL_PARAM();
235    Int_t parerror=glparam->Query_GL_PARAM(1,301,dbc); // parameters stored in DB in GL_PRAM table    Int_t parerror=glparam->Query_GL_PARAM(1,301,dbc); // parameters stored in DB in GL_PRAM table
236        //
237      // Orientation variables
238      //
239      UInt_t evfrom = 0;
240      UInt_t jumped = 0;
241      Int_t itr = -1;    
242      Double_t A1;
243      Double_t A2;
244      Double_t A3;
245      Double_t Px = 0;
246      Double_t Py = 0;      
247      Double_t Pz = 0;  
248      TTree *ttof = 0;
249      ToFLevel2 *tof = new ToFLevel2();
250      OrientationInfo *PO = new OrientationInfo();
251      Int_t nz = 6;
252      Float_t zin[6];
253      Int_t nevtofl2 = 0;
254      //  
255    if ( parerror<0 ) {    if ( parerror<0 ) {
256      code = parerror;      code = parerror;
257      goto closeandexit;      goto closeandexit;
# Line 236  int OrbitalInfoCore(UInt_t run, TFile *f Line 271  int OrbitalInfoCore(UInt_t run, TFile *f
271    //    //
272    // End IGRF stuff//    // End IGRF stuff//
273    //    //
274      for (Int_t ip=0;ip<nz;ip++){
275        zin[ip] = tof->GetZTOF(tof->GetToFPlaneID(ip));
276      };
277      //
278      if ( !standalone ){
279        //
280        // Does it contain the Tracker tree?
281        //
282        ttof = (TTree*)file->Get("ToF");
283        if ( !ttof ) {
284          if ( verbose ) printf(" OrbitalInfo - ERROR: no tof tree\n");
285          code = -900;
286          goto closeandexit;
287        };
288        ttof->SetBranchAddress("ToFLevel2",&tof);  
289        nevtofl2 = ttof->GetEntries();
290      };
291    //    //
292    // Let's start!    // Let's start!
293    //    //
# Line 291  int OrbitalInfoCore(UInt_t run, TFile *f Line 342  int OrbitalInfoCore(UInt_t run, TFile *f
342    // number of run to be processed    // number of run to be processed
343    //    //
344    numbofrun = runinfo->GetNoRun();    numbofrun = runinfo->GetNoRun();
345      UInt_t totnorun = runinfo->GetRunEntries();
346    //    //
347    // Try to access the OrbitalInfo tree in the file, if it exists we are reprocessing data if not we are processing a new run    // Try to access the OrbitalInfo tree in the file, if it exists we are reprocessing data if not we are processing a new run
348    //    //
# Line 317  int OrbitalInfoCore(UInt_t run, TFile *f Line 369  int OrbitalInfoCore(UInt_t run, TFile *f
369      //      //
370      if (verbose) printf("\n Preparing the pre-processing...\n");      if (verbose) printf("\n Preparing the pre-processing...\n");
371      //      //
372      if ( run == 0 ){      if ( run == 0 || totnorun == 1 ){
373        //        //
374        // we are reprocessing all the file        // we are reprocessing all the file
375        // if we are reprocessing everything we don't need to copy any old event and we can just work with the new tree and delete the old one immediately        // if we are reprocessing everything we don't need to copy any old event and we can just work with the new tree and delete the old one immediately
# Line 336  int OrbitalInfoCore(UInt_t run, TFile *f Line 388  int OrbitalInfoCore(UInt_t run, TFile *f
388        //        //
389        // copying old tree to a new file        // copying old tree to a new file
390        //        //
391          gSystem->MakeDirectory(OrbitalInfofolder.str().c_str());
392          myfold = true;
393        tempfile = new TFile(tempname.str().c_str(),"RECREATE");        tempfile = new TFile(tempname.str().c_str(),"RECREATE");
394        tempOrbitalInfo = OrbitalInfotrclone->CloneTree(-1,"fast");        tempOrbitalInfo = OrbitalInfotrclone->CloneTree(-1,"fast");
395        tempOrbitalInfo->SetName("OrbitalInfo-old");        tempOrbitalInfo->SetName("OrbitalInfo-old");
# Line 356  int OrbitalInfoCore(UInt_t run, TFile *f Line 410  int OrbitalInfoCore(UInt_t run, TFile *f
410    file->cd();    file->cd();
411    OrbitalInfotr = new TTree("OrbitalInfo-new","PAMELA OrbitalInfo data");    OrbitalInfotr = new TTree("OrbitalInfo-new","PAMELA OrbitalInfo data");
412    OrbitalInfotr->SetAutoSave(900000000000000LL);    OrbitalInfotr->SetAutoSave(900000000000000LL);
413      orbitalinfo->Set();//ELENA **TEMPORANEO?**
414    OrbitalInfotr->Branch("OrbitalInfo","OrbitalInfo",&orbitalinfo);    OrbitalInfotr->Branch("OrbitalInfo","OrbitalInfo",&orbitalinfo);
415    //    //
416    if ( reproc && !reprocall ){    if ( reproc && !reprocall ){
# Line 398  int OrbitalInfoCore(UInt_t run, TFile *f Line 453  int OrbitalInfoCore(UInt_t run, TFile *f
453    //    //
454    // Loop over the run to be processed    // Loop over the run to be processed
455    //    //
     
456    for (UInt_t irun=0; irun < numbofrun; irun++){    for (UInt_t irun=0; irun < numbofrun; irun++){
457      //      //
458      // retrieve the first run ID to be processed using the RunInfo list      // retrieve the first run ID to be processed using the RunInfo list
# Line 451  int OrbitalInfoCore(UInt_t run, TFile *f Line 505  int OrbitalInfoCore(UInt_t run, TFile *f
505      // print out informations      // print out informations
506      //      //
507      totevent = runinfo->NEVENTS;      totevent = runinfo->NEVENTS;
508        evfrom = runinfo->EV_FROM;
509      //cout<<"totevents = "<<totevent<<"\n";      //cout<<"totevents = "<<totevent<<"\n";
510      if (verbose){      if (verbose){
511        printf("\n LEVEL0 data file: %s \n",fname.Data());        printf("\n LEVEL0 data file: %s \n",fname.Data());
# Line 458  int OrbitalInfoCore(UInt_t run, TFile *f Line 513  int OrbitalInfoCore(UInt_t run, TFile *f
513        printf(" RUN TRAILER absolute time is: %u \n",runinfo->RUNTRAILER_TIME);        printf(" RUN TRAILER absolute time is: %u \n",runinfo->RUNTRAILER_TIME);
514        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);        printf(" %i events to be processed for run %u: from %i to %i \n\n",totevent,idRun,runinfo->EV_FROM+1,runinfo->EV_FROM+totevent);
515      }//      }//
516        //
517        //    if ( !totevent ) goto closeandexit;
518      // Open Level0 file      // Open Level0 file
519      l0File = new TFile(fname.Data());      l0File = new TFile(fname.Data());
520      if ( !l0File ) {      if ( !l0File ) {
# Line 484  int OrbitalInfoCore(UInt_t run, TFile *f Line 541  int OrbitalInfoCore(UInt_t run, TFile *f
541      // end EM      // end EM
542      nevents = l0head->GetEntries();      nevents = l0head->GetEntries();
543      //      //
544      if ( nevents < 1 ) {      if ( nevents < 1 && totevent ) {
545        if ( debug ) printf(" OrbitalInfo - ERROR: Level0 file is empty\n\n");        if ( debug ) printf(" OrbitalInfo - ERROR: Level0 file is empty\n\n");
546        l0File->Close();        l0File->Close();
547        code = -11;        code = -11;
548        goto closeandexit;        goto closeandexit;
549      };      };
550      //      //
551      if ( runinfo->EV_TO > nevents-1 ) {      if ( runinfo->EV_TO > nevents-1 && totevent ) {
552        if ( debug ) printf(" OrbitalInfo - ERROR: too few entries in the registry tree\n");        if ( debug ) printf(" OrbitalInfo - ERROR: too few entries in the registry tree\n");
553        l0File->Close();        l0File->Close();
554        code = -12;        code = -12;
# Line 515  int OrbitalInfoCore(UInt_t run, TFile *f Line 572  int OrbitalInfoCore(UInt_t run, TFile *f
572            
573      l0trm = (TTree*)l0File->Get("Mcmd");      l0trm = (TTree*)l0File->Get("Mcmd");
574      neventsm = l0trm->GetEntries();      neventsm = l0trm->GetEntries();
575        //    neventsm = 0;
576      //      //
577      if (neventsm == 0){      if (neventsm == 0){
578        if ( debug ) printf("InclinationInfo - ERROR: No quaternions in this File");        if ( debug ) printf("InclinationInfo - WARNING: No quaternions in this File");
579        l0File->Close();        //      l0File->Close();
580        code = -13;        code = 900;
581        goto closeandexit;        //      goto closeandexit;
582      }      }
583      //      //
584            
585      l0trm->SetBranchAddress("Mcmd", &mcmdev);      l0trm->SetBranchAddress("Mcmd", &mcmdev);
586      l0trm->SetBranchAddress("Header", &eh);      //    l0trm->SetBranchAddress("Header", &eh);
587      //      //
588      //      //
589      //      //
# Line 549  int OrbitalInfoCore(UInt_t run, TFile *f Line 607  int OrbitalInfoCore(UInt_t run, TFile *f
607      //      //
608      if (verbose) printf("\n Ready to start! \n\n Processed events: \n\n");      if (verbose) printf("\n Ready to start! \n\n Processed events: \n\n");
609      //      //
610         //
611      for ( re = runinfo->EV_FROM; re < (runinfo->EV_FROM+runinfo->NEVENTS); re++){      for ( re = runinfo->EV_FROM; re < (runinfo->EV_FROM+runinfo->NEVENTS); re++){
612            
613        //        //
# Line 564  int OrbitalInfoCore(UInt_t run, TFile *f Line 623  int OrbitalInfoCore(UInt_t run, TFile *f
623        //        //
624        // paranoid check        // paranoid check
625        //        //
626        if ( (atime > runinfo->RUNTRAILER_TIME) || (atime < runinfo->RUNHEADER_TIME)  ) {        if ( (atime > (runinfo->RUNTRAILER_TIME+1)) || (atime < (runinfo->RUNHEADER_TIME-1))  ) {
627          if (verbose) printf(" OrbitalInfo - WARNING: event at time outside the run time window, skipping it\n");          if (verbose) printf(" OrbitalInfo - WARNING: event at time outside the run time window, skipping it\n");
628          debug = true;          jumped++;
629    //      debug = true;
630          continue;          continue;
631        }        }
632    
633          //
634          // retrieve tof informations
635          //
636          if ( !reprocall ){
637            itr = nobefrun + (re - evfrom - jumped);
638            //itr = re-(46438+200241);
639          } else {
640            itr = runinfo->GetFirstEntry() + (re - evfrom - jumped);
641          };
642          //
643          if ( !standalone ){
644            if ( itr > nevtofl2 ){  
645              if ( verbose ) printf(" OrbitalInfo - ERROR: no tof events with entry = %i in Level2 file\n",itr);
646              if ( debug ) printf(" nobefrun %u re %u evfrom %u jumped %u reprocall %i \n",nobefrun,re,evfrom,jumped,reprocall);
647              l0File->Close();
648              code = -901;
649              goto closeandexit;
650            };
651            //
652            tof->Clear();
653            //
654            ttof->GetEntry(itr);
655            //
656          };
657        //        //
658        procev++;        procev++;
659        //        //
# Line 576  int OrbitalInfoCore(UInt_t run, TFile *f Line 661  int OrbitalInfoCore(UInt_t run, TFile *f
661        //        //
662        orbitalinfo->Clear();        orbitalinfo->Clear();
663        //        //
664          OrbitalInfoTrkVar *t_orb = new OrbitalInfoTrkVar();
665          if( !(orbitalinfo->OrbitalInfoTrk) ) orbitalinfo->OrbitalInfoTrk = new TClonesArray("OrbitalInfoTrkVar",2);
666          TClonesArray &tor = *orbitalinfo->OrbitalInfoTrk;
667          //
668        // Fill OBT, pkt_num and absTime        // Fill OBT, pkt_num and absTime
669        //              //      
       ph = eh->GetPscuHeader();  
670        orbitalinfo->pkt_num = ph->GetCounter();        orbitalinfo->pkt_num = ph->GetCounter();
671        orbitalinfo->OBT = ph->GetOrbitalTime();        orbitalinfo->OBT = ph->GetOrbitalTime();
672        orbitalinfo->absTime = atime;        orbitalinfo->absTime = atime;
# Line 612  int OrbitalInfoCore(UInt_t run, TFile *f Line 700  int OrbitalInfoCore(UInt_t run, TFile *f
700        //        //
701        cOrbit orbits(*gltle->GetTle());        cOrbit orbits(*gltle->GetTle());
702        //        //
703          if ( debug ) printf(" I am Here \n");
704          //
705        // synchronize with quaternions data        // synchronize with quaternions data
706        //        //
707        if ( isf ){        if ( isf && neventsm>0 ){
708            if ( debug ) printf(" I am here \n");
709          //          //
710          // First event          // First event
711          //          //
# Line 626  int OrbitalInfoCore(UInt_t run, TFile *f Line 717  int OrbitalInfoCore(UInt_t run, TFile *f
717            tmpSize = mcmdev->Records->GetEntries();            tmpSize = mcmdev->Records->GetEntries();
718            numrec = tmpSize;            numrec = tmpSize;
719            for (Int_t j3 = 0;j3<tmpSize;j3++){            for (Int_t j3 = 0;j3<tmpSize;j3++){
720                if ( debug ) printf(" eh eh eh \n");
721              mcmdrc = (pamela::McmdRecord*)mcmdev->Records->At(j3);              mcmdrc = (pamela::McmdRecord*)mcmdev->Records->At(j3);
722              if ((int)mcmdrc->ID1 == 226){              if ((int)mcmdrc->ID1 == 226){
723                L_QQ_Q_l_upper->fill(mcmdrc->McmdData);                L_QQ_Q_l_upper->fill(mcmdrc->McmdData);
# Line 689  int OrbitalInfoCore(UInt_t run, TFile *f Line 781  int OrbitalInfoCore(UInt_t run, TFile *f
781        };        };
782      closethisloop:      closethisloop:
783        //        //
784        if ((atime>(UInt_t)upperqtime)||(atime<(UInt_t)lowerqtime)){        if ( debug ) printf(" I am There \n");
785          //
786          if (((atime>(UInt_t)upperqtime)||(atime<(UInt_t)lowerqtime)) && neventsm>0 ){
787            if ( debug ) printf(" I am there \n");
788          //          //
789          lowerqtime = upperqtime;          lowerqtime = upperqtime;
790          UInt_t maxloop = 100000000;              Long64_t maxloop = 100000000LL;
791          UInt_t mn = 0;          Long64_t mn = 0;
792          bool gh=false;          bool gh=false;
793          ooi=oi;          ooi=oi;
794          if ( verbose ) printf(" OrbitalInfoCore: sync with quaternions data upperqtime %u lowerqtime %u atime %u \n",(UInt_t)upperqtime,(UInt_t)lowerqtime,atime);          if ( verbose ) printf(" OrbitalInfoCore: sync with quaternions data upperqtime %u lowerqtime %u atime %u \n",(UInt_t)upperqtime,(UInt_t)lowerqtime,atime);
# Line 701  int OrbitalInfoCore(UInt_t run, TFile *f Line 796  int OrbitalInfoCore(UInt_t run, TFile *f
796            if ( mn > maxloop ){            if ( mn > maxloop ){
797              if ( verbose ) printf(" OrbitalInfoCore: quaternions sync out of range! exiting\n");              if ( verbose ) printf(" OrbitalInfoCore: quaternions sync out of range! exiting\n");
798              gh = true;              gh = true;
799                neventsm = 0;
800            };            };
801            mn++;            mn++;
802            if (oi<5) oi++;            if (oi<5) oi++;
803            else oi=0;            else oi=0;
804            if (oi==0){            if (oi==0 && numrec > 0){
805                if ( debug ) printf(" mumble \n");
806              mcreen++;              mcreen++;
807              if (mcreen == numrec){              if (mcreen == numrec){
808                mctren++;                mctren++;
# Line 745  int OrbitalInfoCore(UInt_t run, TFile *f Line 842  int OrbitalInfoCore(UInt_t run, TFile *f
842          if ( verbose ) printf(" OrbitalInfoCore: sync with quaternions data now we have upperqtime %u lowerqtime %u atime %u \n",(UInt_t)upperqtime,(UInt_t)lowerqtime,atime);          if ( verbose ) printf(" OrbitalInfoCore: sync with quaternions data now we have upperqtime %u lowerqtime %u atime %u \n",(UInt_t)upperqtime,(UInt_t)lowerqtime,atime);
843        };        };
844        //        //
845          if ( debug ) printf(" I am THIS \n");
846          //
847        // Fill in quaternions and angles        // Fill in quaternions and angles
848        //        //
849        if ((atime<=(UInt_t)upperqtime)&&(atime>=(UInt_t)lowerqtime)){            if ((atime<=(UInt_t)upperqtime)&&(atime>=(UInt_t)lowerqtime)&& neventsm>0){      
850            if ( debug ) printf(" I am this \n");
851          UInt_t tut = holeq(lowerqtime, upperqtime, L_QQ_Q_l_lower, L_QQ_Q_l_upper, oi);          UInt_t tut = holeq(lowerqtime, upperqtime, L_QQ_Q_l_lower, L_QQ_Q_l_upper, oi);
852          if (oi == 0){          if (oi == 0){
853            if ((tut!=5)||(tut!=6)){            if ((tut!=5)||(tut!=6)){
# Line 830  int OrbitalInfoCore(UInt_t run, TFile *f Line 930  int OrbitalInfoCore(UInt_t run, TFile *f
930              }              }
931            }                        }            
932          }          }
933            //
934          orbitalinfo->mode = holeq(lowerqtime, upperqtime, L_QQ_Q_l_lower, L_QQ_Q_l_upper, oi);          orbitalinfo->mode = holeq(lowerqtime, upperqtime, L_QQ_Q_l_lower, L_QQ_Q_l_upper, oi);
935                    //
936        } else {        } else {
937          orbitalinfo->mode = -1;          if ( debug ) printf(" ops no incl! \n");
938            orbitalinfo->mode = 10;
939        };        };
940            //
941          // ops no inclination information
942          //
943          if ( orbitalinfo->q0< -999 || orbitalinfo->q1 < -999 || orbitalinfo->q2 < -999 || orbitalinfo->q3 < -999 || orbitalinfo->q0 != orbitalinfo->q0 || orbitalinfo->q1 != orbitalinfo->q1 || orbitalinfo->q2 != orbitalinfo->q2 || orbitalinfo->q3 != orbitalinfo->q3 ){
944            orbitalinfo->mode = 10;
945            orbitalinfo->q0 = -1000.;
946            orbitalinfo->q1 = -1000.;
947            orbitalinfo->q2 = -1000.;
948            orbitalinfo->q3 = -1000.;
949            orbitalinfo->etha = -1000.;
950            orbitalinfo->phi = -1000.;
951            orbitalinfo->theta = -1000.;
952          };
953          //
954          // #########################################################################################################################  
955        //        //
956        // fill orbital positions        // fill orbital positions
957        //                //        
# Line 865  int OrbitalInfoCore(UInt_t run, TFile *f Line 981  int OrbitalInfoCore(UInt_t run, TFile *f
981          orbitalinfo->BB0 = babs/bequ;          orbitalinfo->BB0 = babs/bequ;
982          orbitalinfo->L = xl;                orbitalinfo->L = xl;      
983          // Set Stormer vertical cutoff using L shell.          // Set Stormer vertical cutoff using L shell.
984          orbitalinfo->cutoff[0] = 14.9/(xl*xl);          orbitalinfo->cutoffsvl = 14.9/(xl*xl);
985          //          //
986        };              };      
987        //        //
988          // pitch angles
989          //
990          if ( orbitalinfo->mode != 10 && orbitalinfo->mode != 5 && orbitalinfo->mode !=7 && orbitalinfo->mode != 9 ){
991            //
992            Float_t Bx = -orbitalinfo->Bdown;                       //don't need for PamExp ExpOnly for all geography areas
993            Float_t By = orbitalinfo->Beast;                        //don't need for PamExp ExpOnly for all geography areas
994            Float_t Bz = orbitalinfo->Bnorth;                       //don't need for PamExp ExpOnly for all geography areas
995            //
996            TMatrixD Fij = PO->ECItoGreenwich(PO->QuatoECI(orbitalinfo->q0,orbitalinfo->q1,orbitalinfo->q2,orbitalinfo->q3),orbitalinfo->absTime);
997            TMatrixD Dij = PO->GreenwichtoGEO(orbitalinfo->lat,orbitalinfo->lon,Fij);
998            TMatrixD Iij = PO->ColPermutation(Dij);
999            //
1000            orbitalinfo->Iij.ResizeTo(Iij);
1001            orbitalinfo->Iij = Iij;
1002            //
1003            A1 = Iij(0,2);
1004            A2 = Iij(1,2);
1005            A3 = Iij(2,2);
1006            //      
1007            //      orbitalinfo->pamzenitangle = (Float_t)PO->GetPitchAngle(1,0,0,A1,A2,A3);                        // Angle between zenit and Pamela's main axiz
1008            //      orbitalinfo->pamBangle = (Float_t)PO->GetPitchAngle(A1,A2,A3,Bx,By,Bz);                 // Angle between Pamela's main axiz and B
1009            //
1010            if ( !standalone && tof->ntrk() > 0 ){
1011              //
1012              Int_t nn = 0;
1013              for(Int_t nt=0; nt < tof->ntrk(); nt++){  
1014                //
1015                ToFTrkVar *ptt = tof->GetToFTrkVar(nt);
1016                Double_t E11x = ptt->xtr_tof[0]; // tr->x[0];
1017                Double_t E11y = ptt->ytr_tof[0]; //tr->y[0];
1018                Double_t E11z = zin[0];
1019                Double_t E22x = ptt->xtr_tof[3];//tr->x[3];
1020                Double_t E22y = ptt->ytr_tof[3];//tr->y[3];
1021                Double_t E22z = zin[3];
1022                if ( (E11x < 100. && E11y < 100. && E22x < 100. && E22y < 100.) || ptt->trkseqno != -1  ){
1023                  Double_t norm = sqrt(pow(E22x-E11x,2)+pow(E22y-E11y,2)+pow(E22z-E11z,2));
1024                  Double_t MyAzim = TMath::RadToDeg()*atan(TMath::Abs(E22y-E11y)/TMath::Abs(E22x-E11x));
1025                  if(E22x-E11x>=0 && E22y-E11y <0) MyAzim =  360. - MyAzim;
1026                  if(E22x-E11x>=0 && E22y-E11y >=0) MyAzim = MyAzim;
1027                  if(E22x-E11x<0 && E22y-E11y >0) MyAzim = 180. - MyAzim;
1028                  if(E22x-E11x<0 && E22y-E11y <0) MyAzim = 180. + MyAzim;
1029                  Px = (E22x-E11x)/norm;
1030                  Py = (E22y-E11y)/norm;
1031                  Pz = (E22z-E11z)/norm;
1032                  //
1033                  t_orb->trkseqno = ptt->trkseqno;
1034                  //
1035                  TMatrixD Eij = PO->PamelatoGEO(Iij,Px,Py,Pz);
1036                  t_orb->Eij.ResizeTo(Eij);
1037                  t_orb->Eij = Eij;
1038                  //
1039                  TMatrixD Sij = PO->PamelatoGEO(Fij,Px,Py,Pz);
1040                  t_orb->Sij.ResizeTo(Sij);
1041                  t_orb->Sij = Sij;
1042                  //            
1043                  t_orb->pitch = (Float_t)PO->GetPitchAngle(Eij(0,0),Eij(1,0),Eij(2,0),Bx,By,Bz);
1044                  //
1045                  //
1046                  Double_t omega = PO->GetPitchAngle(Eij(0,0),Eij(1,0),Eij(2,0),cos(orbitalinfo->lon+TMath::Pi()/2)-sin(orbitalinfo->lon+TMath::Pi()/2),cos(orbitalinfo->lon+TMath::Pi()/2)+sin(orbitalinfo->lon+TMath::Pi()/2),1);
1047                  //
1048                  t_orb->cutoff = 59.3/(pow(orbitalinfo->L,2)*pow((1+sqrt(1-pow(orbitalinfo->L,-3/2)*cos(omega))),2));
1049                  //
1050                  if ( debug ) printf(" orbitalinfo->cutoffsvl %f vitaly %f \n",orbitalinfo->cutoffsvl,t_orb->cutoff);
1051                  //
1052                  new(tor[nn]) OrbitalInfoTrkVar(*t_orb);
1053                  nn++;
1054                  //
1055                  t_orb->Clear();
1056                  //
1057                };
1058                //
1059              };
1060            } else {
1061              if ( debug ) printf(" mmm... mode %u standalone %i ntrk %i \n",orbitalinfo->mode,standalone,tof->ntrk());
1062            };
1063            //
1064          } else {
1065            if ( !standalone && tof->ntrk() > 0 ){
1066              //
1067              Int_t nn = 0;
1068              for(Int_t nt=0; nt < tof->ntrk(); nt++){  
1069                //
1070                ToFTrkVar *ptt = tof->GetToFTrkVar(nt);
1071                if ( ptt->trkseqno != -1  ){
1072                  //
1073                  t_orb->trkseqno = ptt->trkseqno;
1074                  //
1075                  t_orb->Eij = 0;  
1076                  //
1077                  t_orb->Sij = 0;
1078                  //            
1079                  t_orb->pitch = -1000.;
1080                  //
1081                  t_orb->cutoff = -1000.;
1082                  //
1083                  new(tor[nn]) OrbitalInfoTrkVar(*t_orb);
1084                  nn++;
1085                  //
1086                  t_orb->Clear();
1087                  //
1088                };
1089                //
1090              };    
1091            };
1092          };
1093          //
1094        // Fill the class        // Fill the class
1095        //        //
1096        OrbitalInfotr->Fill();        OrbitalInfotr->Fill();
1097        //        //
1098          delete t_orb;
1099          //
1100      }; // loop over the events in the run      }; // loop over the events in the run
1101      //      //
1102      // Here you may want to clear some variables before processing another run        // Here you may want to clear some variables before processing another run  
1103      //      //
1104      delete dbtime;      delete dbtime;
1105      delete L_QQ_Q_l_upper;      if ( L_QQ_Q_l_upper ) delete L_QQ_Q_l_upper;
1106      delete L_QQ_Q_l_lower;      if ( L_QQ_Q_l_lower ) delete L_QQ_Q_l_lower;
1107      delete RYPang_upper;      if ( RYPang_upper ) delete RYPang_upper;
1108      delete RYPang_lower;      if ( RYPang_lower ) delete RYPang_lower;
1109    }; // process all the runs    }; // process all the runs
1110        
1111    if (verbose) printf("\n Finished processing data \n");    if (verbose) printf("\n Finished processing data \n");
# Line 921  int OrbitalInfoCore(UInt_t run, TFile *f Line 1145  int OrbitalInfoCore(UInt_t run, TFile *f
1145    //    //
1146    if ( l0File ) l0File->Close();    if ( l0File ) l0File->Close();
1147    if ( tempfile ) tempfile->Close();                if ( tempfile ) tempfile->Close();            
1148    gSystem->Unlink(tempname.str().c_str());    if ( myfold ) gSystem->Unlink(tempname.str().c_str());
1149    //    //
1150    if ( runinfo ) runinfo->Close();        if ( runinfo ) runinfo->Close();    
1151    if ( OrbitalInfotr ) OrbitalInfotr->SetName("OrbitalInfo");        if ( OrbitalInfotr ) OrbitalInfotr->SetName("OrbitalInfo");    
1152      if ( tof ) tof->Delete();
1153      if ( ttof ) ttof->Delete();
1154      //
1155    if ( file ){    if ( file ){
1156      file->cd();      file->cd();
1157      file->Write();      file->Write();
1158    };    };
1159    //    //
1160    gSystem->Unlink(OrbitalInfofolder.str().c_str());    if ( myfold ) gSystem->Unlink(OrbitalInfofolder.str().c_str());
1161    //    //
1162    // the end    // the end
1163    //    //
1164      if ( dbc ){
1165        dbc->Close();
1166        delete dbc;
1167      };
1168    if (verbose) printf("\n Exiting...\n");    if (verbose) printf("\n Exiting...\n");
1169    if(OrbitalInfotr)OrbitalInfotr->Delete();    if(OrbitalInfotr)OrbitalInfotr->Delete();
1170    //    //
1171      if ( PO ) delete PO;
1172    if ( orbitalinfo ) delete orbitalinfo;    if ( orbitalinfo ) delete orbitalinfo;
1173    if ( orbitalinfoclone ) delete orbitalinfoclone;    if ( orbitalinfoclone ) delete orbitalinfoclone;
1174    if ( glroot ) delete glroot;    if ( glroot ) delete glroot;
Legend:
Removed from v.1.17  
changed lines
  Added in v.1.35
  ViewVC Help
Powered by ViewVC 1.1.23