/[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.28 by mocchiut, Thu Sep 4 15:46:18 2008 UTC revision 1.30 by mocchiut, Wed Oct 1 15:25:44 2008 UTC
# Line 70  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 = true;
75      //
76    if ( OrbitalInfoargc > 0 ){    if ( OrbitalInfoargc > 0 ){
77      i = 0;      i = 0;
78      while ( i < OrbitalInfoargc ){      while ( i < OrbitalInfoargc ){
# Line 88  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 225  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 245  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 368  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 462  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 581  int OrbitalInfoCore(UInt_t run, TFile *f Line 625  int OrbitalInfoCore(UInt_t run, TFile *f
625        //        //
626        if ( (atime > runinfo->RUNTRAILER_TIME) || (atime < runinfo->RUNHEADER_TIME)  ) {        if ( (atime > runinfo->RUNTRAILER_TIME) || (atime < runinfo->RUNHEADER_TIME)  ) {
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            jumped++;
629  //      debug = true;  //      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 591  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 644  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 727  int OrbitalInfoCore(UInt_t run, TFile *f Line 801  int OrbitalInfoCore(UInt_t run, TFile *f
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 860  int OrbitalInfoCore(UInt_t run, TFile *f Line 935  int OrbitalInfoCore(UInt_t run, TFile *f
935          //          //
936        } else {        } else {
937          if ( debug ) printf(" ops no incl! \n");          if ( debug ) printf(" ops no incl! \n");
938          orbitalinfo->mode = -1;          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 896  int OrbitalInfoCore(UInt_t run, TFile *f Line 985  int OrbitalInfoCore(UInt_t run, TFile *f
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            A1 = Iij(0,2);
1001            A2 = Iij(1,2);
1002            A3 = Iij(2,2);
1003            //      
1004            orbitalinfo->pamzenitangle = (Float_t)PO->GetPitchAngle(1,0,0,A1,A2,A3);                        // Angle between zenit and Pamela's main axiz
1005            orbitalinfo->pamBangle = (Float_t)PO->GetPitchAngle(A1,A2,A3,Bx,By,Bz);                 // Angle between Pamela's main axiz and B
1006            //
1007            if ( !standalone && tof->ntrk() > 0 ){
1008              //
1009              Int_t nn = 0;
1010              for(Int_t nt=0; nt < tof->ntrk(); nt++){  
1011                //
1012                ToFTrkVar *ptt = tof->GetToFTrkVar(nt);
1013                Double_t E11x = ptt->xtr_tof[0]; // tr->x[0];
1014                Double_t E11y = ptt->ytr_tof[0]; //tr->y[0];
1015                Double_t E11z = zin[0];
1016                Double_t E22x = ptt->xtr_tof[3];//tr->x[3];
1017                Double_t E22y = ptt->ytr_tof[3];//tr->y[3];
1018                Double_t E22z = zin[3];
1019                if ( E11x < 100. && E11y < 100. && E22x < 100. && E22y < 100.  ){
1020                  Double_t norm = sqrt(pow(E22x-E11x,2)+pow(E22y-E11y,2)+pow(E22z-E11z,2));
1021                  Double_t MyAzim = TMath::RadToDeg()*atan(TMath::Abs(E22y-E11y)/TMath::Abs(E22x-E11x));
1022                  if(E22x-E11x>=0 && E22y-E11y <0) MyAzim =  360. - MyAzim;
1023                  if(E22x-E11x>=0 && E22y-E11y >=0) MyAzim = MyAzim;
1024                  if(E22x-E11x<0 && E22y-E11y >0) MyAzim = 180. - MyAzim;
1025                  if(E22x-E11x<0 && E22y-E11y <0) MyAzim = 180. + MyAzim;
1026                  Px = (E22x-E11x)/norm;
1027                  Py = (E22y-E11y)/norm;
1028                  Pz = (E22z-E11z)/norm;
1029                  //
1030                  TMatrixD Eij = PO->PamelatoGEO(Iij,Px,Py,Pz);            
1031                  //            
1032                  t_orb->trkseqno = ptt->trkseqno;
1033                  t_orb->pitch = (Float_t)PO->GetPitchAngle(Eij(0,0),Eij(1,0),Eij(2,0),Bx,By,Bz);
1034                  //
1035                  new(tor[nn]) OrbitalInfoTrkVar(*t_orb);
1036                  nn++;
1037                  //
1038                  t_orb->Clear();
1039                  //
1040                };
1041                //
1042              };
1043            } else {
1044              if ( debug ) printf(" mmm... mode %u standalone %i ntrk %i \n",orbitalinfo->mode,standalone,tof->ntrk());
1045            };
1046            //
1047          };
1048          //
1049        // Fill the class        // Fill the class
1050        //        //
1051        OrbitalInfotr->Fill();        OrbitalInfotr->Fill();
1052        //        //
1053          delete t_orb;
1054          //
1055      }; // loop over the events in the run      }; // loop over the events in the run
1056      //      //
1057      // Here you may want to clear some variables before processing another run        // Here you may want to clear some variables before processing another run  
# Line 952  int OrbitalInfoCore(UInt_t run, TFile *f Line 1104  int OrbitalInfoCore(UInt_t run, TFile *f
1104    //    //
1105    if ( runinfo ) runinfo->Close();        if ( runinfo ) runinfo->Close();    
1106    if ( OrbitalInfotr ) OrbitalInfotr->SetName("OrbitalInfo");        if ( OrbitalInfotr ) OrbitalInfotr->SetName("OrbitalInfo");    
1107      if ( tof ) tof->Delete();
1108      if ( ttof ) ttof->Delete();
1109      //
1110    if ( file ){    if ( file ){
1111      file->cd();      file->cd();
1112      file->Write();      file->Write();
# Line 968  int OrbitalInfoCore(UInt_t run, TFile *f Line 1123  int OrbitalInfoCore(UInt_t run, TFile *f
1123    if (verbose) printf("\n Exiting...\n");    if (verbose) printf("\n Exiting...\n");
1124    if(OrbitalInfotr)OrbitalInfotr->Delete();    if(OrbitalInfotr)OrbitalInfotr->Delete();
1125    //    //
1126      if ( PO ) delete PO;
1127    if ( orbitalinfo ) delete orbitalinfo;    if ( orbitalinfo ) delete orbitalinfo;
1128    if ( orbitalinfoclone ) delete orbitalinfoclone;    if ( orbitalinfoclone ) delete orbitalinfoclone;
1129    if ( glroot ) delete glroot;    if ( glroot ) delete glroot;
Legend:
Removed from v.1.28  
changed lines
  Added in v.1.30
  ViewVC Help
Powered by ViewVC 1.1.23