/[PAMELA software]/calo/flight/FUTILITIES/macros/FCaloTRKCALOALIG.cxx
ViewVC logotype

Annotation of /calo/flight/FUTILITIES/macros/FCaloTRKCALOALIG.cxx

Parent Directory Parent Directory | Revision Log Revision Log

Revision 1.1 - (hide annotations) (download)
Wed Mar 22 15:07:46 2006 UTC (18 years, 8 months ago) by mocchiut
Branch point for: FUTILITIES, MAIN 
Initial revision
1 mocchiut 1.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     }
  ViewVC Help
Powered by ViewVC 1.1.23