/[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.7 by mocchiut, Mon Jan 21 10:24:11 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      file->Purge();
315    }
316    
317    void CaloFranzini::WriteFullNMean(TMatrixD *qpl, Int_t bin){
318      file->cd();
319      TString name = Form("fnqplmeann%i",bin);
320      file->WriteObject(&(*qpl),name.Data());
321      file->Purge();
322    }
323    
324    void CaloFranzini::WriteInvertedLongMatrix(TMatrixD mat, Int_t bin){
325      file->cd();
326    TString name = Form("matrixn%i",bin);    TString name = Form("matrixn%i",bin);
327      //  mat.Write(name.Data());
328      file->WriteObject(&mat,name.Data());
329    }
330    
331    void CaloFranzini::WriteInvertedFullMatrix(TMatrixD mat, Int_t bin){
332      file->cd();
333      TString name = Form("fmatrixn%i",bin);
334      //  mat.Write(name.Data());
335      file->WriteObject(&mat,name.Data());
336    }
337    
338    void CaloFranzini::WriteLongMatrix(TMatrixD *mat, Int_t bin){
339      file->cd();
340      TString name = Form("origmatrixn%i",bin);
341      //  mat.Write(name.Data());
342      file->WriteObject(&(*mat),name.Data());
343    }
344    
345    void CaloFranzini::WriteFullMatrix(TMatrixD *mat, Int_t bin){
346      file->cd();
347      TString name = Form("origfmatrixn%i",bin);
348      //  mat.Write(name.Data());
349      file->WriteObject(&(*mat),name.Data());
350      file->Purge();
351    }
352    
353    void CaloFranzini::WriteFullNMatrix(TMatrixF *mat, Int_t bin){
354      file->cd();
355      TString name = Form("origfnmatrixn%i",bin);
356      //  mat.Write(name.Data());
357    file->WriteObject(&(*mat),name.Data());    file->WriteObject(&(*mat),name.Data());
358      file->Purge();
359  }  }
360    
361  void CaloFranzini::CloseMatrixFile(){  void CaloFranzini::CloseMatrixFile(){
362      file->cd();
363    file->Close();    file->Close();
364  }  }
365    
# Line 235  Bool_t CaloFranzini::Open(TString matrix Line 380  Bool_t CaloFranzini::Open(TString matrix
380      return(false);      return(false);
381    };    };
382    //    //
383      if ( !this->LoadBin() ){
384        printf(" %s \n",matrixfile.Data());
385        return(false);
386      };
387      //
388      if ( dolong ){
389        if ( !this->LoadLong() ){
390          printf(" %s \n",matrixfile.Data());
391          return(false);
392        };
393        //
394        if ( !this->LoadMatrices() ){
395          printf(" %s \n",matrixfile.Data());
396          return(false);
397        };
398      };
399      //
400      if ( dofull ){
401        if ( !this->LoadFull() ){
402          printf(" %s \n",matrixfile.Data());
403          return(false);
404        };
405        //
406        if ( !this->LoadFullMatrices() ){
407          printf(" %s \n",matrixfile.Data());
408          return(false);
409        };
410      };
411      //
412      //
413      return(true);
414      //
415    }
416    
417    Bool_t CaloFranzini::LoadBin(){
418      //
419    TArrayI *numbin = (TArrayI*)file->Get("nbinenergy");    TArrayI *numbin = (TArrayI*)file->Get("nbinenergy");
420    if ( !numbin ){    if ( !numbin ){
421      printf(" ERROR: cannot read number of bins from file %s \n",matrixfile.Data());      printf(" ERROR: cannot read number of bins from file ");
422      return(false);      return(false);
423    };    };
424    nbin = (Int_t)numbin->At(0);    nbin = (Int_t)numbin->At(0);
425    if ( nbin <= 0 ){    if ( nbin <= 0 ){
426      printf(" ERROR: cannot work with 0 energy bins (from file %s) \n",matrixfile.Data());      printf(" ERROR: cannot work with 0 energy bins from file ");
427      return(false);      return(false);
428    };    };
429    //    //
430    brig = (TArrayF*)file->Get("binrig");    brig = (TArrayF*)file->Get("binrig");
431    if ( !brig ){    if ( !brig ){
432      printf(" ERROR: cannot read rigidity binning from file %s \n",matrixfile.Data());      printf(" ERROR: cannot read rigidity binning from file ");
433      return(false);      return(false);
434    };    };
435    //    //
436      brigm=(TArrayF*)file->Get("binrigmean");
437      if ( !brigm ){
438        printf(" ERROR: cannot read mean rigidity binning from file ");
439        return(false);
440      };
441      //
442      return(true);
443    }
444    
445    Bool_t CaloFranzini::LoadLong(){
446      //
447      for (Int_t i=0;i<17;i++){
448        TString name = Form("qplmeann%i",i);
449        qplmean[i] = (TArrayF*)file->Get(name.Data());
450        if ( !qplmean[i] ){
451          printf(" ERROR: cannot read average from file ");
452          return(false);
453        };
454      };
455      //
456      return(true);
457    }
458    
459    Bool_t CaloFranzini::LoadFull(){
460      //
461      for (Int_t i=0;i<17;i++){
462        TString name = Form("fqplmeann%i",i);
463        fqplmean[i] = (TMatrixD*)file->Get(name.Data());
464        if ( !fqplmean[i] ){
465          printf(" ERROR: cannot read average from file ");
466          return(false);
467        };
468      };
469      //
470    return(true);    return(true);
471    }
472    
473    Bool_t CaloFranzini::LoadMatrices(){
474      //
475      for (Int_t i=0;i<17;i++){
476            TString name1 = Form("matrixn%i",i);
477            hmat[i] = (TMatrixD*)file->Get(name1.Data());
478      };
479    //    //
480      return(true);
481    }
482    
483    Bool_t CaloFranzini::LoadFullMatrices(){
484      //
485      for (Int_t i=0;i<17;i++){
486            TString name1 = Form("fmatrixn%i",i);
487            hfmat[i] = (TMatrixF*)file->Get(name1.Data());
488      };
489      //
490      return(true);
491  }  }
492    
493  TMatrixD *CaloFranzini::LoadCovarianceMatrix(Float_t rig){  TMatrixD *CaloFranzini::LoadCovarianceMatrix(Float_t rig){
494    //    //
495    TString name;    Int_t mv = 0;
496    for (Int_t i = 0; i<nbin-1; i++){    for (Int_t i = 0; i<nbin-1; i++){
497      if ( rig>=brig->At(i) && rig < brig->At(i+1) ){      if ( rig>=brig->At(i) && rig < brig->At(i+1) ){
498        name = Form("matrixn%i",i);        mv = i;
499        break;        break;
500      };      };
501    };    };
502    if ( rig < brig->At(0) ){    if ( rig < brig->At(0) ){
503      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));
504      name = "matrixn0";      mv = 0;
     printf(" Using matrix %s \n",name.Data());  
505    };    };
506    if ( rig >= brig->At(nbin) ){    if ( rig >= brig->At(nbin-1) ){
507      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));
508      name = Form("matrixn%i",nbin-1);      mv = nbin-2;
     printf(" Using matrix %s \n",name.Data());  
