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

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 1.1.1.1 - (show annotations) (download) (vendor branch)
Tue Dec 6 10:59:28 2005 UTC (19 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 //
2 //- Emiliano Mocchiutti
3 //
4 // CaloADC2MIP.c version 4.01 (2005-11-18)
5 //
6 // The only input needed is
7 //
8 // Changelog:
9 //
10 // 4.00 - 4.01 (2005-11-18): use pathtocalibration library.
11 //
12 // 3.06 - 4.00 (2005-11-17): compiled!
13 //
14 // 3.05 - 3.06 (2005-08-15): 64 bit arch bugs fixed.
15 //
16 // 3.04 - 3.05 (2005-07-21): don't load yodaUtility.c anymore, use clone routines instead.
17 //
18 // 3.03 - 3.04 (2005-07-12): small changes for compilation.
19 //
20 // 3.02 - 3.03 (2005-06-28): changes for the new software architecture.
21 //
22 // 3.01 - 3.02 (2005-06-07): Check the date to see if we can run into trouble with the yoda version.
23 //
24 // 3.00 - 3.01 (2005-05-17): Changes accordingly to changes in calofindcalibs
25 //
26 // 2.02 - 3.00 (2005-04-18): Several changes, added many subroutines...
27 //
28 // 2.01 - 2.02 (2005-04-11): Do not calibrate more than once if there is an OBT jump.
29 //
30 // 2.00 - 2.01 (2005-04-07): Changes accordly to changes in the calofindcalibs subroutine. Changed tree name of calibration data.
31 //
32 // 2.00 (2005-03-14): Several changes, probably this is the first stable version.
33 //
34 // 1.00 (2005-03-04): Created and partially working.
35 //
36 #include <fstream>
37 #include <sstream>
38 #include <string>
39 #include <TTree.h>
40 #include <TClassEdit.h>
41 #include <TObject.h>
42 #include <TList.h>
43 #include <TSystem.h>
44 #include <TSystemDirectory.h>
45 #include <TString.h>
46 #include <TFile.h>
47 #include <TClass.h>
48 #include <TCanvas.h>
49 #include <TH1.h>
50 #include <TH1F.h>
51 #include <TH2D.h>
52 #include <TLatex.h>
53 #include <TPad.h>
54 #include <TPaveLabel.h>
55 #include <TChain.h>
56 #include <TGraph.h>
57 #include <TMultiGraph.h>
58 #include <TLine.h>
59 #if !defined (__CINT__)
60 extern const char *pathtocalibration();
61 #include <event/PamelaRun.h>
62 #include <event/physics/calorimeter/CalorimeterEvent.h>
63 #include <event/physics/trigger/TriggerEvent.h>
64 #include <event/CalibCalPedEvent.h>
65 #include <caloclasses.h>
66 #endif
67 #include <CaloFunctions.h>
68
69 void CaloADC2MIP(TString Filename, TString calcalibfile = "", TString flist=""){
70 gROOT->GetListOfCanvases()->Delete("a");
71 gDirectory->GetList()->Delete();
72 #if defined (__CINT__)
73 emicheckLib();
74 stringstream llib;
75 const char *pamlib = gSystem->ExpandPathName("$PAM_LIB");
76 llib << pamlib << "/caloclasses_h.so";
77 gSystem->Load(llib.str().c_str());
78 stringstream lolib;
79 lolib.str("");
80 lolib << pamlib << "/UTpathtoc_c.so";
81 gSystem->Load(lolib.str().c_str());
82 #endif
83 TString calcalib;
84 struct Evento evento;
85 struct Calib calib;
86 Int_t debug = 0;
87 Int_t b[4];
88 Int_t etime,evno;
89 Float_t sdexy[2][22][96],sdexyc[2][22][96];
90 Float_t edelta[2][22][96];
91 Float_t estrip[2][22][96];
92 Float_t estripnull[2][22][96];
93 Float_t delta;
94 delta = 1.;
95 evento.emin = 0.7;
96 Int_t UPYES = 0;
97 Int_t se =5;
98 Int_t pre;
99 Int_t rdone=0;
100 Int_t done;
101 bool isCOMP = 0;
102 bool isFULL = 0;
103 bool isRAW = 0;
104 TH1F *pfmip[2][22][96];
105 TF1 *pfitsnr[2][22][96];
106 stringstream Fmip;
107 for (Int_t i=0; i<2;i++){
108 for (Int_t j=0; j<22;j++){
109 for (Int_t m=0; m<96;m++){
110 Fmip.str("");
111 Fmip << "Fmip " << i;
112 Fmip << " " << j;
113 Fmip << " " << m;
114 pfmip[i][j][m] = new TH1F(Fmip.str().c_str(),"",56,0.,55.);
115 pfitsnr[i][j][m] = new TF1;
116 };
117 };
118 };
119 //
120 // the file which contains the calorimeter ADC to MIP conversion values is file2 = path to the script dir / CaloADC2MIP.root
121 //
122 const char *pam_calib = pathtocalibration();
123 // const char *pam_calib=gSystem->Getenv("PAM_CALIB");
124 stringstream file2;
125 // char *file2;
126 //file2 = "/misc/wizard/wizard3/pamela/integr/script/caloScripts-2.00/devel/CaloADC2MIP.rootasdalsdjaksjdlkasjdlkajsdklajsdlkjakljal/";
127 file2.str("");
128 file2 << pam_calib << "/CaloADC2MIP.root";
129 // sprintf(file2,"%s/CaloADC2MIP.root",pam_calib);
130 //
131 // this the figures and fit (4224 x 2 objects!)
132 //
133 stringstream file3;
134 file3.str("");
135 file3 << pam_calib << "/CaloADC2MIPf.root";
136 // char *file3;
137 //file3 = "/misc/wizard/wizard3/pamela/integr/script/caloScripts-2.00/devel/CaloADC2MIP.rootasdalsdjaksjdlkasjdlkajsdklajsdlkjakljal/";
138 //sprintf(file3,"%s/CaloADC2MIPf.root",pam_calib);
139 //
140 stringstream file4;
141 file4.str("");
142 // char *file4;
143 //file4 = "/misc/wizard/wizard3/pamela/integr/script/caloScripts-2.00/devel/CaloADC2MIP.rootasdalsdjaksjdlkasjdlkajsdklajsdlkjakljal/";
144 // sprintf(file4,"%s/CaloADC2MIPf",pam_calib);
145 file4 << pam_calib << "/CaloADC2MIPf";
146 //
147 stringstream file5;
148 file5.str("");
149 // char *file5;
150 //file5 = "/misc/wizard/wizard3/pamela/integr/script/caloScripts-2.00/devel/CaloADC2MIP.rootasdalsdjaksjdlkasjdlkajsdklajsdlkjakljal/";
151 // sprintf(file5,"%s/CaloADC2MIPdata.root",pam_calib);
152 file5 << pam_calib << "/CaloADC2MIPdata.root";
153 //
154 Int_t nfile = 0;
155 Int_t fileup = 0;
156 if ( flist != "" ){
157 printf("\n You have given an input file list \n\n Files: \n");
158 ifstream in;
159 in.open(flist, ios::in);
160 char ctime[20];
161 while (1) {
162 in >> ctime;
163 if (!in.good()) break;
164 printf(" File number %i = %s \n",nfile+1,ctime);
165 nfile++;
166 // if ( nfile > 499 ) {
167 // printf(" No more than 500 files please! \n ");
168 // return;
169 //};
170 };
171 in.close();
172 if ( !nfile ) {
173 printf("\n\n The list is empty, exiting... \n\n");
174 return;
175 };
176 } else {
177 nfile = 1;
178 };
179 printf("\n");
180 //
181 // First check if the calibration file exists:
182 //
183 Int_t EFILE = 0;
184 printf(" Check the existence of CaloADC2MIP.root file: \n\n");
185 TFile tfile(file2.str().c_str());
186 // TFile tfile2(file3);
187 // if ( !tfile.IsZombie() && !tfile2.IsZombie()) {
188 if ( !tfile.IsZombie() ) {
189 printf(" The file already exists! Check if an update is needed \n\n");
190 EFILE = 1;
191 tfile.Close();
192 } else {
193 printf("\n OK, I will create it!\n\n");
194 EFILE = 0;
195 };
196 //
197 // the check for an update is done once even if we have a list.
198 //
199 TTree *tree;
200 CalorimeterCalibration *calo = new CalorimeterCalibration();
201 if ( EFILE == 1) {
202 //
203 // The file exists, open it (READ ONLY) and check if any update is needed:
204 //
205 TFile *hfile = 0;
206 hfile = new TFile(file2.str().c_str(),"READ","Calorimeter CALIBRATION data");
207 tree = (TTree*)hfile->Get("CaloADC");
208 // CalorimeterCalibration *
209 calo = new CalorimeterCalibration();
210 tree->SetBranchAddress("Event", &calo);
211 //
212 Long64_t nevents = tree->GetEntries();
213 printf("\n Number of updates: %i\n",(int)nevents);
214 tree->GetEntry(nevents-1);
215 //
216 // Read the status variable of the last event: if it is 0 then we need an update, if it is one print out informations and exit
217 //
218 if ( calo->status == 0 ) {
219 printf("\n An update is needed: \n");
220 Int_t upstrip = 0;
221 for (Int_t i = 0; i<2; i++){
222 for (Int_t j = 0; j<22; j++){
223 for (Int_t l = 0; l<96; l++){
224 if ( calo->mask[i][j][l] == 0. ) upstrip++;
225 };
226 };
227 };
228 printf(" %i strip out of 4224 have no statistic enough.\n\n",upstrip);
229 } else {
230 printf(" No update is needed! \n The data files used to obtain the ADC to MIP conversion are: \n");
231 for (Int_t i = 0; i<nevents; i++){
232 tree->GetEntry(i);
233 const char *fout = calo->fname;
234 printf(" - %s \n",fout);
235 };
236 printf("\n\n");
237 return;
238 };
239 hfile->Close();
240 };
241 //
242 const char *filename2 = Filename;
243 //
244 TString lfile;
245 TString sfile;
246 //
247 Int_t e = 0;
248 ifstream in;
249 in.open(flist, ios::in);
250 char ctime[20];
251 TString file2f;
252 TString pfile;
253 const char *cali;
254 const char *ffile;
255 Int_t calibex = 0;
256 //
257 // TString *file[17];
258 TString filename;
259 const char *nome;
260 TFile *hfile = 0;
261 calo = new CalorimeterCalibration();
262 Long64_t nevents;
263 Int_t used = 0;
264 const char *compare;
265 // TTree *tree = 0;
266 Int_t hf3 = 0;
267 //Int_t done = 0;
268 Int_t ufile = 0;
269 TFile *headerFile = 0;
270 TFile *caloFile = 0;
271 TFile *triggerFile = 0;
272 file2f = "";
273 pfile = "";
274 stringstream fmippa;
275 const char *porca;
276 Int_t S4;
277 Int_t upnn;
278 Int_t notgood = 0;
279 Int_t notgood1 = 0;
280 // char *files = 0;
281 stringstream files;
282 TFile *hfileb = 0;
283 //
284 Int_t pl1, pl2, strmin, strmax;
285 Int_t spre=0;
286 Int_t surstrip = 0;
287 Float_t estrippa = 0.;
288 Int_t surmam = 0;
289 Int_t surman = 0;
290 Int_t surmatrix[3][3];
291 //
292 TString fififile;
293 const char *file;
294 //
295 stringstream nofi;
296 nofi.str("");
297 TFile *hfile3 = 0;
298 //
299 while ( e < nfile ){
300 if ( flist != "" ){
301 in >> ctime;
302 nofi.str("");
303 nofi << filename2 << "/";
304 nofi << ctime;
305 lfile = TString(nofi.str().c_str());
306 // file = "dw_000000_000.dat";
307 // file = getFilename(lfile);
308 fififile = getFilename(lfile);
309 file = fififile;
310 nofi.str("");
311 nofi << file;
312 sfile = TString(nofi.str().c_str());
313 } else {
314 nofi.str("");
315 nofi << filename2;
316 lfile = TString(nofi.str().c_str());
317 // file = "dw_000000_000.dat";
318 //file = getFilename(lfile);
319 fififile = getFilename(lfile);
320 file = fififile;
321 nofi.str("");
322 nofi << file;
323 sfile = TString(nofi.str().c_str());
324 };
325 //
326 filename = lfile;
327 nome = filename;
328 printf("\n Processing file number %i: %s \n\n",e,nome);
329 nome = sfile;
330 //
331 // this must be done for each file.
332 //
333 const char *calname;
334 calname = "";
335 Int_t wused = 0;
336 TTree *otr;
337 pamela::calorimeter::CalorimeterEvent *de = 0;
338 pamela::trigger::TriggerEvent *trigger = 0;
339 pamela::PscuHeader *ph = 0;
340 pamela::EventHeader *eh = 0;
341 //
342 if ( EFILE == 1 ){
343 hfile = new TFile(file2.str().c_str(),"UPDATE","Calorimeter CALIBRATION data");
344 tree = (TTree*)hfile->Get("CaloADC");
345 calo = new CalorimeterCalibration();
346 tree->SetBranchAddress("Event", &calo);
347 //
348 // Check if the input file has already been used
349 //
350 printf(" Check if file %s has already been used... \n",nome);
351 nevents = tree->GetEntries();
352 used = 0;
353 for (Int_t i = 0; i<nevents; i++){
354 tree->GetEntry(i);
355 compare= calo->fname;
356 if ( !strcmp(compare,nome) ) {
357 used = 1;
358 break;
359 };
360 };
361 if ( used ) {
362 printf(" The file %s has already been used!\n You cannot run twice this program with that file!\n\n",nome);
363 hfile->Close();
364 goto jumpfile;
365 } else {
366 printf(" Good, the file has not been used yet. \n\n");
367 nevents = tree->GetEntries() -1;
368 tree->GetEntry(nevents);
369 };
370 } else {
371 //
372 // If the file doesn't exist create it
373 //
374 calo = new CalorimeterCalibration();
375 hfile = new TFile(file2.str().c_str(),"NEW","Calorimeter CALIBRATION data");
376 if ( hfile->IsZombie() ){
377 printf(" ERROR: Cannot create file \n");
378 return;
379 };
380 tree = new TTree("CaloADC","Calorimeter calibration data");
381 tree->Branch("Event","CalorimeterCalibration",&calo);
382 };
383 //
384 //
385 //
386 hf3 = 0;
387 if ( e == 0 ){
388 TFile ttfile(file3.str().c_str());
389 if ( !ttfile.IsZombie() ) {
390 printf(" Update figures file! \n\n");
391 ttfile.Close();
392 ufile = 1;
393 } else {
394 ufile = 0;
395 };
396 };
397 if ( (EFILE == 1 || ufile == 1) && e == 0 ){
398 printf(" Reading the old MIP figures to be updated...\n ");
399 hfile3 = new TFile(file3.str().c_str());
400 hf3 = 1;
401 //
402 stringstream fmino;
403 fmino.str("");
404 for (Int_t l = 0; l < 2; l++){
405 for (Int_t m = 0; m < 22; m++){
406 for (Int_t n =0 ; n < 96; n++){
407 fmino.str("");
408 fmino << "fmip " << l;
409 fmino << " " << m;
410 fmino << " " << n;
411 // char *name = Form("fmip %i %i %i",l,m,n);
412 TH1F *temp = (TH1F*)hfile3->Get(fmino.str().c_str());
413 pfmip[l][m][n] = (TH1F*)temp->Clone();
414 };
415 };
416 };
417 printf(" ...done!\n");
418 };
419 EFILE = 1;
420 //
421 // Check if the input calibration file contains any calibration data
422 //
423 done = 0;
424 printf(" Check for calorimeter calibration: \n");
425 calcalib = filename;
426 findcalibs: printf(" Calibration file: %s \n",calname);
427 calname = calcalib;
428 for (Int_t s=0; s<4;s++){
429 for (Int_t d = 0; d<50; d++){
430 calib.ttime[s][d] = 0 ;
431 if ( d < 49 ) calib.time[s][d] = 0 ;
432 };
433 };
434 for (Int_t m = 0; m < 2 ; m++ ){
435 for (Int_t k = 0; k < 22; k++ ){
436 for (Int_t l = 0; l < 96; l++ ){
437 calib.calped[m][k][l] = 0. ;
438 evento.dexy[m][k][l] = 0. ;
439 evento.dexyc[m][k][l] = 0. ;
440 edelta[m][k][l] = 0.;
441 estrip[m][k][l] = 0.;
442 estripnull[m][k][l] = 0.;
443 };
444 };
445 }
446 wused = 0;
447 CaloFindCalibs(calcalib, calcalib, wused, calib);
448 //if ( wused == 1 ) calcalibfile = filename;
449 //CaloFindCalibs(calcalib, calib);
450 //
451 // print on the screen the results:
452 //
453 ffile = filename;
454 printf(" ------ %s ------- \n \n",ffile);
455 calibex = 0;
456 for (Int_t s=0; s<4;s++){
457 printf(" ** SECTION %i **\n",s);
458 for (Int_t d = 0; d<51; d++){
459 if ( calib.ttime[s][d] != 0 ) {
460 calibex++;
461 if ( calib.fcode[s][d] != 10 ){
462 file2f = "";
463 stringcopy(file2f,calcalib);
464 pfile = (TString)whatnamewith(file2f,calib.fcode[s][d]);
465 } else {
466 pfile = (TString)calcalib;
467 };
468 ffile = pfile;
469 printf(" - from time %i to time %i use calibration at\n time %i, file: %s \n",calib.time[s][d],calib.time[s][d+1],calib.ttime[s][d],ffile);
470 };
471 };
472 printf("\n");
473 };
474 printf(" ----------------------------------------------------------------------- \n \n");
475 //
476 if ( calibex < 1 ) {
477 printf("No calibration data in this file!\n");
478 calcalib = calcalibfile;
479 if ( !done && calcalibfile != "" ) {
480 cali = calcalibfile;
481 printf(" Search in the calibration file %s \n",cali);
482 done = 1;
483 goto findcalibs;
484 } else {
485 printf(" Cannot find any calibration for this file! \n");
486 goto jumpfile;
487 };
488 };
489 if ( calibex < 4 ) {
490 printf(" WARNING! No full calibration data in this file! \n\n");
491 };
492 //
493 // upgrade data...
494 //
495 headerFile = emigetFile(filename, "Physics", "Header");
496 if ( !headerFile ){
497 printf("\n WARNING: No header file!\n");
498 goto jumpfile;
499 };
500 caloFile = emigetFile(filename, "Physics", "Calorimeter");
501 if ( !caloFile ){
502 printf("\n WARNING: No calorimeter file! \n\n");
503 headerFile->Close();
504 goto jumpfile;
505 };
506 triggerFile = emigetFile(filename, "Trigger");
507 if ( !triggerFile ){
508 printf("ERROR: no Trigger file! \n");
509 headerFile->Close();
510 caloFile->Close();
511 goto jumpfile;
512 }
513 //Takes the tree of the header file
514 otr = (TTree*)headerFile->Get("Pscu");
515 otr->AddFriend("Calorimeter",caloFile);
516 otr->AddFriend("Trigger",triggerFile);
517 //
518 otr->SetBranchAddress("Header", &eh);
519 otr->SetBranchAddress("Event", &de);
520 otr->SetBranchAddress("Trigger.Event", &trigger);
521 //
522 nevents = otr->GetEntries();
523 //
524 // calibrate before starting
525 //
526 for (Int_t s = 0; s < 4; s++){
527 b[s]=0;
528 if ( calib.fcode[s][b[s]] != 10 ){
529 stringcopy(file2f,calcalib);
530 pfile = (TString)whatnamewith(file2f,calib.fcode[s][b[s]]);
531 } else {
532 pfile = (TString)calcalib;
533 };
534 porca = pfile;
535 if ( debug ) printf(" porca miseria %s \n",porca);
536 CaloPede(pfile,s,calib.ttime[s][b[s]],calib);
537 };
538 //
539 // run over all the events
540 //
541 printf("\n Processed events: \n\n");
542 //
543 for (Int_t i = 0; i < nevents; i++){
544 evno = i;
545 if ( i%1000 == 0 && i > 0 ) printf(" %iK \n",i/1000);
546 printf(" %i \n",i);
547 //
548 // read from the header of the event the absolute time at which it was recorded
549 //
550 otr->GetEntry(i);
551 //
552 S4 = trigger->patterntrig[1] & (1<<0);
553 //
554 if ( !S4 || i == 0 ) {
555 ph = eh->GetPscuHeader();
556 etime = ph->GetOrbitalTime();
557 //
558 // for each event check that the calibration we are using are still within calibration limits, if not call the next calibration
559 //
560 if ( !calib.obtjump ) {
561 for (Int_t s = 0; s < 4; s++){
562 if ( calib.ttime[s][b[s]] ){
563 while ( etime > calib.time[s][b[s]+1] ){
564 printf(" CALORIMETER: \n" );
565 printf(" - Section %i, event at time %i while old calibration time limit at %i. Use new calibration at time %i -\n",s,etime,calib.time[s][b[s]+1],calib.ttime[s][b[s]+1]);
566 printf(" END CALORIMETER. \n\n" );
567 b[s]++;
568 if ( calib.fcode[s][b[s]] != 10 ){
569 stringcopy(file2f,calcalib);
570 pfile = (TString)whatnamewith(file2f,calib.fcode[s][b[s]]);
571 } else {
572 pfile = (TString)calcalib;
573 };
574 porca = pfile;
575 if ( debug ) printf(" porca miseria %s \n",porca);
576 CaloPede(pfile,s,calib.ttime[s][b[s]],calib);
577 };
578 };
579 };
580 };
581 //
582 // do whatever you want with data
583 //
584 evento.iev = de->iev;
585 //
586 // run over views and planes
587 //
588 for (Int_t l = 0; l < 2; l++){
589 for (Int_t m = 0; m < 22; m++){
590 //
591 // determine the section number
592 //
593 if (l == 0 && m%2 == 0) se = 3;
594 if (l == 0 && m%2 != 0) se = 2;
595 if (l == 1 && m%2 == 0) se = 0;
596 if (l == 1 && m%2 != 0) se = 1;
597 //
598 // determine what kind of event we are going to analyze
599 //
600 isCOMP = 0;
601 isFULL = 0;
602 isRAW = 0;
603 if ( de->stwerr[se] & (1 << 16) ) isCOMP = 1;
604 if ( de->stwerr[se] & (1 << 17) ) isFULL = 1;
605 if ( de->stwerr[se] & (1 << 3) ) isRAW = 1;
606 //
607 // save the prevoius energy deposit and calibration in sbase, sdexy, sdexyc
608 //
609 pre = -1;
610 if ( isRAW ){
611 for (Int_t nn = 0; nn < 96; nn++){
612 if ( nn%16 == 0 ) pre++;
613 evento.base[l][m][pre] = calib.calbase[l][m][pre];
614 sdexy[l][m][nn] = evento.dexy[l][m][nn];
615 evento.dexy[l][m][nn] = de->dexy[l][m][nn] ;
616 sdexyc[l][m][nn] = evento.dexyc[l][m][nn];
617 evento.dexyc[l][m][nn] = de->dexyc[l][m][nn] ;
618 };
619 };
620 //
621 // run over strips
622 //
623 pre = -1;
624 for (Int_t n =0 ; n < 96; n++){
625 if ( n%16 == 0 ) {
626 pre++;
627 rdone = 0;
628 done = 0;
629 };
630 if ( calo->mask[l][m][n] == 0. && ( isRAW || ( de->dexyc[l][m][n] > 0. && de->dexyc[l][m][n] < 32000.)) ){
631 //
632 // baseline check and calculation
633 //
634 if ( isRAW ) {
635 if ( !rdone ){
636 CaloFindBaseRaw(calib,evento,l,m,pre);
637 rdone = 1;
638 };
639 } else {
640 if ( !done ){
641 evento.base[l][m][pre] = de->base[l][m][pre] ;
642 evento.dexyc[l][m][n] = de->dexyc[l][m][n] ;
643 };
644 };
645 //
646 // no suitable new baseline, use old ones and compress data!
647 //
648 if ( !done && (evento.base[l][m][pre] > 30000. || evento.base[l][m][pre] == 0.) ){
649 evento.base[l][m][pre] = calib.sbase[l][m][pre];
650 upnn = n+16;
651 if ( upnn > 96 ) upnn = 96;
652 for ( Int_t nn = n; nn<upnn; nn++ ){
653 evento.dexyc[l][m][nn] = de->dexyc[l][m][nn] ;
654 };
655 CaloCompressData(calib,evento,l,m,pre);
656 done = 1;
657 };
658 //
659 // CALIBRATION ALGORITHM
660 //
661 pl1 = -1;
662 pl2 = -1;
663 pl1 = m - 1;
664 pl2 = m + 1;
665 if ( pl1 < 0 ) {
666 pl1 = pl2;
667 pl2 = m + 2;
668 };
669 if ( pl2 > 21 ) {
670 pl2 = pl1;
671 pl1 = m - 2;
672 };
673 strmin = n - 1;
674 if ( strmin < 0 ) strmin = 0;
675 strmax = n + 1;
676 if ( strmax > 95 ) strmax = 95;
677 surstrip = 0;
678 surmam = 0;
679 surman = 0;
680 for (Int_t im = 0; im<3 ;im++){
681 for (Int_t jm = 0; jm<3 ;jm++){
682 surmatrix[im][jm] = 0;
683 };
684 };
685 for (Int_t spl = pl1; spl<pl2+1; spl++){
686 if ( spl != m ){
687 surmam = spl - pl1;
688 for (Int_t sst = strmin; sst <= strmax ; sst++){
689 surman = sst - strmin;
690 if ( sst < 16 ) spre = 0;
691 if ( 15 < sst && sst < 32 ) spre = 1;
692 if ( 31 < sst && sst < 48 ) spre = 2;
693 if ( 47 < sst && sst < 64 ) spre = 3;
694 if ( 63 < sst && sst < 80 ) spre = 4;
695 if ( 79 < sst ) spre = 5;
696 if ( isRAW ){
697 estrippa = ( evento.dexy[l][spl][sst] - calib.calped[l][spl][sst] - calib.sbase[l][spl][spre])/26. ;
698 } else {
699 estrippa = ( evento.dexyc[l][spl][sst] - calib.calped[l][spl][sst] - calib.sbase[l][spl][spre])/26. ;
700 };
701 if ( estrippa > 0.7 && estrippa < 10. && calib.calgood[l][spl][sst] == 0 ){
702 surstrip++;
703 surmatrix[surmam][surman]++;
704 //printf("surstrip = %i estrip = %f \n",surstrip,estrip);
705 };
706 };
707 };
708 };
709 if ( (surmatrix[0][0] && surmatrix[2][0]) ) surstrip = 0;
710 if ( (surmatrix[0][2] && surmatrix[2][2]) ) surstrip = 0;
711 //
712 //
713 estrip[l][m][n] = ( evento.dexyc[l][m][n] - calib.calped[l][m][n] - evento.base[l][m][pre]) ;
714 if ( estrip[l][m][n] >= 0. && estrip[l][m][n] <= 56. && surstrip >=2 ){
715 pfmip[l][m][n]->Fill(estrip[l][m][n]);
716 UPYES = 1;
717 };
718 calib.sbase[l][m][pre] = evento.base[l][m][pre];
719 };
720 };
721 };
722 };
723 };
724 };
725 //
726 // save filename
727 //
728 calo->fname = sfile;
729 tree->Fill();
730 hfile->Write();
731 fileup++;
732 hfile->Close();
733 //
734 // create a backup copy!
735 //
736
737 if ( nfile > 1 && ( e%10 == 0 || e == nfile-1 ) && e > 0 ){
738 files.str("");
739 files << file4.str().c_str() << (e+1);
740 files << ".bak";
741 // files = Form("%s%i.bak",file4,e+1);
742 printf(" Save the backup file with figures:\n %s \n ",files.str().c_str());
743 hfileb = new TFile(files.str().c_str(),"RECREATE","Calorimeter CALIBRATION figures");
744 if ( hfileb->IsZombie() ){
745 printf(" ERROR: Cannot create file!!! \n");
746 goto oltre;
747 };
748 TH1F *bfmip[2][22][96];
749 stringstream bmippa;
750 bmippa.str("");
751 for (Int_t l = 0; l < 2; l++){
752 for (Int_t m = 0; m < 22; m++){
753 for (Int_t n =0 ; n < 96; n++){
754 bmippa.str("");
755 bmippa << "bmip " << l;
756 bmippa << " " << m;
757 bmippa << " " << n;
758 bfmip[l][m][n] = new TH1F(bmippa.str().c_str(),"",56,0.,55.);
759 bfmip[l][m][n] = (TH1F*)pfmip[l][m][n]->Clone(bmippa.str().c_str());
760 bfmip[l][m][n]->Write();
761 };
762 };
763 };
764 hfileb->Close();
765 };
766 //
767 //
768 //
769 headerFile->Close();
770 caloFile->Close();
771 triggerFile->Close();
772 goto oltre;
773 jumpfile: printf(" File not processed! \n ");
774 hfile->Close();
775 oltre:
776 // gObjectTable->Print();
777 e++;
778 };
779 //
780 if ( !UPYES ) goto end;
781 printf("\n OK, now I will fit the MIP distribution for each strip \n\n");
782 //
783 //TFile *hfile = 0;
784 hfile = new TFile(file2.str().c_str(),"UPDATE","Calorimeter CALIBRATION data");
785 // CalorimeterCalibration *calo = 0;
786 calo = new CalorimeterCalibration();
787 tree = (TTree*)hfile->Get("CaloADC");
788 tree->SetBranchAddress("Event", &calo);
789 nevents = tree->GetEntries();
790 if ( (nevents-1-fileup) > 0 ){
791 tree->GetEntry(nevents-1-fileup);
792 };
793 //
794 // Fit all 4224 distributions and save parameters, fit, figures (all basically).
795 //
796 notgood = 0;
797 notgood1 = 0;
798 //
799 // TCanvas *c1;
800 //c1 = new TCanvas("canvas","canvas",800,800);
801 //c1->Draw();
802 //
803 for (Int_t l = 0; l < 2; l++){
804 for (Int_t m = 0; m < 22; m++){
805 for (Int_t n =0 ; n < 96; n++){
806 Double_t SNRPeak = 0.;
807 if ( calo->mask[l][m][n] == 0. ){
808 //
809 // Setting fit range and start values
810 //
811 Double_t fr[2];
812 Double_t sv[4], pllo[4], plhi[4], fp[4], fpe[4];
813 if ( SNRPeak > 16. && SNRPeak < 35. ){
814 fr[0] = (Float_t)SNRPeak - 7.;
815 sv[0]= fp[0];
816 sv[1]= fp[1];
817 sv[2]= fp[2];
818 sv[3]= fp[3];
819 } else {
820 fr[0] = 19.;
821 sv[0] = 2.8;
822 sv[1] = 21.0;
823 sv[2] = 1000.0;
824 sv[3] = 0.1;
825 };
826 Int_t doitagain = 0;
827 fitting: printf("Fitting strip %i %i %i \n",l,m,n);
828 fr[1] = 60.;
829 pllo[0]=0.5; pllo[1]=5.0; pllo[2]=1.0; pllo[3]=0.2;
830 plhi[0]=8.0; plhi[1]=50.0; plhi[2]=1000000.0; plhi[3]=8.0;
831 Double_t chisqr;
832 Int_t ndf;
833 //
834 // fit the figure pfmip and put fit results in pfitsnr
835 //
836 pfitsnr[l][m][n] = langaufit(pfmip[l][m][n],fr,sv,pllo,plhi,fp,fpe,&chisqr,&ndf);
837 //
838 // search for the maximum of the distribution, it will be the conversion factor between ADC channels and the MIP
839 //
840 Double_t SNRPeak = langaupro(fp);
841 printf("\n Conversion factor: %f \n\n",SNRPeak);
842 //
843 //
844 //
845 if ( ( SNRPeak < 0. || (SNRPeak < 16. || SNRPeak > 35.) || chisqr > 100.) && doitagain < 3 ){
846 printf("\n The fit is not good, I will try again, step %i \n\n",doitagain);
847 doitagain++;
848 if ( chisqr > 100. ) {
849 fr[0] = fr[0] + 5.;
850 sv[0] = fp[0];
851 sv[1] = fp[1];
852 sv[2] = fp[2];
853 sv[3] = fp[3];
854 goto fitting;
855 };
856 if ( doitagain == 2 ) {
857 fr[0] = 19.;
858 sv[0] = 2.8;
859 sv[1] = 21.0;
860 sv[2] = 1000.0;
861 sv[3] = 0.1;
862 goto fitting;
863 };
864 if ( doitagain == 1 ) {
865 fr[0] = 22.;
866 sv[0] = 2.8;
867 sv[1] = 25.0;
868 sv[2] = 1000.0;
869 sv[3] = 0.1;
870 goto fitting;
871 };
872 if ( doitagain == 3 ) {
873 fr[0] = 12.;
874 sv[0] = 2.8;
875 sv[1] = 15.0;
876 sv[2] = 1000.0;
877 sv[3] = 0.1;
878 goto fitting;
879 };
880 };
881
882 //
883 // save data
884 //
885 calo->mip[l][m][n] = (float)SNRPeak;
886 calo->ermip[l][m][n] = (float)fpe[1];
887 for (Int_t a = 0; a < 4 ; a++){
888 calo->fp[a][l][m][n] = (float)fp[a];
889 calo->fpe[a][l][m][n] = (float)fpe[a];
890 };
891 calo->ndf[l][m][n] = (float)ndf;
892 calo->chi2[l][m][n] = (float)chisqr;
893 if ( calo->fp[1][l][m][n] > 16. && calo->fp[1][l][m][n] < 35. ) notgood1++;
894 //
895 if ( ndf!=0 && chisqr/ndf < 2. && fpe[1] < 0.15 && fp[1] > 15. && fp[1] < 40. && SNRPeak > 0. ) {
896 printf(" THIS IS A GOOD FIT, SAVED! \n\n");
897 calo->mask[l][m][n] = 1.;
898 } else {
899 notgood++;
900 calo->mask[l][m][n] = 0.;
901 };
902 printf("Fitting done, strip %i %i %i \n\n\n",l,m,n);
903 //
904 // draw the figure and fit on the screen
905 //
906 //c1->Clear();
907 //c1->cd();
908 //pfmip[l][m][n]->DrawCopy();
909 //pfitsnr[l][m][n]->DrawCopy("lsame");
910 //c1->Modified();
911 //c1->Update();
912 };
913 };
914 };
915 };
916 //
917 calo->status = 0;
918 if ( !notgood ) calo->status = 1;
919 //
920 // save fit data
921 //
922 tree->Fill();
923 hfile->Write();
924 hfile->Close();
925 //
926 // Save histograms and fit figures
927 //
928 TFile *hfile2;
929 hfile2 = new TFile(file3.str().c_str(),"RECREATE","Calorimeter CALIBRATION figures");
930 if ( hfile2->IsZombie() ){
931 printf(" ERROR: Cannot create file \n");
932 return;
933 };
934 TH1F *sfmip[2][22][96];
935 TF1 *sfitsnr[2][22][96];
936 fmippa.str("");
937 for (Int_t l = 0; l < 2; l++){
938 for (Int_t m = 0; m < 22; m++){
939 for (Int_t n =0 ; n < 96; n++){
940 fmippa.str("");
941 fmippa << "fmip " << l;
942 fmippa << " " << m;
943 fmippa << " " << n;
944 sfmip[l][m][n] = new TH1F(fmippa.str().c_str(),"",56,0.,55.);
945 sfmip[l][m][n] = (TH1F*)pfmip[l][m][n]->Clone(fmippa.str().c_str());
946 fmippa.str("");
947 fmippa << "fit " << l;
948 fmippa << " " << m;
949 fmippa << " " << n;
950 sfitsnr[l][m][n] = (TF1*)pfitsnr[l][m][n]->Clone(fmippa.str().c_str());
951 sfmip[l][m][n]->Write();
952 sfitsnr[l][m][n]->Write();
953 };
954 };
955 };
956 hfile2->Close();
957 if ( hf3 ) hfile3->Close();
958 //
959 if ( notgood ) {
960 printf("\n An update will be necessary:\n %i strip with too low statistic (%i of them gave a result anyway)\n\n",notgood,notgood1);
961 };
962 //
963 end: printf("\n");
964 //
965 // char *file4 = 0;
966 //file4 = "/misc/wizard/wizard3/pamela/integr/script/caloScripts-2.00/devel/CaloADC2MIP.rootasdalsdjaksjdlkasjdlkajsdklajsdlkjakljal/";
967 //sprintf(file4,"%s/CaloADC2MIPf*.bak",startingdir);
968 //char *remfile = 0;
969 //remfile = Form("rm -f %s",file4);
970 //gSystem->Exec(remfile);
971 }
972
973 void CaloRAWADC2MIPPLOT(TString Filename, TString calcalibfile = "", TString flist=""){
974 gROOT->GetListOfCanvases()->Delete("a");
975 gROOT->Reset("a");
976 gDirectory->GetList()->Delete();
977 #if defined (__CINT__)
978 emicheckLib();
979 stringstream llib;
980 const char *pamlib = gSystem->ExpandPathName("$PAM_LIB");
981 llib << pamlib << "/caloclasses_h.so";
982 gSystem->Load(llib.str().c_str());
983 stringstream lolib;
984 lolib.str("");
985 lolib << pamlib << "/UTpathtoc_c.so";
986 gSystem->Load(lolib.str().c_str());
987 #endif
988 TString calcalib;
989 struct Evento evento;
990 struct Calib calib;
991 Int_t debug = 0;
992 Int_t b[4];
993 Int_t etime,evno;
994 Float_t sdexy[2][22][96],sdexyc[2][22][96];
995 Float_t edelta[2][22][96];
996 Float_t estrip[2][22][96];
997 Float_t estripnull[2][22][96];
998 Float_t delta;
999 delta = 1.;
1000 evento.emin = 0.7;
1001 Int_t se =5;
1002 Int_t pre;
1003 Int_t rdone = 0;
1004 Int_t done;
1005 bool isCOMP = 0;
1006 bool isFULL = 0;
1007 bool isRAW = 0;
1008 TH1F *pfmip[2][22][96];
1009 TF1 *pfitsnr[2][22][96];
1010 stringstream Fmip;
1011 Fmip.str("");
1012 for (Int_t i=0; i<2;i++){
1013 for (Int_t j=0; j<22;j++){
1014 for (Int_t m=0; m<96;m++){
1015 Fmip.str("");
1016 Fmip << "Fmip " << i;
1017 Fmip << " " << j;
1018 Fmip << " " << m;
1019 pfmip[i][j][m] = new TH1F(Fmip.str().c_str(),"",70,-15.,55.);
1020 pfitsnr[i][j][m] = new TF1;
1021 };
1022 };
1023 };
1024 //
1025 // the file which contains the calorimeter ADC to MIP conversion values is file2 = path to the script dir / CaloADC2MIP.root
1026 //
1027 const char *pam_calib = pathtocalibration();
1028 // const char *pam_calib=gSystem->Getenv("PAM_CALIB");
1029 stringstream file2;
1030 file2.str("");
1031 file2 << pam_calib << "/CaloADC2MIP.raw";
1032 // sprintf(file2,"%s/CaloADC2MIP.raw",pam_calib);
1033 //
1034 // this the figures and fit (4224 x 2 objects!)
1035 //
1036 stringstream file3;
1037 file3 << pam_calib << "/CaloADC2MIPfr.raw";
1038 // char *file3;
1039 //file3 = "/misc/wizard/wizard3/pamela/integr/script/caloScripts-2.00/devel/CaloADC2MIP.rootasdalsdjaksjdlkasjdlkajsdklajsdlkjakljal/";
1040 //sprintf(file3,"%s/CaloADC2MIPfr.raw",pam_calib);
1041 //
1042 stringstream file4;
1043 file4.str("");
1044 // char *file4;
1045 //file4 = "/misc/wizard/wizard3/pamela/integr/script/caloScripts-2.00/devel/CaloADC2MIP.rootasdalsdjaksjdlkasjdlkajsdklajsdlkjakljal/";
1046 // sprintf(file4,"%s/CaloADC2MIPfr",pam_calib);
1047 file4 << pam_calib << "/CaloADC2MIPfr";
1048 //
1049 stringstream file5;
1050 file5.str("");
1051 file5 << pam_calib << "/CaloADC2MIPdata.raw";
1052 //
1053 Int_t nfile = 0;
1054 Int_t fileup = 0;
1055 if ( flist != "" ){
1056 printf("\n You have given an input file list \n\n Files: \n");
1057 ifstream in;
1058 in.open(flist, ios::in);
1059 char ctime[20];
1060 while (1) {
1061 in >> ctime;
1062 if (!in.good()) break;
1063 printf(" File number %i = %s \n",nfile+1,ctime);
1064 nfile++;
1065 // if ( nfile > 499 ) {
1066 // printf(" No more than 500 files please! \n ");
1067 // return;
1068 //};
1069 };
1070 in.close();
1071 if ( !nfile ) {
1072 printf("\n\n The list is empty, exiting... \n\n");
1073 return;
1074 };
1075 } else {
1076 nfile = 1;
1077 };
1078 printf("\n");
1079 //
1080 // First check if the calibration file exists:
1081 //
1082 Int_t EFILE = 0;
1083 printf(" Check the existence of CaloADC2MIP.raw file: \n\n");
1084 TFile tfile(file2.str().c_str());
1085 // TFile tfile2(file3);
1086 // if ( !tfile.IsZombie() && !tfile2.IsZombie()) {
1087 if ( !tfile.IsZombie() ) {
1088 printf(" The file already exists! Check if an update is needed \n\n");
1089 EFILE = 1;
1090 tfile.Close();
1091 } else {
1092 printf("\n OK, I will create it!\n\n");
1093 EFILE = 0;
1094 };
1095 //
1096 // the check for an update is done once even if we have a list.
1097 //
1098 TTree *tree;
1099 if ( EFILE == 1) {
1100 //
1101 // The file exists, open it (READ ONLY) and check if any update is needed:
1102 //
1103 TFile *hfile = 0;
1104 hfile = new TFile(file2.str().c_str(),"READ","Calorimeter CALIBRATION data");
1105 CalorimeterCalibration *calo = 0;
1106 tree = (TTree*)hfile->Get("CaloADC");
1107 tree->SetBranchAddress("Event", &calo);
1108 //
1109 Long64_t nevents = tree->GetEntries();
1110 printf("\n Number of updates: %i\n",(int)nevents);
1111 tree->GetEntry(nevents-1);
1112 //
1113 // Read the status variable of the last event: if it is 0 then we need an update, if it is one print out informations and exit
1114 //
1115 if ( calo->status == 0 ) {
1116 printf("\n An update is needed: \n");
1117 Int_t upstrip = 0;
1118 for (Int_t i = 0; i<2; i++){
1119 for (Int_t j = 0; j<22; j++){
1120 for (Int_t l = 0; l<96; l++){
1121 if ( calo->mask[i][j][l] == 0. ) upstrip++;
1122 };
1123 };
1124 };
1125 printf(" %i strip out of 4224 have no statistic enough.\n\n",upstrip);
1126 } else {
1127 printf(" No update is needed! \n The data files used to obtain the ADC to MIP conversion are: \n");
1128 for (Int_t i = 0; i<nevents; i++){
1129 tree->GetEntry(i);
1130 const char *fnout = calo->fname;
1131 printf(" - %s \n",fnout);
1132 };
1133 printf("\n\n");
1134 return;
1135 };
1136 hfile->Close();
1137 };
1138 //
1139 const char *filename2 = Filename;
1140 //
1141 TString lfile;
1142 TString sfile;
1143 //
1144 Int_t e = 0;
1145 ifstream in;
1146 in.open(flist, ios::in);
1147 char ctime[20];
1148 TString pfile;
1149 const char *cali;
1150 const char *ffile;
1151 Int_t calibex = 0;
1152 //
1153 // TString *file[17];
1154 TString filename;
1155 const char *nome = 0;
1156 TFile *hfile = 0;
1157 CalorimeterCalibration *calo = 0;
1158 Long64_t nevents;
1159 Int_t used = 0;
1160 const char *compare = 0;
1161 Int_t hf3 = 0;
1162 TFile *headerFile = 0;
1163 TFile *caloFile = 0;
1164 TFile *triggerFile = 0;
1165 TString file2f = "";
1166 pfile = "";
1167 const char *porca = 0;
1168 Int_t S4;
1169 Int_t upnn;
1170 // char *files = 0;
1171 stringstream files;
1172 TFile *hfileb = 0;
1173 //
1174 Int_t pl1, pl2, strmin, strmax;
1175 Int_t spre=0;
1176 Int_t surstrip = 0;
1177 Float_t estrippa = 0.;
1178 //
1179 TFile *hfile3 = 0;
1180 TString fififile;
1181 const char *file;
1182 //
1183 stringstream finome;
1184 finome.str("");
1185 while ( e < nfile ){
1186 if ( flist != "" ){
1187 in >> ctime;
1188 finome.str("");
1189 finome << filename2 << "/";
1190 finome << ctime;
1191 lfile = TString(finome.str().c_str());
1192 // file = "dw_000000_000.dat";
1193 //file = getFilename(lfile);
1194 fififile = getFilename(lfile);
1195 file = fififile;
1196 finome.str("");
1197 finome << file;
1198 sfile = TString(finome.str().c_str());
1199 } else {
1200 finome.str("");
1201 finome << filename2;
1202 lfile = TString(finome.str().c_str());
1203 //file = "dw_000000_000.dat";
1204 //file = getFilename(lfile);
1205 fififile = getFilename(lfile);
1206 file = fififile;
1207 finome.str("");
1208 finome << file;
1209 sfile = TString(finome.str().c_str());
1210 };
1211 //
1212 filename = lfile;
1213 nome = filename;
1214 printf("\n Processing file number %i: %s \n\n",e,nome);
1215 nome = sfile;
1216 //
1217 // this must be done for each file.
1218 //
1219 TTree *otr;
1220 pamela::calorimeter::CalorimeterEvent *de = 0;
1221 pamela::trigger::TriggerEvent *trigger = 0;
1222 pamela::PscuHeader *ph = 0;
1223 pamela::EventHeader *eh = 0;
1224 const char *calname;
1225 Int_t wused = 0;
1226 //
1227 TTree *tree;
1228 if ( EFILE == 1 ){
1229 hfile = new TFile(file2.str().c_str(),"UPDATE","Calorimeter CALIBRATION data");
1230 tree = (TTree*)hfile->Get("CaloADC");
1231 tree->SetBranchAddress("Event", &calo);
1232 //
1233 // Check if the input file has already been used
1234 //
1235 printf(" Check if file %s has already been used... \n",nome);
1236 nevents = tree->GetEntries();
1237 used = 0;
1238 for (Int_t i = 0; i<nevents; i++){
1239 tree->GetEntry(i);
1240 compare= calo->fname;
1241 if ( !strcmp(compare,nome) ) {
1242 used = 1;
1243 break;
1244 };
1245 };
1246 if ( used ) {
1247 printf(" The file %s has already been used!\n You cannot run twice this program with that file!\n\n",nome);
1248 hfile->Close();
1249 goto jumpfile;
1250 } else {
1251 printf(" Good, the file has not been used yet. \n\n");
1252 nevents = tree->GetEntries() -1;
1253 tree->GetEntry(nevents);
1254 };
1255 } else {
1256 //
1257 // If the file doesn't exist create it
1258 //
1259 CalorimeterCalibration *calo = new CalorimeterCalibration();
1260 hfile = new TFile(file2.str().c_str(),"NEW","Calorimeter CALIBRATION data");
1261 if ( hfile->IsZombie() ){
1262 printf(" ERROR: Cannot create file \n");
1263 return;
1264 };
1265 tree = new TTree("CaloADC","Calorimeter calibration data");
1266 tree->Branch("Event","CalorimeterCalibration",&calo);
1267 };
1268 //
1269 //
1270 //
1271 hf3 = 0;
1272 if ( EFILE == 1 && e == 0 ){
1273 printf(" Reading the old MIP figures to be updated...\n ");
1274 hfile3 = new TFile(file3.str().c_str());
1275 hf3 = 1;
1276 //
1277 stringstream fmippa;
1278 fmippa.str("");
1279 for (Int_t l = 0; l < 2; l++){
1280 for (Int_t m = 0; m < 22; m++){
1281 for (Int_t n =0 ; n < 96; n++){
1282 //char *name = Form("fmip %i %i %i",l,m,n);
1283 fmippa.str("");
1284 fmippa << "fmip " << l;
1285 fmippa << " " << m;
1286 fmippa << " " << n;
1287 TH1F *temp = (TH1F*)hfile3->Get(fmippa.str().c_str());
1288 pfmip[l][m][n] = (TH1F*)temp->Clone();
1289 };
1290 };
1291 };
1292 printf(" ...done!\n");
1293 };
1294 EFILE = 1;
1295 //
1296 // Check if the input calibration file contains any calibration data
1297 //
1298 done = 0;
1299 printf(" Check for calorimeter calibration: \n");
1300 calcalib = filename;
1301 findcalibs: calname = calcalib;
1302 printf(" Calibration file: %s \n",calname);
1303 for (Int_t s=0; s<4;s++){
1304 for (Int_t d = 0; d<50; d++){
1305 calib.ttime[s][d] = 0 ;
1306 if ( d < 49 ) calib.time[s][d] = 0 ;
1307 };
1308 };
1309 for (Int_t m = 0; m < 2 ; m++ ){
1310 for (Int_t k = 0; k < 22; k++ ){
1311 for (Int_t l = 0; l < 96; l++ ){
1312 calib.calped[m][k][l] = 0. ;
1313 evento.dexy[m][k][l] = 0. ;
1314 evento.dexyc[m][k][l] = 0. ;
1315 edelta[m][k][l] = 0.;
1316 estrip[m][k][l] = 0.;
1317 estripnull[m][k][l] = 0.;
1318 };
1319 };
1320 }
1321 wused = 0;
1322 CaloFindCalibs(calcalib, calcalib, wused, calib);
1323 // CaloFindCalibs(calcalib, calib);
1324 //
1325 // print on the screen the results:
1326 //
1327 ffile = filename;
1328 printf(" ------ %s ------- \n \n",ffile);
1329 calibex = 0;
1330 for (Int_t s=0; s<4;s++){
1331 printf(" ** SECTION %i **\n",s);
1332 for (Int_t d = 0; d<51; d++){
1333 if ( calib.ttime[s][d] != 0 ) {
1334 calibex++;
1335 if ( calib.fcode[s][d] != 10 ){
1336 file2f = "";
1337 stringcopy(file2f,calcalib);
1338 pfile = (TString)whatnamewith(file2f,calib.fcode[s][d]);
1339 } else {
1340 pfile = (TString)calcalib;
1341 };
1342 ffile = pfile;
1343 printf(" - from time %i to time %i use calibration at\n time %i, file: %s \n",calib.time[s][d],calib.time[s][d+1],calib.ttime[s][d],ffile);
1344 };
1345 };
1346 printf("\n");
1347 };
1348 printf(" ----------------------------------------------------------------------- \n \n");
1349 //
1350 if ( calibex < 1 ) {
1351 printf("No calibration data in this file!\n");
1352 calcalib = calcalibfile;
1353 if ( !done && calcalibfile != "" ) {
1354 cali = calcalibfile;
1355 printf(" Search in the calibration file %s \n",cali);
1356 done = 1;
1357 goto findcalibs;
1358 } else {
1359 printf(" Cannot find any calibration for this file! \n");
1360 goto jumpfile;
1361 };
1362 };
1363 if ( calibex < 4 ) {
1364 printf(" WARNING! No full calibration data in this file! \n\n");
1365 };
1366 //
1367 // upgrade data...
1368 //
1369 headerFile = emigetFile(filename, "Physics", "Header");
1370 if ( !headerFile ){
1371 printf("\n WARNING: No header file!\n");
1372 goto jumpfile;
1373 };
1374 caloFile = emigetFile(filename, "Physics", "Calorimeter");
1375 if ( !caloFile ){
1376 printf("\n WARNING: No calorimeter file! \n\n");
1377 headerFile->Close();
1378 goto jumpfile;
1379 };
1380 triggerFile = emigetFile(filename, "Trigger");
1381 if ( !triggerFile ){
1382 printf("ERROR: no Trigger file! \n");
1383 headerFile->Close();
1384 caloFile->Close();
1385 goto jumpfile;
1386 }
1387 //Takes the tree of the header file
1388 otr = (TTree*)headerFile->Get("Pscu");
1389 otr->AddFriend("Calorimeter",caloFile);
1390 otr->AddFriend("Trigger",triggerFile);
1391 //
1392 otr->SetBranchAddress("Header", &eh);
1393 otr->SetBranchAddress("Event", &de);
1394 otr->SetBranchAddress("Trigger.Event", &trigger);
1395 //
1396 nevents = otr->GetEntries();
1397 //
1398 // calibrate before starting
1399 //
1400 for (Int_t s = 0; s < 4; s++){
1401 b[s]=0;
1402 if ( calib.fcode[s][b[s]] != 10 ){
1403 stringcopy(file2f,calcalib);
1404 pfile = (TString)whatnamewith(file2f,calib.fcode[s][b[s]]);
1405 } else {
1406 pfile = (TString)calcalib;
1407 };
1408 porca = pfile;
1409 if ( debug ) printf(" porca miseria %s \n",porca);
1410 CaloPede(pfile,s,calib.ttime[s][b[s]],calib);
1411 };
1412 //
1413 // run over all the events
1414 //
1415 printf("\n Processed events: \n\n");
1416 //
1417 for (Int_t i = 0; i < nevents; i++){
1418 evno = i;
1419 if ( i%1000 == 0 && i > 0 ) printf(" %iK \n",i/1000);
1420 //
1421 // read from the header of the event the absolute time at which it was recorded
1422 //
1423 otr->GetEntry(i);
1424 //
1425 S4 = trigger->patterntrig[1] & (1<<0);
1426 //
1427 if ( !S4 || i == 0 ) {
1428 ph = eh->GetPscuHeader();
1429 etime = ph->GetOrbitalTime();
1430 //
1431 // for each event check that the calibration we are using are still within calibration limits, if not call the next calibration
1432 //
1433 if ( !calib.obtjump ) {
1434 for (Int_t s = 0; s < 4; s++){
1435 if ( calib.ttime[s][b[s]] ){
1436 while ( etime > calib.time[s][b[s]+1] ){
1437 printf(" CALORIMETER: \n" );
1438 printf(" - Section %i, event at time %i while old calibration time limit at %i. Use new calibration at time %i -\n",s,etime,calib.time[s][b[s]+1],calib.ttime[s][b[s]+1]);
1439 printf(" END CALORIMETER. \n\n" );
1440 b[s]++;
1441 if ( calib.fcode[s][b[s]] != 10 ){
1442 stringcopy(file2f,calcalib);
1443 pfile = (TString)whatnamewith(file2f,calib.fcode[s][b[s]]);
1444 } else {
1445 pfile = (TString)calcalib;
1446 };
1447 porca = pfile;
1448 if ( debug ) printf(" porca miseria %s \n",porca);
1449 CaloPede(pfile,s,calib.ttime[s][b[s]],calib);
1450 };
1451 };
1452 };
1453 };
1454 //
1455 // do whatever you want with data
1456 //
1457 evento.iev = de->iev;
1458 //
1459 // run over views and planes
1460 //
1461 for (Int_t l = 0; l < 2; l++){
1462 for (Int_t m = 0; m < 22; m++){
1463 //
1464 // determine the section number
1465 //
1466 if (l == 0 && m%2 == 0) se = 3;
1467 if (l == 0 && m%2 != 0) se = 2;
1468 if (l == 1 && m%2 == 0) se = 0;
1469 if (l == 1 && m%2 != 0) se = 1;
1470 //
1471 // determine what kind of event we are going to analyze
1472 //
1473 isCOMP = 0;
1474 isFULL = 0;
1475 isRAW = 0;
1476 if ( de->stwerr[se] & (1 << 16) ) isCOMP = 1;
1477 if ( de->stwerr[se] & (1 << 17) ) isFULL = 1;
1478 if ( de->stwerr[se] & (1 << 3) ) isRAW = 1;
1479 //
1480 // save the prevoius energy deposit and calibration in sbase, sdexy, sdexyc
1481 //
1482 pre = -1;
1483 if ( isRAW ){
1484 for (Int_t nn = 0; nn < 96; nn++){
1485 if ( nn%16 == 0 ) pre++;
1486 evento.base[l][m][pre] = calib.calbase[l][m][pre];
1487 sdexy[l][m][nn] = evento.dexy[l][m][nn];
1488 evento.dexy[l][m][nn] = de->dexy[l][m][nn] ;
1489 sdexyc[l][m][nn] = evento.dexyc[l][m][nn];
1490 evento.dexyc[l][m][nn] = de->dexyc[l][m][nn] ;
1491 };
1492 //
1493 // run over strips
1494 //
1495 pre = -1;
1496 for (Int_t n =0 ; n < 96; n++){
1497 if ( n%16 == 0 ) {
1498 pre++;
1499 rdone = 0;
1500 done = 0;
1501 };
1502 if ( calo->mask[l][m][n] == 0. && ( isRAW || ( de->dexyc[l][m][n] > 0. && de->dexyc[l][m][n] < 32000.)) ){
1503 //
1504 // baseline check and calculation
1505 //
1506 if ( !rdone ){
1507 CaloFindBaseRawNC(calib,evento,l,m,pre);
1508 rdone = 1;
1509 };
1510 //
1511 // no suitable new baseline, use old ones and compress data!
1512 //
1513 if ( !done && (evento.base[l][m][pre] > 30000. || evento.base[l][m][pre] == 0.) ){
1514 evento.base[l][m][pre] = calib.sbase[l][m][pre];
1515 upnn = n+16;
1516 if ( upnn > 96 ) upnn = 96;
1517 for ( Int_t nn = n; nn<upnn; nn++ ){
1518 evento.dexyc[l][m][nn] = de->dexyc[l][m][nn] ;
1519 };
1520 done = 1;
1521 };
1522 //
1523 // CALIBRATION ALGORITHM
1524 //
1525 pl1 = -1;
1526 pl2 = -1;
1527 pl1 = m - 1;
1528 pl2 = m + 1;
1529 if ( pl1 < 0 ) {
1530 pl1 = pl2;
1531 pl2 = m + 2;
1532 };
1533 if ( pl2 > 21 ) {
1534 pl2 = pl1;
1535 pl1 = m - 2;
1536 };
1537 strmin = n - 2;
1538 if ( strmin < 0 ) strmin = 0;
1539 strmax = n + 2;
1540 if ( strmax > 95 ) strmax = 95;
1541 surstrip = 0;
1542 for (Int_t spl = pl1; spl<pl2+1; spl++){
1543 if ( spl != m ){
1544 for (Int_t sst = strmin; sst<strmax+1; sst++){
1545 if ( sst < 16 ) spre = 0;
1546 if ( 15 < sst && sst < 32 ) spre = 2;
1547 if ( 31 < sst && sst < 48 ) spre = 3;
1548 if ( 47 < sst && sst < 64 ) spre = 4;
1549 if ( 63 < sst && sst < 80 ) spre = 5;
1550 if ( 79 < sst ) spre = 6;
1551 if ( isRAW ){
1552 estrippa = ( evento.dexy[l][spl][sst] - calib.calped[l][spl][sst] - calib.sbase[l][spl][spre])/26. ;
1553 } else {
1554 estrippa = ( evento.dexyc[l][spl][sst] - calib.calped[l][spl][sst] - calib.sbase[l][spl][spre])/26. ;
1555 };
1556 if ( estrippa > 0.7 && estrippa < 10. && calib.calgood[l][spl][sst] == 0 ){
1557 surstrip++;
1558 //printf("surstrip = %i estrip = %f \n",surstrip,estrip);
1559 };
1560 };
1561 };
1562 };
1563 //
1564 //
1565 estrip[l][m][n] = ( evento.dexyc[l][m][n] - calib.calped[l][m][n] - evento.base[l][m][pre]) ;
1566 pfmip[l][m][n]->Fill(estrip[l][m][n]);
1567 calib.sbase[l][m][pre] = evento.base[l][m][pre];
1568 };
1569 };
1570 };
1571 };
1572 };
1573 };
1574 };
1575 //
1576 // save filename
1577 //
1578 calo->fname = sfile;
1579 tree->Fill();
1580 hfile->Write();
1581 fileup++;
1582 hfile->Close();
1583 //
1584 // create a backup copy!
1585 //
1586 if ( nfile > 1 && ( e%10 == 0 || e == nfile-1 ) && e > 0 ){
1587 files.str("");
1588 files << file4.str().c_str() << (e+1);
1589 files << ".bak";
1590 // files = Form("%s%i.bak",file4,e+1);
1591 printf(" Save the backup file with figures:\n %s \n ",files.str().c_str());
1592 hfileb = new TFile(files.str().c_str(),"RECREATE","Calorimeter CALIBRATION figures");
1593 if ( hfileb->IsZombie() ){
1594 printf(" ERROR: Cannot create file!!! \n");
1595 goto oltre;
1596 };
1597 TH1F *bfmip[2][22][96];
1598 stringstream bmippa;
1599 for (Int_t l = 0; l < 2; l++){
1600 for (Int_t m = 0; m < 22; m++){
1601 for (Int_t n =0 ; n < 96; n++){
1602 bmippa.str("");
1603 bmippa << "bmip " << l;
1604 bmippa << " " << m;
1605 bmippa << " " << n;
1606 bfmip[l][m][n] = new TH1F(bmippa.str().c_str(),"",70,-15.,55.);
1607 bfmip[l][m][n] = (TH1F*)pfmip[l][m][n]->Clone(bmippa.str().c_str());
1608 bfmip[l][m][n]->Write();
1609 };
1610 };
1611 };
1612 hfileb->Close();
1613 };
1614 //
1615 //
1616 //
1617 headerFile->Close();
1618 caloFile->Close();
1619 triggerFile->Close();
1620 goto oltre;
1621 jumpfile: printf(" File not processed! \n ");
1622 hfile->Close();
1623 oltre:
1624 // gObjectTable->Print();
1625 e++;
1626 };
1627 //
1628 // Save histograms and fit figures
1629 //
1630 TFile *hfile2;
1631 hfile2 = new TFile(file3.str().c_str(),"RECREATE","Calorimeter CALIBRATION figures");
1632 if ( hfile2->IsZombie() ){
1633 printf(" ERROR: Cannot create file \n");
1634 return;
1635 };
1636 TH1F *sfmip[2][22][96];
1637 stringstream fmippa;
1638 for (Int_t l = 0; l < 2; l++){
1639 for (Int_t m = 0; m < 22; m++){
1640 for (Int_t n =0 ; n < 96; n++){
1641 fmippa.str("");
1642 fmippa << "fmip " << l;
1643 fmippa << " " << m;
1644 fmippa << " " << n;
1645 sfmip[l][m][n] = new TH1F(fmippa.str().c_str(),"",72,-15.,55.);
1646 sfmip[l][m][n] = (TH1F*)pfmip[l][m][n]->Clone(fmippa.str().c_str());
1647 sfmip[l][m][n]->Write();
1648 };
1649 };
1650 };
1651 hfile2->Close();
1652 if ( hf3 ) hfile3->Close();
1653 }
1654
1655 void CaloRAWADC2MIPDATA(TString Filename, TString calcalibfile = "", TString flist=""){
1656 gROOT->GetListOfCanvases()->Delete("a");
1657 gROOT->Reset("a");
1658 gDirectory->GetList()->Delete();
1659 #if defined (__CINT__)
1660 emicheckLib();
1661 stringstream llib;
1662 const char *pamlib = gSystem->ExpandPathName("$PAM_LIB");
1663 llib << pamlib << "/caloclasses_h.so";
1664 gSystem->Load(llib.str().c_str());
1665 stringstream lolib;
1666 lolib.str("");
1667 lolib << pamlib << "/UTpathtoc_c.so";
1668 gSystem->Load(lolib.str().c_str());
1669 #endif
1670 TString calcalib;
1671 struct Evento evento;
1672 struct Calib calib;
1673 Int_t debug = 0;
1674 Int_t b[4];
1675 Int_t etime,evno;
1676 Float_t sdexy[2][22][96],sdexyc[2][22][96];
1677 Float_t edelta[2][22][96];
1678 Float_t estrip[2][22][96];
1679 Float_t estripnull[2][22][96];
1680 Float_t delta;
1681 delta = 1.;
1682 evento.emin = 0.7;
1683 Int_t se =5;
1684 Int_t pre;
1685 Int_t rdone = 0;
1686 Int_t done;
1687 bool isCOMP = 0;
1688 bool isFULL = 0;
1689 bool isRAW = 0;
1690 //
1691 // the file which contains the calorimeter ADC to MIP conversion values is file2 = path to the script dir / CaloADC2MIP.root
1692 //
1693 const char *pam_calib = pathtocalibration();
1694 // const char *pam_calib=gSystem->Getenv("PAM_CALIB");
1695 stringstream file5;
1696 file5.str("");
1697 file5 << pam_calib << "/CaloADC2MIPdata.raw";
1698 //
1699 Int_t nfile = 0;
1700 Int_t fileup = 0;
1701 if ( flist != "" ){
1702 printf("\n You have given an input file list \n\n Files: \n");
1703 ifstream in;
1704 in.open(flist, ios::in);
1705 char ctime[20];
1706 while (1) {
1707 in >> ctime;
1708 if (!in.good()) break;
1709 printf(" File number %i = %s \n",nfile+1,ctime);
1710 nfile++;
1711 // if ( nfile > 499 ) {
1712 // printf(" No more than 500 files please! \n ");
1713 // return;
1714 //};
1715 };
1716 in.close();
1717 if ( !nfile ) {
1718 printf("\n\n The list is empty, exiting... \n\n");
1719 return;
1720 };
1721 } else {
1722 nfile = 1;
1723 };
1724 printf("\n");
1725 //
1726 // First check if the calibration file exists:
1727 //
1728 //
1729 const char *filename2 = Filename;
1730 //
1731 TString lfile;
1732 TString sfile;
1733 //
1734 Int_t e = 0;
1735 ifstream in;
1736 in.open(flist, ios::in);
1737 char ctime[20];
1738 TString file2f;
1739 TString pfile;
1740 const char *cali;
1741 const char *ffile;
1742 Int_t calibex = 0;
1743 //
1744 // TString *file[17];
1745 TString filename;
1746 const char *nome = 0;
1747 Long64_t nevents;
1748 TFile *headerFile = 0;
1749 TFile *caloFile = 0;
1750 file2f = "";
1751 pfile = "";
1752 const char *porca = 0;
1753 Int_t upnn;
1754 //
1755 //
1756 TFile *hfl1 = 0;
1757 hfl1 = new TFile(file5.str().c_str(),"RECREATE","Calorimeter LEVEL1 data");
1758 //
1759 // class CalorimeterLevel1
1760 //
1761 CalorimeterADCRAW *clev1 = new CalorimeterADCRAW();
1762 TTree *trl1 = 0;
1763 trl1 = new TTree("CaloADCRAW","PAMELA Level1 calorimeter data");
1764 trl1->Branch("Rawdata","CalorimeterADCRAW",&clev1);
1765 //
1766 Int_t glbc = 0;
1767 //
1768 TString fififile;
1769 const char *file;
1770 //
1771 stringstream lfillo;
1772 while ( e < nfile ){
1773 if ( flist != "" ){
1774 in >> ctime;
1775 lfillo.str("");
1776 lfillo << filename2 << "/";
1777 lfillo << ctime;
1778 lfile = TString(lfillo.str().c_str());
1779 //file = "dw_000000_000.dat";
1780 //file = getFilename(lfile);
1781 fififile = getFilename(lfile);
1782 file = fififile;
1783 lfillo.str("");
1784 lfillo << file;
1785 sfile = TString(lfillo.str().c_str());
1786 } else {
1787 lfillo.str("");
1788 lfillo << filename2;
1789 lfile = TString(lfillo.str().c_str());
1790 //file = "dw_000000_000.dat";
1791 //file = getFilename(lfile);
1792 fififile = getFilename(lfile);
1793 file = fififile;
1794 lfillo.str("");
1795 lfillo << file;
1796 sfile = TString(lfillo.str().c_str());
1797 };
1798 //
1799 filename = lfile;
1800 nome = filename;
1801 printf("\n Processing file number %i: %s \n\n",e,nome);
1802 nome = sfile;
1803 //
1804 // Check if the input calibration file contains any calibration data
1805 //
1806 TTree *otr;
1807 pamela::calorimeter::CalorimeterEvent *de = 0;
1808 pamela::PscuHeader *ph = 0;
1809 pamela::EventHeader *eh = 0;
1810 //
1811 done = 0;
1812 printf(" Check for calorimeter calibration: \n");
1813 calcalib = filename;
1814 findcalibs: const char *calname = calcalib;
1815 printf(" Calibration file: %s \n",calname);
1816 for (Int_t s=0; s<4;s++){
1817 for (Int_t d = 0; d<50; d++){
1818 calib.ttime[s][d] = 0 ;
1819 if ( d < 49 ) calib.time[s][d] = 0 ;
1820 };
1821 };
1822 for (Int_t m = 0; m < 2 ; m++ ){
1823 for (Int_t k = 0; k < 22; k++ ){
1824 for (Int_t l = 0; l < 96; l++ ){
1825 calib.calped[m][k][l] = 0. ;
1826 evento.dexy[m][k][l] = 0. ;
1827 evento.dexyc[m][k][l] = 0. ;
1828 edelta[m][k][l] = 0.;
1829 estrip[m][k][l] = 0.;
1830 estripnull[m][k][l] = 0.;
1831 };
1832 };
1833 }
1834 Int_t wused = 0;
1835 CaloFindCalibs(filename, calcalibfile, wused, calib);
1836 if ( wused == 1 ) calcalibfile = filename;
1837 //
1838 // print on the screen the results:
1839 //
1840 ffile = filename;
1841 printf(" ------ %s ------- \n \n",ffile);
1842 calibex = 0;
1843 for (Int_t s=0; s<4;s++){
1844 printf(" ** SECTION %i **\n",s);
1845 for (Int_t d = 0; d<51; d++){
1846 if ( calib.ttime[s][d] != 0 ) {
1847 calibex++;
1848 if ( calib.fcode[s][d] != 10 ){
1849 file2f = "";
1850 stringcopy(file2f,calcalib);
1851 pfile = (TString)whatnamewith(file2f,calib.fcode[s][d]);
1852 } else {
1853 pfile = (TString)calcalib;
1854 };
1855 ffile = pfile;
1856 printf(" - from time %i to time %i use calibration at\n time %i, file: %s \n",calib.time[s][d],calib.time[s][d+1],calib.ttime[s][d],ffile);
1857 };
1858 };
1859 printf("\n");
1860 };
1861 printf(" ----------------------------------------------------------------------- \n \n");
1862 //
1863 if ( calibex < 1 ) {
1864 printf("No calibration data in this file!\n");
1865 calcalib = calcalibfile;
1866 if ( !done && calcalibfile != "" ) {
1867 cali = calcalibfile;
1868 printf(" Search in the calibration file %s \n",cali);
1869 done = 1;
1870 goto findcalibs;
1871 } else {
1872 printf(" Cannot find any calibration for this file! \n");
1873 goto jumpfile;
1874 };
1875 };
1876 if ( calibex < 4 ) {
1877 printf(" WARNING! No full calibration data in this file! \n\n");
1878 };
1879 //
1880 // upgrade data...
1881 //
1882 headerFile = emigetFile(filename, "Physics", "Header");
1883 if ( !headerFile ){
1884 printf("\n WARNING: No header file!\n");
1885 goto jumpfile;
1886 };
1887 caloFile = emigetFile(filename, "Physics", "Calorimeter");
1888 if ( !caloFile ){
1889 printf("\n WARNING: No calorimeter file! \n\n");
1890 headerFile->Close();
1891 goto jumpfile;
1892 };
1893 //Takes the tree of the header file
1894 otr = (TTree*)headerFile->Get("Pscu");
1895 otr->AddFriend("Calorimeter",caloFile);
1896 //
1897 otr->SetBranchAddress("Header", &eh);
1898 otr->SetBranchAddress("Event", &de);
1899 //
1900 nevents = otr->GetEntries();
1901 //
1902 // calibrate before starting
1903 //
1904 for (Int_t s = 0; s < 4; s++){
1905 b[s]=0;
1906 if ( calib.fcode[s][b[s]] != 10 ){
1907 stringcopy(file2f,calcalib);
1908 pfile = (TString)whatnamewith(file2f,calib.fcode[s][b[s]]);
1909 } else {
1910 pfile = (TString)calcalib;
1911 };
1912 porca = pfile;
1913 if ( debug ) printf(" porca miseria %s \n",porca);
1914 CaloPede(pfile,s,calib.ttime[s][b[s]],calib);
1915 };
1916 //
1917 // run over all the events
1918 //
1919 printf("\n Processed events: \n\n");
1920 //
1921 Int_t updata;
1922 for (Int_t i = 0; i < nevents; i++){
1923 evno = i;
1924 if ( i%1000 == 0 && i > 0 ) printf(" %iK \n",i/1000);
1925 //
1926 // read from the header of the event the absolute time at which it was recorded
1927 //
1928 otr->GetEntry(i);
1929 //
1930 //
1931 ph = eh->GetPscuHeader();
1932 etime = ph->GetOrbitalTime();
1933 //
1934 // for each event check that the calibration we are using are still within calibration limits, if not call the next calibration
1935 //
1936 if ( !calib.obtjump ) {
1937 for (Int_t s = 0; s < 4; s++){
1938 if ( calib.ttime[s][b[s]] ){
1939 while ( etime > calib.time[s][b[s]+1] ){
1940 printf(" CALORIMETER: \n" );
1941 printf(" - Section %i, event at time %i while old calibration time limit at %i. Use new calibration at time %i -\n",s,etime,calib.time[s][b[s]+1],calib.ttime[s][b[s]+1]);
1942 printf(" END CALORIMETER. \n\n" );
1943 b[s]++;
1944 if ( calib.fcode[s][b[s]] != 10 ){
1945 stringcopy(file2f,calcalib);
1946 pfile = (TString)whatnamewith(file2f,calib.fcode[s][b[s]]);
1947 } else {
1948 pfile = (TString)calcalib;
1949 };
1950 porca = pfile;
1951 if ( debug ) printf(" porca miseria %s \n",porca);
1952 CaloPede(pfile,s,calib.ttime[s][b[s]],calib);
1953 };
1954 };
1955 };
1956 };
1957 //
1958 // do whatever you want with data
1959 //
1960 evento.iev = de->iev;
1961 //
1962 // run over views and planes
1963 //
1964 updata = 0;
1965 glbc++;
1966 for (Int_t l = 0; l < 2; l++){
1967 for (Int_t m = 0; m < 22; m++){
1968 //
1969 // determine the section number
1970 //
1971 if (l == 0 && m%2 == 0) se = 3;
1972 if (l == 0 && m%2 != 0) se = 2;
1973 if (l == 1 && m%2 == 0) se = 0;
1974 if (l == 1 && m%2 != 0) se = 1;
1975 //
1976 // determine what kind of event we are going to analyze
1977 //
1978 isCOMP = 0;
1979 isFULL = 0;
1980 isRAW = 0;
1981 if ( de->stwerr[se] & (1 << 16) ) isCOMP = 1;
1982 if ( de->stwerr[se] & (1 << 17) ) isFULL = 1;
1983 if ( de->stwerr[se] & (1 << 3) ) isRAW = 1;
1984 //
1985 // save the prevoius energy deposit and calibration in sbase, sdexy, sdexyc
1986 //
1987 pre = -1;
1988 if ( isRAW || isFULL ){
1989 //clev1 = new CalorimeterADCRAW();
1990 updata = 1;
1991 for (Int_t nn = 0; nn < 96; nn++){
1992 if ( nn%16 == 0 ) pre++;
1993 evento.base[l][m][pre] = calib.calbase[l][m][pre];
1994 sdexy[l][m][nn] = evento.dexy[l][m][nn];
1995 evento.dexy[l][m][nn] = de->dexy[l][m][nn] ;
1996 sdexyc[l][m][nn] = evento.dexyc[l][m][nn];
1997 // !!!
1998 evento.dexyc[l][m][nn] = de->dexy[l][m][nn] ;
1999 // !!!
2000 };
2001 //
2002 // run over strips
2003 //
2004 pre = -1;
2005 for (Int_t n =0 ; n < 96; n++){
2006 if ( n%16 == 0 ) {
2007 pre++;
2008 rdone = 0;
2009 done = 0;
2010 };
2011 if ( isRAW || ( de->dexyc[l][m][n] > 0. && de->dexyc[l][m][n] < 32000.) || isFULL ){
2012 if ( isFULL && !rdone ) rdone = 1;
2013 //
2014 // baseline check and calculation
2015 //
2016 if ( !rdone ){
2017 CaloFindBaseRawNC(calib,evento,l,m,pre);
2018 rdone = 1;
2019 };
2020 //
2021 // no suitable new baseline, use old ones and compress data!
2022 //
2023 if ( !done && (evento.base[l][m][pre] > 30000. || evento.base[l][m][pre] == 0.) ){
2024 evento.base[l][m][pre] = calib.sbase[l][m][pre];
2025 upnn = n+16;
2026 if ( upnn > 96 ) upnn = 96;
2027 for ( Int_t nn = n; nn<upnn; nn++ ){
2028 evento.dexyc[l][m][nn] = de->dexy[l][m][nn] ;
2029 };
2030 done = 1;
2031 };
2032 //
2033 // CALIBRATION ALGORITHM
2034 //
2035 estrip[l][m][n] = ( evento.dexyc[l][m][n] - calib.calped[l][m][n] - evento.base[l][m][pre]) ;
2036 calib.sbase[l][m][pre] = evento.base[l][m][pre];
2037 clev1->diffbas[l][m][pre] = evento.base[l][m][pre];
2038 clev1->estrip[l][m][n] = estrip[l][m][n];
2039 };
2040 };
2041 };
2042 };
2043 };
2044 // if ( updata ){
2045 clev1->evno = glbc;
2046 clev1->etime = etime;
2047 trl1->Fill();
2048 for (Int_t l =0 ; l < 2; l++){
2049 for (Int_t m =0 ; m < 22; m++){
2050 pre = -1;
2051 for (Int_t n =0 ; n < 96; n++){
2052 // memcpy(estrip, estripnull, sizeof(estrip));
2053 if ( n%16 == 0 ) pre++;
2054 estrip[l][m][n] = 0.;
2055 evento.base[l][m][pre] = 0.;
2056 clev1->diffbas[l][m][pre] = evento.base[l][m][pre];
2057 clev1->estrip[l][m][n] = estrip[l][m][n];
2058 };
2059 };
2060 };
2061 //gObjectTable->Print();
2062 //};
2063 };
2064 //
2065 // save filename
2066 //
2067 fileup++;
2068 //
2069 headerFile->Close();
2070 caloFile->Close();
2071 goto oltre;
2072 jumpfile: printf(" File not processed! \n ");
2073 oltre:
2074 e++;
2075 };
2076 hfl1->Write();
2077 hfl1->Close();
2078 }
2079
2080 void Calo4224STATUS(TString filename = "", TString style = ""){
2081 //
2082 // this routine shows the calorimeter calibration status with a figure. You can choose style between "error" (or "") to show the error on the
2083 // peak maximum determined and "status" to show which strips have converged to a certain value.
2084 //
2085 gROOT->GetListOfCanvases()->Delete("a");
2086 gROOT->Reset("a");
2087 gDirectory->GetList()->Delete();
2088 #if defined (__CINT__)
2089 emicheckLib();
2090 stringstream llib;
2091 const char *pamlib = gSystem->ExpandPathName("$PAM_LIB");
2092 llib << pamlib << "/caloclasses_h.so";
2093 gSystem->Load(llib.str().c_str());
2094 stringstream lolib;
2095 lolib.str("");
2096 lolib << pamlib << "/UTpathtoc_c.so";
2097 gSystem->Load(lolib.str().c_str());
2098 #endif
2099 if ( filename == "" ){
2100 const char *pam_calib = pathtocalibration();
2101 //const char *pam_calib=gSystem->Getenv("PAM_CALIB");
2102 stringstream filenome;
2103 filenome.str("");
2104 filenome << pam_calib << "/CaloADC2MIP.root";
2105 // filename = Form("%s/CaloADC2MIP.root",pam_calib);
2106 filename = (TString)filenome.str().c_str();
2107 };
2108 if ( style == "" ) style = "error";
2109 const char *file2 = filename;
2110 TFile *hfile;
2111 hfile = new TFile(file2,"READ","Calorimeter CALIBRATION data");
2112 CalorimeterCalibration *calo = 0;
2113 TTree *tree;
2114 tree = (TTree*)hfile->Get("CaloADC");
2115 tree->SetBranchAddress("Event", &calo);
2116 //
2117 Long64_t nevents = tree->GetEntries();
2118 tree->GetEntry(nevents-1);
2119 //
2120 // Book the histograms:
2121 //
2122 //
2123 TH2F *Xview = new TH2F("x-view event ","",96,-0.5,95.5,22,-0.5,21.5);
2124 TH2F *Yview = new TH2F("y-view event ","",96,-0.5,95.5,22,-0.5,21.5);
2125 //
2126 // figures:
2127 //
2128 TCanvas *figura2 = new TCanvas("Calorimeter:_strip_RMS", "Calorimeter:_strip_RMS", 750, 650);
2129 figura2->SetFillColor(10);
2130 figura2->Range(0,0,100,100);
2131 Int_t bgcolor = 10;
2132 TPad *pd1 = new TPad("pd1","This is pad1",0.02,0.05,0.88,0.49,bgcolor);
2133 TPad *pd2 = new TPad("pd2","This is pad2",0.02,0.51,0.88,0.95,bgcolor);
2134 TPad *palette = new TPad("palette","This is palette",0.90,0.05,0.98,0.90,bgcolor);
2135 figura2->cd();
2136 TLatex *t=new TLatex();
2137 t->SetTextFont(32);
2138 t->SetTextColor(1);
2139 t->SetTextSize(0.03);
2140 t->SetTextAlign(12);
2141 t->DrawLatex(90.,92.5,"error");
2142 pd1->Range(0,0,100,100);
2143 pd2->Range(0,0,100,100);
2144 palette->Range(0,0,101,100);
2145 pd1->SetTicks();
2146 pd2->SetTicks();
2147 pd1->Draw();
2148 pd2->Draw();
2149 palette->Draw();
2150 palette->cd();
2151 const char *style2 = style;
2152 if ( !strcmp(style2,"error") ){
2153 // palette
2154 TPaveLabel *box1 = new TPaveLabel(1,2,99,14,"no fit","");
2155 box1->SetTextFont(32);
2156 box1->SetTextColor(1);
2157 box1->SetTextSize(0.25);
2158 box1->SetFillColor(10);
2159 box1->Draw();
2160 TPaveLabel *box2 = new TPaveLabel(1,16,99,28,"< 0.15","");
2161 box2->SetTextFont(32);
2162 box2->SetTextColor(1);
2163 box2->SetTextSize(0.25);
2164 box2->SetFillColor(38);
2165 box2->Draw();
2166 TPaveLabel *box3 = new TPaveLabel(1,30,99,42,"0.15 - 0.25","");
2167 box3->SetTextFont(32);
2168 box3->SetTextColor(1);
2169 box3->SetTextSize(0.2);
2170 box3->SetFillColor(4);
2171 box3->Draw();
2172 TPaveLabel *box4 = new TPaveLabel(1,44,99,56,"0.25 - 0.55","");
2173 box4->SetTextFont(32);
2174 box4->SetTextColor(1);
2175 box4->SetTextSize(0.2);
2176 box4->SetFillColor(3);
2177 box4->Draw();
2178 TPaveLabel *box5 = new TPaveLabel(1,58,99,70,"0.55 - 1","");
2179 box5->SetTextFont(32);
2180 box5->SetTextColor(1);
2181 box5->SetTextSize(0.25);
2182 box5->SetFillColor(2);
2183 box5->Draw();
2184 TPaveLabel *box6 = new TPaveLabel(1,72,99,84,"1 - 5","");
2185 box6->SetTextFont(32);
2186 box6->SetTextColor(1);
2187 box6->SetTextSize(0.25);
2188 box6->SetFillColor(6);
2189 box6->Draw();
2190 TPaveLabel *box7 = new TPaveLabel(1,86,99,98,"> 5","");
2191 box7->SetTextFont(32);
2192 box7->SetTextColor(10);
2193 box7->SetTextSize(0.25);
2194 box7->SetFillColor(1);
2195 box7->Draw();
2196 } else {
2197 TPaveLabel *box1 = new TPaveLabel(1,2,99,14,"not done","");
2198 box1->SetTextFont(32);
2199 box1->SetTextColor(1);
2200 box1->SetTextSize(0.25);
2201 box1->SetFillColor(10);
2202 box1->Draw();
2203 TPaveLabel *box2 = new TPaveLabel(1,16,99,28,"done","");
2204 box2->SetTextFont(32);
2205 box2->SetTextColor(1);
2206 box2->SetTextSize(0.25);
2207 box2->SetFillColor(38);
2208 box2->Draw();
2209 };
2210 //
2211 pd1->cd();
2212 gStyle->SetOptStat(0);
2213 Xview->SetXTitle("strip");
2214 Xview->SetYTitle("X - plane");
2215 Xview->GetYaxis()->SetTitleOffset(0.5);
2216 Xview->SetFillColor(bgcolor);
2217 Xview->Fill(1.,1.,1.);
2218 Xview->Draw("box");
2219 pd1->Update();
2220 pd2->cd();
2221 gStyle->SetOptStat(0);
2222 Yview->SetXTitle("strip");
2223 Yview->SetYTitle("Y - plane");
2224 Yview->GetYaxis()->SetTitleOffset(0.5);
2225 Yview->SetFillColor(bgcolor);
2226 Yview->Fill(1.,1.,1.);
2227 Yview->Draw("box");
2228 pd2->Update();
2229 //
2230 // run over views and planes
2231 //
2232 stringstream xview;
2233 stringstream yview;
2234 for (Int_t m = 0; m < 22; m++){
2235 for (Int_t l = 0; l < 2; l++){
2236 for (Int_t n = 0; n < 96; n++){
2237 xview.str("");
2238 xview << "x-view event " << n;
2239 xview << " " << m;
2240 xview << " " << l;
2241 yview.str("");
2242 yview << "y-view event " << n;
2243 yview << " " << m;
2244 yview << " " << l;
2245 gDirectory->Delete(xview.str().c_str());
2246 gDirectory->Delete(yview.str().c_str());
2247 TH2F *Xview = new TH2F(xview.str().c_str(),"",96,-0.5,95.5,22,-0.5,21.5);
2248 TH2F *Yview = new TH2F(yview.str().c_str(),"",96,-0.5,95.5,22,-0.5,21.5);
2249 if ( !strcmp(style2,"error") ){
2250 if ( calo->fpe[1][l][m][n] < 0.15 ){
2251 Xview->SetFillColor(38);
2252 Yview->SetFillColor(38);
2253 };
2254 if ( calo->fpe[1][l][m][n] >= 0.15 ){
2255 Xview->SetFillColor(4);
2256 Yview->SetFillColor(4);
2257 };
2258 if ( calo->fpe[1][l][m][n] >= 0.25){
2259 Xview->SetFillColor(3);
2260 Yview->SetFillColor(3);
2261 };
2262 if ( calo->fpe[1][l][m][n] >= 0.55 ){
2263 Xview->SetFillColor(2);
2264 Yview->SetFillColor(2);
2265 };
2266 if ( calo->fpe[1][l][m][n] >= 1. ){
2267 Xview->SetFillColor(6);
2268 Yview->SetFillColor(6);
2269 };
2270 if ( calo->fpe[1][l][m][n] >= 5. ){
2271 Xview->SetFillColor(1);
2272 Yview->SetFillColor(1);
2273 };
2274 if ( calo->mip[l][m][n] < 15. || calo->mip[l][m][n] > 40. ){
2275 Xview->SetFillColor(10);
2276 Yview->SetFillColor(10);
2277 };
2278 } else {
2279 if ( calo->mask[l][m][n] == 1 ){
2280 Xview->SetFillColor(38);
2281 Yview->SetFillColor(38);
2282 } else {
2283 Xview->SetFillColor(10);
2284 Yview->SetFillColor(10);
2285 };
2286 };
2287 if ( l == 0 ) {
2288 Xview->Fill(n,m,1.);
2289 pd1->cd();
2290 Xview->Draw("box same");
2291 };
2292 if ( l == 1 ) {
2293 Yview->Fill(n,m,1.);
2294 pd2->cd();
2295 Yview->Draw("box same");
2296 };
2297 };
2298 };
2299 };
2300 figura2->cd();
2301 gStyle->SetOptDate(1);
2302 //
2303 t=new TLatex();
2304 t->SetTextFont(32);
2305 t->SetTextColor(8);
2306 t->SetTextAlign(12);
2307 t->SetTextSize(0.02);
2308 figura2->Update();
2309 pd1->Update();
2310 pd2->Update();
2311 //
2312 }
2313
2314 void Calo4224FIT(TString filename = "", TString filevalue = "", TString type=""){
2315 //
2316 // this routine shows the 4224 figures with fit from the main program
2317 //
2318 gROOT->GetListOfCanvases()->Delete("a");
2319 gROOT->Reset("a");
2320 gDirectory->GetList()->Delete();
2321 const char* startingdir = gSystem->WorkingDirectory();
2322 #if defined (__CINT__)
2323 emicheckLib();
2324 stringstream llib;
2325 const char *pamlib = gSystem->ExpandPathName("$PAM_LIB");
2326 llib << pamlib << "/caloclasses_h.so";
2327 gSystem->Load(llib.str().c_str());
2328 stringstream lolib;
2329 lolib.str("");
2330 lolib << pamlib << "/UTpathtoc_c.so";
2331 gSystem->Load(lolib.str().c_str());
2332 #endif
2333 if ( filename == "" ){
2334 filename = "CaloADC2MIPf.root";
2335 };
2336 stringstream filenome;
2337 filenome.str("");
2338 if ( filevalue == "" ){
2339 const char *pam_calib = pathtocalibration();
2340 //const char *pam_calib=gSystem->Getenv("PAM_CALIB");
2341 filenome.str("");
2342 filenome << pam_calib << "/CaloADC2MIP.root";
2343 // filevalue = Form("%s/CaloADC2MIP.root",pam_calib);
2344 //filevalue = (TString)filenome.str();
2345 } else {
2346 filenome << filevalue.Data();
2347 };
2348 // EM
2349 const char *file2;
2350 // file2 = filevalue;
2351 file2 = filenome.str().c_str();
2352 TFile *hfile = 0;
2353 hfile = new TFile(file2,"UPDATE","Calorimeter CALIBRATION data");
2354 CalorimeterCalibration *calo = 0;
2355 TTree *tree;
2356 tree = (TTree*)hfile->Get("CaloADC");
2357 tree->SetBranchAddress("Event", &calo);
2358 Long64_t nevents = tree->GetEntries();
2359 tree->GetEntry(nevents-1);
2360 // EM end
2361 //
2362 const char *file = filename;
2363 const char *pam_calib = pathtocalibration();
2364 // const char *pam_calib=gSystem->Getenv("PAM_CALIB");
2365 stringstream file3;
2366 // file3 = "/misc/wizard/wizard3/pamela/integr/script/caloScripts-2.00/devel/CaloADC2MIP.rootasdalsdjaksjdlkasjdlkajsdklajsdlkjakljal/";
2367 // sprintf(file3,"%s/%s",pam_calib,file);
2368 file3.str("");
2369 file3 << pam_calib << "/";
2370 file3 << file;
2371 TFile hfile3(file3.str().c_str());
2372 Int_t l = 0;
2373 TCanvas *c1;
2374 Int_t ty = 1;
2375 if ( type != "" ) ty = 0;
2376 c1 = new TCanvas("canvas","canvas",800,800);
2377 c1->Draw();
2378 gStyle->SetOptFit(111);
2379 gPad->SetLogy();
2380 gPad->SetTickx();
2381 gPad->SetTicky();
2382 gPad->SetGrid();
2383 Int_t j = -1;
2384 Int_t k = -1;
2385 Int_t h = -1;;
2386 Int_t changed = 0;
2387 beginning:
2388 if ( j >= 0 ) l = j;
2389 stringstream fmippa;
2390 TF1 *tempf = 0;
2391 TH1F *temp;
2392 while ( l < 2){
2393 Int_t m = 0;
2394 if ( k >= 0 ) m = k;
2395 while ( m < 22 ){
2396 Int_t n =0 ;
2397 if ( h >= 0 ) n = h;
2398 while ( n < 96 ){
2399 c1->Clear();
2400 c1->cd();
2401 fmippa.str("");
2402 fmippa << "fmip " << l;
2403 fmippa << " " << m;
2404 fmippa << " " << n;
2405 // char *name = Form("fmip %i %i %i",l,m,n);
2406 temp = (TH1F*)hfile3.Get(fmippa.str().c_str());
2407 if ( ty ) {
2408 fmippa.str("");
2409 fmippa << "fit " << l;
2410 fmippa << " " << m;
2411 fmippa << " " << n;
2412 // char *name = Form("fit %i %i %i",l,m,n);
2413 tempf = (TF1*)hfile3.Get(fmippa.str().c_str());
2414 };
2415 j = -1;
2416 k = -1;
2417 h = -1;
2418 temp->DrawCopy();
2419 if ( ty ) tempf->DrawCopy("lsame");
2420 c1->Update();
2421 //
2422 //
2423 //
2424 printf(" **** STRIP: %i %i %i **** \n",l,m,n);
2425 printf(" MIP: %f \n",calo->mip[l][m][n]);
2426 printf(" ERMIP: %f \n",calo->ermip[l][m][n]);
2427 printf(" FP[0-3]: %f / %f / %f / %f \n",calo->fp[0][l][m][n],calo->fp[1][l][m][n],calo->fp[2][l][m][n],calo->fp[3][l][m][n]);
2428 printf(" FPE[0-3]: %f / %f / %f / %f \n",calo->fpe[0][l][m][n],calo->fpe[1][l][m][n],calo->fpe[2][l][m][n],calo->fpe[3][l][m][n]);
2429 printf(" CHI2: %f \n",calo->chi2[l][m][n]);
2430 printf(" NDF: %f \n",calo->ndf[l][m][n]);
2431 printf(" MASK: %f \n",calo->mask[l][m][n]);
2432 //
2433 // what to do?
2434 //
2435 char tellme[256];
2436 char tellme2[256];
2437 stringstream input;
2438 stringstream input2;
2439 stringstream out;
2440 stringstream stc;
2441 input.str("a");
2442 out.str("a");
2443 while ( !strcmp(input.str().c_str(),out.str().c_str()) ) {
2444 printf("\nPress <enter> to continue, b<enter> to go backward, j<enter> to jump to a\nfigure, p<enter> to save the figure, f<enter> to do the fit, q<enter> to quit: \n");
2445 cin.getline(tellme,256);
2446 input.str("");
2447 input << tellme;
2448 out.str("1");
2449 //
2450 stc.str("f");
2451 if ( !strcmp(input.str().c_str(),stc.str().c_str()) ) {
2452 printf("\n Do you want to give a number by hand, h<enter>, or do you want to try a better fit, f<enter>?");
2453 cin.getline(tellme,256);
2454 input.str("");
2455 input << tellme;
2456 out.str("f");
2457 if ( !strcmp(input.str().c_str(),stc.str().c_str()) ) {
2458 Double_t fr[2],vtemp;
2459 Double_t sv[4], pllo[4], plhi[4], ffp[4], ffpe[4];
2460 for ( Int_t fk = 0; fk < 4; fk++){
2461 sv[fk] = (int)calo->fp[fk][l][m][n];
2462 printf(" Give input parameter [fp[%i]=%f] \n",fk,sv[fk]);
2463 cin.getline(tellme2,256);
2464 input2.str("");
2465 input2 << tellme2;
2466 vtemp = (Double_t)atoi(input2.str().c_str());
2467 if ( vtemp != 0. ) sv[fk] = vtemp;
2468 printf(" -> fp[%i] = %f \n\n",fk,sv[fk]);
2469 };
2470 printf(" Give lower limit of fitting range fr[0] \n");
2471 cin.getline(tellme2,256);
2472 input2.str("");
2473 input2 << tellme2;
2474 fr[0] = (Double_t)atoi(input2.str().c_str());
2475 printf(" -> fr[0] = %f \n\n",fr[0]);
2476 printf(" Give upper limit of fitting range fr[1] \n");
2477 cin.getline(tellme2,256);
2478 input2.str("");
2479 input2 << tellme2;
2480 fr[1] = (Double_t)atoi(input2.str().c_str());
2481 printf(" -> fr[1] = %f \n",fr[1]);
2482 pllo[0]=0.5; pllo[1]=5.0; pllo[2]=1.0; pllo[3]=0.2;
2483 plhi[0]=8.0; plhi[1]=50.0; plhi[2]=1000000.0; plhi[3]=8.0;
2484 Double_t chisqr;
2485 Int_t fndf;
2486 //
2487 // fit the figure pfmip and put fit results in pfitsnr
2488 //
2489 tempf = langaufit(temp,fr,sv,pllo,plhi,ffp,ffpe,&chisqr,&fndf,"R0");
2490 Double_t SNRPeak = langaupro(ffp);
2491 printf("\n Conversion factor: %f \n",SNRPeak);
2492 tempf->DrawCopy("lsame");
2493 c1->Update();
2494 c1->Modified();
2495 c1->Update();
2496 printf(" Do you want to save this result y/n<enter> ?\n");
2497 cin.getline(tellme,256);
2498 input.str("");
2499 input << tellme;
2500 out.str("y");
2501 if ( !strcmp(input.str().c_str(),stc.str().c_str()) ) {
2502 calo->mip[l][m][n] = (float)SNRPeak;
2503 calo->ermip[l][m][n] = (float)ffpe[1];
2504 for (Int_t a = 0; a < 4 ; a++){
2505 calo->fp[a][l][m][n] = (float)ffp[a];
2506 calo->fpe[a][l][m][n] = (float)ffpe[a];
2507 };
2508 calo->ndf[l][m][n] = (float)fndf;
2509 calo->chi2[l][m][n] = (float)chisqr;
2510 //
2511 if ( fndf!=0 && (float)chisqr/(float)fndf < 2. && ffpe[1] < 0.15 && ffp[1] > 15. && ffp[1] < 40. ) {
2512 calo->mask[l][m][n] = 1.;
2513 } else {
2514 calo->mask[l][m][n] = 0.;
2515 };
2516 printf(" Saved! \n");
2517 changed = 1;
2518 };
2519 } else {
2520 printf(" Give the MIP conversion value \n");
2521 cin.getline(tellme2,256);
2522 input2.str("");
2523 input2 << tellme2;
2524 calo->mip[l][m][n] = (Float_t)atoi(input2.str().c_str());
2525 printf(" -> mip = %f \n",calo->mip[l][m][n]);
2526 printf(" Give the mask value \n");
2527 cin.getline(tellme2,256);
2528 input2.str("");
2529 input2 << tellme2;
2530 calo->mask[l][m][n] = (Float_t)atoi(input2.str().c_str());
2531 printf(" -> mask = %f \n",calo->mask[l][m][n]);
2532 changed = 1;
2533 };
2534 out.str("");
2535 out << input.str().c_str();
2536 };
2537 stc.str("b");
2538 if ( !strcmp(input.str().c_str(),stc.str().c_str()) ) {
2539 if ( n > 0 ) {
2540 printf("WARNING: going one figure backward!\n\n");
2541 n -= 2;
2542 } else {
2543 printf("This is the first strip of plane, you can't go backward! \n");
2544 out.str("");
2545 out << input.str().c_str();
2546 };
2547 };
2548 stc.str("j");
2549 if ( !strcmp(input.str().c_str(),stc.str().c_str()) ) {
2550 printf("\n Enter the view number you want to jump to (0 = x, 1 = y): ");
2551 cin.getline(tellme2,256);
2552 input2.str("");
2553 input2 << tellme2;
2554 j = atoi(input2.str().c_str());
2555 if ( j < 0 || j > 1 ) {
2556 printf("\n You can choose between 0 and 1 \n");
2557 out.str("");
2558 out << input.str().c_str();
2559 } else {
2560 printf("\n Enter the plane number you want to jump to (0 to 21): ");
2561 cin.getline(tellme2,256);
2562 input2.str("");
2563 input2 << tellme2;
2564 k = atoi(input2.str().c_str());
2565 if ( k < 0 || k > 21 ) {
2566 printf("\n You can choose between 0 and 21 \n");
2567 out.str("");
2568 out << input.str().c_str();
2569 } else {
2570 printf("\n Enter the strip number you want to jump to (0 to 95): ");
2571 cin.getline(tellme2,256);
2572 input2.str("");
2573 input2 << tellme2;
2574 h = atoi(input2.str().c_str());
2575 if ( h < 0 || h > 95 ) {
2576 printf("\n You can choose between 0 and 95 \n");
2577 out.str("");
2578 out << input.str().c_str();
2579 } else {
2580 printf("\n Jumping to strip %i %i %i \n\n",j,k,h);
2581 goto beginning;
2582 };
2583 };
2584 };
2585 };
2586 //
2587 stc.str("q");
2588 if ( !strcmp(input.str().c_str(),stc.str().c_str()) ) {
2589 printf("Exiting...\n");
2590 goto end;
2591 };
2592 //
2593 stc.str("p");
2594 if ( !strcmp(input.str().c_str(),stc.str().c_str()) ) {
2595 printf("Enter a file extension recognized by ROOT (ps, eps, gif,...):\n");
2596 cin.getline(tellme2,256);
2597 input2.str("");
2598 input2 << tellme2;
2599 out.str("");
2600 out << input.str().c_str();
2601 stringstream figsave;
2602 figsave.str("");
2603 figsave << startingdir << "/mip2adc_";
2604 figsave << l << "_";
2605 figsave << m << "_";
2606 figsave << n << ".";
2607 figsave << input2.str().c_str();
2608 // char *figsave = 0;
2609 //figsave = Form("%s/mip2adc_%i_%i_%i.%s",startingdir,l,m,n,input2);
2610 c1->SaveAs(figsave.str().c_str());
2611 printf("\n");
2612 };
2613 };
2614 n++;
2615 };
2616 m++;
2617 };
2618 l++;
2619 };
2620 end: printf("\n");
2621 if ( changed ){
2622 tree->Fill();
2623 hfile->Write();
2624 };
2625 hfile->Close();
2626 }
2627
2628 void Calo4224MIPVALUES(TString filevalue = "", TString type=""){
2629 //
2630 // this routine shows the 4224 figures with fit from the main program
2631 //
2632 gROOT->GetListOfCanvases()->Delete("a");
2633 gROOT->Reset("a");
2634 gDirectory->GetList()->Delete();
2635 const char* startingdir = gSystem->WorkingDirectory();
2636 #if defined (__CINT__)
2637 emicheckLib();
2638 stringstream llib;
2639 const char *pamlib = gSystem->ExpandPathName("$PAM_LIB");
2640 llib << pamlib << "/caloclasses_h.so";
2641 gSystem->Load(llib.str().c_str());
2642 stringstream lolib;
2643 lolib.str("");
2644 lolib << pamlib << "/UTpathtoc_c.so";
2645 gSystem->Load(lolib.str().c_str());
2646 #endif
2647 stringstream filenome;
2648 filenome.str("");
2649 if ( filevalue == "" ){
2650 const char *pam_calib = pathtocalibration();
2651 //const char *pam_calib=gSystem->Getenv("PAM_CALIB");
2652 filenome << pam_calib << "/CaloADC2MIP.root";
2653 // filevalue = Form("%s/CaloADC2MIP.root",pam_calib);
2654 //filevalue = (TString)filenome.str().c_str();
2655 } else {
2656 filenome << filevalue.Data();
2657 };
2658 // EM
2659 const char *file2;
2660 file2 = filenome.str().c_str();//filevalue;
2661 TFile *hfile = 0;
2662 hfile = new TFile(file2,"READ","Calorimeter CALIBRATION data");
2663 CalorimeterCalibration *calo = 0;
2664 TTree *tree;
2665 tree = (TTree*)hfile->Get("CaloADC");
2666 tree->SetBranchAddress("Event", &calo);
2667 Long64_t nevents = tree->GetEntries();
2668 tree->GetEntry(nevents-1);
2669 // EM end
2670 //
2671 TCanvas *c1;
2672 Int_t ty = 1;
2673 if ( type != "" ) ty = 0;
2674 c1 = new TCanvas("canvas","canvas",1200,800);
2675 c1->Draw();
2676 gStyle->SetOptFit(111);
2677 // gPad->SetLiny();
2678 gPad->SetTickx();
2679 gPad->SetTicky();
2680 // gPad->SetGrid();
2681 Int_t k = -1;
2682 Int_t changed = 0;
2683 Float_t mipf[1]={0};
2684 Float_t val[1]={0};
2685 TLatex *text=new TLatex();
2686 TLine *linea;
2687 stringstream testo;
2688 // char *testo=" looong long line as workaround since root is not working other ways dsf sdf sdd fsd fsd sdfsfdsdkasasd dasdkajsdhkajs da hdkjashdkjashdkjashdkjashdkjashdkasj";
2689 beginning:
2690 Int_t m = 0;
2691 if ( k >= 0 ) m = k;
2692 while ( m < 22 ){
2693 Int_t l = 0;
2694 c1->Clear();
2695 c1->cd();
2696 TMultiGraph *mg = new TMultiGraph();
2697 while ( l < 2){
2698 Int_t n = 0;
2699 while ( n < 96 ){
2700 val[0] = (float)n + 100.*(float)l ;
2701 mipf[0] = (float)calo->mip[l][m][n];// + 100.*(float)l;
2702 TGraph *graph = new TGraph(1, val, mipf);
2703 graph->SetMarkerColor(2);
2704 graph->SetMarkerStyle(22+l);
2705 graph->SetMarkerSize(0.7);
2706 mg->Add(graph);
2707 printf(" **** STRIP: %i %i %i **** \n",l,m,n);
2708 printf(" MIP: %f \n",calo->mip[l][m][n]);
2709 printf(" ERMIP: %f \n",calo->ermip[l][m][n]);
2710 printf(" FP[0-3]: %f / %f / %f / %f \n",calo->fp[0][l][m][n],calo->fp[1][l][m][n],calo->fp[2][l][m][n],calo->fp[3][l][m][n]);
2711 printf(" FPE[0-3]: %f / %f / %f / %f \n",calo->fpe[0][l][m][n],calo->fpe[1][l][m][n],calo->fpe[2][l][m][n],calo->fpe[3][l][m][n]);
2712 printf(" CHI2: %f \n",calo->chi2[l][m][n]);
2713 printf(" NDF: %f \n",calo->ndf[l][m][n]);
2714 printf(" MASK: %f \n",calo->mask[l][m][n]);
2715 n++;
2716 };
2717 l++;
2718 };
2719 mg->SetMaximum(50.);
2720 mg->SetMinimum(0.);
2721 mg->Draw("AP");
2722 for (Int_t e = 0; e<2 ; e++){
2723 for (Int_t ep = 0; ep<96 ; ep++){
2724 if ( ep%16 == 0 ) {
2725 linea = new TLine((float)ep+100.*(float)e,0.,(float)ep+100.*(float)e,50.);
2726 linea->SetLineColor(15);
2727 linea->SetLineStyle(2);
2728 linea->SetLineWidth(1);
2729 linea->Draw("lsame");
2730 if ( ep > 0 ){
2731 linea = new TLine((float)ep-1.+100.*(float)e,0.,(float)ep-1.+100.*(float)e,50.);
2732 linea->SetLineColor(15);
2733 linea->SetLineStyle(2);
2734 linea->SetLineWidth(1);
2735 linea->Draw("lsame");
2736 };
2737 };
2738 if ( ep == 95 ){
2739 linea = new TLine(95.+100.*(float)e,0.,95.+100.*(float)e,50.);
2740 linea->SetLineColor(15);
2741 linea->SetLineStyle(2);
2742 linea->SetLineWidth(1);
2743 linea->Draw("lsame");
2744 };
2745 };
2746 };
2747 text->SetTextAngle(0);
2748 text->SetTextFont(32);
2749 text->SetTextColor(1);
2750 text->SetTextSize(0.050); // 0.02
2751 text->SetTextAlign(12);
2752 testo.str("");
2753 testo << "Plane " << m;
2754 // sprintf(testo,"Plane %i",m);
2755 text->DrawLatex(85.,52.,testo.str().c_str());
2756 c1->Update();
2757 //
2758 // what to do?
2759 //
2760
2761 char tellme[256];
2762 char tellme2[256];
2763 stringstream input;
2764 stringstream input2;
2765 stringstream out;
2766 stringstream stc;
2767 input.str("a");
2768 out.str("a");
2769 while ( !strcmp(input.str().c_str(),out.str().c_str()) ) {
2770 printf("\nPress <enter> to continue, b<enter> to go backward, j<enter> to jump to a\nfigure, p<enter> to save the figure, q<enter> to quit: \n");
2771 cin.getline(tellme,256);
2772 input.str("");
2773 input << tellme;
2774 out.str("1");
2775 //
2776 stc.str("b");
2777 if ( !strcmp(input.str().c_str(),stc.str().c_str()) ) {
2778 if ( m > 0 ) {
2779 printf("WARNING: going one figure backward!\n\n");
2780 m -= 2;
2781 } else {
2782 printf("This is the first plane, you can't go backward! \n");
2783 out.str("");
2784 out << input.str().c_str();
2785 };
2786 };
2787 stc.str("j");
2788 if ( !strcmp(input.str().c_str(),stc.str().c_str()) ) {
2789 printf("\n Enter the plane number you want to jump to (0 to 21): ");
2790 cin.getline(tellme2,256);
2791 input2.str("");
2792 input2 << tellme2;
2793 k = atoi(input2.str().c_str());
2794 if ( k < 0 || k > 21 ) {
2795 printf("\n You can choose between 0 and 21 \n");
2796 out.str("");
2797 out << input.str().c_str();
2798 } else {
2799 printf("\n Jumping to plane %i \n\n",k);
2800 m = k;
2801 goto beginning;
2802 };
2803 };
2804 //
2805 stc.str("q");
2806 if ( !strcmp(input.str().c_str(),stc.str().c_str()) ) {
2807 printf("Exiting...\n");
2808 goto end;
2809 };
2810 //
2811 stc.str("p");
2812 if ( !strcmp(input.str().c_str(),stc.str().c_str()) ) {
2813 printf("Enter a file extension recognized by ROOT (ps, eps, gif,...):\n");
2814 cin.getline(tellme2,256);
2815 input2.str("");
2816 input2 << tellme2;
2817 out.str("");
2818 out << input.str().c_str();
2819 stringstream figsave;
2820 figsave.str("");
2821 figsave << startingdir << "/mip2adc_plane_";
2822 figsave << m << ".";
2823 figsave << input2.str().c_str();
2824 // char *figsave = 0;
2825 //figsave = Form("%s/mip2adc_%i_%i_%i.%s",startingdir,l,m,n,input2);
2826 c1->SaveAs(figsave.str().c_str());
2827 printf("\n");
2828 };
2829 };
2830 m++;
2831 };
2832 end: printf("\n");
2833 if ( changed ){
2834 tree->Fill();
2835 hfile->Write();
2836 };
2837 hfile->Close();
2838 }
2839
2840 void Calo4224BAK( TString filename ){
2841 //
2842 // this shows the 4224 figure from a backup file
2843 //
2844 gROOT->GetListOfCanvases()->Delete("a");
2845 gROOT->Reset("a");
2846 gDirectory->GetList()->Delete();
2847 const char* startingdir = gSystem->WorkingDirectory();
2848 #if defined (__CINT__)
2849 emicheckLib();
2850 stringstream llib;
2851 const char *pamlib = gSystem->ExpandPathName("$PAM_LIB");
2852 llib << pamlib << "/caloclasses_h.so";
2853 gSystem->Load(llib.str().c_str());
2854 stringstream lolib;
2855 lolib.str("");
2856 lolib << pamlib << "/UTpathtoc_c.so";
2857 gSystem->Load(lolib.str().c_str());
2858 #endif
2859 const char *file3 = filename;
2860 TFile hfile3(file3);
2861 Int_t l = 0;
2862 TCanvas *c1;
2863 c1 = new TCanvas("canvas","canvas",800,800);
2864 c1->Draw();
2865 gStyle->SetOptFit(111);
2866 gPad->SetLogy();
2867 Int_t j = 0;
2868 Int_t k = 0;
2869 Int_t h = 0;
2870 beginning:
2871 if ( j ) l = j;
2872 stringstream bmippa;
2873 while ( l < 2){
2874 Int_t m = 0;
2875 if ( k ) m = k;
2876 while ( m < 22 ){
2877 Int_t n =0 ;
2878 if ( h ) n = h;
2879 while ( n < 96 ){
2880 c1->Clear();
2881 c1->cd();
2882 bmippa.str("");
2883 bmippa << "bmip " << l;
2884 bmippa << " " << m;
2885 bmippa << " " << n;
2886 // char *name = Form("bmip %i %i %i",l,m,n);
2887 TH1F *temp = (TH1F*)hfile3.Get(bmippa.str().c_str());
2888 j = 0;
2889 k = 0;
2890 h = 0;
2891 temp->DrawCopy();
2892 c1->Update();
2893 //
2894 // what to do?
2895 //
2896 char tellme[256];
2897 char tellme2[256];
2898 stringstream input;
2899 stringstream input2;
2900 stringstream out;
2901 stringstream stc;
2902 input.str("a");
2903 out.str("a");
2904 while ( !strcmp(input.str().c_str(),out.str().c_str()) ) {
2905 printf("\nPress <enter> to continue, b<enter> to go backward, j<enter> to jump to a\nfigure, p<enter> to save the figure, q<enter> to quit: \n");
2906 cin.getline(tellme,256);
2907 input.str("");
2908 input << tellme;
2909 out.str("1");
2910 //
2911 stc.str("b");
2912 if ( !strcmp(input.str().c_str(),stc.str().c_str()) ) {
2913 if ( n > 0 ) {
2914 printf("WARNING: going one figure backward!\n\n");
2915 n -= 2;
2916 } else {
2917 printf("This is the first strip of plane, you can't go backward! \n");
2918 out.str("");
2919 out << input.str().c_str();
2920 };
2921 };
2922 stc.str("j");
2923 if ( !strcmp(input.str().c_str(),stc.str().c_str()) ) {
2924 printf("\n Enter the view number you want to jump to (0 = x, 1 = y): ");
2925 cin.getline(tellme2,256);
2926 input2.str("");
2927 input2 << tellme2;
2928 j = atoi(input2.str().c_str());
2929 if ( j < 0 || j > 1 ) {
2930 printf("\n You can choose between 0 and 1 \n");
2931 out.str("");
2932 out << input.str().c_str();
2933 } else {
2934 printf("\n Enter the plane number you want to jump to (0 to 21): ");
2935 cin.getline(tellme2,256);
2936 input2.str("");
2937 input2 << tellme2;
2938 k = atoi(input2.str().c_str());
2939 if ( k < 0 || k > 21 ) {
2940 printf("\n You can choose between 0 and 21 \n");
2941 out.str("");
2942 out << input.str().c_str();
2943 } else {
2944 printf("\n Enter the strip number you want to jump to (0 to 95): ");
2945 cin.getline(tellme2,256);
2946 input2.str("");
2947 input2 << tellme2;
2948 h = atoi(input2.str().c_str());
2949 if ( h < 0 || h > 95 ) {
2950 printf("\n You can choose between 0 and 95 \n");
2951 out.str("");
2952 out << input.str().c_str();
2953 } else {
2954 printf("\n Jumping to strip %i %i %i \n\n",j,k,h);
2955 goto beginning;
2956 };
2957 };
2958 };
2959 };
2960 //
2961 stc.str("q");
2962 if ( !strcmp(input.str().c_str(),stc.str().c_str()) ) {
2963 printf("Exiting...\n");
2964 goto end;
2965 };
2966 //
2967 stc.str("p");
2968 if ( !strcmp(input.str().c_str(),stc.str().c_str()) ) {
2969 printf("Enter a file extension recognized by ROOT (ps, eps, gif,...):\n");
2970 cin.getline(tellme2,256);
2971 input2.str("");
2972 input2 << tellme2;
2973 out.str("");
2974 out << input.str().c_str();
2975 stringstream figsave;
2976 figsave.str("");
2977 figsave << startingdir << "/mip2adcbak_";
2978 figsave << l << "_";
2979 figsave << m << "_";
2980 figsave << n << ".";
2981 figsave << input2.str().c_str();
2982 // char *figsave = 0;
2983 //figsave = Form("%s/mip2adc_%i_%i_%i.%s",startingdir,l,m,n,input2);
2984 c1->SaveAs(figsave.str().c_str());
2985 printf("\n");
2986 };
2987 };
2988 n++;
2989 };
2990 m++;
2991 };
2992 l++;
2993 };
2994 end: printf("\n");
2995 }
2996
2997 void CaloBAKFIT( TString filename ){
2998 //
2999 // this routine perform the 4224 fit on a backup figure file and save - with extension .bak - the two output file
3000 // in the correct format, ready to be updated etc. etc.
3001 // NB: you will LOSE information about used files
3002 //
3003 gROOT->GetListOfCanvases()->Delete("a");
3004 gROOT->Reset("a");
3005 gDirectory->GetList()->Delete();
3006 #if defined (__CINT__)
3007 emicheckLib();
3008 stringstream llib;
3009 const char *pamlib = gSystem->ExpandPathName("$PAM_LIB");
3010 llib << pamlib << "/caloclasses_h.so";
3011 gSystem->Load(llib.str().c_str());
3012 stringstream lolib;
3013 lolib.str("");
3014 lolib << pamlib << "/UTpathtoc_c.so";
3015 gSystem->Load(lolib.str().c_str());
3016 #endif
3017 const char *pam_calib = pathtocalibration();
3018 // const char *pam_calib=gSystem->Getenv("PAM_CALIB");
3019 stringstream file2;
3020 //char *file2;
3021 // file2 = "/misc/wizard/wizard3/pamela/integr/script/caloScripts-2.00/devel/CaloADC2MIP.rootasdalsdjaksjdlkasjdlkajsdklajsdlkjakljal/";
3022 file2.str("");
3023 file2 << pam_calib << "/CaloADC2MIP.bak";
3024 // sprintf(file2,"%s/CaloADC2MIP.bak",pam_calib);
3025 stringstream file4;
3026 file4.str("");
3027 // char *file4;
3028 //file4 = "/misc/wizard/wizard3/pamela/integr/script/caloScripts-2.00/devel/CaloADC2MIP.rootasdalsdjaksjdlkasjdlkajsdklajsdlkjakljal/";
3029 // sprintf(file4,"%s/CaloADC2MIPf.bak",pam_calib);
3030 file4 << pam_calib << "/CaloADC2MIPf.bak";
3031 //
3032 Int_t notgood = 0;
3033 //
3034 TFile *hfile = 0;
3035 CalorimeterCalibration *calo = new CalorimeterCalibration();
3036 hfile = new TFile(file2.str().c_str(),"RECREATE","Calorimeter CALIBRATION data");
3037 if ( hfile->IsZombie() ){
3038 printf(" ERROR: Cannot create file \n");
3039 return;
3040 };
3041 TTree *tree = 0;
3042 tree = new TTree("CaloADC","Calorimeter calibration data");
3043 tree->Branch("Event","CalorimeterCalibration",&calo);
3044 //
3045 const char *file3 = filename;
3046 TFile hfile3(file3);
3047 TCanvas *c1;
3048 c1 = new TCanvas("canvas","canvas",800,800);
3049 c1->Draw();
3050 gStyle->SetOptFit(111);
3051 TH1F *pfmip[2][22][96];
3052 TF1 *pfitsnr[2][22][96];
3053 stringstream Fmip;
3054 for (Int_t i=0; i<2;i++){
3055 for (Int_t j=0; j<22;j++){
3056 for (Int_t m=0; m<96;m++){
3057 Fmip.str("");
3058 Fmip << "Fmip " << i;
3059 Fmip << " " << j;
3060 Fmip << " " << m;
3061 pfmip[i][j][m] = new TH1F(Fmip.str().c_str(),"",56,0.,55.);
3062 pfitsnr[i][j][m] = new TF1;
3063 };
3064 };
3065 };
3066 gStyle->SetOptFit(111);
3067 gPad->SetLogy();
3068 //
3069 Int_t l = 0;
3070 stringstream bmippa;
3071 while ( l < 2){
3072 Int_t m = 0;
3073 while ( m < 22 ){
3074 Int_t n =0 ;
3075 Double_t SNRPeak = 0.;
3076 while ( n < 96 ){
3077 c1->Clear();
3078 c1->cd();
3079 bmippa.str("");
3080 bmippa << "bmip " << l;
3081 bmippa << " " << m;
3082 bmippa << " " << n;
3083 // char *name = Form("bmip %i %i %i",l,m,n);
3084 TH1F *temp = (TH1F*)hfile3.Get(bmippa.str().c_str());
3085 bmippa.str("");
3086 bmippa << "fmip " << l;
3087 bmippa << " " << m;
3088 bmippa << " " << n;
3089 pfmip[l][m][n] = (TH1F*)temp->Clone(bmippa.str().c_str());
3090 Int_t doitagain = 0;
3091 Double_t fr[2];
3092 Double_t sv[4], pllo[4], plhi[4], fp[4], fpe[4];
3093 if ( SNRPeak > 16. && SNRPeak < 35. ){
3094 fr[0] = (Float_t)SNRPeak - 7.;
3095 sv[0]= fp[0];
3096 sv[1]= fp[1];
3097 sv[2]= fp[2];
3098 sv[3]= fp[3];
3099 } else {
3100 fr[0] = 19.;
3101 sv[0] = 2.8;
3102 sv[1] = 21.0;
3103 sv[2] = 1000.0;
3104 sv[3] = 0.1;
3105 };
3106 fitting: printf("Fitting strip %i %i %i \n",l,m,n);
3107 fr[1] = 60.;
3108 pllo[0]=0.5; pllo[1]=5.0; pllo[2]=1.0; pllo[3]=0.2;
3109 plhi[0]=8.0; plhi[1]=50.0; plhi[2]=1000000.0; plhi[3]=8.0;
3110 Double_t chisqr;
3111 Int_t ndf;
3112 pfitsnr[l][m][n] = langaufit(pfmip[l][m][n],fr,sv,pllo,plhi,fp,fpe,&chisqr,&ndf);
3113 //
3114 Double_t SNRPeak = langaupro(fp);
3115 printf("\n Conversion factor: %f \n\n",SNRPeak);
3116 //
3117 //
3118 //
3119 if ( (SNRPeak<0. || (SNRPeak < 16. || SNRPeak > 35.) || chisqr > 100.) && doitagain < 3 ){
3120 printf("\n The fit is not good, I will try again, step %i \n\n",doitagain);
3121 doitagain++;
3122 if ( chisqr > 100. ) {
3123 fr[0] = fr[0] + 5.;
3124 sv[0] = fp[0];
3125 sv[1] = fp[1];
3126 sv[2] = fp[2];
3127 sv[3] = fp[3];
3128 goto fitting;
3129 };
3130 if ( doitagain == 1 ) {
3131 fr[0] = 19.;
3132 sv[0] = 2.8;
3133 sv[1] = 21.0;
3134 sv[2] = 1000.0;
3135 sv[3] = 0.1;
3136 goto fitting;
3137 };
3138 if ( doitagain == 2 ) {
3139 fr[0] = 22.;
3140 sv[0] = 2.8;
3141 sv[1] = 25.0;
3142 sv[2] = 1000.0;
3143 sv[3] = 0.1;
3144 goto fitting;
3145 };
3146 if ( doitagain == 3 ) {
3147 fr[0] = 12.;
3148 sv[0] = 2.8;
3149 sv[1] = 15.0;
3150 sv[2] = 1000.0;
3151 sv[3] = 0.1;
3152 goto fitting;
3153 };
3154 };
3155 //
3156 calo->mip[l][m][n] = (float)SNRPeak;
3157 calo->ermip[l][m][n] = (float)fpe[1];
3158 for (Int_t a = 0; a < 4 ; a++){
3159 calo->fp[a][l][m][n] = (float)fp[a];
3160 calo->fpe[a][l][m][n] = (float)fpe[a];
3161 };
3162 calo->ndf[l][m][n] = (float)ndf;
3163 calo->chi2[l][m][n] = (float)chisqr;
3164 //
3165 if ( ndf!=0 && chisqr/ndf < 2. && fpe[1] < 0.15 && fp[1] > 15. && fp[1] < 40. && SNRPeak > 0. ) {
3166 printf("\n THIS IS A GOOD FIT, SAVED! \n\n");
3167 calo->mask[l][m][n] = 1.;
3168 } else {
3169 notgood++;
3170 calo->mask[l][m][n] = 0.;
3171 };
3172 printf("Fitting done, strip %i %i %i \n\n\n",l,m,n);
3173 c1->cd();
3174 pfmip[l][m][n]->DrawCopy();
3175 pfitsnr[l][m][n]->DrawCopy("lsame");
3176 c1->Modified();
3177 c1->Update();
3178 n++;
3179 };
3180 m++;
3181 };
3182 l++;
3183 };
3184 //
3185 calo->status = 0;
3186 if ( !notgood ) calo->status = 1;
3187 //
3188 tree->Fill();
3189 hfile->Write();
3190 hfile->Close();
3191 //
3192 // Save histograms
3193 //
3194 TFile *hfile2;
3195 hfile2 = new TFile(file4.str().c_str(),"RECREATE","Calorimeter CALIBRATION figures");
3196 if ( hfile2->IsZombie() ){
3197 printf(" ERROR: Cannot create file \n");
3198 return;
3199 };
3200 TH1F *sfmip[2][22][96];
3201 TF1 *sfitsnr[2][22][96];
3202 stringstream fmipfi;
3203 for (Int_t l = 0; l < 2; l++){
3204 for (Int_t m = 0; m < 22; m++){
3205 for (Int_t n =0 ; n < 96; n++){
3206 fmipfi.str("");
3207 fmipfi << "fmip " << l;
3208 fmipfi << " " << m;
3209 fmipfi << " " << n;
3210 sfmip[l][m][n] = new TH1F(fmipfi.str().c_str(),"",56,0.,55.);
3211 sfmip[l][m][n] = (TH1F*)pfmip[l][m][n]->Clone(fmipfi.str().c_str());
3212 fmipfi.str("");
3213 fmipfi << "fit " << l;
3214 fmipfi << " " << m;
3215 fmipfi << " " << n;
3216 sfitsnr[l][m][n] = (TF1*)pfitsnr[l][m][n]->Clone(fmipfi.str().c_str());
3217 sfmip[l][m][n]->Write();
3218 sfitsnr[l][m][n]->Write();
3219 };
3220 };
3221 };
3222 hfile2->Close();
3223 printf("\n");
3224 hfile3.Close();
3225 Calo4224STATUS("CaloADC2MIP.bak");
3226 }
3227
3228 void CaloLOOKATSTRIP(Int_t view, Int_t plane, Int_t strip, Int_t fromevno, Int_t toevno,Int_t fromtime=0, Int_t totime=1000000000, Int_t fit = 0){
3229 emicheckLib();
3230 #if defined (__CINT__)
3231 stringstream llib;
3232 const char *pamlib = gSystem->ExpandPathName("$PAM_LIB");
3233 llib << pamlib << "/caloclasses_h.so";
3234 gSystem->Load(llib.str().c_str());
3235 stringstream lolib;
3236 lolib.str("");
3237 lolib << pamlib << "/UTpathtoc_c.so";
3238 gSystem->Load(lolib.str().c_str());
3239 const char *pam_calib = pathtocalibration();
3240 // const char *pam_calib=gSystem->Getenv("PAM_CALIB");
3241 stringstream fname;
3242 fname << pam_calib << "/CaloADC2MIPdata.raw";
3243 // char *fname=Form("%s/CaloADC2MIPdata.raw",pam_calib);
3244 TFile *f = new TFile(fname.str().c_str());
3245 TTree *tr = (TTree*) f->Get("CaloADCRAW");
3246 // printf(" input1 = %s \n input2 = %s \n",input1,input2);
3247 TCanvas *myc = new TCanvas("myc","myc",900,800);
3248 gStyle->SetOptFit(111);
3249 gStyle->SetOptStat(10000001);
3250 gStyle->SetStatH(0.5);
3251 gStyle->SetTitleW(0.6);
3252 stringstream input1;
3253 stringstream input2;
3254 if ( fromtime == 0 ){
3255 input1.str("");
3256 input1 << "estrip[" << view;
3257 input1 << "][" << plane;
3258 input1 << "][" << strip;
3259 input1 << "]:evno";
3260 // input1 = Form("estrip[%i][%i][%i]:evno",view,plane,strip);
3261 input2.str("");
3262 input2 << "evno>" << fromevno;
3263 input2 << " && evno<" << toevno;
3264 input2 << " && estrip[" << view;
3265 input2 << "][" << plane ;
3266 input2 << "][" << strip;
3267 input2 << "]>-50.";
3268 // input2 = Form("evno>%i && evno<%i && estrip[%i][%i][%i]>-50.",fromevno,toevno,view,plane,strip);
3269 tr->Draw(input1.str().c_str(),input2.str().c_str(),"l");
3270 } else {
3271 input1.str("");
3272 input1 << "estrip[" << view;
3273 input1 << "][" << plane;
3274 input1 << "][" << strip;
3275 input1 << "]:evno>>ha2";
3276 // input1 = Form("estrip[%i][%i][%i]:etime>>ha2",view,plane,strip);
3277 input2.str("");
3278 input2 << "etime>" << fromtime;
3279 input2 << " && etime<" << totime;
3280 input2 << " && evno>" << fromevno;
3281 input2 << " && evno<" << toevno;
3282 input2 << " && estrip[" << view;
3283 input2 << "][" << plane ;
3284 input2 << "][" << strip;
3285 input2 << "]>-50.";
3286 // input2 = Form("etime>%i && etime<%i && evno>%i && evno<%i && estrip[%i][%i][%i]>-50.",fromtime,totime,fromevno,toevno,view,plane,strip);
3287 tr->Draw(input1.str().c_str(),input2.str().c_str(),"prof");
3288 };
3289 if ( !fit ) return;
3290 Double_t par[4];
3291 TF1 *func = new TF1("fitraw",fitraw,fromtime,totime,4);
3292 func->SetParameters(15.,fromtime,0.00001,-7.);
3293 func->SetParNames("p0","p1","p2","p3");
3294 ha2->Fit("fitraw","r");
3295 stringstream figname;
3296 figname.str("");
3297 figname << "strip" << view;
3298 figname << plane;
3299 figname << strip;
3300 figname << "t" << fromtime;
3301 figname << ".ps";
3302 // char *figname = 0;
3303 //figname = Form("strip%i%i%it%i.ps",view,plane,strip,fromtime);
3304 myc->SaveAs(figname.str().c_str());
3305 #endif
3306 }
  ViewVC Help
Powered by ViewVC 1.1.23