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

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 1.1.1.1 - (hide annotations) (download) (vendor branch)
Tue Dec 6 10:59:28 2005 UTC (20 years ago) by mocchiut
Branch: MAIN, UTILITIES
CVS Tags: start, v2r01, v2r00, HEAD
Changes since 1.1: +0 -0 lines
File MIME type: text/plain 
Imported sources
1 mocchiut 1.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