/[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.8 by mocchiut, Fri Jan 25 15:09:07 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 58  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    //------------------------------------------------------------------    //------------------------------------------------------------------
# Line 251  bool Select( PamLevel2* event ){ Line 252  bool Select( PamLevel2* event ){
252    };    };
253    //  if ( event->GetCaloLevel1()->qtotpl(0) > 7. ) return false;    //  if ( event->GetCaloLevel1()->qtotpl(0) > 7. ) return false;
254    
255    if ( rigidity > 2.2 || rigidity < 1.5 ) return false;    //  if ( rigidity > 2.2 || rigidity < 1.5 ) return false;
256    //  printf(" rig %f CRIG %i SRIG %i \n",rigidity,CRIG,SRIG);    //  printf(" rig %f CRIG %i SRIG %i \n",rigidity,CRIG,SRIG);
257    //    //
258    return true;    return true;
# Line 264  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());
     cf->LoadBin();  
279      if ( !FULL ){      if ( !FULL ){
280          cf->LoadBin();
281        cf->LoadLong();        cf->LoadLong();
282      } else {      } else {
283          cf->LoadBin(true);
284        cf->LoadFull();        cf->LoadFull();
285      };      };
286    } else {    } else {
# Line 304  void CreateHistos( PamLevel2* event , TS Line 311  void CreateHistos( PamLevel2* event , TS
311    memset(rmean, 0, 17*sizeof(Float_t));    memset(rmean, 0, 17*sizeof(Float_t));
312    memset(ntot, 0, 17*sizeof(Int_t));    memset(ntot, 0, 17*sizeof(Int_t));
313    //  memset(finmat, 0, 43*191*sizeof(Int_t));      //  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      for (Int_t i=3; i < 4 ; i++){      //  for (Int_t i=3; i < 4 ; i++){
318      if ( !FULL ){      if ( !FULL ){
319        matrix[i] = new TMatrixD(43,43);        matrix[i] = new TMatrixD(43,43);
320        qplane[i] = new TArrayF(43);        qplane[i] = new TArrayF(43);
# Line 318  void CreateHistos( PamLevel2* event , TS Line 326  void CreateHistos( PamLevel2* event , TS
326          //      fnmat = new TMatrixF(4128,4128);          //      fnmat = new TMatrixF(4128,4128);
327          //      fmatrix = new TMatrixF(8213,8213);          //      fmatrix = new TMatrixF(8213,8213);
328          //      fnmat = new TMatrixF(8213,8213);          //      fnmat = new TMatrixF(8213,8213);
329          fmatrix = new TMatrixF(MDIM,MDIM);          //      fmatrix = new TMatrixF(MDIM,MDIM);
330          //      fnmat = new TMatrixF(MDIM,MDIM);          //      fnmat = new TMatrixF(MDIM,MDIM);
331          fnmat[i] = new TMatrixF(43,191);          //      fmatrix[i] = new TMatrixF(1849,1849);
332  //      cf->WriteFullMatrix(fmatrix, i);          //      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);          //      cf->WriteFullNMatrix(fnmat, i);
340  //      delete fmatrix;          //      delete fmatrix;
341          //      delete fnmat;          //      delete fnmat;
342          //fnmat[i] = new TMatrixI(8213,8213);          //fnmat[i] = new TMatrixI(8213,8213);
343        } else {        } 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)          //      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);//          //      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          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
348          cf->WriteFullNMean(fnqplane, i);          //
349          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 + ...
350          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  ...
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      };      };
# Line 394  void FindAverage( PamLevel2* L2, int iev Line 415  void FindAverage( PamLevel2* L2, int iev
415      //      //
416      // FULL CALORIMETER      // FULL CALORIMETER
417      //      //
418      fqplane = cf->LoadFullAverage(rbi);      //    fqplane = cf->LoadFullAverage(rbi);
419      fnqplane = cf->LoadFullNAverage(rbi);      //    fnqplane = cf->LoadFullNAverage(rbi);
420      CaloTrkVar *ct = L2->GetTrack(0)->GetCaloTrack();      CaloTrkVar *ct = L2->GetTrack(0)->GetCaloTrack();
421      //      //
422      Int_t nplane = 0;      Int_t nplane = 0;
# Line 409  void FindAverage( PamLevel2* L2, int iev Line 430  void FindAverage( PamLevel2* L2, int iev
430      Int_t mstrip = 0;      Int_t mstrip = 0;
431      //      //
432      for (Int_t j=0; j<2; j++){      for (Int_t j=0; j<2; j++){
433        for (Int_t i=0; i<21; i++){        for (Int_t i=0; i<22; i++){
434          nplane = 1 - j + 2*i;          nplane = 1 - j + 2*i;
435          if ( nplane > 37 ) nplane--;          if ( nplane > 37 ) nplane--;
436          //          //
# Line 417  void FindAverage( PamLevel2* L2, int iev Line 438  void FindAverage( PamLevel2* L2, int iev
438          //          //
439          cd = 95 - cs;            cd = 95 - cs;  
440          //          //
441            Int_t oldstr = -1;
442          for (Int_t k=0; k<191; k++){          for (Int_t k=0; k<191; k++){
443            mstrip = cd + k;            mstrip = cd + k;
444            //      if ( mstrip < (191-cs) ) (*fnqplane[rbi])[nplane][mstrip] += 1.;            //      if ( mstrip < (191-cs) ) (*fnqplane[rbi])[nplane][mstrip] += 1.;
445            if ( mstrip < (191-cs) ) (*fnqplane)[nplane][mstrip] += 1.;            //      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      };      };
# Line 439  void FindAverage( PamLevel2* L2, int iev Line 466  void FindAverage( PamLevel2* L2, int iev
466        //        //
467        mstrip = cd + strip;        mstrip = cd + strip;
468        //        //
469          Int_t lstr = cf->ConvertStrip(mstrip);
470        //      (*fqplane[rbi])[nplane][mstrip] += mip;        //      (*fqplane[rbi])[nplane][mstrip] += mip;
471        (*fqplane)[nplane][mstrip] += mip;        //      (*fqplane)[nplane][mstrip] += mip;
472          (*fqplane[rbi])[nplane][lstr] += mip;
473        //        //
474      };      };
475      //      //
476      cf->WriteFullMean(fqplane, rbi);      //    cf->WriteFullMean(fqplane, rbi);
477      cf->WriteFullNMean(fnqplane, rbi);      //    cf->WriteFullNMean(fnqplane, rbi);
478      cf->UnLoadFullAverage(rbi);      //    cf->UnLoadFullAverage(rbi);
479      cf->UnLoadFullNAverage(rbi);      //    cf->UnLoadFullNAverage(rbi);
480      delete fqplane;      //    delete fqplane;
481      delete fnqplane;      //    delete fnqplane;
482      //      //
483    };    };
484  }  }
# Line 476  void CalculateAverage(){ Line 505  void CalculateAverage(){
505    } else {    } else {
506      //      //
507      for (Int_t i=0; i<nbin-1; i++){      for (Int_t i=0; i<nbin-1; i++){
508        fqplane = cf->LoadFullAverage(i);        //      fqplane = cf->LoadFullAverage(i);
509        fnqplane = cf->LoadFullNAverage(i);        //      fnqplane = cf->LoadFullNAverage(i);
510        if ( ntot[i] > 0 ) rmean[i] /= (Float_t)(ntot[i]);        if ( ntot[i] > 0 ) rmean[i] /= (Float_t)(ntot[i]);
511        //        //
512        for (Int_t j=0; j<43 ; j++){        for (Int_t j=0; j<43 ; j++){
513          for (Int_t k=0; k<191; k++){          //      for (Int_t k=0; k<191; k++){
514  //        if ( (*fnqplane[i])[j][k] > 0 ){            //      for (Int_t k=0; k<43; k++){
515  //          (*fqplane[i])[j][k] /= (Float_t)(*fnqplane[i])[j][k];          //      for (Int_t k=0; k<31; k++){
516  //        } else {          for (Int_t k=0; k<11; k++){
517  //          (*fqplane[i])[j][k] = 0.;            //      if ( (*fnqplane[i])[j][k] > 0 ){  
518  //        };            //        (*fqplane[i])[j][k] /= (Float_t)(*fnqplane[i])[j][k];
519  //        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]);            //      } else {
520            if ( (*fnqplane)[j][k] > 0 ){              //        (*fqplane[i])[j][k] = 0.;
521              if ( (*fqplane)[j][k] == 0. ) (*fqplane)[j][k] = 0.7;            //      };
522              (*fqplane)[j][k] /= (Float_t)(*fnqplane)[j][k];            //      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 {            } else {
527              (*fqplane)[j][k] = 0.;              (*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]);            //      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        cf->WriteFullMean(fqplane, i);        cf->WriteFullMean(fqplane[i], i);
533        cf->WriteFullNMean(fnqplane, i);        cf->WriteFullNMean(fnqplane[i], i);
534        cf->UnLoadFullAverage(i);        //      cf->UnLoadFullAverage(i);
535        cf->UnLoadFullNAverage(i);        //      cf->UnLoadFullNAverage(i);
536        delete fqplane;        //      delete fqplane;
537        delete fnqplane;        //      delete fnqplane;
538      };      };
539      //      //
540  //     for (Int_t i=0; i<nbin-1; i++){      //     for (Int_t i=0; i<nbin-1; i++){
541  //       //      //       //
542  //       cf->WriteFullMean(fqplane[i], i);      //       cf->WriteFullMean(fqplane[i], i);
543  //       //      //       //
544  //     };      //     };
545    };    };
546    //      //  
547    cf->WriteNumBin(nbin);    cf->WriteNumBin(nbin);
# Line 518  void CalculateAverage(){ Line 550  void CalculateAverage(){
550    cf->WriteRigBin(rigbin);    cf->WriteRigBin(rigbin);
551    //    //
552    TArrayF *rmeanbin = new TArrayF(17, rmean);    TArrayF *rmeanbin = new TArrayF(17, rmean);
553    TFile *file =  cf->GetFile();    if ( FULL ){
554    file->cd();      TFile *file =  cf->GetFullFile();
555    file->WriteObject(&(*rmeanbin), "binrigmean");      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  }  }
# Line 540  void FindMatrix( PamLevel2* L2, int iev Line 578  void FindMatrix( PamLevel2* L2, int iev
578      };      };
579    };    };
580    //    //
 if ( rbi != 3 ) return;  
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    //    //
# Line 583  if ( rbi != 3 ) return; Line 620  if ( rbi != 3 ) return;
620      //      //
621      // FULL CALORIMETER      // FULL CALORIMETER
622      //      //
623      printf(" Retrieve matrix %i IEV %i \n",rbi,iev);      //    if ( rbi != 3 ) return;
624        printf(" matrix %i IEV %i \n",rbi,iev);
625      //    fmatrix = cf->LoadFullMatrix(rbi);      //    fmatrix = cf->LoadFullMatrix(rbi);
626  //    cf->LoadFullMatrix(rbi,fmatrix);      //    cf->LoadFullMatrix(rbi,fmatrix);
627      //    fnmat = cf->LoadFullNMatrix(rbi);      //    fnmat = cf->LoadFullNMatrix(rbi);
628      printf(" done \n");      //    printf(" done \n");
629      printf(" start loop \n");      //    printf(" start loop \n");
630      //      //
631      CaloTrkVar *ct = L2->GetTrack(0)->GetCaloTrack();      CaloTrkVar *ct = L2->GetTrack(0)->GetCaloTrack();
632      //      //
# Line 600  if ( rbi != 3 ) return; Line 638  if ( rbi != 3 ) return;
638      //      //
639      //    Int_t mindgf = 48;      //    Int_t mindgf = 48;
640      //    Int_t dgf = 143;      //    Int_t dgf = 143;
641  //       Int_t mindgf = 0; //tutto      //       Int_t mindgf = 0; //tutto
642  //        Int_t dgf = 191; //tutto      //        Int_t dgf = 191; //tutto
643  //    Int_t mindgf = 94;      //    Int_t mindgf = 94;
644  //    Int_t dgf = 96;      //    Int_t dgf = 96;
645      Int_t mindgf = 84;      //    Int_t mindgf = 84;
646      Int_t dgf = 106;      //    Int_t dgf = 106;
647        Int_t mindgf = 0;
648        //    Int_t dgf = 43;
649        Int_t dgf = 11;
650      Int_t cs = 0;      Int_t cs = 0;
651      Int_t cd = 0;      Int_t cd = 0;
652      Int_t mstrip = 0;      Int_t mstrip = 0;
653      //      //
654      Float_t mipv[43][191];      //    Float_t mipv[43][43];
655      memset(mipv,0,43*191*sizeof(Float_t));      //    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++ ){      for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
662        //        //
# Line 626  if ( rbi != 3 ) return; Line 671  if ( rbi != 3 ) return;
671        //        //
672        mstrip = cd + strip;        mstrip = cd + strip;
673        //        //
674        mipv[nplane][mstrip] = mip;        Int_t lstr = cf->ConvertStrip(mstrip);
675          mipv[nplane][lstr] += mip;
676        //        //
677      };      };
678      //      //
# Line 651  if ( rbi != 3 ) return; Line 697  if ( rbi != 3 ) return;
697      //      //
698      Int_t toto = 0;      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++){      for (Int_t nplane1=0; nplane1<43; nplane1++){
712        if ( nplane1 >= 37 ) nn1 = nplane1 + 1;        if ( nplane1 >= 37 ) nn1 = nplane1 + 1;
713        vi1 = 1;        vi1 = 1;
# Line 661  if ( rbi != 3 ) return; Line 718  if ( rbi != 3 ) return;
718        //        //
719        cd1 = 95 - cs1;        cd1 = 95 - cs1;
720        //        //
721        mstrip1min = TMath::Max(mindgf,(cd1+0));        Int_t at1 = TMath::Max(0,(cd1+0));
722        mstrip1max = TMath::Min(dgf,(cd1+95)) + 1;        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 \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);
729        //        //
730        for (Int_t mstrip1=mstrip1min; mstrip1<mstrip1max; mstrip1++){        for (Int_t mstrip1=mstrip1min; mstrip1<mstrip1max; mstrip1++){
731          //      printf(".\n");          //      printf(".\n");
732          //          //
733          mj = -1;          mj = -1;
734          //          //
735          mip1 = mipv[nplane1][mstrip1] - cf->GetFullAverageAt(nplane1,mstrip1,erig,rbi);          mip1 = mipv[nplane1][mstrip1] - cf->GetFullAverageAt(nplane1,mstrip1,erig,therigb);
736          //          //
737          mi = (nplane1 * 191) + mstrip1;          //      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;          //      if ( mstrip1 > mstrip1min ) break;
743          //      if ( mstrip1 > dgf ) break;          //      if ( mstrip1 > dgf ) break;
# Line 697  if ( rbi != 3 ) return; Line 761  if ( rbi != 3 ) return;
761              //              //
762              //    mstrip2min = cd2 + 0;              //    mstrip2min = cd2 + 0;
763              //    mstrip2max = cd2 + 95;              //    mstrip2max = cd2 + 95;
764              mstrip2min = TMath::Max(mindgf,(cd2+0));              Int_t t1 = TMath::Max(0,(cd2+0));
765              mstrip2max = TMath::Min(dgf,(cd2+95)) + 1;              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);              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++){              for (Int_t mstrip2=mstrip2min; mstrip2<mstrip2max; mstrip2++){
772                //                //              
773                mip2 = mipv[nplane2][mstrip2] - cf->GetFullAverageAt(nplane2,mstrip2,erig,rbi);                mip2 = mipv[nplane2][mstrip2] - cf->GetFullAverageAt(nplane2,mstrip2,erig,therigb);
774                //                //
775                if ( mip2 != 0. ){                if ( mip2 != 0. ){
776                  //                  //
777                  mj = (nplane2 * 191) + mstrip2;                  //              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++;                  //            mj++;
805                  //                  //
806                  //            if ( mstrip2 > mstrip2min ) break;                  //            if ( mstrip2 > mstrip2min ) break;
# Line 717  if ( rbi != 3 ) return; Line 809  if ( rbi != 3 ) return;
809                  //              if ( mstrip1 >= mstrip1min && mstrip1 <= mstrip1max && mstrip2 >= mstrip2min && mstrip2 <= mstrip2max){                  //              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));                  //              (*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.;                  //                (*fnmat[rbi])[mi][mj] += 1.;
812                  (*fmatrix)[mi][mj] += (mip1 * mip2); // giusto                  (*fmatrix[rbi])[mi][mj] += (mip1 * mip2); // giusto
813  //              (*fmatrix)[mi][mj] += (mip1 * mip2) * 1000000.;                  //              (*fmatrix)[mi][mj] += (mip1 * mip2) * 1000000.;
814                  toto++;                  toto++;
815                  //            (*fmatrix)[mi][mj] += 1.;                  //            (*fmatrix)[mi][mj] += 1.;
816                  //              cf->GetFullAverageAt(nplane1,mstrip1,erig,rbi);                  //              cf->GetFullAverageAt(nplane1,mstrip1,erig,rbi);
# Line 734  if ( rbi != 3 ) return; Line 826  if ( rbi != 3 ) return;
826      //        //  
827      printf(" toto = %i \n",toto);      printf(" toto = %i \n",toto);
828      printf("\n done \n");      printf("\n done \n");
829      printf(" write matrix \n");      //    printf(" write matrix \n");
830  //    cf->WriteFullMatrix(fmatrix, rbi);      //    cf->WriteFullMatrix(fmatrix, rbi);
831      //    cf->WriteFullNMatrix(fnmat, rbi);      //    cf->WriteFullNMatrix(fnmat, rbi);
832      printf(" done \n");      //    printf(" done \n");
833      printf(" unload matrix \n");      //    printf(" unload matrix \n");
834   //   cf->UnLoadFullMatrix(rbi);      //   cf->UnLoadFullMatrix(rbi);
835      //    cf->UnLoadFullNMatrix(rbi);      //    cf->UnLoadFullNMatrix(rbi);
836      printf(" done \n");      //    printf(" done \n");
837      printf(" delete matrix \n");      //    printf(" delete matrix \n");
838  //    delete fmatrix;      //    delete fmatrix;
839      //    delete fnmat;      //    delete fnmat;
840      printf(" done \n");      //    printf(" done \n");
841    };    };
842  }  }
843    
# Line 801  void SaveHistos(){ Line 893  void SaveHistos(){
893        //        //
894        // FULL        // FULL
895        //        //
896        //      for (Int_t i=0; i<nbin-1; i++){        for (Int_t i=0; i<nbin-1; i++){
897  //      for (Int_t i=3; i<5; i++){          //      for (Int_t i=3; i<5; i++){
898        for (Int_t i=3; i<4; i++){          //      for (Int_t i=3; i<4; i++){
899          //          //
900          // determine the average matrix          // determine the average matrix
901          //              //    
902  //      fmatrix = cf->LoadFullMatrix(i);          //      fmatrix = cf->LoadFullMatrix(i);
903          //      fnmat = cf->LoadFullNMatrix(i);          //      fnmat = cf->LoadFullNMatrix(i);
904          //          //
905  //      for (Int_t ii=0; ii<MDIM; ii++){          //      for (Int_t ii=0; ii<MDIM; ii++){
906  //        for (Int_t j=0; j<MDIM; j++){          //        for (Int_t j=0; j<MDIM; j++){
907  // //       if ( (*fnmat[i])[ii][j] > 0. ){          // //       if ( (*fnmat[i])[ii][j] > 0. ){
908  // //         (*fmatrix[i])[ii][j] /= (*fnmat[i])[ii][j];          // //         (*fmatrix[i])[ii][j] /= (*fnmat[i])[ii][j];
909  // //       } else {          // //       } else {
910  // //         (*fmatrix[i])[ii][j] = 0.;          // //         (*fmatrix[i])[ii][j] = 0.;
911  // //       };          // //       };
912  //          if ( (*fnmat)[ii][j] > 0. ){          //          if ( (*fnmat)[ii][j] > 0. ){
913  //            (*fmatrix)[ii][j] /= (*fnmat)[ii][j];          //            (*fmatrix)[ii][j] /= (*fnmat)[ii][j];
914  //          } else {          //          } else {
915  //            (*fmatrix)[ii][j] = 0.;          //            (*fmatrix)[ii][j] = 0.;
916  //          };          //          };
917  //        };          //        };
918  //      };          //      };
919  //          //
920  //      TMatrixF *mymat = new TMatrixF(129,129);  //      TMatrixD *mymat3 = new TMatrixD(129,129);
921          TMatrixF *mymat = new TMatrixF(989,989);  //      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;          Int_t i1 = -1;
929          Int_t j1 = -1;          Int_t j1 = -1;
930          int mi,mj;          //      int mi,mj;
931          Int_t nonzero = 0;          Int_t nonzero = 0;
932          Int_t nonzero1 = 0;          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<191; j++){            //      for (Int_t j=0; j<191; j++){
935  //          if ( (*fnmat[i])[ii][j] > 0. ){            //      for (Int_t j=0; j<43; j++){
936  //            (*fmatrix[i])[ii][j] /= (*fnmat[i])[ii][j];            //      for (Int_t j=0; j<31; j++){
937  //          } else {            for (Int_t j=0; j<11; j++){
938  //            (*fmatrix[i])[ii][j] = 0.;              //      if ( (*fnmat[i])[ii][j] > 0. ){
939  //          };              //        (*fmatrix[i])[ii][j] /= (*fnmat[i])[ii][j];
940              i1 =  (ii * 191) + j;              //      } 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;              //      j1 = -1;
948              for (Int_t iij=0; iij<43; iij++){              for (Int_t iij=0; iij<43; iij++){
949                for (Int_t jj=0; jj<191; jj++){                //              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;                  //              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++;                  //              j1++;
976                  //              if ( finmat[ii][j] > 0 ){                  //              if ( finmat[ii][j] > 0 ){
977                  //                (*fmatrix)[i1][j1] /= finmat[ii][j];                  //                (*fmatrix)[i1][j1] /= finmat[ii][j];
978                  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]) ){
979                    (*fmatrix)[i1][j1] = 1.;                    (*fmatrix[i])[i1][j1] = 1.;
980                  } else {                  } else {
981                    (*fmatrix)[i1][j1] /= (*fnmat[i])[ii][j];                    (*fmatrix[i])[i1][j1] /= (*fnmat[i])[ii][j];
982                    nonzero++;                    nonzero++;
983                    if ( i1 == 0 ) nonzero1++;                    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 ){                  //              if ( j>=94 && j <=96 && jj >=94 && jj<=96 ){
1019                  //                      mi = (ii*3) + j -94;                  //                      mi = (ii*3) + j -94;
1020                  //                      mj = (iij*3) + jj -94;                  //                      mj = (iij*3) + jj -94;
# Line 862  void SaveHistos(){ Line 1022  void SaveHistos(){
1022                  //              };                  //              };
1023    
1024    
1025                  if ( j>=84 && j <=106 && jj >=84 && jj<=106 ){                  //              if ( j>=84 && j <=106 && jj >=84 && jj<=106 ){
1026                          mi = (ii*3) + j -84;                  //                      mi = (ii*3) + j -84;
1027                          mj = (iij*3) + jj -84;                  //                      mj = (iij*3) + jj -84;
1028                          (*mymat)[mi][mj] = (*fmatrix)[i1][j1];                  //                      (*mymat)[mi][mj] = (*fmatrix)[i1][j1];
1029                  };                  //              };
1030    
1031                };                };
1032              };              };
# Line 874  void SaveHistos(){ Line 1034  void SaveHistos(){
1034          };          };
1035          //          //
1036          printf(" Matrix has %i non-zero elements \n",nonzero);            printf(" Matrix has %i non-zero elements \n",nonzero);  
1037          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);
1038          //          //
1039  //      Bool_t BAD = false;          //      Bool_t BAD = false;
1040  //      for (Int_t ii=0; ii<43; ii++){          //      for (Int_t ii=0; ii<43; ii++){
1041  //        for (Int_t j=0; j<191; j++){          //        for (Int_t j=0; j<191; j++){
1042  //          //          //          //
1043  //          i1 =  (ii * 191) + j;          //          i1 =  (ii * 191) + j;
1044  //          //          //          //
1045  //          for (Int_t iij=0; iij<43; iij++){          //          for (Int_t iij=0; iij<43; iij++){
1046  //            for (Int_t jj=0; jj<191; jj++){          //            for (Int_t jj=0; jj<191; jj++){
1047  //              //          //              //
1048  //              j1 = (iij * 191) + jj;          //              j1 = (iij * 191) + jj;
1049  //              //          //              //
1050  //              //              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]);
1051  //              if ( (*fmatrix)[i1][j1] == 0. || !((*fmatrix)[i1][j1]==(*fmatrix)[i1][j1]) ){          //              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]);          //                printf(" ROW %i COLUMN %i VALUE %f \n",i1,j1,(*fmatrix)[i1][j1]);
1053  //                printf(" che schifo! \n");          //                printf(" che schifo! \n");
1054  //                BAD = true;          //                BAD = true;
1055  //              };          //              };
1056  //              //          //              //
1057  //            };          //            };
1058  //          };          //          };
1059  //        };          //        };
1060  //      };          //      };
1061  //      //          //      //
1062  //      if ( BAD ) printf(" questa matrice fa cagare \n");          //      if ( BAD ) printf(" questa matrice fa cagare \n");
1063          //          //
1064          //          //
1065          //      cf->WriteFullMatrix(fmatrix[i],i);          //      cf->WriteFullMatrix(fmatrix, i);
         cf->WriteFullMatrix(fmatrix, i);  
