/[PAMELA software]/calo/flight/CaloFranzini/src/Calib.cpp
ViewVC logotype

Diff of /calo/flight/CaloFranzini/src/Calib.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.1 by mocchiut, Thu Dec 13 17:08:09 2007 UTC revision 1.2 by mocchiut, Tue Dec 18 09:55:05 2007 UTC
# Line 14  Line 14 
14  //  //
15  using namespace std;  using namespace std;
16  //  //
17    extern Bool_t MATRIX;
18  CaloFranzini *cf;  CaloFranzini *cf;
19  Int_t nbin;  Int_t nbin;
20  Float_t rig[18];  Float_t rig[18];
21  Float_t rmean[18];  Float_t rmean[17];
22  Float_t qqplane[18][43];  //Float_t qqplane[17][43];
23  Int_t nnqplane[18][43];  //Int_t nnqplane[17][43];
24    
25  TArrayF *qplane[18];  TArrayF *qplane[17];
26  TArrayI *nqplane[18];  TArrayI *nqplane[17];
27  TMatrixD *matrix[18];  TMatrixD *matrix[17];
28  TMatrixD *nmat[18];  TMatrixD *nmat[17];
29    
30  //===============================================================================  //===============================================================================
31  bool Select( PamLevel2* event ){  bool Select( PamLevel2* event ){
# Line 134  bool Select( PamLevel2* event ){ Line 135  bool Select( PamLevel2* event ){
135  void CreateHistos( PamLevel2* event , TString file){  void CreateHistos( PamLevel2* event , TString file){
136    
137    cf = new CaloFranzini(event);    cf = new CaloFranzini(event);
138    cf->CreateMatrixFile(file.Data());    //
139      if ( MATRIX ){
140        cf->UpdateMatrixFile(file.Data());
141        cf->LoadBin();
142      } else {
143        cf->CreateMatrixFile(file.Data());
144      };
145      //
146    //    //
147    nbin = 18;    nbin = 18;
148    rig[0] = 0.1;    rig[0] = 0.1;
# Line 156  void CreateHistos( PamLevel2* event , TS Line 164  void CreateHistos( PamLevel2* event , TS
164    rig[16] = 200.;    rig[16] = 200.;
165    rig[17] = 400.;    rig[17] = 400.;
166    //    //
167  //  memset(qplane, 0, 44*18*sizeof(Float_t));    memset(rmean, 0, 17*sizeof(Float_t));
   memset(rmean, 0, 18*sizeof(Float_t));  
   memset(qqplane, 0, 43*18*sizeof(Float_t));  
   memset(nnqplane, 0, 43*18*sizeof(Float_t));  
168    //    //
169    for (Int_t i=0; i < 18 ; i++){    for (Int_t i=0; i < 17 ; i++){
170      matrix[i] = new TMatrixD(43,43);      matrix[i] = new TMatrixD(43,43);
171      qplane[i] = new TArrayF(43);      qplane[i] = new TArrayF(43);
172      nqplane[i] = new TArrayI(43);      nqplane[i] = new TArrayI(43);
173      nmat[i] = new TMatrixD(43,43);      nmat[i] = new TMatrixD(43,43);
174    };    };
 //  for (Int_t i=0; i < 18 ; i++){  
 //      for (Int_t j = 0; j < 43;i++){  
 //              (*qplane[i])[j] = 0.;  
 //              qplane[i]->Reset();  
 //              nqplane[i]->Reset();  
 //              for (Int_t k = 0; k < 43;k++){  
 //                      (*matrix[i])[j][k] = 0.;  
 //                      (*nmat[i])[j][k] = 0.;  
 //              };  
 //      };  
 //  };  
175    //    //
176  }  }
177    
178  //===============================================================  //===============================================================
179  void FindAverage( PamLevel2* L2, int iev ){  void FindAverage( PamLevel2* L2, int iev ){
180    //    //
   //  printf("SELEZIONATO \n");  
181    Float_t erig = L2->GetTrack(0)->GetTrkTrack()->GetRigidity();    Float_t erig = L2->GetTrack(0)->GetTrkTrack()->GetRigidity();
182    //    //
183    Int_t rbi = 0;    Int_t rbi = 0;
184    for (Int_t i = 0; i<nbin-1; i++){ // < 17    for (Int_t i = 0; i<nbin-1; i++){
185      if ( erig>=rig[i] && erig < rig[i+1] ){      if ( erig>=rig[i] && erig < rig[i+1] ){
186        rbi = i;        rbi = i;
187        break;        break;
188      };      };
189    };    };
190    //    //
   //  printf("SELEZIONATO rbi %i erig %f rig[0] %f rig[nbin] %f \n",rbi, erig, rig[0],rig[nbin-1]);  
   //  
191    if ( erig < rig[0] ) return;    if ( erig < rig[0] ) return;
192    if ( erig >= rig[nbin-1] ) return;    if ( erig >= rig[nbin-1] ) return;
193    //    //
# Line 204  void FindAverage( PamLevel2* L2, int iev Line 195  void FindAverage( PamLevel2* L2, int iev
195    //    //
196    Int_t dgf = 43;    Int_t dgf = 43;
197    //    //
 //  for (Int_t i=0; i < 22; i++){  
 //    if ( L2->GetTrack(0)->GetCaloTrack()->tibar[i][1] < 0 ){  
 //      dgf = 2 * i;  
  //     break;  
   //  };  
  //   if ( L2->GetTrack(0)->GetCaloTrack()->tibar[i][0] < 0 ){  
  //     dgf = 1 + 2 * i;  
  //     break;  
  //   };  
  // };  
   //  
 //  if ( dgf < 43 && dgf > 37 ) dgf--;  
   //  
198    for (Int_t i=0; i<dgf; i++){    for (Int_t i=0; i<dgf; i++){
199      (*nqplane[rbi])[i]++;      (*nqplane[rbi])[i]++;
         nnqplane[rbi][i]++;  
200    };    };
201    //    //
   //  printf(" SELEZIONATO rig %f dgf %i \n",erig,dgf);  
   //  
   //  
202    // Fill the estrip matrix    // Fill the estrip matrix
203    //    //
204    Int_t nplane = 0;    Int_t nplane = 0;
# Line 239  void FindAverage( PamLevel2* L2, int iev Line 213  void FindAverage( PamLevel2* L2, int iev
213      //      //
214      nplane = 1 - view + 2 * plane;      nplane = 1 - view + 2 * plane;
215      if ( nplane > 37 ) nplane--;      if ( nplane > 37 ) nplane--;
     //    printf(" mip %f view %i plane %i => nplane %i \n",mip,view,plane,nplane);  
216      if ( nplane < dgf ){      if ( nplane < dgf ){
       //      printf(" PRIMA: qpl %f nqpl %i \n",(*qplane[rbi])[nplane],(*nqplane[rbi])[nplane]);  
217        (*qplane[rbi])[nplane] += mip;        (*qplane[rbi])[nplane] += mip;
         qqplane[rbi][nplane] += mip;  
       //      printf(" DOPO: qpl %f nqpl %i \n",(*qplane[rbi])[nplane],(*nqplane[rbi])[nplane]);  
218      };      };
219      //      //
220    };    };
221  }  }
222    
223    void CalculateAverage(){
224      for (Int_t i=0; i<nbin-1; i++){
225        if ( (*nqplane[i])[0] > 0 ) rmean[i] /= (Float_t)(*nqplane[i])[0];
226        for (Int_t j=0; j<43 ; j++){
227          if ( (*nqplane[i])[j] > 0 ){  
228            (*qplane[i])[j] /= (Float_t)(*nqplane[i])[j];
229          } else {
230            (*qplane[i])[j] = 0.;
231          };
232          printf(" BIN %i plane %i average energy %f qplane %f nqplane %i \n",i,j,rmean[i],(*qplane[i])[j],(*nqplane[i])[j]);
233        };
234      };
235      
236      cf->WriteNumBin(nbin);
237    
238      TArrayF *rigbin = new TArrayF(18, rig);
239      cf->WriteRigBin(rigbin);
240    
241      TArrayF *rmeanbin = new TArrayF(17, rmean);
242      TFile *file =  cf->GetFile();
243      file->cd();
244      file->WriteObject(&(*rmeanbin), "binrigmean");
245      //
246      for (Int_t i=0; i<nbin-1; i++){
247        //
248        cf->WriteLongMean(qplane[i], i);
249        //
250      };
251      //
252    }
253    
254  //===============================================================  //===============================================================
255  void FindMatrix( PamLevel2* L2, int iev ){  void FindMatrix( PamLevel2* L2, int iev ){
256    //    //
   //  printf("SELEZIONATO \n");  
257    Float_t erig = L2->GetTrack(0)->GetTrkTrack()->GetRigidity();    Float_t erig = L2->GetTrack(0)->GetTrkTrack()->GetRigidity();
258    //    //
259    Int_t rbi = 0;    Int_t rbi = 0;
# Line 264  void FindMatrix( PamLevel2* L2, int iev Line 264  void FindMatrix( PamLevel2* L2, int iev
264      };      };
265    };    };
266    //    //
   //  printf("SELEZIONATO rbi %i erig %f rig[0] %f rig[nbin] %f \n",rbi, erig, rig[0],rig[nbin-1]);  
   //  
