TrackerLevel2/src/TrkCore.cpp

/**
 * \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>
#include <TBenchmark.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>

// .........................................................
// 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[])   
{
    // ---------------------------
    // Define some basic variables
    // ---------------------------

    TString  filename     = 0; 
    TString  trkcalibfile = 0;
    Int_t    last_trk_calib_used = 0;
    Long64_t nentries = 0LL;
    
    // LEVEL0 input 
    TFile         *f0 = 0;
    TFile         *f0_c = 0;
    TTree         *physicsTree = 0;
    TBranch       *b_trk       = 0;
    TBranch       *b_header    = 0;
    EventHeader   *header  = 0;
    PscuHeader    *pscu    = 0;
    TrkLevel0     *l0_event = 0;
    // 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);
    p->HandleCustomPar(ncustom,vcustom);
    if( p->DebugMode() )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(p->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 ");
        t_level2->Branch("TrkLevel2","TrkLevel2",&l2_event);
        if(p->get1){
            if(p->VerboseMode())cout << endl << "Requested LEVEL1 output" << endl;
            t_level2->Branch("TrkLevel1","TrkLevel1",&l1_event);
            t_level2->BranchRef();
        }else if(p->geth){
            if(p->VerboseMode())cout << endl << "Requested Hough-Transform output" << endl;
//          t_level2->Branch("TrkHough","TrkHough",&lh_event);
//          t_level2->BranchRef();
        };
    };
        

// -------------------------------------------
// 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
    //
    UInt_t pedsig_time =0;
    UInt_t B_time      =0;
    UInt_t mip_time    =0;
    UInt_t charge_time =0;
    UInt_t eta_time    =0;
    UInt_t align_time  =0;
    UInt_t mask_time   =0;
    //
    // 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(UInt_t idRun = from_run; idRun <= to_run; idRun++){
        
        if(p->VerboseMode()) cout << endl<<" ========================= Run: "<< idRun << endl;
        UInt_t runheadtime       = 0;
        UInt_t runtrailtime      = 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;

        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;
            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;
            evfrom         = runinfo->EV_FROM;
            evto           = runinfo->EV_TO;
            nevents        = runinfo->NEVENTS;
            trk_calib_used = runinfo->TRK_CALIB_USED;
        };
        //
        if(p->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(p->VerboseMode()) cout << "Opening LEVEL0 file: "<< filename << endl;
            FileStat_t t;
            if( gSystem->GetPathInfo(filename.Data(),t) )throw -6;
            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_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",&reg);
            physicsTree->SetBranchAddress("Header",&header);
            
            nentries = physicsTree->GetEntries(); 
            if ( nentries < 1 ) throw -11;
            if ( nentries < (evto+1)) throw -12;
        
        };
        
        GL_TIMESYNC *dbtime = new GL_TIMESYNC(id_root_l0,"ID",dbc);     

        // =============================================================
        // retrieve calibration file needed to reduce data 
        // =============================================================
        // if run OBT is > last calibration "expiration date"
        //  - search for new calibration packet
        //  - load calibration parameters (full + truncated)
        if(p->VerboseMode() )cout << "Full pedestals for cluster finding:";
        if(runheadtime > pedsig_time){
            if( q2.Query_GL_TRK_CALIB(runheadtime,dbc) )throw -53;
            pedsig_time = q2.TO_TIME;
            if(q2.EV_ROOT_CALIBTRK1 != q2.EV_ROOT_CALIBTRK2) 
                printf("WARNING!! ---> EV_ROOT_CALIBTRK1=%d it's different from EV_ROOT_CALIBTRK2=%d \n\n",q2.EV_ROOT_CALIBTRK1,q2.EV_ROOT_CALIBTRK2);
            if( q3.Query_GL_ROOT(q2.ID_ROOT_L0,dbc) )throw -51;
                trkcalibfile = q3.PATH + q3.NAME;
                if(p->VerboseMode()) cout << " >> Loading from LEVEL0 file: "<< trkcalibfile << endl;
            if(trkcalibfile.CompareTo(filename)){
                        FileStat_t t;
                        if( gSystem->GetPathInfo(trkcalibfile.Data(),t) )throw -6;
                        f0_c=new TFile(trkcalibfile);
                        if ( !f0_c ) throw -6;
            }else{
                        f0_c = f0;
            };
            if(p->VerboseMode()) {
                        cout << " calibration entries "<< q2.EV_ROOT_CALIBTRK1 << " " << q2.EV_ROOT_CALIBTRK2; 
                        cout << " (from time "<< q2.FROM_TIME <<" to time "<< q2.TO_TIME <<")"<<endl;
            };
            pedsigbad_.FillACalibFrom(f0_c,q2.EV_ROOT_CALIBTRK1,q2.EV_ROOT_CALIBTRK2);
        }else{
                if(p->VerboseMode() )cout<< " already loaded. "<< endl;
        };
        if( p->DebugMode() ) for(int i=0; i<12; i++) cout << " DSP "<< i << " "<< pedsigbad_.pedestal[64][12][i] << endl;
        // =============================================================
        // retrieve calibration file needed to uncompress data
        // =============================================================
        // if the run was compressed using default calib
        // load truncated pedestals from default
        // otherwise reload them from on-line calibration
        if(p->VerboseMode() )cout << "Truncated pedestals for uncompression:";
        if(trk_calib_used==104 && last_trk_calib_used!=104){
            
            if(p->VerboseMode() )cout << " >> Loading default calibration " << endl; 
            if (q4.Query_GL_PARAM(runheadtime,7,dbc) )throw -52;            
            pedsigbad_.FillTCalibFrom(q4.PATH+q4.NAME);
            
        }else if( (trk_calib_used==1 | trk_calib_used==2) &&  last_trk_calib_used==104  ){

            if(p->VerboseMode()) cout << ">> Re-loading on-line calibration " << endl;
            pedsigbad_.FillTCalibFrom(f0_c,q2.EV_ROOT_CALIBTRK1,q2.EV_ROOT_CALIBTRK2);
        }else{
            if(p->VerboseMode()) cout << ">> Using on-line calibration " << endl;
        };
        if( p->DebugMode() ) for(int i=0; i<12; i++) cout << " DSP "<< i << " "<< pedsigbad_.pedestal_t[64][12][i] << endl;
        last_trk_calib_used = trk_calib_used;
        
        // =============================================================
        // retrieve information about parameters to process LEVEL2
        // =============================================================  
        //
        // magnetic field
        //
        if(p->VerboseMode()) cout << "Loading parameter files : "<< endl;
        
        if(runheadtime > B_time){
            if( q4.Query_GL_PARAM(runheadtime,1,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,2,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,3,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,4,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,5,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,6,dbc) )throw -52;
            mask_time = q4.TO_TIME; 
            path_.FillWith(q4.PATH+q4.NAME);
            readvkmask_();
            if(path_.error) throw -210;
        };

        TBenchmark *reduction = new TBenchmark();
        if(p->VerboseMode())reduction->Start("reduction");
        Int_t ntrk=0;
        // ====================================================
        // start looping on events cointained in the data file  
        // ====================================================         
        for (UInt_t re = evfrom; re < evfrom+nevents; re++){
            
            ev_count++;
            
            if ( p->DebugMode() && re%100 == 0 && re > 0 ) cout << ".";
            
            b_trk->GetEntry(re);
            b_header->GetEntry(re);
            pscu = header->GetPscuHeader();
                                
            if ( dbtime->DBabsTime(pscu->GetOrbitalTime()) > runtrailtime || dbtime->DBabsTime(pscu->GetOrbitalTime()) < runheadtime) {
                        
                if (p->VerboseMode()){ 
                    printf(" TrkCore - WARNING: event outside the run time window, skipping it\n");
                    cout << " OBT "<<pscu->GetOrbitalTime()<<" ABS-time "<<dbtime->DBabsTime(pscu->GetOrbitalTime())<<" RUN "<<runheadtime<<"-"<<runtrailtime<<endl;
                };              
            }else{
                                
                p->ProcessEvent(l0_event);
                                        
                // ----------------
                //  LEVEL1 output
                // ----------------
                if(p->get1){
                    if(p->ifroot1){ // root
                        l1_event->Clear();
                        l1_event->SetFromLevel1Struct(&level1event_);
                        //                              t_level1->Fill();
                    }else{       // hbook
                        throw -299;
                    };
                };
                // ----------------
                //  HOUGH output
                // ----------------
                if(p->geth){
                    if(p->ifrooth){ // root
//                      lh_event->Clear();
//                      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++;
                };
            };
        }; // end loop on events
        if(p->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; 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->get2){
        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
        
    };
    
    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);
}

1	mocchiut	1.1	/**
2			* \file TrkCore.cpp
3			* \author Elena Vannuccini
4			*/
5			// .........................................................
6			// ROOT header files
7			// .........................................................
8			#include <TTree.h>
9			#include <TClassEdit.h>
10			#include <TObject.h>
11			#include <TList.h>
12			#include <TSystem.h>
13			#include <TSystemDirectory.h>
14			#include <TString.h>
15			#include <TFile.h>
16			#include <TClass.h>
17			#include <TCanvas.h>
18			#include <TH1.h>
19			#include <TH1F.h>
20			#include <TH2D.h>
21			#include <TLatex.h>
22			#include <TPad.h>
23			#include <TSQLServer.h>
24			#include <TSQLRow.h>
25			#include <TSQLResult.h>
26	pam-fi	1.11	#include <TBenchmark.h>
27
28	mocchiut	1.1	// .........................................................
29			// other general header files
30			// .........................................................
31			#include <fstream>
32			#include <iostream>
33			// .........................................................
34			// files in the inc directory
35			// .........................................................
36			#include <RunInfo.h>
37			#include <GLTables.h>
38			#include <TrkVerl2.h>
39			#include <TrkProcess.h>
40	pam-fi	1.2	//#include <TrkStruct.h>
41	mocchiut	1.1	#include <TrkLevel2.h>
42	pam-fi	1.6	#include <TrkLevel1.h>
43	mocchiut	1.1	#include <TrkLevel0.h>
44			// .........................................................
45			// YODA header files
46			// .........................................................
47			#include <PamelaRun.h>
48	pam-fi	1.7	//#include <RegistryEvent.h>
49	mocchiut	1.1	#include <physics/tracker/TrackerEvent.h>
50			#include <CalibTrk1Event.h>
51			#include <CalibTrk2Event.h>
52	pam-fi	1.2
53	mocchiut	1.1	// .........................................................
54			// namespaces
55			// .........................................................
56			using namespace std;
57			using namespace pamela;
58			// ================================================================================
59			//
60			//
61			// ================================================================================
62			/**
63			* \brief Tracker data reduction routine.
64			*
65			* It performs data reduction LEVEL0->LEVEL1->LEVEL2, producing LEVEL1 and or LEVEL2 output, in ROOT or HBOOK format.
66			* Input parameters:
67			* @param run id of the run to be processed (if run=0 a whole file is reprocessed)
68			* @param dbc pointer to BD server
69			* @param file1 output LEVEL1 file name (null if no LEVEL1 output is required)
70			* @param file2 output LEVEL2 file name (null if no LEVEL2 output is required)
71			* @param standalone (flag to run the program in standalone mode, that is without reading RunInfo)
72			*/
73			//short int TrkCore(Int_t run, TSQLServer *dbc, TString file1, TString file2, Bool_t standalone)
74	pam-fi	1.7	int TrkCore(UInt_t run, TFile f2, TSQLServer dbc, int ncustom, char*vcustom[])
75	mocchiut	1.1	{
76			// ---------------------------
77			// Define some basic variables
78			// ---------------------------
79
80			TString filename = 0;
81			TString trkcalibfile = 0;
82			Int_t last_trk_calib_used = 0;
83	pam-fi	1.7	Long64_t nentries = 0LL;
84	pam-fi	1.2
85	mocchiut	1.1	// LEVEL0 input
86			TFile *f0 = 0;
87	pam-fi	1.2	TFile *f0_c = 0;
88	mocchiut	1.1	TTree *physicsTree = 0;
89			TBranch *b_trk = 0;
90			TBranch *b_header = 0;
91			EventHeader *header = 0;
92			PscuHeader *pscu = 0;
93			TrkLevel0 *l0_event = 0;
94			// RunInfo
95			ItoRunInfo *runinfo = 0;
96	pam-fi	1.7	TArrayI *runlist = 0;
97			UInt_t from_run;
98			UInt_t to_run;
99	mocchiut	1.1	Bool_t reprocessing = false;
100			//LEVEL2 output - ROOT
101			TTree *t_level2 = 0;
102			TTree *t_clone = 0;
103			TrkLevel2 *l2_event = new TrkLevel2();
104			TrkLevel2 *l2_clone = new TrkLevel2();
105	pam-fi	1.13	//LEVEL1 output - ROOT
106			TrkLevel1 *l1_event = new TrkLevel1();
107			// TrkHough *lh_event = new TrkHough();
108	mocchiut	1.1	// -----------------------
109			// -----------------------
110			// Handle input parameters
111			// (data)
112			//
113			// - open/create output files, determining the environment root/hbook from the estension
114			// - create the level1/level2 tree+branch/nt-uple
115			// -----------------------
116			// -----------------------
117			TrkProcess *p = new TrkProcess(run,f2);
118			p->HandleCustomPar(ncustom,vcustom);
119	pam-fi	1.2	if( p->DebugMode() )p->Dump();
120	mocchiut	1.1
121			// ===================================
122			// Open/Create level2 output file
123			// ===================================
124			if(p->get2){
125	pam-fi	1.13	//-------------------------------------------
126			// read RunInfo
127			//-------------------------------------------
128			if(!(p->standalone)){
129			// Open "RunInfo" tree and get list of run
130			runinfo = new ItoRunInfo(f2);
131			char *trkversion = TrkInfo(false);
132			int runinfo_error = runinfo->Update(run,"TRK",trkversion);
133			if( runinfo_error ) throw runinfo_error;
134			runlist = runinfo->GetRunList();
135			reprocessing = runinfo->IsReprocessing();//????
136			if(p->VerboseMode()){
137			cout << "#events "<< runinfo->GetFileEntries() << endl;// #events in the file
138			cout << "#runs "<< runinfo->GetRunEntries() << endl;// #runs in the file
139			cout << "1-st run "<< runlist->At(0) << endl;
140			cout << "last run "<< runlist->At(runinfo->GetRunEntries()-1) << endl;
141			cout << "1-st event "<< runinfo->GetFirstEntry() << endl;// first event of our run
142			cout << "last event+1 "<< runinfo->GetLastEntry() << endl;// first event AFTER the last event of our run
143			cout << "reprocessing "<< runinfo->IsReprocessing() << endl << endl;
144			};
145			};
146			//-------------------------------------------
147			//
148			//-------------------------------------------
149			// Take (if present) the old Tracker tree
150			t_clone = (TTree*)f2->Get("Tracker");
151			if( t_clone != NULL ) t_clone->SetBranchAddress("TrkLevel2",&l2_clone);
152			// Create NEW Tracker tree
153			t_level2 = new TTree("Tracker","PAMELA tracker level2 data ");
154			t_level2->Branch("TrkLevel2","TrkLevel2",&l2_event);
155			if(p->get1){
156			if(p->VerboseMode())cout << endl << "Requested LEVEL1 output" << endl;
157			t_level2->Branch("TrkLevel1","TrkLevel1",&l1_event);
158			t_level2->BranchRef();
159			}else if(p->geth){
160			if(p->VerboseMode())cout << endl << "Requested Hough-Transform output" << endl;
161			// t_level2->Branch("TrkHough","TrkHough",&lh_event);
162			// t_level2->BranchRef();
163	pam-fi	1.7	};
164	pam-fi	1.13	};
165	mocchiut	1.1
166
167			// -------------------------------------------
168			// define runs to be processed/reprocessed
169			// -------------------------------------------
170			if(run == 0){
171	pam-fi	1.13	// reprocessing ALL runs
172			if(p->standalone)throw -298; // reprocessing not implemented
173			from_run = runlist->At(0);
174			to_run = runlist->At(runinfo->GetRunEntries()-1);
175	mocchiut	1.1	}else{
176			// processing/reprocessing ONE single run
177	pam-fi	1.13	from_run = run;
178			to_run = run;
179	mocchiut	1.1	};
180
181			//
182			// init "expiration date" of calibrations and parameters
183			//
184	pam-fi	1.7	UInt_t pedsig_time =0;
185	pam-fi	1.13	UInt_t B_time =0;
186			UInt_t mip_time =0;
187			UInt_t charge_time =0;
188			UInt_t eta_time =0;
189			UInt_t align_time =0;
190			UInt_t mask_time =0;
191	mocchiut	1.1	//
192			// init event counter
193			//
194			Int_t ev_count =0;
195			//
196			// create query-results objects
197			//
198	pam-fi	1.7	GL_RUN q1 = GL_RUN();
199	mocchiut	1.1	GL_TRK_CALIB q2 = GL_TRK_CALIB();
200			GL_ROOT q3 = GL_ROOT();
201			GL_PARAM q4 = GL_PARAM();
202
203			// ------------------------------------------------------------
204			// if reprocessing one run, copy all the events BEFORE the run
205			// ------------------------------------------------------------
206			if( !(p->standalone) ){
207	pam-fi	1.13	for(UInt_t i=0; i<runinfo->GetFirstEntry(); i++){
208			t_clone->GetEntry(i);
209			l2_event = l2_clone;
210			t_level2->Fill();
211			l2_event->Clear();
212			// COPY COPY COPY
213			};
214	mocchiut	1.1	};
215			// ------------------------------------------------------------
216			// ------------------------------------------------------------
217			// START LOOP OVER RUNS TO PROCESSED/REPROCESSED
218			// ------------------------------------------------------------
219			// ------------------------------------------------------------
220	pam-fi	1.7	for(UInt_t idRun = from_run; idRun <= to_run; idRun++){
221	mocchiut	1.1
222	pam-fi	1.2	if(p->VerboseMode()) cout << endl<<" ========================= Run: "<< idRun << endl;
223	pam-fi	1.7	UInt_t runheadtime = 0;
224			UInt_t runtrailtime = 0;
225	mocchiut	1.1	UInt_t evfrom = 0;
226			UInt_t evto = 0;
227	pam-fi	1.7	UInt_t nevents = 0;
228			Int_t id_root_l0 =-1;
229	mocchiut	1.1	Int_t trk_calib_used = 0;
230	pam-fi	1.7
231	mocchiut	1.1	if(p->standalone){
232			// ==============================
233			// first query: retrieve run info
234			// ==============================
235	pam-fi	1.7	if (q1.Query_GL_RUN(idRun,dbc) )throw -50;
236			id_root_l0 = q1.ID_ROOT_L0;
237	mocchiut	1.1	runheadtime = q1.RUNHEADER_TIME;
238			runtrailtime = q1.RUNTRAILER_TIME;
239	pam-fi	1.7	evfrom = q1.EV_FROM;
240			evto = q1.EV_TO;
241	pam-fi	1.13	nevents = q1.NEVENTS;
242			trk_calib_used = q1.TRK_CALIB_USED;
243	mocchiut	1.1	}else{
244			// ==============================
245			// get run info from RunInfo tree
246			// ==============================
247	pam-fi	1.3	int runinfo_error = runinfo->GetRunInfo(idRun);
248	pam-fi	1.7	if( runinfo_error ) throw runinfo_error;
249			id_root_l0 = runinfo->ID_ROOT_L0;
250	mocchiut	1.1	runheadtime = runinfo->RUNHEADER_TIME;
251			runtrailtime = runinfo->RUNTRAILER_TIME;
252	pam-fi	1.7	evfrom = runinfo->EV_FROM;
253			evto = runinfo->EV_TO;
254	pam-fi	1.13	nevents = runinfo->NEVENTS;
255			trk_calib_used = runinfo->TRK_CALIB_USED;
256	mocchiut	1.1	};
257			//
258	pam-fi	1.2	if(p->VerboseMode()){
259	pam-fi	1.7	cout << "ROOT file ID "<< id_root_l0 << endl;
260	mocchiut	1.1	cout << "RunHeader time "<< runheadtime << endl;
261			cout << "RunTrailer time "<< runtrailtime << endl;
262	pam-fi	1.7	cout << " from event "<< evfrom << endl;
263			cout << " to event "<< evto << endl;
264			cout << " num. events "<< nevents << endl;
265			cout << "trk-calibration used "<< trk_calib_used << endl;
266	mocchiut	1.1	};
267			// ========================================================
268			// second query: search the LEVEL0 file that contains idRun
269			// ========================================================
270			TString lastfilename = filename;
271	pam-fi	1.7	if( q3.Query_GL_ROOT(id_root_l0,dbc) )throw -51;
272	mocchiut	1.1	filename = q3.PATH + q3.NAME;
273			// ========================================================
274			// Open the input LEVEL0 data file
275			// ========================================================
276			if(filename.CompareTo(lastfilename)){
277			if(!lastfilename.IsNull())f0->Close();
278			//if( debug ) cout << "Opening LEVEL0 file: "<< filename << endl;
279	pam-fi	1.2	if(p->VerboseMode()) cout << "Opening LEVEL0 file: "<< filename << endl;
280			FileStat_t t;
281			if( gSystem->GetPathInfo(filename.Data(),t) )throw -6;
282	mocchiut	1.1	f0 = new TFile(filename);
283			if ( !f0 ) throw -6;
284			physicsTree = (TTree*)f0->Get("Physics"); if(!physicsTree) throw -7;
285			b_header = physicsTree ->GetBranch("Header"); if(!b_header ) throw -8;
286	pam-fi	1.7	// b_registry = physicsTree ->GetBranch("Registry"); if(!b_registry ) throw -9;
287	mocchiut	1.1	b_trk = physicsTree ->GetBranch("Tracker"); if(!b_trk ) throw -203;
288			physicsTree->SetBranchAddress("Tracker" ,&l0_event);
289	pam-fi	1.7	// physicsTree->SetBranchAddress("Registry",&reg);
290	mocchiut	1.1	physicsTree->SetBranchAddress("Header",&header);
291
292	pam-fi	1.7	nentries = physicsTree->GetEntries();
293			if ( nentries < 1 ) throw -11;
294			if ( nentries < (evto+1)) throw -12;
295
296	mocchiut	1.1	};
297
298	pam-fi	1.7	GL_TIMESYNC *dbtime = new GL_TIMESYNC(id_root_l0,"ID",dbc);
299
300	mocchiut	1.1	// =============================================================
301	pam-fi	1.2	// retrieve calibration file needed to reduce data
302	mocchiut	1.1	// =============================================================
303	pam-fi	1.2	// if run OBT is > last calibration "expiration date"
304			// - search for new calibration packet
305			// - load calibration parameters (full + truncated)
306	pam-fi	1.5	if(p->VerboseMode() )cout << "Full pedestals for cluster finding:";
307	pam-fi	1.2	if(runheadtime > pedsig_time){
308			if( q2.Query_GL_TRK_CALIB(runheadtime,dbc) )throw -53;
309			pedsig_time = q2.TO_TIME;
310	pam-fi	1.7	if(q2.EV_ROOT_CALIBTRK1 != q2.EV_ROOT_CALIBTRK2)
311			printf("WARNING!! ---> EV_ROOT_CALIBTRK1=%d it's different from EV_ROOT_CALIBTRK2=%d \n\n",q2.EV_ROOT_CALIBTRK1,q2.EV_ROOT_CALIBTRK2);
312			if( q3.Query_GL_ROOT(q2.ID_ROOT_L0,dbc) )throw -51;
313	pam-fi	1.5	trkcalibfile = q3.PATH + q3.NAME;
314			if(p->VerboseMode()) cout << " >> Loading from LEVEL0 file: "<< trkcalibfile << endl;
315	pam-fi	1.2	if(trkcalibfile.CompareTo(filename)){
316	pam-fi	1.5	FileStat_t t;
317			if( gSystem->GetPathInfo(trkcalibfile.Data(),t) )throw -6;
318			f0_c=new TFile(trkcalibfile);
319			if ( !f0_c ) throw -6;
320	pam-fi	1.2	}else{
321	pam-fi	1.5	f0_c = f0;
322	pam-fi	1.2	};
323			if(p->VerboseMode()) {
324	pam-fi	1.7	cout << " calibration entries "<< q2.EV_ROOT_CALIBTRK1 << " " << q2.EV_ROOT_CALIBTRK2;
325			cout << " (from time "<< q2.FROM_TIME <<" to time "<< q2.TO_TIME <<")"<<endl;
326	pam-fi	1.2	};
327	pam-fi	1.7	pedsigbad_.FillACalibFrom(f0_c,q2.EV_ROOT_CALIBTRK1,q2.EV_ROOT_CALIBTRK2);
328	pam-fi	1.5	}else{
329			if(p->VerboseMode() )cout<< " already loaded. "<< endl;
330	pam-fi	1.2	};
331			if( p->DebugMode() ) for(int i=0; i<12; i++) cout << " DSP "<< i << " "<< pedsigbad_.pedestal[64][12][i] << endl;
332			// =============================================================
333			// retrieve calibration file needed to uncompress data
334			// =============================================================
335			// if the run was compressed using default calib
336			// load truncated pedestals from default
337			// otherwise reload them from on-line calibration
338			if(p->VerboseMode() )cout << "Truncated pedestals for uncompression:";
339	mocchiut	1.1	if(trk_calib_used==104 && last_trk_calib_used!=104){
340
341	pam-fi	1.5	if(p->VerboseMode() )cout << " >> Loading default calibration " << endl;
342	pam-fi	1.7	if (q4.Query_GL_PARAM(runheadtime,7,dbc) )throw -52;
343	pam-fi	1.2	pedsigbad_.FillTCalibFrom(q4.PATH+q4.NAME);
344
345			}else if( (trk_calib_used==1 \| trk_calib_used==2) && last_trk_calib_used==104 ){
346
347			if(p->VerboseMode()) cout << ">> Re-loading on-line calibration " << endl;
348	pam-fi	1.7	pedsigbad_.FillTCalibFrom(f0_c,q2.EV_ROOT_CALIBTRK1,q2.EV_ROOT_CALIBTRK2);
349	pam-fi	1.2	}else{
350			if(p->VerboseMode()) cout << ">> Using on-line calibration " << endl;
351			};
352			if( p->DebugMode() ) for(int i=0; i<12; i++) cout << " DSP "<< i << " "<< pedsigbad_.pedestal_t[64][12][i] << endl;
353			last_trk_calib_used = trk_calib_used;
354
355	mocchiut	1.1	// =============================================================
356			// retrieve information about parameters to process LEVEL2
357			// =============================================================
358			//
359			// magnetic field
360			//
361	pam-fi	1.2	if(p->VerboseMode()) cout << "Loading parameter files : "<< endl;
362
363	mocchiut	1.1	if(runheadtime > B_time){
364	pam-fi	1.7	if( q4.Query_GL_PARAM(runheadtime,1,dbc) )throw -52;
365	mocchiut	1.1	B_time = q4.TO_TIME;
366			l2_event->LoadField(q4.PATH+q4.NAME);
367	pam-fi	1.2	// l2_event->LoadField(q4.PATH+q4.NAME);
368	mocchiut	1.1	if(path_.error) throw -215;
369			};
370			//
371			// mip conversion parameters
372			//
373			if(runheadtime > mip_time){
374	pam-fi	1.7	if( q4.Query_GL_PARAM(runheadtime,2,dbc) )throw -52;
375	mocchiut	1.1	mip_time = q4.TO_TIME;
376			path_.FillWith(q4.PATH+q4.NAME);
377			readmipparam_();
378			if(path_.error) throw -212;
379			};
380			//
381			// charge correlation parameters
382			//
383			if(runheadtime > charge_time){
384	pam-fi	1.7	if( q4.Query_GL_PARAM(runheadtime,3,dbc) )throw -52;
385	mocchiut	1.1	charge_time = q4.TO_TIME;
386			path_.FillWith(q4.PATH+q4.NAME);
387			readchargeparam_();
388			if(path_.error) throw -213;
389			};
390			//
391			// eta p.f.a. parameters
392			//
393			if(runheadtime > eta_time){
394	pam-fi	1.7	if( q4.Query_GL_PARAM(runheadtime,4,dbc) )throw -52;
395	mocchiut	1.1	eta_time = q4.TO_TIME;
396			path_.FillWith(q4.PATH+q4.NAME);
397			readetaparam_();
398			if(path_.error) throw -214;
399			};
400			//
401			// alignment parameters
402			//
403			if(runheadtime > align_time){
404	pam-fi	1.7	if( q4.Query_GL_PARAM(runheadtime,5,dbc) )throw -52;
405	mocchiut	1.1	align_time = q4.TO_TIME;
406			path_.FillWith(q4.PATH+q4.NAME);
407			readalignparam_();
408			if(path_.error) throw -211;
409			};
410			//
411			// viking mask
412			//
413			if(runheadtime > mask_time){
414	pam-fi	1.7	if( q4.Query_GL_PARAM(runheadtime,6,dbc) )throw -52;
415	mocchiut	1.1	mask_time = q4.TO_TIME;
416			path_.FillWith(q4.PATH+q4.NAME);
417			readvkmask_();
418			if(path_.error) throw -210;
419			};
420
421	pam-fi	1.11	TBenchmark *reduction = new TBenchmark();
422			if(p->VerboseMode())reduction->Start("reduction");
423			Int_t ntrk=0;
424	mocchiut	1.1	// ====================================================
425			// start looping on events cointained in the data file
426			// ====================================================
427	pam-fi	1.7	for (UInt_t re = evfrom; re < evfrom+nevents; re++){
428	mocchiut	1.1
429			ev_count++;
430
431	pam-fi	1.13	if ( p->DebugMode() && re%100 == 0 && re > 0 ) cout << ".";
432	mocchiut	1.1
433	pam-fi	1.13	b_trk->GetEntry(re);
434			b_header->GetEntry(re);
435			pscu = header->GetPscuHeader();
436	pam-fi	1.7
437	pam-fi	1.13	if ( dbtime->DBabsTime(pscu->GetOrbitalTime()) > runtrailtime \|\| dbtime->DBabsTime(pscu->GetOrbitalTime()) < runheadtime) {
438	pam-fi	1.7
439	pam-fi	1.14	if (p->VerboseMode()){
440			printf(" TrkCore - WARNING: event outside the run time window, skipping it\n");
441			cout << " OBT "<<pscu->GetOrbitalTime()<<" ABS-time "<<dbtime->DBabsTime(pscu->GetOrbitalTime())<<" RUN "<<runheadtime<<"-"<<runtrailtime<<endl;
442			};
443	pam-fi	1.13	}else{
444	pam-fi	1.7
445	pam-fi	1.13	p->ProcessEvent(l0_event);
446	pam-fi	1.7
447	pam-fi	1.13	// ----------------
448			// LEVEL1 output
449			// ----------------
450			if(p->get1){
451			if(p->ifroot1){ // root
452			l1_event->Clear();
453			l1_event->SetFromLevel1Struct(&level1event_);
454			// t_level1->Fill();
455			}else{ // hbook
456			throw -299;
457			};
458			};
459			// ----------------
460			// HOUGH output
461			// ----------------
462			if(p->geth){
463			if(p->ifrooth){ // root
464			// lh_event->Clear();
465			// lh_event->SetFromHoughStruct(&houghevent_);
466			}else{ // hbook
467			throw -299;
468			};
469			};
470			// ----------------
471			// LEVEL2 output
472			// ----------------
473			if(p->get2){
474			l2_event->Clear();
475			if(p->get1) l2_event->SetFromLevel2Struct(&level2event_,l1_event);//set references to level1
476			else l2_event->SetFromLevel2Struct(&level2event_);
477	pam-fi	1.10	// l2_event->Dump();
478	pam-fi	1.13	t_level2->Fill();
479			if( l2_event->ntrk()>0 )ntrk++;
480	mocchiut	1.1	};
481	pam-fi	1.13	};
482	mocchiut	1.1	}; // end loop on events
483	pam-fi	1.11	if(p->VerboseMode()){
484	pam-fi	1.13	cout << " >>> processed "<< ev_count <<" events"<< endl;
485			if(ev_count)cout << ntrk << " events with at least one track ("<<(Int_t)(100*ntrk/ev_count)<<"%)\n";
486			reduction->Show("reduction");
487	pam-fi	1.11	}
488			delete reduction;
489	pam-fi	1.7
490			delete dbtime;
491
492	pam-fi	1.13	}; // end loop on runs
493
494
495	mocchiut	1.1	// ------------------------------------------------------------
496			// if reprocessing one run, copy all the events AFTER the run
497			// ------------------------------------------------------------
498			if( !(p->standalone) ){
499			for(UInt_t i=runinfo->GetLastEntry()+1; i<runinfo->GetFileEntries(); i++){
500			t_clone->GetEntry(i);
501			l2_event = l2_clone;
502			t_level2->Fill();
503			l2_event->Clear();
504			// COPY COPY COPY
505			};
506			};
507			// ---------------
508			// close the files
509			// ---------------
510			if(p->get2){
511	pam-fi	1.13	if( t_clone )t_clone->Delete("all");//delete old tree from file
512			if( !(p->standalone) )runinfo->Close();
513			f2->Write("Tracker");
514			if( t_level2 )t_level2->Delete(); //delete new tree from memory
515
516			};
517	pam-fi	1.6
518	pam-fi	1.13	if( f0->IsOpen() ) f0->Close();
519	pam-fi	1.2	if( f0_c->IsOpen() )f0_c->Close();
520	pam-fi	1.13
521			// lh_event->Delete();
522			l1_event->Delete();
523			l2_event->Delete();
524			l2_clone->Delete();
525	pam-fi	1.2
526			return(p->ostatus);
527	mocchiut	1.1	}
528