/[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.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 14  Line 15 
15  //  //
16  using namespace std;  using namespace std;
17  //  //
18    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[18];  Float_t rmean[17];
27  Float_t qqplane[18][43];  Int_t ntot[17];
28  Int_t nnqplane[18][43];  //Float_t qqplane[17][43];
29    //Int_t nnqplane[17][43];
30  TArrayF *qplane[18];  
31  TArrayI *nqplane[18];  TArrayF *qplane[17];
32  TMatrixD *matrix[18];  TArrayI *nqplane[17];
33  TMatrixD *nmat[18];  TMatrixD *matrix[17];
34    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 ){
# Line 40  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 66  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 87  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 104  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 115  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->GetCaloLevel1()->qtotpl(0) > 7. ) return false;    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      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 134  bool Select( PamLevel2* event ){ Line 255  bool Select( PamLevel2* event ){
255  void CreateHistos( PamLevel2* event , TString file){  void CreateHistos( PamLevel2* event , TString file){
256    
257    cf = new CaloFranzini(event);    cf = new CaloFranzini(event);
258    cf->CreateMatrixFile(file.Data());    //
259      if ( MATRIX ){
260        cf->UpdateMatrixFile(file.Data());
261        cf->LoadBin();
262        if ( !FULL ){
263          cf->LoadLong();
264        } else {
265          cf->LoadFull();
266        };
267      } else {
268        cf->CreateMatrixFile(file.Data());
269      };
270      //
271    //    //
272    nbin = 18;    nbin = 18;
273    rig[0] = 0.1;    rig[0] = 0.1;
# Line 156  void CreateHistos( PamLevel2* event , TS Line 289  void CreateHistos( PamLevel2* event , TS
289    rig[16] = 200.;    rig[16] = 200.;
290    rig[17] = 400.;    rig[17] = 400.;
291    //    //
292  //  memset(qplane, 0, 44*18*sizeof(Float_t));    memset(rmean, 0, 17*sizeof(Float_t));
293    memset(rmean, 0, 18*sizeof(Float_t));    memset(ntot, 0, 17*sizeof(Int_t));
294    memset(qqplane, 0, 43*18*sizeof(Float_t));    //
295    memset(nnqplane, 0, 43*18*sizeof(Float_t));    for (Int_t i=0; i < 17 ; i++){
296    //      if ( !FULL ){
297    for (Int_t i=0; i < 18 ; i++){        matrix[i] = new TMatrixD(43,43);
298      matrix[i] = new TMatrixD(43,43);        qplane[i] = new TArrayF(43);
299      qplane[i] = new TArrayF(43);        nqplane[i] = new TArrayI(43);
300      nqplane[i] = new TArrayI(43);        nmat[i] = new TMatrixD(43,43);
301      nmat[i] = new TMatrixD(43,43);      } else {
302    };        if ( MATRIX ){
303  //  for (Int_t i=0; i < 18 ; i++){          fmatrix[i] = new TMatrixF(4128,4128);
304  //      for (Int_t j = 0; j < 43;i++){          fnmat[i] = new TMatrixF(4128,4128);
305  //              (*qplane[i])[j] = 0.;          //fnmat[i] = new TMatrixI(8213,8213);
306  //              qplane[i]->Reset();        } else {
307  //              nqplane[i]->Reset();          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  //              for (Int_t k = 0; k < 43;k++){          fnqplane[i] = new TMatrixD(43,191);//
309  //                      (*matrix[i])[j][k] = 0.;        };
310  //                      (*nmat[i])[j][k] = 0.;      };
311  //              };    };
 //      };  
 //  };  
312    //    //
313  }  }
314    
315  //===============================================================  //===============================================================
316  void FindAverage( PamLevel2* L2, int iev ){  void FindAverage( PamLevel2* L2, int iev ){
317    //    //
   //  printf("SELEZIONATO \n");  
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++){ // < 17    for (Int_t i = 0; i<nbin-1; i++){
324      if ( erig>=rig[i] && erig < rig[i+1] ){      if ( erig>=rig[i] && erig < rig[i+1] ){
325        rbi = i;        rbi = i;
326        break;        break;
327      };      };
328    };    };
329    //    //
   //  printf("SELEZIONATO rbi %i erig %f rig[0] %f rig[nbin] %f \n",rbi, erig, rig[0],rig[nbin-1]);  
   //  
330    if ( erig < rig[0] ) return;    if ( erig < rig[0] ) return;
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 < 22; i++){      //
339  //    if ( L2->GetTrack(0)->GetCaloTrack()->tibar[i][1] < 0 ){      for (Int_t i=0; i<dgf; i++){
340  //      dgf = 2 * i;        (*nqplane[rbi])[i]++;
341   //     break;      };
342    //  };      //
343   //   if ( L2->GetTrack(0)->GetCaloTrack()->tibar[i][0] < 0 ){      // Fill the estrip matrix
344   //     dgf = 1 + 2 * i;      //
345   //     break;      Int_t nplane = 0;
346   //   };      Int_t view = 0;
347   // };      Int_t plane = 0;
348    //      Int_t strip = 0;
349  //  if ( dgf < 43 && dgf > 37 ) dgf--;      Float_t mip = 0.;
350    //      //
351    for (Int_t i=0; i<dgf; i++){      for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
352      (*nqplane[rbi])[i]++;        //
353          nnqplane[rbi][i]++;        mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
354    };        //
355    //        nplane = 1 - view + 2 * plane;
356    //  printf(" SELEZIONATO rig %f dgf %i \n",erig,dgf);        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    // Fill the estrip matrix        if ( nplane < dgf ){
360    //          (*qplane[rbi])[nplane] += mip;
361    Int_t nplane = 0;        };
362    Int_t view = 0;        //
363    Int_t plane = 0;      };
364    Int_t strip = 0;    } else {
365    Float_t mip = 0.;      //
366    //      // FULL CALORIMETER
367    for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){      //
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      //    printf(" mip %f view %i plane %i => nplane %i \n",mip,view,plane,nplane);        mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
400      if ( nplane < dgf ){        //
401        //      printf(" PRIMA: qpl %f nqpl %i \n",(*qplane[rbi])[nplane],(*nqplane[rbi])[nplane]);        nplane = 1 - view + 2 * plane;
402        (*qplane[rbi])[nplane] += mip;        if ( nplane > 37 ) nplane--;
403          qqplane[rbi][nplane] += mip;        //
404        //      printf(" DOPO: qpl %f nqpl %i \n",(*qplane[rbi])[nplane],(*nqplane[rbi])[nplane]);        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(){
419      //  
420      if ( !FULL ){
421        for (Int_t i=0; i<nbin-1; i++){
422          if ( (*nqplane[i])[0] > 0 ) rmean[i] /= (Float_t)(*nqplane[i])[0];
423          for (Int_t j=0; j<43 ; j++){
424            if ( (*nqplane[i])[j] > 0 ){  
425              (*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        };
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);
461      //
462      TArrayF *rigbin = new TArrayF(18, rig);
463      cf->WriteRigBin(rigbin);
464      //
465      TArrayF *rmeanbin = new TArrayF(17, rmean);
466      TFile *file =  cf->GetFile();
467      file->cd();
468      file->WriteObject(&(*rmeanbin), "binrigmean");
469      //
470      //
471    }
472    
473  //===============================================================  //===============================================================
474  void FindMatrix( PamLevel2* L2, int iev ){  void FindMatrix( PamLevel2* L2, int iev ){
475    //    //
   //  printf("SELEZIONATO \n");  
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 264  void FindMatrix( PamLevel2* L2, int iev Line 485  void FindMatrix( PamLevel2* L2, int iev
485      };      };
486    };    };
487    //    //
   //  printf("SELEZIONATO rbi %i erig %f rig[0] %f rig[nbin] %f \n",rbi, erig, rig[0],rig[nbin-1]);  
   //  
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    //      Int_t dgf = 43;
493  //  for (Int_t i=0; i < 22; i++){      //
494  //    if ( L2->GetTrack(0)->GetCaloTrack()->tibar[i][1] < 0 ){      for (Int_t i=0; i<dgf; i++){
495  //      dgf = 2 * i;        for (Int_t j=0; j<dgf; j++){
496   //     break;          (*nmat[rbi])[i][j] += 1.;
497   //   };        };
  //   if ( L2->GetTrack(0)->GetCaloTrack()->tibar[i][0] < 0 ){  
  //     dgf = 1 + 2 * i;  
  //     break;  
  //   };  
 //  };  
   //  
 //  if ( dgf < 43 && dgf > 37 ) dgf--;  
   //  
   for (Int_t i=0; i<dgf; i++){  
     for (Int_t j=0; j<dgf; j++){  
       (*nmat[rbi])[i][j] += 1.;  
498      };      };
   };  
   //  
   //  printf(" SELEZIONATO rig %f dgf %i \n",erig,dgf);  
   //  
   //  
   // 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    for (Int_t i=0; i<dgf; i++){          (*matrix[rbi])[i][j] += (hpl[i] - cf->GetAverageAt(i,erig)) * (hpl[j] - cf->GetAverageAt(j,erig));
524      for (Int_t j=0; j<dgf; j++){        };
       // 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] - (*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]);  
525      };      };
526    };    } else {
527  }      //
528        // FULL CALORIMETER
529  void CalculateAverage(){      //
530    for (Int_t i=0; i<nbin; i++){      CaloTrkVar *ct = L2->GetTrack(0)->GetCaloTrack();
531      if ( (*nqplane[i])[0] > 0 ) rmean[i] /= (Float_t)(*nqplane[i])[0];      //
532      for (Int_t j=0; j<43 ; j++){      Int_t nplane = 0;
533        if ( (*nqplane[i])[j] > 0 ){        Int_t view = 0;
534          (*qplane[i])[j] /= (Float_t)(*nqplane[i])[j];      Int_t plane = 0;
535        } else {      Int_t strip = 0;
536          (*qplane[i])[j] = 0.;      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        };        };
       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]);  