267    if ( erig < rig[0] ) return;    if ( erig < rig[0] ) return;
268    if ( erig >= rig[nbin-1] ) return;    if ( erig >= rig[nbin-1] ) return;
269    //    //
270    Int_t dgf = 43;    Int_t dgf = 43;
271    //  //  printf("1qui \n");
 //  for (Int_t i=0; i < 22; i++){  
 //    if ( L2->GetTrack(0)->GetCaloTrack()->tibar[i][1] < 0 ){  
 //      dgf = 2 * i;  
  //     break;  
  //   };  
  //   if ( L2->GetTrack(0)->GetCaloTrack()->tibar[i][0] < 0 ){  
  //     dgf = 1 + 2 * i;  
  //     break;  
  //   };  
 //  };  
   //  
 //  if ( dgf < 43 && dgf > 37 ) dgf--;  
272    //    //
273    for (Int_t i=0; i<dgf; i++){    for (Int_t i=0; i<dgf; i++){
274      for (Int_t j=0; j<dgf; j++){      for (Int_t j=0; j<dgf; j++){
275        (*nmat[rbi])[i][j] += 1.;        (*nmat[rbi])[i][j] += 1.;
276      };      };
277    };    };
278    //  //  printf("2qui \n");
   //  printf(" SELEZIONATO rig %f dgf %i \n",erig,dgf);  
   //  
