/[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.3 by mocchiut, Thu Dec 13 17:08:11 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;    sel = true;
40    cont = false;    cont = false;
41    N = 0;    N = 0;
# Line 44  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, 43*sizeof(Float_t));    memset(qplane, 0, 43*sizeof(Float_t));
# Line 80  void CaloFranzini::Delete(){ Line 82  void CaloFranzini::Delete(){
82  }  }
83    
84  void CaloFranzini::SetNoWpreSampler(Int_t n){  void CaloFranzini::SetNoWpreSampler(Int_t n){
85    if ( NC+n < 23 ){    Int_t nc2 = NC/2;
86      if ( NC >= 37 ) nc2 = (NC+1)/2;
87      if ( nc2+n < 23 ){
88      N = n;      N = n;
89    } else {    } else {
90      printf(" ERROR! Calorimeter is made of 22 W planes\n");      printf(" ERROR! Calorimeter is made of 22 W planes\n");
# Line 145  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 162  void CaloFranzini::Process(Int_t itr){ Line 167  void CaloFranzini::Process(Int_t itr){
167      nplane = 1 - view + 2 * (plane - N);      nplane = 1 - view + 2 * (plane - N);
168      if ( nplane > (37-(2*N)) ) nplane--;      if ( nplane > (37-(2*N)) ) nplane--;
169      //      //
170      if ( plane == (18+N) ) mip = 0.;  //    if ( plane == (18+N) ) mip = 0.;
171      if ( nplane > -1 ) qplane[nplane] += mip;      if ( nplane > -1 ) qplane[nplane] += mip;
172      //      //
173    };    };
# Line 208  void CaloFranzini::Process(Int_t itr){ Line 213  void CaloFranzini::Process(Int_t itr){
213      //      //
214      degfre = TMath::Min(dgf,NC);      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          if ( i != mask18b ){          if ( i != mask18b ){
222            for (Int_t j = 0; j < degfre; j++){              for (Int_t j = 0; j < degfre; j++){  
223              if ( j != mask18b ){              if ( j != mask18b ){
224                longtzeta += (qplane[i]-qplmean->At(i)) * (*lmatrix)[i][j] * (qplane[j]-qplmean->At(j));                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 227  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 262  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();    file->cd();
294    TArrayI nbi(1, &numbin);    TArrayI nbi(1, &numbin);
# Line 280  void CaloFranzini::WriteLongMean(TArrayF Line 307  void CaloFranzini::WriteLongMean(TArrayF
307    file->WriteObject(&(*qpl),name.Data());    file->WriteObject(&(*qpl),name.Data());
308  }  }
309    
310    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){  void CaloFranzini::WriteInvertedLongMatrix(TMatrixD mat, Int_t bin){
317    file->cd();    file->cd();
318    TString name = Form("matrixn%i",bin);    TString name = Form("matrixn%i",bin);
# Line 287  void CaloFranzini::WriteInvertedLongMatr Line 320  void CaloFranzini::WriteInvertedLongMatr
320    file->WriteObject(&mat,name.Data());    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){  void CaloFranzini::WriteLongMatrix(TMatrixD *mat, Int_t bin){
331    file->cd();    file->cd();
332    TString name = Form("origmatrixn%i",bin);    TString name = Form("origmatrixn%i",bin);
# Line 294  void CaloFranzini::WriteLongMatrix(TMatr Line 334  void CaloFranzini::WriteLongMatrix(TMatr
334    file->WriteObject(&(*mat),name.Data());    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());
342    }
343    
344  void CaloFranzini::CloseMatrixFile(){  void CaloFranzini::CloseMatrixFile(){
345    file->cd();    file->cd();
346    file->Close();    file->Close();
# Line 316  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-1) ){    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-1));      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-2);      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;
     printf(" Using qplmean %s \n",name.Data());  
511    };    };
512    if ( rig >= brig->At(nbin-1) ){    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));      printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brig->At(nbin-1));
514      name = Form("qplmeann%i",nbin-2);      mv=nbin-2;
     printf(" Using qplmean %s \n",name.Data());  
515    };    };
516    //    //
517    TArrayF *qplmean = (TArrayF*)file->Get(name.Data());    return(qplmean[mv]);
518    //    //
519    return(qplmean);  }
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 >= brigm->At(nbin-2) ){
586        if ( rig >= brig->At(nbin-2) ) {
587    //      printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brigm->At(nbin-2));
588        };
589        therigb = nbin - 3;
590      };
591      //
592      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(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){  void CaloFranzini::DrawLongAverage(Float_t rig){
725    //    //
# Line 412  void CaloFranzini::DrawLongAverage(Float Line 741  void CaloFranzini::DrawLongAverage(Float
741    if ( thf ) thf->Delete();    if ( thf ) thf->Delete();
742    thf = new TH1F(thid,thid,44,-0.5,43.5);    thf = new TH1F(thid,thid,44,-0.5,43.5);
743    tcf->cd();    tcf->cd();
   //  Int_t pp=0;  
744    Float_t qpl[43];    Float_t qpl[43];
745    for (Int_t st=0;st<43;st++){    for (Int_t st=0;st<43;st++){
746      qpl[st] = ll->At(st);      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++){    for (Int_t st=0;st<44;st++){
793      if ( st == 37 ){      if ( st == 37 ){
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.5
  ViewVC Help
Powered by ViewVC 1.1.23