/[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.7 by mocchiut, Mon Jan 21 10:24:10 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];  Int_t MDIM = 8213;
29    //Int_t MDIM = 4128;
30  TArrayF *qplane[18];  //Float_t qqplane[17][43];
31  TArrayI *nqplane[18];  //Int_t nnqplane[17][43];
32  TMatrixD *matrix[18];  
33  TMatrixD *nmat[18];  TArrayF *qplane[17];
34    TArrayI *nqplane[17];
35    TMatrixD *matrix[17];
36    TMatrixD *nmat[17];
37    
38    //TMatrixD *fqplane;
39    //TMatrixD *fnqplane;
40    TMatrixD *fqplane[17];
41    TMatrixD *fnqplane[17];
42    TMatrixD *fmatrix[17];
43    //TMatrixF *fnmat[17];
44    //TMatrixD *fmatrix;
45    //TMatrixD *fnmat;
46    //TMatrixF *fmatrix;
47    //TMatrixF *fnmat;
48    TMatrixF *fnmat[17];
49    //Int_t finmat[43][191];
50    
51  //===============================================================================  //===============================================================================
52  bool Select( PamLevel2* event ){  bool Select( PamLevel2* event ){
# Line 40  bool Select( PamLevel2* event ){ Line 62  bool Select( PamLevel2* event ){
62    // tracker pre-selection    // tracker pre-selection
63    //------------------------------------------------------------------    //------------------------------------------------------------------
64    TrkTrack *trk = track->GetTrkTrack();    TrkTrack *trk = track->GetTrkTrack();
65      float rigidity = trk->GetRigidity();
66      if ( CRIG ) rigidity = event->GetCaloLevel2()->qtot/260.;
67      if ( SRIG ) rigidity = event->GetGPamela()->P0;
68    bool TRACK__OK = false;    bool TRACK__OK = false;
69    if(    if(
70       trk->chi2 >0     &&       trk->chi2 >0     &&
# Line 66  bool Select( PamLevel2* event ){ Line 90  bool Select( PamLevel2* event ){
90       !event->GetAcLevel2()->CARDhit() &&       !event->GetAcLevel2()->CARDhit() &&
91       !event->GetAcLevel2()->CAThit() &&       !event->GetAcLevel2()->CAThit() &&
92       true ) TOF__OK = true;       true ) TOF__OK = true;
93    if( !TOF__OK )return false;    if( !TOF__OK && !SIMU)return false;
94    //------------------------------------------------------    //------------------------------------------------------
95    // no albedo    // no albedo
96    //------------------------------------------------------    //------------------------------------------------------
97    if(    track->GetToFTrack()->beta[12]<=0.2 ||    if(    !SIMU && (track->GetToFTrack()->beta[12]<=0.2 ||
98           track->GetToFTrack()->beta[12] >= 1.5 ) return false;                     track->GetToFTrack()->beta[12] >= 1.5) ) return false;
99    
100    //------------------------------------------------------    //------------------------------------------------------
101    bool CUT1 = false;    bool CUT1 = false;
102    if(    if(
103       trk->nstep<100   &&       trk->nstep<100   &&
104       trk->GetRigidity()<400.   &&       rigidity<400.   &&
105       trk->GetRigidity()>0.1   &&       rigidity>0.1   &&
      trk->GetDeflection()<0.   &&  
106       trk->resx[0]<0.001 &&       trk->resx[0]<0.001 &&
107       trk->resx[5]<0.001 &&       trk->resx[5]<0.001 &&
108       track->IsSolved() &&       track->IsSolved() &&
# Line 87  bool Select( PamLevel2* event ){ Line 110  bool Select( PamLevel2* event ){
110       true ) CUT1 = true;       true ) CUT1 = true;
111    //------------------------------------------------------    //------------------------------------------------------
112    if( !CUT1 )return false;    if( !CUT1 )return false;
113      if ( trk->GetDeflection()>0. && !SIMU ) return false;
114    
115      //
116      //  ELENA'S CUT
117      //
118      //
119      // lever-arm 6
120      //====================================================
121      bool LX6=false;
122      if(
123         track->GetTrkTrack()->GetLeverArmX()==6 &&
124         !track->GetTrkTrack()->IsBad(0,0) &&
125         !track->GetTrkTrack()->IsBad(5,0) &&
126         track->GetTrkTrack()->resx[0]<0.001 &&
127         track->GetTrkTrack()->resx[5]<0.001 &&
128         track->GetTrkTrack()->IsInsideCavity() &&
129         true ) LX6 = true;
130      
131      //====================================================
132      // lever-arm 5
133      //====================================================
134      bool LX5=false;
135      if(
136         track->GetTrkTrack()->GetLeverArmX()==5 &&
137         true ){
138        if(
139           track->GetTrkTrack()->XGood(0) && track->GetTrkTrack()->XGood(4)
140           ){
141          
142          if(
143             !track->GetTrkTrack()->IsBad(0,0) &&
144             !track->GetTrkTrack()->IsBad(4,0) &&
145             track->GetTrkTrack()->resx[0]<0.001 &&
146             track->GetTrkTrack()->resx[4]<0.001 &&
147             track->GetTrkTrack()->IsInsideCavity() &&
148             true) LX5 = true;        
149        }else if (
150                  track->GetTrkTrack()->XGood(1) && track->GetTrkTrack()->XGood(5)
151                  ){
152          
153          if(
154             !track->GetTrkTrack()->IsBad(1,0) &&
155             !track->GetTrkTrack()->IsBad(5,0) &&
156             track->GetTrkTrack()->resx[1]<0.001 &&
157             track->GetTrkTrack()->resx[5]<0.001 &&
158             track->GetTrkTrack()->IsInsideCavity() &&
159             true) LX5 = true;
160        }
161      }
162      if ( !LX5 && !LX6 ) return false;
163    Float_t defl = trk->GetDeflection();    Float_t defl = trk->GetDeflection();
164    float p0 = 1.111588e+00;    float p0 = 1.111588e+00;
165    float p1 = 1.707656e+00;    float p1 = 1.707656e+00;
# Line 104  bool Select( PamLevel2* event ){ Line 177  bool Select( PamLevel2* event ){
177       track->GetTrkTrack()->chi2 < chi2m &&       track->GetTrkTrack()->chi2 < chi2m &&
178       true ) CUT2 = true;         true ) CUT2 = true;  
179    if ( !CUT2 ) return false;    if ( !CUT2 ) return false;
180      float dedxtrk    = trk->GetDEDX();
181      //  float zcutn =  9.   + 20./(rigidity*rigidity);
182      float zcut2 =  3.   + 4.3/(rigidity*rigidity);
183      float zcut1 =  0.52 + 0.455/(rigidity*rigidity);
184      Bool_t Z1 = false;
185      if(dedxtrk > zcut1 && dedxtrk < zcut2){
186        Z1=true;
187      }
188      if ( !Z1 && !SIMU ) return false;
189    //------------------------------------------------------    //------------------------------------------------------
190    //    //
191    // energy momentum match    // energy momentum match
# Line 115  bool Select( PamLevel2* event ){ Line 197  bool Select( PamLevel2* event ){
197    //    //
198    for (Int_t i=0; i < 22; i++){    for (Int_t i=0; i < 22; i++){
199      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 ){
200          return false;        return false;
201        };
202      };
203      //
204      if ( event->GetCaloLevel2()->qtot == 0. ) return false;
205      if ( rigidity>5. && track->GetCaloTrack()->qtrack/event->GetCaloLevel2()->qtot < 0.4 ) return false;
206      if ( rigidity<1. && track->GetToFTrack()->beta[12] < 0.8 ) return false;
207      if ( rigidity>50. ){
208        if ( trk->GetNX()<5  &&
209             trk->GetNY()<4 ) return false;
210        //
211        Bool_t sphit = false;
212        for ( Int_t plane = 0; plane < 6; plane++){
213          if ( !trk->XGood(plane) ){
214            for (Int_t sing = 0; sing < event->GetTrkLevel2()->nclsx(); sing++){
215              TClonesArray &t = *(event->GetTrkLevel2()->SingletX);
216              TrkSinglet *singlet = (TrkSinglet*)t[sing];
217              if ( (singlet->plane-1) == plane ){
218                Float_t x = (singlet->coord[0]+singlet->coord[1])/2.;
219                if ( fabs(track->GetTrkTrack()->xv[plane] - x) < 1. ) sphit = true;
220              };
221            };
222          };
223          if ( !trk->YGood(plane) ){
224            for (Int_t sing = 0; sing < event->GetTrkLevel2()->nclsy(); sing++){
225              TClonesArray &t = *(event->GetTrkLevel2()->SingletY);
226              TrkSinglet *singlet = (TrkSinglet*)t[sing];
227              if ( (singlet->plane-1) == plane ){
228                Float_t x1 = (singlet->coord[0]);
229                Float_t x2 = (singlet->coord[1]);
230                if ( fabs(track->GetTrkTrack()->yv[plane] - x1) < 1. ) sphit = true;
231                if ( fabs(track->GetTrkTrack()->yv[plane] - x2) < 1. ) sphit = true;
232              };
233            };
234          };
235      };      };
236        if ( sphit ) return false; // spurious hit along the track
237    };    };
238    //    //
239    if ( event->GetCaloLevel1()->qtotpl(0) > 7. ) return false;    Int_t ti0 = track->GetCaloTrack()->tibar[0][1]-1;
240      
241      Int_t view = 0;
242      Int_t plane = 0;
243      Int_t strip = 0;
244      Float_t mip = 0.;
245      //
246      for ( Int_t i=0; i<event->GetCaloLevel1()->istrip; i++ ){
247        //
248        mip = event->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
249        if ( view == 1 && plane == 0 && strip == ti0 && mip > 4.) return false;
250        if ( view == 1 && (plane >0 || strip > ti0) ) break;
251      };
252      //  if ( event->GetCaloLevel1()->qtotpl(0) > 7. ) return false;
253    
254      //  if ( rigidity > 2.2 || rigidity < 1.5 ) return false;
255      //  printf(" rig %f CRIG %i SRIG %i \n",rigidity,CRIG,SRIG);
256    //    //
257    return true;    return true;
258  }  }
# Line 134  bool Select( PamLevel2* event ){ Line 267  bool Select( PamLevel2* event ){
267  void CreateHistos( PamLevel2* event , TString file){  void CreateHistos( PamLevel2* event , TString file){
268    
269    cf = new CaloFranzini(event);    cf = new CaloFranzini(event);
270    cf->CreateMatrixFile(file.Data());    //
271      if ( MATRIX ){
272        cf->UpdateMatrixFile(file.Data());
273        cf->LoadBin();
274        if ( !FULL ){
275          cf->LoadLong();
276        } else {
277          cf->LoadFull();
278        };
279      } else {
280        cf->CreateMatrixFile(file.Data());
281      };
282      //
283    //    //
284    nbin = 18;    nbin = 18;
285    rig[0] = 0.1;    rig[0] = 0.1;
# Line 156  void CreateHistos( PamLevel2* event , TS Line 301  void CreateHistos( PamLevel2* event , TS
301    rig[16] = 200.;    rig[16] = 200.;
302    rig[17] = 400.;    rig[17] = 400.;
303    //    //
304  //  memset(qplane, 0, 44*18*sizeof(Float_t));    memset(rmean, 0, 17*sizeof(Float_t));
305    memset(rmean, 0, 18*sizeof(Float_t));    memset(ntot, 0, 17*sizeof(Int_t));
306    memset(qqplane, 0, 43*18*sizeof(Float_t));    //  memset(finmat, 0, 43*191*sizeof(Int_t));  
307    memset(nnqplane, 0, 43*18*sizeof(Float_t));  //  Double_t tol = 1E-20;
308    //    //
309    for (Int_t i=0; i < 18 ; i++){    for (Int_t i=0; i < 17 ; i++){
310      matrix[i] = new TMatrixD(43,43);      //  for (Int_t i=3; i < 4 ; i++){
311      qplane[i] = new TArrayF(43);      if ( !FULL ){
312      nqplane[i] = new TArrayI(43);        matrix[i] = new TMatrixD(43,43);
313      nmat[i] = new TMatrixD(43,43);        qplane[i] = new TArrayF(43);
314    };        nqplane[i] = new TArrayI(43);
315  //  for (Int_t i=0; i < 18 ; i++){        nmat[i] = new TMatrixD(43,43);
316  //      for (Int_t j = 0; j < 43;i++){      } else {
317  //              (*qplane[i])[j] = 0.;        if ( MATRIX ){
318  //              qplane[i]->Reset();          //      fmatrix = new TMatrixF(4128,4128);
319  //              nqplane[i]->Reset();          //      fnmat = new TMatrixF(4128,4128);
320  //              for (Int_t k = 0; k < 43;k++){          //      fmatrix = new TMatrixF(8213,8213);
321  //                      (*matrix[i])[j][k] = 0.;          //      fnmat = new TMatrixF(8213,8213);
322  //                      (*nmat[i])[j][k] = 0.;          //      fmatrix = new TMatrixF(MDIM,MDIM);
323  //              };          //      fnmat = new TMatrixF(MDIM,MDIM);
324  //      };          //      fmatrix[i] = new TMatrixF(1849,1849);
325  //  };          //      fnmat[i] = new TMatrixF(43,43);
326            fmatrix[i] = new TMatrixD(1333,1333);
327    //      fmatrix[i]->SetTol(tol);
328            fnmat[i] = new TMatrixF(43,31);
329            //      cf->WriteFullMatrix(fmatrix, i);
330            //      cf->WriteFullNMatrix(fnmat, i);
331            //      delete fmatrix;
332            //      delete fnmat;
333            //fnmat[i] = new TMatrixI(8213,8213);
334          } else {
335            //      fqplane = 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)
336            //      fnqplane = new TMatrixD(43,191);//
337            //      fqplane[i] = new TMatrixD(43,43); // 43 planes x 43 "strip", where 43 = 50 + 14 + 5 + 2 + 2 + 2 + 2 + 2 + 2 + 2 + 2 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + [1] + ...
338            //      fnqplane[i] = new TMatrixD(43,43);//                                     0    1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20   21
339            //
340            fqplane[i] = new TMatrixD(43,31); // 43 planes x 43 "strip", where 43 = 50 + 14 + 6 + 5 + 3 + 3 + 3 + 2 + 2 + 2 + 1 + 1 + 1 + 1 + 1 + [1] + 1 + 1 + 1 + 1 + ...
341            fnqplane[i] = new TMatrixD(43,31);//                                     0    1   2   3   4   5   6   7   8   9  10  11  12  13  14    15  16  17  18  19  ...
342            //
343            //      cf->WriteFullMean(fqplane, i);
344            //      cf->WriteFullNMean(fnqplane, i);
345            //      delete fqplane;
346            //      delete fnqplane;
347            //
348          };
349        };
350      };
351    //    //
352  }  }
353    
354  //===============================================================  //===============================================================
355  void FindAverage( PamLevel2* L2, int iev ){  void FindAverage( PamLevel2* L2, int iev ){
356    //    //
   //  printf("SELEZIONATO \n");  
357    Float_t erig = L2->GetTrack(0)->GetTrkTrack()->GetRigidity();    Float_t erig = L2->GetTrack(0)->GetTrkTrack()->GetRigidity();
358      if ( SRIG ) erig = L2->GetGPamela()->P0;
359      if ( CRIG ) erig = L2->GetCaloLevel2()->qtot/260.;
360    //    //
361    Int_t rbi = 0;    Int_t rbi = 0;
362    for (Int_t i = 0; i<nbin-1; i++){ // < 17    for (Int_t i = 0; i<nbin-1; i++){
363      if ( erig>=rig[i] && erig < rig[i+1] ){      if ( erig>=rig[i] && erig < rig[i+1] ){
364        rbi = i;        rbi = i;
365        break;        break;
366      };      };
367    };    };
368    //    //
   //  printf("SELEZIONATO rbi %i erig %f rig[0] %f rig[nbin] %f \n",rbi, erig, rig[0],rig[nbin-1]);  
   //  
369    if ( erig < rig[0] ) return;    if ( erig < rig[0] ) return;
370    if ( erig >= rig[nbin-1] ) return;    if ( erig >= rig[nbin-1] ) return;
371    //    //
372    rmean[rbi] += erig;    rmean[rbi] += erig;
373      ntot[rbi]++;
374    //    //
375    Int_t dgf = 43;    if (!FULL ){
376    //      Int_t dgf = 43;
377  //  for (Int_t i=0; i < 22; i++){      //
378  //    if ( L2->GetTrack(0)->GetCaloTrack()->tibar[i][1] < 0 ){      for (Int_t i=0; i<dgf; i++){
379  //      dgf = 2 * i;        (*nqplane[rbi])[i]++;
380   //     break;      };
381    //  };      //
382   //   if ( L2->GetTrack(0)->GetCaloTrack()->tibar[i][0] < 0 ){      // Fill the estrip matrix
383   //     dgf = 1 + 2 * i;      //
384   //     break;      Int_t nplane = 0;
385   //   };      Int_t view = 0;
386   // };      Int_t plane = 0;
387    //      Int_t strip = 0;
388  //  if ( dgf < 43 && dgf > 37 ) dgf--;      Float_t mip = 0.;
389    //      //
390    for (Int_t i=0; i<dgf; i++){      for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
391      (*nqplane[rbi])[i]++;        //
392          nnqplane[rbi][i]++;        mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
393    };        //
394    //        nplane = 1 - view + 2 * plane;
395    //  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());
396    //        //printf(" IEV %i OBT %u pkt %u file %s \n",iev,L2->GetOrbitalInfo()->OBT,L2->GetOrbitalInfo()->pkt_num,L2->GetPamTree()->GetFile()->GetName());
397    //        if ( nplane > 37 ) nplane--;
398    // Fill the estrip matrix        if ( nplane < dgf ){
399    //          (*qplane[rbi])[nplane] += mip;
400    Int_t nplane = 0;        };
401    Int_t view = 0;        //
402    Int_t plane = 0;      };
403    Int_t strip = 0;    } else {
404    Float_t mip = 0.;      //
405    //      // FULL CALORIMETER
   for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){  
406      //      //
407      mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);      //    fqplane = cf->LoadFullAverage(rbi);
408        //    fnqplane = cf->LoadFullNAverage(rbi);
409        CaloTrkVar *ct = L2->GetTrack(0)->GetCaloTrack();
410      //      //
411      nplane = 1 - view + 2 * plane;      Int_t nplane = 0;
412      if ( nplane > 37 ) nplane--;      Int_t view = 0;
413      //    printf(" mip %f view %i plane %i => nplane %i \n",mip,view,plane,nplane);      Int_t plane = 0;
414      if ( nplane < dgf ){      Int_t strip = 0;
415        //      printf(" PRIMA: qpl %f nqpl %i \n",(*qplane[rbi])[nplane],(*nqplane[rbi])[nplane]);      Float_t mip = 0.;
416        (*qplane[rbi])[nplane] += mip;      //
417          qqplane[rbi][nplane] += mip;      Int_t cs = 0;
418        //      printf(" DOPO: qpl %f nqpl %i \n",(*qplane[rbi])[nplane],(*nqplane[rbi])[nplane]);      Int_t cd = 0;
419        Int_t mstrip = 0;
420        //
421        for (Int_t j=0; j<2; j++){
422          for (Int_t i=0; i<22; i++){
423            nplane = 1 - j + 2*i;
424            if ( nplane > 37 ) nplane--;
425            //
426            cs = ct->tibar[i][j] - 1;
427            //
428            cd = 95 - cs;  
429            //
430            Int_t oldstr = -1;
431            for (Int_t k=0; k<191; k++){
432              mstrip = cd + k;
433              //      if ( mstrip < (191-cs) ) (*fnqplane[rbi])[nplane][mstrip] += 1.;
434              //      if ( mstrip < (191-cs) ) (*fnqplane)[nplane][mstrip] += 1.;
435              Int_t lstr = cf->ConvertStrip(mstrip);
436              if ( oldstr != lstr ){
437                (*fnqplane[rbi])[nplane][lstr] += 1.;
438                oldstr = lstr;
439              };    
440            };
441          };
442      };      };
443      //      //
444        //
445        for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
446          //
447          mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
448          //
449          nplane = 1 - view + 2 * plane;
450          if ( nplane > 37 ) nplane--;
451          //
452          cs = ct->tibar[plane][view] - 1;
453          //
454          cd = 95 - cs;
455          //
456          mstrip = cd + strip;
457          //
458          Int_t lstr = cf->ConvertStrip(mstrip);
459          //      (*fqplane[rbi])[nplane][mstrip] += mip;
460          //      (*fqplane)[nplane][mstrip] += mip;
461          (*fqplane[rbi])[nplane][lstr] += mip;
462          //
463        };
464        //
465        //    cf->WriteFullMean(fqplane, rbi);
466        //    cf->WriteFullNMean(fnqplane, rbi);
467        //    cf->UnLoadFullAverage(rbi);
468        //    cf->UnLoadFullNAverage(rbi);
469        //    delete fqplane;
470        //    delete fnqplane;
471        //
472    };    };
473  }  }
474    
475    void CalculateAverage(){
476      //  
477      if ( !FULL ){
478        for (Int_t i=0; i<nbin-1; i++){
479          if ( (*nqplane[i])[0] > 0 ) rmean[i] /= (Float_t)(*nqplane[i])[0];
480          for (Int_t j=0; j<43 ; j++){
481            if ( (*nqplane[i])[j] > 0 ){  
482              (*qplane[i])[j] /= (Float_t)(*nqplane[i])[j];
483            } else {
484              (*qplane[i])[j] = 0.;
485            };
486            printf(" BIN %i plane %i average energy %f qplane %f nqplane %i \n",i,j,rmean[i],(*qplane[i])[j],(*nqplane[i])[j]);
487          };
488        };
489        for (Int_t i=0; i<nbin-1; i++){
490          //
491          cf->WriteLongMean(qplane[i], i);
492          //
493        };
494      } else {
495        //
496        for (Int_t i=0; i<nbin-1; i++){
497          //      fqplane = cf->LoadFullAverage(i);
498          //      fnqplane = cf->LoadFullNAverage(i);
499          if ( ntot[i] > 0 ) rmean[i] /= (Float_t)(ntot[i]);
500          //
501          for (Int_t j=0; j<43 ; j++){
502            //      for (Int_t k=0; k<191; k++){
503            //      for (Int_t k=0; k<43; k++){
504            for (Int_t k=0; k<31; k++){
505              //      if ( (*fnqplane[i])[j][k] > 0 ){  
506              //        (*fqplane[i])[j][k] /= (Float_t)(*fnqplane[i])[j][k];
507              //      } else {
508              //        (*fqplane[i])[j][k] = 0.;
509              //      };
510              //      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]);
511              if ( (*fnqplane[i])[j][k] > 0 ){  
512                if ( (*fqplane[i])[j][k] == 0. ) (*fqplane[i])[j][k] = 0.7;
513                (*fqplane[i])[j][k] /= (Float_t)(*fnqplane[i])[j][k];
514              } else {
515                (*fqplane[i])[j][k] = 0.;
516              };
517              //      printf(" BIN %i plane %i strip %i average energy %f qplane %f nqplane %f \n",i,j,k,rmean[i],(*fqplane)[j][k],(*fnqplane)[j][k]);
518            };
519          };
520          cf->WriteFullMean(fqplane[i], i);
521          cf->WriteFullNMean(fnqplane[i], i);
522          //      cf->UnLoadFullAverage(i);
523          //      cf->UnLoadFullNAverage(i);
524          //      delete fqplane;
525          //      delete fnqplane;
526        };
527        //
528        //     for (Int_t i=0; i<nbin-1; i++){
529        //       //
530        //       cf->WriteFullMean(fqplane[i], i);
531        //       //
532        //     };
533      };
534      //  
535      cf->WriteNumBin(nbin);
536      //
537      TArrayF *rigbin = new TArrayF(18, rig);
538      cf->WriteRigBin(rigbin);
539      //
540      TArrayF *rmeanbin = new TArrayF(17, rmean);
541      TFile *file =  cf->GetFile();
542      file->cd();
543      file->WriteObject(&(*rmeanbin), "binrigmean");
544      //
545      //
546    }
547    
548  //===============================================================  //===============================================================
549  void FindMatrix( PamLevel2* L2, int iev ){  void FindMatrix( PamLevel2* L2, int iev ){
550    //    //
   //  printf("SELEZIONATO \n");  
551    Float_t erig = L2->GetTrack(0)->GetTrkTrack()->GetRigidity();    Float_t erig = L2->GetTrack(0)->GetTrkTrack()->GetRigidity();
552      if ( SRIG ) erig = L2->GetGPamela()->P0;
553      if ( CRIG ) erig = L2->GetCaloLevel2()->qtot/260.;
554    //    //
555    Int_t rbi = 0;    Int_t rbi = 0;
556    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 560  void FindMatrix( PamLevel2* L2, int iev
560      };      };
561    };    };
562    //    //
   //  printf("SELEZIONATO rbi %i erig %f rig[0] %f rig[nbin] %f \n",rbi, erig, rig[0],rig[nbin-1]);  
   //  