279    //    //
280    // Fill the estrip matrix    // Fill the estrip matrix
281    //    //
# Line 314  void FindMatrix( PamLevel2* L2, int iev Line 298  void FindMatrix( PamLevel2* L2, int iev
298      //      //
299    };    };
300    //    //
301    //  printf("3qui \n");
302      //
303    for (Int_t i=0; i<dgf; i++){    for (Int_t i=0; i<dgf; i++){
304      for (Int_t j=0; j<dgf; j++){      for (Int_t j=0; j<dgf; j++){
305        // if ( rbi == 1 ) printf(" PRIMA: matrix %f qpli %f qplj %f hpli %f hplj %f \n",(*matrix[rbi])[i][j],(*qplane[rbi])[i],(*qplane[rbi])[j],hpl[i],hpl[j]);        (*matrix[rbi])[i][j] += (hpl[i] - cf->GetAverageAt(i,erig)) * (hpl[j] - cf->GetAverageAt(j,erig));
       (*matrix[rbi])[i][j] += (hpl[i] - (*qplane[rbi])[i]) * (hpl[j] - (*qplane[rbi])[j]);  
       //      if ( rbi == 1 ) printf(" DOPO: matrix %f qpli %f qplj %f hpli %f hplj %f \n",(*matrix[rbi])[i][j],(*qplane[rbi])[i],(*qplane[rbi])[j],hpl[i],hpl[j]);  
306      };      };
307    };    };
308    //  printf("4qui \n");  
309      //  gObjectTable->Print();
310  }  }
311    
 void CalculateAverage(){  
   for (Int_t i=0; i<nbin; i++){  
     if ( (*nqplane[i])[0] > 0 ) rmean[i] /= (Float_t)(*nqplane[i])[0];  
     for (Int_t j=0; j<43 ; j++){  
       if ( (*nqplane[i])[j] > 0 ){    
         (*qplane[i])[j] /= (Float_t)(*nqplane[i])[j];  
       } else {  
         (*qplane[i])[j] = 0.;  
       };  
       printf(" BIN %i plane %i average energy %f qplane %f nqplane %i \n",i,j,rmean[i],(*qplane[i])[j],(*nqplane[i])[j]);  
 //      if ( nnqplane[i][j] > 0 ){    
 //      qqplane[i][j] /= (Float_t)nnqplane[i][j];  
 //      } else {  
 //      qqplane[i][j] = 0.;  
 //      };  
 //      printf(" BBIN %i plane %i average energy %f qplane %f nqplane %i \n",i,j,rmean[i],qqplane[i][j],nnqplane[i][j]);  
     };  
   };  
     
   cf->WriteNumBin(nbin);  
   
   TArrayF *rigbin = new TArrayF(18, rig);  
   cf->WriteRigBin(rigbin);  
   
   TArrayF *rmeanbin = new TArrayF(18, rmean);  
   TFile *file =  cf->GetFile();  
   file->cd();  
   //  rigbin->Write("binrig");  
   file->WriteObject(&(*rmeanbin), "binrigmean");  
   
 //  cf->WriteRigBin(rigbin);  
     
   for (Int_t i=0; i<nbin; i++){  
       
     cf->WriteLongMean(qplane[i], i);  
   
   };  
 }  
