/[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.4 by mocchiut, Tue Dec 18 09:55:07 2007 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 31  CaloFranzini::CaloFranzini(PamLevel2 *l2 Line 32  CaloFranzini::CaloFranzini(PamLevel2 *l2
32    dofull = false;    dofull = false;
33    sntr = 0;    sntr = 0;
34    //    //
35      sel = true;
36      cont = false;
37      N = 0;
38      NC = 43;
39      //
40      mask18b = -1;
41      //
42    Clear();    Clear();
43    //    //
44  }  }
# Line 44  void CaloFranzini::Clear(){ Line 52  void CaloFranzini::Clear(){
52    fulltzeta = -100.;    fulltzeta = -100.;
53    degfre = 0;    degfre = 0;
54    memset(estrip, 0, 2*22*96*sizeof(Float_t));    memset(estrip, 0, 2*22*96*sizeof(Float_t));
55    memset(qplane, 0, 2*22*sizeof(Float_t));    memset(qplane, 0, 43*sizeof(Float_t));
56    //    //
57  }  }
58    
# Line 55  void CaloFranzini::Print(){ Line 63  void CaloFranzini::Print(){
63    printf("========================================================================\n");    printf("========================================================================\n");
64    printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);    printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);
65    printf(" debug                [debug flag]:.. %i\n",debug);    printf(" debug                [debug flag]:.. %i\n",debug);
66      printf(" degree of freedom                :.. %i\n",this->GetDegreeOfFreedom());  
67    printf(" longtzeta                        :.. %f\n",longtzeta);      printf(" longtzeta                        :.. %f\n",longtzeta);  
68    printf(" fulltzeta                        :.. %f\n",fulltzeta);      printf(" longtzeta normalized             :.. %f\n",this->GetNormLongTZeta());  
69      //  printf(" fulltzeta                        :.. %f\n",fulltzeta);  
70      //  printf(" longtzeta normalized             :.. %f\n",this->GetNormFullTZeta());  
71    printf("========================================================================\n");    printf("========================================================================\n");
72    //    //
73  }  }
# Line 69  void CaloFranzini::Delete(){ Line 80  void CaloFranzini::Delete(){
80    //    //
81  }  }
82    
83    void CaloFranzini::SetNoWpreSampler(Int_t n){
84      if ( NC+n < 23 ){
85        N = n;
86      } else {
87        printf(" ERROR! Calorimeter is made of 22 W planes\n");
88        printf(" you are giving N presampler = %i and N calo = %i \n",n,NC);
89        printf(" WARNING: using default values NWpre = 0, NWcalo = 22\n");
90        NC = 43;
91        N = 0;
92      };
93    }
94    
95    void CaloFranzini::SetNoWcalo(Int_t n){
96      if ( N+n < 23 ){
97        NC = n*2;
98        if ( NC >37 ) NC--;
99      } else {
100        printf(" ERROR! Calorimeter is made of 22 W planes\n");
101        printf(" you are giving N W presampler = %i and N W calo = %i \n",N,n);
102        printf(" WARNING: using default values NWpre = 0, NWcalo = 22\n");
103        NC = 43;
104        N = 0;
105      };
106    }
107    
108  void CaloFranzini::Process(){  void CaloFranzini::Process(){
109    this->Process(0);    this->Process(0);
# Line 124  void CaloFranzini::Process(Int_t itr){ Line 159  void CaloFranzini::Process(Int_t itr){
159      mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);      mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
160      //      //
161      estrip[view][plane][strip] = mip;      estrip[view][plane][strip] = mip;
162      nplane = 1 - view + 2 * plane;      //    
163      qplane[nplane] += mip;      nplane = 1 - view + 2 * (plane - N);
164        if ( nplane > (37-(2*N)) ) nplane--;
165        //
166        if ( plane == (18+N) ) mip = 0.;
167        if ( nplane > -1 ) qplane[nplane] += mip;
168      //      //
169    };    };
170    //    //
# Line 133  void CaloFranzini::Process(Int_t itr){ Line 172  void CaloFranzini::Process(Int_t itr){
172    //    //
173    if ( dolong ){    if ( dolong ){
174      //      //
175      degfre = 44;      if ( cont ){
176          for (Int_t i=0; i<22; i++){
177            if ( i == (18+N) ){
178              mask18b =  1 + 2 * (i - N);
179              break;
180            };
181          };
182        };  
183        //  
184        if ( sel ){
185          for (Int_t i=0; i<22; i++){
186            if ( i == (18-N) ){
187              mask18b =  1 + 2 * (i - N);
188              break;
189            };
190          };
191        };
192        //
193        if ( mask18b == 37 ) mask18b = -1;
194        //
195        Int_t dgf = 43;
196      //      //
197      for (Int_t i=0; i < 22; i++){      for (Int_t i=0; i < 22; i++){
198        if ( L2->GetTrack(itr)->GetCaloTrack()->tibar[i][1] < 0 ){        if ( L2->GetTrack(itr)->GetCaloTrack()->tibar[i][1] < 0 ){
199          degfre = 2 * i;          dgf = 2 * i;
200          break;          break;
201        };        };
202        if ( L2->GetTrack(itr)->GetCaloTrack()->tibar[i][0] < 0 ){        if ( L2->GetTrack(itr)->GetCaloTrack()->tibar[i][0] < 0 ){
203        degfre = 1 + 2 * i;        dgf = 1 + 2 * i;
204        break;        break;
205        };        };
206      };      };
207      //      //
208        if ( dgf < 43 && dgf > 37 ) dgf--;
209        //
210        degfre = TMath::Min(dgf,NC);
211        //
212      if ( degfre > 0 ){      if ( degfre > 0 ){
       TArrayF *qplmean = this->LoadLongAverage(rig);  
       TMatrixD *lmatrix = this->LoadCovarianceMatrix(rig);  
213        for (Int_t i = 0; i < degfre; i++){        for (Int_t i = 0; i < degfre; i++){
214          for (Int_t j = 0; j < degfre; j++){              if ( i != mask18b ){
215            longtzeta += (qplane[i]-qplmean->At(i)) * (*lmatrix)[i][j] * (qplane[j]-qplmean->At(j));            for (Int_t j = 0; j < degfre; j++){  
216                if ( j != mask18b ){
217                  longtzeta += (qplane[i]-this->GetAverageAt(i,rig)) * this->GetHmatrixAt(i,j,rig) * (qplane[j]-this->GetAverageAt(j,rig));
218                };
219              };
220          };          };
221        };        };
222      };      };
# Line 186  Bool_t CaloFranzini::CreateMatrixFile(TS Line 251  Bool_t CaloFranzini::CreateMatrixFile(TS
251    //    //
252    if ( !file || file->IsZombie() ){    if ( !file || file->IsZombie() ){
253      file = new TFile(matrixfile.Data(),"RECREATE");      file = new TFile(matrixfile.Data(),"RECREATE");
254        printf(" Create file %s \n",file->GetName());
255    } else {    } else {
256      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());
257      return(false);      return(false);
# Line 195  Bool_t CaloFranzini::CreateMatrixFile(TS Line 261  Bool_t CaloFranzini::CreateMatrixFile(TS
261    //    //
262  }  }
263    
264    Bool_t CaloFranzini::UpdateMatrixFile(TString matrixfile){  
265      //
266      file = new TFile(matrixfile.Data(),"UPDATE");
267      //
268      if ( !file || file->IsZombie() ){
269        printf(" ERROR: file %s already existing!\n Choose another name or delete the old file\n",matrixfile.Data());
270        return(false);
271      };
272      //
273      return(true);
274      //
275    }
276    
277  void CaloFranzini::WriteNumBin(Int_t numbin){  void CaloFranzini::WriteNumBin(Int_t numbin){
278      file->cd();
279    TArrayI nbi(1, &numbin);    TArrayI nbi(1, &numbin);
280    file->WriteObject(&nbi, "nbinenergy");    file->WriteObject(&nbi, "nbinenergy");
281  }  }
282    
283  void CaloFranzini::WriteRigBin(TArrayF *rigbin){  void CaloFranzini::WriteRigBin(TArrayF *rigbin){
284    file->WriteObject(&rigbin, "binrig");    file->cd();
285      //  rigbin->Write("binrig");
286      file->WriteObject(&(*rigbin), "binrig");
287  }  }
288    
289  void CaloFranzini::WriteLongMean(TArrayF *qpl, Int_t bin){  void CaloFranzini::WriteLongMean(TArrayF *qpl, Int_t bin){
290      file->cd();
291    TString name = Form("qplmeann%i",bin);    TString name = Form("qplmeann%i",bin);
292    file->WriteObject(&qpl,name.Data());    file->WriteObject(&(*qpl),name.Data());
293  }  }
294    
295  void CaloFranzini::WriteLongMatrix(TMatrixD *mat, Int_t bin){  void CaloFranzini::WriteInvertedLongMatrix(TMatrixD mat, Int_t bin){
296      file->cd();
297    TString name = Form("matrixn%i",bin);    TString name = Form("matrixn%i",bin);
298      //  mat.Write(name.Data());
299      file->WriteObject(&mat,name.Data());
300    }
301    
302    void CaloFranzini::WriteLongMatrix(TMatrixD *mat, Int_t bin){
303      file->cd();
304      TString name = Form("origmatrixn%i",bin);
305      //  mat.Write(name.Data());
306    file->WriteObject(&(*mat),name.Data());    file->WriteObject(&(*mat),name.Data());
307  }  }
308    
309  void CaloFranzini::CloseMatrixFile(){  void CaloFranzini::CloseMatrixFile(){
310      file->cd();
311    file->Close();    file->Close();
312  }  }
313    
# Line 235  Bool_t CaloFranzini::Open(TString matrix Line 328  Bool_t CaloFranzini::Open(TString matrix
328      return(false);      return(false);
329    };    };
330    //    //
331      if ( !this->LoadBin() ){
332        printf(" %s \n",matrixfile.Data());
333        return(false);
334      };
335      if ( !this->LoadMatrices() ){
336        printf(" %s \n",matrixfile.Data());
337        return(false);
338      };
339      //
340      //
341      return(true);
342      //
343    }
344    
345    Bool_t CaloFranzini::LoadBin(){
346      //
347    TArrayI *numbin = (TArrayI*)file->Get("nbinenergy");    TArrayI *numbin = (TArrayI*)file->Get("nbinenergy");
348    if ( !numbin ){    if ( !numbin ){
349      printf(" ERROR: cannot read number of bins from file %s \n",matrixfile.Data());      printf(" ERROR: cannot read number of bins from file ");
350      return(false);      return(false);
351    };    };
352    nbin = (Int_t)numbin->At(0);    nbin = (Int_t)numbin->At(0);
353    if ( nbin <= 0 ){    if ( nbin <= 0 ){
354      printf(" ERROR: cannot work with 0 energy bins (from file %s) \n",matrixfile.Data());      printf(" ERROR: cannot work with 0 energy bins from file ");
355      return(false);      return(false);
356    };    };
357    //    //
358    brig = (TArrayF*)file->Get("binrig");    brig = (TArrayF*)file->Get("binrig");
359    if ( !brig ){    if ( !brig ){
360      printf(" ERROR: cannot read rigidity binning from file %s \n",matrixfile.Data());      printf(" ERROR: cannot read rigidity binning from file ");
361        return(false);
362      };
363      //
364      brigm=(TArrayF*)file->Get("binrigmean");
365      if ( !brigm ){
366        printf(" ERROR: cannot read mean rigidity binning from file ");
367      return(false);      return(false);
368    };    };
369    //    //
370      for (Int_t i=0;i<17;i++){
371            TString name = Form("qplmeann%i",i);
372            qplmean[i] = (TArrayF*)file->Get(name.Data());
373            if ( !qplmean[i] ){
374                    printf(" ERROR: cannot read average from file ");
375                    return(false);
376            };
377      };
378      //
379    return(true);    return(true);
380    }
381    
382    Bool_t CaloFranzini::LoadMatrices(){
383      //
384      for (Int_t i=0;i<17;i++){
385            TString name1 = Form("matrixn%i",i);
386            hmat[i] = (TMatrixD*)file->Get(name1.Data());
387      };
388    //    //
389      return(true);
390  }  }
391    
392  TMatrixD *CaloFranzini::LoadCovarianceMatrix(Float_t rig){  TMatrixD *CaloFranzini::LoadCovarianceMatrix(Float_t rig){
393    //    //
394    TString name;    Int_t mv = 0;
395    for (Int_t i = 0; i<nbin-1; i++){    for (Int_t i = 0; i<nbin-1; i++){
396      if ( rig>=brig->At(i) && rig < brig->At(i+1) ){      if ( rig>=brig->At(i) && rig < brig->At(i+1) ){
397        name = Form("matrixn%i",i);        mv = i;
398        break;        break;
399      };      };
400    };    };
401    if ( rig < brig->At(0) ){    if ( rig < brig->At(0) ){
402      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));
403      name = "matrixn0";      mv = 0;
     printf(" Using matrix %s \n",name.Data());  
404    };    };
405    if ( rig >= brig->At(nbin) ){    if ( rig >= brig->At(nbin-1) ){
406      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));
407      name = Form("matrixn%i",nbin-1);      mv = nbin-2;
     printf(" Using matrix %s \n",name.Data());  
