/[PAMELA software]/calo/flight/CaloFranzini/src/CaloFranzini.cpp
ViewVC logotype

Diff of /calo/flight/CaloFranzini/src/CaloFranzini.cpp

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

revision 1.5 by mocchiut, Thu Jan 3 10:02:28 2008 UTC revision 1.10 by mocchiut, Tue Aug 4 13:59:15 2009 UTC
# Line 16  CaloFranzini::CaloFranzini(){ Line 16  CaloFranzini::CaloFranzini(){
16    
17  CaloFranzini::CaloFranzini(PamLevel2 *l2p){    CaloFranzini::CaloFranzini(PamLevel2 *l2p){  
18    //    //
19    file = NULL;    lfile = NULL;
20      ffile = NULL;
21    brig = NULL;    brig = NULL;
22    brigm = NULL;    brigm = NULL;
23    nbin = 0;    nbin = 0;
# Line 52  void CaloFranzini::Clear(){ Line 53  void CaloFranzini::Clear(){
53    longtzeta = 0.;    longtzeta = 0.;
54    fulltzeta = 0.;    fulltzeta = 0.;
55    degfre = 0;    degfre = 0;
56      fdegfre = 0;
57    memset(estrip, 0, 2*22*96*sizeof(Float_t));    memset(estrip, 0, 2*22*96*sizeof(Float_t));
58    memset(qplane, 0, 43*sizeof(Float_t));    memset(qplane, 0, 43*sizeof(Float_t));
59    //    //
60      numneg = 0;
61      numpos = 0;
62      negfulltzeta = 0.;
63      posfulltzeta = 0.;
64      aveposvar = 0.;
65      avenegvar = 0.;
66      minsvalue =  numeric_limits<UInt_t>::max();
67      maxsvalue =  numeric_limits<UInt_t>::min();
68      //
69  }  }
70    
71  void CaloFranzini::Print(){  void CaloFranzini::Print(){
# Line 64  void CaloFranzini::Print(){ Line 75  void CaloFranzini::Print(){
75    printf("========================================================================\n");    printf("========================================================================\n");
76    printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);    printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);
77    printf(" debug                [debug flag]:.. %i\n",debug);    printf(" debug                [debug flag]:.. %i\n",debug);
78    printf(" degree of freedom                :.. %i\n",this->GetDegreeOfFreedom());      printf(" long  degree of freedom          :.. %i\n",this->GetLongDegreeOfFreedom());  
79    printf(" longtzeta                        :.. %f\n",longtzeta);      printf(" longtzeta                        :.. %f\n",longtzeta);  
80    printf(" longtzeta normalized             :.. %f\n",this->GetNormLongTZeta());      printf(" longtzeta normalized             :.. %f\n",this->GetNormLongTZeta());  
81    //  printf(" fulltzeta                        :.. %f\n",fulltzeta);      printf(" full degree of freedom           :.. %i\n",this->GetFullDegreeOfFreedom());  
82    //  printf(" longtzeta normalized             :.. %f\n",this->GetNormFullTZeta());      printf(" fulltzeta                        :.. %f\n",fulltzeta);  
83      printf(" fulltzeta normalized             :.. %f\n",this->GetNormFullTZeta());  
84      printf(" fulltzeta negative contribute    :.. %f\n",negfulltzeta);  
85      printf(" fulltzeta positive contribute    :.. %f\n",posfulltzeta);  
86      printf(" fulltzeta number of negatives    :.. %i\n",numneg);  
87      printf(" fulltzeta number of positives    :.. %i\n",numpos);  
88      printf(" fulltzeta minimum variance       :.. %f\n",minsvalue);  
89      printf(" fulltzeta maximum variance       :.. %f\n",maxsvalue);  
90      printf(" fulltzeta average positive var   :.. %f\n",aveposvar);  
91      printf(" fulltzeta average negative var   :.. %f\n",avenegvar);  
92    printf("========================================================================\n");    printf("========================================================================\n");
93    //    //
94  }  }
95    
96  void CaloFranzini::Delete(){  void CaloFranzini::Delete(){
97    //    //
98    if ( file ) file->Close();    if ( ffile ) ffile->Close();
99      if ( lfile ) lfile->Close();
100    //    //
101    Clear();    Clear();
102    //    //
# Line 108  void CaloFranzini::SetNoWcalo(Int_t n){ Line 129  void CaloFranzini::SetNoWcalo(Int_t n){
129    };    };
130  }  }
131    
132    void CaloFranzini::SplitInto(Int_t NoWpreSampler, Int_t NoWcalo){
133      this->SetNoWpreSampler(0);
134      this->SetNoWcalo(0);
135      if ( NoWpreSampler < NoWcalo ){
136              this->SetNoWpreSampler(NoWpreSampler);
137              this->SetNoWcalo(NoWcalo);
138      } else {
139              this->SetNoWcalo(NoWcalo);            
140              this->SetNoWpreSampler(NoWpreSampler);
141      };
142    }
143    
144    
145  void CaloFranzini::Process(){  void CaloFranzini::Process(){
146    this->Process(0);    this->Process(0);
147  }  }
# Line 120  void CaloFranzini::Process(Int_t itr){ Line 154  void CaloFranzini::Process(Int_t itr){
154      return;      return;
155    };    };
156    //    //
157    if ( !nbin || !file || !brig ){    if ( !nbin || !brig || (!lfile && !ffile) ){
158      printf(" ERROR: it seems covariance matrix file has not been opened (CaloFranzini::Open()) \n");      printf(" ERROR: it seems covariance matrix file has not been opened (CaloFranzini::Open()) \n");
159      printf(" ERROR: CaloFranzini variables _NOT_ filled \n");      printf(" ERROR: CaloFranzini variables _NOT_ filled \n");
160      return;      return;
# Line 158  void CaloFranzini::Process(Int_t itr){ Line 192  void CaloFranzini::Process(Int_t itr){
192    Int_t plane = 0;    Int_t plane = 0;
193    Int_t strip = 0;    Int_t strip = 0;
194    Float_t mip = 0.;    Float_t mip = 0.;
   for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){  
     //  
     mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);  
     //  
     estrip[view][plane][strip] = mip;  
     //      
     nplane = 1 - view + 2 * (plane - N);  
     if ( nplane > (37-(2*N)) ) nplane--;  
     //  
 //    if ( plane == (18+N) ) mip = 0.;  
     if ( nplane > -1 ) qplane[nplane] += mip;  
     //  
   };  
