/[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.1 by mocchiut, Tue Dec 4 12:05:29 2007 UTC revision 1.3 by mocchiut, Thu Dec 13 17:08:11 2007 UTC
# Line 29  CaloFranzini::CaloFranzini(PamLevel2 *l2 Line 29  CaloFranzini::CaloFranzini(PamLevel2 *l2
29    debug = false;    debug = false;
30    dolong = true;    dolong = true;
31    dofull = false;    dofull = false;
32      sntr = 0;
33      //
34      sel = true;
35      cont = false;
36      N = 0;
37      NC = 43;
38      //
39      mask18b = -1;
40    //    //
41    Clear();    Clear();
42    //    //
# Line 39  void CaloFranzini::Clear(){ Line 47  void CaloFranzini::Clear(){
47    OBT = 0;    OBT = 0;
48    PKT = 0;    PKT = 0;
49    atime = 0;    atime = 0;
50    longtzeta = 0.;    longtzeta = -100.;
51    fulltzeta = 0.;    fulltzeta = -100.;
52      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 53  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  }  }
73    
74  void CaloFranzini::Delete(){  void CaloFranzini::Delete(){
75      //
76    if ( file ) file->Close();    if ( file ) file->Close();
77      //
78    Clear();    Clear();
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(){
108      this->Process(0);
109    }
110    
111    void CaloFranzini::Process(Int_t itr){
112      //  
113      if ( !L2 ){
114        printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n");
115        printf(" ERROR: CaloFranzini variables _NOT_ filled \n");
116        return;
117      };
118      //
119      if ( !nbin || !file || !brig ){
120        printf(" ERROR: it seems covariance matrix file has not been opened (CaloFranzini::Open()) \n");
121        printf(" ERROR: CaloFranzini variables _NOT_ filled \n");
122        return;
123      };
124      //
125      Bool_t newentry = false;
126      //
127      if ( L2->IsORB() ){
128        if ( L2->GetOrbitalInfo()->pkt_num != PKT || L2->GetOrbitalInfo()->OBT != OBT || L2->GetOrbitalInfo()->absTime != atime || itr != sntr ){
129          newentry = true;
130          OBT = L2->GetOrbitalInfo()->OBT;
131          PKT = L2->GetOrbitalInfo()->pkt_num;
132          atime = L2->GetOrbitalInfo()->absTime;
133          sntr = itr;
134        };
135      } else {
136        newentry = true;
137      };
138      //
139      if ( !newentry ) return;
140      //
141      // Some variables
142      //
143      if ( debug ) printf(" Processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime);
144      //
145      this->Clear();
146      //
147      Float_t rig = L2->GetTrack(itr)->GetTrkTrack()->GetRigidity();
148      //
149      // Fill the estrip matrix
150      //
151      Int_t nplane = 0;
152      Int_t view = 0;
153      Int_t plane = 0;
154      Int_t strip = 0;
155      Float_t mip = 0.;
156      for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
157        //
158        mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
159        //
160        estrip[view][plane][strip] = mip;
161        //    
162        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      //
170      //
171      //
172      if ( dolong ){
173        //
174        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++){
197          if ( L2->GetTrack(itr)->GetCaloTrack()->tibar[i][1] < 0 ){
198            dgf = 2 * i;
199            break;
200          };
201          if ( L2->GetTrack(itr)->GetCaloTrack()->tibar[i][0] < 0 ){
202          dgf = 1 + 2 * i;
203          break;
204          };
205        };
206        //
207        if ( dgf < 43 && dgf > 37 ) dgf--;
208        //
209        degfre = TMath::Min(dgf,NC);
210        //
211        if ( degfre > 0 ){
212          TArrayF *qplmean = this->LoadLongAverage(rig);
213          TMatrixD *lmatrix = this->LoadCovarianceMatrix(rig);
214          for (Int_t i = 0; i < degfre; i++){
215            if ( i != mask18b ){
216              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        };
224      };
225      if ( dofull ){
226        printf(" ERROR: NOT IMPLEMENTED YET\n");
227      printf(" ERROR: CaloFranzini variables _NOT_ filled \n");      printf(" ERROR: CaloFranzini variables _NOT_ filled \n");
228    };      };  
229    //    //
230    if ( debug ) this->Print();    if ( debug ) this->Print();
231    if ( debug ) printf(" exit \n");    if ( debug ) printf(" exit \n");
232    //  }
233    
234    
235    Float_t CaloFranzini::GetNormLongTZeta(){
236      Process();
237      Float_t normz = 0.;
238      if ( degfre > 0 ) normz = longtzeta/(Float_t)degfre;
239      return normz;
240    }
241    
242    Float_t CaloFranzini::GetNormFullTZeta(){
243      Process();
244      Float_t normz = 0.;
245      if ( degfre > 0 ) normz = fulltzeta/(Float_t)degfre;
246      return normz;
247  }  }
248    
249  Bool_t CaloFranzini::CreateMatrixFile(TString matrixfile){    Bool_t CaloFranzini::CreateMatrixFile(TString matrixfile){  
# Line 76  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 86  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 160  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 187  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 198  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.1  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.23