/[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.6 by mocchiut, Tue Jan 15 12:41:38 2008 UTC revision 1.7 by mocchiut, Mon Jan 21 10:24:10 2008 UTC
# Line 35  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;  //TMatrixD *fqplane;
39  TMatrixD *fnqplane;  //TMatrixD *fnqplane;
40  //TMatrixD *fqplane[17];  TMatrixD *fqplane[17];
41  //TMatrixD *fnqplane[17];  TMatrixD *fnqplane[17];
42  //TMatrixF *fmatrix[17];  TMatrixD *fmatrix[17];
43  //TMatrixF *fnmat[17];  //TMatrixF *fnmat[17];
44  //TMatrixD *fmatrix;  //TMatrixD *fmatrix;
45  //TMatrixD *fnmat;  //TMatrixD *fnmat;
46  TMatrixF *fmatrix;  //TMatrixF *fmatrix;
47  //TMatrixF *fnmat;  //TMatrixF *fnmat;
48  TMatrixF *fnmat[17];  TMatrixF *fnmat[17];
49  //Int_t finmat[43][191];  //Int_t finmat[43][191];
# Line 251  bool Select( PamLevel2* event ){ Line 251  bool Select( PamLevel2* event ){
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;    //  if ( rigidity > 2.2 || rigidity < 1.5 ) return false;
255    //  printf(" rig %f CRIG %i SRIG %i \n",rigidity,CRIG,SRIG);    //  printf(" rig %f CRIG %i SRIG %i \n",rigidity,CRIG,SRIG);
256    //    //
257    return true;    return true;
# Line 304  void CreateHistos( PamLevel2* event , TS Line 304  void CreateHistos( PamLevel2* event , TS
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));      //  memset(finmat, 0, 43*191*sizeof(Int_t));  
307    //  Double_t tol = 1E-20;
308    //    //
309  //  for (Int_t i=0; i < 17 ; i++){    for (Int_t i=0; i < 17 ; i++){
310      for (Int_t i=3; i < 4 ; i++){      //  for (Int_t i=3; i < 4 ; i++){
311      if ( !FULL ){      if ( !FULL ){
312        matrix[i] = new TMatrixD(43,43);        matrix[i] = new TMatrixD(43,43);
313        qplane[i] = new TArrayF(43);        qplane[i] = new TArrayF(43);
# Line 318  void CreateHistos( PamLevel2* event , TS Line 319  void CreateHistos( PamLevel2* event , TS
319          //      fnmat = new TMatrixF(4128,4128);          //      fnmat = new TMatrixF(4128,4128);
320          //      fmatrix = new TMatrixF(8213,8213);          //      fmatrix = new TMatrixF(8213,8213);
321          //      fnmat = new TMatrixF(8213,8213);          //      fnmat = new TMatrixF(8213,8213);
322          fmatrix = new TMatrixF(MDIM,MDIM);          //      fmatrix = new TMatrixF(MDIM,MDIM);
323          //      fnmat = new TMatrixF(MDIM,MDIM);          //      fnmat = new TMatrixF(MDIM,MDIM);
324          fnmat[i] = new TMatrixF(43,191);          //      fmatrix[i] = new TMatrixF(1849,1849);
325  //      cf->WriteFullMatrix(fmatrix, i);          //      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);          //      cf->WriteFullNMatrix(fnmat, i);
331  //      delete fmatrix;          //      delete fmatrix;
332          //      delete fnmat;          //      delete fnmat;
333          //fnmat[i] = new TMatrixI(8213,8213);          //fnmat[i] = new TMatrixI(8213,8213);
334        } else {        } 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)          //      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);//          //      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          cf->WriteFullMean(fqplane, i);          //      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          cf->WriteFullNMean(fnqplane, i);          //
340          delete fqplane;          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          delete fnqplane;          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      };      };
# Line 394  void FindAverage( PamLevel2* L2, int iev Line 404  void FindAverage( PamLevel2* L2, int iev
404      //      //
405      // FULL CALORIMETER      // FULL CALORIMETER
406      //      //
407      fqplane = cf->LoadFullAverage(rbi);      //    fqplane = cf->LoadFullAverage(rbi);
408      fnqplane = cf->LoadFullNAverage(rbi);      //    fnqplane = cf->LoadFullNAverage(rbi);
409      CaloTrkVar *ct = L2->GetTrack(0)->GetCaloTrack();      CaloTrkVar *ct = L2->GetTrack(0)->GetCaloTrack();
410      //      //
411      Int_t nplane = 0;      Int_t nplane = 0;
# Line 409  void FindAverage( PamLevel2* L2, int iev Line 419  void FindAverage( PamLevel2* L2, int iev
419      Int_t mstrip = 0;      Int_t mstrip = 0;
420      //      //
421      for (Int_t j=0; j<2; j++){      for (Int_t j=0; j<2; j++){
422        for (Int_t i=0; i<21; i++){        for (Int_t i=0; i<22; i++){
423          nplane = 1 - j + 2*i;          nplane = 1 - j + 2*i;
424          if ( nplane > 37 ) nplane--;          if ( nplane > 37 ) nplane--;
425          //          //
# Line 417  void FindAverage( PamLevel2* L2, int iev Line 427  void FindAverage( PamLevel2* L2, int iev
427          //          //
428          cd = 95 - cs;            cd = 95 - cs;  
429          //          //
430            Int_t oldstr = -1;
431          for (Int_t k=0; k<191; k++){          for (Int_t k=0; k<191; k++){
432            mstrip = cd + k;            mstrip = cd + k;
433            //      if ( mstrip < (191-cs) ) (*fnqplane[rbi])[nplane][mstrip] += 1.;            //      if ( mstrip < (191-cs) ) (*fnqplane[rbi])[nplane][mstrip] += 1.;
434            if ( mstrip < (191-cs) ) (*fnqplane)[nplane][mstrip] += 1.;            //      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      };      };
# Line 439  void FindAverage( PamLevel2* L2, int iev Line 455  void FindAverage( PamLevel2* L2, int iev
455        //        //
456        mstrip = cd + strip;        mstrip = cd + strip;
457        //        //
458          Int_t lstr = cf->ConvertStrip(mstrip);
459        //      (*fqplane[rbi])[nplane][mstrip] += mip;        //      (*fqplane[rbi])[nplane][mstrip] += mip;
460        (*fqplane)[nplane][mstrip] += mip;        //      (*fqplane)[nplane][mstrip] += mip;
461          (*fqplane[rbi])[nplane][lstr] += mip;
462        //        //
463      };      };
464      //      //
465      cf->WriteFullMean(fqplane, rbi);      //    cf->WriteFullMean(fqplane, rbi);
466      cf->WriteFullNMean(fnqplane, rbi);      //    cf->WriteFullNMean(fnqplane, rbi);
467      cf->UnLoadFullAverage(rbi);      //    cf->UnLoadFullAverage(rbi);
468      cf->UnLoadFullNAverage(rbi);      //    cf->UnLoadFullNAverage(rbi);
469      delete fqplane;      //    delete fqplane;
470      delete fnqplane;      //    delete fnqplane;
471      //      //
472    };    };
473  }  }
# Line 476  void CalculateAverage(){ Line 494  void CalculateAverage(){
494    } else {    } else {
495      //      //
496      for (Int_t i=0; i<nbin-1; i++){      for (Int_t i=0; i<nbin-1; i++){
497        fqplane = cf->LoadFullAverage(i);        //      fqplane = cf->LoadFullAverage(i);
498        fnqplane = cf->LoadFullNAverage(i);        //      fnqplane = cf->LoadFullNAverage(i);
499        if ( ntot[i] > 0 ) rmean[i] /= (Float_t)(ntot[i]);        if ( ntot[i] > 0 ) rmean[i] /= (Float_t)(ntot[i]);
500        //        //
501        for (Int_t j=0; j<43 ; j++){        for (Int_t j=0; j<43 ; j++){
502          for (Int_t k=0; k<191; k++){          //      for (Int_t k=0; k<191; k++){
503  //        if ( (*fnqplane[i])[j][k] > 0 ){            //      for (Int_t k=0; k<43; k++){
504  //          (*fqplane[i])[j][k] /= (Float_t)(*fnqplane[i])[j][k];          for (Int_t k=0; k<31; k++){
505  //        } else {            //      if ( (*fnqplane[i])[j][k] > 0 ){  
506  //          (*fqplane[i])[j][k] = 0.;            //        (*fqplane[i])[j][k] /= (Float_t)(*fnqplane[i])[j][k];
507  //        };            //      } else {
508  //        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]);            //        (*fqplane[i])[j][k] = 0.;
509            if ( (*fnqplane)[j][k] > 0 ){              //      };
510              if ( (*fqplane)[j][k] == 0. ) (*fqplane)[j][k] = 0.7;            //      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              (*fqplane)[j][k] /= (Float_t)(*fnqplane)[j][k];            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 {            } else {
515              (*fqplane)[j][k] = 0.;              (*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]);            //      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);        cf->WriteFullMean(fqplane[i], i);
521        cf->WriteFullNMean(fnqplane, i);        cf->WriteFullNMean(fnqplane[i], i);
522        cf->UnLoadFullAverage(i);        //      cf->UnLoadFullAverage(i);
523        cf->UnLoadFullNAverage(i);        //      cf->UnLoadFullNAverage(i);
524        delete fqplane;        //      delete fqplane;
525        delete fnqplane;        //      delete fnqplane;
526      };      };
527      //      //
528  //     for (Int_t i=0; i<nbin-1; i++){      //     for (Int_t i=0; i<nbin-1; i++){
529  //       //      //       //
530  //       cf->WriteFullMean(fqplane[i], i);      //       cf->WriteFullMean(fqplane[i], i);
531  //       //      //       //
532  //     };      //     };
533    };    };
534    //      //  
535    cf->WriteNumBin(nbin);    cf->WriteNumBin(nbin);
# Line 540  void FindMatrix( PamLevel2* L2, int iev Line 560  void FindMatrix( PamLevel2* L2, int iev
560      };      };
561    };    };
562    //    //
 if ( rbi != 3 ) return;  
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    //    //
# Line 583  if ( rbi != 3 ) return; Line 602  if ( rbi != 3 ) return;
602      //      //
603      // FULL CALORIMETER      // FULL CALORIMETER
604      //      //
605      printf(" Retrieve matrix %i IEV %i \n",rbi,iev);      //    if ( rbi != 3 ) return;
606        printf(" matrix %i IEV %i \n",rbi,iev);
607      //    fmatrix = cf->LoadFullMatrix(rbi);      //    fmatrix = cf->LoadFullMatrix(rbi);
608  //    cf->LoadFullMatrix(rbi,fmatrix);      //    cf->LoadFullMatrix(rbi,fmatrix);
609      //    fnmat = cf->LoadFullNMatrix(rbi);      //    fnmat = cf->LoadFullNMatrix(rbi);
610      printf(" done \n");      //    printf(" done \n");
611      printf(" start loop \n");      //    printf(" start loop \n");
612      //      //
613      CaloTrkVar *ct = L2->GetTrack(0)->GetCaloTrack();      CaloTrkVar *ct = L2->GetTrack(0)->GetCaloTrack();
614      //      //
# Line 600  if ( rbi != 3 ) return; Line 620  if ( rbi != 3 ) return;
620      //      //
621      //    Int_t mindgf = 48;      //    Int_t mindgf = 48;
622      //    Int_t dgf = 143;      //    Int_t dgf = 143;
623  //       Int_t mindgf = 0; //tutto      //       Int_t mindgf = 0; //tutto
624  //        Int_t dgf = 191; //tutto      //        Int_t dgf = 191; //tutto
625  //    Int_t mindgf = 94;      //    Int_t mindgf = 94;
626  //    Int_t dgf = 96;      //    Int_t dgf = 96;
627      Int_t mindgf = 84;      //    Int_t mindgf = 84;
628      Int_t dgf = 106;      //    Int_t dgf = 106;
629        Int_t mindgf = 0;
630        Int_t dgf = 43;
631      Int_t cs = 0;      Int_t cs = 0;
632      Int_t cd = 0;      Int_t cd = 0;
633      Int_t mstrip = 0;      Int_t mstrip = 0;
634      //      //
635      Float_t mipv[43][191];      //    Float_t mipv[43][43];
636      memset(mipv,0,43*191*sizeof(Float_t));      //    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++ ){      for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
641        //        //
# Line 626  if ( rbi != 3 ) return; Line 650  if ( rbi != 3 ) return;
650        //        //
651        mstrip = cd + strip;        mstrip = cd + strip;
652        //        //
653        mipv[nplane][mstrip] = mip;        Int_t lstr = cf->ConvertStrip(mstrip);
654          mipv[nplane][lstr] += mip;
655        //        //
656      };      };
657      //      //
# Line 661  if ( rbi != 3 ) return; Line 686  if ( rbi != 3 ) return;
686        //        //
687        cd1 = 95 - cs1;        cd1 = 95 - cs1;
688        //        //
689        mstrip1min = TMath::Max(mindgf,(cd1+0));        Int_t at1 = TMath::Max(0,(cd1+0));
690        mstrip1max = TMath::Min(dgf,(cd1+95)) + 1;        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 \n",nplane1,mstrip1min,mstrip1max);        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++){        for (Int_t mstrip1=mstrip1min; mstrip1<mstrip1max; mstrip1++){
699          //      printf(".\n");          //      printf(".\n");
# Line 673  if ( rbi != 3 ) return; Line 702  if ( rbi != 3 ) return;
702          //          //
703          mip1 = mipv[nplane1][mstrip1] - cf->GetFullAverageAt(nplane1,mstrip1,erig,rbi);          mip1 = mipv[nplane1][mstrip1] - cf->GetFullAverageAt(nplane1,mstrip1,erig,rbi);
704          //          //
705          mi = (nplane1 * 191) + mstrip1;          //      mi = (nplane1 * 191) + mstrip1;
706            //      mi = (nplane1 * 43) + mstrip1;
707            mi = (nplane1 * 31) + mstrip1;
708          //          //
709          //      if ( mstrip1 > mstrip1min ) break;          //      if ( mstrip1 > mstrip1min ) break;
710          //      if ( mstrip1 > dgf ) break;          //      if ( mstrip1 > dgf ) break;
# Line 697  if ( rbi != 3 ) return; Line 728  if ( rbi != 3 ) return;
728              //              //
729              //    mstrip2min = cd2 + 0;              //    mstrip2min = cd2 + 0;
730              //    mstrip2max = cd2 + 95;              //    mstrip2max = cd2 + 95;
731              mstrip2min = TMath::Max(mindgf,(cd2+0));              Int_t t1 = TMath::Max(0,(cd2+0));
732              mstrip2max = TMath::Min(dgf,(cd2+95)) + 1;              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);              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++){              for (Int_t mstrip2=mstrip2min; mstrip2<mstrip2max; mstrip2++){
739                //                //              
740                mip2 = mipv[nplane2][mstrip2] - cf->GetFullAverageAt(nplane2,mstrip2,erig,rbi);                mip2 = mipv[nplane2][mstrip2] - cf->GetFullAverageAt(nplane2,mstrip2,erig,rbi);
741                //                //
742                if ( mip2 != 0. ){                if ( mip2 != 0. ){
743                  //                  //
744                  mj = (nplane2 * 191) + mstrip2;                  //              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++;                  //            mj++;
757                  //                  //
758                  //            if ( mstrip2 > mstrip2min ) break;                  //            if ( mstrip2 > mstrip2min ) break;
# Line 717  if ( rbi != 3 ) return; Line 761  if ( rbi != 3 ) return;
761                  //              if ( mstrip1 >= mstrip1min && mstrip1 <= mstrip1max && mstrip2 >= mstrip2min && mstrip2 <= mstrip2max){                  //              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));                  //              (*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.;                  //                (*fnmat[rbi])[mi][mj] += 1.;
764                  (*fmatrix)[mi][mj] += (mip1 * mip2); // giusto                  (*fmatrix[rbi])[mi][mj] += (mip1 * mip2); // giusto
765  //              (*fmatrix)[mi][mj] += (mip1 * mip2) * 1000000.;                  //              (*fmatrix)[mi][mj] += (mip1 * mip2) * 1000000.;
766                  toto++;                  toto++;
767                  //            (*fmatrix)[mi][mj] += 1.;                  //            (*fmatrix)[mi][mj] += 1.;
768                  //              cf->GetFullAverageAt(nplane1,mstrip1,erig,rbi);                  //              cf->GetFullAverageAt(nplane1,mstrip1,erig,rbi);
# Line 734  if ( rbi != 3 ) return; Line 778  if ( rbi != 3 ) return;
778      //        //  
779      printf(" toto = %i \n",toto);      printf(" toto = %i \n",toto);
780      printf("\n done \n");      printf("\n done \n");
781      printf(" write matrix \n");      //    printf(" write matrix \n");
782  //    cf->WriteFullMatrix(fmatrix, rbi);      //    cf->WriteFullMatrix(fmatrix, rbi);
783      //    cf->WriteFullNMatrix(fnmat, rbi);      //    cf->WriteFullNMatrix(fnmat, rbi);
784      printf(" done \n");      //    printf(" done \n");
785      printf(" unload matrix \n");      //    printf(" unload matrix \n");
786   //   cf->UnLoadFullMatrix(rbi);      //   cf->UnLoadFullMatrix(rbi);
787      //    cf->UnLoadFullNMatrix(rbi);      //    cf->UnLoadFullNMatrix(rbi);
788      printf(" done \n");      //    printf(" done \n");
789      printf(" delete matrix \n");      //    printf(" delete matrix \n");
790  //    delete fmatrix;      //    delete fmatrix;
791      //    delete fnmat;      //    delete fnmat;
792      printf(" done \n");      //    printf(" done \n");
793    };    };
794  }  }
795    
# Line 801  void SaveHistos(){ Line 845  void SaveHistos(){
845        //        //
846        // FULL        // FULL
847        //        //
848        //      for (Int_t i=0; i<nbin-1; i++){        for (Int_t i=0; i<nbin-1; i++){
849  //      for (Int_t i=3; i<5; i++){          //      for (Int_t i=3; i<5; i++){
850        for (Int_t i=3; i<4; i++){          //      for (Int_t i=3; i<4; i++){
851          //          //
852          // determine the average matrix          // determine the average matrix
853          //              //    
854  //      fmatrix = cf->LoadFullMatrix(i);          //      fmatrix = cf->LoadFullMatrix(i);
855          //      fnmat = cf->LoadFullNMatrix(i);          //      fnmat = cf->LoadFullNMatrix(i);
856          //          //
857  //      for (Int_t ii=0; ii<MDIM; ii++){          //      for (Int_t ii=0; ii<MDIM; ii++){
858  //        for (Int_t j=0; j<MDIM; j++){          //        for (Int_t j=0; j<MDIM; j++){
859  // //       if ( (*fnmat[i])[ii][j] > 0. ){          // //       if ( (*fnmat[i])[ii][j] > 0. ){
860  // //         (*fmatrix[i])[ii][j] /= (*fnmat[i])[ii][j];          // //         (*fmatrix[i])[ii][j] /= (*fnmat[i])[ii][j];
861  // //       } else {          // //       } else {
862  // //         (*fmatrix[i])[ii][j] = 0.;          // //         (*fmatrix[i])[ii][j] = 0.;
863  // //       };          // //       };
864  //          if ( (*fnmat)[ii][j] > 0. ){          //          if ( (*fnmat)[ii][j] > 0. ){
865  //            (*fmatrix)[ii][j] /= (*fnmat)[ii][j];          //            (*fmatrix)[ii][j] /= (*fnmat)[ii][j];
866  //          } else {          //          } else {
867  //            (*fmatrix)[ii][j] = 0.;          //            (*fmatrix)[ii][j] = 0.;
868  //          };          //          };
869  //        };          //        };
870  //      };          //      };
871  //          //
872  //      TMatrixF *mymat = new TMatrixF(129,129);  //      TMatrixD *mymat3 = new TMatrixD(129,129);
873          TMatrixF *mymat = new TMatrixF(989,989);  //      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;          Int_t i1 = -1;
881          Int_t j1 = -1;          Int_t j1 = -1;
882          int mi,mj;          //      int mi,mj;
883          Int_t nonzero = 0;          Int_t nonzero = 0;
884          Int_t nonzero1 = 0;          Int_t nonzero1 = 0;
885          for (Int_t ii=0; ii<43; ii++){          for (Int_t ii=0; ii<43; ii++){
886            for (Int_t j=0; j<191; j++){            //      for (Int_t j=0; j<191; j++){
887  //          if ( (*fnmat[i])[ii][j] > 0. ){            //      for (Int_t j=0; j<43; j++){
888  //            (*fmatrix[i])[ii][j] /= (*fnmat[i])[ii][j];            for (Int_t j=0; j<31; j++){
889  //          } else {              //      if ( (*fnmat[i])[ii][j] > 0. ){
890  //            (*fmatrix[i])[ii][j] = 0.;              //        (*fmatrix[i])[ii][j] /= (*fnmat[i])[ii][j];
891  //          };              //      } else {
892              i1 =  (ii * 191) + j;              //        (*fmatrix[i])[ii][j] = 0.;
893                //      };
894                //      i1 =  (ii * 191) + j;
895                //      i1 =  (ii * 43) + j;
896                i1 =  (ii * 31) + j;
897              //      j1 = -1;              //      j1 = -1;
898              for (Int_t iij=0; iij<43; iij++){              for (Int_t iij=0; iij<43; iij++){
899                for (Int_t jj=0; jj<191; jj++){                //              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                  j1 = (iij * 191) + jj;                  if ( j1 < 0 ) j1 += 1333;
911                    if ( j1 >= 1333 ) j1 -= 1333;
912    
913    //              j1 = (iij * 31) + jj;
914                  //              j1++;                  //              j1++;
915                  //              if ( finmat[ii][j] > 0 ){                  //              if ( finmat[ii][j] > 0 ){
916                  //                (*fmatrix)[i1][j1] /= finmat[ii][j];                  //                (*fmatrix)[i1][j1] /= finmat[ii][j];
917                  if ( (*fnmat[i])[ii][j] == 0. || (*fmatrix)[i1][j1] == 0. || !((*fmatrix)[i1][j1] == (*fmatrix)[i1][j1]) ){                  if ( (*fnmat[i])[ii][j] == 0. || (*fmatrix[i])[i1][j1] == 0. || !((*fmatrix[i])[i1][j1] == (*fmatrix[i])[i1][j1]) ){
918                    (*fmatrix)[i1][j1] = 1.;                    (*fmatrix[i])[i1][j1] = 1.;
919                  } else {                  } else {
920                    (*fmatrix)[i1][j1] /= (*fnmat[i])[ii][j];                    (*fmatrix[i])[i1][j1] /= (*fnmat[i])[ii][j];
921                    nonzero++;                    nonzero++;
922                    if ( i1 == 0 ) nonzero1++;                    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 ){                  //              if ( j>=94 && j <=96 && jj >=94 && jj<=96 ){
958                  //                      mi = (ii*3) + j -94;                  //                      mi = (ii*3) + j -94;
959                  //                      mj = (iij*3) + jj -94;                  //                      mj = (iij*3) + jj -94;
# Line 862  void SaveHistos(){ Line 961  void SaveHistos(){
961                  //              };                  //              };
962    
963    
964                  if ( j>=84 && j <=106 && jj >=84 && jj<=106 ){                  //              if ( j>=84 && j <=106 && jj >=84 && jj<=106 ){
965                          mi = (ii*3) + j -84;                  //                      mi = (ii*3) + j -84;
966                          mj = (iij*3) + jj -84;                  //                      mj = (iij*3) + jj -84;
967                          (*mymat)[mi][mj] = (*fmatrix)[i1][j1];                  //                      (*mymat)[mi][mj] = (*fmatrix)[i1][j1];
968                  };                  //              };
969    
970                };                };
971              };              };
# Line 874  void SaveHistos(){ Line 973  void SaveHistos(){
973          };          };
974          //          //
975          printf(" Matrix has %i non-zero elements \n",nonzero);            printf(" Matrix has %i non-zero elements \n",nonzero);  
976          printf(" Matrix has %i non-zero elements on the first row\n",nonzero1);          //      printf(" Matrix has %i non-zero elements on the first row\n",nonzero1);
977          //          //
978  //      Bool_t BAD = false;          //      Bool_t BAD = false;
979  //      for (Int_t ii=0; ii<43; ii++){          //      for (Int_t ii=0; ii<43; ii++){
980  //        for (Int_t j=0; j<191; j++){          //        for (Int_t j=0; j<191; j++){
981  //          //          //          //
982  //          i1 =  (ii * 191) + j;          //          i1 =  (ii * 191) + j;
983  //          //          //          //
984  //          for (Int_t iij=0; iij<43; iij++){          //          for (Int_t iij=0; iij<43; iij++){
985  //            for (Int_t jj=0; jj<191; jj++){          //            for (Int_t jj=0; jj<191; jj++){
986  //              //          //              //
987  //              j1 = (iij * 191) + jj;          //              j1 = (iij * 191) + jj;
988  //              //          //              //
989  //              //              printf(" ROW %i COLUMN %i VALUE %f \n",i1,j1,(*fmatrix)[i1][j1]);          //              //              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]) ){          //              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]);          //                printf(" ROW %i COLUMN %i VALUE %f \n",i1,j1,(*fmatrix)[i1][j1]);
992  //                printf(" che schifo! \n");          //                printf(" che schifo! \n");
993  //                BAD = true;          //                BAD = true;
994  //              };          //              };
995  //              //          //              //
996  //            };          //            };
997  //          };          //          };
998  //        };          //        };
999  //      };          //      };
1000  //      //          //      //
1001  //      if ( BAD ) printf(" questa matrice fa cagare \n");          //      if ( BAD ) printf(" questa matrice fa cagare \n");
1002          //          //
1003          //          //
1004          //      cf->WriteFullMatrix(fmatrix[i],i);          cf->WriteFullMatrix(fmatrix[i],i);
1005          cf->WriteFullMatrix(fmatrix, i);          //      cf->WriteFullMatrix(fmatrix, i);
1006          //      cf->WriteFullNMatrix(fnmat, i);          //      cf->WriteFullNMatrix(fnmat, i);
1007          cf->WriteFullNMatrix(fnmat[i], i);          cf->WriteFullNMatrix(fnmat[i], i);
1008          //          //
1009          //      if ( fmatrix[i]->Determinant() == 0. ){          //      TDecompSVD svd(*fmatrix[i]);
1010          if ( fmatrix->Determinant() == 0. ){          //      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);            printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i);
1017          } else {          } else {
1018            Double_t det = 0.;            //    };
1019            //      TMatrixF invmatrix = (TMatrixF)(fmatrix[i]->Invert(&det));            //    if ( i == 3 ){
1020            //      printf(" Bin %i determinant is %f \n",i,det);            //    if ( ok ){
1021            //      cf->WriteInvertedFullMatrix((TMatrixF)invmatrix,i);            //      Double_t tol = 1E-20;
1022            TMatrixF invmatrix = (TMatrixF)(fmatrix->Invert(&det));            //      TDecompSVD svd((*fmatrix)[i],tol);
1023            printf(" Bin %i determinant is %f \n",i,det);            //      svd.Decompose();
1024            cf->WriteInvertedFullMatrix((TMatrixF)invmatrix,i);            //      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 ( mymat->Determinant() == 0. ){  //      if ( mymat3->Determinant() == 0. ){
1036            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);
1037          } else {  //      } else {
1038            Double_t det = 0.;  //        Double_t det = 0.;
1039            TMatrixF invmatrix = (TMatrixF)(mymat->Invert(&det));  //        TMatrixD invmatrix = (TMatrixD)(mymat3->Invert(&det));
1040            printf(" Bin %i determinant is %f \n",i,det);  //        printf(" Mymat3 determinant is %f \n",det);
1041            cf->WriteInvertedFullMatrix((TMatrixF)invmatrix,i);  //        cf->WriteInvertedFullMatrix((TMatrixD)invmatrix,1103);
1042          };  //      };
1043          cf->WriteFullMatrix(mymat, 99);  //      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    //      cf->WriteFullMatrix(mymat5, 105);
1053    //      if ( mymat7->Determinant() == 0. ){
1054    //        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    //        printf(" Mymat7 determinant is %f \n",det);
1059    //        cf->WriteInvertedFullMatrix((TMatrixD)invmatrix,1107);
1060    //      };
1061    //      cf->WriteFullMatrix(mymat7, 107);
1062    //      if ( mymat9->Determinant() == 0. ){
1063    //        printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i);
1064    //      } else {
1065    //        Double_t det = 0.;
1066    //        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);          //      cf->UnLoadFullMatrix(i);
1093          //      cf->UnLoadFullNMatrix(i);          //      cf->UnLoadFullNMatrix(i);
1094          delete fmatrix;          //      delete fmatrix;
1095          //      delete fnmat;          //      delete fnmat;
1096          //          //
1097        };        };
Legend:
Removed from v.1.6  
changed lines
  Added in v.1.7
  ViewVC Help
Powered by ViewVC 1.1.23