/[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.3 by mocchiut, Thu Jan 3 10:02:26 2008 UTC revision 1.6 by mocchiut, Tue Jan 15 12:41:38 2008 UTC
# Line 25  Int_t nbin; Line 25  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];  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 33  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[17];  TMatrixD *fqplane;
39  TMatrixD *fnqplane[17];  TMatrixD *fnqplane;
40  TMatrixF *fmatrix[17];  //TMatrixD *fqplane[17];
41    //TMatrixD *fnqplane[17];
42    //TMatrixF *fmatrix[17];
43    //TMatrixF *fnmat[17];
44    //TMatrixD *fmatrix;
45    //TMatrixD *fnmat;
46    TMatrixF *fmatrix;
47    //TMatrixF *fnmat;
48  TMatrixF *fnmat[17];  TMatrixF *fnmat[17];
49    //Int_t finmat[43][191];
50    
51  //===============================================================================  //===============================================================================
52  bool Select( PamLevel2* event ){  bool Select( PamLevel2* event ){
# Line 241  bool Select( PamLevel2* event ){ Line 250  bool Select( PamLevel2* event ){
250      if ( view == 1 && (plane >0 || strip > ti0) ) break;      if ( view == 1 && (plane >0 || strip > ti0) ) break;
251    };    };
252    //  if ( event->GetCaloLevel1()->qtotpl(0) > 7. ) return false;    //  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 291  void CreateHistos( PamLevel2* event , TS Line 303  void CreateHistos( PamLevel2* event , TS
303    //    //
304    memset(rmean, 0, 17*sizeof(Float_t));    memset(rmean, 0, 17*sizeof(Float_t));
305    memset(ntot, 0, 17*sizeof(Int_t));    memset(ntot, 0, 17*sizeof(Int_t));
306      //  memset(finmat, 0, 43*191*sizeof(Int_t));  
307    //    //
308    for (Int_t i=0; i < 17 ; i++){  //  for (Int_t i=0; i < 17 ; i++){
309        for (Int_t i=3; i < 4 ; i++){
310      if ( !FULL ){      if ( !FULL ){
311        matrix[i] = new TMatrixD(43,43);        matrix[i] = new TMatrixD(43,43);
312        qplane[i] = new TArrayF(43);        qplane[i] = new TArrayF(43);
# Line 300  void CreateHistos( PamLevel2* event , TS Line 314  void CreateHistos( PamLevel2* event , TS
314        nmat[i] = new TMatrixD(43,43);        nmat[i] = new TMatrixD(43,43);
315      } else {      } else {
316        if ( MATRIX ){        if ( MATRIX ){
317          fmatrix[i] = new TMatrixF(4128,4128);          //      fmatrix = new TMatrixF(4128,4128);
318          fnmat[i] = new TMatrixF(4128,4128);          //      fnmat = new TMatrixF(4128,4128);
319            //      fmatrix = new TMatrixF(8213,8213);
320            //      fnmat = new TMatrixF(8213,8213);
321            fmatrix = new TMatrixF(MDIM,MDIM);
322            //      fnmat = new TMatrixF(MDIM,MDIM);
323            fnmat[i] = new TMatrixF(43,191);
324    //      cf->WriteFullMatrix(fmatrix, i);
325            //      cf->WriteFullNMatrix(fnmat, i);
326    //      delete fmatrix;
327            //      delete fnmat;
328          //fnmat[i] = new TMatrixI(8213,8213);          //fnmat[i] = new TMatrixI(8213,8213);
329        } else {        } else {
330          fqplane[i] = new TMatrixD(43,191); // 43 planes x 191 strip (= 1 + 95 x 2, one strip is the one transversed by the track that could be on the extreme right or left)          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)
331          fnqplane[i] = new TMatrixD(43,191);//          fnqplane = new TMatrixD(43,191);//
332            //
333            cf->WriteFullMean(fqplane, i);
334            cf->WriteFullNMean(fnqplane, i);
335            delete fqplane;
336            delete fnqplane;
337            //
338        };        };
339      };      };
340    };    };
# Line 365  void FindAverage( PamLevel2* L2, int iev Line 394  void FindAverage( PamLevel2* L2, int iev
394      //      //
395      // FULL CALORIMETER      // FULL CALORIMETER
396      //      //
397        fqplane = cf->LoadFullAverage(rbi);
398        fnqplane = cf->LoadFullNAverage(rbi);
399      CaloTrkVar *ct = L2->GetTrack(0)->GetCaloTrack();      CaloTrkVar *ct = L2->GetTrack(0)->GetCaloTrack();
400      //      //
401      Int_t nplane = 0;      Int_t nplane = 0;
# Line 388  void FindAverage( PamLevel2* L2, int iev Line 419  void FindAverage( PamLevel2* L2, int iev
419          //          //
420          for (Int_t k=0; k<191; k++){          for (Int_t k=0; k<191; k++){
421            mstrip = cd + k;            mstrip = cd + k;
422            if ( mstrip < (191-cs) ) (*fnqplane[rbi])[nplane][mstrip] += 1.;            //      if ( mstrip < (191-cs) ) (*fnqplane[rbi])[nplane][mstrip] += 1.;
423              if ( mstrip < (191-cs) ) (*fnqplane)[nplane][mstrip] += 1.;
424          };          };
425        };        };
426      };      };
# Line 407  void FindAverage( PamLevel2* L2, int iev Line 439  void FindAverage( PamLevel2* L2, int iev
439        //        //
440        mstrip = cd + strip;        mstrip = cd + strip;
441        //        //
442        (*fqplane[rbi])[nplane][mstrip] += mip;        //      (*fqplane[rbi])[nplane][mstrip] += mip;
443          (*fqplane)[nplane][mstrip] += mip;
444        //        //
445      };      };
446            //
447        cf->WriteFullMean(fqplane, rbi);
448        cf->WriteFullNMean(fnqplane, rbi);
449        cf->UnLoadFullAverage(rbi);
450        cf->UnLoadFullNAverage(rbi);
451        delete fqplane;
452        delete fnqplane;
453      //      //
454    };    };
455  }  }
# Line 435  void CalculateAverage(){ Line 474  void CalculateAverage(){
474        //        //
475      };      };
476    } else {    } else {
477        //
478      for (Int_t i=0; i<nbin-1; i++){      for (Int_t i=0; i<nbin-1; i++){
479          fqplane = cf->LoadFullAverage(i);
480          fnqplane = cf->LoadFullNAverage(i);
481        if ( ntot[i] > 0 ) rmean[i] /= (Float_t)(ntot[i]);        if ( ntot[i] > 0 ) rmean[i] /= (Float_t)(ntot[i]);
482        //        //
483        for (Int_t j=0; j<43 ; j++){        for (Int_t j=0; j<43 ; j++){
484          for (Int_t k=0; k<191; k++){          for (Int_t k=0; k<191; k++){
485            if ( (*fnqplane[i])[j][k] > 0 ){    //        if ( (*fnqplane[i])[j][k] > 0 ){  
486              (*fqplane[i])[j][k] /= (Float_t)(*fnqplane[i])[j][k];  //          (*fqplane[i])[j][k] /= (Float_t)(*fnqplane[i])[j][k];
487    //        } else {
488    //          (*fqplane[i])[j][k] = 0.;
489    //        };
490    //        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]);
491              if ( (*fnqplane)[j][k] > 0 ){  
492                if ( (*fqplane)[j][k] == 0. ) (*fqplane)[j][k] = 0.7;
493                (*fqplane)[j][k] /= (Float_t)(*fnqplane)[j][k];
494            } else {            } else {
495              (*fqplane[i])[j][k] = 0.;              (*fqplane)[j][k] = 0.;
496            };            };
497            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]);            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]);
498          };          };
499        };        };
500          cf->WriteFullMean(fqplane, i);
501          cf->WriteFullNMean(fnqplane, i);
502          cf->UnLoadFullAverage(i);
503          cf->UnLoadFullNAverage(i);
504          delete fqplane;
505          delete fnqplane;
506      };      };
507      for (Int_t i=0; i<nbin-1; i++){      //
508        //  //     for (Int_t i=0; i<nbin-1; i++){
509        cf->WriteFullMean(fqplane[i], i);  //       //
510        //  //       cf->WriteFullMean(fqplane[i], i);
511      };  //       //
512    //     };
513    };    };
514    //      //  
515    cf->WriteNumBin(nbin);    cf->WriteNumBin(nbin);
# Line 485  void FindMatrix( PamLevel2* L2, int iev Line 540  void FindMatrix( PamLevel2* L2, int iev
540      };      };
541    };    };
542    //    //
543    if ( rbi != 3 ) return;
544    if ( erig < rig[0] ) return;    if ( erig < rig[0] ) return;
545    if ( erig >= rig[nbin-1] ) return;    if ( erig >= rig[nbin-1] ) return;
546    //    //
# Line 527  void FindMatrix( PamLevel2* L2, int iev Line 583  void FindMatrix( PamLevel2* L2, int iev
583      //      //
584      // FULL CALORIMETER      // FULL CALORIMETER
585      //      //
586        printf(" Retrieve matrix %i IEV %i \n",rbi,iev);
587        //    fmatrix = cf->LoadFullMatrix(rbi);
588    //    cf->LoadFullMatrix(rbi,fmatrix);
589        //    fnmat = cf->LoadFullNMatrix(rbi);
590        printf(" done \n");
591        printf(" start loop \n");
592        //
593      CaloTrkVar *ct = L2->GetTrack(0)->GetCaloTrack();      CaloTrkVar *ct = L2->GetTrack(0)->GetCaloTrack();
594      //      //
595      Int_t nplane = 0;      Int_t nplane = 0;
# Line 535  void FindMatrix( PamLevel2* L2, int iev Line 598  void FindMatrix( PamLevel2* L2, int iev
598      Int_t strip = 0;      Int_t strip = 0;
599      Float_t mip = 0.;      Float_t mip = 0.;
600      //      //
601      Int_t mindgf = 48;      //    Int_t mindgf = 48;
602      Int_t dgf = 143;      //    Int_t dgf = 143;
603    //       Int_t mindgf = 0; //tutto
604    //        Int_t dgf = 191; //tutto
605    //    Int_t mindgf = 94;
606    //    Int_t dgf = 96;
607        Int_t mindgf = 84;
608        Int_t dgf = 106;
609      Int_t cs = 0;      Int_t cs = 0;
610      Int_t cd = 0;      Int_t cd = 0;
611      Int_t mstrip = 0;      Int_t mstrip = 0;
# Line 561  void FindMatrix( PamLevel2* L2, int iev Line 630  void FindMatrix( PamLevel2* L2, int iev
630        //        //
631      };      };
632      //      //
633      Float_t mip1;      Float_t mip1 = 1.;
634      Int_t cs1;      Int_t cs1;
635      Int_t cd1;      Int_t cd1;
636      Float_t mip2;      Float_t mip2 = 1.;
637      Int_t cs2;      Int_t cs2;
638      Int_t cd2;      Int_t cd2;
639      Int_t mi = -1;      Int_t mi = -1;
# Line 580  void FindMatrix( PamLevel2* L2, int iev Line 649  void FindMatrix( PamLevel2* L2, int iev
649      Int_t mstrip2min = 0;      Int_t mstrip2min = 0;
650      Int_t mstrip2max = 0;      Int_t mstrip2max = 0;
651      //      //
652        Int_t toto = 0;
653        //
654      for (Int_t nplane1=0; nplane1<43; nplane1++){      for (Int_t nplane1=0; nplane1<43; nplane1++){
655        for (Int_t mstrip1=0; mstrip1<191; mstrip1++){        if ( nplane1 >= 37 ) nn1 = nplane1 + 1;
656          vi1 = 1;
657          if ( nn1%2 ) vi1 = 0;
658          pl1 = (nn1 - 1 + vi1)/2;
659          //
660          cs1 = ct->tibar[pl1][vi1] - 1; // convertire nplane in pl1 e vi1
661          //
662          cd1 = 95 - cs1;
663          //
664          mstrip1min = TMath::Max(mindgf,(cd1+0));
665          mstrip1max = TMath::Min(dgf,(cd1+95)) + 1;
666          //
667          if ( nplane1 == 0 || nplane1 == 42 ) printf(" pl %i mstrip1min %i mstrip1max %i \n",nplane1,mstrip1min,mstrip1max);
668          //
669          for (Int_t mstrip1=mstrip1min; mstrip1<mstrip1max; mstrip1++){
670            //      printf(".\n");
671            //
672          mj = -1;          mj = -1;
673          //          //
674          if ( nplane1 >= 37 ) nn1 = nplane1 + 1;          mip1 = mipv[nplane1][mstrip1] - cf->GetFullAverageAt(nplane1,mstrip1,erig,rbi);
675          vi1 = 1;          //
676          if ( nn1%2 ) vi1 = 0;          mi = (nplane1 * 191) + mstrip1;
677          pl1 = (nn1 - 1 + vi1)/2;          //
678          //          //      if ( mstrip1 > mstrip1min ) break;
679          cs1 = ct->tibar[pl1][vi1] - 1; // convertire nplane in pl1 e vi1          //      if ( mstrip1 > dgf ) break;
680          //          //      if ( mstrip1 >= mindgf && mstrip1 <= dgf && mstrip1 >= mstrip1min && mstrip1 <= mstrip1max ){
681          cd1 = 95 - cs1;          //
682          //          //      finmat[nplane1][mstrip1]++;
683          mstrip1min = cd1 + 0;          (*fnmat[rbi])[nplane1][mstrip1] += 1.;
684          mstrip1max = cd1 + 95;          //
685          mip1 = mipv[nplane1][mstrip1];          if ( mip1 != 0. ){
686          //            //
         if ( mstrip1 >= mindgf && mstrip1 <= dgf ){  
           mi++;  
           //  
687            for (Int_t nplane2=0; nplane2<43; nplane2++){            for (Int_t nplane2=0; nplane2<43; nplane2++){
688              for (Int_t mstrip2=0; mstrip2<191; mstrip2++){              //
689                //              if ( nplane2 >= 37 ) nn2 = nplane2 + 1;
690                if ( nplane2 >= 37 ) nn2 = nplane2 + 1;              vi2 = 1;
691                vi2 = 1;              if ( nn2%2 ) vi2 = 0;
692                if ( nn2%2 ) vi2 = 0;              pl1 = (nn2 - 1 + vi2)/2;
693                pl1 = (nn2 - 1 + vi2)/2;              //
694                //              cs2 = ct->tibar[pl2][vi2] - 1;
695                cs2 = ct->tibar[pl2][vi2] - 1;              //
696                cd2 = 95 - cs2;
697                //
698                //    mstrip2min = cd2 + 0;
699                //    mstrip2max = cd2 + 95;
700                mstrip2min = TMath::Max(mindgf,(cd2+0));
701                mstrip2max = TMath::Min(dgf,(cd2+95)) + 1;
702                //
703                if ( nplane1 == 0 && nplane2 == 0 && mstrip1==mstrip1min ) printf(" mstrip2min %i mstrip2max %i \n",mstrip2min,mstrip2max);
704                //
705                for (Int_t mstrip2=mstrip2min; mstrip2<mstrip2max; mstrip2++){
706                //                //
707                cd2 = 95 - cs2;                mip2 = mipv[nplane2][mstrip2] - cf->GetFullAverageAt(nplane2,mstrip2,erig,rbi);
708                //                //
709                mstrip2min = cd2 + 0;                if ( mip2 != 0. ){
710                mstrip2max = cd2 + 95;                  //
711                //                  mj = (nplane2 * 191) + mstrip2;
712                mip2 = mipv[nplane2][mstrip2];                  //            mj++;
713                //                  //
714                if ( mstrip2 >= mindgf && mstrip2 <= dgf ){                  //            if ( mstrip2 > mstrip2min ) break;
715                  mj++;                  //            if ( mstrip2 > dgf ) break;
716                  (*fmatrix[rbi])[mi][mj] += (mipv[nplane1][mstrip1] - cf->GetFullAverageAt(nplane1,mstrip1,erig)) * (mipv[nplane2][mstrip2] - cf->GetFullAverageAt(nplane2,mstrip2,erig));                  //            if ( mstrip2 >= mindgf && mstrip2 <= dgf && mstrip2 >= mstrip2min && mstrip2 <= mstrip2max ){
717                  if ( mstrip1 >= mstrip1min && mstrip1 <= mstrip1max && mstrip2 >= mstrip2min && mstrip2 <= mstrip2max){                  //              if ( mstrip1 >= mstrip1min && mstrip1 <= mstrip1max && mstrip2 >= mstrip2min && mstrip2 <= mstrip2max){
718                    (*fnmat[rbi])[mi][mj] += 1.;                  //              (*fmatrix[rbi])[mi][mj] += (mipv[nplane1][mstrip1] - cf->GetFullAverageAt(nplane1,mstrip1,erig)) * (mipv[nplane2][mstrip2] - cf->GetFullAverageAt(nplane2,mstrip2,erig));
719                  };                  //                (*fnmat[rbi])[mi][mj] += 1.;
720                    (*fmatrix)[mi][mj] += (mip1 * mip2); // giusto
721    //              (*fmatrix)[mi][mj] += (mip1 * mip2) * 1000000.;
722                    toto++;
723                    //            (*fmatrix)[mi][mj] += 1.;
724                    //              cf->GetFullAverageAt(nplane1,mstrip1,erig,rbi);
725                    //              cf->GetFullAverageAt(nplane2,mstrip2,erig,rbi);
726                    //                (*fnmat)[mi][mj] += 1.;
727                    //              };
728                };                };
729              };              };
730            };            };
# Line 630  void FindMatrix( PamLevel2* L2, int iev Line 732  void FindMatrix( PamLevel2* L2, int iev
732        };        };
733      };      };
734      //        //  
735        printf(" toto = %i \n",toto);
736        printf("\n done \n");
737        printf(" write matrix \n");
738    //    cf->WriteFullMatrix(fmatrix, rbi);
739        //    cf->WriteFullNMatrix(fnmat, rbi);
740        printf(" done \n");
741        printf(" unload matrix \n");
742     //   cf->UnLoadFullMatrix(rbi);
743        //    cf->UnLoadFullNMatrix(rbi);
744        printf(" done \n");
745        printf(" delete matrix \n");
746    //    delete fmatrix;
747        //    delete fnmat;
748        printf(" done \n");
749    };    };
750  }  }
751    
# Line 685  void SaveHistos(){ Line 801  void SaveHistos(){
801        //        //
802        // FULL        // FULL
803        //        //
804        for (Int_t i=0; i<nbin-1; i++){        //      for (Int_t i=0; i<nbin-1; i++){
805    //      for (Int_t i=3; i<5; i++){
806          for (Int_t i=3; i<4; i++){
807          //          //
808          // determine the average matrix          // determine the average matrix
809          //              //    
810          for (Int_t ii=0; ii<4171; ii++){  //      fmatrix = cf->LoadFullMatrix(i);
811            for (Int_t j=0; j<4171; j++){          //      fnmat = cf->LoadFullNMatrix(i);
812              if ( (*fnmat[i])[ii][j] > 0. ){          //
813                (*fmatrix[i])[ii][j] /= (*fnmat[i])[ii][j];  //      for (Int_t ii=0; ii<MDIM; ii++){
814              } else {  //        for (Int_t j=0; j<MDIM; j++){
815                (*fmatrix[i])[ii][j] = 0.;  // //       if ( (*fnmat[i])[ii][j] > 0. ){
816    // //         (*fmatrix[i])[ii][j] /= (*fnmat[i])[ii][j];
817    // //       } else {
818    // //         (*fmatrix[i])[ii][j] = 0.;
819    // //       };
820    //          if ( (*fnmat)[ii][j] > 0. ){
821    //            (*fmatrix)[ii][j] /= (*fnmat)[ii][j];
822    //          } else {
823    //            (*fmatrix)[ii][j] = 0.;
824    //          };
825    //        };
826    //      };
827    //
828    //      TMatrixF *mymat = new TMatrixF(129,129);
829            TMatrixF *mymat = new TMatrixF(989,989);
830            Int_t i1 = -1;
831            Int_t j1 = -1;
832            int mi,mj;
833            Int_t nonzero = 0;
834            Int_t nonzero1 = 0;
835            for (Int_t ii=0; ii<43; ii++){
836              for (Int_t j=0; j<191; j++){
837    //          if ( (*fnmat[i])[ii][j] > 0. ){
838    //            (*fmatrix[i])[ii][j] /= (*fnmat[i])[ii][j];
839    //          } else {
840    //            (*fmatrix[i])[ii][j] = 0.;
841    //          };
842                i1 =  (ii * 191) + j;
843                //      j1 = -1;
844                for (Int_t iij=0; iij<43; iij++){
845                  for (Int_t jj=0; jj<191; jj++){
846                    //
847                    j1 = (iij * 191) + jj;
848                    //              j1++;
849                    //              if ( finmat[ii][j] > 0 ){
850                    //                (*fmatrix)[i1][j1] /= finmat[ii][j];
851                    if ( (*fnmat[i])[ii][j] == 0. || (*fmatrix)[i1][j1] == 0. || !((*fmatrix)[i1][j1] == (*fmatrix)[i1][j1]) ){
852                      (*fmatrix)[i1][j1] = 1.;
853                    } else {
854                      (*fmatrix)[i1][j1] /= (*fnmat[i])[ii][j];
855                      nonzero++;
856                      if ( i1 == 0 ) nonzero1++;
857                    };
858                    //              if ( j>=94 && j <=96 && jj >=94 && jj<=96 ){
859                    //                      mi = (ii*3) + j -94;
860                    //                      mj = (iij*3) + jj -94;
861                    //                      (*mymat)[mi][mj] = (*fmatrix)[i1][j1];
862                    //              };
863    
864    
865                    if ( j>=84 && j <=106 && jj >=84 && jj<=106 ){
866                            mi = (ii*3) + j -84;
867                            mj = (iij*3) + jj -84;
868                            (*mymat)[mi][mj] = (*fmatrix)[i1][j1];
869                    };
870    
871                  };
872              };              };
873            };            };
874          };          };
875          //          //
876          cf->WriteFullMatrix(fmatrix[i],i);          printf(" Matrix has %i non-zero elements \n",nonzero);  
877            printf(" Matrix has %i non-zero elements on the first row\n",nonzero1);
878            //
879    //      Bool_t BAD = false;
880    //      for (Int_t ii=0; ii<43; ii++){
881    //        for (Int_t j=0; j<191; j++){
882    //          //
883    //          i1 =  (ii * 191) + j;
884    //          //
885    //          for (Int_t iij=0; iij<43; iij++){
886    //            for (Int_t jj=0; jj<191; jj++){
887    //              //
888    //              j1 = (iij * 191) + jj;
889    //              //
890    //              //              printf(" ROW %i COLUMN %i VALUE %f \n",i1,j1,(*fmatrix)[i1][j1]);
891    //              if ( (*fmatrix)[i1][j1] == 0. || !((*fmatrix)[i1][j1]==(*fmatrix)[i1][j1]) ){
892    //                printf(" ROW %i COLUMN %i VALUE %f \n",i1,j1,(*fmatrix)[i1][j1]);
893    //                printf(" che schifo! \n");
894    //                BAD = true;
895    //              };
896    //              //
897    //            };
898    //          };
899    //        };
900    //      };
901    //      //
902    //      if ( BAD ) printf(" questa matrice fa cagare \n");
903            //
904            //
905            //      cf->WriteFullMatrix(fmatrix[i],i);
906            cf->WriteFullMatrix(fmatrix, i);
907            //      cf->WriteFullNMatrix(fnmat, i);
908            cf->WriteFullNMatrix(fnmat[i], i);
909          //          //
910          if ( fmatrix[i]->Determinant() == 0. ){          //      if ( fmatrix[i]->Determinant() == 0. ){
911            if ( fmatrix->Determinant() == 0. ){
912            printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i);            printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i);
913          } else {          } else {
914            Double_t det = 0.;            Double_t det = 0.;
915            TMatrixF invmatrix = (TMatrixF)(fmatrix[i]->Invert(&det));            //      TMatrixF invmatrix = (TMatrixF)(fmatrix[i]->Invert(&det));
916              //      printf(" Bin %i determinant is %f \n",i,det);
917              //      cf->WriteInvertedFullMatrix((TMatrixF)invmatrix,i);
918              TMatrixF invmatrix = (TMatrixF)(fmatrix->Invert(&det));
919            printf(" Bin %i determinant is %f \n",i,det);            printf(" Bin %i determinant is %f \n",i,det);
920            cf->WriteInvertedFullMatrix((TMatrixF)invmatrix,i);            cf->WriteInvertedFullMatrix((TMatrixF)invmatrix,i);
921          };          };
922    
923            if ( mymat->Determinant() == 0. ){
924              printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i);
925            } else {
926              Double_t det = 0.;
927              TMatrixF invmatrix = (TMatrixF)(mymat->Invert(&det));
928              printf(" Bin %i determinant is %f \n",i,det);
929              cf->WriteInvertedFullMatrix((TMatrixF)invmatrix,i);
930            };
931            cf->WriteFullMatrix(mymat, 99);
932    
933    
934            //
935            cf->UnLoadFullMatrix(i);
936            //      cf->UnLoadFullNMatrix(i);
937            delete fmatrix;
938            //      delete fnmat;
939            //
940        };        };
941      };      };
942      //      //
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.6
  ViewVC Help
Powered by ViewVC 1.1.23