/** * \file TrkCore.cpp * \author Elena Vannuccini */ // ......................................................... // ROOT header files // ......................................................... #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include // ......................................................... // other general header files // ......................................................... #include #include // ......................................................... // files in the inc directory // ......................................................... #include #include #include #include //#include #include #include #include #include // ......................................................... // YODA header files // ......................................................... #include //#include #include #include #include // ......................................................... // namespaces // ......................................................... using namespace std; using namespace pamela; // ================================================================================ // // // ================================================================================ /** * \brief Tracker data reduction routine. * * It performs data reduction LEVEL0->LEVEL1->LEVEL2, producing LEVEL1 and or LEVEL2 output, in ROOT or HBOOK format. * Input parameters: * @param run id of the run to be processed (if run=0 a whole file is reprocessed) * @param dbc pointer to BD server * @param file1 output LEVEL1 file name (null if no LEVEL1 output is required) * @param file2 output LEVEL2 file name (null if no LEVEL2 output is required) * @param standalone (flag to run the program in standalone mode, that is without reading RunInfo) */ //short int TrkCore(Int_t run, TSQLServer *dbc, TString file1, TString file2, Bool_t standalone) //int TrkCore(UInt_t run, TFile *f2, TSQLServer *dbc, int ncustom, char*vcustom[]) // EMILIANO int TrkCore(UInt_t run, TFile *f2, GL_TABLES *glt, int ncustom, char*vcustom[]) // EMILIANO { // --------------------------- // Define some basic variables // --------------------------- TString filename = ""; Long64_t nentries = 0LL; // LEVEL0 input TFile *f0 = 0; TTree *physicsTree = 0; TBranch *b_trk = 0; TBranch *b_header = 0; EventHeader *header = 0; PscuHeader *pscu = 0; TrkLevel0 *l0_event = new TrkLevel0(); // RunInfo ItoRunInfo *runinfo = 0; TArrayI *runlist = 0; UInt_t from_run; UInt_t to_run; Bool_t reprocessing = false; //LEVEL2 output - ROOT TTree *t_level2 = 0; TTree *t_clone = 0; TrkLevel2 *l2_event = new TrkLevel2(); TrkLevel2 *l2_clone = new TrkLevel2(); //LEVEL1 output - ROOT TrkLevel1 *l1_event = new TrkLevel1(); TrkHough *lh_event = new TrkHough(); // ----------------------- // ----------------------- // Handle input parameters // (data) // // - open/create output files, determining the environment root/hbook from the estension // - create the level1/level2 tree+branch/nt-uple // ----------------------- // ----------------------- TrkProcess *p = new TrkProcess(run,f2); if( p->HandleCustomPar(ncustom,vcustom) )return(p->ostatus);; if( TrkParams::VerboseMode() )p->Dump(); // =================================== // Open/Create level2 output file // =================================== if(p->get2){ //------------------------------------------- // read RunInfo //------------------------------------------- if(!(p->standalone)){ // Open "RunInfo" tree and get list of run runinfo = new ItoRunInfo(f2); char *trkversion = TrkInfo(false); int runinfo_error = runinfo->Update(run,"TRK",trkversion); if( runinfo_error ) throw runinfo_error; runlist = runinfo->GetRunList(); reprocessing = runinfo->IsReprocessing();//???? if(TrkParams::VerboseMode()){ cout << "#events "<< runinfo->GetFileEntries() << endl;// #events in the file cout << "#runs "<< runinfo->GetRunEntries() << endl;// #runs in the file cout << "1-st run "<< runlist->At(0) << endl; cout << "last run "<< runlist->At(runinfo->GetRunEntries()-1) << endl; cout << "1-st event "<< runinfo->GetFirstEntry() << endl;// first event of our run cout << "last event+1 "<< runinfo->GetLastEntry() << endl;// first event AFTER the last event of our run cout << "reprocessing "<< runinfo->IsReprocessing() << endl << endl; }; }; //------------------------------------------- // //------------------------------------------- // Take (if present) the old Tracker tree t_clone = (TTree*)f2->Get("Tracker"); if( t_clone != NULL ) t_clone->SetBranchAddress("TrkLevel2",&l2_clone); // Create NEW Tracker tree t_level2 = new TTree("Tracker","PAMELA tracker level2 data "); l2_event->Set(); // ****NBNBNBN***** t_level2->Branch("TrkLevel2","TrkLevel2",&l2_event); if(p->get1){ if(TrkParams::VerboseMode())cout << endl << "Requested LEVEL1 output" << endl; l1_event->Set(); // ****NBNBNBN***** t_level2->Branch("TrkLevel1","TrkLevel1",&l1_event); } if(p->geth){ if(TrkParams::VerboseMode())cout << endl << "Requested Hough-Transform output" << endl; t_level2->Branch("TrkHough","TrkHough",&lh_event); }; }; // ------------------------------------------- // define runs to be processed/reprocessed // ------------------------------------------- if(run == 0){ // reprocessing ALL runs if(p->standalone)throw -298; // reprocessing not implemented from_run = runlist->At(0); to_run = runlist->At(runinfo->GetRunEntries()-1); }else{ // processing/reprocessing ONE single run from_run = run; to_run = run; }; // // init event counter // Int_t ev_count =0; // // create query-results objects // GL_RUN q1 = GL_RUN(); GL_TRK_CALIB q2 = GL_TRK_CALIB(); GL_ROOT q3 = GL_ROOT(); GL_PARAM q4 = GL_PARAM(); // ------------------------------------------------------------ // if reprocessing one run, copy all the events BEFORE the run // ------------------------------------------------------------ if( !(p->standalone) ){ for(UInt_t i=0; iGetFirstEntry(); i++){ if ( t_clone->GetEntry(i) <= 0 ) throw -36;// EM *l2_event = *l2_clone; t_level2->Fill(); l2_event->Clear(); // COPY COPY COPY }; }; // ------------------------------------------------------------ // ------------------------------------------------------------ // START LOOP OVER RUNS TO PROCESS/REPROCESS // ------------------------------------------------------------ // ------------------------------------------------------------ for(UInt_t idRun = from_run; idRun <= to_run; idRun++){ if(TrkParams::VerboseMode()) cout << endl<<" ========================= Run: "<< idRun << endl; UInt_t runheadtime = 0; UInt_t runtrailtime = 0; UInt_t runheadobt = 0; UInt_t runtrailobt = 0; UInt_t evfrom = 0; UInt_t evto = 0; UInt_t nevents = 0; Int_t id_root_l0 =-1; Int_t trk_calib_used = 0; TString host = glt->CGetHost(); // EMILIANO TString user = glt->CGetUser(); // EMILIANO TString psw = glt->CGetPsw(); // EMILIANO // TString host = "mysql://localhost/pamelaprod"; // TString user = "anonymous"; // TString psw = ""; // const char *pamdbhost=gSystem->Getenv("PAM_DBHOST"); // const char *pamdbuser=gSystem->Getenv("PAM_DBUSER"); // const char *pamdbpsw=gSystem->Getenv("PAM_DBPSW"); // if ( !pamdbhost ) pamdbhost = ""; // if ( !pamdbuser ) pamdbuser = ""; // if ( !pamdbpsw ) pamdbpsw = ""; // if ( strcmp(pamdbhost,"") ) host = pamdbhost; // if ( strcmp(pamdbuser,"") ) user = pamdbuser; // if ( strcmp(pamdbpsw,"") ) psw = pamdbpsw; TSQLServer *dbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data()); stringstream myquery; // EMILIANO myquery.str(""); // EMILIANO myquery << "SET time_zone='+0:00'"; // EMILIANO delete dbc->Query(myquery.str().c_str()); // EMILIANO if(p->standalone){ // ============================== // first query: retrieve run info // ============================== if (q1.Query_GL_RUN(idRun,dbc) )throw -50; id_root_l0 = q1.ID_ROOT_L0; runheadtime = q1.RUNHEADER_TIME; runtrailtime = q1.RUNTRAILER_TIME; runheadobt = q1.RUNHEADER_OBT; runtrailobt = q1.RUNTRAILER_OBT; evfrom = q1.EV_FROM; evto = q1.EV_TO; nevents = q1.NEVENTS; trk_calib_used = q1.TRK_CALIB_USED; }else{ // ============================== // get run info from RunInfo tree // ============================== int runinfo_error = runinfo->GetRunInfo(idRun); if( runinfo_error ) throw runinfo_error; id_root_l0 = runinfo->ID_ROOT_L0; runheadtime = runinfo->RUNHEADER_TIME; runtrailtime = runinfo->RUNTRAILER_TIME; runheadobt = runinfo->RUNHEADER_OBT; runtrailobt = runinfo->RUNTRAILER_OBT; evfrom = runinfo->EV_FROM; evto = runinfo->EV_TO; nevents = runinfo->NEVENTS; trk_calib_used = runinfo->TRK_CALIB_USED; }; // if(TrkParams::VerboseMode()){ cout << "ROOT file ID "<< id_root_l0 << endl; cout << "RunHeader time "<< runheadtime << endl; cout << "RunTrailer time "<< runtrailtime << endl; cout << " from event "<< evfrom << endl; cout << " to event "<< evto << endl; cout << " num. events "<< nevents << endl; cout << "trk-calibration used "<< trk_calib_used << endl; }; // ======================================================== // second query: search the LEVEL0 file that contains idRun // ======================================================== TString lastfilename = filename; if( q3.Query_GL_ROOT(id_root_l0,dbc) )throw -51; filename = q3.PATH + q3.NAME; // ======================================================== // Open the input LEVEL0 data file // ======================================================== if(filename.CompareTo(lastfilename)){ if(!lastfilename.IsNull())f0->Close(); //if( debug ) cout << "Opening LEVEL0 file: "<< filename << endl; if(TrkParams::VerboseMode()) cout << "Opening LEVEL0 file: "<< filename << endl; FileStat_t t; if( gSystem->GetPathInfo(filename.Data(),t) )throw -6; if ( f0 ) f0->Close(); f0 = new TFile(filename); if ( !f0 ) throw -6; physicsTree = (TTree*)f0->Get("Physics"); if(!physicsTree) throw -7; b_header = physicsTree ->GetBranch("Header"); if(!b_header ) throw -8; b_trk = physicsTree ->GetBranch("Tracker"); if(!b_trk ) throw -203; l0_event->Set(); physicsTree->SetBranchAddress("Tracker" ,l0_event->GetPointerToTrackerEvent()); physicsTree->SetBranchAddress("Header",&header); nentries = physicsTree->GetEntries(); if ( nentries < 1 && nevents ) throw -11; // EMILIANO: go on if the file is empty, why not? in principle we should not have any event to process but if the case the next line will take care of exiting; exit only if the file is empty but nevents is not zero if ( nentries < (evto+1) && nevents > 0 ) throw -12; // EMILIANO: if NEVENTS = 0 and the file is empty everything is ok but due to a mistake in the // variable EV_TO type (UInt_t instead of Int_t) that we don't want to change to avoid changing a lot of code, evto becomes inf hence // the condition ( nentries < (evto+1) ) is satisfied and DV exit with error even if the error is only in the DB. }; GL_TIMESYNC *dbtime = new GL_TIMESYNC(id_root_l0,"ID",dbc); // ============================================================= // retrieve information about parameters to process LEVEL2 // ============================================================= TrkParams::Set(runinfo->GetGL_RUN(),dbc); for(int i=0; inpar; i++){ if(TrkParams::VerboseMode())cout<<" ((( force parameters from input path )))"<parpath[i],p->partype[i]); } TrkParams::Load(); if( !TrkParams::IsLoaded() )throw -52; // ============================================================= // retrieve calibration file needed to reduce data // ============================================================= TrkParams::SetCalib(runinfo->GetGL_RUN(),dbc); TrkParams::LoadCalib( ); if( !TrkParams::CalibIsLoaded() )throw -53; TBenchmark *reduction = new TBenchmark(); if(TrkParams::VerboseMode())reduction->Start("reduction"); Int_t ntrk=0; // ==================================================== // start looping on events cointained in the data file // ==================================================== if(dbc){ dbc->Close(); delete dbc;}; // for (UInt_t re = evfrom+min(p->nskip,nevents); re < evfrom+nevents; re++){ ev_count++; // if ( TrkParams::DebugMode() && re%100 == 0 && re > 0 ) cout << "."; if ( b_trk->GetEntry(re) <= 0 ) throw -36;//EM if ( b_header->GetEntry(re) <= 0 ) throw -36;//EM pscu = header->GetPscuHeader(); if( TrkParams::DebugMode() )cout << ">>> "<GetOrbitalTime()<DBabsTime(pscu->GetOrbitalTime()) > (runtrailtime+1) || dbtime->DBabsTime(pscu->GetOrbitalTime()) < (runheadtime-1)) { if (TrkParams::VerboseMode()){ printf(" TrkCore - WARNING: event outside the run time window, skipping it\n"); cout << " OBT "<GetOrbitalTime()<<" RUN "<DBabsTime(pscu->GetOrbitalTime())<<" RUN "<ProcessEvent(l0_event); // ---------------- // LEVEL1 output // ---------------- if(p->get1){ if(p->ifroot1){ // root l1_event->Clear(); // l1_event->SetFromLevel1Struct(&level1event_,p->full1); l1_event->SetFromLevel1Struct(p->full1); // t_level1->Fill(); }else{ // hbook throw -299; }; }; // ---------------- // HOUGH output // ---------------- if(p->geth){ if(p->ifrooth){ // root lh_event->Delete(); lh_event->SetFromHoughStruct(&houghevent_); }else{ // hbook throw -299; }; }; // ---------------- // LEVEL2 output // ---------------- if(p->get2){ l2_event->Clear(); if(p->get1){ l2_event->SetFromLevel2Struct(&level2event_,l1_event);//set references to level1 }else{ l2_event->SetFromLevel2Struct(&level2event_); } // l2_event->Dump(); t_level2->Fill(); if( l2_event->ntrk()>0 )ntrk++; if(TrkParams::VerboseMode())l2_event->Dump(); }; }; }; // end loop on events if(TrkParams::VerboseMode()){ cout << " >>> processed "<< ev_count <<" events"<< endl; if(ev_count)cout << ntrk << " events with at least one track ("<<(Int_t)(100*ntrk/ev_count)<<"%)\n"; reduction->Show("reduction"); } delete reduction; delete dbtime; }; // end loop on runs // ------------------------------------------------------------ // if reprocessing one run, copy all the events AFTER the run // ------------------------------------------------------------ if( !(p->standalone) ){ for(UInt_t i=runinfo->GetLastEntry()+1; iGetFileEntries(); i++){ if ( t_clone->GetEntry(i) <= 0 ) throw -36;//EM *l2_event = *l2_clone; t_level2->Fill(); l2_event->Clear(); // COPY COPY COPY }; }; // --------------- // close the files // --------------- if(p->get2){ if( t_clone )t_clone->Delete("all");//delete old tree from file if( !(p->standalone) )runinfo->Close(); // t_level2->Write("Tracker", TObject::kOverwrite); f2->Write("Tracker", TObject::kOverwrite); if( t_level2 )t_level2->Delete(); //delete new tree from memory }; if(f0) if(f0->IsOpen()) f0->Close(); // if( f0_c->IsOpen() )f0_c->Close(); // lh_event->Delete(); l1_event->Delete(); l2_event->Delete(); l2_clone->Delete(); return(p->ostatus); }