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