195    //    //
196    //    //
197    //    //
198    if ( dolong ){    if ( dolong ){
199      //      //
200        for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
201          //
202          mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
203          //
204          estrip[view][plane][strip] = mip;
205          //    
206          nplane = 1 - view + 2 * (plane - N);
207          if ( nplane > (37-(2*N)) ) nplane--;
208          //
209          //    if ( plane == (18+N) ) mip = 0.;
210          if ( nplane > -1 ) qplane[nplane] += mip;
211          //
212        };
213        //
214      if ( cont ){      if ( cont ){
215        for (Int_t i=0; i<22; i++){        for (Int_t i=0; i<22; i++){
216          if ( i == (18+N) ){          if ( i == (18+N) ){
# Line 213  void CaloFranzini::Process(Int_t itr){ Line 248  void CaloFranzini::Process(Int_t itr){
248      //      //
249      degfre = TMath::Min(dgf,NC);      degfre = TMath::Min(dgf,NC);
250      //      //
251      Float_t longzdiag = 0.;      //    Float_t longzdiag = 0.;
252      Float_t longzout = 0.;      //    Float_t longzout = 0.;
253      //      //
254      if ( degfre > 0 ){      if ( degfre > 0 ){
255        for (Int_t i = 0; i < degfre; i++){        for (Int_t i = 0; i < degfre; i++){
256          if ( i != mask18b ){          if ( i != mask18b ){
257            for (Int_t j = 0; j < degfre; j++){              for (Int_t j = 0; j < degfre; j++){  
258              if ( j != mask18b ){              if ( j != mask18b ){
259                if ( i == j ){  //            if ( i == j ){
260                  longzdiag += (qplane[i]-this->GetAverageAt(i,rig)) * this->GetHmatrixAt(i,j,rig) * (qplane[j]-this->GetAverageAt(j,rig));  //              longzdiag += (qplane[i]-this->GetAverageAt(i,rig)) * this->GetHmatrixAt(i,j,rig) * (qplane[j]-this->GetAverageAt(j,rig));
261                  if ( debug ) printf(" %i %i %f %f %f %f %f\n",i,j,qplane[i],this->GetAverageAt(i,rig),this->GetHmatrixAt(i,j,rig),qplane[j],this->GetAverageAt(j,rig));  //              if ( debug ) printf(" %i %i %f %f %f %f %f\n",i,j,qplane[i],this->GetAverageAt(i,rig),this->GetHmatrixAt(i,j,rig),qplane[j],this->GetAverageAt(j,rig));
262                } else {  //            } else {
263                  longzout += (qplane[i]-this->GetAverageAt(i,rig)) * this->GetHmatrixAt(i,j,rig) * (qplane[j]-this->GetAverageAt(j,rig));  //              longzout += (qplane[i]-this->GetAverageAt(i,rig)) * this->GetHmatrixAt(i,j,rig) * (qplane[j]-this->GetAverageAt(j,rig));
264                  if ( debug && i == (j+1) ) printf(" %i %i %f %f %f %f %f\n",i,j,qplane[i],this->GetAverageAt(i,rig),this->GetHmatrixAt(i,j,rig),qplane[j],this->GetAverageAt(j,rig));  //              if ( debug && i == (j+1) ) printf(" %i %i %f %f %f %f %f\n",i,j,qplane[i],this->GetAverageAt(i,rig),this->GetHmatrixAt(i,j,rig),qplane[j],this->GetAverageAt(j,rig));
265                };  //            };
266                  //
267                longtzeta += (qplane[i]-this->GetAverageAt(i,rig)) * this->GetHmatrixAt(i,j,rig) * (qplane[j]-this->GetAverageAt(j,rig));                longtzeta += (qplane[i]-this->GetAverageAt(i,rig)) * this->GetHmatrixAt(i,j,rig) * (qplane[j]-this->GetAverageAt(j,rig));
268                  //
269              };              };
270            };            };
271          };          };
272        };        };
273        if ( debug ) printf(" diagonal total %f out of diagonal total %f \n",longzdiag,longzout);        //      if ( debug ) printf(" diagonal total %f out of diagonal total %f \n",longzdiag,longzout);
274      };      };
275    };    };
276    if ( dofull ){    if ( dofull ){
277      printf(" ERROR: NOT IMPLEMENTED YET\n");      //    printf(" ERROR: NOT IMPLEMENTED YET\n");
278      printf(" ERROR: CaloFranzini variables _NOT_ filled \n");      //    printf(" ERROR: CaloFranzini variables _NOT_ filled \n");
279    };        //
280        // FULL CALORIMETER
281        //
282        CaloTrkVar *ct = L2->GetTrack(0)->GetCaloTrack();
283        //
284        Int_t dgf = 0;
285        Int_t cs = 0;
286        Int_t cd = 0;
287        Int_t mstrip = 0;
288        //
289        Int_t maxnpl = 0;
290        Float_t mipv[43][11];
291        memset(mipv,0,43*11*sizeof(Float_t));
292        //
293        for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
294          //
295          mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
296          //
297          //      nplane = 1 - view + 2 * plane;
298          //      if ( nplane > 37 ) nplane--;
299          nplane = 1 - view + 2 * (plane - N);
300          if ( nplane > (37-(2*N)) ) nplane--;
301          //
302          cs = ct->tibar[plane][view] - 1;
303          //
304          if ( ct->tibar[plane][view] > -1 ){
305            //
306            cd = 95 - cs;
307            //
308            mstrip = cd + strip;
309            //
310            Int_t lstr = this->ConvertStrip(mstrip);
311            //
312            if ( nplane > maxnpl ) maxnpl = nplane;
313            if ( nplane > -1 ) mipv[nplane][lstr] += mip;
314            //
315          };
316          //      mipv[nplane][lstr] += mip;
317          //
318        };
319        //
320        Float_t mip1 = 1.;
321        Int_t cs1;
322        Int_t cd1;
323        Float_t mip2 = 1.;
324        Int_t cs2;
325        Int_t cd2;
326        Int_t mi = -1;
327        Int_t mj = -1;
328        Int_t nn1 = 0;
329        Int_t pl1 = 0;
330        Int_t vi1 = 0;
331        Int_t nn2 = 0;
332        Int_t pl2 = 0;
333        Int_t vi2 = 0;
334        Int_t mstrip1min = 0;
335        Int_t mstrip1max = 0;
336        Int_t mstrip2min = 0;
337        Int_t mstrip2max = 0;
338        //
339        Int_t tonpl = TMath::Min(maxnpl,NC);
340        //
341        Int_t rbi = 0;
342        for (Int_t i = 0; i<nbin-1; i++){
343          if ( rig>=brigm->At(i) && rig < brigm->At(i+1) ){
344            rbi = i;
345            break;
346          };
347        };
348        //
349        //
350        Int_t therigb = rbi;
351        //
352        if ( rig < brigm->At(2) ){
353    //      therigb = 0;
354          therigb = 2;
355        };
356        //    if ( rig >= brigm->At(nbin-5) ){
357        if ( rig >= brigm->At(nbin-2) ){
358          therigb = nbin - 3;
359          //      therigb = nbin - 5;
360        };
361        Int_t mtherig = therigb;
362        if ( mtherig >= 13 ) mtherig = 12;
363        //
364        if ( debug ) printf(" rig %f brigm0 %f brigm n2 %f nbin %i \n",rig,brigm->At(0),brigm->At(nbin-2),nbin);
365        //
366        if ( cont ){
367          for (Int_t i=0; i<22; i++){
368            if ( i == (18+N) ){
369              mask18b =  1 + 2 * (i - N);
370              break;
371            };
372          };
373        };  
374        //  
375        if ( sel ){
376          for (Int_t i=0; i<22; i++){
377            if ( i == (18-N) ){
378              mask18b =  1 + 2 * (i - N);
379              break;
380            };
381          };
382        };
383        //
384        if ( mask18b == 37 ) mask18b = -1;
385        //
386        if ( debug ) printf("Start main loop \n");
387        //
388        for (Int_t nplane1=0; nplane1<tonpl; nplane1++){
389          //
390          if ( nplane1 != mask18b ){
391            //
392            if ( nplane1 >= 37 ) nn1 = nplane1 + 1;
393            vi1 = 1;
394            if ( nn1%2 ) vi1 = 0;
395            pl1 = (nn1 - 1 + vi1)/2;
396            //
397            cs1 = ct->tibar[pl1][vi1] - 1;
398            //
399            if ( ct->tibar[pl1][vi1] > -1 ){
400              //
401              dgf++;
402              //
403              cd1 = 95 - cs1;
404              //
405              Int_t at1 = TMath::Max(0,(cd1+0));
406              Int_t at2 = TMath::Min(190,(cd1+95));
407              mstrip1min = this->ConvertStrip(at1);
408              mstrip1max = this->ConvertStrip(at2) + 1;
409              //
410              for (Int_t mstrip1=mstrip1min; mstrip1<mstrip1max; mstrip1++){
411                //
412                mj = -1;
413                //
414                //      if ( debug ) printf(" get mip1 %i %i %f %i \n",nplane1,mstrip1,rig,therigb);
415                mip1 = mipv[nplane1][mstrip1] - this->GetFullAverageAt(nplane1,mstrip1,rig,therigb);
416                if ( mip1 < 0 ){
417                  numneg++;
418                  avenegvar += mip1;
419                };
420                if ( mip1 >= 0 ){
421                  numpos++;
422                  aveposvar += mip1;
423                };
424                if ( mip1 < minsvalue ) minsvalue = mip1;
425                if ( mip1 >= maxsvalue ) maxsvalue = mip1;
426                //
427                mi = (nplane1 * 11) + mstrip1;
428                //
429                if ( mip1 != 0. ){
430                  //
431                  for (Int_t nplane2=0; nplane2<tonpl; nplane2++){
432                    //
433                    if ( nplane2 != mask18b ){
434                      //
435                      if ( nplane2 >= 37 ) nn2 = nplane2 + 1;
436                      vi2 = 1;
437                      if ( nn2%2 ) vi2 = 0;
438                      pl1 = (nn2 - 1 + vi2)/2;
439                      //
440                      cs2 = ct->tibar[pl2][vi2] - 1;
441                      //
442                      if ( ct->tibar[pl2][vi2] > -1 ){
443                        //
444                        cd2 = 95 - cs2;
445                        //
446                        Int_t t1 = TMath::Max(0,(cd2+0));
447                        Int_t t2 = TMath::Min(190,(cd2+95));
448                        mstrip2min = this->ConvertStrip(t1);
449                        mstrip2max = this->ConvertStrip(t2) + 1;
450                        //
451                        for (Int_t mstrip2=mstrip2min; mstrip2<mstrip2max; mstrip2++){
452                          //              
453                          mip2 = mipv[nplane2][mstrip2] - this->GetFullAverageAt(nplane2,mstrip2,rig,therigb);
454                          //
455                          if ( mip2 != 0. ){
456                            //
457    //                      Int_t sh = -5 + nplane1;
458    //                      if ( sh > 5 ) sh -= 11*nplane1;
459    //                      //
460    //                      mj = (nplane2 * 11) + mstrip2 + sh;
461    //                      //
462    //                      if ( mj < 0 ) mj += 473;
463    //                      if ( mj >= 473 ) mj -= 473;
464    //                      //
465                            mj = (nplane2 * 11) + mstrip2;
466                            //
467                            Float_t svalue =  mip1 * this->GetFullHmatrixAt(mi,mj,rig,mtherig,therigb) * mip2;
468                            //                      fulltzeta += mip1 * this->GetFullHmatrixAt(mi,mj,rig,therigb) * mip2;
469                            fulltzeta += svalue;
470                            if ( svalue < 0. ) negfulltzeta += svalue;
471                            if ( svalue > 0. ) posfulltzeta += svalue;
472                            //              (*fmatrix[rbi])[mi][mj] += (mip1 * mip2);
473                            //                      if ( mstrip1 == mstrip2 ) printf("\n\n=> nplane1 %i nplane2 %i mstrip1 %i mstrip2 %i \n => mip1 %f H %f mip2 %f \n => mipv1 %f ave1 %f mipv2 %f ave2 %f \n => rig %f rigbin %i mtherigb %i\n",nplane1,nplane2,mstrip1,mstrip2,mip1,this->GetFullHmatrixAt(mi,mj,rig,mtherig),mip2,mipv[nplane1][mstrip1],this->GetFullAverageAt(nplane1,mstrip1,rig,therigb),mipv[nplane2][mstrip2],this->GetFullAverageAt(nplane2,mstrip2,rig,therigb),rig,therigb,mtherig);
474                            //
475                          };
476                        };
477                      };
478                    };
479                  };
480                };
481              };
482            };
483          };
484        };
485        //
486        fdegfre = dgf*11;
487        if ( numpos ) aveposvar /= numpos;
488        if ( numneg ) avenegvar /= numneg;
489        //
490      };
491    //    //
492  //  if ( debug ) this->Print();    //  if ( debug ) this->Print();
493    if ( debug ) printf(" exit \n");    if ( debug ) printf(" exit \n");
494  }  }
495    
# Line 249  void CaloFranzini::Process(Int_t itr){ Line 497  void CaloFranzini::Process(Int_t itr){
497  Float_t CaloFranzini::GetNormLongTZeta(){  Float_t CaloFranzini::GetNormLongTZeta(){
498    Process();    Process();
499    Float_t normz = 0.;    Float_t normz = 0.;
500    if ( degfre > 0 ) normz = longtzeta/(Float_t)degfre;    if ( degfre != 0 ) normz = longtzeta/(Float_t)degfre;
501    return normz;    return normz;
502  }  }
503    
504  Float_t CaloFranzini::GetNormFullTZeta(){  Float_t CaloFranzini::GetNormFullTZeta(){
505    Process();    Process();
506    Float_t normz = 0.;    Float_t normz = 0.;
507    if ( degfre > 0 ) normz = fulltzeta/(Float_t)degfre;    if ( fdegfre != 0 ) normz = fulltzeta/(Float_t)fdegfre;
508    return normz;    return normz;
509  }  }
510    
511  Bool_t CaloFranzini::CreateMatrixFile(TString matrixfile){    Bool_t CaloFranzini::CreateMatrixFile(TString matrixfile){  
512    //    //
513    file = new TFile(matrixfile.Data(),"READ");    lfile = new TFile(matrixfile.Data(),"READ");
514    //    //
515    if ( !file || file->IsZombie() ){    if ( !lfile || lfile->IsZombie() ){
516      file = new TFile(matrixfile.Data(),"RECREATE");      if ( dolong ){
517      printf(" Create file %s \n",file->GetName());        lfile = new TFile(matrixfile.Data(),"RECREATE");
518          printf(" Create file %s \n",lfile->GetName());
519        };
520        if ( dofull ){
521          ffile = new TFile(matrixfile.Data(),"RECREATE");
522          printf(" Create file %s \n",ffile->GetName());
523        };
524    } else {    } else {
525        lfile->Close();
526      printf(" ERROR: file %s already existing!\n Choose another name or delete the old file\n",matrixfile.Data());      printf(" ERROR: file %s already existing!\n Choose another name or delete the old file\n",matrixfile.Data());
527      return(false);      return(false);
528    };    };
# Line 278  Bool_t CaloFranzini::CreateMatrixFile(TS Line 533  Bool_t CaloFranzini::CreateMatrixFile(TS
533    
534  Bool_t CaloFranzini::UpdateMatrixFile(TString matrixfile){    Bool_t CaloFranzini::UpdateMatrixFile(TString matrixfile){  
535    //    //
536    file = new TFile(matrixfile.Data(),"UPDATE");    if ( dolong ){
537        lfile = new TFile(matrixfile.Data(),"UPDATE");
538        //
539        if ( !lfile || lfile->IsZombie() ){
540          printf(" ERROR: file %s already existing!\n Choose another name or delete the old file\n",matrixfile.Data());
541          return(false);
542        };
543      };
544    //    //
545    if ( !file || file->IsZombie() ){    if ( dofull ){
546      printf(" ERROR: file %s already existing!\n Choose another name or delete the old file\n",matrixfile.Data());      ffile = new TFile(matrixfile.Data(),"UPDATE");
547      return(false);      //
548        if ( !ffile || ffile->IsZombie() ){
549          printf(" ERROR: file %s already existing!\n Choose another name or delete the old file\n",matrixfile.Data());
550          return(false);
551        };
552    };    };
553    //    //
554    return(true);    return(true);
# Line 290  Bool_t CaloFranzini::UpdateMatrixFile(TS Line 556  Bool_t CaloFranzini::UpdateMatrixFile(TS
556  }  }
557    
558  void CaloFranzini::WriteNumBin(Int_t numbin){  void CaloFranzini::WriteNumBin(Int_t numbin){
559    file->cd();    if ( dolong ){
560    TArrayI nbi(1, &numbin);      lfile->cd();
561    file->WriteObject(&nbi, "nbinenergy");      TArrayI nbi(1, &numbin);
562        lfile->WriteObject(&nbi, "nbinenergy");
563      };
564      if ( dofull ){
565        ffile->cd();
566        TArrayI nbi(1, &numbin);
567        ffile->WriteObject(&nbi, "nbinenergy");
568      };
569  }  }
570    
571  void CaloFranzini::WriteRigBin(TArrayF *rigbin){  void CaloFranzini::WriteRigBin(TArrayF *rigbin){
572    file->cd();    if ( dolong ){
573    //  rigbin->Write("binrig");      lfile->cd();
574    file->WriteObject(&(*rigbin), "binrig");      //  rigbin->Write("binrig");
575        lfile->WriteObject(&(*rigbin), "binrig");
576      };
577      if ( dofull ){
578        ffile->cd();
579        //  rigbin->Write("binrig");
580        ffile->WriteObject(&(*rigbin), "binrig");
581      };
582  }  }
583    
584  void CaloFranzini::WriteLongMean(TArrayF *qpl, Int_t bin){  void CaloFranzini::WriteLongMean(TArrayF *qpl, Int_t bin){
585    file->cd();    lfile->cd();
586    TString name = Form("qplmeann%i",bin);    TString name = Form("qplmeann%i",bin);
587    file->WriteObject(&(*qpl),name.Data());    lfile->WriteObject(&(*qpl),name.Data());
588  }  }
589    
590  void CaloFranzini::WriteFullMean(TMatrixD *qpl, Int_t bin){  void CaloFranzini::WriteFullMean(TMatrixD *qpl, Int_t bin){
591    file->cd();    ffile->cd();
592    TString name = Form("fqplmeann%i",bin);    TString name = Form("fqplmeann%i",bin);
593    file->WriteObject(&(*qpl),name.Data());    ffile->WriteObject(&(*qpl),name.Data());
594      ffile->Purge();
595    }
596    
597    void CaloFranzini::WriteFullNMean(TMatrixD *qpl, Int_t bin){
598      ffile->cd();
599      TString name = Form("fnqplmeann%i",bin);
600      ffile->WriteObject(&(*qpl),name.Data());
601      ffile->Purge();
602  }  }
603    
604  void CaloFranzini::WriteInvertedLongMatrix(TMatrixD mat, Int_t bin){  void CaloFranzini::WriteInvertedLongMatrix(TMatrixD mat, Int_t bin){
605    file->cd();    lfile->cd();
606    TString name = Form("matrixn%i",bin);    TString name = Form("matrixn%i",bin);
607    //  mat.Write(name.Data());    //  mat.Write(name.Data());
608    file->WriteObject(&mat,name.Data());    lfile->WriteObject(&mat,name.Data());
609  }  }
610    
611  void CaloFranzini::WriteInvertedFullMatrix(TMatrixF mat, Int_t bin){  void CaloFranzini::WriteInvertedFullMatrix(TMatrixD mat, Int_t bin){
612    file->cd();    ffile->cd();
613    TString name = Form("fmatrixn%i",bin);    TString name = Form("fmatrixn%i",bin);
614    //  mat.Write(name.Data());    //  mat.Write(name.Data());
615    file->WriteObject(&mat,name.Data());    ffile->WriteObject(&mat,name.Data());
616  }  }
617    
618  void CaloFranzini::WriteLongMatrix(TMatrixD *mat, Int_t bin){  void CaloFranzini::WriteLongMatrix(TMatrixD *mat, Int_t bin){
619    file->cd();    lfile->cd();
620    TString name = Form("origmatrixn%i",bin);    TString name = Form("origmatrixn%i",bin);
621    //  mat.Write(name.Data());    //  mat.Write(name.Data());
622    file->WriteObject(&(*mat),name.Data());    lfile->WriteObject(&(*mat),name.Data());
623  }  }
624    
625  void CaloFranzini::WriteFullMatrix(TMatrixF *mat, Int_t bin){  void CaloFranzini::WriteFullMatrix(TMatrixD *mat, Int_t bin){
626    file->cd();    ffile->cd();
627    TString name = Form("origfmatrixn%i",bin);    TString name = Form("origfmatrixn%i",bin);
628    //  mat.Write(name.Data());    //  mat.Write(name.Data());
629    file->WriteObject(&(*mat),name.Data());    ffile->WriteObject(&(*mat),name.Data());
630      ffile->Purge();
631  }  }
632    
633  void CaloFranzini::CloseMatrixFile(){  void CaloFranzini::WriteFullNMatrix(TMatrixF *mat, Int_t bin){
634    file->cd();    ffile->cd();
635    file->Close();    TString name = Form("origfnmatrixn%i",bin);
636      //  mat.Write(name.Data());
637      ffile->WriteObject(&(*mat),name.Data());
638      ffile->Purge();
639  }  }
640    
641    void CaloFranzini::CloseMatrixFile(){
642      if ( dolong ){
643        lfile->cd();
644        lfile->Close();
645      };
646      if ( dofull ){
647        ffile->cd();
648        ffile->Close();
649      };
650    }
651    
652  Bool_t CaloFranzini::Open(TString matrixfile){    Bool_t CaloFranzini::Open(TString matrixfile){  
653      return(this->Open(matrixfile,""));
654    }
655    
656    Bool_t CaloFranzini::Open(TString longmatrixfile, TString fullmatrixfile){  
657    //    //
658    // find matrix file    // find matrix file
659    //    //
660    if ( !strcmp(matrixfile.Data(),"") ){    if ( !strcmp(longmatrixfile.Data(),"") ){
661      if (dolong) matrixfile = (TString)gSystem->ExpandPathName("$PAM_CALIB")+"/cal-param/covmatrix_longel.root";          if (dolong) longmatrixfile = (TString)gSystem->ExpandPathName("$PAM_CALIB")+"/cal-param/covmatrix_longel.root";    
     if (dofull) matrixfile = (TString)gSystem->ExpandPathName("$PAM_CALIB")+"/cal-param/covmatrix_fullel.root";      
662    };    };
663    //    if ( !strcmp(fullmatrixfile.Data(),"") ){
664    file = new TFile(matrixfile.Data(),"READ");      if (dofull) fullmatrixfile = (TString)gSystem->ExpandPathName("$PAM_CALIB")+"/cal-param/covmatrix_fullel.root";    
   //  
   if ( !file || file->IsZombie() ){  
     printf(" ERROR: cannot open file %s \n",matrixfile.Data());  
     return(false);  
   };  
   //  
   if ( !this->LoadBin() ){  
     printf(" %s \n",matrixfile.Data());  
     return(false);  
665    };    };
666    //    //
667    if ( dolong ){    if ( dolong ){
668        //
669        lfile = new TFile(longmatrixfile.Data(),"READ");
670        //
671        if ( !lfile || lfile->IsZombie() ){
672          printf(" ERROR: cannot open file %s \n",longmatrixfile.Data());
673          return(false);
674        };
675        //
676        if ( !this->LoadBin() ){
677          printf(" %s \n",longmatrixfile.Data());
678          return(false);
679        };
680        //
681      if ( !this->LoadLong() ){      if ( !this->LoadLong() ){
682        printf(" %s \n",matrixfile.Data());        printf(" %s \n",longmatrixfile.Data());
683        return(false);        return(false);
684      };      };
685      //      //
686      if ( !this->LoadMatrices() ){      if ( !this->LoadMatrices() ){
687        printf(" %s \n",matrixfile.Data());        printf(" %s \n",longmatrixfile.Data());
688        return(false);        return(false);
689      };      };
690    };    };
691    //    //
692    if ( dofull ){    if ( dofull ){
693        //
694        ffile = new TFile(fullmatrixfile.Data(),"READ");
695        //
696        if ( !ffile || ffile->IsZombie() ){
697          printf(" ERROR: cannot open file %s \n",fullmatrixfile.Data());
698          return(false);
699        };
700        //
701        if ( !dolong ){
702          //
703          if ( !this->LoadBin(true) ){
704            printf(" %s \n",fullmatrixfile.Data());
705            return(false);
706          };      
707          //
708        };
709      if ( !this->LoadFull() ){      if ( !this->LoadFull() ){
710        printf(" %s \n",matrixfile.Data());        printf(" %s \n",fullmatrixfile.Data());
711        return(false);        return(false);
712      };      };
713      //      //
714      if ( !this->LoadFullMatrices() ){      if ( !this->LoadFullMatrices() ){
715        printf(" %s \n",matrixfile.Data());        printf(" %s \n",fullmatrixfile.Data());
716        return(false);        return(false);
717      };      };
718    };    };
# Line 397  Bool_t CaloFranzini::Open(TString matrix Line 722  Bool_t CaloFranzini::Open(TString matrix
722    //    //
723  }  }
724    
725    
726  Bool_t CaloFranzini::LoadBin(){  Bool_t CaloFranzini::LoadBin(){
727      return(this->LoadBin(false));
728    }
729    
730    Bool_t CaloFranzini::LoadBin(Bool_t full){
731    //    //
732    TArrayI *numbin = (TArrayI*)file->Get("nbinenergy");    if ( !full ) {
733    if ( !numbin ){      //
734      printf(" ERROR: cannot read number of bins from file ");      TArrayI *numbin = (TArrayI*)lfile->Get("nbinenergy");
735      return(false);      if ( !numbin ){
736    };        printf(" ERROR: cannot read number of bins from file ");
737    nbin = (Int_t)numbin->At(0);        return(false);
738    if ( nbin <= 0 ){      };
739      printf(" ERROR: cannot work with 0 energy bins from file ");      nbin = (Int_t)numbin->At(0);
740      return(false);      if ( nbin <= 0 ){
741    };        printf(" ERROR: cannot work with 0 energy bins from file ");
742    //        return(false);
743    brig = (TArrayF*)file->Get("binrig");      };
744    if ( !brig ){      //
745      printf(" ERROR: cannot read rigidity binning from file ");      brig = (TArrayF*)lfile->Get("binrig");
746      return(false);      if ( !brig ){
747    };        printf(" ERROR: cannot read rigidity binning from file ");
748    //        return(false);
749    brigm=(TArrayF*)file->Get("binrigmean");      };
750    if ( !brigm ){      //
751      printf(" ERROR: cannot read mean rigidity binning from file ");      brigm=(TArrayF*)lfile->Get("binrigmean");
752      return(false);      if ( !brigm ){
753          printf(" ERROR: cannot read mean rigidity binning from file ");
754          return(false);
755        };
756        //
757      } else {
758        //
759        TArrayI *numbin = (TArrayI*)ffile->Get("nbinenergy");
760        if ( !numbin ){
761          printf(" ERROR: cannot read number of bins from file ");
762          return(false);
763        };
764        nbin = (Int_t)numbin->At(0);
765        if ( nbin <= 0 ){
766          printf(" ERROR: cannot work with 0 energy bins from file ");
767          return(false);
768        };
769        //
770        brig = (TArrayF*)ffile->Get("binrig");
771        if ( !brig ){
772          printf(" ERROR: cannot read rigidity binning from file ");
773          return(false);
774        };
775        //
776        brigm=(TArrayF*)ffile->Get("binrigmean");
777        if ( !brigm ){
778          printf(" ERROR: cannot read mean rigidity binning from file ");
779          return(false);
780        };
781    };    };
782    //    //
783    return(true);    return(true);
# Line 429  Bool_t CaloFranzini::LoadLong(){ Line 787  Bool_t CaloFranzini::LoadLong(){
787    //    //
788    for (Int_t i=0;i<17;i++){    for (Int_t i=0;i<17;i++){
789      TString name = Form("qplmeann%i",i);      TString name = Form("qplmeann%i",i);
790      qplmean[i] = (TArrayF*)file->Get(name.Data());      qplmean[i] = (TArrayF*)lfile->Get(name.Data());
791      if ( !qplmean[i] ){      if ( !qplmean[i] ){
792        printf(" ERROR: cannot read average from file ");        printf(" ERROR: cannot read average from file ");
793        return(false);        return(false);
# Line 443  Bool_t CaloFranzini::LoadFull(){ Line 801  Bool_t CaloFranzini::LoadFull(){
801    //    //
802    for (Int_t i=0;i<17;i++){    for (Int_t i=0;i<17;i++){
803      TString name = Form("fqplmeann%i",i);      TString name = Form("fqplmeann%i",i);
804      fqplmean[i] = (TMatrixD*)file->Get(name.Data());      fqplmean[i] = (TMatrixD*)ffile->Get(name.Data());
805      if ( !fqplmean[i] ){      if ( !fqplmean[i] ){
806        printf(" ERROR: cannot read average from file ");        printf(" ERROR: cannot read average from file ");
807        return(false);        return(false);
# Line 457  Bool_t CaloFranzini::LoadMatrices(){ Line 815  Bool_t CaloFranzini::LoadMatrices(){
815    //    //
816    for (Int_t i=0;i<17;i++){    for (Int_t i=0;i<17;i++){
817          TString name1 = Form("matrixn%i",i);          TString name1 = Form("matrixn%i",i);
818          hmat[i] = (TMatrixD*)file->Get(name1.Data());          hmat[i] = (TMatrixD*)lfile->Get(name1.Data());
819    };    };
820    //    //
821    return(true);    return(true);
# Line 467  Bool_t CaloFranzini::LoadFullMatrices(){ Line 825  Bool_t CaloFranzini::LoadFullMatrices(){
825    //    //
826    for (Int_t i=0;i<17;i++){    for (Int_t i=0;i<17;i++){
827          TString name1 = Form("fmatrixn%i",i);          TString name1 = Form("fmatrixn%i",i);
828          hfmat[i] = (TMatrixF*)file->Get(name1.Data());          hfmat[i] = (TMatrixD*)ffile->Get(name1.Data());
829    };    };
830    //    //
831    return(true);    return(true);
# Line 495  TMatrixD *CaloFranzini::LoadCovarianceMa Line 853  TMatrixD *CaloFranzini::LoadCovarianceMa
853    //    //
854  }  }
855    
856    TMatrixD *CaloFranzini::LoadFullAverage(Int_t rigbin){
857      //
858      TString name = Form("fqplmeann%i",rigbin);
859      TMatrixD *fmean=(TMatrixD*)ffile->Get(name.Data());
860      //
861      return(fmean);
862      //
863    }
864    
865    TMatrixF *CaloFranzini::LoadFullMatrix(Int_t rigbin){
866      //
867      TString name = Form("origfmatrixn%i",rigbin);
868      TMatrixF *fmatri=(TMatrixF*)ffile->Get(name.Data());
869      //
870      return(fmatri);
871      //
872    }
873    
874    void CaloFranzini::LoadFullMatrix(Int_t rigbin, TMatrixF *fmatri){
875      //
876      TString name = Form("origfmatrixn%i",rigbin);
877      fmatri=(TMatrixF*)ffile->Get(name.Data());
878      //
879    }
880    
881    void CaloFranzini::UnLoadFullAverage(Int_t rigbin){
882      //
883      TString name = Form("fqplmeann%i",rigbin);
884      ffile->Delete(name.Data());
885      //
886    }
887    
888    void CaloFranzini::UnLoadFullMatrix(Int_t rigbin){
889      //
890      TString name = Form("origfmatrixn%i",rigbin);
891      ffile->Delete(name.Data());
892      //
893    }
894    
895    TMatrixF *CaloFranzini::LoadFullNMatrix(Int_t rigbin){
896      //
897      TString name = Form("origfnmatrixn%i",rigbin);
898      TMatrixF *fnmatri=(TMatrixF*)ffile->Get(name.Data());
899      //
900      return(fnmatri);
901      //
902    }
903    
904    void CaloFranzini::UnLoadFullNMatrix(Int_t rigbin){
905      //
906      TString name = Form("origfnmatrixn%i",rigbin);
907      ffile->Delete(name.Data());
908      //
909    }
910    
911    TMatrixD *CaloFranzini::LoadFullNAverage(Int_t rigbin){
912      //
913      TString name = Form("fnqplmeann%i",rigbin);
914      TMatrixD *fnmean=(TMatrixD*)ffile->Get(name.Data());
915      //
916      return(fnmean);
917      //
918    }
919    
920    void CaloFranzini::UnLoadFullNAverage(Int_t rigbin){
921      //
922      TString name = Form("fnqplmeann%i",rigbin);
923      ffile->Delete(name.Data());
924      //
925    }
926    
927    
928  TArrayF *CaloFranzini::LoadLongAverage(Float_t rig){  TArrayF *CaloFranzini::LoadLongAverage(Float_t rig){
929    //    //
# Line 558  Float_t CaloFranzini::GetAverageAt(Int_t Line 987  Float_t CaloFranzini::GetAverageAt(Int_t
987    Float_t b = log(maxene/minene)/(maxrig-minrig);    Float_t b = log(maxene/minene)/(maxrig-minrig);
988    Float_t a = minene/exp(minrig*b);    Float_t a = minene/exp(minrig*b);
989    Float_t ave = a*exp(b*rig);    Float_t ave = a*exp(b*rig);
990      if ( a == 0. || minene == 0. || ave != ave ){
991        //  if ( a == 0. || minene == 0. ){
992        Float_t m = (maxene-minene)/(maxrig-minrig);
993        Float_t q = minene - m * minrig;
994        ave = rig * m + q;
995      };
996    //    //
997    return(ave);    return(ave);
998    //    //
999  }  }
1000    
1001  Float_t CaloFranzini::GetFullAverageAt(Int_t plane, Int_t strip, Float_t rig){  Float_t CaloFranzini::GetFullAverageAt(Int_t plane, Int_t strip, Float_t rig){
1002    //    //
1003    Int_t therigb = 0;    Int_t therigb = 0;
# Line 572  Float_t CaloFranzini::GetFullAverageAt(I Line 1008  Float_t CaloFranzini::GetFullAverageAt(I
1008      };      };
1009    };    };
1010    //    //
   Float_t minrig;  
   Float_t maxrig;  
   //  
1011    //    //
1012    if ( rig < brigm->At(0) ){    if ( rig < brigm->At(0) ){
1013      if ( rig < brig->At(0) ){      if ( rig < brig->At(0) ){
# Line 589  Float_t CaloFranzini::GetFullAverageAt(I Line 1022  Float_t CaloFranzini::GetFullAverageAt(I
1022      therigb = nbin - 3;      therigb = nbin - 3;
1023    };    };
1024    //    //
1025      return(this->GetFullAverageAt(plane,strip,rig,therigb));
1026      //
1027    }
1028    
1029    Float_t CaloFranzini::GetFullAverageAt(Int_t plane, Int_t strip, Float_t rig, Int_t therigb){
1030      //
1031      Float_t minrig;
1032      Float_t maxrig;
1033      //
1034    minrig = brigm->At(therigb);    minrig = brigm->At(therigb);
1035    maxrig = brigm->At(therigb+1);    maxrig = brigm->At(therigb+1);
1036    //    //
# Line 602  Float_t CaloFranzini::GetFullAverageAt(I Line 1044  Float_t CaloFranzini::GetFullAverageAt(I
1044    Float_t b = log(maxene/minene)/(maxrig-minrig);    Float_t b = log(maxene/minene)/(maxrig-minrig);
1045    Float_t a = minene/exp(minrig*b);    Float_t a = minene/exp(minrig*b);
1046    Float_t ave = a*exp(b*rig);    Float_t ave = a*exp(b*rig);
1047      if ( a == 0. || minene == 0. || ave != ave ){
1048        Float_t m = (maxene-minene)/(maxrig-minrig);
1049        Float_t q = minene - m * minrig;
1050        ave = rig * m + q;
1051      };
1052      //  
1053      //  ave += (44.-plane)*strip;
1054      //if ( a == 0. ) ave = 0.;
1055      if ( !(ave == ave) ) printf("a %f b %f ave %f maxene %f minene %f maxrig %f minrig %f \n",a,b,ave,maxene,minene,maxrig,minrig);
1056    //    //
1057    return(ave);    return(ave);
1058    //    //
# Line 621  Float_t CaloFranzini::GetHmatrixAt(Int_t Line 1072  Float_t CaloFranzini::GetHmatrixAt(Int_t
1072    //    //
1073    if ( rig < brigm->At(0) ){    if ( rig < brigm->At(0) ){
1074      if ( rig < brig->At(0) ){      if ( rig < brig->At(0) ){
1075  //      printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brigm->At(0));        //      printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brigm->At(0));
1076      };      };
1077      therigb = 0;      therigb = 0;
1078    };    };
1079  //  if ( rig >= brigm->At(nbin-4) ){ // -2    //  if ( rig >= brigm->At(nbin-4) ){ // -2
1080    if ( rig >= brigm->At(nbin-2) ){    if ( rig >= brigm->At(nbin-2) ){
1081      if ( rig >= brig->At(nbin-2) ) {      if ( rig >= brig->At(nbin-2) ) {
1082  //      printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brigm->At(nbin-2));        //      printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brigm->At(nbin-2));
1083      };      };
1084  //    therigb = nbin - 5;// -3      //    therigb = nbin - 5;// -3
1085      therigb = nbin - 3;      therigb = nbin - 3;
1086    };    };
1087    //    //
# Line 673  Float_t CaloFranzini::GetFullHmatrixAt(I Line 1124  Float_t CaloFranzini::GetFullHmatrixAt(I
1124      };      };
1125    };    };
1126    //    //
   Float_t minrig;  
   Float_t maxrig;  
   //  
1127    if ( rig < brigm->At(0) ){    if ( rig < brigm->At(0) ){
1128      if ( rig < brig->At(0) ){      if ( rig < brig->At(0) ){
1129  //      printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brigm->At(0));        //      printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brigm->At(0));
1130      };      };
1131      therigb = 0;      therigb = 0;
1132    };    };
1133  //  if ( rig >= brigm->At(nbin-4) ){ // -2  //  if ( rig >= brigm->At(nbin-4) ){ // -2
1134    if ( rig >= brigm->At(nbin-2) ){    if ( rig >= brigm->At(nbin-2) ){
1135        //if ( rig >= brigm->At(nbin-5) ){
1136      if ( rig >= brig->At(nbin-2) ) {      if ( rig >= brig->At(nbin-2) ) {
1137  //      printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brigm->At(nbin-2));        //      printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brigm->At(nbin-2));
1138      };      };
1139  //    therigb = nbin - 5;// -3      //    therigb = nbin - 5;// -3
1140        //  therigb = nbin - 5;// -3
1141      therigb = nbin - 3;      therigb = nbin - 3;
1142    };    };
1143    //    //
1144    if ( therigb < 2 ) therigb = 2;    if ( therigb < 2 ) therigb = 2;
1145      if ( therigb > 13 ) therigb = 13;
1146      //
1147      return(this->GetFullHmatrixAt(iindex,jindex,rig,therigb));
1148      //
1149    }
1150    
1151    
1152    Float_t CaloFranzini::GetFullHmatrixAt(Int_t iindex, Int_t jindex, Float_t rig, Int_t therigb){
1153      //
1154      return(this->GetFullHmatrixAt(iindex,jindex,rig,therigb,0));
1155      //
1156    }
1157    
1158    Float_t CaloFranzini::GetFullHmatrixAt(Int_t iindex, Int_t jindex, Float_t rig, Int_t therigb, Int_t mtherig){  
1159      //
1160      //  Int_t lofit = 12;
1161      Int_t lofit = 3;
1162      //
1163      Float_t lowrig = 0.;
1164      Float_t minrig;
1165      Float_t maxrig;
1166      //
1167      if ( mtherig > lofit ) lowrig = brigm->At(therigb-1);
1168    minrig = brigm->At(therigb);    minrig = brigm->At(therigb);
1169    maxrig = brigm->At(therigb+1);    maxrig = brigm->At(therigb+1);
1170    //  printf(" therigb %i minrig %f maxrig %f rig %f %i i %i j \n",therigb,minrig,maxrig,rig,iindex,jindex);    //if ( therigb > 10 )  printf(" therigb %i minrig %f maxrig %f rig %f %i i %i j \n",therigb,minrig,maxrig,rig,iindex,jindex);
1171    //    //
1172      Float_t lowene = 0.;
1173      if ( mtherig > lofit ) lowene = (*hfmat[therigb-1])[iindex][jindex];
1174    Float_t minene = (*hfmat[therigb])[iindex][jindex];    Float_t minene = (*hfmat[therigb])[iindex][jindex];
1175    Float_t maxene = (*hfmat[therigb+1])[iindex][jindex];    Float_t maxene = (*hfmat[therigb+1])[iindex][jindex];
1176    //  printf(" therigb %i minrig %f maxrig %f minene %f maxene %f a %f b %f rig %f ave %f \n",therigb,minrig,maxrig,minene,maxene,a,b,rig,ave);    //  printf(" therigb %i minrig %f maxrig %f minene %f maxene %f a %f b %f rig %f ave %f \n",therigb,minrig,maxrig,minene,maxene,a,b,rig,ave);
1177    //    //
1178  //  Float_t a = 0.;    Float_t ave = 0.;
 //  Float_t b = 0.;  
 //  Float_t ave = 0.;  
 //  if ( minene == 0. ){  
 //        
 //  } else {  
 //      b = log(maxene/minene)/(maxrig-minrig);  
 //      a = minene/exp(minrig*b);  
 //      ave = a*exp(b*rig);  
 //  };  
1179    //    //
1180    Float_t m = (maxene-minene)/(maxrig-minrig);    //  Float_t a = 0.;
1181    Float_t q = minene - m * minrig;    //  Float_t b = 0.;
1182    Float_t ave = rig * m + q;    //  Float_t ave = 0.;
1183    if ( debug ) printf(" therigb %i minrig %f maxrig %f minene %f maxene %f a %f b %f rig %f ave %f \n",therigb,minrig,maxrig,minene,maxene,m,q,rig,ave);    //  if ( minene == 0. ){
1184      //    
1185      //  } else {
1186      //    b = log(maxene/minene)/(maxrig-minrig);
1187      //    a = minene/exp(minrig*b);
1188      //    ave = a*exp(b*rig);
1189      //  };
1190      //
1191      if ( mtherig > lofit ){
1192        //
1193        // FIT
1194        //
1195    //     Float_t x[3]={lowrig,minrig,maxrig};
1196    //     Float_t y[3]={lowene,minene,maxene};
1197    //     //
1198    //     TGraph *tg= new TGraph(3,x,y);
1199    //     //
1200    //     //     gStyle->SetLabelSize(0.04);
1201    //     //     gStyle->SetNdivisions(510,"XY");
1202    //     //     TCanvas *c = new TCanvas();
1203    //     //     c->Draw();
1204    //     //     c->cd();
1205    //     //     tg->Draw("AP");
1206    //     //
1207    //     TF1 *fun = new TF1("fun","pol2");
1208    //     tg->Fit("fun","QNC");
1209    //     //
1210    //     ave = fun->GetParameter(0) + rig * fun->GetParameter(1) + rig * rig * fun->GetParameter(2);
1211    //     //
1212    //     //    printf(" therigb %i ave %f rig %f lowrig %f minrig %f maxrig %f lowene %f minene %f maxene %f \n",therigb,ave,rig,lowrig,minrig,maxrig,lowene,minene,maxene);
1213    //     //
1214    //     //
1215    //     //     c->Modified();
1216    //     //     c->Update();
1217    //     //     gSystem->ProcessEvents();
1218    //     //     gSystem->Sleep(6000);
1219    //     //     c->Close();
1220    //     //     gStyle->SetLabelSize(0);
1221    //     //     gStyle->SetNdivisions(1,"XY");
1222    //     //
1223    //    tg->Delete();
1224    //    fun->Delete();
1225        Float_t mrigl = (lowrig+minrig)/2.;
1226        Float_t mrigu = (minrig+maxrig)/2.;
1227        Float_t menel = (lowene+minene)/2.;
1228        Float_t meneu = (minene+maxene)/2.;
1229        Float_t m = (meneu-menel)/(mrigu-mrigl);
1230        Float_t q = menel - m * mrigl;
1231        ave = rig * m + q;
1232        //
1233      } else {
1234        Float_t m = (maxene-minene)/(maxrig-minrig);
1235        Float_t q = minene - m * minrig;
1236        ave = rig * m + q;
1237        if ( debug ) printf(" therigb %i minrig %f maxrig %f minene %f maxene %f a %f b %f rig %f ave %f \n",therigb,minrig,maxrig,minene,maxene,m,q,rig,ave);
1238      };
1239    //    //
1240    //    //
1241    return(ave);    return(ave);
# Line 806  void CaloFranzini::DrawLongAverage(Int_t Line 1327  void CaloFranzini::DrawLongAverage(Int_t
1327    gStyle->SetNdivisions(1,"XY");    gStyle->SetNdivisions(1,"XY");
1328    //    //
1329  };  };
1330    
1331    Int_t CaloFranzini::ConvertStrip(Int_t mstrip){  
1332      //
1333      Int_t lastrip = 0;
1334    //   //
1335    //   if ( mstrip < 50 ) lastrip = 0;
1336    //   //
1337    //   if ( mstrip >= 50 && mstrip < 64 ) lastrip = 1;
1338    //   //
1339    //   if ( mstrip >= 64 && mstrip < 70 ) lastrip = 2;
1340    //   //
1341    //   if ( mstrip >= 70 && mstrip < 75 ) lastrip = 3;
1342    //   //
1343    //   if ( mstrip >= 75 && mstrip < 84 ){    
1344    //     lastrip = (int)trunc((mstrip - 75)/3.) + 4;
1345    //   };
1346    //   //
1347    //   if ( mstrip >= 84 && mstrip < 90 ){    
1348    //     lastrip = (int)trunc((mstrip - 84)/2.) + 7;
1349    //   };
1350    //   //
1351    //   if ( mstrip >= 90 && mstrip < 101 ){    
1352    //     lastrip = mstrip - 90 + 10;
1353    //   };
1354    //   //
1355    //   if ( mstrip >= 101 && mstrip < 107 ){    
1356    //     lastrip = (int)trunc((mstrip - 101)/2.) + 21;
1357    //   };
1358    //   //
1359    //   //
1360    //   if ( mstrip >= 107 && mstrip < 116 ){    
1361    //     lastrip = (int)trunc((mstrip - 107)/3.) + 24;
1362    //   };
1363    //   //
1364    //   if ( mstrip >= 116 && mstrip < 121 ) lastrip = 27;
1365    //   //
1366    //   if ( mstrip >= 121 && mstrip < 127 ) lastrip = 28;
1367    //   //
1368    //   if ( mstrip >= 127 && mstrip < 141 ) lastrip = 29;
1369    //   //
1370    //   if ( mstrip >= 141 ) lastrip = 30;
1371    //   //
1372      //
1373      if ( mstrip < 83 ) lastrip = 0;
1374      //
1375      if ( mstrip >= 83 && mstrip < 90 ) lastrip = 1;
1376      //
1377      if ( mstrip >= 90 && mstrip < 93 ) lastrip = 2;
1378      //
1379      if ( mstrip >= 93 && mstrip < 98 ){    
1380        lastrip = mstrip - 93 + 3;
1381      };
1382      //
1383      if ( mstrip >= 98 && mstrip < 101 ) lastrip = 8;
1384      //
1385      if ( mstrip >= 101 && mstrip < 107 ) lastrip = 9;
1386      //
1387      if ( mstrip >= 107 ) lastrip = 10;
1388      //
1389      return(lastrip);
1390    }
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.10
  ViewVC Help
Powered by ViewVC 1.1.23