/[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.2 by mocchiut, Tue Dec 18 09:55:05 2007 UTC revision 1.3 by mocchiut, Thu Jan 3 10:02:26 2008 UTC
# Line 6  Line 6 
6  #include <TH1F.h>  #include <TH1F.h>
7  #include <TH2F.h>  #include <TH2F.h>
8  #include <TMatrixD.h>  #include <TMatrixD.h>
9    #include <TMatrixF.h>
10  #include <TArrayF.h>  #include <TArrayF.h>
11  #include <TArrayI.h>  #include <TArrayI.h>
12  //  //
# Line 15  Line 16 
16  using namespace std;  using namespace std;
17  //  //
18  extern Bool_t MATRIX;  extern Bool_t MATRIX;
19    extern Bool_t FULL;
20    extern Bool_t SIMU;
21    extern Bool_t CRIG;
22    extern Bool_t SRIG;
23  CaloFranzini *cf;  CaloFranzini *cf;
24  Int_t nbin;  Int_t nbin;
25  Float_t rig[18];  Float_t rig[18];
26  Float_t rmean[17];  Float_t rmean[17];
27    Int_t ntot[17];
28  //Float_t qqplane[17][43];  //Float_t qqplane[17][43];
29  //Int_t nnqplane[17][43];  //Int_t nnqplane[17][43];
30    
# Line 27  TArrayI *nqplane[17]; Line 33  TArrayI *nqplane[17];
33  TMatrixD *matrix[17];  TMatrixD *matrix[17];
34  TMatrixD *nmat[17];  TMatrixD *nmat[17];
35    
36    TMatrixD *fqplane[17];
37    TMatrixD *fnqplane[17];
38    TMatrixF *fmatrix[17];
39    TMatrixF *fnmat[17];
40    
41    
42  //===============================================================================  //===============================================================================
43  bool Select( PamLevel2* event ){  bool Select( PamLevel2* event ){
44    
# Line 41  bool Select( PamLevel2* event ){ Line 53  bool Select( PamLevel2* event ){
53    // tracker pre-selection    // tracker pre-selection
54    //------------------------------------------------------------------    //------------------------------------------------------------------
55    TrkTrack *trk = track->GetTrkTrack();    TrkTrack *trk = track->GetTrkTrack();
56      float rigidity = trk->GetRigidity();
57      if ( CRIG ) rigidity = event->GetCaloLevel2()->qtot/260.;
58      if ( SRIG ) rigidity = event->GetGPamela()->P0;
59    bool TRACK__OK = false;    bool TRACK__OK = false;
60    if(    if(
61       trk->chi2 >0     &&       trk->chi2 >0     &&
# Line 67  bool Select( PamLevel2* event ){ Line 81  bool Select( PamLevel2* event ){
81       !event->GetAcLevel2()->CARDhit() &&       !event->GetAcLevel2()->CARDhit() &&
82       !event->GetAcLevel2()->CAThit() &&       !event->GetAcLevel2()->CAThit() &&
83       true ) TOF__OK = true;       true ) TOF__OK = true;
84    if( !TOF__OK )return false;    if( !TOF__OK && !SIMU)return false;
85    //------------------------------------------------------    //------------------------------------------------------
86    // no albedo    // no albedo
87    //------------------------------------------------------    //------------------------------------------------------
88    if(    track->GetToFTrack()->beta[12]<=0.2 ||    if(    !SIMU && (track->GetToFTrack()->beta[12]<=0.2 ||
89           track->GetToFTrack()->beta[12] >= 1.5 ) return false;                     track->GetToFTrack()->beta[12] >= 1.5) ) return false;
90    
91    //------------------------------------------------------    //------------------------------------------------------
92    bool CUT1 = false;    bool CUT1 = false;
93    if(    if(
94       trk->nstep<100   &&       trk->nstep<100   &&
95       trk->GetRigidity()<400.   &&       rigidity<400.   &&
96       trk->GetRigidity()>0.1   &&       rigidity>0.1   &&
      trk->GetDeflection()<0.   &&  
97       trk->resx[0]<0.001 &&       trk->resx[0]<0.001 &&
98       trk->resx[5]<0.001 &&       trk->resx[5]<0.001 &&
99       track->IsSolved() &&       track->IsSolved() &&
# Line 88  bool Select( PamLevel2* event ){ Line 101  bool Select( PamLevel2* event ){
101       true ) CUT1 = true;       true ) CUT1 = true;
102    //------------------------------------------------------    //------------------------------------------------------
103    if( !CUT1 )return false;    if( !CUT1 )return false;
104      if ( trk->GetDeflection()>0. && !SIMU ) return false;
105    
106      //
107      //  ELENA'S CUT
108      //
109      //
110      // lever-arm 6
111      //====================================================
112      bool LX6=false;
113      if(
114         track->GetTrkTrack()->GetLeverArmX()==6 &&
115         !track->GetTrkTrack()->IsBad(0,0) &&
116         !track->GetTrkTrack()->IsBad(5,0) &&
117         track->GetTrkTrack()->resx[0]<0.001 &&
118         track->GetTrkTrack()->resx[5]<0.001 &&
119         track->GetTrkTrack()->IsInsideCavity() &&
120         true ) LX6 = true;
121      
122      //====================================================
123      // lever-arm 5
124      //====================================================
125      bool LX5=false;
126      if(
127         track->GetTrkTrack()->GetLeverArmX()==5 &&
128         true ){
129        if(
130           track->GetTrkTrack()->XGood(0) && track->GetTrkTrack()->XGood(4)
131           ){
132          
133          if(
134             !track->GetTrkTrack()->IsBad(0,0) &&
135             !track->GetTrkTrack()->IsBad(4,0) &&
136             track->GetTrkTrack()->resx[0]<0.001 &&
137             track->GetTrkTrack()->resx[4]<0.001 &&
138             track->GetTrkTrack()->IsInsideCavity() &&
139             true) LX5 = true;        
140        }else if (
141                  track->GetTrkTrack()->XGood(1) && track->GetTrkTrack()->XGood(5)
142                  ){
143          
144          if(
145             !track->GetTrkTrack()->IsBad(1,0) &&
146             !track->GetTrkTrack()->IsBad(5,0) &&
147             track->GetTrkTrack()->resx[1]<0.001 &&
148             track->GetTrkTrack()->resx[5]<0.001 &&
149             track->GetTrkTrack()->IsInsideCavity() &&
150             true) LX5 = true;
151        }
152      }
153      if ( !LX5 && !LX6 ) return false;
154    Float_t defl = trk->GetDeflection();    Float_t defl = trk->GetDeflection();
155    float p0 = 1.111588e+00;    float p0 = 1.111588e+00;
156    float p1 = 1.707656e+00;    float p1 = 1.707656e+00;
# Line 105  bool Select( PamLevel2* event ){ Line 168  bool Select( PamLevel2* event ){
168       track->GetTrkTrack()->chi2 < chi2m &&       track->GetTrkTrack()->chi2 < chi2m &&
169       true ) CUT2 = true;         true ) CUT2 = true;  
170    if ( !CUT2 ) return false;    if ( !CUT2 ) return false;
171      float dedxtrk    = trk->GetDEDX();
172      //  float zcutn =  9.   + 20./(rigidity*rigidity);
173      float zcut2 =  3.   + 4.3/(rigidity*rigidity);
174      float zcut1 =  0.52 + 0.455/(rigidity*rigidity);
175      Bool_t Z1 = false;
176      if(dedxtrk > zcut1 && dedxtrk < zcut2){
177        Z1=true;
178      }
179      if ( !Z1 && !SIMU ) return false;
180    //------------------------------------------------------    //------------------------------------------------------
181    //    //
182    // energy momentum match    // energy momentum match
# Line 116  bool Select( PamLevel2* event ){ Line 188  bool Select( PamLevel2* event ){
188    //    //
189    for (Int_t i=0; i < 22; i++){    for (Int_t i=0; i < 22; i++){
190      if ( track->GetCaloTrack()->tibar[i][1] < 0 || track->GetCaloTrack()->tibar[i][0] < 0 ){      if ( track->GetCaloTrack()->tibar[i][1] < 0 || track->GetCaloTrack()->tibar[i][0] < 0 ){
191          return false;        return false;
192        };
193      };
194      //
195      if ( event->GetCaloLevel2()->qtot == 0. ) return false;
196      if ( rigidity>5. && track->GetCaloTrack()->qtrack/event->GetCaloLevel2()->qtot < 0.4 ) return false;
197      if ( rigidity<1. && track->GetToFTrack()->beta[12] < 0.8 ) return false;
198      if ( rigidity>50. ){
199        if ( trk->GetNX()<5  &&
200             trk->GetNY()<4 ) return false;
201        //
202        Bool_t sphit = false;
203        for ( Int_t plane = 0; plane < 6; plane++){
204          if ( !trk->XGood(plane) ){
205            for (Int_t sing = 0; sing < event->GetTrkLevel2()->nclsx(); sing++){
206              TClonesArray &t = *(event->GetTrkLevel2()->SingletX);
207              TrkSinglet *singlet = (TrkSinglet*)t[sing];
208              if ( (singlet->plane-1) == plane ){
209                Float_t x = (singlet->coord[0]+singlet->coord[1])/2.;
210                if ( fabs(track->GetTrkTrack()->xv[plane] - x) < 1. ) sphit = true;
211              };
212            };
213          };
214          if ( !trk->YGood(plane) ){
215            for (Int_t sing = 0; sing < event->GetTrkLevel2()->nclsy(); sing++){
216              TClonesArray &t = *(event->GetTrkLevel2()->SingletY);
217              TrkSinglet *singlet = (TrkSinglet*)t[sing];
218              if ( (singlet->plane-1) == plane ){
219                Float_t x1 = (singlet->coord[0]);
220                Float_t x2 = (singlet->coord[1]);
221                if ( fabs(track->GetTrkTrack()->yv[plane] - x1) < 1. ) sphit = true;
222                if ( fabs(track->GetTrkTrack()->yv[plane] - x2) < 1. ) sphit = true;
223              };
224            };
225          };
226      };      };
227        if ( sphit ) return false; // spurious hit along the track
228    };    };
229    //    //
230    if ( event->GetCaloLevel1()->qtotpl(0) > 7. ) return false;    Int_t ti0 = track->GetCaloTrack()->tibar[0][1]-1;
231      
232      Int_t view = 0;
233      Int_t plane = 0;
234      Int_t strip = 0;
235      Float_t mip = 0.;
236      //
237      for ( Int_t i=0; i<event->GetCaloLevel1()->istrip; i++ ){
238        //
239        mip = event->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
240        if ( view == 1 && plane == 0 && strip == ti0 && mip > 4.) return false;
241        if ( view == 1 && (plane >0 || strip > ti0) ) break;
242      };
243      //  if ( event->GetCaloLevel1()->qtotpl(0) > 7. ) return false;
244    //    //
245    return true;    return true;
246  }  }
# Line 139  void CreateHistos( PamLevel2* event , TS Line 259  void CreateHistos( PamLevel2* event , TS
259    if ( MATRIX ){    if ( MATRIX ){
260      cf->UpdateMatrixFile(file.Data());      cf->UpdateMatrixFile(file.Data());
261      cf->LoadBin();      cf->LoadBin();
262        if ( !FULL ){
263          cf->LoadLong();
264        } else {
265          cf->LoadFull();
266        };
267    } else {    } else {
268      cf->CreateMatrixFile(file.Data());      cf->CreateMatrixFile(file.Data());
269    };    };
# Line 165  void CreateHistos( PamLevel2* event , TS Line 290  void CreateHistos( PamLevel2* event , TS
290    rig[17] = 400.;    rig[17] = 400.;
291    //    //
292    memset(rmean, 0, 17*sizeof(Float_t));    memset(rmean, 0, 17*sizeof(Float_t));
293      memset(ntot, 0, 17*sizeof(Int_t));
294    //    //
295    for (Int_t i=0; i < 17 ; i++){    for (Int_t i=0; i < 17 ; i++){
296      matrix[i] = new TMatrixD(43,43);      if ( !FULL ){
297      qplane[i] = new TArrayF(43);        matrix[i] = new TMatrixD(43,43);
298      nqplane[i] = new TArrayI(43);        qplane[i] = new TArrayF(43);
299      nmat[i] = new TMatrixD(43,43);        nqplane[i] = new TArrayI(43);
300          nmat[i] = new TMatrixD(43,43);
301        } else {
302          if ( MATRIX ){
303            fmatrix[i] = new TMatrixF(4128,4128);
304            fnmat[i] = new TMatrixF(4128,4128);
305            //fnmat[i] = new TMatrixI(8213,8213);
306          } else {
307            fqplane[i] = new TMatrixD(43,191); // 43 planes x 191 strip (= 1 + 95 x 2, one strip is the one transversed by the track that could be on the extreme right or left)
308            fnqplane[i] = new TMatrixD(43,191);//
309          };
310        };
311    };    };
312    //    //
313  }  }
# Line 179  void CreateHistos( PamLevel2* event , TS Line 316  void CreateHistos( PamLevel2* event , TS
316  void FindAverage( PamLevel2* L2, int iev ){  void FindAverage( PamLevel2* L2, int iev ){
317    //    //
318    Float_t erig = L2->GetTrack(0)->GetTrkTrack()->GetRigidity();    Float_t erig = L2->GetTrack(0)->GetTrkTrack()->GetRigidity();
319      if ( SRIG ) erig = L2->GetGPamela()->P0;
320      if ( CRIG ) erig = L2->GetCaloLevel2()->qtot/260.;
321    //    //
322    Int_t rbi = 0;    Int_t rbi = 0;
323    for (Int_t i = 0; i<nbin-1; i++){    for (Int_t i = 0; i<nbin-1; i++){
# Line 192  void FindAverage( PamLevel2* L2, int iev Line 331  void FindAverage( PamLevel2* L2, int iev
331    if ( erig >= rig[nbin-1] ) return;    if ( erig >= rig[nbin-1] ) return;
332    //    //
333    rmean[rbi] += erig;    rmean[rbi] += erig;
334      ntot[rbi]++;
335    //    //
336    Int_t dgf = 43;    if (!FULL ){
337    //      Int_t dgf = 43;
338    for (Int_t i=0; i<dgf; i++){      //
339      (*nqplane[rbi])[i]++;      for (Int_t i=0; i<dgf; i++){
340    };        (*nqplane[rbi])[i]++;
341    //      };
342    // Fill the estrip matrix      //
343    //      // Fill the estrip matrix
344    Int_t nplane = 0;      //
345    Int_t view = 0;      Int_t nplane = 0;
346    Int_t plane = 0;      Int_t view = 0;
347    Int_t strip = 0;      Int_t plane = 0;
348    Float_t mip = 0.;      Int_t strip = 0;
349    //      Float_t mip = 0.;
350    for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){      //
351        for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
352          //
353          mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
354          //
355          nplane = 1 - view + 2 * plane;
356          if ( erig > 4. && nplane == 0 && mip > 15. ) printf(" IEV %i erig %f OBT %u pkt %u file %s \n",iev,erig,L2->GetOrbitalInfo()->OBT,L2->GetOrbitalInfo()->pkt_num,L2->GetPamTree()->GetFile()->GetName());
357          //printf(" IEV %i OBT %u pkt %u file %s \n",iev,L2->GetOrbitalInfo()->OBT,L2->GetOrbitalInfo()->pkt_num,L2->GetPamTree()->GetFile()->GetName());
358          if ( nplane > 37 ) nplane--;
359          if ( nplane < dgf ){
360            (*qplane[rbi])[nplane] += mip;
361          };
362          //
363        };
364      } else {
365        //
366        // FULL CALORIMETER
367        //
368        CaloTrkVar *ct = L2->GetTrack(0)->GetCaloTrack();
369        //
370        Int_t nplane = 0;
371        Int_t view = 0;
372        Int_t plane = 0;
373        Int_t strip = 0;
374        Float_t mip = 0.;
375        //
376        Int_t cs = 0;
377        Int_t cd = 0;
378        Int_t mstrip = 0;
379        //
380        for (Int_t j=0; j<2; j++){
381          for (Int_t i=0; i<21; i++){
382            nplane = 1 - j + 2*i;
383            if ( nplane > 37 ) nplane--;
384            //
385            cs = ct->tibar[i][j] - 1;
386            //
387            cd = 95 - cs;  
388            //
389            for (Int_t k=0; k<191; k++){
390              mstrip = cd + k;
391              if ( mstrip < (191-cs) ) (*fnqplane[rbi])[nplane][mstrip] += 1.;
392            };
393          };
394        };
395      //      //
     mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);  
396      //      //
397      nplane = 1 - view + 2 * plane;      for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
398      if ( nplane > 37 ) nplane--;        //
399      if ( nplane < dgf ){        mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
400        (*qplane[rbi])[nplane] += mip;        //
401          nplane = 1 - view + 2 * plane;
402          if ( nplane > 37 ) nplane--;
403          //
404          cs = ct->tibar[plane][view] - 1;
405          //
406          cd = 95 - cs;
407          //
408          mstrip = cd + strip;
409          //
410          (*fqplane[rbi])[nplane][mstrip] += mip;
411          //
412      };      };
413        
414      //      //
415    };    };
416  }  }
417    
418  void CalculateAverage(){  void CalculateAverage(){
419    for (Int_t i=0; i<nbin-1; i++){    //  
420      if ( (*nqplane[i])[0] > 0 ) rmean[i] /= (Float_t)(*nqplane[i])[0];    if ( !FULL ){
421      for (Int_t j=0; j<43 ; j++){      for (Int_t i=0; i<nbin-1; i++){
422        if ( (*nqplane[i])[j] > 0 ){          if ( (*nqplane[i])[0] > 0 ) rmean[i] /= (Float_t)(*nqplane[i])[0];
423          (*qplane[i])[j] /= (Float_t)(*nqplane[i])[j];        for (Int_t j=0; j<43 ; j++){
424        } else {          if ( (*nqplane[i])[j] > 0 ){  
425          (*qplane[i])[j] = 0.;            (*qplane[i])[j] /= (Float_t)(*nqplane[i])[j];
426            } else {
427              (*qplane[i])[j] = 0.;
428            };
429            printf(" BIN %i plane %i average energy %f qplane %f nqplane %i \n",i,j,rmean[i],(*qplane[i])[j],(*nqplane[i])[j]);
430        };        };
431        printf(" BIN %i plane %i average energy %f qplane %f nqplane %i \n",i,j,rmean[i],(*qplane[i])[j],(*nqplane[i])[j]);      };
432        for (Int_t i=0; i<nbin-1; i++){
433          //
434          cf->WriteLongMean(qplane[i], i);
435          //
436        };
437      } else {
438    
439        for (Int_t i=0; i<nbin-1; i++){
440          if ( ntot[i] > 0 ) rmean[i] /= (Float_t)(ntot[i]);
441          //
442          for (Int_t j=0; j<43 ; j++){
443            for (Int_t k=0; k<191; k++){
444              if ( (*fnqplane[i])[j][k] > 0 ){  
445                (*fqplane[i])[j][k] /= (Float_t)(*fnqplane[i])[j][k];
446              } else {
447                (*fqplane[i])[j][k] = 0.;
448              };
449              printf(" BIN %i plane %i strip %i average energy %f qplane %f nqplane %f \n",i,j,k,rmean[i],(*fqplane[i])[j][k],(*fnqplane[i])[j][k]);
450            };
451          };
452        };
453        for (Int_t i=0; i<nbin-1; i++){
454          //
455          cf->WriteFullMean(fqplane[i], i);
456          //
457      };      };
458    };    };
459        //  
460    cf->WriteNumBin(nbin);    cf->WriteNumBin(nbin);
461      //
462    TArrayF *rigbin = new TArrayF(18, rig);    TArrayF *rigbin = new TArrayF(18, rig);
463    cf->WriteRigBin(rigbin);    cf->WriteRigBin(rigbin);
464      //
465    TArrayF *rmeanbin = new TArrayF(17, rmean);    TArrayF *rmeanbin = new TArrayF(17, rmean);
466    TFile *file =  cf->GetFile();    TFile *file =  cf->GetFile();
467    file->cd();    file->cd();
468    file->WriteObject(&(*rmeanbin), "binrigmean");    file->WriteObject(&(*rmeanbin), "binrigmean");
469    //    //
   for (Int_t i=0; i<nbin-1; i++){  
     //  
     cf->WriteLongMean(qplane[i], i);  
     //  
   };  
470    //    //
471  }  }
472    
# Line 255  void CalculateAverage(){ Line 474  void CalculateAverage(){
474  void FindMatrix( PamLevel2* L2, int iev ){  void FindMatrix( PamLevel2* L2, int iev ){
475    //    //
476    Float_t erig = L2->GetTrack(0)->GetTrkTrack()->GetRigidity();    Float_t erig = L2->GetTrack(0)->GetTrkTrack()->GetRigidity();
477      if ( SRIG ) erig = L2->GetGPamela()->P0;
478      if ( CRIG ) erig = L2->GetCaloLevel2()->qtot/260.;
479    //    //
480    Int_t rbi = 0;    Int_t rbi = 0;
481    for (Int_t i = 0; i<nbin-1; i++){    for (Int_t i = 0; i<nbin-1; i++){
# Line 267  void FindMatrix( PamLevel2* L2, int iev Line 488  void FindMatrix( PamLevel2* L2, int iev
488    if ( erig < rig[0] ) return;    if ( erig < rig[0] ) return;
489    if ( erig >= rig[nbin-1] ) return;    if ( erig >= rig[nbin-1] ) return;
490    //    //
491    Int_t dgf = 43;    if ( !FULL ){
492  //  printf("1qui \n");      Int_t dgf = 43;
493    //      //
494    for (Int_t i=0; i<dgf; i++){      for (Int_t i=0; i<dgf; i++){
495      for (Int_t j=0; j<dgf; j++){        for (Int_t j=0; j<dgf; j++){
496        (*nmat[rbi])[i][j] += 1.;          (*nmat[rbi])[i][j] += 1.;
497          };
498      };      };
   };  
 //  printf("2qui \n");  
   //  
   // Fill the estrip matrix  
   //  
   Int_t nplane = 0;  
   Int_t view = 0;  
   Int_t plane = 0;  
   Int_t strip = 0;  
   Float_t mip = 0.;  
   Float_t hpl[43];  
   memset(hpl,0,43*sizeof(Float_t));  
   for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){  
499      //      //
500      mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);      // Fill the estrip matrix
501      //      //
502      nplane = 1 - view + 2 * plane;      Int_t nplane = 0;
503      if ( nplane > 37 ) nplane--;      Int_t view = 0;
504      if ( nplane < dgf ){      Int_t plane = 0;
505        hpl[nplane] += mip;      Int_t strip = 0;
506        Float_t mip = 0.;
507        Float_t hpl[43];
508        memset(hpl,0,43*sizeof(Float_t));
509        for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
510          //
511          mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
512          //
513          nplane = 1 - view + 2 * plane;
514          if ( nplane > 37 ) nplane--;
515          if ( nplane < dgf ){
516            hpl[nplane] += mip;
517          };
518          //
519      };      };
520      //      //
521    };      for (Int_t i=0; i<dgf; i++){
522    //        for (Int_t j=0; j<dgf; j++){
523  //  printf("3qui \n");          (*matrix[rbi])[i][j] += (hpl[i] - cf->GetAverageAt(i,erig)) * (hpl[j] - cf->GetAverageAt(j,erig));
524    //        };
525    for (Int_t i=0; i<dgf; i++){      };
526      for (Int_t j=0; j<dgf; j++){    } else {
527        (*matrix[rbi])[i][j] += (hpl[i] - cf->GetAverageAt(i,erig)) * (hpl[j] - cf->GetAverageAt(j,erig));      //
528        // FULL CALORIMETER
529        //
530        CaloTrkVar *ct = L2->GetTrack(0)->GetCaloTrack();
531        //
532        Int_t nplane = 0;
533        Int_t view = 0;
534        Int_t plane = 0;
535        Int_t strip = 0;
536        Float_t mip = 0.;
537        //
538        Int_t mindgf = 48;
539        Int_t dgf = 143;
540        Int_t cs = 0;
541        Int_t cd = 0;
542        Int_t mstrip = 0;
543        //
544        Float_t mipv[43][191];
545        memset(mipv,0,43*191*sizeof(Float_t));
546        //
547        for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
548          //
549          mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
550          //
551          nplane = 1 - view + 2 * plane;
552          if ( nplane > 37 ) nplane--;
553          //
554          cs = ct->tibar[plane][view] - 1;
555          //
556          cd = 95 - cs;
557          //
558          mstrip = cd + strip;
559          //
560          mipv[nplane][mstrip] = mip;
561          //
562      };      };
563        //
564        Float_t mip1;
565        Int_t cs1;
566        Int_t cd1;
567        Float_t mip2;
568        Int_t cs2;
569        Int_t cd2;
570        Int_t mi = -1;
571        Int_t mj = -1;
572        Int_t nn1 = 0;
573        Int_t pl1 = 0;
574        Int_t vi1 = 0;
575        Int_t nn2 = 0;
576        Int_t pl2 = 0;
577        Int_t vi2 = 0;
578        Int_t mstrip1min = 0;
579        Int_t mstrip1max = 0;
580        Int_t mstrip2min = 0;
581        Int_t mstrip2max = 0;
582        //
583        for (Int_t nplane1=0; nplane1<43; nplane1++){
584          for (Int_t mstrip1=0; mstrip1<191; mstrip1++){
585            mj = -1;
586            //
587            if ( nplane1 >= 37 ) nn1 = nplane1 + 1;
588            vi1 = 1;
589            if ( nn1%2 ) vi1 = 0;
590            pl1 = (nn1 - 1 + vi1)/2;
591            //
592            cs1 = ct->tibar[pl1][vi1] - 1; // convertire nplane in pl1 e vi1
593            //
594            cd1 = 95 - cs1;
595            //
596            mstrip1min = cd1 + 0;
597            mstrip1max = cd1 + 95;
598            mip1 = mipv[nplane1][mstrip1];
599            //
600            if ( mstrip1 >= mindgf && mstrip1 <= dgf ){
601              mi++;
602              //
603              for (Int_t nplane2=0; nplane2<43; nplane2++){
604                for (Int_t mstrip2=0; mstrip2<191; mstrip2++){
605                  //
606                  if ( nplane2 >= 37 ) nn2 = nplane2 + 1;
607                  vi2 = 1;
608                  if ( nn2%2 ) vi2 = 0;
609                  pl1 = (nn2 - 1 + vi2)/2;
610                  //
611                  cs2 = ct->tibar[pl2][vi2] - 1;
612                  //
613                  cd2 = 95 - cs2;
614                  //
615                  mstrip2min = cd2 + 0;
616                  mstrip2max = cd2 + 95;
617                  //
618                  mip2 = mipv[nplane2][mstrip2];
619                  //
620                  if ( mstrip2 >= mindgf && mstrip2 <= dgf ){
621                    mj++;
622                    (*fmatrix[rbi])[mi][mj] += (mipv[nplane1][mstrip1] - cf->GetFullAverageAt(nplane1,mstrip1,erig)) * (mipv[nplane2][mstrip2] - cf->GetFullAverageAt(nplane2,mstrip2,erig));
623                    if ( mstrip1 >= mstrip1min && mstrip1 <= mstrip1max && mstrip2 >= mstrip2min && mstrip2 <= mstrip2max){
624                      (*fnmat[rbi])[mi][mj] += 1.;
625                    };
626                  };
627                };
628              };
629            };
630          };
631        };
632        //  
633    };    };
 //  printf("4qui \n");    
   //  gObjectTable->Print();  
