/[PAMELA software]/DarthVader/CalorimeterLevel2/src/CaloLevel0.cpp
ViewVC logotype

Diff of /DarthVader/CalorimeterLevel2/src/CaloLevel0.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.1 by mocchiut, Sun Mar 18 20:22:54 2007 UTC revision 1.23 by mocchiut, Thu Dec 4 13:50:42 2008 UTC
# Line 64  CaloLevel0::CaloLevel0(){ Line 64  CaloLevel0::CaloLevel0(){
64    memset(mip, 0, 2*22*96*sizeof(Float_t));    memset(mip, 0, 2*22*96*sizeof(Float_t));
65    memset(base, 0, 2*22*6*sizeof(Float_t));    memset(base, 0, 2*22*6*sizeof(Float_t));
66    memset(sbase, 0, 2*22*6*sizeof(Float_t));    memset(sbase, 0, 2*22*6*sizeof(Float_t));
67      memset(obadmask, 0, 2*22*96*sizeof(Int_t));
68      memset(obadpulsemask, 0, 2*22*6*sizeof(Int_t));
69      memset(ctprecor, 0, 2*22*6*sizeof(Float_t));
70      memset(ctsicor, 0, 2*22*9*sizeof(Float_t));
71      memset(ctneigcor, 0, 2*22*6*sizeof(Float_t));
72    calopar1 = true;    calopar1 = true;
73    calopar2 = true;    calopar2 = true;
74    calopar3 = true;    calopar3 = true;
75      calopar4 = true;
76      calopar5 = true;
77    crosst = true;    crosst = true;
78    ftcalopar1 = 0;    ftcalopar1 = 0;
79    ttcalopar1 = 0;    ttcalopar1 = 0;
# Line 74  CaloLevel0::CaloLevel0(){ Line 81  CaloLevel0::CaloLevel0(){
81    ttcalopar2 = 0;    ttcalopar2 = 0;
82    ftcalopar3 = 0;    ftcalopar3 = 0;
83    ttcalopar3 = 0;    ttcalopar3 = 0;
84      ftcalopar4 = 0;
85      ttcalopar4 = 0;
86      ftcalopar5 = 0;
87      ttcalopar5 = 0;
88  }  }
89    
90  void CaloLevel0::SetCrossTalk(Bool_t ct){  void CaloLevel0::SetCrossTalk(Bool_t ct){
91    crosst = ct;    crosst = ct;
92  };  }
93    
94    void CaloLevel0::SetCrossTalkType(Bool_t ct){
95      ctground = ct;
96    }
97    
98    void CaloLevel0::SetCrossTalkType(Int_t ct){
99      if ( ct == 0 ) ctground = true;
100      if ( ct == 1 ){
101        ctground = false;
102        noselfct = false;
103      };
104      if ( ct == 2 ){
105        ctground = false;
106        noselfct = true;
107      };
108    }
109    
110  void CaloLevel0::SetVerbose(Bool_t ct){  void CaloLevel0::SetVerbose(Bool_t ct){
111    verbose = ct;    verbose = ct;
112  };  }
113    
114  /**  /**
115   * Initialize CaloLevel0 object   * Initialize CaloLevel0 object
116  **/  **/
117  void CaloLevel0::ProcessingInit(TSQLServer *dbc, UInt_t hs, Int_t &sgnl, TTree *l0tree, Bool_t isdeb, Bool_t isverb){  void CaloLevel0::ProcessingInit(TSQLServer *dbc, UInt_t hs, Int_t &sgnl, TTree *l0tree, Bool_t isdeb, Bool_t isverb){
118      if ( !dbc->IsConnected() ) throw -116;  
119      this->InitDo(dbc,hs,sgnl,l0tree,isdeb,isverb);
120    }
121    
122    /**
123     * Initialize CaloLevel0 object
124    **/
125    void CaloLevel0::ProcessingInit(GL_TABLES *glt, UInt_t hs, Int_t &sgnl, TTree *l0tree, Bool_t isdeb, Bool_t isverb){
126      //
127      const TString host = glt->CGetHost();
128      const TString user = glt->CGetUser();
129      const TString psw = glt->CGetPsw();
130      TSQLServer *dbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data());
131      if ( !dbc->IsConnected() ) throw -116;  
132      this->InitDo(dbc,hs,sgnl,l0tree,isdeb,isverb);
133      dbc->Close();
134      delete dbc;
135    }
136    
137    
138    void CaloLevel0::InitDo(TSQLServer *dbc, UInt_t hs, Int_t &sgnl, TTree *l0tree, Bool_t isdeb, Bool_t isverb){
139      stringstream myquery;
140      myquery.str("");
141      myquery << "SET time_zone='+0:00'";
142      dbc->Query(myquery.str().c_str());
143    //    //
144    debug = isdeb;    debug = isdeb;
145    verbose = isverb;    verbose = isverb;
# Line 169  void CaloLevel0::ProcessingInit(TSQLServ Line 221  void CaloLevel0::ProcessingInit(TSQLServ
221    //    //
222  }  }
223    
224  Int_t CaloLevel0::ChkCalib(TSQLServer *dbc, UInt_t atime){  Int_t CaloLevel0::ChkCalib(GL_TABLES *glt, UInt_t atime){
225    Int_t sgnl = 0;    Int_t sgnl = 0;
226    for ( Int_t s = 0; s < 4; s++){    for ( Int_t s = 0; s < 4; s++){
227      if ( atime > totime[s] ){      if ( atime > totime[s] ){
228        if ( !dbc->IsConnected() ) throw -116;        sgnl = Update(glt,atime,s);
       sgnl = Update(dbc,atime,s);  
229        if ( sgnl < 0 ) return(sgnl);        if ( sgnl < 0 ) return(sgnl);
230      };      };
231    };    };
232    return(sgnl);    return(sgnl);
233  }  }
234    
235  Int_t CaloLevel0::ChkParam(TSQLServer *dbc, UInt_t runheader){  Int_t CaloLevel0::ChkParam(TSQLServer *dbc, UInt_t runheader, Bool_t mechal){
236      Int_t sig = this->ChkParamDo(dbc,runheader,mechal);
237      return(sig);
238    }
239    
240    Int_t CaloLevel0::ChkParam(GL_TABLES *glt, UInt_t runheader, Bool_t mechal){
241      const TString host = glt->CGetHost();
242      const TString user = glt->CGetUser();
243      const TString psw = glt->CGetPsw();
244      TSQLServer *dbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data());
245      if ( !dbc->IsConnected() ) throw -116;
246      stringstream myquery;
247      myquery.str("");
248      myquery << "SET time_zone='+0:00'";
249      dbc->Query(myquery.str().c_str());
250      //
251      Int_t sig = this->ChkParamDo(dbc,runheader,mechal);
252      dbc->Close();
253      delete dbc;
254      return(sig);
255    }
256    
257    Int_t CaloLevel0::ChkParamDo(TSQLServer *dbc, UInt_t runheader, Bool_t mechal){
258      //
259    stringstream calfile;    stringstream calfile;
260    stringstream bmfile;    stringstream bmfile;
261    stringstream aligfile;    stringstream aligfile;
# Line 191  Int_t CaloLevel0::ChkParam(TSQLServer *d Line 265  Int_t CaloLevel0::ChkParam(TSQLServer *d
265    GL_PARAM *glparam = new GL_PARAM();    GL_PARAM *glparam = new GL_PARAM();
266    //    //
267    if ( calopar1 || ( ttcalopar1 != 0 && ttcalopar1 < runheader ) ){    if ( calopar1 || ( ttcalopar1 != 0 && ttcalopar1 < runheader ) ){
268        //
269        if ( debug ) printf(" calopar1 %i ftcalopar1 %u ttcalopar1 %u runheader %u \n",calopar1,ftcalopar1,ttcalopar1,runheader);
270        //
271        if ( calopar1 ){
272          //
273          // determine where I can find calorimeter ADC to MIP conversion file  
274          //
275          if ( verbose ) printf(" Querying DB for calorimeter parameters files...\n");
276          //
277          error = 0;
278          error = glparam->Query_GL_PARAM(runheader,101,dbc);
279          if ( error < 0 ) return(error);
280          //
281          calfile.str("");
282          calfile << glparam->PATH.Data() << "/";
283          calfile << glparam->NAME.Data();
284          //
285          if ( verbose ) printf("\n Using ADC to MIP conversion file: \n %s \n",calfile.str().c_str());
286          f = fopen(calfile.str().c_str(),"rb");
287          if ( !f ){
288            if ( verbose ) printf(" CALORIMETER - ERROR: no ADC to MIP file!\n");
289            return(-105);
290          };
291          //
292          for (Int_t m = 0; m < 2 ; m++ ){
293            for (Int_t k = 0; k < 22; k++ ){
294              for (Int_t l = 0; l < 96; l++ ){
295                fread(&mip[m][k][l],sizeof(mip[m][k][l]),1,f);
296                if ( debug ) printf(" %f \n",mip[m][k][l]);
297              };
298            };
299          };
300          fclose(f);      
301        };
302        //
303      calopar1 = false;      calopar1 = false;
304      //      //
305      // determine where I can find calorimeter ADC to MIP conversion file        // flight extra corrections:
306      //      //
307      if ( verbose ) printf(" Querying DB for calorimeter parameters files...\n");      if ( verbose ) printf(" Querying DB for calorimeter flight ADC to MIP files...\n");
308      //      //
309      error = 0;      error = 0;
310      error = glparam->Query_GL_PARAM(runheader,101,dbc);      error = glparam->Query_GL_PARAM(runheader,110,dbc);
311      if ( error < 0 ) return(error);      if ( error < 0 ) return(error);
312      //      //
313      calfile.str("");      calfile.str("");
# Line 207  Int_t CaloLevel0::ChkParam(TSQLServer *d Line 316  Int_t CaloLevel0::ChkParam(TSQLServer *d
316      ftcalopar1 = glparam->FROM_TIME;      ftcalopar1 = glparam->FROM_TIME;
317      ttcalopar1 = glparam->TO_TIME;      ttcalopar1 = glparam->TO_TIME;
318      //      //
319      if ( verbose ) printf("\n Using ADC to MIP conversion file: \n %s \n",calfile.str().c_str());      if ( verbose ) printf("\n Using ADC to MIP special conversion file: \n %s \n",calfile.str().c_str());
320      f = fopen(calfile.str().c_str(),"rb");      ifstream spfile;
321      if ( !f ){      spfile.open(calfile.str().c_str());
322        if ( verbose ) printf(" CALORIMETER - ERROR: no ADC to MIP file!\n");      if ( !spfile ){
323        return(-105);        if ( verbose ) printf(" CALORIMETER - ERROR: no special calibration file!\n");
324          return(-123);
325      };      };
326        //  
327        Int_t vview = 0;
328        Int_t vplane = 0;
329        Int_t vstrip = 0;
330        Float_t vval = 0.;
331        while ( spfile >> vview && spfile >> vplane && spfile >> vstrip && spfile >> vval){
332          if ( debug ) printf(" Setting ADC to MIP conversion factor: view %i plane %i strip %i mip %f  \n",vview,vplane,vstrip,vval);
333          mip[vview][vplane][vstrip] = vval;
334        };
335      //      //
     for (Int_t m = 0; m < 2 ; m++ ){  
       for (Int_t k = 0; k < 22; k++ ){  
         for (Int_t l = 0; l < 96; l++ ){  
           fread(&mip[m][k][l],sizeof(mip[m][k][l]),1,f);  
         };  
       };  
     };  
     fclose(f);  
336    };    };
337    //    //
338      //
339    if ( calopar2 || ( ttcalopar2 != 0 && ttcalopar2 < runheader ) ){    if ( calopar2 || ( ttcalopar2 != 0 && ttcalopar2 < runheader ) ){
340        //
341        if ( debug ) printf(" calopar2 %i ftcalopar2 %u ttcalopar2 %u runheader %u \n",calopar2,ftcalopar2,ttcalopar2,runheader);
342      calopar2 = false;      calopar2 = false;
343      //      //
344      // determine where I can find calorimeter alignment file        // determine where I can find calorimeter alignment file  
# Line 240  Int_t CaloLevel0::ChkParam(TSQLServer *d Line 354  Int_t CaloLevel0::ChkParam(TSQLServer *d
354      ftcalopar2 = glparam->FROM_TIME;      ftcalopar2 = glparam->FROM_TIME;
355      ttcalopar2 = glparam->TO_TIME;      ttcalopar2 = glparam->TO_TIME;
356      //      //
357      if ( verbose ) printf("\n Using alignment file: \n %s \n",aligfile.str().c_str());      if ( verbose ) printf("\n Using parameter file: \n %s \n",aligfile.str().c_str());
358      f = fopen(aligfile.str().c_str(),"rb");      f = fopen(aligfile.str().c_str(),"rb");
359      if ( !f ){      if ( !f ){
360        if ( verbose ) printf(" CALORIMETER - ERROR: no alignement file!\n");        if ( verbose ) printf(" CALORIMETER - ERROR: no parameter file!\n");
361        return(-106);        return(-106);
362      };      };
363      //      //
364      fread(&clevel1->xalig,sizeof(clevel1->xalig),1,f);      if ( !mechal ){
365      if ( debug ) printf(" xalig = %f \n",clevel1->xalig);        //
366      fread(&clevel1->yalig,sizeof(clevel1->yalig),1,f);        fread(&clevel1->xalig,sizeof(clevel1->xalig),1,f);
367      if ( debug ) printf(" yalig = %f \n",clevel1->yalig);        if ( debug ) printf(" xalig = %f \n",clevel1->xalig);
368      fread(&clevel1->zalig,sizeof(clevel1->zalig),1,f);        fread(&clevel1->yalig,sizeof(clevel1->yalig),1,f);
369      if ( debug ) printf(" zalig = %f \n",clevel1->zalig);        if ( debug ) printf(" yalig = %f \n",clevel1->yalig);
370          fread(&clevel1->zalig,sizeof(clevel1->zalig),1,f);
371          if ( debug ) printf(" zalig = %f \n",clevel1->zalig);
372        } else {
373          if ( verbose ) printf("\n Using MECHANICAL alignement parameters \n");
374          //
375          CaloStrip cs = CaloStrip();
376          cs.UseMechanicalAlig();
377          clevel1->xalig = cs.GetXalig();
378          if ( debug ) printf(" xalig = %f \n",clevel1->xalig);
379          clevel1->yalig = cs.GetYalig();
380          if ( debug ) printf(" yalig = %f \n",clevel1->yalig);
381          clevel1->zalig = cs.GetZalig();
382          if ( debug ) printf(" zalig = %f \n",clevel1->zalig);
383          //
384          Float_t tmp = 0;
385          fread(&tmp,sizeof(clevel1->xalig),1,f);
386          fread(&tmp,sizeof(clevel1->yalig),1,f);
387          fread(&tmp,sizeof(clevel1->zalig),1,f);
388          //      clevel1->zalig = -265.82;
389          //    
390        };
391      fread(&clevel1->emin,sizeof(clevel1->emin),1,f);        fread(&clevel1->emin,sizeof(clevel1->emin),1,f);  
392      if ( debug ) printf(" signal threshold = %f \n",clevel1->emin);      if ( debug ) printf(" signal threshold = %f \n",clevel1->emin);
393      //      //
# Line 262  Int_t CaloLevel0::ChkParam(TSQLServer *d Line 397  Int_t CaloLevel0::ChkParam(TSQLServer *d
397    // Load offline bad strip mask    // Load offline bad strip mask
398    //    //
399    if ( calopar3 || ( ttcalopar3 != 0 && ttcalopar3 < runheader ) ){    if ( calopar3 || ( ttcalopar3 != 0 && ttcalopar3 < runheader ) ){
400        if ( debug ) printf(" calopar3 %i ftcalopar3 %u ttcalopar3 %u runheader %u \n",calopar3,ftcalopar3,ttcalopar3,runheader);
401      calopar3 = false;      calopar3 = false;
402      //      //
403      // determine where I can find calorimeter alignment file        // determine where I can find calorimeter alignment file  
# Line 308  Int_t CaloLevel0::ChkParam(TSQLServer *d Line 444  Int_t CaloLevel0::ChkParam(TSQLServer *d
444      badfile.close();      badfile.close();
445    };    };
446    //    //
447      // calopar4
448      //
449      if ( calopar4 || ( ttcalopar4 != 0 && ttcalopar4 < runheader ) ){
450        //
451        if ( debug ) printf(" calopar4 %i ftcalopar4 %u ttcalopar4 %u runheader %u \n",calopar4,ftcalopar4,ttcalopar4,runheader);
452        //
453        calopar4 = false;
454        //
455        // flight extra corrections:
456        //
457        if ( verbose ) printf(" Querying DB for calorimeter max rms file...\n");
458        //
459        error = 0;
460        error = glparam->Query_GL_PARAM(runheader,109,dbc);
461        if ( error < 0 ) return(error);
462        //
463        calfile.str("");
464        calfile << glparam->PATH.Data() << "/";
465        calfile << glparam->NAME.Data();
466        ftcalopar4 = glparam->FROM_TIME;
467        ttcalopar4 = glparam->TO_TIME;
468        //
469        if ( verbose ) printf("\n Using calorimeter max rms file: \n %s \n",calfile.str().c_str());
470        ifstream spfile;
471        spfile.open(calfile.str().c_str());
472        if ( !spfile ){
473          if ( verbose ) printf(" CALORIMETER - ERROR: no max rms file!\n");
474          return(-124);
475        };
476        //  
477        Int_t vview = 0;
478        Int_t vplane = 0;
479        Int_t vval = 0;
480        for (Int_t l=0; l<2; l++){
481          for (Int_t m=0; m<22; m++){
482            maxrms[l][m] = 26;
483          };
484        };
485        while ( spfile >> vview && spfile >> vplane && spfile >> vval){
486          if ( debug ) printf(" Setting view %i plane %i max rms %i  \n",vview,vplane,vval);
487          maxrms[vview][vplane] = vval;
488        };
489        spfile.close();
490        //
491      };
492      //
493      // calopar5
494      //
495      if ( calopar5 || ( ttcalopar5 != 0 && ttcalopar5 < runheader ) ){
496        //
497        if ( debug ) printf(" calopar5 %i ftcalopar5 %u ttcalopar5 %u runheader %u \n",calopar5,ftcalopar5,ttcalopar5,runheader);
498        //
499        calopar5 = false;
500        //
501        // flight extra corrections:
502        //
503        if ( verbose ) printf(" Querying DB for calorimeter noise to signal threshold file...\n");
504        //
505        error = 0;
506        error = glparam->Query_GL_PARAM(runheader,111,dbc);
507        if ( error < 0 ) return(error);
508        //
509        calfile.str("");
510        calfile << glparam->PATH.Data() << "/";
511        calfile << glparam->NAME.Data();
512        ftcalopar5 = glparam->FROM_TIME;
513        ttcalopar5 = glparam->TO_TIME;
514        //
515        if ( verbose ) printf("\n Using calorimeter noise to signal threshold file: \n %s \n",calfile.str().c_str());
516        ifstream spfile;
517        spfile.open(calfile.str().c_str());
518        if ( !spfile ){
519          if ( verbose ) printf(" CALORIMETER - ERROR: no noise to signal threshold file!\n");
520          return(-125);
521        };
522        //  
523        Int_t vview = 0;
524        Int_t vplane = 0;
525        Int_t vstrip = 0;
526        Float_t vval = 0.;
527        for (Int_t l=0; l<2; l++){
528          for (Int_t m=0; m<22; m++){
529            for (Int_t n=0; n<96; n++){
530              memin[l][m][n] = 0.7;
531            };
532          };
533        };
534        while ( spfile >> vview && spfile >> vplane && spfile >> vstrip && spfile >> vval){
535          if ( vstrip == -1 ){
536            for (Int_t ll=0; ll<96; ll++){
537              if ( debug ) printf(" Setting view %i plane %i strip %i noise to signal ratio %f  \n",vview,vplane,ll,vval);
538              memin[vview][vplane][ll] = vval;
539            };
540          } else {
541            if ( debug ) printf(" Setting view %i plane %i strip %i noise to signal ratio %f  \n",vview,vplane,vstrip,vval);
542            memin[vview][vplane][vstrip] = vval;
543          };
544        };
545        spfile.close();
546        //
547      };
548      //
549      //
550    delete glparam;    delete glparam;
551    //    //
552    return(0);    return(0);
553  }  }
554    
555    Int_t CaloLevel0::CalcCrossTalkCorr(TSQLServer *dbc, UInt_t runheader, Bool_t ctusetable){
556      Int_t sig = CalcCrossTalkCorrDo(dbc,runheader,ctusetable);
557      return(sig);  
558    };
559    
560    Int_t CaloLevel0::CalcCrossTalkCorr(TSQLServer *dbc, UInt_t runheader){
561      Int_t sig = CalcCrossTalkCorrDo(dbc,runheader,true);
562      return(sig);
563    }
564    
565  void CaloLevel0::FindBaseRaw(Int_t l, Int_t m, Int_t pre){  Int_t CaloLevel0::CalcCrossTalkCorr(GL_TABLES *glt, UInt_t runheader, Bool_t usetable){
566      Float_t minstrip = 100000.;    const TString host = glt->CGetHost();
567      Float_t rms = 0.;    const TString user = glt->CGetUser();
568      base[l][m][pre] = 0.;    const TString psw = glt->CGetPsw();
569      for (Int_t e = pre*16; e < (pre+1)*16 ; e++){    TSQLServer *dbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data());
570          if ( calgood[l][m][e] == 0. && obadmask[l][m][e] == 0 &&  dexy[l][m][e]-calped[l][m][e] < minstrip &&  dexy[l][m][e] > 0.) {    if ( !dbc->IsConnected() ) throw -116;
571              minstrip = dexy[l][m][e]-calped[l][m][e];    stringstream myquery;
572              rms = calthr[l][m][pre];    myquery.str("");
573          };    myquery << "SET time_zone='+0:00'";
574      dbc->Query(myquery.str().c_str());
575      //
576      Int_t sig = CalcCrossTalkCorrDo(dbc,runheader,usetable);
577      dbc->Close();
578      delete dbc;
579      //
580      return(sig);
581      //  
582    };
583    
584    Int_t CaloLevel0::CalcCrossTalkCorr(GL_TABLES *glt, UInt_t runheader){
585      const TString host = glt->CGetHost();
586      const TString user = glt->CGetUser();
587      const TString psw = glt->CGetPsw();
588      TSQLServer *dbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data());
589      if ( !dbc->IsConnected() ) throw -116;
590      stringstream myquery;
591      myquery.str("");
592      myquery << "SET time_zone='+0:00'";
593      dbc->Query(myquery.str().c_str());
594      //
595      Int_t sig = CalcCrossTalkCorrDo(dbc,runheader,true);
596      dbc->Close();
597      delete dbc;
598      //
599      return(sig);
600      //
601    }
602    
603    Int_t CaloLevel0::CalcCrossTalkCorrDo(TSQLServer *dbc, UInt_t runheader, Bool_t usetable){
604      //
605      if ( ctground ) return(0);
606      //
607      Int_t error = 0;
608      GL_PARAM *glparam = new GL_PARAM();  
609      //
610      // determine where I can find file with offline bad pulser mask
611      //
612      stringstream bmfile;
613      error = 0;
614      error = glparam->Query_GL_PARAM(runheader,105,dbc);
615      if ( error < 0 ) return(error);
616      //
617      bmfile.str("");
618      bmfile << glparam->PATH.Data() << "/";
619      bmfile << glparam->NAME.Data();
620      //
621      ifstream badfile;
622      if ( verbose ) printf("\n Using bad pulser offline mask file: \n %s \n\n",bmfile.str().c_str());
623      badfile.open(bmfile.str().c_str());
624      if ( !badfile ){
625        if ( verbose ) printf(" CALORIMETER - ERROR: no bad pulser offline mask file!\n");
626        return(-115);
627      };
628      //
629      Bool_t isdone = false;
630      Int_t bad = 0;
631      Int_t view = 1;
632      Int_t pre = 0;
633      Int_t plane = 21;
634      while ( !isdone ) {
635        badfile >> bad;
636        obadpulsemask[view][plane][pre] = bad;
637        if ( debug && bad ) printf(" SETTING view %i plane %i pre %i  BAD = %i \n",view,plane,pre,bad);
638        pre++;
639        if ( pre > 5 ){
640          pre = 0;
641          plane--;
642          if ( plane < 0 ){
643            plane = 21;
644            view--;
645          };
646          if ( view < 0 ) isdone = true;
647      };      };
648      if ( minstrip != 100000. ) {    };
649          Float_t strip6s = 0.;    //
650          for (Int_t e = pre*16; e < (pre+1)*16 ; e++){    badfile.close();
651              if ( (dexy[l][m][e]-calped[l][m][e]) >= minstrip && (dexy[l][m][e]-calped[l][m][e]) <= (minstrip+rms) ) {    if ( !usetable ){
652                  strip6s += 1.;      //
653                  base[l][m][pre] += (dexy[l][m][e] - calped[l][m][e]);      // Let's start with cross-talk correction calculation
654        //
655        GL_CALOPULSE_CALIB *glp = new GL_CALOPULSE_CALIB();
656        Float_t adcp[2][22][96];
657        Float_t adcpcal[2][22][96];
658        memset(adcp , 0, 2*22*96*sizeof(Float_t));
659        memset(adcpcal , 0, 2*22*96*sizeof(Float_t));
660        //
661        UInt_t pampli = 0;  
662        for (Int_t s=0; s<4; s++){
663          //
664          // Save into matrix adcp the values of the highest pulse calibration (pulse amplitude = 2)
665          //
666          pampli = 2;
667          error = 0;
668          error = glp->Query_GL_CALOPULSE_CALIB(runheader,s,pampli,dbc);
669          if ( error < 0 ){
670            if ( verbose ) printf(" CALORIMETER - ERROR: error from GLTables\n");
671            return(error);
672          };
673          //
674          UInt_t idcalib = glp->ID_ROOT_L0;
675          UInt_t fromtime = glp->FROM_TIME;
676          UInt_t calibno = glp->EV_ROOT;
677          //
678          // determine path and name and entry of the calibration file
679          //
680          GL_ROOT *glroot = new GL_ROOT();  
681          if ( verbose ) printf("\n");
682          if ( verbose ) printf(" ** SECTION %i **\n",s);
683          //
684          error = 0;
685          error = glroot->Query_GL_ROOT(idcalib,dbc);
686          if ( error < 0 ){
687            if ( verbose ) printf(" CALORIMETER - ERROR: error from GLTables\n");
688            return(error);
689          };
690          //        
691          stringstream name;
692          name.str("");
693          name << glroot->PATH.Data() << "/";
694          name << glroot->NAME.Data();
695          //
696          TString fcalname = (TString)name.str().c_str();
697          ifstream myfile;
698          myfile.open(fcalname.Data());
699          if ( !myfile ){    
700            return(-107);
701          };
702          myfile.close();
703          //
704          TFile *File = new TFile(fcalname.Data());
705          if ( !File ) return(-108);
706          TTree *tr = (TTree*)File->Get("CalibCalPulse2");
707          if ( !tr ) return(-119);
708          //
709          TBranch *calo = tr->GetBranch("CalibCalPulse2");
710          //
711          pamela::CalibCalPulse2Event *ce = 0;
712          tr->SetBranchAddress("CalibCalPulse2", &ce);
713          //
714          Long64_t ncalibs = calo->GetEntries();
715          //
716          if ( !ncalibs ) return(-110);
717          //
718          calo->GetEntry(calibno);
719          if ( verbose ) printf(" PULSE2 using entry %u from file %s",calibno,fcalname.Data());
720          //
721          // retrieve calibration table
722          //
723          if ( ce->pstwerr[s] && ce->pperror[s] == 0 && ce->unpackError == 0 ){
724            for ( Int_t d=0 ; d<11 ;d++  ){
725              for ( Int_t j=0; j<96 ;j++){
726                if ( s == 2 ){
727                  adcp[0][2*d+1][j] = ce->calpuls[3][d][j];
728              };              };
729              //              if ( s == 3 ){
730              //  compression                adcp[0][2*d][j] = ce->calpuls[1][d][j];
             //  
             if ( abs((int)(dexy[l][m][e]-calped[l][m][e])) <= (minstrip+rms) ) {  
                 dexyc[l][m][e] = 0.;  
             } else {  
                 dexyc[l][m][e] = dexy[l][m][e];  
731              };              };
732                if ( s == 0 ){
733                  adcp[1][2*d][j] = ce->calpuls[0][d][j];
734                };
735                if ( s == 1 ){
736                  adcp[1][2*d+1][j] = ce->calpuls[2][d][j];
737                };
738              };
739          };          };
740          if ( strip6s >= 9. ){              } else {
741              Double_t arro = base[l][m][pre]/strip6s;          if ( verbose ) printf(" CALORIMETER - ERROR: problems finding a good calibration in this file! \n\n ");
742              Float_t deci = 1000.*((float)arro - float(int(arro)));                              return(-111);
743              if ( deci < 500. ) {        };
744                  arro = double(int(arro));        //
745          File->Close();
746          delete glroot;
747          //
748          // Save into matrix adcpcal the calibrated values of the pulse calibration (subtraction of pulse amplitude = 0 relative to the pulse2 calibration used)
749          //
750          pampli = 0;
751          error = 0;
752          error = glp->Query_GL_CALOPULSE_CALIB(fromtime,s,pampli,dbc);
753          if ( error < 0 ){
754            if ( verbose ) printf(" CALORIMETER - ERROR: error from GLTables\n");
755            return(error);
756          };
757          //
758          idcalib = glp->ID_ROOT_L0;
759          calibno = glp->EV_ROOT;
760          //
761          // determine path and name and entry of the calibration file
762          //
763          glroot = new GL_ROOT();  
764          if ( verbose ) printf("\n");
765          if ( verbose ) printf(" ** SECTION %i **\n",s);
766          //
767          error = 0;
768          error = glroot->Query_GL_ROOT(idcalib,dbc);
769          if ( error < 0 ){
770            if ( verbose ) printf(" CALORIMETER - ERROR: error from GLTables\n");
771            return(error);
772          };
773          //        
774          name.str("");
775          name << glroot->PATH.Data() << "/";
776          name << glroot->NAME.Data();
777          //
778          fcalname = (TString)name.str().c_str();
779          myfile.open(fcalname.Data());
780          if ( !myfile ){    
781            return(-107);
782          };
783          myfile.close();
784          //
785          TFile *File1 = new TFile(fcalname.Data());
786          if ( !File1 ) return(-108);
787          TTree *tr1 = (TTree*)File1->Get("CalibCalPulse1");
788          if ( !tr1 ) return(-120);
789          //
790          TBranch *calo1 = tr1->GetBranch("CalibCalPulse1");
791          //
792          pamela::CalibCalPulse1Event *ce1 = 0;
793          tr1->SetBranchAddress("CalibCalPulse1", &ce1);
794          //
795          ncalibs = calo1->GetEntries();
796          //
797          if ( !ncalibs ) return(-110);
798          //
799          calo1->GetEntry(calibno);
800          if ( verbose ) printf(" PULSE1 using entry %u from file %s",calibno,fcalname.Data());
801          //
802          // retrieve calibration table
803          //
804          if ( ce1->pstwerr[s] && ce1->pperror[s] == 0 && ce1->unpackError == 0 ){
805            for ( Int_t d=0 ; d<11 ;d++  ){
806              for ( Int_t j=0; j<96 ;j++){
807                if ( s == 2 ){
808                  adcpcal[0][2*d+1][j] = adcp[0][2*d+1][j] - ce1->calpuls[3][d][j];
809                };
810                if ( s == 3 ){
811                  adcpcal[0][2*d][j] = adcp[0][2*d][j] - ce1->calpuls[1][d][j];
812                };
813                if ( s == 0 ){
814                  adcpcal[1][2*d][j] = adcp[1][2*d][j] - ce1->calpuls[0][d][j];
815                };
816                if ( s == 1 ){
817                  adcpcal[1][2*d+1][j] = adcp[1][2*d+1][j] - ce1->calpuls[2][d][j];
818                };
819              };
820            };
821          } else {
822            if ( verbose ) printf(" CALORIMETER - ERROR: problems finding a good calibration in this file! \n\n ");
823            return(-111);
824          };
825          //
826          File1->Close();
827          //
828          delete glroot;
829          //
830        };// loop on the four sections
831        //
832        //
833        delete glp;
834        //
835        // Ok, now we can try to calculate the cross-talk correction for each pre-amplifier
836        //
837        for ( Int_t v=0; v<2; v++){
838          if ( debug ) printf(" \n\n NEW VIEW \n");
839          for ( Int_t p=0; p<22; p++){
840            for ( Int_t npre=0; npre<6; npre++){
841              ctprecor[v][p][npre] = 1000.;
842              ctneigcor[v][p][npre] = 1000.;
843              Int_t str0=npre*16;
844              Int_t str16= -1 + (1+npre)*16;
845              //
846              UInt_t neigc = 0;
847              UInt_t farc = 0;
848              UInt_t pulsc = 0;
849              Float_t sigpulsed = 0.;
850              Float_t neigbase = 0.;
851              Float_t farbase = 0.;
852              //
853              // Loop over the strip of the pre and sum all signal far away from pulsed strip, signal in the neighbour(s) strip(s) and save the pulsed signal
854              // moreover count the number of strips in each case
855              //
856              for (Int_t s=str0; s<=str16; s++){
857                if ( adcpcal[v][p][s] > 10000.){
858                  sigpulsed = adcpcal[v][p][s];
859                  pulsc++;
860                  if ( s > str0 ){
861                    neigbase += adcpcal[v][p][s-1];
862                    neigc++;
863                    farbase -= adcpcal[v][p][s-1];
864                    farc--;
865                  };
866                  if ( s < str16 ){
867                    neigbase += adcpcal[v][p][s+1];
868                    neigc++;
869                    farbase -= adcpcal[v][p][s+1];
870                    farc--;
871                  };
872              } else {              } else {
873                  arro = 1. + double(int(arro));                farc++;
874                  farbase += adcpcal[v][p][s];
875              };              };
876              base[l][m][pre] = arro;            };
877          } else {            //
878              base[l][m][pre] = 31000.;            // Now calculate the corrections
879              for (Int_t e = pre*16; e < (pre+1)*16 ; e++){            //
880                  dexyc[l][m][e] = dexy[l][m][e];            Float_t avefarbase = 0.;
881              if ( farc ) avefarbase = farbase/(Float_t)farc;
882              Float_t aveneigbase = 0.;
883              if ( neigc ) aveneigbase = neigbase/(Float_t)neigc;
884              //
885              if ( pulsc == 1 && farc && neigc ){
886                ctprecor[v][p][npre] = -avefarbase/(sigpulsed+fabs(avefarbase));
887                ctneigcor[v][p][npre] = fabs(aveneigbase-avefarbase)/(sigpulsed+fabs(avefarbase));    
888                if ( debug ) printf(" Cross-talk correction View %i Plane %i Pre %i : pre-correction: %f neighbour strips correction %f \n",v,p,npre,ctprecor[v][p][npre],ctneigcor[v][p][npre]);
889              } else {
890                //
891                // did not find the pulsed strip or more than one pulsed strip found!
892                //
893                if ( debug ) printf(" Problems finding the cross-talk corrections: \n View %i Plane %i Pre %i number of pulsed strip %i \n Average faraway baseline %f number of strips %i Average neighbour baseline %f number of neighbour strips %i \n",v,p,npre,pulsc,avefarbase,farc,aveneigbase,neigc);
894                //
895              };
896            };
897            if ( debug ) printf(" \n ==================== \n");
898          };
899        };
900      } else {
901        //
902        // use pre-amply table
903        //    
904        //
905        // determine where I can find file with offline neighbour correction table
906        //
907        stringstream bmfile2;
908        error = 0;
909        error = glparam->Query_GL_PARAM(runheader,106,dbc);
910        if ( error < 0 ) return(error);
911        //
912        bmfile2.str("");
913        bmfile2 << glparam->PATH.Data() << "/";
914        bmfile2 << glparam->NAME.Data();
915        //
916        ifstream badfile2;
917        if ( verbose ) printf("\n Using pre-amply neighbour crosstalk table file: \n %s \n\n",bmfile2.str().c_str());
918        badfile2.open(bmfile2.str().c_str());
919        if ( !badfile2 ){
920          if ( verbose ) printf(" CALORIMETER - ERROR: no pre-amply neighbour crosstalk table file!\n");
921          return(-121);
922        };
923        //  
924        Int_t vview = 0;
925        Int_t vplane = 0;
926        Int_t vpre = 0;
927        Float_t vcorr = 0.;
928        while ( badfile2 >> vview && badfile2 >> vplane && badfile2 >> vpre && badfile2 >> vcorr){
929          if ( debug ) printf(" Pre-amply neighbour correction: view %i plane %i pre %i correction %f  \n",vview,vplane,vpre,vcorr);
930          ctneigcor[vview][vplane][vpre] = vcorr;
931        };
932        //
933        // determine where I can find file with offline SECOND neighbour correction table
934        //
935        stringstream bmfile3;
936        error = 0;
937        error = glparam->Query_GL_PARAM(runheader,107,dbc);
938        if ( error < 0 ) return(error);
939        //
940        bmfile3.str("");
941        bmfile3 << glparam->PATH.Data() << "/";
942        bmfile3 << glparam->NAME.Data();
943        //
944        ifstream badfile3;
945        if ( verbose ) printf("\n Using pre-amply second neighbour crosstalk table file: \n %s \n\n",bmfile3.str().c_str());
946        badfile3.open(bmfile3.str().c_str());
947        if ( !badfile3 ){
948          if ( verbose ) printf(" CALORIMETER - ERROR: no pre-amply second neighbour crosstalk table file!\n");
949          return(-122);
950        };
951        //  
952        Int_t pview = 0;
953        Int_t pplane = 0;
954        Int_t ppre = 0;
955        Float_t pcorr = 0.;
956        while ( badfile3 >> pview && badfile3 >> pplane && badfile3 >> ppre && badfile3 >> pcorr){
957          if ( debug ) printf(" Pre-amply second neighbour correction: view %i plane %i pre %i correction %f  \n",pview,pplane,ppre,-pcorr);
958          ctprecor[pview][pplane][ppre] = -pcorr; // data are saved as negatives in the file
959        };
960        //
961        // determine where to find the file containing the Silicon crosstalk correction table
962        //
963        stringstream bmfile4;
964        error = 0;
965        error = glparam->Query_GL_PARAM(runheader,108,dbc);
966        if ( error < 0 ) return(error);
967        //
968        bmfile4.str("");
969        bmfile4 << glparam->PATH.Data() << "/";
970        bmfile4 << glparam->NAME.Data();
971        //
972        ifstream badfile4;
973        if ( verbose ) printf("\n Using Silicon crosstalk table file: \n %s \n\n",bmfile4.str().c_str());
974        badfile4.open(bmfile4.str().c_str());
975        if ( !badfile4 ){
976          if ( verbose ) printf(" CALORIMETER - ERROR: no Silicon crosstalk table file!\n");
977          return(-125);
978        };
979        //  
980        Int_t spview = 0;
981        Int_t spplane = 0;
982        Int_t psil = 0;
983        Float_t spcorr = 0.;
984        memset(ctsicor, 0, 2*22*9*sizeof(Float_t));
985        while ( badfile4 >> spview && badfile4 >> spplane && badfile4 >> psil && badfile4 >> spcorr){
986          if ( debug ) printf(" Silicon correction: view %i plane %i silicon %i correction %f  \n",spview,spplane,psil,-spcorr);
987          ctsicor[spview][spplane][psil] = -spcorr; // data are saved as negatives in the file
988        };
989        //
990      };
991      //
992      delete glparam;
993      //
994      // Check the calculated corrections
995      //
996      Int_t opre=0;
997      Int_t ppre=0;
998      Bool_t found = false;
999      for ( Int_t v=0; v<2; v++){
1000        for ( Int_t p=0; p<22; p++){
1001          for ( Int_t npre=0; npre<6; npre++){
1002            if ( ctprecor[v][p][npre] == 1000. || ctneigcor[v][p][npre] == 1000. || obadpulsemask[v][p][npre] != 0 ){
1003              if ( debug ) printf(" Cross-talk correction CHANGED for view %i Plane %i Pre %i\n    BEFORE: pre-correction: %f neighbour strips correction %f \n",v,p,npre,ctprecor[v][p][npre],ctneigcor[v][p][npre]);
1004              if ( npre%2 ){
1005                opre = npre-1;
1006              } else {
1007                opre = npre+1;          
1008              };
1009              if ( ctprecor[v][p][opre] == 1000. || ctneigcor[v][p][opre] == 1000. || obadpulsemask[v][p][opre] != 0 ){
1010                ppre=0;
1011                found = false;
1012                while ( ppre < 6 ){
1013                  if ( ctprecor[v][p][ppre] != 1000. && ctneigcor[v][p][ppre] != 1000. && !obadpulsemask[v][p][ppre] ){
1014                    found = true;
1015                    ctprecor[v][p][npre] = ctprecor[v][p][ppre];
1016                    ctneigcor[v][p][npre] = ctneigcor[v][p][ppre];
1017                    break;
1018                  };
1019                  ppre++;
1020                };
1021                if ( !found ){
1022                  if ( verbose ) printf(" WARNING: cannot find a good cross-talk correction for view %i plane %i pre %i \n Setting to default values 0.002 0.002\n",v,p,npre);
1023                  ctprecor[v][p][npre] = 0.002;
1024                  ctneigcor[v][p][npre] = 0.002;
1025              };              };
1026              } else {
1027                ctprecor[v][p][npre] = ctprecor[v][p][opre];
1028                ctneigcor[v][p][npre] = ctneigcor[v][p][opre];
1029              };
1030              if ( debug ) printf("    AFTER: pre-correction: %f neighbour strips correction %f \n",ctprecor[v][p][npre],ctneigcor[v][p][npre]);
1031          };          };
1032      } else {        };
1033        };
1034      };  
1035      //
1036      return(0);
1037    }
1038    
1039    void CaloLevel0::FindBaseCompress(Int_t l, Int_t m, Int_t pre){
1040      Int_t n = 0;
1041      Float_t q = 0;
1042      this->FindBaseCompress(l,m,pre,n,q);
1043    }
1044    
1045    void CaloLevel0::FindBaseCompress(Int_t l, Int_t m, Int_t pre, Int_t &nst, Float_t &qp){
1046      for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
1047        dexy[l][m][e] = dexyc[l][m][e];
1048      };  
1049      this->FindBaseRaw(l,m,pre,nst,qp);
1050    }
1051    
1052    void CaloLevel0::FindBaseRaw(Int_t l, Int_t m, Int_t pre){
1053      Int_t n = 0;
1054      Float_t q = 0;
1055      this->FindBaseRaw(l,m,pre,n,q);
1056    }
1057    
1058    void CaloLevel0::FindBaseRaw(Int_t l, Int_t m, Int_t pre, Int_t &nst, Float_t &qp){
1059      //
1060      Float_t minstrip = 100000.;
1061      Float_t rms = 0.;
1062      Int_t process = 0;
1063      Int_t onlmask[16];
1064      memset(onlmask, 0, 16*sizeof(Int_t));
1065      //
1066      while ( process < 2 ){
1067        //
1068        minstrip = 100000.;
1069        rms = 0.;
1070        base[l][m][pre] = 0.;
1071        qp = 0.;
1072        //
1073        Int_t spos = -1;
1074        Int_t ee = 0;
1075        for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
1076          if ( calgood[l][m][e] == 0. && obadmask[l][m][e] == 0 &&  dexy[l][m][e]-calped[l][m][e] < minstrip &&  dexy[l][m][e] > 0. && onlmask[ee] == 0 ) {
1077            minstrip = dexy[l][m][e]-calped[l][m][e];
1078            rms = calthr[l][m][pre];
1079            spos = ee;
1080          };
1081          ee++;
1082          qp += (dexy[l][m][e]-calped[l][m][e]-sbase[l][m][e]);
1083        };
1084        //
1085        if ( debug && l==0 ){
1086          printf("\n BASELINE CALCULATION for view %i pl %i pre %i: \n => minstrip %f rms %f \n => qp = %f \n",l,m,pre,minstrip,rms,qp);
1087        };
1088        if ( minstrip != 100000. ) {
1089          Float_t strip6s = 0.;
1090          for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
1091            if ( (dexy[l][m][e]-calped[l][m][e]) >= minstrip && (dexy[l][m][e]-calped[l][m][e]) <= (minstrip+rms) ) {
1092              strip6s += 1.;
1093              base[l][m][pre] += (dexy[l][m][e] - calped[l][m][e]);
1094            };
1095            //
1096            //  compression
1097            //
1098            //      if ( abs((int)(dexy[l][m][e]-calped[l][m][e])) <= (minstrip+rms) ) {
1099            //        dexyc[l][m][e] = 0.;
1100            //      } else {
1101            dexyc[l][m][e] = dexy[l][m][e];
1102              //    };
1103          };
1104          //
1105          if ( strip6s == 1. && process < 1 ){
1106            onlmask[spos] = 1;
1107            process++;
1108            if ( debug ) printf(" Warning, only one strip to calculate baseline: minstrip %f rms %f spos %i l %i m %i pre %i \n",minstrip,rms,spos,l,m,pre);
1109            continue;
1110          };
1111          process += 2;
1112          nst = (Int_t)strip6s;
1113          //
1114          if ( debug ){
1115            printf(" strip6s %f \n",strip6s);
1116          };
1117          //        if ( strip6s >= 9. ){      
1118          if ( (strip6s >= 2. && process == 2) || (strip6s >= 9. && process > 2) ){    
1119            //if ( (strip6s >= 4. && process == 2) || (strip6s >= 9. && process > 2) ){        
1120            Double_t arro = base[l][m][pre]/strip6s;
1121            Float_t deci = 1000.*((float)arro - float(int(arro)));                
1122            if ( deci < 500. ) {
1123              arro = double(int(arro));
1124            } else {
1125              arro = 1. + double(int(arro));
1126            };
1127            base[l][m][pre] = arro;
1128            //
1129            // if too few strips were used to determine the baseline check if it is comparable with the previous event, if not mark it as bad
1130            //
1131            if ( debug && process > 2 ) printf(" AGH low strip value was discarded process %i strip6s %f minstrip %f rms %f spos %i\n",process,strip6s,minstrip,rms,spos);
1132            if ( debug ) printf(" Calculated baseline: base %f sbase-0.02*qp %f \n",base[l][m][pre],(-qp*0.02+sbase[l][m][pre]));
1133            //
1134            if ( strip6s < 4 && base[l][m][pre] > (-0.015*qp+sbase[l][m][pre]) && sbase[l][m][pre] > 0. ){
1135              if ( debug ) printf(" Suspicious calculated baseline: base %f sbase-0.02*qp %f strip6s %i \n",base[l][m][pre],(-qp*0.02+sbase[l][m][pre]),(Int_t)strip6s);
1136              base[l][m][pre] = 31000.;    
1137              for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
1138                dexyc[l][m][e] = dexy[l][m][e];
1139              };      
1140            };
1141          } else {
1142          base[l][m][pre] = 31000.;          base[l][m][pre] = 31000.;
1143            for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
1144              dexyc[l][m][e] = dexy[l][m][e];
1145            };
1146          };
1147        } else {
1148          process += 2;
1149          base[l][m][pre] = 31000.;
1150          for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
1151            dexyc[l][m][e] = dexy[l][m][e];
1152          };
1153      };      };
1154      };
1155  }  }
1156    
1157  Int_t CaloLevel0::Calibrate(Int_t ei){  Int_t CaloLevel0::Calibrate(Int_t ei){
# Line 369  Int_t CaloLevel0::Calibrate(Int_t ei){ Line 1162  Int_t CaloLevel0::Calibrate(Int_t ei){
1162    //    //
1163    // if it was not a selftrigger event, could it ever been a selftrigger event? if so trigty = 3.    // if it was not a selftrigger event, could it ever been a selftrigger event? if so trigty = 3.
1164    //    //
1165      clevel2->nsatstrip = 0.;
1166    Int_t val = 0;    Int_t val = 0;
1167    Int_t del = 1100;    Int_t del = 1100;
1168      for (Int_t sec = 0; sec < 4; sec++){
1169        for (Int_t dsec = 0; dsec < 7; dsec++){
1170          val = (Int_t)de->calselftrig[sec][dsec];
1171          del = delay(val);
1172          clevel2->selfdelay[sec][dsec] = del;
1173        };
1174      };
1175      val = 0;
1176      del = 1100;
1177    if ( clevel2->trigty != 2. ){    if ( clevel2->trigty != 2. ){
1178      Bool_t ck = false;      Bool_t ck = false;
1179      for (Int_t sec = 0; sec < 4; sec++){      for (Int_t sec = 0; sec < 4; sec++){
# Line 406  Int_t CaloLevel0::Calibrate(Int_t ei){ Line 1209  Int_t CaloLevel0::Calibrate(Int_t ei){
1209    Float_t ener;    Float_t ener;
1210    Int_t doneb = 0;    Int_t doneb = 0;
1211    Int_t donec = 0;    Int_t donec = 0;
1212    Int_t ck = 0;    Int_t ck[2][22][6];
1213      memset(ck, 0, 2*22*6*sizeof(Int_t));
1214    Int_t ipre = 0;    Int_t ipre = 0;
1215    Int_t ip[3] = {0};    //  Int_t ip[3] = {0};
1216      Int_t ip[3] = {0,0,0};
1217      Int_t ipp = 0;
1218    Float_t base0, base1, base2;    Float_t base0, base1, base2;
1219    base0 = 0.;    base0 = 0.;
1220    base1 = 0.;    base1 = 0.;
1221    base2 = 0.;    base2 = 0.;
1222    Float_t qpre[6] = {0.,0.,0.,0.,0.,0.};    Float_t qpre[2][22][6];
1223      memset(qpre, 0, 2*22*6*sizeof(Float_t));
1224    Float_t ene[96];    Float_t ene[96];
1225    Int_t chdone[4] = {0,0,0,0};    Int_t chdone[4] = {0,0,0,0};
1226    Int_t pe = 0;    Int_t pe = 0;
# Line 421  Int_t CaloLevel0::Calibrate(Int_t ei){ Line 1228  Int_t CaloLevel0::Calibrate(Int_t ei){
1228    Float_t ener0 = 0.;    Float_t ener0 = 0.;
1229    Float_t cbase0 = 0.;    Float_t cbase0 = 0.;
1230    Bool_t pproblem = false;    Bool_t pproblem = false;
1231      Bool_t negbase = false;
1232    //    //
1233    Float_t tim = 0.;    Float_t tim = 0.;
1234    Int_t plo = 0;    Int_t plo = 0;
# Line 432  Int_t CaloLevel0::Calibrate(Int_t ei){ Line 1240  Int_t CaloLevel0::Calibrate(Int_t ei){
1240    for (Int_t l = 0; l < 2; l++){    for (Int_t l = 0; l < 2; l++){
1241      for (Int_t m = 0; m < 22; m++){      for (Int_t m = 0; m < 22; m++){
1242        //        //
1243        // determine the section number        // determine the section number  
1244        //        //
1245          negbase = false;
1246        se = 5;        se = 5;
1247        if (l == 0 && m%2 == 0) se = 3;        if (l == 0 && m%2 == 0) se = 3;
1248        if (l == 0 && m%2 != 0) se = 2;        if (l == 0 && m%2 != 0) se = 2;
1249        if (l == 1 && m%2 == 0) se = 1;        if (l == 1 && m%2 != 0) se = 1;
1250        if (l == 1 && m%2 != 0) se = 0;                  if (l == 1 && m%2 == 0) se = 0;          
1251        //        //
1252        // determine what kind of event we are going to analyze        // determine what kind of event we are going to analyze
1253        //        //
# Line 459  Int_t CaloLevel0::Calibrate(Int_t ei){ Line 1268  Int_t CaloLevel0::Calibrate(Int_t ei){
1268          };          };
1269          clevel2->perr[se] = 0;          clevel2->perr[se] = 0;
1270          if ( de->perror[se] != 0 ){          if ( de->perror[se] != 0 ){
1271            clevel2->perr[se] = 1;            clevel2->perr[se] = de->perror[se];
1272            pe++;            pe++;
1273          };          };
1274          clevel2->swerr[se] = 0;          clevel2->swerr[se] = 0;
# Line 475  Int_t CaloLevel0::Calibrate(Int_t ei){ Line 1284  Int_t CaloLevel0::Calibrate(Int_t ei){
1284          pre = -1;          pre = -1;
1285          //          //
1286          for (Int_t nn = 0; nn < 96; nn++){                            for (Int_t nn = 0; nn < 96; nn++){                  
1287            ene[nn] = 0.;            //      ene[nn] = 0.;
1288            dexy[l][m][nn] = de->dexy[l][m][nn] ;            dexy[l][m][nn] = de->dexy[l][m][nn] ;
1289            dexyc[l][m][nn] = de->dexyc[l][m][nn] ;            dexyc[l][m][nn] = de->dexyc[l][m][nn] ;
1290          };          };
# Line 484  Int_t CaloLevel0::Calibrate(Int_t ei){ Line 1293  Int_t CaloLevel0::Calibrate(Int_t ei){
1293          //          //
1294          pre = -1;          pre = -1;
1295          cbase0 = 0.;          cbase0 = 0.;
1296            Int_t nstt[2];
1297            Float_t rqp[2];
1298          for (Int_t i = 0; i < 3; i++){          for (Int_t i = 0; i < 3; i++){
1299              nstt[0] = 1000;
1300              nstt[1] = 1000;
1301              rqp[0] = 0.;
1302              rqp[1] = 0.;
1303            for (Int_t j = 0; j < 2; j++){            for (Int_t j = 0; j < 2; j++){
1304              pre = j + i*2;              pre = j + i*2;
1305              //              //
1306              // baseline check and calculation              // baseline check and calculation
1307              //              //
1308              if ( !isRAW ) {              if ( !isRAW ){
1309                base[l][m][pre] = de->base[l][m][pre] ;                  //
1310                  // if it is a compress event with fully transmitted pre try to calculate the baseline
1311                  //
1312                  if ( de->base[l][m][pre] != 0. && de->base[l][m][pre]<31000. ) {
1313                    base[l][m][pre] = de->base[l][m][pre] ;  
1314                  } else {
1315                    FindBaseCompress(l,m,pre,nstt[j],rqp[j]);
1316                  };
1317                cbase0 += base[l][m][pre];                cbase0 += base[l][m][pre];
1318              } else {              } else {
1319                //                //
1320                // if it is a raw event and we haven't checked                // if it is a raw event calculate the baseline.
               // yet, calculate the baseline.  
1321                //                //
1322                FindBaseRaw(l,m,pre);                FindBaseRaw(l,m,pre,nstt[j],rqp[j]);
1323                cbase0 += base[l][m][pre];                cbase0 += base[l][m][pre];
1324              };              };
1325            };            };
1326              //
1327              // if we are able to calculate the baseline with more than 3 strips on one pre and not in the other one choose the pre with more calculated strips
1328              //
1329              if ( nstt[0] < 4 && nstt[1] >= 4 && nstt[0] != 1000 && nstt[1] != 1000 ) base[l][m][pre-1] = 31000.;
1330              if ( nstt[0] >= 4 && nstt[1] < 4 && nstt[0] != 1000 && nstt[1] != 1000 ) base[l][m][pre] = 31000.;
1331              //      //
1332              //      // if we are NOT able to calculate the baseline with more than 3 strips on both pres take the baseline (if any) of the one which has less energy
1333              //      //
1334              //      if ( nstt[0] < 4 && nstt[1] < 4 ){
1335              //        if ( rqp[0] >= rqp[1] ) base[l][m][pre-1] = 31000.;
1336              //        if ( rqp[0] < rqp[1] ) base[l][m][pre] = 31000.;
1337              //      };
1338          };          };
1339          //          //
1340          // run over strips          // run over strips
# Line 512  Int_t CaloLevel0::Calibrate(Int_t ei){ Line 1345  Int_t CaloLevel0::Calibrate(Int_t ei){
1345            ip[i] = 0;            ip[i] = 0;
1346            for (Int_t n = i*32 ; n < (i+1)*32 ; n++){                            for (Int_t n = i*32 ; n < (i+1)*32 ; n++){                
1347              if (n%16 == 0) {              if (n%16 == 0) {
               ck = 0;  
1348                done = 0;                done = 0;
1349                doneb = 0;                doneb = 0;
1350                donec = 0;                donec = 0;
1351                pre++;                pre++;
1352                qpre[pre] = 0.;                ck[l][m][pre] = 0;
1353                  qpre[l][m][pre] = 0.;
1354              };              };
1355              //              //
1356              // baseline check and calculation              // baseline check and calculation
# Line 526  Int_t CaloLevel0::Calibrate(Int_t ei){ Line 1359  Int_t CaloLevel0::Calibrate(Int_t ei){
1359              //              //
1360              if ( !done ){              if ( !done ){
1361                if ( (base[l][m][pre] == 31000. || base[l][m][pre] == 0.) ){                if ( (base[l][m][pre] == 31000. || base[l][m][pre] == 0.) ){
1362                  ck = 1;                  ck[l][m][pre] = 1;
1363                  if (pre%2 == 0) {                  if (pre%2 == 0) {
1364                    ip[i] = pre + 1;                    ip[i] = pre + 1;
1365                  } else {                  } else {
1366                    ip[i] = pre - 1;                    ip[i] = pre - 1;
1367                  };                  };
1368                  if ( (base[l][m][ip[i]] == 31000. || base[l][m][ip[i]] == 0.) ){                  if ( (base[l][m][ip[i]] == 31000. || base[l][m][ip[i]] == 0. || !crosst ) ){
1369                    //                    //
1370                    ck = 2;                    ck[l][m][pre] = 2;
1371                    if ( sbase[l][m][pre] == 31000. || sbase[l][m][pre] == 0. ) {                    if ( sbase[l][m][pre] == 31000. || sbase[l][m][pre] == 0. ) {
1372                      ck = 3;                      ck[l][m][pre] = 3;
1373                    };                    };
1374                  };                  };
                 done = 1;  
1375                };                };
1376                  done = 1;
1377              };                                };                  
1378              //              //
1379              // CALIBRATION ALGORITHM              // CALIBRATION ALGORITHM
1380              //              //
1381              if ( !doneb ){              if ( !doneb ){
1382                switch (ck) {                if ( debug ) printf(" ck[l][m][pre] is %i \n",ck[l][m][pre]);
1383                  switch (ck[l][m][pre]) {
1384                case 0:                case 0:
1385                  base0 = base[l][m][pre];                  base0 = base[l][m][pre];
1386                  base2 = calbase[l][m][pre];                  base2 = calbase[l][m][pre];
1387                    if ( debug ) printf(" base0 = base l%i m%i pre%i = %f base2 = calbase l m pre = %f \n",l,m,pre,base[l][m][pre],calbase[l][m][pre]);
1388                  break;                  break;
1389                case 1:                case 1:
1390                  base0 = base[l][m][ip[i]];                  base0 = base[l][m][ip[i]];
1391                  base2 = calbase[l][m][ip[i]];                  base2 = calbase[l][m][ip[i]];
1392                    if ( debug ) printf(" base0 = base l%i m%i ip(i)%i = %f base2 = calbase l m ip(i) = %f \n",l,m,ip[i],base[l][m][ip[i]],calbase[l][m][ip[i]]);
1393                  break;                  break;
1394                case 2:                case 2:
1395                  base0 = sbase[l][m][pre];                  base0 = sbase[l][m][pre];
1396                  base2 = calbase[l][m][pre];                  base2 = calbase[l][m][pre];    
1397                    if ( debug ) printf(" base0 = sbase l%i m%i pre%i = %f base2 = calbase l m pre = %f \n",l,m,pre,sbase[l][m][pre],calbase[l][m][pre]);
1398                  break;                  break;
1399                case 3:                case 3:
1400                  base0 = calbase[l][m][pre];                  base0 = calbase[l][m][pre];
1401                  base2 = calbase[l][m][pre];                  base2 = calbase[l][m][pre];
1402                    if ( debug ) printf(" base0 = calbase l%i m%i pre%i = %f base2 = calbase l m pre = %f \n",l,m,pre,calbase[l][m][pre],calbase[l][m][pre]);
1403                  break;                  break;
1404                };                };
1405                base1 = calbase[l][m][pre];                base1 = calbase[l][m][pre];
# Line 570  Int_t CaloLevel0::Calibrate(Int_t ei){ Line 1408  Int_t CaloLevel0::Calibrate(Int_t ei){
1408              ener = dexyc[l][m][n];              ener = dexyc[l][m][n];
1409              ener0 += ener;              ener0 += ener;
1410              clevel1->estrip[n][m][l] = 0.;              clevel1->estrip[n][m][l] = 0.;
1411                if ( de->base[l][m][pre] < 0 ) negbase = true;
1412              if ( base0>0 && base0 < 30000. ){              if ( base0>0 && base0 < 30000. ){
1413                if ( !donec && (base0 - base1 + base2) != 0. ){                //
1414                  sbase[l][m][pre] = base0 - base1 + base2;                // save the baseline only if the energy release is "small"
1415                  //
1416                  if ( !donec && (base0 + base1 - base2) != 0. && (n+1)%16==0 ){
1417                    if ( qpre[l][m][pre] < 200. ) sbase[l][m][pre] = base0 + base1 - base2;
1418                  donec = 1;                  donec = 1;
1419                };                };
1420                if ( ener > 0. ){                if ( ener > 0. ){
# Line 580  Int_t CaloLevel0::Calibrate(Int_t ei){ Line 1422  Int_t CaloLevel0::Calibrate(Int_t ei){
1422                  //                  //
1423                  // OK, now in estrip we have the energy deposit in MIP of all the strips for this event (at the end of loops of course)                  // OK, now in estrip we have the energy deposit in MIP of all the strips for this event (at the end of loops of course)
1424                  //                  //
1425                  qpre[pre] += clevel1->estrip[n][m][l];                  if ( clevel1->estrip[n][m][l] > 0. ) qpre[l][m][pre] += clevel1->estrip[n][m][l];
1426                };                  //
1427              };                  //
           };  
           if ( crosst ){  
             if (ck == 1){  
               if (ip[i]%2 == 0) {  
                 ipre = ip[i] + 1;  
               } else {  
                 ipre = ip[i] - 1;  
               };  
               for (Int_t j = ipre*16 ; j < (ipre+1)*16 ; j++){  
                 clevel1->estrip[j][m][l] += (qpre[ipre] - qpre[ip[i]]) * 0.00478;  
               };  
             };  
             if (ck == 2){  
               for (Int_t j = i*32 ; j < (i+1)*32 ; j++){  
                 ipre = j/16 + 1;  
                 clevel1->estrip[j][m][l] +=  qpre[ipre] * 0.00478;  
1428                };                };
1429              };              };
1430            };            };
1431          };          };
1432          //          //
1433          if ( ener0 == 0. && cbase0 == 0. && !pproblem && clevel2->perr[se] == 0){          // check if there were problems with 5.7 or glitches in the power supply
1434            //
1435            if ( ((ener0 == 0. && cbase0 == 0.) || negbase ) && !pproblem && clevel2->perr[se] == 0){
1436            if ( verbose ) printf(" L0 entry %i : calorimeter power problems! event marked as bad perr %f swerr %X view %i plane %i \n",ei,de->perror[se],de->stwerr[se],l,m);            if ( verbose ) printf(" L0 entry %i : calorimeter power problems! event marked as bad perr %f swerr %X view %i plane %i \n",ei,de->perror[se],de->stwerr[se],l,m);
1437            pproblem = true;            pproblem = true;
1438            pe++;            pe++;
1439          };          };
1440          //          //
1441          Int_t j4 = -4;        } else {
1442          Int_t jjj = -3;          for (Int_t nn = 0; nn < 96; nn++){                  
1443          Int_t jj = -2;            clevel1->estrip[nn][m][l] = 0.;
1444          for (Int_t j = 0 ; j < 100 ; j++){          };
1445            jj++;        };
1446            jjj++;      };
1447            j4++;    };
1448            if ( j < 96 ) ene[j] = clevel1->estrip[j][m][l];    //
1449            if ( crosst ){    // run over views and planes to apply crosstalk corrections
1450              if ( jj >= 0 && jj < 96 ){    //
1451                if ( jj != 0 && jj != 32 && jj != 64 ) ene[jj-1] += -clevel1->estrip[jj][m][l] * 0.01581;    for (Int_t l = 0; l < 2; l++){
1452                if ( jj != 31 && jj != 63 && jj != 95 ) ene[jj+1] += -clevel1->estrip[jj][m][l] * 0.01581;                              for (Int_t m = 0; m < 22; m++){
1453              };        //
1454              if ( jjj >= 0 && jjj < 96 ){        // determine the section number  
1455                if ( jjj != 0 && jjj != 32 && jjj != 64 ) clevel1->estrip[jjj-1][m][l] +=  -ene[jjj] * 0.01581;        //
1456                if ( jjj != 31 && jjj != 63 && jjj != 95 ) clevel1->estrip[jjj+1][m][l] +=  -ene[jjj] * 0.01581;        se = 5;
1457              };        if (l == 0 && m%2 == 0) se = 3;
1458            };        if (l == 0 && m%2 != 0) se = 2;
1459            if ( j4 >= 0 && j4 < 96 ){        if (l == 1 && m%2 != 0) se = 1;
1460          if (l == 1 && m%2 == 0) se = 0;          
1461          //
1462          // check for any error in the event
1463          //
1464          if ( clevel2->crc[se] == 0 && (clevel1->good2 == 1 || clevel2->trigty >= 2) ){
1465            //
1466            // Cross-talk corrections  
1467            //
1468            if ( crosst ){
1469              //
1470              // energy on silicon ladders
1471              //
1472              Float_t qsi[3];
1473              qsi[0] = qpre[l][m][0]+qpre[l][m][1];
1474              qsi[1] = qpre[l][m][2]+qpre[l][m][3];
1475              qsi[2] = qpre[l][m][4]+qpre[l][m][5];      
1476              //
1477              for ( pre = 1; pre < 6; pre += 2 ){  
1478                Int_t ladder = (pre - 1)/2;
1479              //              //
1480              // NOTICE: THE FOLLOWING LINE EXCLUDE ALL STRIPS FOR WHICH THE RMS*4 IS GREATER THAN 26 !!! <=============== IMPORTANT! =================>              // If the noselfct flag is set the strip doesn't suffer the self crosstalk due to electronics so we must subtract some energy
1481              //              //
1482              if ( obadmask[l][m][j4] == 1 || clevel1->estrip[j4][m][l] <= clevel1->emin || calrms[l][m][j4] > 26 ){              if ( noselfct ){
1483                clevel1->estrip[j4][m][l] = 0.;                for (Int_t j = ladder*32 ; j < (ladder+1)*32 ; j++){
1484                    ipre = j/16 ;  
1485                    if ( clevel1->estrip[j][m][l] != 0. ) clevel1->estrip[j][m][l] -=  clevel1->estrip[j][m][l] * ctprecor[l][m][ipre];
1486                  };
1487              };              };
1488              //              //
1489              // code and save the energy for each strip in svstrip              // Using the neighbour pre baseline
1490              //              //
1491              if ( clevel1->estrip[j4][m][l] > clevel1->emin ){              if (ck[l][m][pre] == 1 || ck[l][m][pre-1] == 1){
1492                //                //
1493                tim = 100000.;                // pre-amplifier effect on baseline when using the neighbour pre (ck=1)
1494                plo = m;                //
1495                fbi = 0;                if (ck[l][m][pre] == 1){
1496                if ( clevel1->estrip[j4][m][l] > 0.99995 ){                  ipre = pre;
1497                  tim = 10000.;                  ipp = pre - 1;
1498                  plo = m;                } else {
1499                  fbi = 1;                  ipre = pre - 1;
1500                };                  ipp = pre;
               if ( clevel1->estrip[j4][m][l] > 9.9995 ){  
                 tim = 1000.;  
                 plo = 22 + m;  
                 fbi = 1;  
               };  
               if ( clevel1->estrip[j4][m][l] > 99.995 ){  
                 tim = 100.;  
                 plo = 22 + m;  
                 fbi = 0;  
1501                };                };
1502                if ( clevel1->estrip[j4][m][l] > 999.95 ){                Int_t it = 0;
1503                  tim = 10.;                Float_t nqpre = 0.;
1504                  plo = 44 + m;                //
1505                  fbi = 0;                if ( debug ) printf(" CK1 Limit for while: 0.07 \n");
1506                  for (Int_t j = ipre*16 ; j < (ipre+1)*16 ; j++){
1507                    if ( !ctground ){
1508                      if ( clevel1->estrip[j][m][l] != 0. ) clevel1->estrip[j][m][l] += - qpre[l][m][ipp] * ctprecor[l][m][ipp];
1509                    } else {
1510                      if ( clevel1->estrip[j][m][l] != 0. ) clevel1->estrip[j][m][l] += - qpre[l][m][ipp] * 0.00478;
1511                    };
1512                    if ( clevel1->estrip[j][m][l] > 0. ) nqpre += clevel1->estrip[j][m][l] ;
1513                };                };
1514                if ( clevel1->estrip[j4][m][l] > 9999.5 ){                qpre[l][m][ipre] = nqpre;
1515                  tim = 1.;                nqpre = 0.;
1516                  plo = 66 + m;                Float_t deltaqpre = qpre[l][m][ipre];          
1517                  fbi = 0;                //
1518                  // these values are empirically determined, usually the routine converge due to deltaqsi and the latest applied correction is based on less than 1 mip
1519                  //
1520                  while ( it < 10 && deltaqpre > 0.07 ){
1521                    nqpre = 0.;
1522                    for (Int_t j = ipre*16 ; j < (ipre+1)*16 ; j++){
1523                      if ( !ctground ){
1524                        if ( debug ) printf(" CK1 pre correction: iteration %i deltaqpre %f  ctprecor %f  TOTAL CORRECTION %f \n",it,deltaqpre,ctprecor[l][m][ipre],deltaqpre * ctprecor[l][m][ipre]);
1525                        if ( clevel1->estrip[j][m][l] != 0. ) clevel1->estrip[j][m][l] += deltaqpre * ctprecor[l][m][ipre];
1526                      } else {
1527                        if ( clevel1->estrip[j][m][l] != 0. ) clevel1->estrip[j][m][l] += deltaqpre * 0.00478;
1528                      };
1529                      if ( clevel1->estrip[j][m][l] > 0. ) nqpre += clevel1->estrip[j][m][l] ;
1530                    };
1531                    if ( ctground ) it = 100;
1532                    it++;
1533                    deltaqpre = nqpre - qpre[l][m][ipre];
1534                    if ( debug ) printf(" CK1 BEFORE: qpre %f \n",qpre[l][m][ipre]);  
1535                    qpre[l][m][ipre] = nqpre;
1536                    if ( debug ) printf(" CK1 AFTER: qpre %f \n",qpre[l][m][ipre]);  
1537                };                };
1538                //                //
1539                cle = (Int_t)lroundf(tim*clevel1->estrip[j4][m][l]);              };
1540                //
1541                // No baseline calculation due to high energy release
1542                //
1543                if (ck[l][m][pre] == 2 && ck[l][m][pre-1] == 2){
1544                  //
1545                  // y^
1546                  //  |
1547                  //  |    6 7 8
1548                  //  |    3 4 5  
1549                  //  |    0 1 2
1550                  //  | --------------------------------------> x
1551                //                //
1552                  Int_t si1 = 0;
1553                  Int_t si2 = 0;
1554                  Int_t si3 = 0;
1555                if ( l == 0 ){                if ( l == 0 ){
1556                  //                  if ( ladder == 0 ){
1557                  // +-PPSSmmmm.mmmm                    si1 = 0;
1558                  //                    si2 = 3;
1559                  svstrip[istrip] = fbi*1000000000 + plo*10000000 + j4*100000 + cle;                    si3 = 6;
1560                    };
1561                    if ( ladder == 1 ){
1562                      si1 = 1;
1563                      si2 = 4;
1564                      si3 = 7;
1565                    };
1566                    if ( ladder == 2 ){
1567                      si1 = 2;
1568                      si2 = 5;
1569                      si3 = 8;
1570                    };
1571                } else {                } else {
1572                  svstrip[istrip] = -(fbi*1000000000 + plo*10000000 + j4*100000 + cle);                  if ( ladder == 0 ){
1573                      si1 = 0;
1574                      si2 = 1;
1575                      si3 = 2;
1576                    };
1577                    if ( ladder == 1 ){
1578                      si1 = 3;
1579                      si2 = 4;
1580                      si3 = 5;
1581                    };
1582                    if ( ladder == 2 ){
1583                      si1 = 6;
1584                      si2 = 7;
1585                      si3 = 8;
1586                    };
1587                };                };
1588                //                //
1589                //              if ( ei >= -770 ) printf(" j %i l %i m %i estrip %f \n",j4,l,m,clevel1->estrip[j4][m][l]);                // Find the energy distribution along the considered plane looking at the two sandwiching plane of the other view.
               //              if ( ei >= -770 ) printf(" num lim %i fbi %i tim %f plo %i cle %i \n",numeric_limits<Int_t>::max(),fbi,tim,plo,cle);  
               //              if ( ei >= -770 ) printf(" svstrip %i \n",svstrip[istrip]);  
1590                //                //
1591                istrip++;                Float_t sied[3] = {0.,0.,0.};
1592                  Int_t othv = !l;
1593                  Int_t othpl1 = m - 1;
1594                  Int_t othpl2 = m + 1;
1595                  Float_t oprof[3] = {0.,0.,0.};
1596                  for(Int_t s=0; s<3; s++){
1597                    for(Int_t t=(s*32); t<32*(s + 1); t++){      
1598                      if ( othpl1 > -1 ) {
1599                        oprof[s] += clevel1->estrip[othv][othpl1][t];
1600                      };
1601                      if ( othpl2 < 22 ) {
1602                        oprof[s] += clevel1->estrip[othv][othpl2][t];
1603                      };          
1604                    };
1605                  };
1606                  Float_t otote = fabs(oprof[0]) + fabs(oprof[1]) + fabs(oprof[2]);
1607                  for(Int_t g=0; g<3; g++){
1608                    if ( otote > 0. ){
1609                      sied[g] = fabs(oprof[g])/otote;
1610                    } else {
1611                      sied[g] = 1./3.;
1612                    };
1613                  };
1614                  //
1615                  //
1616                  //
1617                  Int_t it = 0;
1618                  Int_t jpre = 0;
1619                  Float_t nqsi = 0.;
1620                  Float_t snqsi = qsi[ladder];
1621                  Float_t nqpre[2] = {0.,0.};
1622                  Float_t deltaqsi = qsi[ladder];
1623                  Float_t deltaqpre[2];
1624                  deltaqpre[0] = qpre[l][m][pre-1];
1625                  deltaqpre[1] = qpre[l][m][pre];
1626                  //
1627                  if ( debug ) printf(" Limit for while: 0.07 it < 10 \n");
1628                  //
1629                  // these values are empirically determined, usually the routine converge due to deltaqsi and the latest applied correction is based on less than 1 mip
1630                  //
1631                  while ( it < 10 && (deltaqsi > 0.07 || deltaqpre[0] > 0.07 || deltaqpre[1] > 0.07) ){
1632                    nqsi = 0.;
1633                    nqpre[0] = 0.;
1634                    nqpre[1] = 0.;
1635                    for (Int_t j = ladder*32 ; j < (ladder+1)*32 ; j++){
1636                      ipre = 0;
1637                      if ( j > (ladder*32)+15 ) ipre = 1;
1638                      jpre = j/16 ;
1639                      //
1640                      // Silicon effect on the baseline when using the same pre previous baseline (ck = 2) + pre-amply effect
1641                      //          
1642                      if ( !ctground ){
1643                        if ( debug ) printf(" silicon correction: iteration %i deltaqsi[%i] %f ctsicor %f %f %f sied %f %f %f si %i %i %i TOTAL CORRECTION %f \n",it,ladder,deltaqsi,ctsicor[l][m][si1],ctsicor[l][m][si2],ctsicor[l][m][si3],sied[0],sied[1],sied[2],si1,si2,si3,deltaqsi * (ctsicor[l][m][si1] * sied[0] + ctsicor[l][m][si2] * sied[1] + ctsicor[l][m][si3] * sied[2]));
1644                        if ( debug ) printf(" pre correction: iteration %i deltaqpre[0] %f deltaqpre[1] %f ctprecor %f  TOTAL CORRECTION %f \n",it,deltaqpre[0],deltaqpre[1],ctprecor[l][m][jpre],deltaqpre[ipre] * ctprecor[l][m][jpre]);
1645                        if ( clevel1->estrip[j][m][l] != 0. ) clevel1->estrip[j][m][l] +=  (deltaqsi * (ctsicor[l][m][si1] * sied[0] + ctsicor[l][m][si2] * sied[1] + ctsicor[l][m][si3] * sied[2])/mip[l][m][j]) + deltaqpre[ipre] * ctprecor[l][m][jpre];
1646                      } else {
1647                        if ( clevel1->estrip[j][m][l] != 0. ) clevel1->estrip[j][m][l] +=  0. + qpre[l][m][jpre] * 0.00478; // no correction
1648                      };
1649                      if ( clevel1->estrip[j][m][l] > 0. ) nqsi += clevel1->estrip[j][m][l] ;
1650                      if ( clevel1->estrip[j][m][l] > 0. ) nqpre[ipre] += clevel1->estrip[j][m][l] ;
1651                    };      
1652                    if ( ctground ) it = 100;
1653                    deltaqsi = nqsi-snqsi;
1654                    snqsi = nqsi;
1655                    it++;
1656                    deltaqpre[0] = nqpre[0] - qpre[l][m][pre-1];
1657                    deltaqpre[1] = nqpre[1] - qpre[l][m][pre];
1658                    if ( debug ) printf(" BEFORE: qpre 0 %f qpre 1 %f \n",qpre[l][m][pre-1],qpre[l][m][pre]);  
1659                    qpre[l][m][pre-1] = nqpre[0];
1660                    qpre[l][m][pre] = nqpre[1];    
1661                    if ( debug ) printf(" AFTER: qpre 0 %f qpre 1 %f \n",qpre[l][m][pre-1],qpre[l][m][pre]);  
1662                  };
1663                  //
1664                  //
1665                  //
1666    //            for (Int_t j = ladder*32 ; j < (ladder+1)*32 ; j++){
1667    //              ipre = j/16 ;          
1668    //              //
1669    //              // pre-amplifier effect on baseline when using the same pre previous event baseline (ck=2)
1670    //              //
1671    //              if ( !ctground ){
1672    //                if ( clevel1->estrip[j][m][l] != 0. ) clevel1->estrip[j][m][l] +=  qpre[l][m][ipre] * ctprecor[l][m][ipre];
1673    //              } else {
1674    //                if ( clevel1->estrip[j][m][l] != 0. ) clevel1->estrip[j][m][l] +=  qpre[l][m][ipre] * 0.00478;
1675    //              };
1676    //            };
1677              };              };
1678            };            };
1679          };          };
1680          //        };
1681        } else {        //
1682          for (Int_t nn = 0; nn < 96; nn++){                          Int_t j4 = -4;
1683            clevel1->estrip[nn][m][l] = 0.;        Int_t jjj = -3;
1684          Int_t jj = -2;
1685          Int_t jjpre = -1;
1686          Int_t jjjpre = -1;
1687          memset(ene, 0, 96*sizeof(Float_t));
1688          for (Int_t j = 0 ; j < 100 ; j++){          
1689            jj++;
1690            jjj++;
1691            j4++;
1692            if ( j < 96 ) ene[j] = clevel1->estrip[j][m][l];
1693            if ( crosst ){
1694              //
1695              // "Real" crosstalk effect on the neighbour strips respect to the one which have seen the energy deposit
1696              //
1697              if ( jj >= 0 && jj < 96 ){
1698                if ( !ctground ){
1699                  if ( jj%16 == 0 ) jjpre++;              
1700                  if ( jj != 0 && jj != 32 && jj != 64 && ene[jj-1] != 0. ) ene[jj-1] += -clevel1->estrip[jj][m][l] * ctneigcor[l][m][jjpre];
1701                  if ( jj != 31 && jj != 63 && jj != 95 && ene[jj+1] != 0. ) ene[jj+1] += -clevel1->estrip[jj][m][l] * ctneigcor[l][m][jjpre];                      
1702                } else {
1703                  if ( jj != 0 && jj != 32 && jj != 64 && ene[jj-1] != 0. ) ene[jj-1] += -clevel1->estrip[jj][m][l] * 0.01581;
1704                  if ( jj != 31 && jj != 63 && jj != 95 && ene[jj+1] != 0. ) ene[jj+1] += -clevel1->estrip[jj][m][l] * 0.01581;                    
1705                };
1706              };
1707              if ( jjj >= 0 && jjj < 96 ){
1708                if ( !ctground ){
1709                  if ( jjj%16 == 0 ) jjjpre++;            
1710                  if ( jjj != 0 && jjj != 32 && jjj != 64 && clevel1->estrip[jjj-1][m][l] != 0. ) clevel1->estrip[jjj-1][m][l] +=  -ene[jjj] * ctneigcor[l][m][jjjpre];
1711                  if ( jjj != 31 && jjj != 63 && jjj != 95 && clevel1->estrip[jjj+1][m][l] !=0. ) clevel1->estrip[jjj+1][m][l] +=  -ene[jjj] * ctneigcor[l][m][jjjpre];
1712                } else {
1713                  if ( jjj != 0 && jjj != 32 && jjj != 64 && clevel1->estrip[jjj-1][m][l] != 0. ) clevel1->estrip[jjj-1][m][l] +=  -ene[jjj] * 0.01581;
1714                  if ( jjj != 31 && jjj != 63 && jjj != 95 && clevel1->estrip[jjj+1][m][l] != 0. ) clevel1->estrip[jjj+1][m][l] +=  -ene[jjj] * 0.01581;
1715                };
1716              };
1717            };
1718            if ( j4 >= 0 && j4 < 96 ){
1719              //
1720              // CALOLEVEL1 CODING AND FILLING
1721              //
1722              //
1723              // NOTICE: THE FOLLOWING LINE EXCLUDE ALL STRIPS FOR WHICH THE RMS*4 IS GREATER THAN 26 !!! <=============== IMPORTANT! =================>
1724              //
1725              if ( obadmask[l][m][j4] == 1 || clevel1->estrip[j4][m][l] <= clevel1->emin || clevel1->estrip[j4][m][l] <= memin[l][m][j4] || calrms[l][m][j4] > maxrms[l][m] ){
1726                clevel1->estrip[j4][m][l] = 0.;
1727              };
1728              //
1729              // code and save the energy for each strip in svstrip    
1730              //
1731              if ( clevel1->estrip[j4][m][l] > clevel1->emin ){
1732                //
1733                Float_t savel1 = clevel1->estrip[j4][m][l];
1734                //      if ( dexyc[l][m][j4] == 32767. ){
1735                if ( dexyc[l][m][j4] > 32000. ){
1736                  savel1 += 5000.;
1737                  clevel2->nsatstrip += 1.;
1738                };
1739                //
1740                tim = 100000.;
1741                plo = m;
1742                fbi = 0;
1743                if ( savel1 > 0.99995 ){
1744                  tim = 10000.;
1745                  plo = m;
1746                  fbi = 1;
1747                };
1748                if ( savel1 > 9.9995 ){
1749                  tim = 1000.;
1750                  plo = 22 + m;
1751                  fbi = 1;
1752                };
1753                if ( savel1 > 99.995 ){
1754                  tim = 100.;
1755                  plo = 22 + m;
1756                  fbi = 0;
1757                };
1758                if ( savel1 > 999.95 ){
1759                  tim = 10.;
1760                  plo = 44 + m;
1761                  fbi = 0;
1762                };
1763                if ( savel1 > 9999.5 ){
1764                  tim = 1.;
1765                  plo = 66 + m;
1766                  fbi = 0;
1767                };
1768                //
1769                cle = (Int_t)lroundf(tim*savel1);
1770                //
1771                if ( l == 0 ){
1772                  //
1773                  // +-PPSSmmmm.mmmm
1774                  //
1775                  svstrip[istrip] = fbi*1000000000 + plo*10000000 + j4*100000 + cle;
1776                } else {
1777                  svstrip[istrip] = -(fbi*1000000000 + plo*10000000 + j4*100000 + cle);
1778                };
1779                //
1780                istrip++;
1781              };
1782          };          };
1783        };        };
1784          //
1785      };      };
1786    };    };
1787      //
1788      // store goodness flag
1789      //
1790    if ( !pe  ){    if ( !pe  ){
1791      clevel2->good = 1;      clevel2->good = 1;
1792    } else {    } else {
1793      clevel2->good = 0;      clevel2->good = 0;
1794    };    };
1795      //
1796      // done
1797      //
1798    return(0);    return(0);
1799  }  }
1800    
# Line 749  void CaloLevel0::FillTrkVar(CaloLevel2 * Line 1843  void CaloLevel0::FillTrkVar(CaloLevel2 *
1843      memcpy(t_ca->tibar,clevel2->cibar,sizeof(clevel2->cibar));      memcpy(t_ca->tibar,clevel2->cibar,sizeof(clevel2->cibar));
1844      memcpy(t_ca->tbar,clevel2->cbar,sizeof(clevel2->cbar));      memcpy(t_ca->tbar,clevel2->cbar,sizeof(clevel2->cbar));
1845      memcpy(ca->planemax,clevel2->planemax,sizeof(clevel2->planemax));      memcpy(ca->planemax,clevel2->planemax,sizeof(clevel2->planemax));
1846        memcpy(ca->selfdelay,clevel2->selfdelay,sizeof(clevel2->selfdelay));
1847      ca->varcfit[2] = clevel2->varcfit[0];      ca->varcfit[2] = clevel2->varcfit[0];
1848      ca->varcfit[3] = clevel2->varcfit[1];      ca->varcfit[3] = clevel2->varcfit[1];
1849      ca->npcfit[2] = clevel2->npcfit[0];      ca->npcfit[2] = clevel2->npcfit[0];
# Line 789  void CaloLevel0::FillCommonVar(CaloLevel Line 1884  void CaloLevel0::FillCommonVar(CaloLevel
1884    memcpy(ca->swerr,clevel2->swerr,sizeof(clevel2->swerr));    memcpy(ca->swerr,clevel2->swerr,sizeof(clevel2->swerr));
1885    memcpy(ca->crc,clevel2->crc,sizeof(clevel2->crc));    memcpy(ca->crc,clevel2->crc,sizeof(clevel2->crc));
1886    ca->nstrip = (Int_t)clevel2->nstrip;    ca->nstrip = (Int_t)clevel2->nstrip;
1887      ca->nsatstrip = (Int_t)clevel2->nsatstrip;
1888    ca->qtot = clevel2->qtot;    ca->qtot = clevel2->qtot;
1889    //  ca->impx = clevel2->impx;    //  ca->impx = clevel2->impx;
1890    //  ca->impy = clevel2->impy;    //  ca->impy = clevel2->impy;
# Line 801  void CaloLevel0::FillCommonVar(CaloLevel Line 1897  void CaloLevel0::FillCommonVar(CaloLevel
1897    ca->selen = clevel2->selen;    ca->selen = clevel2->selen;
1898    memcpy(ca->qq,clevel2->qq,sizeof(clevel2->qq));    memcpy(ca->qq,clevel2->qq,sizeof(clevel2->qq));
1899    memcpy(ca->planemax,clevel2->planemax,sizeof(clevel2->planemax));    memcpy(ca->planemax,clevel2->planemax,sizeof(clevel2->planemax));
1900      memcpy(ca->selfdelay,clevel2->selfdelay,sizeof(clevel2->selfdelay));
1901    ca->varcfit[0] = clevel2->varcfit[0];    ca->varcfit[0] = clevel2->varcfit[0];
1902    ca->varcfit[1] = clevel2->varcfit[1];    ca->varcfit[1] = clevel2->varcfit[1];
1903    ca->npcfit[0] = clevel2->npcfit[0];    ca->npcfit[0] = clevel2->npcfit[0];
# Line 824  void CaloLevel0::ClearStructs(){ Line 1921  void CaloLevel0::ClearStructs(){
1921    ClearCommonVar();    ClearCommonVar();
1922  }  }
1923    
1924    void CaloLevel0::Delete(Option_t *t){
1925      if ( de ) delete de;
1926      delete this;
1927    }
1928    
1929    
1930  void CaloLevel0::RunClose(){  void CaloLevel0::RunClose(){
1931    l0tr->Delete();    l0tr->Delete();
1932    ClearStructs();    ClearStructs();
# Line 832  void CaloLevel0::RunClose(){ Line 1935  void CaloLevel0::RunClose(){
1935    memset(dexyc, 0, 2*22*96*sizeof(Float_t));    memset(dexyc, 0, 2*22*96*sizeof(Float_t));
1936    memset(base, 0, 2*22*6*sizeof(Float_t));    memset(base, 0, 2*22*6*sizeof(Float_t));
1937    memset(sbase, 0, 2*22*6*sizeof(Float_t));    memset(sbase, 0, 2*22*6*sizeof(Float_t));
1938      memset(ctprecor, 0, 2*22*6*sizeof(Float_t));
1939      memset(ctsicor, 0, 2*22*9*sizeof(Float_t));
1940      memset(ctneigcor, 0, 2*22*6*sizeof(Float_t));
1941    //    //
1942  }  }
1943    
# Line 876  void CaloLevel0::ClearCommonVar(){ Line 1982  void CaloLevel0::ClearCommonVar(){
1982    clevel2->wartrig = 0.;    clevel2->wartrig = 0.;
1983    clevel2->good = 0;    clevel2->good = 0;
1984    clevel2->nstrip = 0.;    clevel2->nstrip = 0.;
1985      clevel2->nsatstrip = 0.;
1986    clevel2->qtot = 0.;    clevel2->qtot = 0.;
1987    //  clevel2->impx = 0.;    //  clevel2->impx = 0.;
1988    //  clevel2->impy = 0.;    //  clevel2->impy = 0.;
# Line 891  void CaloLevel0::ClearCommonVar(){ Line 1998  void CaloLevel0::ClearCommonVar(){
1998    memset(clevel2->varcfit, 0, 4*sizeof(Float_t));    memset(clevel2->varcfit, 0, 4*sizeof(Float_t));
1999    memset(clevel2->npcfit, 0, 4*sizeof(Int_t));    memset(clevel2->npcfit, 0, 4*sizeof(Int_t));
2000    memset(clevel2->planemax, 0, 2*sizeof(Int_t));    memset(clevel2->planemax, 0, 2*sizeof(Int_t));
2001      memset(clevel2->selfdelay, 0, 4*7*sizeof(Int_t));
2002    memset(clevel2->fmode, 0, 2*sizeof(Int_t));    memset(clevel2->fmode, 0, 2*sizeof(Int_t));
2003    memset(clevel2->cibar, 0, 2*22*sizeof(Int_t));    memset(clevel2->cibar, 0, 2*22*sizeof(Int_t));
2004    memset(clevel2->cbar, 0, 2*22*sizeof(Float_t));    memset(clevel2->cbar, 0, 2*22*sizeof(Float_t));
# Line 947  void CaloLevel0::ClearCalibVals(Int_t s) Line 2055  void CaloLevel0::ClearCalibVals(Int_t s)
2055    return;    return;
2056  }  }
2057    
2058  Int_t CaloLevel0::Update(TSQLServer *dbc, UInt_t atime, Int_t s){  Int_t CaloLevel0::Update(GL_TABLES *glt, UInt_t atime, Int_t s){
2059    //    //
2060      const TString host = glt->CGetHost();
2061      const TString user = glt->CGetUser();
2062      const TString psw = glt->CGetPsw();
2063      TSQLServer *dbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data());
2064      if ( !dbc->IsConnected() ) throw -116;
2065      stringstream myquery;
2066      myquery.str("");
2067      myquery << "SET time_zone='+0:00'";
2068      dbc->Query(myquery.str().c_str());
2069    Int_t sgnl = 0;    Int_t sgnl = 0;
2070    //    //
2071    GL_CALO_CALIB *glcalo = new GL_CALO_CALIB();    GL_CALO_CALIB *glcalo = new GL_CALO_CALIB();
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.23
  ViewVC Help
Powered by ViewVC 1.1.23