312  //===============================================================  //===============================================================
313  // Save histograms  // Save histograms
314  //  //
# Line 370  void CalculateAverage(){ Line 318  void CalculateAverage(){
318  //  //
319  //===============================================================  //===============================================================
320  void SaveHistos(){  void SaveHistos(){
321        //
322    printf("Finished, calculating average and inverting matrices\n");    if ( MATRIX ){
323        //
324        printf("Finished, calculating average and inverting matrices\n");
     
 //  cf->WriteNumBin(nbin);  
   
 //  TArrayF *rigbin = new TArrayF(18, rig);  
 //  cf->WriteRigBin(rigbin);  
     
   for (Int_t i=0; i<nbin; i++){  
       
 //    cf->WriteLongMean(qplane[i], i);  
   
325      //      //
326      // determine the average matrix      for (Int_t i=0; i<nbin-1; i++){
327      //            //
328      for (Int_t ii=0; ii<43; ii++){        // determine the average matrix
329        for (Int_t j=0; j<43; j++){        //    
330          //      if ( i == 1 ) printf(" num %f \n",(*nmat[i])[ii][j]);        for (Int_t ii=0; ii<43; ii++){
331          if ( (*nmat[i])[ii][j] > 0. ){          for (Int_t j=0; j<43; j++){
332            //      if ( i == 1 ) printf(" value %f \n",(*matrix[i])[ii][j]);            if ( (*nmat[i])[ii][j] > 0. ){
333            (*matrix[i])[ii][j] /= (*nmat[i])[ii][j];              (*matrix[i])[ii][j] /= (*nmat[i])[ii][j];
334          } else {            } else {
335            (*matrix[i])[ii][j] = 0.;              (*matrix[i])[ii][j] = 0.;
336              };
337          };          };
338        };        };
339      };        //
340              cf->WriteLongMatrix(matrix[i],i);
341  //     if ( i == 2 ){        //
342  //       printf("\n");        if ( matrix[i]->Determinant() == 0. ){
343  //       for (Int_t ii=0; ii<43; ii++){          printf("\n");
344  //      for (Int_t j=0; j<43; j++){          for (Int_t ii=0; ii<43; ii++){
345  //        printf(" %.f",(*matrix[i])[ii][j]);                for (Int_t j=0; j<43; j++){
346  //      };              printf(" %.f",(*matrix[i])[ii][j]);  
347  //      printf("\n");            };
348  //       };            printf("\n");
349  //       printf("\n");          };
350  //     };          printf("\n");
351                printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i);
352          } else {
353      cf->WriteLongMatrix(matrix[i],i);          Double_t det = 0.;
354                TMatrixD invmatrix = (TMatrixD)(matrix[i]->Invert(&det));
355      if ( matrix[i]->Determinant() == 0. ){          printf(" Bin %i determinant is %f \n",i,det);
356        printf("\n");          cf->WriteInvertedLongMatrix((TMatrixD)invmatrix,i);
       for (Int_t ii=0; ii<43; ii++){  
         for (Int_t j=0; j<43; j++){  
           printf(" %.f",(*matrix[i])[ii][j]);      
         };  
         printf("\n");  
357        };        };
       printf("\n");  
       printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i);  
     } else {  
       Double_t det = 0.;  
       TMatrixD invmatrix = (TMatrixD)(matrix[i]->Invert(&det));  
       printf(" Bin %i determinant is %f \n",i,det);  
       cf->WriteInvertedLongMatrix((TMatrixD)invmatrix,i);  
358      };      };
359        //
360        printf(" done, closing file and exiting\n");
361        //
362    };    };
363        //
   printf(" done, closing file and exiting\n");  
   
364    cf->CloseMatrixFile();    cf->CloseMatrixFile();
365      //
366    cf->Delete();    cf->Delete();
367      //
368  }  }
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.23