/[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.3 by mocchiut, Thu Dec 13 17:08:11 2007 UTC
# Line 31  CaloFranzini::CaloFranzini(PamLevel2 *l2 Line 31  CaloFranzini::CaloFranzini(PamLevel2 *l2
31    dofull = false;    dofull = false;
32    sntr = 0;    sntr = 0;
33    //    //
34      sel = true;
35      cont = false;
36      N = 0;
37      NC = 43;
38      //
39      mask18b = -1;
40      //
41    Clear();    Clear();
42    //    //
43  }  }
# Line 44  void CaloFranzini::Clear(){ Line 51  void CaloFranzini::Clear(){
51    fulltzeta = -100.;    fulltzeta = -100.;
52    degfre = 0;    degfre = 0;
53    memset(estrip, 0, 2*22*96*sizeof(Float_t));    memset(estrip, 0, 2*22*96*sizeof(Float_t));
54    memset(qplane, 0, 2*22*sizeof(Float_t));    memset(qplane, 0, 43*sizeof(Float_t));
55    //    //
56  }  }
57    
# Line 55  void CaloFranzini::Print(){ Line 62  void CaloFranzini::Print(){
62    printf("========================================================================\n");    printf("========================================================================\n");
63    printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);    printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);
64    printf(" debug                [debug flag]:.. %i\n",debug);    printf(" debug                [debug flag]:.. %i\n",debug);
65      printf(" degree of freedom                :.. %i\n",this->GetDegreeOfFreedom());  
66    printf(" longtzeta                        :.. %f\n",longtzeta);      printf(" longtzeta                        :.. %f\n",longtzeta);  
67    printf(" fulltzeta                        :.. %f\n",fulltzeta);      printf(" longtzeta normalized             :.. %f\n",this->GetNormLongTZeta());  
68      //  printf(" fulltzeta                        :.. %f\n",fulltzeta);  
69      //  printf(" longtzeta normalized             :.. %f\n",this->GetNormFullTZeta());  
70    printf("========================================================================\n");    printf("========================================================================\n");
71    //    //
72  }  }
# Line 69  void CaloFranzini::Delete(){ Line 79  void CaloFranzini::Delete(){
79    //    //
80  }  }
81    
82    void CaloFranzini::SetNoWpreSampler(Int_t n){
83      if ( NC+n < 23 ){
84        N = n;
85      } else {
86        printf(" ERROR! Calorimeter is made of 22 W planes\n");
87        printf(" you are giving N presampler = %i and N calo = %i \n",n,NC);
88        printf(" WARNING: using default values NWpre = 0, NWcalo = 22\n");
89        NC = 43;
90        N = 0;
91      };
92    }
93    
94    void CaloFranzini::SetNoWcalo(Int_t n){
95      if ( N+n < 23 ){
96        NC = n*2;
97        if ( NC >37 ) NC--;
98      } else {
99        printf(" ERROR! Calorimeter is made of 22 W planes\n");
100        printf(" you are giving N W presampler = %i and N W calo = %i \n",N,n);
101        printf(" WARNING: using default values NWpre = 0, NWcalo = 22\n");
102        NC = 43;
103        N = 0;
104      };
105    }
106    
107  void CaloFranzini::Process(){  void CaloFranzini::Process(){
108    this->Process(0);    this->Process(0);
# Line 124  void CaloFranzini::Process(Int_t itr){ Line 158  void CaloFranzini::Process(Int_t itr){
158      mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);      mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
159      //      //
160      estrip[view][plane][strip] = mip;      estrip[view][plane][strip] = mip;
161      nplane = 1 - view + 2 * plane;      //    
162      qplane[nplane] += mip;      nplane = 1 - view + 2 * (plane - N);
163        if ( nplane > (37-(2*N)) ) nplane--;
164        //
165        if ( plane == (18+N) ) mip = 0.;
166        if ( nplane > -1 ) qplane[nplane] += mip;
167      //      //
168    };    };
169    //    //
# Line 133  void CaloFranzini::Process(Int_t itr){ Line 171  void CaloFranzini::Process(Int_t itr){
171    //    //
172    if ( dolong ){    if ( dolong ){
173      //      //
174      degfre = 44;      if ( cont ){
175          for (Int_t i=0; i<22; i++){
176            if ( i == (18+N) ){
177              mask18b =  1 + 2 * (i - N);
178              break;
179            };
180          };
181        };  
182        //  
183        if ( sel ){
184          for (Int_t i=0; i<22; i++){
185            if ( i == (18-N) ){
186              mask18b =  1 + 2 * (i - N);
187              break;
188            };
189          };
190        };
191        //
192        if ( mask18b == 37 ) mask18b = -1;
193        //
194        Int_t dgf = 43;
195      //      //
196      for (Int_t i=0; i < 22; i++){      for (Int_t i=0; i < 22; i++){
197        if ( L2->GetTrack(itr)->GetCaloTrack()->tibar[i][1] < 0 ){        if ( L2->GetTrack(itr)->GetCaloTrack()->tibar[i][1] < 0 ){
198          degfre = 2 * i;          dgf = 2 * i;
199          break;          break;
200        };        };
201        if ( L2->GetTrack(itr)->GetCaloTrack()->tibar[i][0] < 0 ){        if ( L2->GetTrack(itr)->GetCaloTrack()->tibar[i][0] < 0 ){
202        degfre = 1 + 2 * i;        dgf = 1 + 2 * i;
203        break;        break;
204        };        };
205      };      };
206      //      //
207        if ( dgf < 43 && dgf > 37 ) dgf--;
208        //
209        degfre = TMath::Min(dgf,NC);
210        //
211      if ( degfre > 0 ){      if ( degfre > 0 ){
212        TArrayF *qplmean = this->LoadLongAverage(rig);        TArrayF *qplmean = this->LoadLongAverage(rig);
213        TMatrixD *lmatrix = this->LoadCovarianceMatrix(rig);        TMatrixD *lmatrix = this->LoadCovarianceMatrix(rig);
214        for (Int_t i = 0; i < degfre; i++){        for (Int_t i = 0; i < degfre; i++){
215          for (Int_t j = 0; j < degfre; j++){              if ( i != mask18b ){
216            longtzeta += (qplane[i]-qplmean->At(i)) * (*lmatrix)[i][j] * (qplane[j]-qplmean->At(j));            for (Int_t j = 0; j < degfre; j++){  
217                if ( j != mask18b ){
218                  longtzeta += (qplane[i]-qplmean->At(i)) * (*lmatrix)[i][j] * (qplane[j]-qplmean->At(j));
219                };
220              };
221          };          };
222        };        };
223      };      };
# Line 186  Bool_t CaloFranzini::CreateMatrixFile(TS Line 252  Bool_t CaloFranzini::CreateMatrixFile(TS
252    //    //
253    if ( !file || file->IsZombie() ){    if ( !file || file->IsZombie() ){
254      file = new TFile(matrixfile.Data(),"RECREATE");      file = new TFile(matrixfile.Data(),"RECREATE");
255        printf(" Create file %s \n",file->GetName());
256    } else {    } else {
257      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());
258      return(false);      return(false);
# Line 196  Bool_t CaloFranzini::CreateMatrixFile(TS Line 263  Bool_t CaloFranzini::CreateMatrixFile(TS
263  }  }
264    
265  void CaloFranzini::WriteNumBin(Int_t numbin){  void CaloFranzini::WriteNumBin(Int_t numbin){
266      file->cd();
267    TArrayI nbi(1, &numbin);    TArrayI nbi(1, &numbin);
268    file->WriteObject(&nbi, "nbinenergy");    file->WriteObject(&nbi, "nbinenergy");
269  }  }
270    
271  void CaloFranzini::WriteRigBin(TArrayF *rigbin){  void CaloFranzini::WriteRigBin(TArrayF *rigbin){
272    file->WriteObject(&rigbin, "binrig");    file->cd();
273      //  rigbin->Write("binrig");
274      file->WriteObject(&(*rigbin), "binrig");
275  }  }
276    
277  void CaloFranzini::WriteLongMean(TArrayF *qpl, Int_t bin){  void CaloFranzini::WriteLongMean(TArrayF *qpl, Int_t bin){
278      file->cd();
279    TString name = Form("qplmeann%i",bin);    TString name = Form("qplmeann%i",bin);
280    file->WriteObject(&qpl,name.Data());    file->WriteObject(&(*qpl),name.Data());
281  }  }
282    
283  void CaloFranzini::WriteLongMatrix(TMatrixD *mat, Int_t bin){  void CaloFranzini::WriteInvertedLongMatrix(TMatrixD mat, Int_t bin){
284      file->cd();
285    TString name = Form("matrixn%i",bin);    TString name = Form("matrixn%i",bin);
286      //  mat.Write(name.Data());
287      file->WriteObject(&mat,name.Data());
288    }
289    
290    void CaloFranzini::WriteLongMatrix(TMatrixD *mat, Int_t bin){
291      file->cd();
292      TString name = Form("origmatrixn%i",bin);
293      //  mat.Write(name.Data());
294    file->WriteObject(&(*mat),name.Data());    file->WriteObject(&(*mat),name.Data());
295  }  }
296    
297  void CaloFranzini::CloseMatrixFile(){  void CaloFranzini::CloseMatrixFile(){
298      file->cd();
299    file->Close();    file->Close();
300  }  }
301    
# Line 270  TMatrixD *CaloFranzini::LoadCovarianceMa Line 351  TMatrixD *CaloFranzini::LoadCovarianceMa
351      name = "matrixn0";      name = "matrixn0";
352      printf(" Using matrix %s \n",name.Data());      printf(" Using matrix %s \n",name.Data());
353    };    };
354    if ( rig >= brig->At(nbin) ){    if ( rig >= brig->At(nbin-1) ){
355      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));
356      name = Form("matrixn%i",nbin-1);      name = Form("matrixn%i",nbin-2);
357      printf(" Using matrix %s \n",name.Data());      printf(" Using matrix %s \n",name.Data());
358    };    };
359    //    //
# Line 297  TArrayF *CaloFranzini::LoadLongAverage(F Line 378  TArrayF *CaloFranzini::LoadLongAverage(F
378      name = "qplmeann0";      name = "qplmeann0";
379      printf(" Using qplmean %s \n",name.Data());      printf(" Using qplmean %s \n",name.Data());
380    };    };
381    if ( rig >= brig->At(nbin) ){    if ( rig >= brig->At(nbin-1) ){
382      printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brig->At(nbin));      printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brig->At(nbin-1));
383      name = Form("qplmeann%i",nbin-1);      name = Form("qplmeann%i",nbin-2);
384      printf(" Using qplmean %s \n",name.Data());      printf(" Using qplmean %s \n",name.Data());
385    };    };
386    //    //
# Line 308  TArrayF *CaloFranzini::LoadLongAverage(F Line 389  TArrayF *CaloFranzini::LoadLongAverage(F
389    return(qplmean);    return(qplmean);
390    //    //
391  }  }
392    
393    
394    
395    void CaloFranzini::DrawLongAverage(Float_t rig){
396      //
397      TArrayF *ll = this->LoadLongAverage(rig);
398      //
399      gStyle->SetLabelSize(0.04);
400      gStyle->SetNdivisions(510,"XY");
401      //
402      TString hid = Form("cfralongvyvx");  
403      TCanvas *tcf  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
404      if ( tcf ){
405        tcf->cd();
406      } else {
407        tcf = new TCanvas(hid,hid);
408      };
409      //
410      TString thid = Form("hfralongvyvx");  
411      TH1F *thf  = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
412      if ( thf ) thf->Delete();
413      thf = new TH1F(thid,thid,44,-0.5,43.5);
414      tcf->cd();
415      //  Int_t pp=0;
416      Float_t qpl[43];
417      for (Int_t st=0;st<43;st++){
418        qpl[st] = ll->At(st);
419      };
420      for (Int_t st=0;st<44;st++){
421        if ( st == 37 ){
422          thf->Fill(st,0.);
423        } else {
424          Int_t ss = st;
425          if ( st > 37 ) ss--;
426          thf->Fill(st,qpl[ss]);
427        };
428      };
429      thf->Draw();
430      tcf->Modified();
431      tcf->Update();
432      //
433      gStyle->SetLabelSize(0);
434      gStyle->SetNdivisions(1,"XY");
435      //
436    };
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.23