/**
* \file TrkCore.cpp
* \author Elena Vannuccini
*/
// .........................................................
// ROOT header files
// .........................................................
#include <TTree.h>
#include <TClassEdit.h>
#include <TObject.h>
#include <TList.h>
#include <TSystem.h>
#include <TSystemDirectory.h>
#include <TString.h>
#include <TFile.h>
#include <TClass.h>
#include <TCanvas.h>
#include <TH1.h>
#include <TH1F.h>
#include <TH2D.h>
#include <TLatex.h>
#include <TPad.h>
#include <TSQLServer.h>
#include <TSQLRow.h>
#include <TSQLResult.h>
// .........................................................
// other general header files
// .........................................................
#include <fstream>
#include <iostream>
// .........................................................
// files in the inc directory
// .........................................................
#include <RunInfo.h>
#include <GLTables.h>
#include <TrkVerl2.h>
#include <TrkProcess.h>
#include <TrkStruct.h>
#include <TrkLevel2.h>
//#include <TrkLevel1.h>
#include <TrkLevel0.h>
// .........................................................
// YODA header files
// .........................................................
#include <PamelaRun.h>
#include <RegistryEvent.h>
#include <physics/tracker/TrackerEvent.h>
#include <CalibTrk1Event.h>
#include <CalibTrk2Event.h>
// .........................................................
// tracker analysis FORTRAN77 routines
// .........................................................
extern "C" {
extern struct cTrkCalib pedsigbad_;
extern struct cTrkLevel0 level0event_;
extern struct cTrkLevel1 level1event_;
extern struct cTrkLevel2 level2event_;
extern struct cPath path_;
extern struct cBPath bpath_;
int fillpedsigfromdefault_();
int readmipparam_();
int readchargeparam_();
int readvkmask_();
int readalignparam_();
int readetaparam_();
int reductionflight_();
int analysisflight_();
}
// .........................................................
// 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(ULong64_t run, TFile *f2, TSQLServer *dbc, int ncustom, char*vcustom[])
{
// ---------------------------
// Define some basic variables
// ---------------------------
TString Frameworklev1="paw";
TString Frameworklev2="root";
TString filename = 0;
TString trkcalibfile = 0;
Int_t last_trk_calib_used = 0;
Long64_t nevents = 0LL;
// LEVEL0 input
TFile *f0 = 0;
TTree *physicsTree = 0;
TBranch *b_registry = 0;
TBranch *b_trk = 0;
TBranch *b_header = 0;
RegistryEvent *reg = 0;
EventHeader *header = 0;
PscuHeader *pscu = 0;
TrkLevel0 *l0_event = 0;
// RunInfo
ItoRunInfo *runinfo = 0;
TArrayL *runlist = 0;
ULong64_t from_run;
ULong64_t to_run;
Bool_t reprocessing = false;
//LEVEL1 output - ROOT
// TFile *f1;
// TTree *t_level1;
// TrkLevel1 *l1_event;
//LEVEL2 output - ROOT
TTree *t_level2 = 0;
TTree *t_clone = 0;
TrkLevel2 *l2_event = new TrkLevel2();
TrkLevel2 *l2_clone = new TrkLevel2();
// -----------------------
// -----------------------
// 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);
p->HandleCustomPar(ncustom,vcustom);
path_.debug = p->DEBUG; // temporaneo.. vedro` di fare di meglio
if( p->DEBUG )p->Dump();
// p->Dump();
// ===================================
// Open/Create level1 output file
// (NB output files other than LEVEL2 are
// created in a separate folder)
// ===================================
if(p->get1){
if(p->ifroot1){
throw -299;
}else{
throw -299;
};
};
// ===================================
// Open/Create level2 output file
// ===================================
if(p->get2){
// if(p->ifroot2){ // root
//-------------------------------------------
// read RunInfo
//-------------------------------------------
if(!(p->standalone)){
// Open "RunInfo" tree and get list of run
runinfo = new ItoRunInfo(f2);
char *trkversion = TrkInfo(false);
if( runinfo->Update(run,"TRK",trkversion) ) throw -205;
runlist = runinfo->GetRunList();
reprocessing = runinfo->IsReprocessing();//????
if(p->DEBUG){
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 ");
t_level2->Branch("TrkLevel2","TrkLevel2",&l2_event);
// }else{ // hbook
// throw -299;
// }
};
// -------------------------------------------
// 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 "expiration date" of calibrations and parameters
//
ULong64_t pedsig_time =0ULL;
ULong64_t B_time =0ULL;
ULong64_t mip_time =0ULL;
ULong64_t charge_time =0ULL;
ULong64_t eta_time =0ULL;
ULong64_t align_time =0ULL;
ULong64_t mask_time =0ULL;
//
// 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; i<runinfo->GetFirstEntry(); i++){
t_clone->GetEntry(i);
*l2_event = *l2_clone;
t_level2->Fill();
l2_event->Clear();
// COPY COPY COPY
};
};
// ------------------------------------------------------------
// ------------------------------------------------------------
// START LOOP OVER RUNS TO PROCESSED/REPROCESSED
// ------------------------------------------------------------
// ------------------------------------------------------------
for(ULong64_t idRun=from_run; idRun <= to_run; idRun++){
if(p->DEBUG) cout << " ========================= Run: "<< idRun << endl;
ULong64_t runheadtime = 0ULL;
ULong64_t runtrailtime = 0ULL;
UInt_t evfrom = 0;
UInt_t evto = 0;
Int_t id_reg_run =-1;
Int_t trk_calib_used = 0;
if(p->standalone){
// ==============================
// first query: retrieve run info
// ==============================
if (q1.Query_GL_RUN(idRun,dbc) )throw -50;
id_reg_run = q1.ID_REG_RUN;
runheadtime = q1.RUNHEADER_TIME;
runtrailtime = q1.RUNTRAILER_TIME;
evfrom = q1.EV_REG_PHYS_FROM;
evto = q1.EV_REG_PHYS_TO;
trk_calib_used = q1.TRK_CALIB_USED;
}else{
// ==============================
// get run info from RunInfo tree
// ==============================
if( runinfo->GetRunInfo(idRun) ) throw -205;
id_reg_run = runinfo->ID_REG_RUN;
runheadtime = runinfo->RUNHEADER_TIME;
runtrailtime = runinfo->RUNTRAILER_TIME;
evfrom = runinfo->EV_REG_PHYS_FROM;
evto = runinfo->EV_REG_PHYS_TO;
trk_calib_used = runinfo->TRK_CALIB_USED;
};
//
if(p->DEBUG){
cout << "ROOT file ID "<< id_reg_run << endl;
cout << "RunHeader time "<< runheadtime << endl;
cout << "RunTrailer time "<< runtrailtime << endl;
cout << " from event "<< evfrom << endl;
cout << " to event "<< evto << 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_reg_run,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(p->DEBUG) cout << "Opening LEVEL0 file: "<< filename << endl;
f0 = new TFile(filename);
/* if ( !f0 ){
printf("\n\n ERROR: no data file, exiting...\n\n");
return(103);
};*/
if ( !f0 ) throw -6;
physicsTree = (TTree*)f0->Get("Physics"); if(!physicsTree) throw -7;
b_header = physicsTree ->GetBranch("Header"); if(!b_header ) throw -8;
b_registry = physicsTree ->GetBranch("Registry"); if(!b_registry ) throw -9;
b_trk = physicsTree ->GetBranch("Tracker"); if(!b_trk ) throw -203;
physicsTree->SetBranchAddress("Tracker" ,&l0_event);
physicsTree->SetBranchAddress("Registry",®);
physicsTree->SetBranchAddress("Header",&header);
nevents = physicsTree->GetEntries();
if ( nevents < 1 ) throw -11;
};
// =============================================================
// retrieve calibration file needed to process the run
// =============================================================
if(trk_calib_used==104 && last_trk_calib_used!=104){
if(p->DEBUG )cout << "** Using default calibration **" << endl;
if (q4.Query_GL_PARAM(runheadtime,"default calibration",dbc) )throw -52;
path_.FillWith(q4.PATH+q4.NAME);
fillpedsigfromdefault_();
if(path_.error) throw -216;
}else if( (trk_calib_used==1 | trk_calib_used==2) && runheadtime > pedsig_time){
if( q2.Query_GL_TRK_CALIB(runheadtime,dbc) )throw -53;
pedsig_time = q2.TO_TIME;
if(q2.EV_REG_CALIBTRK1 != q2.EV_REG_CALIBTRK2)
printf("WARNING!! ---> EV_REG_CALIBTRK1=%d it's different from EV_REG_CALIBTRK2=%d \n\n",q2.EV_REG_CALIBTRK1,q2.EV_REG_CALIBTRK2);
if( q3.Query_GL_ROOT(q2.ID_REG_CALIBTRK,dbc) )throw -51;
trkcalibfile = q3.PATH + q3.NAME;
// store calibration informations
if(p->DEBUG) cout << "Online calibration: "<< trkcalibfile << endl;
// pedsigbad_ = TrkFlightPede(trkcalibfile,q2.EV_REG_CALIBTRK1,q2.EV_REG_CALIBTRK2);
TFile* f0_c= f0;
if(trkcalibfile.CompareTo(filename))f0_c = new TFile(trkcalibfile);
pedsigbad_.FillFrom(f0_c,q2.EV_REG_CALIBTRK1,q2.EV_REG_CALIBTRK2);
if(trkcalibfile.CompareTo(filename))f0_c->Close();
};
last_trk_calib_used = trk_calib_used;
// =============================================================
// retrieve information about parameters to process LEVEL2
// =============================================================
//
// magnetic field
//
if(p->DEBUG) cout << "Loading parameter files : "<< endl;
if(runheadtime > B_time){
if( q4.Query_GL_PARAM(runheadtime,"field",dbc) )throw -52;
B_time = q4.TO_TIME;
l2_event->LoadField(q4.PATH+q4.NAME);
l2_event->LoadField(q4.PATH+q4.NAME);
if(path_.error) throw -215;
};
//
// mip conversion parameters
//
if(runheadtime > mip_time){
if( q4.Query_GL_PARAM(runheadtime,"trk mip",dbc) )throw -52;
mip_time = q4.TO_TIME;
path_.FillWith(q4.PATH+q4.NAME);
readmipparam_();
if(path_.error) throw -212;
};
//
// charge correlation parameters
//
if(runheadtime > charge_time){
if( q4.Query_GL_PARAM(runheadtime,"trk charge",dbc) )throw -52;
charge_time = q4.TO_TIME;
path_.FillWith(q4.PATH+q4.NAME);
readchargeparam_();
if(path_.error) throw -213;
};
//
// eta p.f.a. parameters
//
if(runheadtime > eta_time){
if( q4.Query_GL_PARAM(runheadtime,"trk pfa",dbc) )throw -52;
eta_time = q4.TO_TIME;
path_.FillWith(q4.PATH+q4.NAME);
readetaparam_();
if(path_.error) throw -214;
};
//
// alignment parameters
//
if(runheadtime > align_time){
if( q4.Query_GL_PARAM(runheadtime,"trk alignment",dbc) )throw -52;
align_time = q4.TO_TIME;
path_.FillWith(q4.PATH+q4.NAME);
readalignparam_();
if(path_.error) throw -211;
};
//
// viking mask
//
if(runheadtime > mask_time){
if( q4.Query_GL_PARAM(runheadtime,"trk mask",dbc) )throw -52;
mask_time = q4.TO_TIME;
path_.FillWith(q4.PATH+q4.NAME);
readvkmask_();
if(path_.error) throw -210;
};
//
//
if( p->DEBUG ){
cout<<"Total number of events : "<<nevents<<endl;
cout<<"Events to analize for run "<<run<<" : "<<evto+1-evfrom;
cout<<" (from "<<evfrom<<" to "<<evto<<")"<<endl;
};
// ====================================================
// start looping on events cointained in the data file
// ====================================================
for (UInt_t re = evfrom; re <= evto; re++){
ev_count++;
// if ( p->DEBUG && re%100 == 0 && re > 0 ) printf(" %iK \n",re/1000);
if ( p->DEBUG && re%100 == 0 && re > 0 ) cout << ".";
b_registry->GetEntry(re);
b_trk->GetEntry(reg->event);
b_header->GetEntry(reg->event);
pscu = header->GetPscuHeader();
// pscu->Print();
// cout <<endl<< " *************** Pkt ID "<< (UINT32)pscu->Counter << endl;
l0_event->GetCommonVar(&level0event_);
// =============================
if(p->get1|p->get2) reductionflight_();
if(p->get2) analysisflight_();
// =============================
// ----------------
// LEVEL1 output
// ----------------
if(p->get1){
if(p->ifroot1){ // root
}else{ // hbook
};
};
// ----------------
// LEVEL2 output
// ----------------
if(p->get2){
// if(p->ifroot2){ // root
l2_event->FillCommonVar(&level2event_);
// l2_event->Dump();
t_level2->Fill();
l2_event->Clear();
// }else{ // hbook
// }
};
}; // end loop on events
if(p->DEBUG)cout << " >>> processed "<< ev_count <<" events"<< endl;
}; // 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; i<runinfo->GetFileEntries(); i++){
t_clone->GetEntry(i);
*l2_event = *l2_clone;
t_level2->Fill();
l2_event->Clear();
// COPY COPY COPY
};
};
// ---------------
// close the files
// ---------------
if(p->get1){
if(p->ifroot1){ // root
}else{ // hbook
}
};
if(p->get2){
// if(p->ifroot2){ // root
if( t_clone )t_clone->Delete("all");//delete old tree from file
if( !(p->standalone) )runinfo->Close();
f2->Write("Tracker");
if( t_level2 )t_level2->Delete(); //delete new tree from memory
// }else{ // hbook
// }
};
f0->Close();
return(0);
}