408    };    };
409    //    //
410    TMatrixD *matrix = (TMatrixD*)file->Get(name.Data());    return(hmat[mv]);
   //  
   return(matrix);  
411    //    //
412  }  }
413    
414    
415  TArrayF *CaloFranzini::LoadLongAverage(Float_t rig){  TArrayF *CaloFranzini::LoadLongAverage(Float_t rig){
416    //    //
417    TString name;    Int_t mv=0;
418    for (Int_t i = 0; i<nbin-1; i++){    for (Int_t i = 0; i<nbin-1; i++){
419      if ( rig>=brig->At(i) && rig < brig->At(i+1) ){      if ( rig>=brig->At(i) && rig < brig->At(i+1) ){
420        name = Form("qplmeann%i",i);          mv = i;
421        break;        break;
422      };      };
423    };    };
424    if ( rig < brig->At(0) ){    if ( rig < brig->At(0) ){
425      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));
426      name = "qplmeann0";      mv = 0;
427      printf(" Using qplmean %s \n",name.Data());    };
428      if ( rig >= brig->At(nbin-1) ){
429        printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brig->At(nbin-1));
430        mv=nbin-2;
431      };
432      //
433      return(qplmean[mv]);
434      //
435    }
436    
437    Float_t CaloFranzini::GetAverageAt(Int_t plane, Float_t rig){
438      //
439      Int_t therigb = 0;
440      for (Int_t i = 0; i<nbin-2; i++){
441        if ( rig>=brigm->At(i) && rig < brigm->At(i+1) ){
442          therigb = i;
443          break;
444        };
445      };
446      //
447      Float_t minrig;
448      Float_t maxrig;
449      //
450      //
451      if ( rig < brigm->At(0) ){
452        if ( rig < brig->At(0) ){
453          printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brigm->At(0));
454        };
455        therigb = 0;
456      };
457      if ( rig >= brigm->At(nbin-2) ){
458        if ( rig >= brig->At(nbin-2) ) {
459          printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brigm->At(nbin-2));
460        };
461        therigb = nbin - 3;
462      };
463      //
464      minrig = brigm->At(therigb);
465      maxrig = brigm->At(therigb+1);
466      //
467      Float_t minene = (*qplmean[therigb])[plane];
468      Float_t maxene = (*qplmean[therigb+1])[plane];
469      //
470      if ( maxrig == minrig ){
471       printf("Unrecoverable ERROR! Matrix will be screwed up... \n");
472       return(0.);
473    };    };
474    if ( rig >= brig->At(nbin) ){    Float_t b = log(maxene/minene)/(maxrig-minrig);
475      printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brig->At(nbin));    Float_t a = minene/exp(minrig*b);
476      name = Form("qplmeann%i",nbin-1);    Float_t ave = a*exp(b*rig);
477      printf(" Using qplmean %s \n",name.Data());    //
478      return(ave);
479      //
480    }
481    
482    Float_t CaloFranzini::GetHmatrixAt(Int_t iindex, Int_t jindex, Float_t rig){
483      Int_t therigb = 0;
484      for (Int_t i = 0; i<nbin-2; i++){
485        if ( rig>=brigm->At(i) && rig < brigm->At(i+1) ){
486          therigb = i;
487          break;
488        };
489    };    };
490    //    //
491    TArrayF *qplmean = (TArrayF*)file->Get(name.Data());    Float_t minrig;
492      Float_t maxrig;
493      //
494      if ( rig < brigm->At(0) ){
495        if ( rig < brig->At(0) ){
496          printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brigm->At(0));
497        };
498        therigb = 0;
499      };
500      if ( rig >= brigm->At(nbin-2) ){
501        if ( rig >= brig->At(nbin-2) ) {
502          printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brigm->At(nbin-2));
503        };
504        therigb = nbin - 3;
505      };
506      //
507      minrig = brigm->At(therigb);
508      maxrig = brigm->At(therigb+1);
509      //
510      Float_t minene = (*hmat[therigb])[iindex][jindex];
511      Float_t maxene = (*hmat[therigb+1])[iindex][jindex];
512    //    //
513    return(qplmean);    Float_t b = log(maxene/minene)/(maxrig-minrig);
514      Float_t a = minene/exp(minrig*b);
515      Float_t ave = a*exp(b*rig);
516      //
517      return(ave);
518    //    //
519  }  }
520    
521    void CaloFranzini::DrawLongAverage(Float_t rig){
522      //
523      TArrayF *ll = this->LoadLongAverage(rig);
524      //
525      gStyle->SetLabelSize(0.04);
526      gStyle->SetNdivisions(510,"XY");
527      //
528      TString hid = Form("cfralongvyvx");  
529      TCanvas *tcf  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
530      if ( tcf ){
531        tcf->cd();
532      } else {
533        tcf = new TCanvas(hid,hid);
534      };
535      //
536      TString thid = Form("hfralongvyvx");  
537      TH1F *thf  = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
538      if ( thf ) thf->Delete();
539      thf = new TH1F(thid,thid,44,-0.5,43.5);
540      tcf->cd();
541      Float_t qpl[43];
542      for (Int_t st=0;st<43;st++){
543        qpl[st] = ll->At(st);
544        printf("st %i qpl %f\n",st,qpl[st]);
545      };
546      for (Int_t st=0;st<44;st++){
547        if ( st == 37 ){
548          thf->Fill(st,0.);
549        } else {
550          Int_t ss = st;
551          if ( st > 37 ) ss--;
552          thf->Fill(st,qpl[ss]);
553        };
554      };
555      thf->Draw();
556      tcf->Modified();
557      tcf->Update();
558      //
559      gStyle->SetLabelSize(0);
560      gStyle->SetNdivisions(1,"XY");
561      //
562    };
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.4
  ViewVC Help
Powered by ViewVC 1.1.23