634  }  }
635    
636  //===============================================================  //===============================================================
# Line 323  void SaveHistos(){ Line 647  void SaveHistos(){
647      //      //
648      printf("Finished, calculating average and inverting matrices\n");      printf("Finished, calculating average and inverting matrices\n");
649      //      //
650      for (Int_t i=0; i<nbin-1; i++){      if ( !FULL ){
651        //        for (Int_t i=0; i<nbin-1; i++){
652        // determine the average matrix          //
653        //              // determine the average matrix
654        for (Int_t ii=0; ii<43; ii++){          //    
655          for (Int_t j=0; j<43; j++){          for (Int_t ii=0; ii<43; ii++){
656            if ( (*nmat[i])[ii][j] > 0. ){            for (Int_t j=0; j<43; j++){
657              (*matrix[i])[ii][j] /= (*nmat[i])[ii][j];              if ( (*nmat[i])[ii][j] > 0. ){
658            } else {                (*matrix[i])[ii][j] /= (*nmat[i])[ii][j];
659              (*matrix[i])[ii][j] = 0.;              } else {
660                  (*matrix[i])[ii][j] = 0.;
661                };
662            };            };
663          };          };
664            //
665            cf->WriteLongMatrix(matrix[i],i);
666            //
667            if ( matrix[i]->Determinant() == 0. ){
668              printf("\n");
669              for (Int_t ii=0; ii<43; ii++){
670                for (Int_t j=0; j<43; j++){
671                  printf(" %.f",(*matrix[i])[ii][j]);        
672                };
673                printf("\n");
674              };
675              printf("\n");
676              printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i);
677            } else {
678              Double_t det = 0.;
679              TMatrixD invmatrix = (TMatrixD)(matrix[i]->Invert(&det));
680              printf(" Bin %i determinant is %f \n",i,det);
681              cf->WriteInvertedLongMatrix((TMatrixD)invmatrix,i);
682            };
683        };        };
684        } else {
685        //        //
686        cf->WriteLongMatrix(matrix[i],i);        // FULL
687        //        //
688        if ( matrix[i]->Determinant() == 0. ){        for (Int_t i=0; i<nbin-1; i++){
689          printf("\n");          //
690          for (Int_t ii=0; ii<43; ii++){          // determine the average matrix
691            for (Int_t j=0; j<43; j++){          //    
692              printf(" %.f",(*matrix[i])[ii][j]);            for (Int_t ii=0; ii<4171; ii++){
693              for (Int_t j=0; j<4171; j++){
694                if ( (*fnmat[i])[ii][j] > 0. ){
695                  (*fmatrix[i])[ii][j] /= (*fnmat[i])[ii][j];
696                } else {
697                  (*fmatrix[i])[ii][j] = 0.;
698                };
699            };            };
           printf("\n");  
700          };          };
701          printf("\n");          //
702          printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i);          cf->WriteFullMatrix(fmatrix[i],i);
703        } else {          //
704          Double_t det = 0.;          if ( fmatrix[i]->Determinant() == 0. ){
705          TMatrixD invmatrix = (TMatrixD)(matrix[i]->Invert(&det));            printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i);
706          printf(" Bin %i determinant is %f \n",i,det);          } else {
707          cf->WriteInvertedLongMatrix((TMatrixD)invmatrix,i);            Double_t det = 0.;
708              TMatrixF invmatrix = (TMatrixF)(fmatrix[i]->Invert(&det));
709              printf(" Bin %i determinant is %f \n",i,det);
710              cf->WriteInvertedFullMatrix((TMatrixF)invmatrix,i);
711            };
712        };        };
713      };      };
714      //      //
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.23