/[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.2 by mocchiut, Tue Dec 4 12:42:54 2007 UTC revision 1.5 by mocchiut, Thu Jan 3 10:02:28 2008 UTC
# Line 18  CaloFranzini::CaloFranzini(PamLevel2 *l2 Line 18  CaloFranzini::CaloFranzini(PamLevel2 *l2
18    //    //
19    file = NULL;    file = NULL;
20    brig = NULL;    brig = NULL;
21      brigm = NULL;
22    nbin = 0;    nbin = 0;
23    //    //
24    L2 = l2p;    L2 = l2p;
# Line 30  CaloFranzini::CaloFranzini(PamLevel2 *l2 Line 31  CaloFranzini::CaloFranzini(PamLevel2 *l2
31    dolong = true;    dolong = true;
32    dofull = false;    dofull = false;
33    sntr = 0;    sntr = 0;
34      OBT = 0;
35      PKT = 0;
36      atime = 0;
37      //
38      crig = false;
39      sel = true;
40      cont = false;
41      N = 0;
42      NC = 43;
43      //
44      mask18b = -1;
45    //    //
46    Clear();    Clear();
47    //    //
# Line 37  CaloFranzini::CaloFranzini(PamLevel2 *l2 Line 49  CaloFranzini::CaloFranzini(PamLevel2 *l2
49    
50  void CaloFranzini::Clear(){  void CaloFranzini::Clear(){
51    //    //
52    OBT = 0;    longtzeta = 0.;
53    PKT = 0;    fulltzeta = 0.;
   atime = 0;  
   longtzeta = -100.;  
   fulltzeta = -100.;  
54    degfre = 0;    degfre = 0;
55    memset(estrip, 0, 2*22*96*sizeof(Float_t));    memset(estrip, 0, 2*22*96*sizeof(Float_t));
56    memset(qplane, 0, 2*22*sizeof(Float_t));    memset(qplane, 0, 43*sizeof(Float_t));
57    //    //
58  }  }
59    
# Line 55  void CaloFranzini::Print(){ Line 64  void CaloFranzini::Print(){
64    printf("========================================================================\n");    printf("========================================================================\n");
65    printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);    printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);
66    printf(" debug                [debug flag]:.. %i\n",debug);    printf(" debug                [debug flag]:.. %i\n",debug);
67      printf(" degree of freedom                :.. %i\n",this->GetDegreeOfFreedom());  
68    printf(" longtzeta                        :.. %f\n",longtzeta);      printf(" longtzeta                        :.. %f\n",longtzeta);  
69    printf(" fulltzeta                        :.. %f\n",fulltzeta);      printf(" longtzeta normalized             :.. %f\n",this->GetNormLongTZeta());  
70      //  printf(" fulltzeta                        :.. %f\n",fulltzeta);  
71      //  printf(" longtzeta normalized             :.. %f\n",this->GetNormFullTZeta());  
72    printf("========================================================================\n");    printf("========================================================================\n");
73    //    //
74  }  }
# Line 69  void CaloFranzini::Delete(){ Line 81  void CaloFranzini::Delete(){
81    //    //
82  }  }
83    
84    void CaloFranzini::SetNoWpreSampler(Int_t n){
85      Int_t nc2 = NC/2;
86      if ( NC >= 37 ) nc2 = (NC+1)/2;
87      if ( nc2+n < 23 ){
88        N = n;
89      } else {
90        printf(" ERROR! Calorimeter is made of 22 W planes\n");
91        printf(" you are giving N presampler = %i and N calo = %i \n",n,NC);
92        printf(" WARNING: using default values NWpre = 0, NWcalo = 22\n");
93        NC = 43;
94        N = 0;
95      };
96    }
97    
98    void CaloFranzini::SetNoWcalo(Int_t n){
99      if ( N+n < 23 ){
100        NC = n*2;
101        if ( NC >37 ) NC--;
102      } else {
103        printf(" ERROR! Calorimeter is made of 22 W planes\n");
104        printf(" you are giving N W presampler = %i and N W calo = %i \n",N,n);
105        printf(" WARNING: using default values NWpre = 0, NWcalo = 22\n");
106        NC = 43;
107        N = 0;
108      };
109    }
110    
111  void CaloFranzini::Process(){  void CaloFranzini::Process(){
112    this->Process(0);    this->Process(0);
# Line 111  void CaloFranzini::Process(Int_t itr){ Line 149  void CaloFranzini::Process(Int_t itr){
149    this->Clear();    this->Clear();
150    //    //
151    Float_t rig = L2->GetTrack(itr)->GetTrkTrack()->GetRigidity();    Float_t rig = L2->GetTrack(itr)->GetTrkTrack()->GetRigidity();
152      if ( crig ) rig = L2->GetCaloLevel2()->qtot/260.;
153    //    //
154    // Fill the estrip matrix    // Fill the estrip matrix
155    //    //
# Line 124  void CaloFranzini::Process(Int_t itr){ Line 163  void CaloFranzini::Process(Int_t itr){
163      mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);      mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
164      //      //
165      estrip[view][plane][strip] = mip;      estrip[view][plane][strip] = mip;
166      nplane = 1 - view + 2 * plane;      //    
167      qplane[nplane] += mip;      nplane = 1 - view + 2 * (plane - N);
168        if ( nplane > (37-(2*N)) ) nplane--;
169        //
170    //    if ( plane == (18+N) ) mip = 0.;
171        if ( nplane > -1 ) qplane[nplane] += mip;
172      //      //
173    };    };
174    //    //
# Line 133  void CaloFranzini::Process(Int_t itr){ Line 176  void CaloFranzini::Process(Int_t itr){
176    //    //
177    if ( dolong ){    if ( dolong ){
178      //      //
179      degfre = 44;      if ( cont ){
180          for (Int_t i=0; i<22; i++){
181            if ( i == (18+N) ){
182              mask18b =  1 + 2 * (i - N);
183              break;
184            };
185          };
186        };  
187        //  
188        if ( sel ){
189          for (Int_t i=0; i<22; i++){
190            if ( i == (18-N) ){
191              mask18b =  1 + 2 * (i - N);
192              break;
193            };
194          };
195        };
196        //
197        if ( mask18b == 37 ) mask18b = -1;
198        //
199        Int_t dgf = 43;
200      //      //
201      for (Int_t i=0; i < 22; i++){      for (Int_t i=0; i < 22; i++){
202        if ( L2->GetTrack(itr)->GetCaloTrack()->tibar[i][1] < 0 ){        if ( L2->GetTrack(itr)->GetCaloTrack()->tibar[i][1] < 0 ){
203          degfre = 2 * i;          dgf = 2 * i;
204          break;          break;
205        };        };
206        if ( L2->GetTrack(itr)->GetCaloTrack()->tibar[i][0] < 0 ){        if ( L2->GetTrack(itr)->GetCaloTrack()->tibar[i][0] < 0 ){
207        degfre = 1 + 2 * i;        dgf = 1 + 2 * i;
208        break;        break;
209        };        };
210      };      };
211      //      //
212        if ( dgf < 43 && dgf > 37 ) dgf--;
213        //
214        degfre = TMath::Min(dgf,NC);
215        //
216        Float_t longzdiag = 0.;
217        Float_t longzout = 0.;
218        //
219      if ( degfre > 0 ){      if ( degfre > 0 ){
       TArrayF *qplmean = this->LoadLongAverage(rig);  
       TMatrixD *lmatrix = this->LoadCovarianceMatrix(rig);  
220        for (Int_t i = 0; i < degfre; i++){        for (Int_t i = 0; i < degfre; i++){
221          for (Int_t j = 0; j < degfre; j++){              if ( i != mask18b ){
222            longtzeta += (qplane[i]-qplmean->At(i)) * (*lmatrix)[i][j] * (qplane[j]-qplmean->At(j));            for (Int_t j = 0; j < degfre; j++){  
223                if ( j != mask18b ){
224                  if ( i == j ){
225                    longzdiag += (qplane[i]-this->GetAverageAt(i,rig)) * this->GetHmatrixAt(i,j,rig) * (qplane[j]-this->GetAverageAt(j,rig));
226                    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));
227                  } else {
228                    longzout += (qplane[i]-this->GetAverageAt(i,rig)) * this->GetHmatrixAt(i,j,rig) * (qplane[j]-this->GetAverageAt(j,rig));
229                    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));
230                  };
231                  longtzeta += (qplane[i]-this->GetAverageAt(i,rig)) * this->GetHmatrixAt(i,j,rig) * (qplane[j]-this->GetAverageAt(j,rig));
232                };
233              };
234          };          };
235        };        };
236          if ( debug ) printf(" diagonal total %f out of diagonal total %f \n",longzdiag,longzout);
237      };      };
238    };    };
239    if ( dofull ){    if ( dofull ){
# Line 161  void CaloFranzini::Process(Int_t itr){ Line 241  void CaloFranzini::Process(Int_t itr){
241      printf(" ERROR: CaloFranzini variables _NOT_ filled \n");      printf(" ERROR: CaloFranzini variables _NOT_ filled \n");
242    };      };  
243    //    //
244    if ( debug ) this->Print();  //  if ( debug ) this->Print();
245    if ( debug ) printf(" exit \n");    if ( debug ) printf(" exit \n");
246  }  }
247    
# Line 186  Bool_t CaloFranzini::CreateMatrixFile(TS Line 266  Bool_t CaloFranzini::CreateMatrixFile(TS
266    //    //
267    if ( !file || file->IsZombie() ){    if ( !file || file->IsZombie() ){
268      file = new TFile(matrixfile.Data(),"RECREATE");      file = new TFile(matrixfile.Data(),"RECREATE");
269        printf(" Create file %s \n",file->GetName());
270    } else {    } else {
271      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());
272      return(false);      return(false);
# Line 195  Bool_t CaloFranzini::CreateMatrixFile(TS Line 276  Bool_t CaloFranzini::CreateMatrixFile(TS
276    //    //
277  }  }
278    
279    Bool_t CaloFranzini::UpdateMatrixFile(TString matrixfile){  
280      //
281      file = new TFile(matrixfile.Data(),"UPDATE");
282      //
283      if ( !file || file->IsZombie() ){
284        printf(" ERROR: file %s already existing!\n Choose another name or delete the old file\n",matrixfile.Data());
285        return(false);
286      };
287      //
288      return(true);
289      //
290    }
291    
292  void CaloFranzini::WriteNumBin(Int_t numbin){  void CaloFranzini::WriteNumBin(Int_t numbin){
293      file->cd();
294    TArrayI nbi(1, &numbin);    TArrayI nbi(1, &numbin);
295    file->WriteObject(&nbi, "nbinenergy");    file->WriteObject(&nbi, "nbinenergy");
296  }  }
297    
298  void CaloFranzini::WriteRigBin(TArrayF *rigbin){  void CaloFranzini::WriteRigBin(TArrayF *rigbin){
299    file->WriteObject(&rigbin, "binrig");    file->cd();
300      //  rigbin->Write("binrig");
301      file->WriteObject(&(*rigbin), "binrig");
302  }  }
303    
304  void CaloFranzini::WriteLongMean(TArrayF *qpl, Int_t bin){  void CaloFranzini::WriteLongMean(TArrayF *qpl, Int_t bin){
305      file->cd();
306    TString name = Form("qplmeann%i",bin);    TString name = Form("qplmeann%i",bin);
307    file->WriteObject(&qpl,name.Data());    file->WriteObject(&(*qpl),name.Data());
308  }  }
309    
310  void CaloFranzini::WriteLongMatrix(TMatrixD *mat, Int_t bin){  void CaloFranzini::WriteFullMean(TMatrixD *qpl, Int_t bin){
311      file->cd();
312      TString name = Form("fqplmeann%i",bin);
313      file->WriteObject(&(*qpl),name.Data());
314    }
315    
316    void CaloFranzini::WriteInvertedLongMatrix(TMatrixD mat, Int_t bin){
317      file->cd();
318    TString name = Form("matrixn%i",bin);    TString name = Form("matrixn%i",bin);
319      //  mat.Write(name.Data());
320      file->WriteObject(&mat,name.Data());
321    }
322    
323    void CaloFranzini::WriteInvertedFullMatrix(TMatrixF mat, Int_t bin){
324      file->cd();
325      TString name = Form("fmatrixn%i",bin);
326      //  mat.Write(name.Data());
327      file->WriteObject(&mat,name.Data());
328    }
329    
330    void CaloFranzini::WriteLongMatrix(TMatrixD *mat, Int_t bin){
331      file->cd();
332      TString name = Form("origmatrixn%i",bin);
333      //  mat.Write(name.Data());
334      file->WriteObject(&(*mat),name.Data());
335    }
336    
337    void CaloFranzini::WriteFullMatrix(TMatrixF *mat, Int_t bin){
338      file->cd();
339      TString name = Form("origfmatrixn%i",bin);
340      //  mat.Write(name.Data());
341    file->WriteObject(&(*mat),name.Data());    file->WriteObject(&(*mat),name.Data());
342  }  }
343    
344  void CaloFranzini::CloseMatrixFile(){  void CaloFranzini::CloseMatrixFile(){
345      file->cd();
346    file->Close();    file->Close();
347  }  }
348    
# Line 235  Bool_t CaloFranzini::Open(TString matrix Line 363  Bool_t CaloFranzini::Open(TString matrix
363      return(false);      return(false);
364    };    };
365    //    //
366      if ( !this->LoadBin() ){
367        printf(" %s \n",matrixfile.Data());
368        return(false);
369      };
370      //
371      if ( dolong ){
372        if ( !this->LoadLong() ){
373          printf(" %s \n",matrixfile.Data());
374          return(false);
375        };
376        //
377        if ( !this->LoadMatrices() ){
378          printf(" %s \n",matrixfile.Data());
379          return(false);
380        };
381      };
382      //
383      if ( dofull ){
384        if ( !this->LoadFull() ){
385          printf(" %s \n",matrixfile.Data());
386          return(false);
387        };
388        //
389        if ( !this->LoadFullMatrices() ){
390          printf(" %s \n",matrixfile.Data());
391          return(false);
392        };
393      };
394      //
395      //
396      return(true);
397      //
398    }
399    
400    Bool_t CaloFranzini::LoadBin(){
401      //
402    TArrayI *numbin = (TArrayI*)file->Get("nbinenergy");    TArrayI *numbin = (TArrayI*)file->Get("nbinenergy");
403    if ( !numbin ){    if ( !numbin ){
404      printf(" ERROR: cannot read number of bins from file %s \n",matrixfile.Data());      printf(" ERROR: cannot read number of bins from file ");
405      return(false);      return(false);
406    };    };
407    nbin = (Int_t)numbin->At(0);    nbin = (Int_t)numbin->At(0);
408    if ( nbin <= 0 ){    if ( nbin <= 0 ){
409      printf(" ERROR: cannot work with 0 energy bins (from file %s) \n",matrixfile.Data());      printf(" ERROR: cannot work with 0 energy bins from file ");
410      return(false);      return(false);
411    };    };
412    //    //
413    brig = (TArrayF*)file->Get("binrig");    brig = (TArrayF*)file->Get("binrig");
414    if ( !brig ){    if ( !brig ){
415      printf(" ERROR: cannot read rigidity binning from file %s \n",matrixfile.Data());      printf(" ERROR: cannot read rigidity binning from file ");
416        return(false);
417      };
418      //
419      brigm=(TArrayF*)file->Get("binrigmean");
420      if ( !brigm ){
421        printf(" ERROR: cannot read mean rigidity binning from file ");
422      return(false);      return(false);
423    };    };
424    //    //
425    return(true);    return(true);
426    }
427    
428    Bool_t CaloFranzini::LoadLong(){
429      //
430      for (Int_t i=0;i<17;i++){
431        TString name = Form("qplmeann%i",i);
432        qplmean[i] = (TArrayF*)file->Get(name.Data());
433        if ( !qplmean[i] ){
434          printf(" ERROR: cannot read average from file ");
435          return(false);
436        };
437      };
438      //
439      return(true);
440    }
441    
442    Bool_t CaloFranzini::LoadFull(){
443      //
444      for (Int_t i=0;i<17;i++){
445        TString name = Form("fqplmeann%i",i);
446        fqplmean[i] = (TMatrixD*)file->Get(name.Data());
447        if ( !fqplmean[i] ){
448          printf(" ERROR: cannot read average from file ");
449          return(false);
450        };
451      };
452      //
453      return(true);
454    }
455    
456    Bool_t CaloFranzini::LoadMatrices(){
457      //
458      for (Int_t i=0;i<17;i++){
459            TString name1 = Form("matrixn%i",i);
460            hmat[i] = (TMatrixD*)file->Get(name1.Data());
461      };
462      //
463      return(true);
464    }
465    
466    Bool_t CaloFranzini::LoadFullMatrices(){
467    //    //
468      for (Int_t i=0;i<17;i++){
469            TString name1 = Form("fmatrixn%i",i);
470            hfmat[i] = (TMatrixF*)file->Get(name1.Data());
471      };
472      //
473      return(true);
474  }  }
475    
476  TMatrixD *CaloFranzini::LoadCovarianceMatrix(Float_t rig){  TMatrixD *CaloFranzini::LoadCovarianceMatrix(Float_t rig){
477    //    //
478    TString name;    Int_t mv = 0;
479    for (Int_t i = 0; i<nbin-1; i++){    for (Int_t i = 0; i<nbin-1; i++){
480      if ( rig>=brig->At(i) && rig < brig->At(i+1) ){      if ( rig>=brig->At(i) && rig < brig->At(i+1) ){
481        name = Form("matrixn%i",i);        mv = i;
482        break;        break;
483      };      };
484    };    };
485    if ( rig < brig->At(0) ){    if ( rig < brig->At(0) ){
486      printf(" WARNING: Event with rigidity lower than the first covariance matrix bin (rig = %f, lower limit = %f)\n",rig,brig->At(0));      printf(" WARNING: Event with rigidity lower than the first covariance matrix bin (rig = %f, lower limit = %f)\n",rig,brig->At(0));
487      name = "matrixn0";      mv = 0;
     printf(" Using matrix %s \n",name.Data());  
488    };    };
489    if ( rig >= brig->At(nbin) ){    if ( rig >= brig->At(nbin-1) ){
490      printf(" WARNING: Event with rigidity higher than the last covariance matrix bin (rig = %f, upper limit = %f)\n",rig,brig->At(nbin));      printf(" WARNING: Event with rigidity higher than the last covariance matrix bin (rig = %f, upper limit = %f)\n",rig,brig->At(nbin-1));
491      name = Form("matrixn%i",nbin-1);      mv = nbin-2;
     printf(" Using matrix %s \n",name.Data());  
492    };    };
493    //    //
494    TMatrixD *matrix = (TMatrixD*)file->Get(name.Data());    return(hmat[mv]);
   //  
   return(matrix);  
495    //    //
496  }  }
497    
498    
499  TArrayF *CaloFranzini::LoadLongAverage(Float_t rig){  TArrayF *CaloFranzini::LoadLongAverage(Float_t rig){
500    //    //
501    TString name;    Int_t mv=0;
502    for (Int_t i = 0; i<nbin-1; i++){    for (Int_t i = 0; i<nbin-1; i++){
503      if ( rig>=brig->At(i) && rig < brig->At(i+1) ){      if ( rig>=brig->At(i) && rig < brig->At(i+1) ){
504        name = Form("qplmeann%i",i);          mv = i;
505        break;        break;
506      };      };
507    };    };
508    if ( rig < brig->At(0) ){    if ( rig < brig->At(0) ){
509      printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brig->At(0));      printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brig->At(0));
510      name = "qplmeann0";      mv = 0;
511      printf(" Using qplmean %s \n",name.Data());    };
512      if ( rig >= brig->At(nbin-1) ){
513        printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brig->At(nbin-1));
514        mv=nbin-2;
515      };
516      //
517      return(qplmean[mv]);
518      //
519    }
520    
521    Float_t CaloFranzini::GetAverageAt(Int_t plane, Float_t rig){
522      //
523      Int_t therigb = 0;
524      for (Int_t i = 0; i<nbin-2; i++){
525        if ( rig>=brigm->At(i) && rig < brigm->At(i+1) ){
526          therigb = i;
527          break;
528        };
529      };
530      //
531      Float_t minrig;
532      Float_t maxrig;
533      //
534      //
535      if ( rig < brigm->At(0) ){
536        if ( rig < brig->At(0) ){
537    //      printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brigm->At(0));
538        };
539        therigb = 0;
540      };
541      if ( rig >= brigm->At(nbin-2) ){
542        if ( rig >= brig->At(nbin-2) ) {
543    //      printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brigm->At(nbin-2));
544        };
545        therigb = nbin - 3;
546      };
547      //
548      minrig = brigm->At(therigb);
549      maxrig = brigm->At(therigb+1);
550      //
551      Float_t minene = (*qplmean[therigb])[plane];
552      Float_t maxene = (*qplmean[therigb+1])[plane];
553      //
554      if ( maxrig == minrig ){
555       printf("Unrecoverable ERROR! Matrix will be screwed up... \n");
556       return(0.);
557      };
558      Float_t b = log(maxene/minene)/(maxrig-minrig);
559      Float_t a = minene/exp(minrig*b);
560      Float_t ave = a*exp(b*rig);
561      //
562      return(ave);
563      //
564    }
565    Float_t CaloFranzini::GetFullAverageAt(Int_t plane, Int_t strip, Float_t rig){
566      //
567      Int_t therigb = 0;
568      for (Int_t i = 0; i<nbin-2; i++){
569        if ( rig>=brigm->At(i) && rig < brigm->At(i+1) ){
570          therigb = i;
571          break;
572        };
573      };
574      //
575      Float_t minrig;
576      Float_t maxrig;
577      //
578      //
579      if ( rig < brigm->At(0) ){
580        if ( rig < brig->At(0) ){
581    //      printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brigm->At(0));
582        };
583        therigb = 0;
584    };    };
585    if ( rig >= brig->At(nbin) ){    if ( rig >= brigm->At(nbin-2) ){
586      printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brig->At(nbin));      if ( rig >= brig->At(nbin-2) ) {
587      name = Form("qplmeann%i",nbin-1);  //      printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brigm->At(nbin-2));
588      printf(" Using qplmean %s \n",name.Data());      };
589        therigb = nbin - 3;
590    };    };
591    //    //
592    TArrayF *qplmean = (TArrayF*)file->Get(name.Data());    minrig = brigm->At(therigb);
593      maxrig = brigm->At(therigb+1);
594      //
595      Float_t minene = (*fqplmean[therigb])[plane][strip];
596      Float_t maxene = (*fqplmean[therigb+1])[plane][strip];
597      //
598      if ( maxrig == minrig ){
599       printf("Unrecoverable ERROR! Matrix will be screwed up... \n");
600       return(0.);
601      };
602      Float_t b = log(maxene/minene)/(maxrig-minrig);
603      Float_t a = minene/exp(minrig*b);
604      Float_t ave = a*exp(b*rig);
605    //    //
606    return(qplmean);    return(ave);
607    //    //
608  }  }
609    
610    Float_t CaloFranzini::GetHmatrixAt(Int_t iindex, Int_t jindex, Float_t rig){
611      Int_t therigb = 0;
612      for (Int_t i = 0; i<nbin-2; i++){
613        if ( rig>=brigm->At(i) && rig < brigm->At(i+1) ){
614          therigb = i;
615          break;
616        };
617      };
618      //
619      Float_t minrig;
620      Float_t maxrig;
621      //
622      if ( rig < brigm->At(0) ){
623        if ( rig < brig->At(0) ){
624    //      printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brigm->At(0));
625        };
626        therigb = 0;
627      };
628    //  if ( rig >= brigm->At(nbin-4) ){ // -2
629      if ( rig >= brigm->At(nbin-2) ){
630        if ( rig >= brig->At(nbin-2) ) {
631    //      printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brigm->At(nbin-2));
632        };
633    //    therigb = nbin - 5;// -3
634        therigb = nbin - 3;
635      };
636      //
637      if ( therigb < 2 ) therigb = 2;
638      minrig = brigm->At(therigb);
639      maxrig = brigm->At(therigb+1);
640      //  printf(" therigb %i minrig %f maxrig %f rig %f %i i %i j \n",therigb,minrig,maxrig,rig,iindex,jindex);
641      //
642      Float_t minene = (*hmat[therigb])[iindex][jindex];
643      Float_t maxene = (*hmat[therigb+1])[iindex][jindex];
644      //  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);
645      //
646    //  Float_t a = 0.;
647    //  Float_t b = 0.;
648    //  Float_t ave = 0.;
649    //  if ( minene == 0. ){
650    //      
651    //  } else {
652    //      b = log(maxene/minene)/(maxrig-minrig);
653    //      a = minene/exp(minrig*b);
654    //      ave = a*exp(b*rig);
655    //  };
656      //
657      Float_t m = (maxene-minene)/(maxrig-minrig);
658      Float_t q = minene - m * minrig;
659      Float_t ave = rig * m + q;
660      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);
661      //
662      //
663      return(ave);
664      //
665    }
666    
667    Float_t CaloFranzini::GetFullHmatrixAt(Int_t iindex, Int_t jindex, Float_t rig){
668      Int_t therigb = 0;
669      for (Int_t i = 0; i<nbin-2; i++){
670        if ( rig>=brigm->At(i) && rig < brigm->At(i+1) ){
671          therigb = i;
672          break;
673        };
674      };
675      //
676      Float_t minrig;
677      Float_t maxrig;
678      //
679      if ( rig < brigm->At(0) ){
680        if ( rig < brig->At(0) ){
681    //      printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brigm->At(0));
682        };
683        therigb = 0;
684      };
685    //  if ( rig >= brigm->At(nbin-4) ){ // -2
686      if ( rig >= brigm->At(nbin-2) ){
687        if ( rig >= brig->At(nbin-2) ) {
688    //      printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brigm->At(nbin-2));
689        };
690    //    therigb = nbin - 5;// -3
691        therigb = nbin - 3;
692      };
693      //
694      if ( therigb < 2 ) therigb = 2;
695      minrig = brigm->At(therigb);
696      maxrig = brigm->At(therigb+1);
697      //  printf(" therigb %i minrig %f maxrig %f rig %f %i i %i j \n",therigb,minrig,maxrig,rig,iindex,jindex);
698      //
699      Float_t minene = (*hfmat[therigb])[iindex][jindex];
700      Float_t maxene = (*hfmat[therigb+1])[iindex][jindex];
701      //  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);
702      //
703    //  Float_t a = 0.;
704    //  Float_t b = 0.;
705    //  Float_t ave = 0.;
706    //  if ( minene == 0. ){
707    //      
708    //  } else {
709    //      b = log(maxene/minene)/(maxrig-minrig);
710    //      a = minene/exp(minrig*b);
711    //      ave = a*exp(b*rig);
712    //  };
713      //
714      Float_t m = (maxene-minene)/(maxrig-minrig);
715      Float_t q = minene - m * minrig;
716      Float_t ave = rig * m + q;
717      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);
718      //
719      //
720      return(ave);
721      //
722    }
723    
724    void CaloFranzini::DrawLongAverage(Float_t rig){
725      //
726      TArrayF *ll = this->LoadLongAverage(rig);
727      //
728      gStyle->SetLabelSize(0.04);
729      gStyle->SetNdivisions(510,"XY");
730      //
731      TString hid = Form("cfralongvyvx");  
732      TCanvas *tcf  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
733      if ( tcf ){
734        tcf->cd();
735      } else {
736        tcf = new TCanvas(hid,hid);
737      };
738      //
739      TString thid = Form("hfralongvyvx");  
740      TH1F *thf  = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
741      if ( thf ) thf->Delete();
742      thf = new TH1F(thid,thid,44,-0.5,43.5);
743      tcf->cd();
744      Float_t qpl[43];
745      for (Int_t st=0;st<43;st++){
746        qpl[st] = ll->At(st);
747        printf("st %i qpl %f\n",st,qpl[st]);
748      };
749      for (Int_t st=0;st<44;st++){
750        if ( st == 37 ){
751          thf->Fill(st,0.);
752        } else {
753          Int_t ss = st;
754          if ( st > 37 ) ss--;
755          thf->Fill(st,qpl[ss]);
756        };
757      };
758      thf->Draw();
759      tcf->Modified();
760      tcf->Update();
761      //
762      gStyle->SetLabelSize(0);
763      gStyle->SetNdivisions(1,"XY");
764      //
765    };
766    
767    void CaloFranzini::DrawLongAverage(Int_t bin){
768      //
769      TArrayF *ll = this->LoadLongAverage(brigm->At(bin));
770      //
771      gStyle->SetLabelSize(0.04);
772      gStyle->SetNdivisions(510,"XY");
773      //
774      TString hid = Form("cfralongvyvx");  
775      TCanvas *tcf  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
776      if ( tcf ){
777        tcf->cd();
778      } else {
779        tcf = new TCanvas(hid,hid);
780      };
781      //
782      TString thid = Form("hfralongvyvx");  
783      TH1F *thf  = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
784      if ( thf ) thf->Delete();
785      thf = new TH1F(thid,thid,44,-0.5,43.5);
786      tcf->cd();
787      Float_t qpl[43];
788      for (Int_t st=0;st<43;st++){
789        qpl[st] = ll->At(st);
790        printf("st %i qpl %f\n",st,qpl[st]);
791      };
792      for (Int_t st=0;st<44;st++){
793        if ( st == 37 ){
794          thf->Fill(st,0.);
795        } else {
796          Int_t ss = st;
797          if ( st > 37 ) ss--;
798          thf->Fill(st,qpl[ss]);
799        };
800      };
801      thf->Draw();
802      tcf->Modified();
803      tcf->Update();
804      //
805      gStyle->SetLabelSize(0);
806      gStyle->SetNdivisions(1,"XY");
807      //
808    };
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.5
  ViewVC Help
Powered by ViewVC 1.1.23