/[PAMELA software]/PamelaLevel2/doc/examples/Histo-set-0.cpp
ViewVC logotype

Annotation of /PamelaLevel2/doc/examples/Histo-set-0.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.2 - (hide annotations) (download)
Mon Dec 11 18:29:01 2006 UTC (17 years, 11 months ago) by pam-fi
Branch: MAIN
CVS Tags: HEAD
Changes since 1.1: +0 -0 lines
FILE REMOVED
*** empty log message ***

1 pam-fi 1.1 ////////////////////////////////////////////////////////////////
2     //
3     // Histogram set for tracker-efficiency study.
4     //
5     // It contains the functions:
6     //
7     // - CreateHistos()
8     // - FillHistos(PamLevel2* event)
9     // - Save Histos()
10     //
11     // Usage (example):
12     // [] .L My-histo-set-macro.C+
13     // [] .L My-selection.C+
14     // [] .L Loop.C+
15     // [] Loop("data-dir/","file-list.txt",100000,"+ALL fillTree fillHistos","output-file.root")
16     //
17     ////////////////////////////////////////////////////////////////
18    
19     #if !defined(__CINT__) || defined(__MAKECINT__)
20    
21     #include <TString.h>
22     #include <TH1F.h>
23     #include <TH2F.h>
24     #include <TMath.h>
25     #include <TLine.h>
26     #include <TPolyMarker.h>
27     #include <TSelector.h>
28     #include <TFile.h>
29    
30     #include <stdlib.h>
31     #include <iostream>
32     using namespace std;
33    
34     #include <PamLevel2.h>
35     #include <CaloAxis.h>
36    
37     #endif
38     //======================
39     // HITOGRAMS declaration
40     //======================
41     TH1F *htofresx;
42     TH1F *htofresy;
43     TH1F *hcalresx;
44     TH1F *hcalresy;
45     TH1F *htrkresx[10];
46     TH1F *htrkresy[10];
47     TH1F *hxcalotrk[6];
48     TH1F *hycalotrk[6];
49     TH1F *hresangx;
50     TH1F *hresangy;
51     //===============================================================
52     // Create histograms
53     //===============================================================
54     void CreateHistos( ){
55    
56     gROOT->cd();//create histos in memory
57    
58     gDirectory->Delete("htofresx");
59     gDirectory->Delete("htofresy");
60     TH1F *htofresx = new TH1F("htofresx","Spatial residuals (S12+S32) - S21 (x)",200,-10.,10.);
61     TH1F *htofresy = new TH1F("htofresy","Spatial residuals (S11+S31) - S22 (y)",200,-10.,10.);
62    
63    
64     gDirectory->Delete("hcalresx");
65     gDirectory->Delete("hcalresy");
66     TH1F *hcalresx = new TH1F("hcalresx","Spatial residuals on CAL - S32 (x)",200,-10.,10.);
67     TH1F *hcalresy = new TH1F("hcalresy","Spatial residuals on CAL - S31 (y)",200,-10.,10.);
68    
69     TString hid;
70     for(Int_t i=0; i<3; i++){
71     hid.Form("htrkresx%i",i);
72     htrkresx[i] = new TH1F(hid.Data(),"Spatial residuals on TRK - TOF (x)",200,-10.,10.);
73     hid.Form("htrkresy%i",i);
74     htrkresy[i] = new TH1F(hid.Data(),"Spatial residuals on TRK - TOF (y)",200,-10.,10.);
75     }
76     for(Int_t i=3; i<9; i++){
77     hid.Form("htrkresx%i",i);
78     htrkresx[i] = new TH1F(hid.Data(),"Spatial residuals on TRK (x)",5000,-.2,.2);
79     hid.Form("htrkresy%i",i);
80     htrkresy[i] = new TH1F(hid.Data(),"Spatial residuals on TRK (y)",5000,-.2,.2);
81     }
82     for(Int_t i=9; i<10; i++){
83     hid.Form("htrkresx%i",i);
84     htrkresx[i] = new TH1F(hid.Data(),"Spatial residuals on TRK - CAL (x)",800,-10.,10.);
85     hid.Form("htrkresy%i",i);
86     htrkresy[i] = new TH1F(hid.Data(),"Spatial residuals on TRK - CAL (y)",800,-10.,10.);
87     }
88    
89    
90     for(Int_t ip=0; ip<6; ip++){
91     hid.Form("hxcalotrk%i",ip);
92     hxcalotrk[ip] = new TH1F(hid.Data(),"Calo axis intersaction on Trk plane (x)",500,-25.,25.);
93     hid.Form("hycalotrk%i",ip);
94     hycalotrk[ip] = new TH1F(hid.Data(),"Calo axis intersaction on Trk plane (y)",500,-25.,25.);
95     }
96    
97     hresangx = new TH1F("hresangx","Angular residual on the top calo plane (x)",400,-40.,40.);
98     hresangy = new TH1F("hresangy","Angular residual on the top calo plane (y)",400,-40.,40.);
99    
100     }
101    
102     //===============================================================
103     // Fill histograms (called event-by-event)
104     //===============================================================
105    
106     void FillHistos( PamLevel2* event ){
107     // ----------------------------------
108     // retrieve histos
109     // ----------------------------------
110     htofresx = dynamic_cast<TH1F*>(gDirectory->FindObject("htofresx"));
111     htofresy = dynamic_cast<TH1F*>(gDirectory->FindObject("htofresy"));
112     // cout << hex << htofresx << dec << endl;
113     hcalresx = dynamic_cast<TH1F*>(gDirectory->FindObject("hcalresx"));
114     hcalresy = dynamic_cast<TH1F*>(gDirectory->FindObject("hcalresy"));
115     TString hid;
116     for(Int_t i=0; i<10; i++){
117     hid.Form("htrkresx%i",i);
118     htrkresx[i] = dynamic_cast<TH1F*>(gDirectory->FindObject(hid.Data()));
119     hid.Form("htrkresy%i",i);
120     htrkresy[i] = dynamic_cast<TH1F*>(gDirectory->FindObject(hid.Data()));
121     }
122     for(Int_t ip=0; ip<6; ip++){
123     hid.Form("hxcalotrk%i",ip);
124     hxcalotrk[ip] = dynamic_cast<TH1F*>(gDirectory->FindObject(hid.Data()));
125     hid.Form("hycalotrk%i",ip);
126     hycalotrk[ip] = dynamic_cast<TH1F*>(gDirectory->FindObject(hid.Data()));
127     }
128     hresangx = dynamic_cast<TH1F*>(gDirectory->FindObject("hresangx"));
129     hresangy = dynamic_cast<TH1F*>(gDirectory->FindObject("hresangy"));
130    
131     // ----------------------------------
132     // calo variables
133     // ----------------------------------
134     CaloAxis *x_axis = new CaloAxis();
135     CaloAxis *y_axis = new CaloAxis();
136    
137     float rcil = 1.;// tolerance (cm)
138     x_axis->FitAxis(event->GetCaloLevel1(),0,rcil);
139     y_axis->FitAxis(event->GetCaloLevel1(),1,rcil);
140    
141     float xs31 = x_axis->par[0]+x_axis->par[1]*event->GetZTOF(31);
142     float xs32 = x_axis->par[0]+x_axis->par[1]*event->GetZTOF(32);
143     float ys31 = y_axis->par[0]+y_axis->par[1]*event->GetZTOF(31);
144     float ys32 = y_axis->par[0]+y_axis->par[1]*event->GetZTOF(32);
145     float qtrack = x_axis->GetQaxis()+y_axis->GetQaxis();
146     int ntrack = x_axis->GetN()+y_axis->GetN();
147    
148     CaloStrip top = CaloStrip();
149     top.Set(0,0,48);
150     float zcalotop = top.GetZ();
151    
152     float xcalotop = x_axis->par[0]+x_axis->par[1]*zcalotop;
153     float ycalotop = y_axis->par[0]+y_axis->par[1]*zcalotop;
154    
155    
156     float xtrk[6];
157     float ytrk[6];
158     for(Int_t ip=0; ip<6; ip++){
159     xtrk[ip] = x_axis->par[0]+x_axis->par[1]*event->GetZTrk(ip+1);
160     ytrk[ip] = y_axis->par[0]+y_axis->par[1]*event->GetZTrk(ip+1);
161     }
162    
163     float caloangx = (float)atan((double)x_axis->par[1])*180./3.1415026;
164     float caloangy = (float)atan((double)y_axis->par[1])*180./3.1415026;
165    
166     float xtrktop = x_axis->par[0]+x_axis->par[1]*event->GetZTrk(1);
167     float ytrktop = y_axis->par[0]+y_axis->par[1]*event->GetZTrk(1);
168    
169     x_axis->Delete();
170     y_axis->Delete();
171    
172     // tolleranza sul contenimento della traccia nel tracker
173     float dztrk = 44.51;
174     float dxtrk = 16.14;
175     float dytrk = 13.14;
176     float axtoll = 10.; //deg
177     float aytoll = 5.; //deg
178     float tollx = tan(atan(dxtrk/dztrk)+axtoll*3.1415026/180.)*dztrk-dxtrk;
179     float tolly = tan(atan(dytrk/dztrk)+aytoll*3.1415026/180.)*dztrk-dytrk;
180     // cout << tollx << " "<<tolly << endl;
181    
182     // ----------------------------------
183     // tof variables
184     // ----------------------------------
185     // get track stored by the tof
186     ToFTrkVar *tof = event->GetToFStoredTrack(-1);
187     if(!tof){
188     cout << " no ToF-track stored "<<endl;
189     return;
190     }
191     float resxs21 = (event->GetZTOF(21)-event->GetZTOF(32))*(tof->xtofpos[0]-tof->xtofpos[2])/(event->GetZTOF(12)-event->GetZTOF(32))+tof->xtofpos[2];
192     float resys22 = (event->GetZTOF(22)-event->GetZTOF(31))*(tof->ytofpos[0]-tof->ytofpos[2])/(event->GetZTOF(11)-event->GetZTOF(31))+tof->ytofpos[2];
193     resxs21 -= tof->xtofpos[1];
194     resys22 -= tof->ytofpos[1];
195    
196     float resxs32 = xs32 - tof->xtofpos[2];
197     float resys31 = ys31 - tof->ytofpos[2];
198    
199     // ----------------------------------
200     // tracker variables
201     // ----------------------------------
202     float resxtrk[]={-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.};
203     float resytrk[]={-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.};
204    
205     float resangx = -999;
206     float resangy = -999;
207     if( event->GetNTracks()==1 ){
208    
209     PamTrack *track = event->GetTrack(0);
210     if( track->chi2 >0 ){
211    
212     Float_t ztraj[13];
213    
214     Int_t i=0;
215     ztraj[i++] = event->GetZTOF(11);
216     ztraj[i++] = event->GetZTOF(12);
217     ztraj[i++] = event->GetZTOF(21);
218     ztraj[i++] = event->GetZTOF(22);
219     ztraj[i++] = event->GetZTrk(1);
220     ztraj[i++] = event->GetZTrk(2);
221     ztraj[i++] = event->GetZTrk(3);
222     ztraj[i++] = event->GetZTrk(4);
223     ztraj[i++] = event->GetZTrk(5);
224     ztraj[i++] = event->GetZTrk(6);
225     ztraj[i++] = event->GetZTOF(31);
226     ztraj[i++] = event->GetZTOF(32);
227     ztraj[i++] = zcalotop;
228     Trajectory traj = Trajectory(13,ztraj);
229     traj.DoTrack2( track->al );
230    
231     // --- S1
232     resxtrk[0] = track->xtofpos[0] - traj.x[1];
233     resytrk[0] = track->ytofpos[0] - traj.y[0];
234     // --- S2
235     resxtrk[1] = track->xtofpos[1] - traj.x[2];
236     resytrk[1] = track->ytofpos[1] - traj.y[3];
237     // --- S3
238     resxtrk[2] = track->xtofpos[2] - traj.x[11];
239     resytrk[2] = track->ytofpos[2] - traj.y[10];
240     // --- trk
241     for(Int_t ip=0; ip<6; ip++){
242     if(track->XGood(ip))resxtrk[3+ip] = track->xm[ip] - traj.x[4+ip];
243     if(track->YGood(ip))resytrk[3+ip] = track->ym[ip] - traj.y[4+ip];
244     }
245     // --- calo
246     resxtrk[9] = xcalotop - traj.x[12];
247     resytrk[9] = ycalotop - traj.y[12];
248    
249     resangx = caloangx - traj.thx[12];
250     resangy = caloangy - traj.thy[12];
251     // cout << caloangx << " -x- " << traj.thx[12]<<endl;
252     // cout << caloangy << " -y- " << traj.thy[12]<<endl;
253    
254     }
255     if(track) delete track;
256     }
257    
258     // ----------------------------------
259     // fill spatial residuals
260     // ----------------------------------
261     htofresx->Fill(resxs21);
262     htofresy->Fill(resys22);
263     hcalresx->Fill(resxs32);
264     hcalresy->Fill(resys31);
265     for(Int_t ih=0; ih<10; ih++){
266     htrkresx[ih]->Fill(resxtrk[ih]);
267     htrkresy[ih]->Fill(resytrk[ih]);
268     }
269     for(Int_t ip=0; ip<6; ip++){
270     hxcalotrk[ip]->Fill(xtrk[ip]);
271     hycalotrk[ip]->Fill(ytrk[ip]);
272     }
273     hresangx->Fill(resangx);
274     hresangy->Fill(resangy);
275    
276    
277     }
278     //===============================================================
279     // Save histograms
280     //===============================================================
281     void SaveHistos( TFile * outf ){
282    
283     TString hid;
284    
285     gROOT->cd();
286     // ----------------------------------
287     // retrieve histos
288     // ----------------------------------
289     htofresx = dynamic_cast<TH1F*>(gDirectory->FindObject("htofresx"));
290     htofresy = dynamic_cast<TH1F*>(gDirectory->FindObject("htofresy"));
291     // cout << hex << htofresx << dec << endl;
292     hcalresx = dynamic_cast<TH1F*>(gDirectory->FindObject("hcalresx"));
293     hcalresy = dynamic_cast<TH1F*>(gDirectory->FindObject("hcalresy"));
294     for(Int_t i=0; i<10; i++){
295     hid.Form("htrkresx%i",i);
296     htrkresx[i] = dynamic_cast<TH1F*>(gDirectory->FindObject(hid.Data()));
297     hid.Form("htrkresy%i",i);
298     htrkresy[i] = dynamic_cast<TH1F*>(gDirectory->FindObject(hid.Data()));
299     }
300     for(Int_t ip=0; ip<6; ip++){
301     hid.Form("hxcalotrk%i",ip);
302     hxcalotrk[ip] = dynamic_cast<TH1F*>(gDirectory->FindObject(hid.Data()));
303     hid.Form("hycalotrk%i",ip);
304     hycalotrk[ip] = dynamic_cast<TH1F*>(gDirectory->FindObject(hid.Data()));
305     }
306     hresangx = dynamic_cast<TH1F*>(gDirectory->FindObject("hresangx"));
307     hresangy = dynamic_cast<TH1F*>(gDirectory->FindObject("hresangy"));
308    
309     // ----------------------------------
310     // save on file
311     // ----------------------------------
312    
313     outf->cd();
314    
315     if(htofresx)htofresx->Write();
316     if(htofresy)htofresy->Write();
317     if(hcalresx)hcalresx->Write();
318     if(hcalresy)hcalresy->Write();
319    
320     for(Int_t i=0; i<10; i++){
321     if(htrkresx[i])htrkresx[i]->Write();
322     if(htrkresy[i])htrkresy[i]->Write();
323     }
324    
325     for(Int_t ip=0; ip<6; ip++){
326     if(hxcalotrk[ip])hxcalotrk[ip]->Write();
327     if(hycalotrk[ip])hycalotrk[ip]->Write();
328     }
329    
330     if(hresangx) hresangx->Write();
331     if(hresangy) hresangy->Write();
332     }
333    

  ViewVC Help
Powered by ViewVC 1.1.23