doc/examples-1/GP.cpp


#if !defined(__CINT__) || defined(__MAKECINT__)


// ROOT includes //
#include <TString.h>
#include <TH1F.h>
#include <TH2F.h>

// standard libraries includes //
#include <stdlib.h>
#include <iostream>
//#include <fstream>
#include <iomanip>
using namespace std;

// PAMELA includes //
#include <PamLevel2.h>

#endif

//======================
// global variables
//======================
Int_t tot=0;
Int_t sel[]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};


//======================
// HISTOGRAMS declaration
//======================
TH1F* hp0gen;
TH1F* hp0trk;
TH1F* hpmtrk;


//===============================================================================
//
//
//
//
//===============================================================================
bool Select( PamLevel2* event ){
//    return false;
    TrkParams::SetVerboseMode();

    //-------------------------------------------
    // check if the event contains selection info
    //-------------------------------------------
    if( event->GetTrkLevel2()==0 ){
        cout << "Missing TrkLevel2 object "<<endl;
        return false;
    }
    if( event->GetCaloLevel2()==0 ){
        cout << "Missing CaloLevel2 object "<<endl;
        return false;
    }
    if( event->GetRunInfo()==0 ){
        cout << "Missing RunInfo object "<<endl;
        return false;
    }
    if( event->GetGPamela()==0 ){
        cout << "Missing GPamela object "<<endl;
        return false;
    }

    tot++;


    // -------------------------------------------------
    // no hits on CAT and CARD
    // -------------------------------------------------
    if( event->GetGPamela()->Nthcat > 0 ||  event->GetGPamela()->Nthcard > 0) return false;
    sel[0]++;

    // --------------------------------------------------
    // select simulated events where the primary particle
    // traverse S2 and S3
    // --------------------------------------------------
    bool  S2OK = false; // the primary particle hits S3
    bool  S3OK = false; // the primary particle hits S3

    for(int ih=0; ih<event->GetGPamela()->Nthtof; ih++){ //loop over tof hits

        // S2
        if( 
            (Int_t) event->GetGPamela()->Ipartof[ih] == event->GetGPamela()->Ipa &&                                //primary type
            ( (Int_t) event->GetGPamela()->Ipltof[ih] == 4 || (Int_t) event->GetGPamela()->Ipltof[ih] == 3) &&     //S21 or S22
            fabs( (event->GetGPamela()->P0 - event->GetGPamela()->P0tof[ih]) / event->GetGPamela()->P0 ) < 0.05 && //energy consistent with the primary (1%)
            true){
            S2OK = true;
        }
        // S3
        if( 
            (Int_t) event->GetGPamela()->Ipartof[ih] == event->GetGPamela()->Ipa &&                                //primary type
            ( (Int_t) event->GetGPamela()->Ipltof[ih] == 5 || (Int_t) event->GetGPamela()->Ipltof[ih] == 6) &&     //S31 or S32
            fabs( (event->GetGPamela()->P0 - event->GetGPamela()->P0tof[ih]) / event->GetGPamela()->P0 ) < 0.05 && //energy consistent with the primary (1%)
            true){
            S3OK = true;
        }
    }
    if(!S2OK || !S3OK)return false;
    sel[1]++;
    

    return true;
    

}
//===============================================================
// Create histograms
//
//
//
//
//
//===============================================================
void CreateHistos( PamLevel2* event , TFile* outf ){
    cout << "Creating histos..."<<endl;
    gROOT->cd();//create histos in memory
    TString hid;
    TString title;

    //===========================================
    // rigidity binning
    Int_t   nbins = 10;
    Float_t rmin  = 1.;
    Float_t rmax  = 100.;
    Float_t dr = (log10(rmax)-log10(rmin))/nbins;
    TArrayD* arr = new TArrayD(nbins+1);
    arr->AddAt((Double_t)rmin,0);
    for(Int_t i=1; i<=nbins; i++){
        Double_t previous = arr->At(i-1);
        Double_t value=0;
        value = (Double_t) pow(10,log10((Float_t)previous)+dr);
        arr->AddAt(value,i);
    };
    Double_t *rbinpt = arr->GetArray();
    //===========================================

    //--------------------------------
    // histos
    //--------------------------------
    hp0gen = new TH1F("hp0gen","Generated momentum (all events)",nbins,rbinpt);
    hp0trk = new TH1F("hp0trk","Generated momentum ",nbins,rbinpt);
    hpmtrk = new TH1F("hpmtrk","Measured  momentum ",nbins,rbinpt);


    //--------------------------------
    // initialization
    //--------------------------------
    cout << "Opening DB connection..."<<endl;
    event->SetDBConnection();

}
//===============================================================
// Retrieve histograms
//
//
//
//
//
//===============================================================
void RetrieveHistos(){

    gROOT->cd();
    TString hid;
    hp0gen=dynamic_cast<TH1F*>(gDirectory->FindObject("hp0gen"));
    hp0trk=dynamic_cast<TH1F*>(gDirectory->FindObject("hp0trk"));
    hpmtrk=dynamic_cast<TH1F*>(gDirectory->FindObject("hpmtrk"));

}
//===============================================================
// Get histograms (... non ho capito la differenza con retrieve)
//
//
//
//
//
//===============================================================
void GetHistos(){

    TString hid;
}
//===============================================================
// Fill histograms (called event-by-event)
//
//
//
//
//
//
//=============================================================== 
void FillHistos( PamLevel2* event ){
//  ----------------------------------
//  retrieve histos
//  ----------------------------------
    RetrieveHistos();


    cout << "================================================="<<endl;
    cout << "part        "<<event->GetGPamela()->Ipa <<endl;
    cout << "P0 (GeV/c)  "<<event->GetGPamela()->P0 <<endl;
    cout << "X0 (cm)     "<<event->GetGPamela()->X0 << endl;
    cout << "Y0 (cm)     "<<event->GetGPamela()->Y0    << endl; 
    cout << "Z0 (cm)     "<<event->GetGPamela()->Z0    << endl; 
    cout << "Theta (rad) "<<event->GetGPamela()->Theta  << endl;
    cout << "Phi (rad)   "<<event->GetGPamela()->Phi    << endl;
    cout <<"Spectrometer hits: "<< endl;
    cout << "hit --";
    cout << "part  ";
    cout << "plane ";
    cout << "X           ";
    cout << "Y           ";
    cout << "Z           ";
    cout << "En.release  ";
    cout << "P           ";
    cout << endl;
    for(int ih=0; ih<event->GetGPamela()->Nthspe; ih++){
        cout << setw(4) << ih << "--";
        cout << setw(6) <<   (Int_t) event->GetGPamela()->Iparspe[ih];//part
        cout << setw(6) <<   (Int_t) event->GetGPamela()->Itrpb[ih];//0-5 o 1-6 piano
        cout << setw(12)<< event->GetGPamela()->Xavspe[ih];
        cout << setw(12)<< event->GetGPamela()->Yavspe[ih];
        cout << setw(12)<< event->GetGPamela()->Zavspe[ih];
        cout << setw(12)<< event->GetGPamela()->Erelspe[ih];//en.release
        cout << setw(12)<< event->GetGPamela()->P0spe[ih];//momentum
        cout << endl;
    }
    

    hp0gen->Fill(event->GetGPamela()->P0);
    
    // one fitted track
    if( event->GetTrkLevel2()->GetNTracks()>0 ){
        sel[2]++;
        hp0trk->Fill(event->GetGPamela()->P0);

//  ---------------------------------
//  Retrieve the track
//  ---------------------------------
        event->SetSortingMethod("+CAL");        // use only calorimeter for image rejection
        PamTrack  *track = event->GetTrack(0);  // get the whole track 
        TrkTrack  *trk   = track->GetTrkTrack();// get tracker track variables
        
        hpmtrk->Fill(trk->GetRigidity());
        
    }
        
}
//===============================================================
// Save histograms
//
//
//
//
//
//===============================================================
void SaveHistos( TFile * outf ){

TString hid; 

gROOT->cd();
//  ----------------------------------
//  retrieve histos
//  ----------------------------------
    RetrieveHistos();

//  ----------------------------------
//  save on file
//  ----------------------------------
    outf->cd();

    hp0gen->Write();
    hp0trk->Write();
    hpmtrk->Write();

    cout << "####------------------------------------------------------------------" <<endl;
    cout << "#### total :"<<tot<<endl;
    if(!tot)return;
    cout << "#### no GP hits on CAT/CARD   :"<<sel[0]<< "  (" << ((float)sel[0]/(float)tot) <<")"<<endl;
    cout << "#### GP hits on S2/S3         :"<<sel[1]<< "  (" << ((float)sel[1]/(float)tot) <<")"<<endl;
    cout << "#### 1 fitted track           :"<<sel[2]<< "  (" << ((float)sel[2]/(float)tot) <<")"<<endl;
    cout << "####------------------------------------------------------------------" <<endl;

}
1	pam-fi	1.1
2			#if !defined(__CINT__) \|\| defined(__MAKECINT__)
3
4
5			// ROOT includes //
6			#include <TString.h>
7			#include <TH1F.h>
8			#include <TH2F.h>
9
10			// standard libraries includes //
11			#include <stdlib.h>
12			#include <iostream>
13			//#include <fstream>
14			#include <iomanip>
15			using namespace std;
16
17			// PAMELA includes //
18			#include <PamLevel2.h>
19
20			#endif
21
22			//======================
23			// global variables
24			//======================
25			Int_t tot=0;
26			Int_t sel[]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
27
28
29			//======================
30			// HISTOGRAMS declaration
31			//======================
32			TH1F* hp0gen;
33			TH1F* hp0trk;
34			TH1F* hpmtrk;
35
36
37
38			//===============================================================================
39			//
40			//
41			//
42			//
43			//===============================================================================
44			bool Select( PamLevel2* event ){
45			// return false;
46			TrkParams::SetVerboseMode();
47
48			//-------------------------------------------
49			// check if the event contains selection info
50			//-------------------------------------------
51			if( event->GetTrkLevel2()==0 ){
52			cout << "Missing TrkLevel2 object "<<endl;
53			return false;
54			}
55			if( event->GetCaloLevel2()==0 ){
56			cout << "Missing CaloLevel2 object "<<endl;
57			return false;
58			}
59			if( event->GetRunInfo()==0 ){
60			cout << "Missing RunInfo object "<<endl;
61			return false;
62			}
63			if( event->GetGPamela()==0 ){
64			cout << "Missing GPamela object "<<endl;
65			return false;
66			}
67
68			tot++;
69
70
71			// -------------------------------------------------
72			// no hits on CAT and CARD
73			// -------------------------------------------------
74			if( event->GetGPamela()->Nthcat > 0 \|\| event->GetGPamela()->Nthcard > 0) return false;
75			sel[0]++;
76
77			// --------------------------------------------------
78			// select simulated events where the primary particle
79			// traverse S2 and S3
80			// --------------------------------------------------
81			bool S2OK = false; // the primary particle hits S3
82			bool S3OK = false; // the primary particle hits S3
83
84			for(int ih=0; ih<event->GetGPamela()->Nthtof; ih++){ //loop over tof hits
85
86			// S2
87			if(
88			(Int_t) event->GetGPamela()->Ipartof[ih] == event->GetGPamela()->Ipa && //primary type
89			( (Int_t) event->GetGPamela()->Ipltof[ih] == 4 \|\| (Int_t) event->GetGPamela()->Ipltof[ih] == 3) && //S21 or S22
90			fabs( (event->GetGPamela()->P0 - event->GetGPamela()->P0tof[ih]) / event->GetGPamela()->P0 ) < 0.05 && //energy consistent with the primary (1%)
91			true){
92			S2OK = true;
93			}
94			// S3
95			if(
96			(Int_t) event->GetGPamela()->Ipartof[ih] == event->GetGPamela()->Ipa && //primary type
97			( (Int_t) event->GetGPamela()->Ipltof[ih] == 5 \|\| (Int_t) event->GetGPamela()->Ipltof[ih] == 6) && //S31 or S32
98			fabs( (event->GetGPamela()->P0 - event->GetGPamela()->P0tof[ih]) / event->GetGPamela()->P0 ) < 0.05 && //energy consistent with the primary (1%)
99			true){
100			S3OK = true;
101			}
102			}
103			if(!S2OK \|\| !S3OK)return false;
104			sel[1]++;
105
106
107			return true;
108
109
110			}
111			//===============================================================
112			// Create histograms
113			//
114			//
115			//
116			//
117			//
118			//===============================================================
119			void CreateHistos( PamLevel2* event , TFile* outf ){
120			cout << "Creating histos..."<<endl;
121			gROOT->cd();//create histos in memory
122			TString hid;
123			TString title;
124
125			//===========================================
126			// rigidity binning
127			Int_t nbins = 10;
128			Float_t rmin = 1.;
129			Float_t rmax = 100.;
130			Float_t dr = (log10(rmax)-log10(rmin))/nbins;
131			TArrayD* arr = new TArrayD(nbins+1);
132			arr->AddAt((Double_t)rmin,0);
133			for(Int_t i=1; i<=nbins; i++){
134			Double_t previous = arr->At(i-1);
135			Double_t value=0;
136			value = (Double_t) pow(10,log10((Float_t)previous)+dr);
137			arr->AddAt(value,i);
138			};
139			Double_t *rbinpt = arr->GetArray();
140			//===========================================
141
142			//--------------------------------
143			// histos
144			//--------------------------------
145			hp0gen = new TH1F("hp0gen","Generated momentum (all events)",nbins,rbinpt);
146			hp0trk = new TH1F("hp0trk","Generated momentum ",nbins,rbinpt);
147			hpmtrk = new TH1F("hpmtrk","Measured momentum ",nbins,rbinpt);
148
149
150
151			//--------------------------------
152			// initialization
153			//--------------------------------
154			cout << "Opening DB connection..."<<endl;
155			event->SetDBConnection();
156
157			}
158			//===============================================================
159			// Retrieve histograms
160			//
161			//
162			//
163			//
164			//
165			//===============================================================
166			void RetrieveHistos(){
167
168			gROOT->cd();
169			TString hid;
170			hp0gen=dynamic_cast<TH1F*>(gDirectory->FindObject("hp0gen"));
171			hp0trk=dynamic_cast<TH1F*>(gDirectory->FindObject("hp0trk"));
172			hpmtrk=dynamic_cast<TH1F*>(gDirectory->FindObject("hpmtrk"));
173
174			}
175			//===============================================================
176			// Get histograms (... non ho capito la differenza con retrieve)
177			//
178			//
179			//
180			//
181			//
182			//===============================================================
183			void GetHistos(){
184
185			TString hid;
186			}
187			//===============================================================
188			// Fill histograms (called event-by-event)
189			//
190			//
191			//
192			//
193			//
194			//
195			//===============================================================
196			void FillHistos( PamLevel2* event ){
197			// ----------------------------------
198			// retrieve histos
199			// ----------------------------------
200			RetrieveHistos();
201
202
203
204
205			cout << "================================================="<<endl;
206			cout << "part "<<event->GetGPamela()->Ipa <<endl;
207			cout << "P0 (GeV/c) "<<event->GetGPamela()->P0 <<endl;
208			cout << "X0 (cm) "<<event->GetGPamela()->X0 << endl;
209			cout << "Y0 (cm) "<<event->GetGPamela()->Y0 << endl;
210			cout << "Z0 (cm) "<<event->GetGPamela()->Z0 << endl;
211			cout << "Theta (rad) "<<event->GetGPamela()->Theta << endl;
212			cout << "Phi (rad) "<<event->GetGPamela()->Phi << endl;
213			cout <<"Spectrometer hits: "<< endl;
214			cout << "hit --";
215			cout << "part ";
216			cout << "plane ";
217			cout << "X ";
218			cout << "Y ";
219			cout << "Z ";
220			cout << "En.release ";
221			cout << "P ";
222			cout << endl;
223			for(int ih=0; ih<event->GetGPamela()->Nthspe; ih++){
224			cout << setw(4) << ih << "--";
225			cout << setw(6) << (Int_t) event->GetGPamela()->Iparspe[ih];//part
226			cout << setw(6) << (Int_t) event->GetGPamela()->Itrpb[ih];//0-5 o 1-6 piano
227			cout << setw(12)<< event->GetGPamela()->Xavspe[ih];
228			cout << setw(12)<< event->GetGPamela()->Yavspe[ih];
229			cout << setw(12)<< event->GetGPamela()->Zavspe[ih];
230			cout << setw(12)<< event->GetGPamela()->Erelspe[ih];//en.release
231			cout << setw(12)<< event->GetGPamela()->P0spe[ih];//momentum
232			cout << endl;
233			}
234
235
236
237
238
239			hp0gen->Fill(event->GetGPamela()->P0);
240
241			// one fitted track
242			if( event->GetTrkLevel2()->GetNTracks()>0 ){
243			sel[2]++;
244			hp0trk->Fill(event->GetGPamela()->P0);
245
246			// ---------------------------------
247			// Retrieve the track
248			// ---------------------------------
249			event->SetSortingMethod("+CAL"); // use only calorimeter for image rejection
250			PamTrack *track = event->GetTrack(0); // get the whole track
251			TrkTrack *trk = track->GetTrkTrack();// get tracker track variables
252
253			hpmtrk->Fill(trk->GetRigidity());
254
255			}
256
257			}
258			//===============================================================
259			// Save histograms
260			//
261			//
262			//
263			//
264			//
265			//===============================================================
266			void SaveHistos( TFile * outf ){
267
268			TString hid;
269
270			gROOT->cd();
271			// ----------------------------------
272			// retrieve histos
273			// ----------------------------------
274			RetrieveHistos();
275
276			// ----------------------------------
277			// save on file
278			// ----------------------------------
279			outf->cd();
280
281			hp0gen->Write();
282			hp0trk->Write();
283			hpmtrk->Write();
284
285			cout << "####------------------------------------------------------------------" <<endl;
286			cout << "#### total :"<<tot<<endl;
287			if(!tot)return;
288			cout << "#### no GP hits on CAT/CARD :"<<sel[0]<< " (" << ((float)sel[0]/(float)tot) <<")"<<endl;
289			cout << "#### GP hits on S2/S3 :"<<sel[1]<< " (" << ((float)sel[1]/(float)tot) <<")"<<endl;
290			cout << "#### 1 fitted track :"<<sel[2]<< " (" << ((float)sel[2]/(float)tot) <<")"<<endl;
291			cout << "####------------------------------------------------------------------" <<endl;
292
293			}