/[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.7 by mocchiut, Mon Jan 21 10:24:10 2008 UTC revision 1.8 by mocchiut, Fri Jan 25 15:09:07 2008 UTC
# 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 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 323  void CreateHistos( PamLevel2* event , TS Line 330  void CreateHistos( PamLevel2* event , TS
330          //      fnmat = new TMatrixF(MDIM,MDIM);          //      fnmat = new TMatrixF(MDIM,MDIM);
331          //      fmatrix[i] = new TMatrixF(1849,1849);          //      fmatrix[i] = new TMatrixF(1849,1849);
332          //      fnmat[i] = new TMatrixF(43,43);          //      fnmat[i] = new TMatrixF(43,43);
333          fmatrix[i] = new TMatrixD(1333,1333);          //      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);  //      fmatrix[i]->SetTol(tol);
         fnmat[i] = new TMatrixF(43,31);  
338          //      cf->WriteFullMatrix(fmatrix, i);          //      cf->WriteFullMatrix(fmatrix, i);
339          //      cf->WriteFullNMatrix(fnmat, i);          //      cf->WriteFullNMatrix(fnmat, i);
340          //      delete fmatrix;          //      delete fmatrix;
# Line 337  void CreateHistos( PamLevel2* event , TS Line 346  void CreateHistos( PamLevel2* event , TS
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] + ...          //      fqplane[i] = new TMatrixD(43,43); // 43 planes x 43 "strip", where 43 = 50 + 14 + 5 + 2 + 2 + 2 + 2 + 2 + 2 + 2 + 2 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + [1] + ...
347          //      fnqplane[i] = new TMatrixD(43,43);//                                     0    1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20   21          //      fnqplane[i] = new TMatrixD(43,43);//                                     0    1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20   21
348          //          //
349          fqplane[i] = new TMatrixD(43,31); // 43 planes x 43 "strip", where 43 = 50 + 14 + 6 + 5 + 3 + 3 + 3 + 2 + 2 + 2 + 1 + 1 + 1 + 1 + 1 + [1] + 1 + 1 + 1 + 1 + ...          //      fqplane[i] = new TMatrixD(43,31); // 43 planes x 43 "strip", where 43 = 50 + 14 + 6 + 5 + 3 + 3 + 3 + 2 + 2 + 2 + 1 + 1 + 1 + 1 + 1 + [1] + 1 + 1 + 1 + 1 + ...
350          fnqplane[i] = new TMatrixD(43,31);//                                     0    1   2   3   4   5   6   7   8   9  10  11  12  13  14    15  16  17  18  19  ...          //      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);          //      cf->WriteFullMean(fqplane, i);
355          //      cf->WriteFullNMean(fnqplane, i);          //      cf->WriteFullNMean(fnqplane, i);
# Line 501  void CalculateAverage(){ Line 512  void CalculateAverage(){
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          //      for (Int_t k=0; k<43; k++){          //      for (Int_t k=0; k<43; k++){
515          for (Int_t k=0; k<31; k++){          //      for (Int_t k=0; k<31; k++){
516            for (Int_t k=0; k<11; k++){
517            //      if ( (*fnqplane[i])[j][k] > 0 ){              //      if ( (*fnqplane[i])[j][k] > 0 ){  
518            //        (*fqplane[i])[j][k] /= (Float_t)(*fnqplane[i])[j][k];            //        (*fqplane[i])[j][k] /= (Float_t)(*fnqplane[i])[j][k];
519            //      } else {            //      } else {
# Line 538  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 627  void FindMatrix( PamLevel2* L2, int iev Line 645  void FindMatrix( PamLevel2* L2, int iev
645      //    Int_t mindgf = 84;      //    Int_t mindgf = 84;
646      //    Int_t dgf = 106;      //    Int_t dgf = 106;
647      Int_t mindgf = 0;      Int_t mindgf = 0;
648      Int_t dgf = 43;      //    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][43];      //    Float_t mipv[43][43];
655      //    memset(mipv,0,43*43*sizeof(Float_t));      //    memset(mipv,0,43*43*sizeof(Float_t));
656      Float_t mipv[43][31];      //    Float_t mipv[43][31];
657      memset(mipv,0,43*31*sizeof(Float_t));      //    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 676  void FindMatrix( PamLevel2* L2, int iev Line 697  void FindMatrix( PamLevel2* L2, int iev
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 700  void FindMatrix( PamLevel2* L2, int iev Line 732  void FindMatrix( PamLevel2* L2, int iev
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;          //      mi = (nplane1 * 43) + mstrip1;
739          mi = (nplane1 * 31) + mstrip1;          //      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 737  void FindMatrix( PamLevel2* L2, int iev Line 770  void FindMatrix( PamLevel2* L2, int iev
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;                  //              mj = (nplane2 * 43) + mstrip2;
779  //              mj = (nplane2 * 31) + mstrip2;  //              mj = (nplane2 * 31) + mstrip2;
780                  Int_t sh = -15 + nplane1;  //              Int_t sh = -15 + nplane1;
781                  if ( sh > 15 ) sh -= 31*nplane1;  //              if ( sh > 15 ) sh -= 31*nplane1;
782                  //  //              //
783                  mj = (nplane2 * 31) + mstrip2 + sh;  //              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                  //                  //
                 if ( mj < 0 ) mj += 1333;  
                 if ( mj >= 1333 ) mj -= 1333;  
802  //              printf(" mi %i mj %i sh %i \n",mi,mj,sh);  //              printf(" mi %i mj %i sh %i \n",mi,mj,sh);
803                  //                  //
804                  //            mj++;                  //            mj++;
# Line 885  void SaveHistos(){ Line 933  void SaveHistos(){
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            //      for (Int_t j=0; j<43; j++){            //      for (Int_t j=0; j<43; j++){
936            for (Int_t j=0; j<31; j++){            //      for (Int_t j=0; j<31; j++){
937              for (Int_t j=0; j<11; j++){
938              //      if ( (*fnmat[i])[ii][j] > 0. ){              //      if ( (*fnmat[i])[ii][j] > 0. ){
939              //        (*fmatrix[i])[ii][j] /= (*fnmat[i])[ii][j];              //        (*fmatrix[i])[ii][j] /= (*fnmat[i])[ii][j];
940              //      } else {              //      } else {
# Line 893  void SaveHistos(){ Line 942  void SaveHistos(){
942              //      };              //      };
943              //      i1 =  (ii * 191) + j;              //      i1 =  (ii * 191) + j;
944              //      i1 =  (ii * 43) + j;              //      i1 =  (ii * 43) + j;
945              i1 =  (ii * 31) + j;              //      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++){                //              for (Int_t jj=0; jj<43; jj++){
951                for (Int_t jj=0; jj<31; jj++){                //              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;                  //              j1 = (iij * 43) + jj;
956                  Int_t sh = -15 + ii;  //              Int_t sh = -15 + ii;
957                  if ( sh > 15 ) sh -= 31*ii;  //              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                  //                  //
                 j1 = (iij * 31) + jj + sh;  
964                  //                  //
965                  if ( j1 < 0 ) j1 += 1333;                  //
966                  if ( j1 >= 1333 ) j1 -= 1333;  //              Int_t sh = -5 + ii;
967    //              if ( sh > 5 ) sh -= 11*ii;
968  //              j1 = (iij * 31) + jj;  //              //
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];
# Line 1001  void SaveHistos(){ Line 1062  void SaveHistos(){
1062          //      if ( BAD ) printf(" questa matrice fa cagare \n");          //      if ( BAD ) printf(" questa matrice fa cagare \n");
1063          //          //
1064          //          //
         cf->WriteFullMatrix(fmatrix[i],i);  
1065          //      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]);          //      TDecompSVD svd(*fmatrix[i]);
1071          //      Bool_t ok = svd.Decompose();          //      Bool_t ok = svd.Decompose();
# Line 1014  void SaveHistos(){ Line 1075  void SaveHistos(){
1075          if ( fmatrix[i]->Determinant() == zero ){          if ( fmatrix[i]->Determinant() == zero ){
1076            //if ( fmatrix->Determinant() == 0. ){            //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            //    };            //    };
1082            //    if ( i == 3 ){            //    if ( i == 3 ){
Legend:
Removed from v.1.7  
changed lines
  Added in v.1.8
  ViewVC Help
Powered by ViewVC 1.1.23