| 1 |
// |
| 2 |
//- Emiliano Mocchiutti |
| 3 |
// |
| 4 |
// Version 1.01 (2005-11-18) |
| 5 |
// |
| 6 |
// questo e` tutto da risistemare una volta che si ha l'output finale. |
| 7 |
// per ora ho sistemato tutto ma mai provato a compilare, si deve comunque cambiare l'apertura |
| 8 |
// del file (ora ne apre due...) attenzione alle strutture che non esistono piu`! |
| 9 |
// |
| 10 |
// Changelog: |
| 11 |
// |
| 12 |
// 1.00 - 1.01 (2005-11-18): Use pathtocalibration routine. |
| 13 |
// |
| 14 |
// 0.00 - 1.00 (2005-09-09): Working. |
| 15 |
// |
| 16 |
// 0.00 (2005-07-25): Created. |
| 17 |
// |
| 18 |
#include <fstream> |
| 19 |
#include <TTree.h> |
| 20 |
#include <TClassEdit.h> |
| 21 |
#include <TObject.h> |
| 22 |
#include <TList.h> |
| 23 |
#include <TSystem.h> |
| 24 |
#include <TSystemDirectory.h> |
| 25 |
#include <TString.h> |
| 26 |
#include <TFile.h> |
| 27 |
#include <TClass.h> |
| 28 |
#include <TCanvas.h> |
| 29 |
#include <TH1.h> |
| 30 |
#include <TH1F.h> |
| 31 |
#include <TH2D.h> |
| 32 |
#include <TLatex.h> |
| 33 |
#include <TPad.h> |
| 34 |
#include <TGraph.h> |
| 35 |
#include <TF1.h> |
| 36 |
#include <TStyle.h> |
| 37 |
// |
| 38 |
#include <string> |
| 39 |
extern const char *pathtocalibration(); |
| 40 |
extern void creadB(const char []); |
| 41 |
extern void ctrack(int , Double_t [], Double_t [], Double_t [], Double_t [], int &); |
| 42 |
// |
| 43 |
extern void stringcopy(TString&, const TString&, Int_t, Int_t); |
| 44 |
extern bool existfile(TString); |
| 45 |
// |
| 46 |
// |
| 47 |
void FCaloTrackerAlignmentRT(TString filename){ |
| 48 |
if ( !existfile(filename) ){ |
| 49 |
printf(" %s :no such file or directory \n",filename.Data()); |
| 50 |
return; |
| 51 |
}; |
| 52 |
const char* startingdir = gSystem->WorkingDirectory(); |
| 53 |
TString path; |
| 54 |
stringcopy(path,startingdir); |
| 55 |
// |
| 56 |
// |
| 57 |
Long64_t nevents = 0; |
| 58 |
Long64_t trnevents = 0; |
| 59 |
// |
| 60 |
// Read the magnetic field map (two files), here it goes the path AND name of one of them. |
| 61 |
// |
| 62 |
printf("\n\n TRACKER: loading the magnetic field maps...\n\n\n"); |
| 63 |
const char *pam_calib = pathtocalibration(); |
| 64 |
stringstream bdir; |
| 65 |
bdir.str(""); |
| 66 |
bdir << pam_calib << "."; |
| 67 |
creadB((char *)bdir.str().c_str()); |
| 68 |
// |
| 69 |
printf(" ...done! \n"); |
| 70 |
// |
| 71 |
Float_t rig; |
| 72 |
// |
| 73 |
Float_t sigmacalo = 10.; |
| 74 |
// |
| 75 |
Float_t piano[22]; |
| 76 |
piano[0] = 0.; |
| 77 |
for ( Int_t ee = 1; ee<22; ee++){ |
| 78 |
if ( ee%2 ){ |
| 79 |
piano[ee] = piano[ee-1] - 8.09; |
| 80 |
} else { |
| 81 |
piano[ee] = piano[ee-1] - 10.09; |
| 82 |
}; |
| 83 |
}; |
| 84 |
// |
| 85 |
Float_t zoff = 0.; // 260 |
| 86 |
// |
| 87 |
Float_t chi2 = 0.; |
| 88 |
Float_t chi2nev = 0.; |
| 89 |
// |
| 90 |
Float_t xalign0 = -120.; |
| 91 |
Float_t yalign0 = -120.; |
| 92 |
Float_t zalign0 = -260.; // 0 |
| 93 |
Float_t phi0 = 0.; |
| 94 |
Float_t the0 = 0.; |
| 95 |
Float_t psi0 = 0.; |
| 96 |
Float_t dxalig = 100.; //50 |
| 97 |
Float_t dyalig = 100.; //50 |
| 98 |
Float_t dzalig = 140.; //80 |
| 99 |
// |
| 100 |
Float_t dpsi = 0.20; |
| 101 |
Float_t dthe = 0.20; |
| 102 |
Float_t dphi = 0.20; |
| 103 |
// |
| 104 |
TString lfile; |
| 105 |
TString sfile; |
| 106 |
TString file2f; |
| 107 |
TString pfile; |
| 108 |
TString fififile; |
| 109 |
// |
| 110 |
Float_t minxalign; |
| 111 |
Float_t minyalign; |
| 112 |
Float_t minzalign; |
| 113 |
Float_t minpsialign; |
| 114 |
Float_t minthealign; |
| 115 |
Float_t minphialign; |
| 116 |
Float_t minxalign0; |
| 117 |
Float_t minyalign0; |
| 118 |
Float_t minzalign0; |
| 119 |
Float_t minpsi0; |
| 120 |
Float_t minthe0; |
| 121 |
Float_t minphi0; |
| 122 |
Int_t minnxa; |
| 123 |
Int_t minnza; |
| 124 |
Int_t minnya; |
| 125 |
Int_t minphi; |
| 126 |
Int_t minpsi; |
| 127 |
Int_t minthe; |
| 128 |
// |
| 129 |
const int npia = 0; |
| 130 |
TFile *hfile = 0; |
| 131 |
// |
| 132 |
// |
| 133 |
const string myfil = (const char *)filename; |
| 134 |
Int_t myposiz = myfil.find("dw_"); |
| 135 |
if ( myposiz == -1 ) { |
| 136 |
myposiz = myfil.find("DW_"); |
| 137 |
}; |
| 138 |
TString fileid; |
| 139 |
stringcopy(fileid,filename,myposiz+3,myposiz+13); |
| 140 |
// |
| 141 |
// |
| 142 |
stringstream namefile; |
| 143 |
TTree *tree; |
| 144 |
const char *pattho = path; |
| 145 |
namefile.str(""); |
| 146 |
namefile << pattho << "/CaloTrkAlignData_"; |
| 147 |
namefile << fileid.Data() << ".root"; |
| 148 |
printf(" Save file %s\n",namefile.str().c_str()); |
| 149 |
hfile = new TFile(namefile.str().c_str(),"RECREATE","Calorimeter-Tracker alignment data"); |
| 150 |
tree = new TTree("CaloTrkAlign","Calorimeter-Tracker alignment data"); |
| 151 |
tree->Branch("chi2",&chi2,"chi2/F"); |
| 152 |
tree->Branch("chi2nev",&chi2nev,"chi2nev/F"); |
| 153 |
tree->Branch("minnxa",&minnxa,"minnxa/I"); |
| 154 |
tree->Branch("minnya",&minnya,"minnya/I"); |
| 155 |
tree->Branch("minnza",&minnza,"minnza/I"); |
| 156 |
tree->Branch("minpsi",&minpsi,"minpsi/I"); |
| 157 |
tree->Branch("minphi",&minphi,"minphi/I"); |
| 158 |
tree->Branch("minthe",&minthe,"minthe/I"); |
| 159 |
tree->Branch("minxalign",&minxalign,"minxalign/F"); |
| 160 |
tree->Branch("minyalign",&minyalign,"minyalign/F"); |
| 161 |
tree->Branch("minzalign",&minzalign,"minzalign/F"); |
| 162 |
tree->Branch("minpsialign",&minpsialign,"minpsialign/F"); |
| 163 |
tree->Branch("minphialign",&minphialign,"minphialign/F"); |
| 164 |
tree->Branch("minthealign",&minthealign,"minthealign/F"); |
| 165 |
// |
| 166 |
const int ngrid = 11; |
| 167 |
Float_t togl = ((float)ngrid-1.)/2.; |
| 168 |
// |
| 169 |
struct Tracklev2 trk; |
| 170 |
struct CaLevel2 calo; |
| 171 |
// |
| 172 |
Int_t i = -1; |
| 173 |
Int_t itr = -1; |
| 174 |
Int_t syncro = 1; |
| 175 |
Int_t pktnum = 0; |
| 176 |
Int_t obt = 0; |
| 177 |
// |
| 178 |
// open files |
| 179 |
// |
| 180 |
TFile *caloFile = 0; |
| 181 |
TFile *trackerFile = 0; |
| 182 |
TTree *tr = 0; |
| 183 |
TTree *otr = 0; |
| 184 |
// |
| 185 |
caloFile = emigetFile(filename, "Physics.Level2", "Calorimeter"); |
| 186 |
if ( !caloFile ){ |
| 187 |
printf("\n ERROR: no calorimeter file! \n\n"); |
| 188 |
return; |
| 189 |
}; |
| 190 |
trackerFile = emigetFile(filename, "Physics.Level2", "Tracker"); |
| 191 |
if ( !trackerFile ){ |
| 192 |
printf("\n ERROR: no tracker file! \n\n"); |
| 193 |
caloFile->Close(); |
| 194 |
return; |
| 195 |
} |
| 196 |
// |
| 197 |
//Takes the tree of the header file |
| 198 |
// |
| 199 |
tr = (TTree*)trackerFile->Get("TrkLevel2"); |
| 200 |
settrklev2(tr,trk); |
| 201 |
// |
| 202 |
otr = (TTree*)caloFile->Get("CaloLevel2"); |
| 203 |
setcalolevel2(otr,calo); |
| 204 |
// |
| 205 |
trnevents = tr->GetEntries(); |
| 206 |
nevents = otr->GetEntries(); |
| 207 |
const int lung = nevents; |
| 208 |
Int_t goodev[lung]; |
| 209 |
Float_t mxtb[11][lung]; |
| 210 |
Float_t mytb[11][lung]; |
| 211 |
// |
| 212 |
i = -1; |
| 213 |
itr = -1; |
| 214 |
syncro = 1; |
| 215 |
pktnum = 0; |
| 216 |
obt = 0; |
| 217 |
printf(" Initializing... \n"); |
| 218 |
Int_t npoint = 44; |
| 219 |
Int_t goodevno = 0; |
| 220 |
while ( i < (nevents-1) ){ |
| 221 |
// |
| 222 |
if ( i%1000 == 0 && i > 0 ) printf(" %iK \n",i/1000); |
| 223 |
// |
| 224 |
// look for tracker data |
| 225 |
// |
| 226 |
itr++; |
| 227 |
i++; |
| 228 |
syncro = 1; |
| 229 |
// |
| 230 |
trkcalosync0: |
| 231 |
// |
| 232 |
// check if we have tracker data |
| 233 |
// |
| 234 |
if ( i >= nevents ){ |
| 235 |
printf(" WARNING: no more calorimeter data.\n"); |
| 236 |
goto skipthisev0; |
| 237 |
}; |
| 238 |
if ( itr > trnevents ){ |
| 239 |
printf(" WARNING: no more tracker data.\n"); |
| 240 |
goto skipthisev0; |
| 241 |
}; |
| 242 |
// |
| 243 |
// retrieve tracker informations |
| 244 |
// |
| 245 |
tr->GetEntry(itr); |
| 246 |
// |
| 247 |
// check synchronization tracker and calorimeter informations: |
| 248 |
// |
| 249 |
otr->GetEntry(i); |
| 250 |
// |
| 251 |
pktnum = calo.pkt_num; |
| 252 |
obt = calo.OBT; |
| 253 |
if ( pktnum != trk.pkt_num || obt != trk.obt ){ |
| 254 |
if ( pktnum > trk.pkt_num || obt > trk.obt ){ |
| 255 |
itr++; |
| 256 |
syncro = 0; |
| 257 |
goto trkcalosync0; |
| 258 |
}; |
| 259 |
if ( pktnum < trk.pkt_num || obt < trk.obt ){ |
| 260 |
i++; |
| 261 |
syncro = 0; |
| 262 |
goto trkcalosync0; |
| 263 |
}; |
| 264 |
}; |
| 265 |
// |
| 266 |
// here we have synchronized data |
| 267 |
// |
| 268 |
if ( !syncro ) { |
| 269 |
syncro = 1; |
| 270 |
}; |
| 271 |
// |
| 272 |
// here we go! |
| 273 |
// |
| 274 |
goodev[i] = 0; |
| 275 |
// |
| 276 |
// select events with one track... |
| 277 |
// |
| 278 |
if ( trk.ntrk == 1 && trk.good2 ){ |
| 279 |
// |
| 280 |
// a good chi2 from tracker and rigidity greater than 1 GV (notice correction in deflection) |
| 281 |
// |
| 282 |
rig = 1./(trk.al[0][4]); |
| 283 |
if ( trk.chi2[0] < 25 && (rig > 0.5 || rig < -0.5) ){ |
| 284 |
// |
| 285 |
// a good fit of the track from calorimeter information |
| 286 |
// |
| 287 |
if ( calo.npcfit[0]>15 && calo.npcfit[1]>15 && calo.varcfit[0]<20. && calo.varcfit[1]<20. && calo.good == 1 ){ //era 19! |
| 288 |
// |
| 289 |
// determine the position of the track according to tracker: |
| 290 |
// |
| 291 |
Double_t zin[22]; |
| 292 |
Double_t xtb[22]; |
| 293 |
Double_t ytb[22]; |
| 294 |
Double_t al_p[5]; |
| 295 |
Int_t ifail; |
| 296 |
Int_t ee; |
| 297 |
Int_t eb; |
| 298 |
// |
| 299 |
minzalign0 = zalign0; |
| 300 |
minzalign = 0.; |
| 301 |
// |
| 302 |
ee = 0; |
| 303 |
eb = 11; |
| 304 |
// |
| 305 |
for ( Int_t ea = 0; ea < 11; ea++) { |
| 306 |
// |
| 307 |
eb--; |
| 308 |
minzalign = minzalign0 + ((float)eb - togl)*(dzalig/((float)ngrid+1.)); |
| 309 |
// |
| 310 |
if ( ea < 5 ) al_p[ea] = (Double_t)trk.al[0][ea]; |
| 311 |
xtb[ee] = 0.; |
| 312 |
ytb[ee] = 0.; |
| 313 |
zin[ee] = (piano[npia] + minzalign)/10.; |
| 314 |
ee++; |
| 315 |
xtb[ee] = 0.; |
| 316 |
ytb[ee] = 0.; |
| 317 |
zin[ee] = (piano[npia] - 5.1 + minzalign)/10.; |
| 318 |
ee++; |
| 319 |
}; |
| 320 |
// |
| 321 |
npoint = 22; |
| 322 |
ifail = 0; |
| 323 |
// |
| 324 |
ctrack(npoint,zin,xtb,ytb,al_p,ifail); |
| 325 |
// |
| 326 |
// here alignment factors for x and y |
| 327 |
// |
| 328 |
if ( !ifail ){ |
| 329 |
goodev[i] = 1; |
| 330 |
goodevno++; |
| 331 |
// |
| 332 |
minzalign0 = zalign0; |
| 333 |
minzalign = 0.; |
| 334 |
ee = 21; |
| 335 |
for ( Int_t ea = 0; ea < 11 ; ea++) { |
| 336 |
mxtb[ea][i] = (float)xtb[ee] * 10.; |
| 337 |
ee--; |
| 338 |
mytb[ea][i] = (float)ytb[ee] * 10.; |
| 339 |
ee--; |
| 340 |
}; |
| 341 |
// |
| 342 |
}; |
| 343 |
}; |
| 344 |
}; |
| 345 |
}; |
| 346 |
// |
| 347 |
skipthisev0: |
| 348 |
continue; |
| 349 |
}; |
| 350 |
printf(" ...done! \n"); |
| 351 |
printf(" %i good events \n",goodevno); |
| 352 |
// |
| 353 |
Int_t ee = 0; |
| 354 |
Float_t the; |
| 355 |
Float_t psi; |
| 356 |
Float_t phi; |
| 357 |
// |
| 358 |
Float_t xcbar; |
| 359 |
Float_t ycbar; |
| 360 |
Float_t xtbar; |
| 361 |
Float_t ytbar; |
| 362 |
// |
| 363 |
// for (Int_t it = 0; it < 3; it++){ |
| 364 |
for (Int_t it = 0; it < 1; it++){ |
| 365 |
// printf(" iteraction %i \n",it); |
| 366 |
// |
| 367 |
if ( it == 0 ){ |
| 368 |
minxalign0 = xalign0; |
| 369 |
minyalign0 = yalign0; |
| 370 |
minzalign0 = zalign0; |
| 371 |
minpsi0 = psi0; |
| 372 |
minthe0 = the0; |
| 373 |
minphi0 = phi0; |
| 374 |
} else { |
| 375 |
minxalign0 = minxalign; |
| 376 |
minyalign0 = minyalign; |
| 377 |
minzalign0 = minzalign; |
| 378 |
minpsi0 = minpsialign; |
| 379 |
minthe0 = minthealign; |
| 380 |
minphi0 = minphialign; |
| 381 |
}; |
| 382 |
// |
| 383 |
for (Int_t nthe = 0; nthe<ngrid; nthe++){ |
| 384 |
//for (Int_t nthe = 5; nthe<6; nthe++){ |
| 385 |
//for (Int_t nthe = 4; nthe<7; nthe++){ |
| 386 |
minthealign = minthe0 + ((float)nthe - togl)*(dthe/((float)ngrid+1.)); |
| 387 |
// |
| 388 |
for (Int_t nphi = 0; nphi<ngrid; nphi++){ |
| 389 |
//for (Int_t nphi = 4; nphi<5; nphi++){ |
| 390 |
//for (Int_t nphi = 5; nphi<6; nphi++){ |
| 391 |
//for (Int_t nphi = 4; nphi<7; nphi++){ |
| 392 |
minphialign = minphi0 + ((float)nphi - togl)*(dphi/((float)ngrid+1.)); |
| 393 |
// |
| 394 |
for (Int_t npsi = 0; npsi<ngrid; npsi++){ |
| 395 |
//for (Int_t npsi = 5; npsi<6; npsi++){ |
| 396 |
//for (Int_t npsi = 4; npsi<7; npsi++){ |
| 397 |
minpsialign = minpsi0 + ((float)npsi - togl)*(dpsi/((float)ngrid+1.)); |
| 398 |
// |
| 399 |
for (Int_t nxa = 0; nxa<ngrid; nxa++){ |
| 400 |
//for (Int_t nxa = 4; nxa<5; nxa++){ |
| 401 |
// printf(" nxa %i \n",nxa); |
| 402 |
minxalign = minxalign0 + ((float)nxa - togl)*(dxalig/((float)ngrid+1.)); |
| 403 |
//printf("minxalign[%i][%i] = %f\n",g,h,minxalign[g][h]); |
| 404 |
//for (Int_t nya = 6; nya<7; nya++){ |
| 405 |
for (Int_t nya = 0; nya<ngrid; nya++){ |
| 406 |
// printf(" nya %i \n",nya); |
| 407 |
minyalign = minyalign0 + ((float)nya - togl)*(dyalig/((float)ngrid+1.)); |
| 408 |
//printf("minyalign[%i][%i] = %f\n",g,h,minyalign[g][h]); |
| 409 |
//for (Int_t nza = 5; nza<6; nza++){ |
| 410 |
//for (Int_t nza = ngrid ; nza>=0; nza--){ |
| 411 |
for (Int_t nza = 0 ; nza < ngrid; nza++){ |
| 412 |
// |
| 413 |
minzalign = minzalign0 + ((float)nza - togl)*(dzalig/((float)ngrid+1.)); |
| 414 |
printf(" Step %i xa %f ya %f za %f phi %f the %f psi %f \n",it,minxalign,minyalign,minzalign,minphialign,minthealign,minpsialign); |
| 415 |
// |
| 416 |
i = -1; |
| 417 |
itr = -1; |
| 418 |
syncro = 1; |
| 419 |
pktnum = 0; |
| 420 |
obt = 0; |
| 421 |
// |
| 422 |
while ( i < (nevents-1) ){ |
| 423 |
// |
| 424 |
// look for tracker data |
| 425 |
// |
| 426 |
Int_t ea = npia; |
| 427 |
itr++; |
| 428 |
i++; |
| 429 |
syncro = 1; |
| 430 |
if ( goodev[i] == 0 ) { |
| 431 |
goto skipthisev; |
| 432 |
}; |
| 433 |
ee = 0; |
| 434 |
the = minthealign; |
| 435 |
psi = minpsialign; |
| 436 |
phi = minphialign; |
| 437 |
trkcalosync: |
| 438 |
// |
| 439 |
// check if we have tracker data |
| 440 |
// |
| 441 |
if ( i >= nevents ){ |
| 442 |
printf(" WARNING: no more calorimeter data.\n"); |
| 443 |
goto closeandexit; |
| 444 |
}; |
| 445 |
if ( itr > trnevents ){ |
| 446 |
printf(" WARNING: no more tracker data.\n"); |
| 447 |
goto closeandexit; |
| 448 |
}; |
| 449 |
// |
| 450 |
// retrieve tracker informations |
| 451 |
// |
| 452 |
tr->GetEntry(itr); |
| 453 |
// |
| 454 |
// check synchronization tracker and calorimeter informations: |
| 455 |
// |
| 456 |
otr->GetEntry(i); |
| 457 |
// |
| 458 |
pktnum = calo.pkt_num; |
| 459 |
obt = calo.OBT; |
| 460 |
if ( pktnum != trk.pkt_num || obt != trk.obt ){ |
| 461 |
if ( pktnum > trk.pkt_num || obt > trk.obt ){ |
| 462 |
itr++; |
| 463 |
syncro = 0; |
| 464 |
goto trkcalosync; |
| 465 |
}; |
| 466 |
if ( pktnum < trk.pkt_num || obt < trk.obt ){ |
| 467 |
i++; |
| 468 |
syncro = 0; |
| 469 |
goto trkcalosync; |
| 470 |
}; |
| 471 |
}; |
| 472 |
// |
| 473 |
// here we have synchronized data |
| 474 |
// |
| 475 |
if ( !syncro ) { |
| 476 |
syncro = 1; |
| 477 |
}; |
| 478 |
// |
| 479 |
if ( goodev[i] == 0 ) { |
| 480 |
goto skipthisev; |
| 481 |
}; |
| 482 |
ea = npia; |
| 483 |
// |
| 484 |
xcbar = (calo.cbar[ea][0] + minxalign) * cos(phi) * cos(the) + (calo.cbar[ea][1] + minyalign) * (cos(phi) * sin(the) * sin(psi) - sin(phi) * cos(psi)) + (piano[ea] - zoff + minzalign) * (cos(phi) * sin(the) * cos(psi) + sin(phi) * sin (psi)); |
| 485 |
ycbar = (calo.cbar[ea][0] + minxalign) * sin(phi) * cos(the) + (calo.cbar[ea][1] + minyalign) * (sin(phi) * sin(the) * sin(phi) + cos(the) * cos(psi) ) + (piano[ea] - 5.1 - zoff + minzalign) * (sin(phi) * sin(the) * cos(psi) - cos(the) * sin(psi)); |
| 486 |
// |
| 487 |
xtbar = mxtb[nza][i]; |
| 488 |
ytbar = -mytb[nza][i]; |
| 489 |
// |
| 490 |
// |
| 491 |
chi2nev += 1.; |
| 492 |
// |
| 493 |
chi2 += ((xcbar-xtbar)*(xcbar-xtbar)+(ycbar-ytbar)*(ycbar-ytbar))/(sigmacalo*sigmacalo); |
| 494 |
// |
| 495 |
skipthisev: |
| 496 |
continue; |
| 497 |
}; |
| 498 |
closeandexit: |
| 499 |
// |
| 500 |
minnxa = nxa; |
| 501 |
minnya = nya; |
| 502 |
minnza = nza; |
| 503 |
minphi = nphi; |
| 504 |
minpsi = npsi; |
| 505 |
minthe = nthe; |
| 506 |
tree->Fill(); |
| 507 |
chi2nev = 0.; |
| 508 |
chi2 = 0.; |
| 509 |
}; |
| 510 |
}; |
| 511 |
}; |
| 512 |
}; |
| 513 |
}; |
| 514 |
}; |
| 515 |
// |
| 516 |
// run over different alignment parameters |
| 517 |
// |
| 518 |
printf("\n\n"); |
| 519 |
dxalig -= 7.; |
| 520 |
dyalig -= 7.; |
| 521 |
dzalig -= 7.; |
| 522 |
dphi -= 0.02; |
| 523 |
dpsi -= 0.02; |
| 524 |
dthe -= 0.02; |
| 525 |
}; |
| 526 |
hfile->Write(); |
| 527 |
hfile->Close(); |
| 528 |
}; |
| 529 |
|
| 530 |
|
| 531 |
void Ffindminimum(TString filename){ |
| 532 |
const char* startingdir = gSystem->WorkingDirectory(); |
| 533 |
TString path; |
| 534 |
stringcopy(path,startingdir); |
| 535 |
Int_t minnxa; |
| 536 |
Int_t minnza; |
| 537 |
Int_t minnya; |
| 538 |
Int_t minphi; |
| 539 |
Int_t minpsi; |
| 540 |
Int_t minthe; |
| 541 |
Float_t minxalign; |
| 542 |
Float_t minyalign; |
| 543 |
Float_t minzalign; |
| 544 |
Float_t minpsialign; |
| 545 |
Float_t minthealign; |
| 546 |
Float_t minphialign; |
| 547 |
Int_t nxa = 0; |
| 548 |
Int_t nza = 0; |
| 549 |
Int_t nya = 0; |
| 550 |
Int_t phi = 0; |
| 551 |
Int_t psi = 0; |
| 552 |
Int_t the = 0; |
| 553 |
Float_t chi2; |
| 554 |
Float_t chi2nev; |
| 555 |
TFile *f = new TFile(filename); |
| 556 |
TTree *tr = (TTree*) f->Get("CaloTrkAlign"); |
| 557 |
tr->SetBranchAddress("chi2",&chi2); |
| 558 |
tr->SetBranchAddress("chi2nev",&chi2nev); |
| 559 |
tr->SetBranchAddress("minnxa",&minnxa); |
| 560 |
tr->SetBranchAddress("minnya",&minnya); |
| 561 |
tr->SetBranchAddress("minnza",&minnza); |
| 562 |
tr->SetBranchAddress("minphi",&minphi); |
| 563 |
tr->SetBranchAddress("minpsi",&minpsi); |
| 564 |
tr->SetBranchAddress("minthe",&minthe); |
| 565 |
tr->SetBranchAddress("minxalign",&minxalign); |
| 566 |
tr->SetBranchAddress("minyalign",&minyalign); |
| 567 |
tr->SetBranchAddress("minzalign",&minzalign); |
| 568 |
tr->SetBranchAddress("minpsialign",&minpsialign); |
| 569 |
tr->SetBranchAddress("minphialign",&minphialign); |
| 570 |
tr->SetBranchAddress("minthealign",&minthealign); |
| 571 |
// |
| 572 |
Long64_t nevents = tr->GetEntries(); |
| 573 |
// Float_t min = 1000000000.; |
| 574 |
Float_t min2 = 1000000000.; |
| 575 |
for (Int_t i = nevents-1; i>=0; i--){ |
| 576 |
tr->GetEntry(i); |
| 577 |
// if ( chi2/chi2nev <= min ){ |
| 578 |
// min = chi2/chi2nev; |
| 579 |
// printf(" New minimum found:\n chi2 = %f minnxa = %i minnya = %i minnza = %i minthe = %i minpsi = %i minphi = %i \n\n",chi2/chi2nev,minnxa,minnya,minnza,minthe,minpsi,minphi); |
| 580 |
// }; |
| 581 |
// if ( minnxa==4 && minnya==6 && minnza==6 && chi2/chi2nev <= min2 ){ |
| 582 |
if ( minnxa<6 && minnxa>2 && minnya<8 && minnya>4 && chi2/chi2nev <= min2 ){ |
| 583 |
min2 = chi2/chi2nev; |
| 584 |
// printf(" 2- New minimum found:\n chi2 = %f minnxa = %i minnya = %i minnza = %i minthe = %i minpsi = %i minphi = %i \n\n",chi2/chi2nev,minnxa,minnya,minnza,minthe,minpsi,minphi); |
| 585 |
nxa = minnxa; |
| 586 |
nya = minnya; |
| 587 |
nza = minnza; |
| 588 |
the = minthe; |
| 589 |
phi = minphi; |
| 590 |
psi = minpsi; |
| 591 |
}; |
| 592 |
}; |
| 593 |
printf("\n The best minimum found is for chi2 = %f with:\n minnxa = %i minnya = %i minnza = %i minthe = %i minpsi = %i minphi = %i \n\n",min2,nxa,nya,nza,the,psi,phi); |
| 594 |
Int_t nps = 0; |
| 595 |
Int_t nph = 0; |
| 596 |
Int_t nth = 0; |
| 597 |
Int_t nx = 0; |
| 598 |
Int_t ny = 0; |
| 599 |
Int_t nz = 0; |
| 600 |
Float_t vnpsi[11]; |
| 601 |
Float_t vnphi[11]; |
| 602 |
Float_t vnthe[11]; |
| 603 |
Float_t vnx[11]; |
| 604 |
Float_t vny[11]; |
| 605 |
Float_t vnz[11]; |
| 606 |
Float_t vrpsi[11]; |
| 607 |
Float_t vrphi[11]; |
| 608 |
Float_t vrthe[11]; |
| 609 |
Float_t vrx[11]; |
| 610 |
Float_t vry[11]; |
| 611 |
Float_t vrz[11]; |
| 612 |
for (Int_t i = nevents-1; i>=0; i--){ |
| 613 |
tr->GetEntry(i); |
| 614 |
if ( minnxa==nxa && minnya==nya && minnza==nza && minthe==the && minphi==phi ){ |
| 615 |
vnpsi[nps] = chi2/chi2nev; |
| 616 |
vrpsi[nps] = minpsialign; |
| 617 |
nps++; |
| 618 |
}; |
| 619 |
if ( minnxa==nxa && minnya==nya && minnza==nza && minthe==the && minpsi==psi ){ |
| 620 |
vnphi[nph] = chi2/chi2nev; |
| 621 |
vrphi[nph] = minphialign; |
| 622 |
nph++; |
| 623 |
}; |
| 624 |
if ( minnxa==nxa && minnya==nya && minnza==nza && minphi==phi && minpsi==psi ){ |
| 625 |
vnthe[nth] = chi2/chi2nev; |
| 626 |
vrthe[nth] = minthealign; |
| 627 |
nth++; |
| 628 |
}; |
| 629 |
if ( minnxa==nxa && minnya==nya && minthe==the && minphi==phi && minpsi==psi ){ |
| 630 |
vnz[nz] = chi2/chi2nev; |
| 631 |
vrz[nz] = minzalign; |
| 632 |
nz++; |
| 633 |
}; |
| 634 |
if ( minnxa==nxa && minnza==nza && minthe==the && minphi==phi && minpsi==psi ){ |
| 635 |
vny[ny] = chi2/chi2nev; |
| 636 |
vry[ny] = -minyalign; |
| 637 |
ny++; |
| 638 |
}; |
| 639 |
if ( minnya==nya && minnza==nza && minthe==the && minphi==phi && minpsi==psi ){ |
| 640 |
vnx[nx] = chi2/chi2nev; |
| 641 |
vrx[nx] = -minxalign; |
| 642 |
nx++; |
| 643 |
}; |
| 644 |
}; |
| 645 |
// |
| 646 |
// |
| 647 |
// |
| 648 |
const string myfil = (const char *)filename; |
| 649 |
Int_t myposiz = myfil.find("dw_"); |
| 650 |
if ( myposiz == -1 ) { |
| 651 |
myposiz = myfil.find("DW_"); |
| 652 |
}; |
| 653 |
TString fileid; |
| 654 |
stringcopy(fileid,filename,myposiz+18,myposiz+28); |
| 655 |
// |
| 656 |
// |
| 657 |
stringstream namefile; |
| 658 |
const char *pattho = path; |
| 659 |
// |
| 660 |
TCanvas *figura1 = new TCanvas("figura1","figura1", 1100, 900); |
| 661 |
TGraph *gr; |
| 662 |
TF1 *fit; |
| 663 |
Float_t par2; |
| 664 |
Float_t par3; |
| 665 |
stringstream titolo; |
| 666 |
// |
| 667 |
figura1->Clear(); |
| 668 |
figura1->cd(); |
| 669 |
gPad->SetTickx(); |
| 670 |
gPad->SetTicky(); |
| 671 |
gr = new TGraph(11,vrx,vnx); |
| 672 |
gr->SetLineWidth(1); |
| 673 |
gr->SetMarkerStyle(20); |
| 674 |
gStyle->SetOptFit(0); |
| 675 |
gr->GetXaxis()->SetTitle("x offset [mm]"); |
| 676 |
gr->GetYaxis()->SetTitle("#chi^{2}"); |
| 677 |
gr->Fit("pol2","q"); |
| 678 |
fit = gr->GetFunction("pol2"); |
| 679 |
par2 = (float)fit->GetParameter(1); |
| 680 |
par3 = (float)fit->GetParameter(2); |
| 681 |
titolo.str(""); |
| 682 |
titolo << " Minimum at x = " << -par2/(2.*par3); |
| 683 |
titolo << " mm"; |
| 684 |
gr->SetTitle(titolo.str().c_str()); |
| 685 |
gr->Draw("AP"); |
| 686 |
namefile.str(""); |
| 687 |
namefile << pattho << "/x_chi2_"; |
| 688 |
namefile << fileid.Data() << ".ps"; |
| 689 |
printf(" Save figure %s\n",namefile.str().c_str()); |
| 690 |
figura1->SaveAs(namefile.str().c_str()); |
| 691 |
// |
| 692 |
figura1->Clear(); |
| 693 |
figura1->cd(); |
| 694 |
gPad->SetTickx(); |
| 695 |
gPad->SetTicky(); |
| 696 |
gr = new TGraph(11,vry,vny); |
| 697 |
gr->SetLineWidth(1); |
| 698 |
gr->SetMarkerStyle(20); |
| 699 |
gStyle->SetOptFit(0); |
| 700 |
gr->GetXaxis()->SetTitle("y offset [mm]"); |
| 701 |
gr->GetYaxis()->SetTitle("#chi^{2}"); |
| 702 |
gr->Fit("pol2","q"); |
| 703 |
fit = gr->GetFunction("pol2"); |
| 704 |
par2 = (float)fit->GetParameter(1); |
| 705 |
par3 = (float)fit->GetParameter(2); |
| 706 |
titolo.str(""); |
| 707 |
titolo << " Minimum at y = " << -par2/(2.*par3); |
| 708 |
titolo << " mm"; |
| 709 |
gr->SetTitle(titolo.str().c_str()); |
| 710 |
gr->Draw("AP"); |
| 711 |
namefile.str(""); |
| 712 |
namefile << pattho << "/y_chi2_"; |
| 713 |
namefile << fileid.Data() << ".ps"; |
| 714 |
printf(" Save figure %s\n",namefile.str().c_str()); |
| 715 |
figura1->SaveAs(namefile.str().c_str()); |
| 716 |
// |
| 717 |
figura1->Clear(); |
| 718 |
figura1->cd(); |
| 719 |
gPad->SetTickx(); |
| 720 |
gPad->SetTicky(); |
| 721 |
gr = new TGraph(11,vrz,vnz); |
| 722 |
gr->SetLineWidth(1); |
| 723 |
gr->SetMarkerStyle(20); |
| 724 |
gStyle->SetOptFit(0); |
| 725 |
gr->GetXaxis()->SetTitle("z offset [mm]"); |
| 726 |
gr->GetYaxis()->SetTitle("#chi^{2}"); |
| 727 |
gr->Fit("pol2","q"); |
| 728 |
fit = gr->GetFunction("pol2"); |
| 729 |
par2 = (float)fit->GetParameter(1); |
| 730 |
par3 = (float)fit->GetParameter(2); |
| 731 |
titolo.str(""); |
| 732 |
titolo << " Minimum at z = " << -par2/(2.*par3); |
| 733 |
titolo << " mm"; |
| 734 |
gr->SetTitle(titolo.str().c_str()); |
| 735 |
gr->Draw("AP"); |
| 736 |
namefile.str(""); |
| 737 |
namefile << pattho << "/z_chi2_"; |
| 738 |
namefile << fileid.Data() << ".ps"; |
| 739 |
printf(" Save figure %s\n",namefile.str().c_str()); |
| 740 |
figura1->SaveAs(namefile.str().c_str()); |
| 741 |
// |
| 742 |
figura1->Clear(); |
| 743 |
figura1->cd(); |
| 744 |
gPad->SetTickx(); |
| 745 |
gPad->SetTicky(); |
| 746 |
gr = new TGraph(11,vrphi,vnphi); |
| 747 |
gr->SetLineWidth(1); |
| 748 |
gr->SetMarkerStyle(20); |
| 749 |
gStyle->SetOptFit(0); |
| 750 |
gr->GetXaxis()->SetTitle("#phi_{z} offset [rad]"); |
| 751 |
gr->GetYaxis()->SetTitle("#chi^{2}"); |
| 752 |
gr->Fit("pol2","q"); |
| 753 |
fit = gr->GetFunction("pol2"); |
| 754 |
par2 = (float)fit->GetParameter(1); |
| 755 |
par3 = (float)fit->GetParameter(2); |
| 756 |
titolo.str(""); |
| 757 |
titolo << " Minimum at #phi_{z} = " << -par2/(2.*par3); |
| 758 |
titolo << " rad"; |
| 759 |
gr->SetTitle(titolo.str().c_str()); |
| 760 |
gr->Draw("AP"); |
| 761 |
namefile.str(""); |
| 762 |
namefile << pattho << "/phi_chi2_"; |
| 763 |
namefile << fileid.Data() << ".ps"; |
| 764 |
printf(" Save figure %s\n",namefile.str().c_str()); |
| 765 |
figura1->SaveAs(namefile.str().c_str()); |
| 766 |
// |
| 767 |
figura1->Clear(); |
| 768 |
figura1->cd(); |
| 769 |
gPad->SetTickx(); |
| 770 |
gPad->SetTicky(); |
| 771 |
gr = new TGraph(11,vrthe,vnthe); |
| 772 |
gr->SetLineWidth(1); |
| 773 |
gr->SetMarkerStyle(20); |
| 774 |
gStyle->SetOptFit(0); |
| 775 |
gr->GetXaxis()->SetTitle("#theta_{y} offset [rad]"); |
| 776 |
gr->GetYaxis()->SetTitle("#chi^{2}"); |
| 777 |
gr->Fit("pol2","q"); |
| 778 |
fit = gr->GetFunction("pol2"); |
| 779 |
par2 = (float)fit->GetParameter(1); |
| 780 |
par3 = (float)fit->GetParameter(2); |
| 781 |
titolo.str(""); |
| 782 |
titolo << " Minimum at #theta_{y} = " << -par2/(2.*par3); |
| 783 |
titolo << " rad"; |
| 784 |
gr->SetTitle(titolo.str().c_str()); |
| 785 |
gr->Draw("AP"); |
| 786 |
namefile.str(""); |
| 787 |
namefile << pattho << "/the_chi2_"; |
| 788 |
namefile << fileid.Data() << ".ps"; |
| 789 |
printf(" Save figure %s\n",namefile.str().c_str()); |
| 790 |
figura1->SaveAs(namefile.str().c_str()); |
| 791 |
// |
| 792 |
figura1->Clear(); |
| 793 |
figura1->cd(); |
| 794 |
gPad->SetTickx(); |
| 795 |
gPad->SetTicky(); |
| 796 |
gr = new TGraph(11,vrpsi,vnpsi); |
| 797 |
gr->SetLineWidth(1); |
| 798 |
gr->SetMarkerStyle(20); |
| 799 |
gStyle->SetOptFit(0); |
| 800 |
gr->GetXaxis()->SetTitle("#psi_{x} offset [rad]"); |
| 801 |
gr->GetYaxis()->SetTitle("#chi^{2}"); |
| 802 |
gr->Fit("pol2","q"); |
| 803 |
fit = gr->GetFunction("pol2"); |
| 804 |
par2 = (float)fit->GetParameter(1); |
| 805 |
par3 = (float)fit->GetParameter(2); |
| 806 |
titolo.str(""); |
| 807 |
titolo << " Minimum at #psi_{x} = " << -par2/(2.*par3); |
| 808 |
titolo << " rad"; |
| 809 |
gr->SetTitle(titolo.str().c_str()); |
| 810 |
gr->Draw("AP"); |
| 811 |
namefile.str(""); |
| 812 |
namefile << pattho << "/psi_chi2_"; |
| 813 |
namefile << fileid.Data() << ".ps"; |
| 814 |
printf(" Save figure %s\n",namefile.str().c_str()); |
| 815 |
figura1->SaveAs(namefile.str().c_str()); |
| 816 |
// |
| 817 |
} |
| 818 |
|
| 819 |
void Ffindminimumnorotation(TString filename){ |
| 820 |
Int_t minnxa; |
| 821 |
Int_t minnza; |
| 822 |
Int_t minnya; |
| 823 |
Int_t minphi; |
| 824 |
Int_t minpsi; |
| 825 |
Int_t minthe; |
| 826 |
Float_t chi2; |
| 827 |
Float_t chi2nev; |
| 828 |
TFile *f = new TFile(filename); |
| 829 |
TTree *tr = (TTree*) f->Get("CaloTrkAlign"); |
| 830 |
tr->SetBranchAddress("chi2",&chi2); |
| 831 |
tr->SetBranchAddress("chi2nev",&chi2nev); |
| 832 |
tr->SetBranchAddress("minnxa",&minnxa); |
| 833 |
tr->SetBranchAddress("minnya",&minnya); |
| 834 |
tr->SetBranchAddress("minnza",&minnza); |
| 835 |
tr->SetBranchAddress("minphi",&minphi); |
| 836 |
tr->SetBranchAddress("minpsi",&minpsi); |
| 837 |
tr->SetBranchAddress("minthe",&minthe); |
| 838 |
Long64_t nevents = tr->GetEntries(); |
| 839 |
Float_t min = 1000000000.; |
| 840 |
for (Int_t i = nevents-1; i>=0; i--){ |
| 841 |
tr->GetEntry(i); |
| 842 |
if ( chi2/chi2nev <= min && minthe==5 && minpsi==5 && minphi == 5 ){ |
| 843 |
min = chi2/chi2nev; |
| 844 |
printf(" New minimum found:\n chi2 = %f minnxa = %i minnya = %i minnza = %i minthe = %i minpsi = %i minphi = %i \n\n",chi2/chi2nev,minnxa,minnya,minnza,minthe,minpsi,minphi); |
| 845 |
}; |
| 846 |
}; |
| 847 |
|
| 848 |
|
| 849 |
} |
Parent Directory
| 
Revision Log