/[PAMELA software]/calo/ground/UTILITIES/macros/CaloTRKCALOALIG.c
ViewVC logotype

Contents of /calo/ground/UTILITIES/macros/CaloTRKCALOALIG.c

Parent Directory Parent Directory | Revision Log Revision Log

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