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

Contents of /calo/flight/FUTILITIES/macros/FCaloFUNCTIONS.cxx

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1.1.1 - (show annotations) (download) (vendor branch)
Wed Mar 22 15:07:47 2006 UTC (18 years, 8 months ago) by mocchiut
Branch: FUTILITIES
CVS Tags: start, v1r01
Changes since 1.1: +0 -0 lines
Flight Utilities calorimeter package 1st release

1 //
2 // Calorimeter useful functions and subroutines - Emiliano Mocchiutti
3 //
4 // FCaloFUNCTIONS.cxx version 1.00 (2006-03-03)
5 //
6 // Programs in this file are called by other programs and cannot be run by hand.
7 //
8 // Changelog:
9 //
10 // 1.00 - 1.01 (2006-03-03): read unique YODA file and compile correctly with ROOT classes etc. Small changes in the fetchpreviousfile/whatnamewith routines.
11 //
12 // 0.00 - 1.00 (2006-03-03) Clone of CaloFunctions.h
13 //
14 #include <iostream>
15 #include <sstream>
16 #include <fstream>
17 #include <TSystem.h>
18 #include <TH1.h>
19 #include <TF1.h>
20 #include <TROOT.h>
21 #include <TStyle.h>
22 #include <TClassEdit.h>
23 #include <TObject.h>
24 #include <TList.h>
25 #include <TSystem.h>
26 #include <TSystemDirectory.h>
27 #include <TString.h>
28 #include <TTree.h>
29 #include <TChain.h>
30 #include <TFile.h>
31 #include <TClass.h>
32 #include <TCanvas.h>
33 #include <TKey.h>
34 #include <TClassTable.h>
35 #include <TThread.h>
36 //
37 #include <PamelaRun.h>
38 #include <physics/calorimeter/CalorimeterEvent.h>
39 #include <CalibCalPedEvent.h>
40 //
41 #include <fcalostructs.h>
42 #include <FCaloFUNCTIONSfun.h>
43 #include <caloclassesfun.h>
44 //
45 using namespace std;
46
47 void stringcopy(TString& s1, const TString& s2, Int_t from=0, Int_t to=0){
48 if ( to == 0 ){
49 Int_t t2length = s2.Length();
50 s1 = "";
51 to = t2length;
52 };
53 for (Int_t i = from; i<to; i++){
54 s1.Append(s2[i],1);
55 };
56 }
57
58 void stringappend(TString& s1, const TString& s2){
59 Int_t t2length = s2.Length();
60 for (Int_t i = 0; i<t2length; i++){
61 s1.Append(s2[i],1);
62 };
63 }
64
65 void CaloCompressData(Calib & calib, Evento & evento, Int_t l, Int_t m, Int_t pre){
66 Float_t minstrip = 100000.;
67 Float_t rms = 0.;
68 for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
69 if ( calib.calgood[l][m][e] == 0. && evento.dexyc[l][m][e]-calib.calped[l][m][e] < minstrip && evento.dexyc[l][m][e] > 0.) {
70 minstrip = evento.dexyc[l][m][e]-calib.calped[l][m][e];
71 rms = calib.calthr[l][m][pre];
72 };
73 };
74 //
75 // compression
76 //
77 if ( minstrip < evento.base[l][m][pre] && minstrip != 100000.){
78 for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
79 if ( evento.dexyc[l][m][e]-calib.calped[l][m][e] <= (minstrip+rms) ) evento.dexyc[l][m][e] = 0.;
80 };
81 };
82 }
83
84 void CaloFindBase(Calib & calib, Evento & evento, Int_t l, Int_t m, Int_t pre){
85 Float_t ominstrip = 100000.;
86 Float_t orms = 0.;
87 Float_t minstrip = 100000.;
88 Float_t rms = 0.;
89 evento.base[l][m][pre] = 0.;
90 for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
91 if ( calib.calgood[l][m][e] == 0. && evento.dexyc[l][m][e]-calib.calped[l][m][e] < minstrip && evento.dexyc[l][m][e] > 0.) {
92 ominstrip = minstrip;
93 minstrip = evento.dexyc[l][m][e]-calib.calped[l][m][e];
94 orms = rms;
95 rms = calib.calthr[l][m][pre];
96 };
97 if ( calib.calgood[l][m][e] == 0. && evento.dexyc[l][m][e]-calib.calped[l][m][e] < ominstrip && evento.dexyc[l][m][e]-calib.calped[l][m][e] >minstrip ) {
98 ominstrip = evento.dexyc[l][m][e]-calib.calped[l][m][e];
99 orms = calib.calthr[l][m][pre];
100 };
101 };
102 if ( minstrip != 100000. ) {
103 if ( ominstrip != 10000.) {
104 minstrip = ominstrip;
105 rms = orms;
106 };
107 Float_t strip6s = 0.;
108 for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
109 if ( (evento.dexyc[l][m][e]-calib.calped[l][m][e]) >= minstrip && (evento.dexyc[l][m][e]-calib.calped[l][m][e]) <= (minstrip+rms) ) {
110 strip6s += 1.;
111 evento.base[l][m][pre] += (evento.dexyc[l][m][e] - calib.calped[l][m][e]);
112 };
113 //
114 // compression
115 //
116 if ( evento.dexyc[l][m][e]-calib.calped[l][m][e] <= (minstrip+rms) ) evento.dexyc[l][m][e] = 0.;
117 };
118 if ( strip6s >= 9. ){
119 Double_t arro = evento.base[l][m][pre]/strip6s;
120 Float_t deci = 1000.*((float)arro - float(int(arro)));
121 if ( deci < 500. ) {
122 arro = double(int(arro));
123 } else {
124 arro = 1. + double(int(arro));
125 };
126 evento.base[l][m][pre] = arro;
127 } else {
128 evento.base[l][m][pre] = 31000.;
129 };
130 } else {
131 evento.base[l][m][pre] = 31000.;
132 };
133 }
134
135 void CaloFindBaseRaw(Calib & calib, Evento & evento, Int_t l, Int_t m, Int_t pre){
136 Float_t minstrip = 100000.;
137 Float_t rms = 0.;
138 evento.base[l][m][pre] = 0.;
139 for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
140 if ( calib.calgood[l][m][e] == 0. && evento.dexy[l][m][e]-calib.calped[l][m][e] < minstrip && evento.dexy[l][m][e] > 0.) {
141 minstrip = evento.dexy[l][m][e]-calib.calped[l][m][e];
142 rms = calib.calthr[l][m][pre];
143 };
144 };
145 if ( minstrip != 100000. ) {
146 Float_t strip6s = 0.;
147 for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
148 if ( (evento.dexy[l][m][e]-calib.calped[l][m][e]) >= minstrip && (evento.dexy[l][m][e]-calib.calped[l][m][e]) <= (minstrip+rms) ) {
149 strip6s += 1.;
150 evento.base[l][m][pre] += (evento.dexy[l][m][e] - calib.calped[l][m][e]);
151 };
152 //
153 // compression
154 //
155 if ( abs((int)(evento.dexy[l][m][e]-calib.calped[l][m][e])) <= (minstrip+rms) ) {
156 evento.dexyc[l][m][e] = 0.;
157 } else {
158 evento.dexyc[l][m][e] = evento.dexy[l][m][e];
159 };
160 };
161 if ( strip6s >= 9. ){
162 Double_t arro = evento.base[l][m][pre]/strip6s;
163 Float_t deci = 1000.*((float)arro - float(int(arro)));
164 if ( deci < 500. ) {
165 arro = double(int(arro));
166 } else {
167 arro = 1. + double(int(arro));
168 };
169 evento.base[l][m][pre] = arro;
170 } else {
171 evento.base[l][m][pre] = 31000.;
172 for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
173 evento.dexyc[l][m][e] = evento.dexy[l][m][e];
174 };
175 };
176 } else {
177 evento.base[l][m][pre] = 31000.;
178 };
179 }
180
181 void CaloFindBaseRawNC(Calib & calib, Evento & evento, Int_t l, Int_t m, Int_t pre){
182 Float_t minstrip = 100000.;
183 Float_t rms = 0.;
184 evento.base[l][m][pre] = 0.;
185 for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
186 if ( calib.calgood[l][m][e] == 0. && evento.dexy[l][m][e]-calib.calped[l][m][e] < minstrip && evento.dexy[l][m][e] > 0.) {
187 minstrip = evento.dexy[l][m][e]-calib.calped[l][m][e];
188 rms = calib.calthr[l][m][pre];
189 };
190 };
191 if ( minstrip != 100000. ) {
192 Float_t strip6s = 0.;
193 for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
194 if ( (evento.dexy[l][m][e]-calib.calped[l][m][e]) >= minstrip && (evento.dexy[l][m][e]-calib.calped[l][m][e]) <= (minstrip+rms) ) {
195 strip6s += 1.;
196 evento.base[l][m][pre] += (evento.dexy[l][m][e] - calib.calped[l][m][e]);
197 };
198 //
199 // compression
200 //
201 evento.dexyc[l][m][e] = evento.dexy[l][m][e];
202 };
203 if ( strip6s >= 9. ){
204 Double_t arro = evento.base[l][m][pre]/strip6s;
205 Float_t deci = 1000.*((float)arro - float(int(arro)));
206 if ( deci < 500. ) {
207 arro = double(int(arro));
208 } else {
209 arro = 1. + double(int(arro));
210 };
211 evento.base[l][m][pre] = arro;
212 } else {
213 evento.base[l][m][pre] = 31000.;
214 for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
215 evento.dexyc[l][m][e] = evento.dexy[l][m][e];
216 };
217 };
218 } else {
219 evento.base[l][m][pre] = 31000.;
220 };
221 }
222
223 Bool_t existfile(TString filename){
224 ifstream myfile;
225 myfile.open(filename.Data());
226 if ( !myfile ){
227 // printf("\n\nERROR: not able to open file: %s \n",filename.Data());
228 return(false);
229 };
230 myfile.close();
231 return(true);
232 }
233
234 int CaloPede(TString filename, Int_t s, Int_t atime, Calib & calib){
235 if ( !existfile(filename) ) return(3);
236 //
237 TFile *File = new TFile(filename.Data());
238 TTree *tr = (TTree*)File->Get("CalibCalPed");
239 //
240 pamela::CalibCalPedEvent *ce = 0;
241 pamela::PscuHeader *cph = 0;
242 pamela::EventHeader *ceh = 0;
243 tr->SetBranchAddress("Header", &ceh);
244 tr->SetBranchAddress("CalibCalPed", &ce);
245 //
246 Long64_t ncalibs = tr->GetEntries();
247 for (Int_t ci = 0; ci < ncalibs ; ci++){
248 tr->GetEntry(ci);
249 cph = ceh->GetPscuHeader();
250 if ( atime == cph->GetOrbitalTime()){
251 calib.iev = ce->iev;
252 if (ce->cstwerr[s] != 0 && ce->cperror[s] == 0 ) {
253 for ( Int_t d=0 ; d<11 ;d++ ){
254 Int_t pre = -1;
255 for ( Int_t j=0; j<96 ;j++){
256 if ( j%16 == 0 ) pre++;
257 if ( s == 2 ){
258 calib.calped[0][2*d+1][j] = ce->calped[3][d][j];
259 calib.cstwerr[3] = ce->cstwerr[3];
260 calib.cperror[3] = ce->cperror[3];
261 calib.calgood[0][2*d+1][j] = ce->calgood[3][d][j];
262 calib.calthr[0][2*d+1][pre] = ce->calthr[3][d][pre];
263 calib.calrms[0][2*d+1][j] = ce->calrms[3][d][j];
264 calib.calbase[0][2*d+1][pre] = ce->calbase[3][d][pre];
265 calib.calvar[0][2*d+1][pre] = ce->calvar[3][d][pre];
266 };
267 if ( s == 3 ){
268 calib.calped[0][2*d][j] = ce->calped[1][d][j];
269 calib.cstwerr[1] = ce->cstwerr[1];
270 calib.cperror[1] = ce->cperror[1];
271 calib.calgood[0][2*d][j] = ce->calgood[1][d][j];
272 calib.calthr[0][2*d][pre] = ce->calthr[1][d][pre];
273 calib.calrms[0][2*d][j] = ce->calrms[1][d][j];
274 calib.calbase[0][2*d][pre] = ce->calbase[1][d][pre];
275 calib.calvar[0][2*d][pre] = ce->calvar[1][d][pre];
276 };
277 if ( s == 0 ){
278 calib.calped[1][2*d][j] = ce->calped[0][d][j];
279 calib.cstwerr[0] = ce->cstwerr[0];
280 calib.cperror[0] = ce->cperror[0];
281 calib.calgood[1][2*d][j] = ce->calgood[0][d][j];
282 calib.calthr[1][2*d][pre] = ce->calthr[0][d][pre];
283 calib.calrms[1][2*d][j] = ce->calrms[0][d][j];
284 calib.calbase[1][2*d][pre] = ce->calbase[0][d][pre];
285 calib.calvar[1][2*d][pre] = ce->calvar[0][d][pre];
286 };
287 if ( s == 1 ){
288 calib.calped[1][2*d+1][j] = ce->calped[2][d][j];
289 calib.cstwerr[2] = ce->cstwerr[2];
290 calib.cperror[2] = ce->cperror[2];
291 calib.calgood[1][2*d+1][j] = ce->calgood[2][d][j];
292 calib.calthr[1][2*d+1][pre] = ce->calthr[2][d][pre];
293 calib.calrms[1][2*d+1][j] = ce->calrms[2][d][j];
294 calib.calbase[1][2*d+1][pre] = ce->calbase[2][d][pre];
295 calib.calvar[1][2*d+1][pre] = ce->calvar[2][d][pre];
296 };
297 };
298 };
299 };
300 };
301 };
302 tr->Delete();
303 return(0);
304 }
305
306 TString getFilename(const TString filename){
307 const string fil = (const char*)filename;
308 Int_t posiz = fil.find("dw_");
309 if ( posiz == -1 ) posiz = fil.find("DW_");
310 if ( posiz == -1 ) return 0;
311 Int_t posiz2 = posiz+13;
312 TString file2;
313 stringcopy(file2,filename,posiz,posiz2);
314 TString pdat(".dat");
315 stringappend(file2,pdat);
316 return file2;
317 }
318
319 Int_t fetchpreviousfile(TString ffile, Int_t s){
320 const char *file = ffile;
321 stringstream fil;
322 fil.str("");
323 fil << file ;
324 Int_t posiz = fil.str().find("dw_");
325 Int_t upper = 0;
326 if ( posiz == -1 ) {
327 posiz = fil.str().find("DW_");
328 upper = 1;
329 };
330 if ( posiz == -1 ) return(1);
331 //
332 TString nomefile = getFilename(ffile);
333 const char *ufnome = nomefile;
334 stringstream nyear;
335 stringstream nmonth;
336 stringstream nday;
337 stringstream nfno;
338 TString year;
339 TString month;
340 TString day;
341 TString fno;
342 TString yno;
343 stringcopy(year, ufnome, 3, 5);
344 stringcopy(month, ufnome, 5, 7);
345 stringcopy(day, ufnome, 7, 9);
346 stringcopy(fno, ufnome, 10, 13);
347 stringcopy(yno, ffile.Data(), posiz+14, posiz+16);
348 const char *cyear = year;
349 const char *cmonth = month;
350 const char *cday = day;
351 const char *cfno = fno;
352 const char *cyno = yno;
353 Int_t jy = atoi(cyear);
354 Int_t jm = atoi(cmonth);
355 Int_t jd = atoi(cday);
356 Int_t jn = atoi(cfno);
357 Int_t jynst = atoi(cyno);
358 Int_t jyn = 0;
359 Int_t inter = 0;
360 char *zy;
361 char *zm;
362 char *zd;
363 char *zn = "";
364 char *zyn;
365 Int_t goodthis = 0;
366 Int_t onegood = 0;
367 Long64_t ncalibs;
368 if ( jn > 1 ) {
369 jn--;
370 } else {
371 jn = 999;
372 if ( jd > 1 ) {
373 jd--;
374 } else {
375 if ( jm > 1 ) {
376 jm--;
377 } else {
378 jm = 12;
379 if ( jy > 3 ) {
380 jy--;
381 } else {
382 return(1);
383 };
384 };
385 if ( jm==4 || jm==6 || jm==9 || jm==11 ){
386 jd = 30;
387 } else {
388 jd = 31;
389 };
390 if ( jm==2 ) jd = 29;
391 };
392 };
393 Int_t gday = 0;
394 Int_t retvar = 1;
395 stringstream nyd;
396 while ( jy>3 ) {
397 while ( jm>0 ){
398 while ( jd > 0 ){
399 while ( jn > 0 ){
400 retvar--;
401 nyear.str("");
402 nyear << jy;
403 if ( jy<10 ){
404 zy = "0";
405 } else {
406 zy = "";
407 };
408 nmonth.str("");
409 nmonth << jm;
410 if ( jm<10 ){
411 zm = "0";
412 } else {
413 zm = "";
414 };
415 nday.str("");
416 nday << jd;
417 if ( jd<10 ){
418 zd = "0";
419 } else {
420 zd = "";
421 };
422 nfno.str("");
423 nfno << jn;
424 if ( jn>99 ) zn = "";
425 if ( jn<100 && jn>9 ) zn = "0";
426 if ( jn<10 ) zn = "00";
427 jyn = jynst;
428 goodthis = 0;
429 onegood = 0;
430 while ( jyn < 100 ){
431 begwhile:
432 nyd.str("");
433 nyd << jyn;
434 if ( jyn<10 ){
435 zyn = "0";
436 } else {
437 zyn = "";
438 };
439 stringstream nfile;
440 nfile.str("");
441 if ( upper ){
442 nfile << "DW_" << zy;
443 } else {
444 nfile << "dw_" << zy;
445 };
446 nfile << nyear.str().c_str();
447 nfile << zm;
448 nfile << nmonth.str().c_str();
449 nfile << zd;
450 nfile << nday.str().c_str() << "_";
451 nfile << zn;
452 nfile << nfno.str().c_str();
453 nfile << zyn;
454 nfile << nyd.str().c_str();
455 nfile << ".root"; // addded
456 //
457 TString fgfile = "";
458 stringcopy(fgfile,ffile,0,posiz);
459 stringappend(fgfile,nfile.str().c_str());
460 //
461 TFile *File = 0;
462 TTree *tr = 0;
463 pamela::CalibCalPedEvent *ce = 0;
464 pamela::PscuHeader *cph = 0;
465 pamela::EventHeader *ceh = 0;
466 if ( !existfile(fgfile) ){
467 if ( jyn > jynst && onegood ) {
468 goodthis = 1;
469 jyn--;
470 goto begwhile;
471 };
472 if ( jyn == jynst ) {
473 onegood = 0;
474 jyn = 100;
475 goto salta;
476 };
477 if ( jyn > jynst && !onegood ) {
478 onegood = 0;
479 jyn = 100;
480 goto salta;
481 };
482 goto salta;
483 } else {
484 onegood = 1;
485 if ( !goodthis ) {
486 goto salta;
487 };
488 };
489 File = new TFile(fgfile.Data());
490 tr = (TTree*)File->Get("CalibCalPed");
491 //Takes the tree of the header file
492 tr->SetBranchAddress("Header", &ceh);
493 tr->SetBranchAddress("CalibCalPed", &ce);
494 ncalibs = tr->GetEntries();
495 if ( ncalibs == 0 ) {
496 jyn = 100;
497 onegood = 0;
498 goodthis = 0;
499 File->Close();
500 goto salta;
501 };
502 inter = 0;
503 for (Int_t c = 0; c < ncalibs; c++){
504 tr->GetEntry(c);
505 cph = ceh->GetPscuHeader();
506 if ( ce->cstwerr[s] != 0 && ce->cperror[s] == 0 ) {
507 inter++;
508 };
509 };
510 if ( inter ){
511 jyn = 100;
512 File->Close();
513 return(retvar);
514 };
515 File->Close();
516 salta:
517 jyn++;
518 };
519 jn--;
520 };
521 jn = 999;
522 jd--;
523 gday++;
524 if ( gday>1 ) return(1);
525 };
526 jm--;
527 if ( jm==4 || jm==6 || jm==9 || jm==11 ){
528 jd = 30;
529 } else {
530 jd = 31;
531 };
532 if ( jm==2 ) jd = 29;
533 };
534 jm = 12;
535 jy--;
536 };
537 return(retvar);
538 }
539
540 TString whatnamewith(TString ffile, Int_t retvar){
541 Int_t debug = 0;
542 const char *file = ffile;
543 stringstream fil;
544 fil.str("");
545 fil << file ;
546 if ( debug ) printf("ffile %s retvar %i \n",file,retvar);
547 Int_t posiz = fil.str().find("dw_");
548 Int_t upper = 0;
549 if ( posiz == -1 ) {
550 posiz = fil.str().find("DW_");
551 upper = 1;
552 };
553 if ( posiz == -1 ) return("wrong");
554 //
555 TString nomefile = getFilename(ffile);
556 const char *ufnome = nomefile;
557 stringstream nyear;
558 stringstream nmonth;
559 stringstream nday;
560 stringstream nfno;
561 TString year;
562 TString month;
563 TString day;
564 TString fno;
565 TString yno;
566 stringcopy(year, ufnome, 3, 5);
567 stringcopy(month, ufnome, 5, 7);
568 stringcopy(day, ufnome, 7, 9);
569 stringcopy(fno, ufnome, 10, 13);
570 stringcopy(yno, ffile.Data(), posiz+14, posiz+16);
571 const char *cyear = year;
572 const char *cmonth = month;
573 const char *cday = day;
574 const char *cfno = fno;
575 const char *cyno = yno;
576 Int_t jy = atoi(cyear);
577 Int_t jm = atoi(cmonth);
578 Int_t jd = atoi(cday);
579 Int_t jn = atoi(cfno);
580 Int_t jynst = atoi(cyno);
581 Int_t jyn = 0;
582 Int_t goodthis = 0;
583 Int_t onegood = 0;
584 char *zy;
585 char *zm;
586 char *zd;
587 char *zn = "";
588 char *zyn;
589 Long64_t ncalibs;
590 TString fgfile;
591 if ( jn > 1 ) {
592 jn--;
593 } else {
594 jn = 999;
595 if ( jd > 1 ) {
596 jd--;
597 } else {
598 if ( jm > 1 ) {
599 jm--;
600 } else {
601 jm = 12;
602 if ( jy > 3 ) {
603 jy--;
604 } else {
605 return("wrong");
606 };
607 };
608 if ( jm==4 || jm==6 || jm==9 || jm==11 ){
609 jd = 30;
610 } else {
611 jd = 31;
612 };
613 if ( jm==2 ) jd = 29;
614 };
615 };
616 Int_t retval = 1;
617 Int_t gday = 0;
618 stringstream nyd;
619 while ( jy>3 ) {
620 while ( jm>0 ){
621 while ( jd > 0 ){
622 while ( jn > 0 ){
623 retval--;
624 if ( retval < retvar ) return("wrong");
625 nyear.str("");
626 nyear << jy;
627 if ( jy<10 ){
628 zy = "0";
629 } else {
630 zy = "";
631 };
632 nmonth.str("");
633 nmonth << jm;
634 if ( jm<10 ){
635 zm = "0";
636 } else {
637 zm = "";
638 };
639 nday.str("");
640 nday << jd;
641 if ( jd<10 ){
642 zd = "0";
643 } else {
644 zd = "";
645 };
646 nfno.str("");
647 nfno << jn;
648 if ( jn>99 ) zn = "";
649 if ( jn<100 && jn>9 ) zn = "0";
650 if ( jn<10 ) zn = "00";
651 jyn = jynst;
652 goodthis = 0;
653 onegood = 0;
654 while ( jyn < 100 ){
655 begwhile:
656 nyd.str("");
657 nyd << jyn;
658 if ( jyn<10 ){
659 zyn = "0";
660 } else {
661 zyn = "";
662 };
663 stringstream nfile;
664 nfile.str("");
665 if ( upper ){
666 nfile << "DW_" << zy;
667 } else {
668 nfile << "dw_" << zy;
669 };
670 nfile << nyear.str().c_str();
671 nfile << zm;
672 nfile << nmonth.str().c_str();
673 nfile << zd;
674 nfile << nday.str().c_str() << "_";
675 nfile << zn;
676 nfile << nfno.str().c_str();
677 nfile << zyn;
678 nfile << nyd.str().c_str();
679 nfile << ".root"; // added
680 //
681 fgfile = "";
682 stringcopy(fgfile,ffile,0,posiz);
683 stringappend(fgfile,nfile.str().c_str());
684 const char *tosee = fgfile;
685 if ( debug ) printf("file: %s retval %i retvar %i \n",tosee,retval,retvar);
686 TFile *File = 0;
687 TTree *tr = 0;
688 //
689 if ( !existfile(fgfile) ){
690 if ( jyn > jynst && onegood ) {
691 goodthis = 1;
692 jyn--;
693 goto begwhile;
694 };
695 if ( jyn == jynst ) {
696 onegood = 0;
697 jyn = 100;
698 goto salta;
699 };
700 if ( jyn > jynst && !onegood ) {
701 onegood = 0;
702 jyn = 100;
703 goto salta;
704 };
705 goto salta;
706 } else {
707 onegood = 1;
708 if ( !goodthis ){
709 goto salta;
710 };
711 };
712 File = new TFile(fgfile.Data());
713 tr = (TTree*)File->Get("CalibCalPed");
714 ncalibs = tr->GetEntries();
715 if ( ncalibs == 0 ) {
716 jyn = 100;
717 onegood = 0;
718 goodthis = 0;
719 File->Close();
720 goto salta;
721 };
722 File->Close();
723 if ( retval == retvar ) return(fgfile);
724 salta:
725 jyn++;
726 };
727 jn--;
728 };
729 jn = 999;
730 jd--;
731 gday++;
732 if ( gday>7 ) return("wrong");
733 };
734 jm--;
735 if ( jm==4 || jm==6 || jm==9 || jm==11 ){
736 jd = 30;
737 } else {
738 jd = 31;
739 };
740 if ( jm==2 ) jd = 29;
741 };
742 jm = 12;
743 jy--;
744 };
745 return("wrong");
746 }
747
748 char *getLEVname(TString filename, TString detc, TString numb, TString frame="root"){
749 // char *file;
750 stringstream file;
751 file.str("");
752 const char *num = numb;
753 const char *det = detc;
754 const char *fra = frame;
755 string fil = (const char *)filename;
756 Int_t posiz = fil.find("dw_");
757 if ( posiz == -1 ) posiz = fil.find("DW_");
758 if ( posiz == -1 ) return(0);
759 Int_t spos = posiz;
760 Int_t epos = posiz+15;
761 TString tmptempf;
762 stringcopy(tmptempf,filename,spos,epos);
763 const char *tempf = tmptempf;
764 file << tempf << ".Physics.Level";
765 file << num << ".";
766 file << det << ".Event.";
767 file << fra;
768 const char *rfile = file.str().c_str();
769 return((char*)rfile);
770 };
771
772 char *getfileLEVname(TString filename, TString detc, TString numb, TString frame="root"){
773 // char *file;
774 stringstream file;
775 const char *num = numb;
776 const char *det = detc;
777 const char *fra = frame;
778 string fil = (const char *)filename;
779 Int_t posiz = fil.find("dw_");
780 if ( posiz == -1 ) posiz = fil.find("DW_");
781 if ( posiz == -1 ) return(0);
782 Int_t spos = posiz;
783 Int_t epos = posiz+13;
784 TString tmptempf;
785 stringcopy(tmptempf,filename,spos,epos);
786 const char *tempf = tmptempf;
787 file.str("");
788 file << tempf << "00.Physics.Level";
789 // file << "00.Physics.Level";
790 file << num << ".";
791 file << det << ".Event.";
792 file << fra;
793 const char *rfile = file.str().c_str();
794 return((char*)rfile);
795 };
796
797 int OLDCaloFindCalibs(TString &filename, Calib & calib){
798 for (Int_t s = 0; s < 4; s++){
799 for (Int_t d = 1; d<50; d++){
800 calib.ttime[s][d] = 0;
801 if ( d < 49 ) calib.time[s][d] = 0;
802 };
803 };
804 if ( !existfile(filename) ) return(1);
805 //
806 TFile *File = new TFile(filename.Data());
807 TTree *tr = (TTree*)File->Get("CalibCalPed");
808 pamela::CalibCalPedEvent *ce = 0;
809 pamela::PscuHeader *cph = 0;
810 pamela::EventHeader *ceh = 0;
811 tr->SetBranchAddress("Header", &ceh);
812 tr->SetBranchAddress("CalibCalPed", &ce);
813 Long64_t ncalibs = tr->GetEntries();
814 Int_t inter;
815 for (Int_t s = 0; s < 4; s++){
816 for (Int_t d = 1; d<50; d++){
817 calib.ttime[s][d] = 0;
818 if ( d < 49 ) calib.time[s][d] = 0;
819 };
820 };
821 for (Int_t s = 0; s < 4; s++){
822 inter = 0;
823 for (Int_t c = 0; c < ncalibs; c++){
824 tr->GetEntry(c);
825 cph = ceh->GetPscuHeader();
826 calib.ttime[s][inter] = 0;
827 if ( ce->cstwerr[s] != 0 && ce->cperror[s] == 0 ) {
828 calib.ttime[s][inter] = cph->GetOrbitalTime();
829 inter++;
830 } else {
831 if ( ce->cstwerr[s] != 0 && ce->cperror[s] != 0 ) {
832 printf(" ERROR: entry %i stwerr %X perror %f \n",c,ce->cstwerr[s],ce->cperror[s]);
833 };
834 };
835 };
836 if ( inter == 0 ){
837 printf(" ERROR: no suitable calibration for section %i in this file!\n",s);
838 };
839 for (Int_t d = 1; d<50; d++){
840 if ( calib.ttime[s][d] != 0 ) {
841 calib.time[s][d-1] = calib.ttime[s][d-1] + (int)((calib.ttime[s][d] - calib.ttime[s][d-1])/2.);
842 } else {
843 if ( d == 1 ) {
844 calib.time[s][d-1] = 0;
845 };
846 };
847 };
848 };
849 File->Close();
850 return(0);
851 }
852
853 int CaloFindCalibs(TString &filename, TString &calibfile, Int_t &wused, Calib & calib){
854 //
855 // First of all check if the file has a monotone growing OBT time function.
856 //
857 Int_t debug = 0;
858 Int_t obtjump = 0;
859 Int_t evjump[50];
860 Int_t OBT = 0;
861 Int_t OOBT = 0;
862 Int_t OBT1 = 0;
863 Int_t LOBT = 0;
864 pamela::PscuHeader *ph = 0;
865 pamela::EventHeader *eh = 0;
866 //
867 if ( !existfile(filename) ) return(1);
868 TFile *File = new TFile(filename.Data());
869 TTree *ctr = (TTree*)File->Get("Physics");
870 ctr->SetBranchStatus("*",0);
871 ctr->SetBranchStatus("Header*",1);
872 ctr->SetBranchStatus("Pscu*",1);
873 Long64_t nevents = ctr->GetEntries();
874 if ( nevents == 0 ) {
875 printf(" The file is empty! \n ");
876 return(1);
877 };
878 ctr->SetBranchAddress("Header", &eh);
879 for (Int_t i = 0; i < nevents; i++){
880 ctr->GetEntry(i);
881 ph = eh->GetPscuHeader();
882 OBT = ph->GetOrbitalTime();
883 if ( !i ) OBT1 = OBT;
884 if ( i == nevents-1 ) LOBT = OBT;
885 if ( OOBT > OBT ) {
886 evjump[obtjump] = i;
887 obtjump++;
888 };
889 OOBT = OBT;
890 };
891 printf("\n Obtjump = %i - FIRST OBT %i - LAST OBT %i \n\n",obtjump,OBT1,LOBT);
892 //
893 TTree *tr;
894 pamela::CalibCalPedEvent *ce = 0;
895 pamela::PscuHeader *cph = 0;
896 pamela::EventHeader *ceh = 0;
897 Long64_t ncalibs;
898 //
899 for (Int_t s = 0; s < 4; s++){
900 for (Int_t d = 1; d<50; d++){
901 calib.ttime[s][d] = 0;
902 if ( d < 49 ) calib.time[s][d] = 0;
903 };
904 };
905 for (Int_t s = 0; s < 4; s++){
906 //
907 Int_t firstlook = 0;
908 repeatsearch:
909 //
910 Int_t inter = 0;
911 //
912 tr = (TTree*)File->Get("CalibCalPed");
913 tr->SetBranchAddress("Header", &ceh);
914 tr->SetBranchAddress("CalibCalPed", &ce);
915 if ( tr->GetEntries() == 0 ){
916 printf(" No calibration in this file! \n");
917 if ( !firstlook ){
918 wused = 2;
919 filename = calibfile;
920 firstlook = 1;
921 goto repeatsearch;
922 } else {
923 wused = 0;
924 goto jte;
925 };
926 };
927 wused = 1;
928 //
929 ncalibs = tr->GetEntries();
930 calib.obtjump = 0;
931 inter = 0;
932 for (Int_t c = 0; c < ncalibs; c++){
933 tr->GetEntry(c);
934 cph = ceh->GetPscuHeader();
935 calib.ttime[s][inter] = 0;
936 if ( ce->cstwerr[s] != 0 && ce->cperror[s] == 0 ) {
937 calib.ttime[s][inter] = cph->GetOrbitalTime();
938 inter++;
939 } else {
940 if ( ce->cstwerr[s] != 0 && ce->cperror[s] != 0 ) {
941 printf(" ERROR: entry %i stwerr %X perror %f \n",c,ce->cstwerr[s],ce->cperror[s]);
942 };
943 };
944 };
945 jte:
946 if ( inter == 0 ){
947 printf(" WARNING: no suitable calibration for section %i in this file!\n",s);
948 printf(" I WILL SEARCH IN PREVIOUS FILES\n");
949 if ( !firstlook ){
950 wused = 2;
951 filename = calibfile;
952 firstlook = 1;
953 goto repeatsearch;
954 };
955 };
956 if ( inter > 50 ){
957 printf(" WARNING: cannot handle more than 50 calibrations for file!\n");
958 printf(" I WILL SEARCH IN PREVIOUS FILES\n");
959 inter = 0;
960 if ( !firstlook ){
961 wused = 2;
962 filename = calibfile;
963 firstlook = 1;
964 goto repeatsearch;
965 };
966 };
967 if ( obtjump ){
968 calib.obtjump = 1;
969 printf(" ERROR: jump in OBT! CPU restarded? \n");
970 printf(" WARNING: I will use only the first calibration (if any)\n");
971 };
972 //
973 Int_t lastcal = 0;
974 Int_t of = 0;
975 if ( OBT1 < calib.ttime[s][0] || !inter ) {
976 //
977 // here we must look for a calibration in the previous file...
978 //
979 TString ffile = "";
980 stringcopy(ffile,filename,0,0);
981 Int_t fcodep = fetchpreviousfile(ffile,s);
982 if ( fcodep == 1 ) {
983 if ( inter ){
984 of = -3;
985 printf(" WARNING: Problems to find a good calibration for section %i \n",s);
986 } else {
987 of = -2;
988 printf(" ERROR: I was not able to find any good calibration for section %i \n",s);
989 };
990 } else {
991 stringcopy(ffile,filename,0,0);
992 TString pfile = whatnamewith(ffile,fcodep);
993 struct Calib tcal;
994 OLDCaloFindCalibs(pfile,tcal);
995 for ( Int_t d = 0; d<51; d++ ){
996 if ( tcal.ttime[s][d] == 0 ){
997 lastcal = d-1;
998 break;
999 };
1000 };
1001 if ( !inter ){
1002 calib.ttime[s][0] = tcal.ttime[s][lastcal];
1003 calib.time[s][0] = OBT1;
1004 calib.time[s][1] = LOBT;
1005 calib.fcode[s][0] = fcodep;
1006 if ( debug ) printf("boh! %i \n",fcodep);
1007 if ( debug ) getchar();
1008 of = -1;
1009 } else {
1010 if ( obtjump ){
1011 if ( debug ) printf("boh2! %i \n",fcodep);
1012 if ( debug ) getchar();
1013 inter = 1;
1014 calib.time[s][0] = OBT1;
1015 calib.time[s][1] = LOBT;
1016 calib.fcode[s][0] = 10;
1017 of = -5;
1018 for (Int_t d = 1; d<50; d++){
1019 calib.ttime[s][d] = 0;
1020 if ( d>1 && d<49 ) calib.time[s][d] = 0;
1021 };
1022 } else {
1023 if ( lastcal >= 0 && tcal.ttime[s][lastcal]<OBT1 && (OBT1-tcal.ttime[s][lastcal])<7200000 && (OBT1-tcal.ttime[s][lastcal])>0 ){
1024 for (Int_t i = inter+1; i>0; i--){
1025 calib.ttime[s][i] = calib.ttime[s][i-1];
1026 };
1027 calib.time[s][0] = OBT1;
1028 calib.ttime[s][0] = tcal.ttime[s][lastcal];
1029 calib.fcode[s][0] = fcodep;
1030 of = 2;
1031 } else {
1032 calib.time[s][0] = OBT1;
1033 calib.fcode[s][0] = 10;
1034 of = 1;
1035 };
1036 };
1037 };
1038 };
1039 } else {
1040 of = 0;
1041 };
1042 if ( debug ) printf("boh3! %i \n",of);
1043 Int_t of1 = 0;
1044 if ( inter && of != 0 && of != 2 && of != -2 && of != -5) {
1045 for (Int_t i = inter+1; i>0; i--){
1046 calib.ttime[s][i] = calib.ttime[s][i-1];
1047 };
1048 of1 = 0;
1049 };
1050 if ( of == -2 ) of = -1;
1051 if ( of == -3 ) {
1052 of = 1;
1053 calib.fcode[s][0] = 10;
1054 };
1055 //
1056 if ( of != -1 && of != -5 ){
1057 if ( of == 2 ) of = 1;
1058 for (Int_t d = 0; d<inter+1; d++){
1059 calib.time[s][d+of] = calib.ttime[s][d+of];
1060 calib.fcode[s][d+of] = 10;
1061 if ( d == inter ){
1062 if ( debug ) printf("boh4! %i \n",of);
1063 calib.time[s][d+of1+of] = LOBT;
1064 calib.fcode[s][d+of] = 0;
1065 };
1066 };
1067 };
1068 // if ( calheadFile ) calheadFile->Close();
1069 // if ( calcalibFile ) calcalibFile->Close();
1070 //if ( File ) File->Close();
1071 };
1072 return(0);
1073 }
1074
1075 void Calo1stcalib(TString filename, Calib & calib, Int_t b[4]){
1076 Float_t estrip[2][22][96];
1077 //
1078 // this is the value of the mip for each strip. To be changed when we will have the real values
1079 //
1080 for (Int_t s=0; s<4;s++){
1081 for (Int_t d = 0; d<50; d++){
1082 calib.ttime[s][d] = 0 ;
1083 if ( d < 49 ) calib.time[s][d] = 0 ;
1084 };
1085 };
1086 for (Int_t m = 0; m < 2 ; m++ ){
1087 for (Int_t k = 0; k < 22; k++ ){
1088 for (Int_t l = 0; l < 96; l++ ){
1089 calib.calped[m][k][l] = 0. ;
1090 estrip[m][k][l] = 0.;
1091 };
1092 };
1093 }
1094 //
1095 // first of all find the calibrations in the file
1096 //
1097 OLDCaloFindCalibs(filename, calib);
1098 //
1099 // print on the screen the results:
1100 //
1101 printf(" ---------------------------------------------------------- \n \n");
1102 Int_t calibex = 0;
1103 for (Int_t s=0; s<4;s++){
1104 Int_t stop = 0;
1105 for (Int_t d = 0; d<48; d++){
1106 if ( calib.ttime[s][d] != 0 ) calibex++;
1107 if ( calib.time[s][0] != 0 ){
1108 if ( d == 0 ) printf(" Section %i from time 0 to time %i use calibration at time %i \n",s,calib.time[s][d],calib.ttime[s][d]);
1109 if ( calib.time[s][d+1] != 0 ) {
1110 printf(" Section %i from time %i to time %i use calibration at time %i \n",s,calib.time[s][d],calib.time[s][d+1],calib.ttime[s][d+1]);
1111 } else {
1112 if ( !stop ){
1113 printf(" Section %i from time %i use calibration at time %i \n",s,calib.time[s][d],calib.ttime[s][d+1]);
1114 stop = 1;
1115 };
1116 };
1117 } else {
1118 if ( calib.ttime[s][d] != 0 ) printf(" Section %i from time 0 use calibration at time %i \n",s,calib.ttime[s][d]);
1119 };
1120 };
1121 printf("\n");
1122 };
1123 printf(" ---------------------------------------------------------- \n");
1124 if ( calibex < 1 ) {
1125 printf("No full calibration data in this file, sorry!\n");
1126 } else {
1127 //
1128 // calibrate before starting
1129 //
1130 for (Int_t s = 0; s < 4; s++){
1131 b[s]=0;
1132 CaloPede(filename,s,calib.ttime[s][b[s]],calib);
1133 };
1134 };
1135 }
1136
1137
1138 void ColorMIP(Float_t mip, int& colo){
1139 // printf("mip = %f \n",mip);
1140 if ( colo > 0 ){
1141 colo = 10;
1142 if ( mip > 0.7 ) colo = 38;
1143 if ( mip > 2. ) colo = 4;
1144 if ( mip > 10. ) colo = 3;
1145 if ( mip > 100. ) colo = 2;
1146 if ( mip > 500. ) colo = 6;
1147 } else {
1148 colo = 10;
1149 if ( mip > 0.7 ) colo = 17;
1150 if ( mip > 2. ) colo = 15;
1151 if ( mip > 10. ) colo = 14;
1152 if ( mip > 100. ) colo = 13;
1153 if ( mip > 500. ) colo = 12;
1154 };
1155 }
1156
1157 void ColorTOFMIP(Float_t mip, int& colo){
1158 // printf("mip = %f \n",mip);
1159 if ( colo > 0 ){
1160 colo = 10;
1161 if ( mip > 0. ) colo = 38;
1162 if ( mip > 2. ) colo = 4;
1163 if ( mip > 10. ) colo = 3;
1164 if ( mip > 100. ) colo = 2;
1165 if ( mip > 500. ) colo = 6;
1166 } else {
1167 colo = 10;
1168 if ( mip > 0. ) colo = 17;
1169 if ( mip > 2. ) colo = 15;
1170 if ( mip > 10. ) colo = 14;
1171 if ( mip > 100. ) colo = 13;
1172 if ( mip > 500. ) colo = 12;
1173 };
1174 }
1175
1176 void *processevents(void *ptr){
1177 char *tellme = (char *)ptr;
1178 cin.getline(tellme,256);
1179 delete TThread::Self();
1180 return(NULL);
1181 }
1182
1183 Int_t WhatToDo(int& i, int& doflag, int& jumpto, Long64_t nevents, Float_t lastevno, Float_t firstevno, const char *file, TString outDir, TString figty, TCanvas &figure) {
1184 char *bw;
1185 if ( jumpto == -1 ){
1186 bw = "bw";
1187 } else {
1188 bw = "";
1189 };
1190 jumpto = 0;
1191 char tellme[256];
1192 stringstream input;
1193 char tellme2[256];
1194 stringstream input2;
1195 char tellme3[256];
1196 stringstream input3;
1197 stringstream out;
1198 out.str("x");
1199 input2.str("zzzzzzzzzzzzzz");
1200 input3.str("z");
1201 input << out.str().c_str();
1202 stringstream stc;
1203 stringstream stc2;
1204 while ( !strcmp(input.str().c_str(),out.str().c_str()) ) {
1205 printf("\nPress <enter> to continue, b<enter> to go backward, j<enter> to jump to a\nparticular event number, p<enter> to save the figure, q<enter> to quit: \n");
1206 //
1207 input.str("");
1208 if ( !gROOT->GetListOfClasses()->FindObject("TRint") ){
1209 TThread *thread = new TThread((TThread::VoidFunc_t)processevents,(void *)tellme);
1210 thread->Run();
1211 while ( thread->GetState() == TThread::kRunningState ){
1212 gSystem->ProcessEvents();
1213 };
1214 thread->Kill();
1215 } else {
1216 cin.getline(tellme,256);
1217 };
1218 input << tellme;
1219 //
1220 out << "y";
1221 //
1222 stc.str("b");
1223 //
1224 if ( !strcmp(input.str().c_str(),stc.str().c_str()) ) {
1225 if ( i > 0 ) {
1226 printf("WARNING: going backward!\n\n");
1227 doflag = 2;
1228 } else {
1229 printf("This is the first event, you can't go backward! \n");
1230 out.str("");
1231 out << input.str().c_str();
1232 };
1233 };
1234 //
1235 stc.str("j");
1236 if ( !strcmp(input.str().c_str(),stc.str().c_str()) ) {
1237 printf("\n Do you want to jump to a progressive number [p] or to an event number [e]? ");
1238 cin.getline(tellme3,256);
1239 input3.str("");
1240 input3 << tellme3;
1241 stc.str("p");
1242 if ( !strcmp(input3.str().c_str(),stc.str().c_str()) ) {
1243 printf("\n Enter the progressive number you want to jump to: ");
1244 cin.getline(tellme2,256);
1245 input2.str("");
1246 input2 << tellme2;
1247 Int_t j;
1248 j = atoi(input2.str().c_str());
1249 if ( j < 1 || j > nevents+1 ) {
1250 printf("\n You can choose between 1 and %i \n",(int)nevents+1);
1251 out.str("");
1252 out << input.str().c_str();
1253 } else {
1254 printf("\n Jumping to progressive number %i\n\n",j);
1255 i = j-2;
1256 };
1257 };
1258 stc.str("e");
1259 if ( !strcmp(input3.str().c_str(),stc.str().c_str()) ) {
1260 printf("\n Enter the event number you want to jump to: ");
1261 cin.getline(tellme2,256);
1262 input2.str("");
1263 input2 << tellme2;
1264 Int_t j;
1265 j = atoi(input2.str().c_str());
1266 if ( j < firstevno || j > lastevno ) {
1267 printf("\n You can choose between %i and %i \n",(int)firstevno,(int)lastevno);
1268 out.str("");
1269 out << input.str().c_str();
1270 } else {
1271 printf("\n Jumping to event number %i\n\n",j);
1272 jumpto = j;
1273 i = -1;
1274 };
1275 };
1276 stc.str("e");
1277 stc2.str("p");
1278 if ( strcmp(input3.str().c_str(),stc.str().c_str()) && strcmp(input3.str().c_str(),stc2.str().c_str()) ) {
1279 printf(" You must type or \"p\" or \"e\"\n");
1280 out.str("");
1281 out << input.str().c_str();
1282 };
1283 };
1284 //
1285 stc.str("q");
1286 if ( !strcmp(input.str().c_str(),stc.str().c_str()) ) {
1287 printf("Exiting...\n");
1288 return(1);
1289 };
1290 //
1291 stc.str("p");
1292 if ( !strcmp(input.str().c_str(),stc.str().c_str()) ) {
1293 printf("Enter a file extension recognized by ROOT (ps, eps, gif,...):\n");
1294 cin.getline(tellme2,256);
1295 input2.str("");
1296 input2 << tellme2;
1297 out.str("");
1298 out << input.str().c_str();
1299 //
1300 TString filename = file;
1301 const string fil = (const char*)filename;
1302 Int_t posiz = fil.find("dw_");
1303 if ( posiz == -1 ) posiz = fil.find("DW_");
1304 Int_t posiz2 = posiz+13;
1305 TString file2;
1306 stringcopy(file2,filename,posiz,posiz2);
1307 const char *figrec = file2;
1308 const char *outdir = outDir;
1309 stringstream figsave;
1310 const char *ty = figty;
1311 figsave.str("");
1312 figsave << outdir << "/";
1313 figsave << ty << "_";
1314 figsave << (i+1) << "_";
1315 figsave << figrec;
1316 figsave << bw << ".";
1317 figsave << input2.str().c_str();
1318 figure.SaveAs(figsave.str().c_str());
1319 printf("\n");
1320 };
1321 };
1322 return(0);
1323 }
1324
1325 void PrintFigure(TString filename, TString outDir, TString figty, TString saveas, TCanvas& figure) {
1326 const string fil = (const char*)filename;
1327 Int_t posiz = fil.find("dw_");
1328 if ( posiz == -1 ) posiz = fil.find("DW_");
1329 Int_t posiz2 = posiz+13;
1330 TString file2;
1331 stringcopy(file2,filename,posiz,posiz2);
1332 //
1333 const char *figrec = file2;
1334 stringstream figsave;
1335 figsave.str("");
1336 figsave << outDir.Data() << "/";
1337 figsave << figrec << "_";
1338 figsave << figty.Data() << ".";
1339 figsave << saveas.Data();
1340 figure.SaveAs(figsave.str().c_str());
1341 return;
1342 }
1343
1344 Double_t langaufun(Double_t *x, Double_t *par) {
1345 // Numeric constants
1346 Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
1347 Double_t mpshift = -0.22278298; // Landau maximum location
1348 // Control constants
1349 Double_t np = 100.0; // number of convolution steps
1350 Double_t sc = 3.; // convolution extends to +-sc Gaussian sigmas
1351 // Variables
1352 Double_t xx;
1353 Double_t mpc;
1354 Double_t fland;
1355 Double_t sum = 0.0;
1356 Double_t xlow,xupp;
1357 Double_t step;
1358 Double_t i;
1359 // MP shift correction
1360 mpc = par[1] - mpshift * par[0];
1361 // Range of convolution integral
1362 xlow = x[0] - sc * par[3];
1363 xupp = x[0] + sc * par[3];
1364 step = (xupp-xlow) / np;
1365 // Convolution integral of Landau and Gaussian by sum
1366 for(i=1.0; i<=np/2; i++) {
1367 xx = xlow + (i-.5) * step;
1368 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
1369 sum += fland * TMath::Gaus(x[0],xx,par[3]);
1370 xx = xupp - (i-.5) * step;
1371 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
1372 sum += fland * TMath::Gaus(x[0],xx,par[3]);
1373 };
1374 return (par[2] * step * sum * invsq2pi / par[3]);
1375 }
1376
1377
1378
1379 TF1 *langaufit(TH1F *his, Double_t *fitrange, Double_t *startvalues, Double_t *parlimitslo, Double_t *parlimitshi, Double_t *fitparams, Double_t *fiterrors, Double_t *ChiSqr, Int_t *NDF, TString drw="QR0"){
1380 // Once again, here are the Landau * Gaussian parameters:
1381 // par[0]=Width (scale) parameter of Landau density
1382 // par[1]=Most Probable (MP, location) parameter of Landau density
1383 // par[2]=Total area (integral -inf to inf, normalization constant)
1384 // par[3]=Width (sigma) of convoluted Gaussian function
1385 //
1386 // Variables for langaufit call:
1387 // his histogram to fit
1388 // fitrange[2] lo and hi boundaries of fit range
1389 // startvalues[4] reasonable start values for the fit
1390 // parlimitslo[4] lower parameter limits
1391 // parlimitshi[4] upper parameter limits
1392 // fitparams[4] returns the final fit parameters
1393 // fiterrors[4] returns the final fit errors
1394 // ChiSqr returns the chi square
1395 // NDF returns ndf
1396 Int_t i;
1397 Char_t FunName[100];
1398 sprintf(FunName,"Fitfcn_%s",his->GetName());
1399 //
1400 TF1 *ffitold = (TF1*)gROOT->GetListOfFunctions()->FindObject(FunName);
1401 if (ffitold) delete ffitold;
1402 //
1403 TF1 *ffit = new TF1(FunName,langaufun,fitrange[0],fitrange[1],4);
1404 ffit->SetParameters(startvalues);
1405 ffit->SetParNames("Width","MP","Area","GSigma");
1406 //
1407 for (i=0; i<4; i++) {
1408 ffit->SetParLimits(i, parlimitslo[i], parlimitshi[i]);
1409 };
1410 his->Fit(FunName,drw); // fit within specified range, use ParLimits, do not plot
1411 //
1412 ffit->GetParameters(fitparams); // obtain fit parameters
1413 for (i=0; i<4; i++) {
1414 fiterrors[i] = ffit->GetParError(i); // obtain fit parameter errors
1415 }
1416 ChiSqr[0] = ffit->GetChisquare(); // obtain chi^2
1417 NDF[0] = ffit->GetNDF(); // obtain ndf
1418 //
1419 return (ffit); // return fit function
1420 }
1421
1422
1423 //Double_t langaupro(Double_t *params, Double_t &maxx) {
1424 Double_t langaupro(Double_t *params) {
1425 // Seaches for the location (x value) at the maximum of the
1426 // Landau-Gaussian convolute and its full width at half-maximum.
1427 //
1428 Double_t p,x;
1429 Double_t step;
1430 Double_t l,lold;
1431 Int_t i = 0;
1432 Int_t MAXCALLS = 5000;
1433 // Search for maximum
1434 p = params[1] - 0.1 * params[0];
1435 step = 0.05 * params[0];
1436 lold = -2.0;
1437 l = -1.0;
1438 // while ( (l != lold) && (i < MAXCALLS) ) {
1439 while ( fabs(l-lold)>10e-12 && (i < MAXCALLS) ) {
1440 i++;
1441 lold = l;
1442 x = p + step;
1443 l = langaufun(&x,params);
1444 if ( l < lold ) step = -step/10.;
1445 p += step;
1446 // printf(" i %i: l %.16g lold %.16g x %.16g p %.16g step %.16g \n",i,l,lold,x,p,step);
1447 // printf(" l %.16g \n",l);
1448 };
1449 // printf(" i %i: l %.16g lold %.16g x %.16g p %.16g step %.16g \n",i,l,lold,x,p,step);
1450 //printf(" l %.16g \n",l);
1451 if (i == MAXCALLS) return (-1.);
1452 //
1453 return(x);
1454 }
1455
1456
1457 void delay(int selftrig, int &del){
1458 int i;
1459 i = 0;
1460 del = 300;
1461 //
1462 while ((i < 16) && ((selftrig & (0x8000>>i)) == 0)) {
1463 del += 50;
1464 i++;
1465 };
1466 if (i > 15) del = 1100;
1467 }
1468
1469 Double_t fitraw(Double_t *x, Double_t *par){
1470 //
1471 Double_t fitval = par[3] + par[0]*TMath::Exp(-par[2]*(x[0]-par[1]));
1472 //
1473 return fitval;
1474 }
1475
1476
1477
1478 int checkmerging(Chmerge & chmerge){
1479 for (Int_t i = chmerge.fromflag; i<1000; i++){
1480 if ( chmerge.pscuOBT == 0 && chmerge.pscucount == 0 ) return(1);
1481 if ( chmerge.pscuOBT == chmerge.lastOBT[i] && chmerge.pscucount == chmerge.lastcount[i] ){
1482 chmerge.fromflag = i;
1483 printf("\n WARNING! %i event number %i at OBT %i already processed, skipped!\n",i,chmerge.pscucount,chmerge.pscuOBT);
1484 return(1);
1485 };
1486 };
1487 return(0);
1488 }
1489
1490 int fillmerging(Chmerge & chmerge){
1491 chmerge.lastOBT[chmerge.scanflag] = chmerge.pscuOBT;
1492 chmerge.lastcount[chmerge.scanflag] = chmerge.pscucount;
1493 chmerge.scanflag++;
1494 if ( chmerge.scanflag > 999 ) return(0);
1495 chmerge.lastOBT[chmerge.scanflag] = 0;
1496 chmerge.lastcount[chmerge.scanflag] = 0;
1497 return(0);
1498 }
1499
1500 //short int checkifpretyodaversion(TString filename){
1501 // string fil = (const char *)filename;
1502 // Int_t posiz = fil.find("dw_");
1503 // if ( posiz == -1 ) posiz = fil.find("DW_");
1504 // if ( posiz == -1 ) return (-1);
1505 // TString year;
1506 // stringcopy(year,filename,posiz+3,posiz+5);
1507 // const char *ye = year;
1508 // Int_t y = atoi(ye);
1509 // TString month;
1510 // stringcopy(month,filename,posiz+5,posiz+7);
1511 // const char *mo = month;
1512 // Int_t m = atoi(mo);
1513 // TString day;
1514 // stringcopy(day,filename,posiz+7,posiz+9);
1515 // const char *da = day;
1516 // Int_t d = atoi(da);
1517 // TString fno;
1518 // stringcopy(fno,filename,posiz+10,posiz+13);
1519 // const char *fn = fno;
1520 // Int_t f = atoi(fn);
1521 // if ( y >= 5 && m >= 5 && d < 15 ) return(1);
1522 // if ( y == 5 && m == 5 && d == 15 && f < 7 ) return(1);
1523 // return(0);
1524 //}
1525
1526

  ViewVC Help
Powered by ViewVC 1.1.23