631      };      };
632    };      //  
     
   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);  
   
633    };    };
634  }  }
635    
636  //===============================================================  //===============================================================
637  // Save histograms  // Save histograms
638  //  //
# Line 370  void CalculateAverage(){ Line 642  void CalculateAverage(){
642  //  //
643  //===============================================================  //===============================================================
644  void SaveHistos(){  void SaveHistos(){
645        //
646    printf("Finished, calculating average and inverting matrices\n");    if ( MATRIX ){
647        //
648        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);  
   
649      //      //
650      // determine the average matrix      if ( !FULL ){
651      //            for (Int_t i=0; i<nbin-1; i++){
652      for (Int_t ii=0; ii<43; ii++){          //
653        for (Int_t j=0; j<43; j++){          // determine the average matrix
654          //      if ( i == 1 ) printf(" num %f \n",(*nmat[i])[ii][j]);          //    
655          if ( (*nmat[i])[ii][j] > 0. ){          for (Int_t ii=0; ii<43; ii++){
656            //      if ( i == 1 ) printf(" value %f \n",(*matrix[i])[ii][j]);            for (Int_t j=0; j<43; j++){
657            (*matrix[i])[ii][j] /= (*nmat[i])[ii][j];              if ( (*nmat[i])[ii][j] > 0. ){
658                  (*matrix[i])[ii][j] /= (*nmat[i])[ii][j];
659                } 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 {          } else {
678            (*matrix[i])[ii][j] = 0.;            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        };        };
     };  
       
 //     if ( i == 2 ){  
 //       printf("\n");  
 //       for (Int_t ii=0; ii<43; ii++){  
 //      for (Int_t j=0; j<43; j++){  
 //        printf(" %.f",(*matrix[i])[ii][j]);      
 //      };  
 //      printf("\n");  
 //       };  
 //       printf("\n");  
 //     };  
       
   
     cf->WriteLongMatrix(matrix[i],i);  
       
     if ( matrix[i]->Determinant() == 0. ){  
       printf("\n");  
       for (Int_t ii=0; ii<43; ii++){  
         for (Int_t j=0; j<43; j++){  
           printf(" %.f",(*matrix[i])[ii][j]);      
         };  
         printf("\n");  
       };  
       printf("\n");  
       printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i);  
684      } else {      } else {
685        Double_t det = 0.;        //
686        TMatrixD invmatrix = (TMatrixD)(matrix[i]->Invert(&det));        // FULL
687        printf(" Bin %i determinant is %f \n",i,det);        //
688        cf->WriteInvertedLongMatrix((TMatrixD)invmatrix,i);        for (Int_t i=0; i<nbin-1; i++){
689            //
690            // determine the average matrix
691            //    
692            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              };
700            };
701            //
702            cf->WriteFullMatrix(fmatrix[i],i);
703            //
704            if ( fmatrix[i]->Determinant() == 0. ){
705              printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i);
706            } else {
707              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        //
715        printf(" done, closing file and exiting\n");
716        //
717    };    };
718        //
   printf(" done, closing file and exiting\n");  
   
719    cf->CloseMatrixFile();    cf->CloseMatrixFile();
720      //
721    cf->Delete();    cf->Delete();
722      //
723  }  }
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.23