1066          //      cf->WriteFullNMatrix(fnmat, i);          //      cf->WriteFullNMatrix(fnmat, i);
1067          cf->WriteFullNMatrix(fnmat[i], i);          //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          //      if ( fmatrix[i]->Determinant() == 0. ){          Double_t zero = (Double_t)0.0;
1074          if ( fmatrix->Determinant() == 0. ){          //
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);            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 {          } else {
1081            Double_t det = 0.;            //    };
1082            //      TMatrixF invmatrix = (TMatrixF)(fmatrix[i]->Invert(&det));            //    if ( i == 3 ){
1083            //      printf(" Bin %i determinant is %f \n",i,det);            //    if ( ok ){
1084            //      cf->WriteInvertedFullMatrix((TMatrixF)invmatrix,i);            //      Double_t tol = 1E-20;
1085            TMatrixF invmatrix = (TMatrixF)(fmatrix->Invert(&det));            //      TDecompSVD svd((*fmatrix)[i],tol);
1086            printf(" Bin %i determinant is %f \n",i,det);            //      svd.Decompose();
1087            cf->WriteInvertedFullMatrix((TMatrixF)invmatrix,i);            //      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 ( mymat->Determinant() == 0. ){  //      if ( mymat3->Determinant() == 0. ){
1099            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);
1100          } else {  //      } else {
1101            Double_t det = 0.;  //        Double_t det = 0.;
1102            TMatrixF invmatrix = (TMatrixF)(mymat->Invert(&det));  //        TMatrixD invmatrix = (TMatrixD)(mymat3->Invert(&det));
1103            printf(" Bin %i determinant is %f \n",i,det);  //        printf(" Mymat3 determinant is %f \n",det);
1104            cf->WriteInvertedFullMatrix((TMatrixF)invmatrix,i);  //        cf->WriteInvertedFullMatrix((TMatrixD)invmatrix,1103);
1105          };  //      };
1106          cf->WriteFullMatrix(mymat, 99);  //      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);          //      cf->UnLoadFullMatrix(i);
1156          //      cf->UnLoadFullNMatrix(i);          //      cf->UnLoadFullNMatrix(i);
1157          delete fmatrix;          //      delete fmatrix;
1158          //      delete fnmat;          //      delete fnmat;
1159          //          //
1160        };        };
Legend:
Removed from v.1.6  
changed lines
  Added in v.1.8
  ViewVC Help
Powered by ViewVC 1.1.23