/[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.6 by mocchiut, Fri Jan 11 15:27:13 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      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){  void CaloFranzini::WriteInvertedLongMatrix(TMatrixD mat, Int_t bin){
325    file->cd();    file->cd();
326    TString name = Form("matrixn%i",bin);    TString name = Form("matrixn%i",bin);
# Line 287  void CaloFranzini::WriteInvertedLongMatr Line 328  void CaloFranzini::WriteInvertedLongMatr
328    file->WriteObject(&mat,name.Data());    file->WriteObject(&mat,name.Data());
329  }  }
330    
331    void CaloFranzini::WriteInvertedFullMatrix(TMatrixF 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){  void CaloFranzini::WriteLongMatrix(TMatrixD *mat, Int_t bin){
339    file->cd();    file->cd();
340    TString name = Form("origmatrixn%i",bin);    TString name = Form("origmatrixn%i",bin);
# Line 294  void CaloFranzini::WriteLongMatrix(TMatr Line 342  void CaloFranzini::WriteLongMatrix(TMatr
342    file->WriteObject(&(*mat),name.Data());    file->WriteObject(&(*mat),name.Data());
343  }  }
344    
345    void CaloFranzini::WriteFullMatrix(TMatrixF *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());
358      file->Purge();
359    }
360    
361  void CaloFranzini::CloseMatrixFile(){  void CaloFranzini::CloseMatrixFile(){
362    file->cd();    file->cd();
363    file->Close();    file->Close();
# Line 316  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);
434      };
435      //
436      brigm=(TArrayF*)file->Get("binrigmean");
437      if ( !brigm ){
438        printf(" ERROR: cannot read mean rigidity binning from file ");
439      return(false);      return(false);
440    };    };
441    //    //
442    return(true);    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);
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-1) ){    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-1));      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-2);      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    return(matrix);  }
578    
579    void CaloFranzini::UnLoadFullNAverage(Int_t rigbin){
580      //
581      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;
     printf(" Using qplmean %s \n",name.Data());  
599    };    };
600    if ( rig >= brig->At(nbin-1) ){    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));      printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brig->At(nbin-1));
602      name = Form("qplmeann%i",nbin-2);      mv=nbin-2;
603      printf(" Using qplmean %s \n",name.Data());    };
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      //
670      //
671      if ( rig < brigm->At(0) ){
672        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      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    TArrayF *qplmean = (TArrayF*)file->Get(name.Data());  //  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    return(qplmean);    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){  void CaloFranzini::DrawLongAverage(Float_t rig){
842    //    //
# Line 412  void CaloFranzini::DrawLongAverage(Float Line 858  void CaloFranzini::DrawLongAverage(Float
858    if ( thf ) thf->Delete();    if ( thf ) thf->Delete();
859    thf = new TH1F(thid,thid,44,-0.5,43.5);    thf = new TH1F(thid,thid,44,-0.5,43.5);
860    tcf->cd();    tcf->cd();
   //  Int_t pp=0;  
861    Float_t qpl[43];    Float_t qpl[43];
862    for (Int_t st=0;st<43;st++){    for (Int_t st=0;st<43;st++){
863      qpl[st] = ll->At(st);      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++){    for (Int_t st=0;st<44;st++){
910      if ( st == 37 ){      if ( st == 37 ){
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.6
  ViewVC Help
Powered by ViewVC 1.1.23