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
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	}