/[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.25 by mocchiut, Sun Sep 9 18:57:29 2007 UTC revision 1.35 by mocchiut, Thu Dec 11 10:08:19 2008 UTC
# Line 58  int OrbitalInfoCore(UInt_t run, TFile *f Line 58  int OrbitalInfoCore(UInt_t run, TFile *f
58    TString psw = glt->CGetPsw();    TString psw = glt->CGetPsw();
59    TSQLServer *dbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data());    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    //    //
68    // Set these to true to have a very verbose output.    // Set these to true to have a very verbose output.
# Line 65  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 83  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 220  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 240  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 363  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 457  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 464  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 490  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 572  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            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 584  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 637  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 720  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 853  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 885  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  
# Line 945  int OrbitalInfoCore(UInt_t run, TFile *f Line 1149  int OrbitalInfoCore(UInt_t run, TFile *f
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();
# Line 961  int OrbitalInfoCore(UInt_t run, TFile *f Line 1168  int OrbitalInfoCore(UInt_t run, TFile *f
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.25  
changed lines
  Added in v.1.35
  ViewVC Help
Powered by ViewVC 1.1.23