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

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 1.1 - (hide annotations) (download)
Tue Dec 6 10:59:28 2005 UTC (19 years, 10 months ago) by mocchiut
Branch point for: MAIN, UTILITIES
File MIME type: text/plain 
Initial revision
1 mocchiut 1.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