509    };    };
510    //    //
511    TMatrixD *matrix = (TMatrixD*)file->Get(name.Data());    return(hmat[mv]);
512      //
513    }
514    
515    TMatrixD *CaloFranzini::LoadFullAverage(Int_t rigbin){
516      //
517      TString name = Form("fqplmeann%i",rigbin);
518      TMatrixD *fmean=(TMatrixD*)file->Get(name.Data());
519      //
520      return(fmean);
521      //
522    }
523    
524    TMatrixF *CaloFranzini::LoadFullMatrix(Int_t rigbin){
525      //
526      TString name = Form("origfmatrixn%i",rigbin);
527      TMatrixF *fmatri=(TMatrixF*)file->Get(name.Data());
528      //
529      return(fmatri);
530      //
531    }
532    
533    void CaloFranzini::LoadFullMatrix(Int_t rigbin, TMatrixF *&fmatri){
534      //
535      TString name = Form("origfmatrixn%i",rigbin);
536      fmatri=(TMatrixF*)file->Get(name.Data());
537      //
538    }
539    
540    void CaloFranzini::UnLoadFullAverage(Int_t rigbin){
541      //
542      TString name = Form("fqplmeann%i",rigbin);
543      file->Delete(name.Data());
544      //
545    }
546    
547    void CaloFranzini::UnLoadFullMatrix(Int_t rigbin){
548      //
549      TString name = Form("origfmatrixn%i",rigbin);
550      file->Delete(name.Data());
551      //
552    }
553    
554    TMatrixF *CaloFranzini::LoadFullNMatrix(Int_t rigbin){
555      //
556      TString name = Form("origfnmatrixn%i",rigbin);
557      TMatrixF *fnmatri=(TMatrixF*)file->Get(name.Data());
558      //
559      return(fnmatri);
560      //
561    }
562    
563    void CaloFranzini::UnLoadFullNMatrix(Int_t rigbin){
564      //
565      TString name = Form("origfnmatrixn%i",rigbin);
566      file->Delete(name.Data());
567      //
568    }
569    
570    TMatrixD *CaloFranzini::LoadFullNAverage(Int_t rigbin){
571      //
572      TString name = Form("fnqplmeann%i",rigbin);
573      TMatrixD *fnmean=(TMatrixD*)file->Get(name.Data());
574      //
575      return(fnmean);
576      //
577    }
578    
579    void CaloFranzini::UnLoadFullNAverage(Int_t rigbin){
580    //    //
581    return(matrix);    TString name = Form("fnqplmeann%i",rigbin);
582      file->Delete(name.Data());
583    //    //
584  }  }
585    
586    
587  TArrayF *CaloFranzini::LoadLongAverage(Float_t rig){  TArrayF *CaloFranzini::LoadLongAverage(Float_t rig){
588    //    //
589    TString name;    Int_t mv=0;
590    for (Int_t i = 0; i<nbin-1; i++){    for (Int_t i = 0; i<nbin-1; i++){
591      if ( rig>=brig->At(i) && rig < brig->At(i+1) ){      if ( rig>=brig->At(i) && rig < brig->At(i+1) ){
592        name = Form("qplmeann%i",i);          mv = i;
593        break;        break;
594      };      };
595    };    };
596    if ( rig < brig->At(0) ){    if ( rig < brig->At(0) ){
597      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));
598      name = "qplmeann0";      mv = 0;
599      printf(" Using qplmean %s \n",name.Data());    };
600      if ( rig >= brig->At(nbin-1) ){
601        printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brig->At(nbin-1));
602        mv=nbin-2;
603      };
604      //
605      return(qplmean[mv]);
606      //
607    }
608    
609    Float_t CaloFranzini::GetAverageAt(Int_t plane, Float_t rig){
610      //
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      //
623      if ( rig < brigm->At(0) ){
624        if ( rig < brig->At(0) ){
625    //      printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brigm->At(0));
626        };
627        therigb = 0;
628      };
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 - 3;
634      };
635      //
636      minrig = brigm->At(therigb);
637      maxrig = brigm->At(therigb+1);
638      //
639      Float_t minene = (*qplmean[therigb])[plane];
640      Float_t maxene = (*qplmean[therigb+1])[plane];
641      //
642      if ( maxrig == minrig ){
643       printf("Unrecoverable ERROR! Matrix will be screwed up... \n");
644       return(0.);
645      };
646      Float_t b = log(maxene/minene)/(maxrig-minrig);
647      Float_t a = minene/exp(minrig*b);
648      Float_t ave = a*exp(b*rig);
649      if ( a == 0. || minene == 0. || ave != ave ){
650        //  if ( a == 0. || minene == 0. ){
651        Float_t m = (maxene-minene)/(maxrig-minrig);
652        Float_t q = minene - m * minrig;
653        ave = rig * m + q;
654      };
655      //
656      return(ave);
657      //
658    }
659    
660    Float_t CaloFranzini::GetFullAverageAt(Int_t plane, Int_t strip, Float_t rig){
661      //
662      Int_t therigb = 0;
663      for (Int_t i = 0; i<nbin-2; i++){
664        if ( rig>=brigm->At(i) && rig < brigm->At(i+1) ){
665          therigb = i;
666          break;
667        };
668    };    };
669    if ( rig >= brig->At(nbin) ){    //
670      printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brig->At(nbin));    //
671      name = Form("qplmeann%i",nbin-1);    if ( rig < brigm->At(0) ){
672      printf(" Using qplmean %s \n",name.Data());      if ( rig < brig->At(0) ){
673    //      printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brigm->At(0));
674        };
675        therigb = 0;
676      };
677      if ( rig >= brigm->At(nbin-2) ){
678        if ( rig >= brig->At(nbin-2) ) {
679    //      printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brigm->At(nbin-2));
680        };
681        therigb = nbin - 3;
682      };
683      //
684      return(this->GetFullAverageAt(plane,strip,rig,therigb));
685      //
686    }
687    
688    Float_t CaloFranzini::GetFullAverageAt(Int_t plane, Int_t strip, Float_t rig, Int_t therigb){
689      //
690      Float_t minrig;
691      Float_t maxrig;
692      //
693      minrig = brigm->At(therigb);
694      maxrig = brigm->At(therigb+1);
695      //
696      Float_t minene = (*fqplmean[therigb])[plane][strip];
697      Float_t maxene = (*fqplmean[therigb+1])[plane][strip];
698      //
699      if ( maxrig == minrig ){
700       printf("Unrecoverable ERROR! Matrix will be screwed up... \n");
701       return(0.);
702      };
703      Float_t b = log(maxene/minene)/(maxrig-minrig);
704      Float_t a = minene/exp(minrig*b);
705      Float_t ave = a*exp(b*rig);
706      if ( a == 0. || minene == 0. || ave != ave ){
707        Float_t m = (maxene-minene)/(maxrig-minrig);
708        Float_t q = minene - m * minrig;
709        ave = rig * m + q;
710      };
711      //  
712      //  ave += (44.-plane)*strip;
713      //if ( a == 0. ) ave = 0.;
714      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);
715      //
716      return(ave);
717      //
718    }
719    
720    Float_t CaloFranzini::GetHmatrixAt(Int_t iindex, Int_t jindex, Float_t rig){
721      Int_t therigb = 0;
722      for (Int_t i = 0; i<nbin-2; i++){
723        if ( rig>=brigm->At(i) && rig < brigm->At(i+1) ){
724          therigb = i;
725          break;
726        };
727    };    };
728    //    //
729    TArrayF *qplmean = (TArrayF*)file->Get(name.Data());    Float_t minrig;
730      Float_t maxrig;
731      //
732      if ( rig < brigm->At(0) ){
733        if ( rig < brig->At(0) ){
734    //      printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brigm->At(0));
735        };
736        therigb = 0;
737      };
738    //  if ( rig >= brigm->At(nbin-4) ){ // -2
739      if ( rig >= brigm->At(nbin-2) ){
740        if ( rig >= brig->At(nbin-2) ) {
741    //      printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brigm->At(nbin-2));
742        };
743    //    therigb = nbin - 5;// -3
744        therigb = nbin - 3;
745      };
746      //
747      if ( therigb < 2 ) therigb = 2;
748      minrig = brigm->At(therigb);
749      maxrig = brigm->At(therigb+1);
750      //  printf(" therigb %i minrig %f maxrig %f rig %f %i i %i j \n",therigb,minrig,maxrig,rig,iindex,jindex);
751      //
752      Float_t minene = (*hmat[therigb])[iindex][jindex];
753      Float_t maxene = (*hmat[therigb+1])[iindex][jindex];
754      //  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);
755      //
756    //  Float_t a = 0.;
757    //  Float_t b = 0.;
758    //  Float_t ave = 0.;
759    //  if ( minene == 0. ){
760    //      
761    //  } else {
762    //      b = log(maxene/minene)/(maxrig-minrig);
763    //      a = minene/exp(minrig*b);
764    //      ave = a*exp(b*rig);
765    //  };
766      //
767      Float_t m = (maxene-minene)/(maxrig-minrig);
768      Float_t q = minene - m * minrig;
769      Float_t ave = rig * m + q;
770      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);
771      //
772      //
773      return(ave);
774      //
775    }
776    
777    Float_t CaloFranzini::GetFullHmatrixAt(Int_t iindex, Int_t jindex, Float_t rig){
778      Int_t therigb = 0;
779      for (Int_t i = 0; i<nbin-2; i++){
780        if ( rig>=brigm->At(i) && rig < brigm->At(i+1) ){
781          therigb = i;
782          break;
783        };
784      };
785      //
786      if ( rig < brigm->At(0) ){
787        if ( rig < brig->At(0) ){
788    //      printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brigm->At(0));
789        };
790        therigb = 0;
791      };
792    //  if ( rig >= brigm->At(nbin-4) ){ // -2
793      if ( rig >= brigm->At(nbin-2) ){
794        if ( rig >= brig->At(nbin-2) ) {
795    //      printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brigm->At(nbin-2));
796        };
797    //    therigb = nbin - 5;// -3
798        therigb = nbin - 3;
799      };
800      //
801      if ( therigb < 2 ) therigb = 2;
802      //
803      return(this->GetFullHmatrixAt(iindex,jindex,rig,therigb));
804      //
805    }
806    
807    Float_t CaloFranzini::GetFullHmatrixAt(Int_t iindex, Int_t jindex, Float_t rig, Int_t therigb){
808      //
809      Float_t minrig;
810      Float_t maxrig;
811      //
812      minrig = brigm->At(therigb);
813      maxrig = brigm->At(therigb+1);
814      //  printf(" therigb %i minrig %f maxrig %f rig %f %i i %i j \n",therigb,minrig,maxrig,rig,iindex,jindex);
815      //
816      Float_t minene = (*hfmat[therigb])[iindex][jindex];
817      Float_t maxene = (*hfmat[therigb+1])[iindex][jindex];
818      //  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);
819      //
820    //  Float_t a = 0.;
821    //  Float_t b = 0.;
822    //  Float_t ave = 0.;
823    //  if ( minene == 0. ){
824    //      
825    //  } else {
826    //      b = log(maxene/minene)/(maxrig-minrig);
827    //      a = minene/exp(minrig*b);
828    //      ave = a*exp(b*rig);
829    //  };
830      //
831      Float_t m = (maxene-minene)/(maxrig-minrig);
832      Float_t q = minene - m * minrig;
833      Float_t ave = rig * m + q;
834      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);
835      //
836      //
837      return(ave);
838      //
839    }
840    
841    void CaloFranzini::DrawLongAverage(Float_t rig){
842      //
843      TArrayF *ll = this->LoadLongAverage(rig);
844      //
845      gStyle->SetLabelSize(0.04);
846      gStyle->SetNdivisions(510,"XY");
847      //
848      TString hid = Form("cfralongvyvx");  
849      TCanvas *tcf  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
850      if ( tcf ){
851        tcf->cd();
852      } else {
853        tcf = new TCanvas(hid,hid);
854      };
855      //
856      TString thid = Form("hfralongvyvx");  
857      TH1F *thf  = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
858      if ( thf ) thf->Delete();
859      thf = new TH1F(thid,thid,44,-0.5,43.5);
860      tcf->cd();
861      Float_t qpl[43];
862      for (Int_t st=0;st<43;st++){
863        qpl[st] = ll->At(st);
864        printf("st %i qpl %f\n",st,qpl[st]);
865      };
866      for (Int_t st=0;st<44;st++){
867        if ( st == 37 ){
868          thf->Fill(st,0.);
869        } else {
870          Int_t ss = st;
871          if ( st > 37 ) ss--;
872          thf->Fill(st,qpl[ss]);
873        };
874      };
875      thf->Draw();
876      tcf->Modified();
877      tcf->Update();
878      //
879      gStyle->SetLabelSize(0);
880      gStyle->SetNdivisions(1,"XY");
881      //
882    };
883    
884    void CaloFranzini::DrawLongAverage(Int_t bin){
885      //
886      TArrayF *ll = this->LoadLongAverage(brigm->At(bin));
887      //
888      gStyle->SetLabelSize(0.04);
889      gStyle->SetNdivisions(510,"XY");
890      //
891      TString hid = Form("cfralongvyvx");  
892      TCanvas *tcf  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
893      if ( tcf ){
894        tcf->cd();
895      } else {
896        tcf = new TCanvas(hid,hid);
897      };
898      //
899      TString thid = Form("hfralongvyvx");  
900      TH1F *thf  = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
901      if ( thf ) thf->Delete();
902      thf = new TH1F(thid,thid,44,-0.5,43.5);
903      tcf->cd();
904      Float_t qpl[43];
905      for (Int_t st=0;st<43;st++){
906        qpl[st] = ll->At(st);
907        printf("st %i qpl %f\n",st,qpl[st]);
908      };
909      for (Int_t st=0;st<44;st++){
910        if ( st == 37 ){
911          thf->Fill(st,0.);
912        } else {
913          Int_t ss = st;
914          if ( st > 37 ) ss--;
915          thf->Fill(st,qpl[ss]);
916        };
917      };
918      thf->Draw();
919      tcf->Modified();
920      tcf->Update();
921      //
922      gStyle->SetLabelSize(0);
923      gStyle->SetNdivisions(1,"XY");
924      //
925    };
926    
927    Int_t CaloFranzini::ConvertStrip(Int_t mstrip){  
928      //
929      Int_t lastrip = 0;
930      //
931      if ( mstrip < 50 ) lastrip = 0;
932      //
933      if ( mstrip >= 50 && mstrip < 64 ) lastrip = 1;
934      //
935      if ( mstrip >= 64 && mstrip < 70 ) lastrip = 2;
936      //
937      if ( mstrip >= 70 && mstrip < 75 ) lastrip = 3;
938      //
939      if ( mstrip >= 75 && mstrip < 84 ){    
940        lastrip = (int)trunc((mstrip - 75)/3.) + 4;
941      };
942      //
943      if ( mstrip >= 84 && mstrip < 90 ){    
944        lastrip = (int)trunc((mstrip - 84)/2.) + 7;
945      };
946      //
947      if ( mstrip >= 90 && mstrip < 101 ){    
948        lastrip = mstrip - 90 + 10;
949      };
950      //
951      if ( mstrip >= 101 && mstrip < 107 ){    
952        lastrip = (int)trunc((mstrip - 101)/2.) + 21;
953      };
954      //
955      //
956      if ( mstrip >= 107 && mstrip < 116 ){    
957        lastrip = (int)trunc((mstrip - 107)/3.) + 24;
958      };
959      //
960      if ( mstrip >= 116 && mstrip < 121 ) lastrip = 27;
961      //
962      if ( mstrip >= 121 && mstrip < 127 ) lastrip = 28;
963      //
964      if ( mstrip >= 127 && mstrip < 141 ) lastrip = 29;
965    //    //
966    return(qplmean);    if ( mstrip >= 141 ) lastrip = 30;
967    //    //
968      return(lastrip);
969  }  }
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.7
  ViewVC Help
Powered by ViewVC 1.1.23