563    if ( erig < rig[0] ) return;    if ( erig < rig[0] ) return;
564    if ( erig >= rig[nbin-1] ) return;    if ( erig >= rig[nbin-1] ) return;
565    //    //
566    Int_t dgf = 43;    if ( !FULL ){
567    //      Int_t dgf = 43;
568  //  for (Int_t i=0; i < 22; i++){      //
569  //    if ( L2->GetTrack(0)->GetCaloTrack()->tibar[i][1] < 0 ){      for (Int_t i=0; i<dgf; i++){
570  //      dgf = 2 * i;        for (Int_t j=0; j<dgf; j++){
571   //     break;          (*nmat[rbi])[i][j] += 1.;
572   //   };        };
  //   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.;  
573      };      };
   };  
   //  
   //  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++ ){  
574      //      //
575      mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);      // Fill the estrip matrix
576      //      //
577      nplane = 1 - view + 2 * plane;      Int_t nplane = 0;
578      if ( nplane > 37 ) nplane--;      Int_t view = 0;
579      if ( nplane < dgf ){      Int_t plane = 0;
580        hpl[nplane] += mip;      Int_t strip = 0;
581        Float_t mip = 0.;
582        Float_t hpl[43];
583        memset(hpl,0,43*sizeof(Float_t));
584        for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
585          //
586          mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
587          //
588          nplane = 1 - view + 2 * plane;
589          if ( nplane > 37 ) nplane--;
590          if ( nplane < dgf ){
591            hpl[nplane] += mip;
592          };
593          //
594      };      };
595      //      //
596    };      for (Int_t i=0; i<dgf; i++){
597    //        for (Int_t j=0; j<dgf; j++){
598    for (Int_t i=0; i<dgf; i++){          (*matrix[rbi])[i][j] += (hpl[i] - cf->GetAverageAt(i,erig)) * (hpl[j] - cf->GetAverageAt(j,erig));
599      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]);  
600      };      };
601    };    } else {
602  }      //
603        // FULL CALORIMETER
604  void CalculateAverage(){      //
605    for (Int_t i=0; i<nbin; i++){      //    if ( rbi != 3 ) return;
606      if ( (*nqplane[i])[0] > 0 ) rmean[i] /= (Float_t)(*nqplane[i])[0];      printf(" matrix %i IEV %i \n",rbi,iev);
607      for (Int_t j=0; j<43 ; j++){      //    fmatrix = cf->LoadFullMatrix(rbi);
608        if ( (*nqplane[i])[j] > 0 ){        //    cf->LoadFullMatrix(rbi,fmatrix);
609          (*qplane[i])[j] /= (Float_t)(*nqplane[i])[j];      //    fnmat = cf->LoadFullNMatrix(rbi);
610        } else {      //    printf(" done \n");
611          (*qplane[i])[j] = 0.;      //    printf(" start loop \n");
612        //
613        CaloTrkVar *ct = L2->GetTrack(0)->GetCaloTrack();
614        //
615        Int_t nplane = 0;
616        Int_t view = 0;
617        Int_t plane = 0;
618        Int_t strip = 0;
619        Float_t mip = 0.;
620        //
621        //    Int_t mindgf = 48;
622        //    Int_t dgf = 143;
623        //       Int_t mindgf = 0; //tutto
624        //        Int_t dgf = 191; //tutto
625        //    Int_t mindgf = 94;
626        //    Int_t dgf = 96;
627        //    Int_t mindgf = 84;
628        //    Int_t dgf = 106;
629        Int_t mindgf = 0;
630        Int_t dgf = 43;
631        Int_t cs = 0;
632        Int_t cd = 0;
633        Int_t mstrip = 0;
634        //
635        //    Float_t mipv[43][43];
636        //    memset(mipv,0,43*43*sizeof(Float_t));
637        Float_t mipv[43][31];
638        memset(mipv,0,43*31*sizeof(Float_t));
639        //
640        for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
641          //
642          mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
643          //
644          nplane = 1 - view + 2 * plane;
645          if ( nplane > 37 ) nplane--;
646          //
647          cs = ct->tibar[plane][view] - 1;
648          //
649          cd = 95 - cs;
650          //
651          mstrip = cd + strip;
652          //
653          Int_t lstr = cf->ConvertStrip(mstrip);
654          mipv[nplane][lstr] += mip;
655          //
656        };
657        //
658        Float_t mip1 = 1.;
659        Int_t cs1;
660        Int_t cd1;
661        Float_t mip2 = 1.;
662        Int_t cs2;
663        Int_t cd2;
664        Int_t mi = -1;
665        Int_t mj = -1;
666        Int_t nn1 = 0;
667        Int_t pl1 = 0;
668        Int_t vi1 = 0;
669        Int_t nn2 = 0;
670        Int_t pl2 = 0;
671        Int_t vi2 = 0;
672        Int_t mstrip1min = 0;
673        Int_t mstrip1max = 0;
674        Int_t mstrip2min = 0;
675        Int_t mstrip2max = 0;
676        //
677        Int_t toto = 0;
678        //
679        for (Int_t nplane1=0; nplane1<43; nplane1++){
680          if ( nplane1 >= 37 ) nn1 = nplane1 + 1;
681          vi1 = 1;
682          if ( nn1%2 ) vi1 = 0;
683          pl1 = (nn1 - 1 + vi1)/2;
684          //
685          cs1 = ct->tibar[pl1][vi1] - 1; // convertire nplane in pl1 e vi1
686          //
687          cd1 = 95 - cs1;
688          //
689          Int_t at1 = TMath::Max(0,(cd1+0));
690          Int_t at2 = TMath::Min(190,(cd1+95));
691          mstrip1min = cf->ConvertStrip(at1);
692          mstrip1max = cf->ConvertStrip(at2) + 1;
693          //      mstrip1min = cf->ConvertStrip(TMath::Max(mindgf,(cd1+0)));
694          //      mstrip1max = cf->ConvertStrip(TMath::Min(dgf,(cd1+95))) + 1;
695          //
696          if ( nplane1 == 0 || nplane1 == 42 ) printf(" pl %i mstrip1min %i mstrip1max %i mindgf %i dgf %i cd1 %i\n",nplane1,mstrip1min,mstrip1max,mindgf,dgf,cd1);
697          //
698          for (Int_t mstrip1=mstrip1min; mstrip1<mstrip1max; mstrip1++){
699            //      printf(".\n");
700            //
701            mj = -1;
702            //
703            mip1 = mipv[nplane1][mstrip1] - cf->GetFullAverageAt(nplane1,mstrip1,erig,rbi);
704            //
705            //      mi = (nplane1 * 191) + mstrip1;
706            //      mi = (nplane1 * 43) + mstrip1;
707            mi = (nplane1 * 31) + mstrip1;
708            //
709            //      if ( mstrip1 > mstrip1min ) break;
710            //      if ( mstrip1 > dgf ) break;
711            //      if ( mstrip1 >= mindgf && mstrip1 <= dgf && mstrip1 >= mstrip1min && mstrip1 <= mstrip1max ){
712            //
713            //      finmat[nplane1][mstrip1]++;
714            (*fnmat[rbi])[nplane1][mstrip1] += 1.;
715            //
716            if ( mip1 != 0. ){
717              //
718              for (Int_t nplane2=0; nplane2<43; nplane2++){
719                //
720                if ( nplane2 >= 37 ) nn2 = nplane2 + 1;
721                vi2 = 1;
722                if ( nn2%2 ) vi2 = 0;
723                pl1 = (nn2 - 1 + vi2)/2;
724                //
725                cs2 = ct->tibar[pl2][vi2] - 1;
726                //
727                cd2 = 95 - cs2;
728                //
729                //    mstrip2min = cd2 + 0;
730                //    mstrip2max = cd2 + 95;
731                Int_t t1 = TMath::Max(0,(cd2+0));
732                Int_t t2 = TMath::Min(190,(cd2+95));
733                mstrip2min = cf->ConvertStrip(t1);
734                mstrip2max = cf->ConvertStrip(t2) + 1;
735                //
736                if ( nplane1 == 0 && nplane2 == 0 && mstrip1==mstrip1min ) printf(" mstrip2min %i mstrip2max %i \n",mstrip2min,mstrip2max);
737                //
738                for (Int_t mstrip2=mstrip2min; mstrip2<mstrip2max; mstrip2++){
739                  //              
740                  mip2 = mipv[nplane2][mstrip2] - cf->GetFullAverageAt(nplane2,mstrip2,erig,rbi);
741                  //
742                  if ( mip2 != 0. ){
743                    //
744                    //              mj = (nplane2 * 191) + mstrip2;
745                    //              mj = (nplane2 * 43) + mstrip2;
746    //              mj = (nplane2 * 31) + mstrip2;
747                    Int_t sh = -15 + nplane1;
748                    if ( sh > 15 ) sh -= 31*nplane1;
749                    //
750                    mj = (nplane2 * 31) + mstrip2 + sh;
751                    //
752                    if ( mj < 0 ) mj += 1333;
753                    if ( mj >= 1333 ) mj -= 1333;
754    //              printf(" mi %i mj %i sh %i \n",mi,mj,sh);
755                    //
756                    //            mj++;
757                    //
758                    //            if ( mstrip2 > mstrip2min ) break;
759                    //            if ( mstrip2 > dgf ) break;
760                    //            if ( mstrip2 >= mindgf && mstrip2 <= dgf && mstrip2 >= mstrip2min && mstrip2 <= mstrip2max ){
761                    //              if ( mstrip1 >= mstrip1min && mstrip1 <= mstrip1max && mstrip2 >= mstrip2min && mstrip2 <= mstrip2max){
762                    //              (*fmatrix[rbi])[mi][mj] += (mipv[nplane1][mstrip1] - cf->GetFullAverageAt(nplane1,mstrip1,erig)) * (mipv[nplane2][mstrip2] - cf->GetFullAverageAt(nplane2,mstrip2,erig));
763                    //                (*fnmat[rbi])[mi][mj] += 1.;
764                    (*fmatrix[rbi])[mi][mj] += (mip1 * mip2); // giusto
765                    //              (*fmatrix)[mi][mj] += (mip1 * mip2) * 1000000.;
766                    toto++;
767                    //            (*fmatrix)[mi][mj] += 1.;
768                    //              cf->GetFullAverageAt(nplane1,mstrip1,erig,rbi);
769                    //              cf->GetFullAverageAt(nplane2,mstrip2,erig,rbi);
770                    //                (*fnmat)[mi][mj] += 1.;
771                    //              };
772                  };
773                };
774              };
775            };
776        };        };
       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]);  
777      };      };
778    };      //  
779          printf(" toto = %i \n",toto);
780    cf->WriteNumBin(nbin);      printf("\n done \n");
781        //    printf(" write matrix \n");
782    TArrayF *rigbin = new TArrayF(18, rig);      //    cf->WriteFullMatrix(fmatrix, rbi);
783    cf->WriteRigBin(rigbin);      //    cf->WriteFullNMatrix(fnmat, rbi);
784        //    printf(" done \n");
785    TArrayF *rmeanbin = new TArrayF(18, rmean);      //    printf(" unload matrix \n");
786    TFile *file =  cf->GetFile();      //   cf->UnLoadFullMatrix(rbi);
787    file->cd();      //    cf->UnLoadFullNMatrix(rbi);
788    //  rigbin->Write("binrig");      //    printf(" done \n");
789    file->WriteObject(&(*rmeanbin), "binrigmean");      //    printf(" delete matrix \n");
790        //    delete fmatrix;
791  //  cf->WriteRigBin(rigbin);      //    delete fnmat;
792          //    printf(" done \n");
   for (Int_t i=0; i<nbin; i++){  
       
     cf->WriteLongMean(qplane[i], i);  
   
793    };    };
794  }  }
795    
796  //===============================================================  //===============================================================
797  // Save histograms  // Save histograms
798  //  //
# Line 370  void CalculateAverage(){ Line 802  void CalculateAverage(){
802  //  //
803  //===============================================================  //===============================================================
804  void SaveHistos(){  void SaveHistos(){
805        //
806    printf("Finished, calculating average and inverting matrices\n");    if ( MATRIX ){
807        //
808        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);  
   
809      //      //
810      // determine the average matrix      if ( !FULL ){
811      //            for (Int_t i=0; i<nbin-1; i++){
812      for (Int_t ii=0; ii<43; ii++){          //
813        for (Int_t j=0; j<43; j++){          // determine the average matrix
814          //      if ( i == 1 ) printf(" num %f \n",(*nmat[i])[ii][j]);          //    
815          if ( (*nmat[i])[ii][j] > 0. ){          for (Int_t ii=0; ii<43; ii++){
816            //      if ( i == 1 ) printf(" value %f \n",(*matrix[i])[ii][j]);            for (Int_t j=0; j<43; j++){
817            (*matrix[i])[ii][j] /= (*nmat[i])[ii][j];              if ( (*nmat[i])[ii][j] > 0. ){
818                  (*matrix[i])[ii][j] /= (*nmat[i])[ii][j];
819                } else {
820                  (*matrix[i])[ii][j] = 0.;
821                };
822              };
823            };
824            //
825            cf->WriteLongMatrix(matrix[i],i);
826            //
827            if ( matrix[i]->Determinant() == 0. ){
828              printf("\n");
829              for (Int_t ii=0; ii<43; ii++){
830                for (Int_t j=0; j<43; j++){
831                  printf(" %.f",(*matrix[i])[ii][j]);        
832                };
833                printf("\n");
834              };
835              printf("\n");
836              printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i);
837          } else {          } else {
838            (*matrix[i])[ii][j] = 0.;            Double_t det = 0.;
839              TMatrixD invmatrix = (TMatrixD)(matrix[i]->Invert(&det));
840              printf(" Bin %i determinant is %f \n",i,det);
841              cf->WriteInvertedLongMatrix((TMatrixD)invmatrix,i);
842          };          };
843        };        };
844      };      } else {
845              //
846  //     if ( i == 2 ){        // FULL
847  //       printf("\n");        //
848  //       for (Int_t ii=0; ii<43; ii++){        for (Int_t i=0; i<nbin-1; i++){
849  //      for (Int_t j=0; j<43; j++){          //      for (Int_t i=3; i<5; i++){
850  //        printf(" %.f",(*matrix[i])[ii][j]);              //      for (Int_t i=3; i<4; i++){
851            //
852            // determine the average matrix
853            //    
854            //      fmatrix = cf->LoadFullMatrix(i);
855            //      fnmat = cf->LoadFullNMatrix(i);
856            //
857            //      for (Int_t ii=0; ii<MDIM; ii++){
858            //        for (Int_t j=0; j<MDIM; j++){
859            // //       if ( (*fnmat[i])[ii][j] > 0. ){
860            // //         (*fmatrix[i])[ii][j] /= (*fnmat[i])[ii][j];
861            // //       } else {
862            // //         (*fmatrix[i])[ii][j] = 0.;
863            // //       };
864            //          if ( (*fnmat)[ii][j] > 0. ){
865            //            (*fmatrix)[ii][j] /= (*fnmat)[ii][j];
866            //          } else {
867            //            (*fmatrix)[ii][j] = 0.;
868            //          };
869            //        };
870            //      };
871            //
872    //      TMatrixD *mymat3 = new TMatrixD(129,129);
873    //      TMatrixD *mymat5 = new TMatrixD(215,215);
874    //      TMatrixD *mymat7 = new TMatrixD(301,301);
875    //      TMatrixD *mymat9 = new TMatrixD(387,387);
876    //      TMatrixD *mymat11 = new TMatrixD(473,473);
877    //      TMatrixD *mymat17 = new TMatrixD(731,731);
878            //      TMatrixF *mymat = new TMatrixF(129,129);
879            //      TMatrixF *mymat = new TMatrixF(989,989);
880            Int_t i1 = -1;
881            Int_t j1 = -1;
882            //      int mi,mj;
883            Int_t nonzero = 0;
884            Int_t nonzero1 = 0;
885            for (Int_t ii=0; ii<43; ii++){
886              //      for (Int_t j=0; j<191; j++){
887              //      for (Int_t j=0; j<43; j++){
888              for (Int_t j=0; j<31; j++){
889                //      if ( (*fnmat[i])[ii][j] > 0. ){
890                //        (*fmatrix[i])[ii][j] /= (*fnmat[i])[ii][j];
891                //      } else {
892                //        (*fmatrix[i])[ii][j] = 0.;
893                //      };
894                //      i1 =  (ii * 191) + j;
895                //      i1 =  (ii * 43) + j;
896                i1 =  (ii * 31) + j;
897                //      j1 = -1;
898                for (Int_t iij=0; iij<43; iij++){
899                  //              for (Int_t jj=0; jj<191; jj++){
900                  //              for (Int_t jj=0; jj<43; jj++){
901                  for (Int_t jj=0; jj<31; jj++){
902                    //
903                    //              j1 = (iij * 191) + jj;
904                    //              j1 = (iij * 43) + jj;
905                    Int_t sh = -15 + ii;
906                    if ( sh > 15 ) sh -= 31*ii;
907                    //
908                    j1 = (iij * 31) + jj + sh;
909                    //
910                    if ( j1 < 0 ) j1 += 1333;
911                    if ( j1 >= 1333 ) j1 -= 1333;
912    
913    //              j1 = (iij * 31) + jj;
914                    //              j1++;
915                    //              if ( finmat[ii][j] > 0 ){
916                    //                (*fmatrix)[i1][j1] /= finmat[ii][j];
917                    if ( (*fnmat[i])[ii][j] == 0. || (*fmatrix[i])[i1][j1] == 0. || !((*fmatrix[i])[i1][j1] == (*fmatrix[i])[i1][j1]) ){
918                      (*fmatrix[i])[i1][j1] = 1.;
919                    } else {
920                      (*fmatrix[i])[i1][j1] /= (*fnmat[i])[ii][j];
921                      nonzero++;
922                      if ( i1 == 0 ) nonzero1++;
923                    };
924                    //
925    //              if ( j>=7 && j <=23 && jj >=7 && jj<=23 ){
926    //                Int_t mi17 = (ii*3) + j -7;
927    //                Int_t mj17 = (iij*3) + jj -7;
928    //                (*mymat17)[mi17][mj17] = (*fmatrix[i])[i1][j1];
929    //              };
930    //              if ( j>=10 && j <=20 && jj >=10 && jj<=20 ){
931    //                Int_t mi11 = (ii*3) + j -10;
932    //                Int_t mj11 = (iij*3) + jj -10;
933    //                (*mymat11)[mi11][mj11] = (*fmatrix[i])[i1][j1];
934    //              };
935    //              if ( j>=11 && j <=19 && jj >=11 && jj<=19 ){
936    //                Int_t mi9 = (ii*3) + j -11;
937    //                Int_t mj9 = (iij*3) + jj -11;
938    //                (*mymat9)[mi9][mj9] = (*fmatrix[i])[i1][j1];
939    //              };
940    //              if ( j>=12 && j <=18 && jj >=12 && jj<=18 ){
941    //                Int_t mi7 = (ii*3) + j -12;
942    //                Int_t mj7 = (iij*3) + jj -12;
943    //                (*mymat7)[mi7][mj7] = (*fmatrix[i])[i1][j1];
944    //              };
945    //              if ( j>=13 && j <=17 && jj >=13 && jj<=17 ){
946    //                Int_t mi5 = (ii*3) + j -13;
947    //                Int_t mj5 = (iij*3) + jj -13;
948    //                (*mymat5)[mi5][mj5] = (*fmatrix[i])[i1][j1];
949    //              };
950    //              if ( j>=14 && j <=16 && jj >=14 && jj<=16 ){
951    //                Int_t mi3 = (ii*3) + j -14;
952    //                Int_t mj3 = (iij*3) + jj -14;
953    //                (*mymat3)[mi3][mj3] = (*fmatrix[i])[i1][j1];
954    //              };
955    
956    
957                    //              if ( j>=94 && j <=96 && jj >=94 && jj<=96 ){
958                    //                      mi = (ii*3) + j -94;
959                    //                      mj = (iij*3) + jj -94;
960                    //                      (*mymat)[mi][mj] = (*fmatrix)[i1][j1];
961                    //              };
962    
963    
964                    //              if ( j>=84 && j <=106 && jj >=84 && jj<=106 ){
965                    //                      mi = (ii*3) + j -84;
966                    //                      mj = (iij*3) + jj -84;
967                    //                      (*mymat)[mi][mj] = (*fmatrix)[i1][j1];
968                    //              };
969    
970                  };
971                };
972              };
973            };
974            //
975            printf(" Matrix has %i non-zero elements \n",nonzero);  
976            //      printf(" Matrix has %i non-zero elements on the first row\n",nonzero1);
977            //
978            //      Bool_t BAD = false;
979            //      for (Int_t ii=0; ii<43; ii++){
980            //        for (Int_t j=0; j<191; j++){
981            //          //
982            //          i1 =  (ii * 191) + j;
983            //          //
984            //          for (Int_t iij=0; iij<43; iij++){
985            //            for (Int_t jj=0; jj<191; jj++){
986            //              //
987            //              j1 = (iij * 191) + jj;
988            //              //
989            //              //              printf(" ROW %i COLUMN %i VALUE %f \n",i1,j1,(*fmatrix)[i1][j1]);
990            //              if ( (*fmatrix)[i1][j1] == 0. || !((*fmatrix)[i1][j1]==(*fmatrix)[i1][j1]) ){
991            //                printf(" ROW %i COLUMN %i VALUE %f \n",i1,j1,(*fmatrix)[i1][j1]);
992            //                printf(" che schifo! \n");
993            //                BAD = true;
994            //              };
995            //              //
996            //            };
997            //          };
998            //        };
999            //      };
1000            //      //
1001            //      if ( BAD ) printf(" questa matrice fa cagare \n");
1002            //
1003            //
1004            cf->WriteFullMatrix(fmatrix[i],i);
1005            //      cf->WriteFullMatrix(fmatrix, i);
1006            //      cf->WriteFullNMatrix(fnmat, i);
1007            cf->WriteFullNMatrix(fnmat[i], i);
1008            //
1009            //      TDecompSVD svd(*fmatrix[i]);
1010            //      Bool_t ok = svd.Decompose();
1011            //
1012            Double_t zero = (Double_t)0.0;
1013            //
1014            if ( fmatrix[i]->Determinant() == zero ){
1015              //if ( fmatrix->Determinant() == 0. ){
1016              printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i);
1017            } else {
1018              //    };
1019              //    if ( i == 3 ){
1020              //    if ( ok ){
1021              //      Double_t tol = 1E-20;
1022              //      TDecompSVD svd((*fmatrix)[i],tol);
1023              //      svd.Decompose();
1024              //      TMatrixD svdInv = svd.Invert();
1025              //      svdInv.Print("svdInv");
1026              //      cout << "condition: " << svd.Condition() << endl;
1027              //      cf->WriteInvertedFullMatrix((TMatrixD)svdInv,999);
1028              
1029              Double_t det = 0.;
1030              TMatrixD invmatrix = (TMatrixD)(fmatrix[i]->Invert(&det));
1031              printf(" Bin %i determinant is %f \n",i,det);
1032              cf->WriteInvertedFullMatrix((TMatrixD)invmatrix,i);
1033            };
1034    
1035    //      if ( mymat3->Determinant() == 0. ){
1036    //        printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i);
1037    //      } else {
1038    //        Double_t det = 0.;
1039    //        TMatrixD invmatrix = (TMatrixD)(mymat3->Invert(&det));
1040    //        printf(" Mymat3 determinant is %f \n",det);
1041    //        cf->WriteInvertedFullMatrix((TMatrixD)invmatrix,1103);
1042    //      };
1043    //      cf->WriteFullMatrix(mymat3, 103);
1044    //      if ( mymat5->Determinant() == 0. ){
1045    //        printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i);
1046    //      } else {
1047    //        Double_t det = 0.;
1048    //        TMatrixD invmatrix = (TMatrixD)(mymat5->Invert(&det));
1049    //        printf(" Mymat5 determinant is %f \n",det);
1050    //        cf->WriteInvertedFullMatrix((TMatrixD)invmatrix,1105);
1051  //      };  //      };
1052  //      printf("\n");  //      cf->WriteFullMatrix(mymat5, 105);
1053  //       };  //      if ( mymat7->Determinant() == 0. ){
1054  //       printf("\n");  //        printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i);
1055  //     };  //      } else {
1056        //        Double_t det = 0.;
1057    //        TMatrixD invmatrix = (TMatrixD)(mymat7->Invert(&det));
1058      cf->WriteLongMatrix(matrix[i],i);  //        printf(" Mymat7 determinant is %f \n",det);
1059        //        cf->WriteInvertedFullMatrix((TMatrixD)invmatrix,1107);
1060      if ( matrix[i]->Determinant() == 0. ){  //      };
1061        printf("\n");  //      cf->WriteFullMatrix(mymat7, 107);
1062        for (Int_t ii=0; ii<43; ii++){  //      if ( mymat9->Determinant() == 0. ){
1063          for (Int_t j=0; j<43; j++){  //        printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i);
1064            printf(" %.f",(*matrix[i])[ii][j]);      //      } else {
1065          };  //        Double_t det = 0.;
1066          printf("\n");  //        TMatrixD invmatrix = (TMatrixD)(mymat9->Invert(&det));
1067    //        printf(" Mymat3 determinant is %f \n",det);
1068    //        cf->WriteInvertedFullMatrix((TMatrixD)invmatrix,1109);
1069    //      };
1070    //      cf->WriteFullMatrix(mymat9, 109);
1071    //      if ( mymat11->Determinant() == 0. ){
1072    //        printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i);
1073    //      } else {
1074    //        Double_t det = 0.;
1075    //        TMatrixD invmatrix = (TMatrixD)(mymat11->Invert(&det));
1076    //        printf(" Mymat11 determinant is %f \n",det);
1077    //        cf->WriteInvertedFullMatrix((TMatrixD)invmatrix,1111);
1078    //      };
1079    //      cf->WriteFullMatrix(mymat11, 111);
1080    //      if ( mymat17->Determinant() == 0. ){
1081    //        printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i);
1082    //      } else {
1083    //        Double_t det = 0.;
1084    //        TMatrixD invmatrix = (TMatrixD)(mymat17->Invert(&det));
1085    //        printf(" Mymat3 determinant is %f \n",det);
1086    //        cf->WriteInvertedFullMatrix((TMatrixD)invmatrix,1117);
1087    //      };
1088    //      cf->WriteFullMatrix(mymat17, 117);
1089    
1090    
1091            //
1092            //      cf->UnLoadFullMatrix(i);
1093            //      cf->UnLoadFullNMatrix(i);
1094            //      delete fmatrix;
1095            //      delete fnmat;
1096            //
1097        };        };
       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);  
1098      };      };
1099        //
1100        printf(" done, closing file and exiting\n");
1101        //
1102    };    };
1103        //
   printf(" done, closing file and exiting\n");  
   
1104    cf->CloseMatrixFile();    cf->CloseMatrixFile();
1105      //
1106    cf->Delete();    cf->Delete();
1107      //
1108  }  }
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.7
  ViewVC Help
Powered by ViewVC 1.1.23