| 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 |
|
Parent Directory
| 
Revision Log