/[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.4 by mocchiut, Tue Dec 18 09:55:07 2007 UTC revision 1.5 by mocchiut, Thu Jan 3 10:02:28 2008 UTC
# Line 31  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 45  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 81  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 146  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 163  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 209  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 ){
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                  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));                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 226  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 292  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 299  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 306  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 332  Bool_t CaloFranzini::Open(TString matrix Line 367  Bool_t CaloFranzini::Open(TString matrix
367      printf(" %s \n",matrixfile.Data());      printf(" %s \n",matrixfile.Data());
368      return(false);      return(false);
369    };    };
370    if ( !this->LoadMatrices() ){    //
371      printf(" %s \n",matrixfile.Data());    if ( dolong ){
372      return(false);      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    //    //
# Line 367  Bool_t CaloFranzini::LoadBin(){ Line 422  Bool_t CaloFranzini::LoadBin(){
422      return(false);      return(false);
423    };    };
424    //    //
425      return(true);
426    }
427    
428    Bool_t CaloFranzini::LoadLong(){
429      //
430    for (Int_t i=0;i<17;i++){    for (Int_t i=0;i<17;i++){
431          TString name = Form("qplmeann%i",i);      TString name = Form("qplmeann%i",i);
432          qplmean[i] = (TArrayF*)file->Get(name.Data());      qplmean[i] = (TArrayF*)file->Get(name.Data());
433          if ( !qplmean[i] ){      if ( !qplmean[i] ){
434                  printf(" ERROR: cannot read average from file ");        printf(" ERROR: cannot read average from file ");
435                  return(false);        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);    return(true);
# Line 389  Bool_t CaloFranzini::LoadMatrices(){ Line 463  Bool_t CaloFranzini::LoadMatrices(){
463    return(true);    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    Int_t mv = 0;    Int_t mv = 0;
# Line 450  Float_t CaloFranzini::GetAverageAt(Int_t Line 534  Float_t CaloFranzini::GetAverageAt(Int_t
534    //    //
535    if ( rig < brigm->At(0) ){    if ( rig < brigm->At(0) ){
536      if ( rig < brig->At(0) ){      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));  //      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;      therigb = 0;
540    };    };
541    if ( rig >= brigm->At(nbin-2) ){    if ( rig >= brigm->At(nbin-2) ){
542      if ( rig >= brig->At(nbin-2) ) {      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));  //      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;      therigb = nbin - 3;
546    };    };
# Line 478  Float_t CaloFranzini::GetAverageAt(Int_t Line 562  Float_t CaloFranzini::GetAverageAt(Int_t
562    return(ave);    return(ave);
563    //    //
564  }  }
565    Float_t CaloFranzini::GetFullAverageAt(Int_t plane, Int_t strip, Float_t rig){
566  Float_t CaloFranzini::GetHmatrixAt(Int_t iindex, Int_t jindex, Float_t rig){    //
567    Int_t therigb = 0;    Int_t therigb = 0;
568    for (Int_t i = 0; i<nbin-2; i++){    for (Int_t i = 0; i<nbin-2; i++){
569      if ( rig>=brigm->At(i) && rig < brigm->At(i+1) ){      if ( rig>=brigm->At(i) && rig < brigm->At(i+1) ){
# Line 491  Float_t CaloFranzini::GetHmatrixAt(Int_t Line 575  Float_t CaloFranzini::GetHmatrixAt(Int_t
575    Float_t minrig;    Float_t minrig;
576    Float_t maxrig;    Float_t maxrig;
577    //    //
578      //
579    if ( rig < brigm->At(0) ){    if ( rig < brigm->At(0) ){
580      if ( rig < brig->At(0) ){      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));  //      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;      therigb = 0;
584    };    };
585    if ( rig >= brigm->At(nbin-2) ){    if ( rig >= brigm->At(nbin-2) ){
586      if ( rig >= brig->At(nbin-2) ) {      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));  //      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;      therigb = nbin - 3;
590    };    };
# Line 507  Float_t CaloFranzini::GetHmatrixAt(Int_t Line 592  Float_t CaloFranzini::GetHmatrixAt(Int_t
592    minrig = brigm->At(therigb);    minrig = brigm->At(therigb);
593    maxrig = brigm->At(therigb+1);    maxrig = brigm->At(therigb+1);
594    //    //
595    Float_t minene = (*hmat[therigb])[iindex][jindex];    Float_t minene = (*fqplmean[therigb])[plane][strip];
596    Float_t maxene = (*hmat[therigb+1])[iindex][jindex];    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);    Float_t b = log(maxene/minene)/(maxrig-minrig);
603    Float_t a = minene/exp(minrig*b);    Float_t a = minene/exp(minrig*b);
604    Float_t ave = a*exp(b*rig);    Float_t ave = a*exp(b*rig);
# Line 518  Float_t CaloFranzini::GetHmatrixAt(Int_t Line 607  Float_t CaloFranzini::GetHmatrixAt(Int_t
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    //    //
726    TArrayF *ll = this->LoadLongAverage(rig);    TArrayF *ll = this->LoadLongAverage(rig);
# Line 560  void CaloFranzini::DrawLongAverage(Float Line 763  void CaloFranzini::DrawLongAverage(Float
763    gStyle->SetNdivisions(1,"XY");    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.4  
changed lines
  Added in v.1.5
  ViewVC Help
Powered by ViewVC 1.1.23