--- calo/flight/CaloFranzini/src/Calib.cpp 2007/12/18 09:55:05 1.2 +++ calo/flight/CaloFranzini/src/Calib.cpp 2008/01/25 15:09:07 1.8 @@ -6,6 +6,7 @@ #include #include #include +#include #include #include // @@ -15,10 +16,17 @@ using namespace std; // extern Bool_t MATRIX; +extern Bool_t FULL; +extern Bool_t SIMU; +extern Bool_t CRIG; +extern Bool_t SRIG; CaloFranzini *cf; Int_t nbin; Float_t rig[18]; Float_t rmean[17]; +Int_t ntot[17]; +Int_t MDIM = 8213; +//Int_t MDIM = 4128; //Float_t qqplane[17][43]; //Int_t nnqplane[17][43]; @@ -27,6 +35,19 @@ TMatrixD *matrix[17]; TMatrixD *nmat[17]; +//TMatrixD *fqplane; +//TMatrixD *fnqplane; +TMatrixD *fqplane[17]; +TMatrixD *fnqplane[17]; +TMatrixD *fmatrix[17]; +//TMatrixF *fnmat[17]; +//TMatrixD *fmatrix; +//TMatrixD *fnmat; +//TMatrixF *fmatrix; +//TMatrixF *fnmat; +TMatrixF *fnmat[17]; +//Int_t finmat[43][191]; + //=============================================================================== bool Select( PamLevel2* event ){ @@ -37,11 +58,14 @@ PamTrack *track = event->GetTrack(0); if(!track)return false; + if ( SIMU ) return true; //------------------------------------------------------------------ // tracker pre-selection //------------------------------------------------------------------ TrkTrack *trk = track->GetTrkTrack(); - + float rigidity = trk->GetRigidity(); + if ( CRIG ) rigidity = event->GetCaloLevel2()->qtot/260.; + if ( SRIG ) rigidity = event->GetGPamela()->P0; bool TRACK__OK = false; if( trk->chi2 >0 && @@ -67,20 +91,19 @@ !event->GetAcLevel2()->CARDhit() && !event->GetAcLevel2()->CAThit() && true ) TOF__OK = true; - if( !TOF__OK )return false; + if( !TOF__OK && !SIMU)return false; //------------------------------------------------------ // no albedo //------------------------------------------------------ - if( track->GetToFTrack()->beta[12]<=0.2 || - track->GetToFTrack()->beta[12] >= 1.5 ) return false; + if( !SIMU && (track->GetToFTrack()->beta[12]<=0.2 || + track->GetToFTrack()->beta[12] >= 1.5) ) return false; //------------------------------------------------------ bool CUT1 = false; if( trk->nstep<100 && - trk->GetRigidity()<400. && - trk->GetRigidity()>0.1 && - trk->GetDeflection()<0. && + rigidity<400. && + rigidity>0.1 && trk->resx[0]<0.001 && trk->resx[5]<0.001 && track->IsSolved() && @@ -88,6 +111,56 @@ true ) CUT1 = true; //------------------------------------------------------ if( !CUT1 )return false; + if ( trk->GetDeflection()>0. && !SIMU ) return false; + + // + // ELENA'S CUT + // + // + // lever-arm 6 + //==================================================== + bool LX6=false; + if( + track->GetTrkTrack()->GetLeverArmX()==6 && + !track->GetTrkTrack()->IsBad(0,0) && + !track->GetTrkTrack()->IsBad(5,0) && + track->GetTrkTrack()->resx[0]<0.001 && + track->GetTrkTrack()->resx[5]<0.001 && + track->GetTrkTrack()->IsInsideCavity() && + true ) LX6 = true; + + //==================================================== + // lever-arm 5 + //==================================================== + bool LX5=false; + if( + track->GetTrkTrack()->GetLeverArmX()==5 && + true ){ + if( + track->GetTrkTrack()->XGood(0) && track->GetTrkTrack()->XGood(4) + ){ + + if( + !track->GetTrkTrack()->IsBad(0,0) && + !track->GetTrkTrack()->IsBad(4,0) && + track->GetTrkTrack()->resx[0]<0.001 && + track->GetTrkTrack()->resx[4]<0.001 && + track->GetTrkTrack()->IsInsideCavity() && + true) LX5 = true; + }else if ( + track->GetTrkTrack()->XGood(1) && track->GetTrkTrack()->XGood(5) + ){ + + if( + !track->GetTrkTrack()->IsBad(1,0) && + !track->GetTrkTrack()->IsBad(5,0) && + track->GetTrkTrack()->resx[1]<0.001 && + track->GetTrkTrack()->resx[5]<0.001 && + track->GetTrkTrack()->IsInsideCavity() && + true) LX5 = true; + } + } + if ( !LX5 && !LX6 ) return false; Float_t defl = trk->GetDeflection(); float p0 = 1.111588e+00; float p1 = 1.707656e+00; @@ -105,6 +178,15 @@ track->GetTrkTrack()->chi2 < chi2m && true ) CUT2 = true; if ( !CUT2 ) return false; + float dedxtrk = trk->GetDEDX(); + // float zcutn = 9. + 20./(rigidity*rigidity); + float zcut2 = 3. + 4.3/(rigidity*rigidity); + float zcut1 = 0.52 + 0.455/(rigidity*rigidity); + Bool_t Z1 = false; + if(dedxtrk > zcut1 && dedxtrk < zcut2){ + Z1=true; + } + if ( !Z1 && !SIMU ) return false; //------------------------------------------------------ // // energy momentum match @@ -116,11 +198,62 @@ // for (Int_t i=0; i < 22; i++){ if ( track->GetCaloTrack()->tibar[i][1] < 0 || track->GetCaloTrack()->tibar[i][0] < 0 ){ - return false; + return false; + }; + }; + // + if ( event->GetCaloLevel2()->qtot == 0. ) return false; + if ( rigidity>5. && track->GetCaloTrack()->qtrack/event->GetCaloLevel2()->qtot < 0.4 ) return false; + if ( rigidity<1. && track->GetToFTrack()->beta[12] < 0.8 ) return false; + if ( rigidity>50. ){ + if ( trk->GetNX()<5 && + trk->GetNY()<4 ) return false; + // + Bool_t sphit = false; + for ( Int_t plane = 0; plane < 6; plane++){ + if ( !trk->XGood(plane) ){ + for (Int_t sing = 0; sing < event->GetTrkLevel2()->nclsx(); sing++){ + TClonesArray &t = *(event->GetTrkLevel2()->SingletX); + TrkSinglet *singlet = (TrkSinglet*)t[sing]; + if ( (singlet->plane-1) == plane ){ + Float_t x = (singlet->coord[0]+singlet->coord[1])/2.; + if ( fabs(track->GetTrkTrack()->xv[plane] - x) < 1. ) sphit = true; + }; + }; + }; + if ( !trk->YGood(plane) ){ + for (Int_t sing = 0; sing < event->GetTrkLevel2()->nclsy(); sing++){ + TClonesArray &t = *(event->GetTrkLevel2()->SingletY); + TrkSinglet *singlet = (TrkSinglet*)t[sing]; + if ( (singlet->plane-1) == plane ){ + Float_t x1 = (singlet->coord[0]); + Float_t x2 = (singlet->coord[1]); + if ( fabs(track->GetTrkTrack()->yv[plane] - x1) < 1. ) sphit = true; + if ( fabs(track->GetTrkTrack()->yv[plane] - x2) < 1. ) sphit = true; + }; + }; + }; }; + if ( sphit ) return false; // spurious hit along the track }; // - if ( event->GetCaloLevel1()->qtotpl(0) > 7. ) return false; + Int_t ti0 = track->GetCaloTrack()->tibar[0][1]-1; + + Int_t view = 0; + Int_t plane = 0; + Int_t strip = 0; + Float_t mip = 0.; + // + for ( Int_t i=0; iGetCaloLevel1()->istrip; i++ ){ + // + mip = event->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip); + if ( view == 1 && plane == 0 && strip == ti0 && mip > 4.) return false; + if ( view == 1 && (plane >0 || strip > ti0) ) break; + }; + // if ( event->GetCaloLevel1()->qtotpl(0) > 7. ) return false; + + // if ( rigidity > 2.2 || rigidity < 1.5 ) return false; + // printf(" rig %f CRIG %i SRIG %i \n",rigidity,CRIG,SRIG); // return true; } @@ -132,13 +265,24 @@ // // //=============================================================== -void CreateHistos( PamLevel2* event , TString file){ +void CreateHistos( PamLevel2* event , TString file){ cf = new CaloFranzini(event); // + if ( FULL ){ + cf->CalculateLongTZeta(false); + cf->CalculateFullTZeta(); + }; + // if ( MATRIX ){ cf->UpdateMatrixFile(file.Data()); - cf->LoadBin(); + if ( !FULL ){ + cf->LoadBin(); + cf->LoadLong(); + } else { + cf->LoadBin(true); + cf->LoadFull(); + }; } else { cf->CreateMatrixFile(file.Data()); }; @@ -165,12 +309,55 @@ rig[17] = 400.; // memset(rmean, 0, 17*sizeof(Float_t)); + memset(ntot, 0, 17*sizeof(Int_t)); + // memset(finmat, 0, 43*191*sizeof(Int_t)); +// Double_t tol = 1E-20; // for (Int_t i=0; i < 17 ; i++){ - matrix[i] = new TMatrixD(43,43); - qplane[i] = new TArrayF(43); - nqplane[i] = new TArrayI(43); - nmat[i] = new TMatrixD(43,43); + // for (Int_t i=3; i < 4 ; i++){ + if ( !FULL ){ + matrix[i] = new TMatrixD(43,43); + qplane[i] = new TArrayF(43); + nqplane[i] = new TArrayI(43); + nmat[i] = new TMatrixD(43,43); + } else { + if ( MATRIX ){ + // fmatrix = new TMatrixF(4128,4128); + // fnmat = new TMatrixF(4128,4128); + // fmatrix = new TMatrixF(8213,8213); + // fnmat = new TMatrixF(8213,8213); + // fmatrix = new TMatrixF(MDIM,MDIM); + // fnmat = new TMatrixF(MDIM,MDIM); + // fmatrix[i] = new TMatrixF(1849,1849); + // fnmat[i] = new TMatrixF(43,43); + // fmatrix[i] = new TMatrixD(1333,1333); + // fnmat[i] = new TMatrixF(43,31); + fmatrix[i] = new TMatrixD(473,473); + fnmat[i] = new TMatrixF(43,11); +// fmatrix[i]->SetTol(tol); + // cf->WriteFullMatrix(fmatrix, i); + // cf->WriteFullNMatrix(fnmat, i); + // delete fmatrix; + // delete fnmat; + //fnmat[i] = new TMatrixI(8213,8213); + } else { + // 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) + // fnqplane = new TMatrixD(43,191);// + // 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] + ... + // 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 + // + // 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 + ... + // 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 ... + fqplane[i] = new TMatrixD(43,11); // 43 planes x 11 "strip", where 43 = 84 + 6 + 3 + 1 + 1 +[1]+ 1 + 1 + 3 + 6 + 84 + fnqplane[i] = new TMatrixD(43,11);// 0 1 2 3 4 5 6 7 8 9 10 + // + // cf->WriteFullMean(fqplane, i); + // cf->WriteFullNMean(fnqplane, i); + // delete fqplane; + // delete fnqplane; + // + }; + }; }; // } @@ -179,6 +366,8 @@ void FindAverage( PamLevel2* L2, int iev ){ // Float_t erig = L2->GetTrack(0)->GetTrkTrack()->GetRigidity(); + if ( SRIG ) erig = L2->GetGPamela()->P0; + if ( CRIG ) erig = L2->GetCaloLevel2()->qtot/260.; // Int_t rbi = 0; for (Int_t i = 0; i= rig[nbin-1] ) return; // rmean[rbi] += erig; + ntot[rbi]++; // - Int_t dgf = 43; - // - for (Int_t i=0; iGetCaloLevel1()->istrip; i++ ){ + if (!FULL ){ + Int_t dgf = 43; + // + for (Int_t i=0; iGetCaloLevel1()->istrip; i++ ){ + // + mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip); + // + nplane = 1 - view + 2 * plane; + if ( erig > 4. && nplane == 0 && mip > 15. ) printf(" IEV %i erig %f OBT %u pkt %u file %s \n",iev,erig,L2->GetOrbitalInfo()->OBT,L2->GetOrbitalInfo()->pkt_num,L2->GetPamTree()->GetFile()->GetName()); + //printf(" IEV %i OBT %u pkt %u file %s \n",iev,L2->GetOrbitalInfo()->OBT,L2->GetOrbitalInfo()->pkt_num,L2->GetPamTree()->GetFile()->GetName()); + if ( nplane > 37 ) nplane--; + if ( nplane < dgf ){ + (*qplane[rbi])[nplane] += mip; + }; + // + }; + } else { + // + // FULL CALORIMETER // - mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip); + // fqplane = cf->LoadFullAverage(rbi); + // fnqplane = cf->LoadFullNAverage(rbi); + CaloTrkVar *ct = L2->GetTrack(0)->GetCaloTrack(); // - nplane = 1 - view + 2 * plane; - if ( nplane > 37 ) nplane--; - if ( nplane < dgf ){ - (*qplane[rbi])[nplane] += mip; + Int_t nplane = 0; + Int_t view = 0; + Int_t plane = 0; + Int_t strip = 0; + Float_t mip = 0.; + // + Int_t cs = 0; + Int_t cd = 0; + Int_t mstrip = 0; + // + for (Int_t j=0; j<2; j++){ + for (Int_t i=0; i<22; i++){ + nplane = 1 - j + 2*i; + if ( nplane > 37 ) nplane--; + // + cs = ct->tibar[i][j] - 1; + // + cd = 95 - cs; + // + Int_t oldstr = -1; + for (Int_t k=0; k<191; k++){ + mstrip = cd + k; + // if ( mstrip < (191-cs) ) (*fnqplane[rbi])[nplane][mstrip] += 1.; + // if ( mstrip < (191-cs) ) (*fnqplane)[nplane][mstrip] += 1.; + Int_t lstr = cf->ConvertStrip(mstrip); + if ( oldstr != lstr ){ + (*fnqplane[rbi])[nplane][lstr] += 1.; + oldstr = lstr; + }; + }; + }; + }; + // + // + for ( Int_t i=0; iGetCaloLevel1()->istrip; i++ ){ + // + mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip); + // + nplane = 1 - view + 2 * plane; + if ( nplane > 37 ) nplane--; + // + cs = ct->tibar[plane][view] - 1; + // + cd = 95 - cs; + // + mstrip = cd + strip; + // + Int_t lstr = cf->ConvertStrip(mstrip); + // (*fqplane[rbi])[nplane][mstrip] += mip; + // (*fqplane)[nplane][mstrip] += mip; + (*fqplane[rbi])[nplane][lstr] += mip; + // }; // + // cf->WriteFullMean(fqplane, rbi); + // cf->WriteFullNMean(fnqplane, rbi); + // cf->UnLoadFullAverage(rbi); + // cf->UnLoadFullNAverage(rbi); + // delete fqplane; + // delete fnqplane; + // }; } void CalculateAverage(){ - for (Int_t i=0; i 0 ) rmean[i] /= (Float_t)(*nqplane[i])[0]; - for (Int_t j=0; j<43 ; j++){ - if ( (*nqplane[i])[j] > 0 ){ - (*qplane[i])[j] /= (Float_t)(*nqplane[i])[j]; - } else { - (*qplane[i])[j] = 0.; + // + if ( !FULL ){ + for (Int_t i=0; i 0 ) rmean[i] /= (Float_t)(*nqplane[i])[0]; + for (Int_t j=0; j<43 ; j++){ + if ( (*nqplane[i])[j] > 0 ){ + (*qplane[i])[j] /= (Float_t)(*nqplane[i])[j]; + } else { + (*qplane[i])[j] = 0.; + }; + printf(" BIN %i plane %i average energy %f qplane %f nqplane %i \n",i,j,rmean[i],(*qplane[i])[j],(*nqplane[i])[j]); + }; + }; + for (Int_t i=0; iWriteLongMean(qplane[i], i); + // + }; + } else { + // + for (Int_t i=0; iLoadFullAverage(i); + // fnqplane = cf->LoadFullNAverage(i); + if ( ntot[i] > 0 ) rmean[i] /= (Float_t)(ntot[i]); + // + for (Int_t j=0; j<43 ; j++){ + // for (Int_t k=0; k<191; k++){ + // for (Int_t k=0; k<43; k++){ + // for (Int_t k=0; k<31; k++){ + for (Int_t k=0; k<11; k++){ + // if ( (*fnqplane[i])[j][k] > 0 ){ + // (*fqplane[i])[j][k] /= (Float_t)(*fnqplane[i])[j][k]; + // } else { + // (*fqplane[i])[j][k] = 0.; + // }; + // 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]); + if ( (*fnqplane[i])[j][k] > 0 ){ + if ( (*fqplane[i])[j][k] == 0. ) (*fqplane[i])[j][k] = 0.7; + (*fqplane[i])[j][k] /= (Float_t)(*fnqplane[i])[j][k]; + } else { + (*fqplane[i])[j][k] = 0.; + }; + // 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 average energy %f qplane %f nqplane %i \n",i,j,rmean[i],(*qplane[i])[j],(*nqplane[i])[j]); + cf->WriteFullMean(fqplane[i], i); + cf->WriteFullNMean(fnqplane[i], i); + // cf->UnLoadFullAverage(i); + // cf->UnLoadFullNAverage(i); + // delete fqplane; + // delete fnqplane; }; + // + // for (Int_t i=0; iWriteFullMean(fqplane[i], i); + // // + // }; }; - + // cf->WriteNumBin(nbin); - + // TArrayF *rigbin = new TArrayF(18, rig); cf->WriteRigBin(rigbin); - - TArrayF *rmeanbin = new TArrayF(17, rmean); - TFile *file = cf->GetFile(); - file->cd(); - file->WriteObject(&(*rmeanbin), "binrigmean"); // - for (Int_t i=0; iWriteLongMean(qplane[i], i); - // + TArrayF *rmeanbin = new TArrayF(17, rmean); + if ( FULL ){ + TFile *file = cf->GetFullFile(); + file->cd(); + file->WriteObject(&(*rmeanbin), "binrigmean"); + } else { + TFile *file = cf->GetFile(); + file->cd(); + file->WriteObject(&(*rmeanbin), "binrigmean"); }; // + // } //=============================================================== void FindMatrix( PamLevel2* L2, int iev ){ // Float_t erig = L2->GetTrack(0)->GetTrkTrack()->GetRigidity(); + if ( SRIG ) erig = L2->GetGPamela()->P0; + if ( CRIG ) erig = L2->GetCaloLevel2()->qtot/260.; // Int_t rbi = 0; for (Int_t i = 0; i= rig[nbin-1] ) return; // - Int_t dgf = 43; -// printf("1qui \n"); - // - for (Int_t i=0; iGetCaloLevel1()->istrip; i++ ){ // - mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip); + // Fill the estrip matrix + // + Int_t nplane = 0; + Int_t view = 0; + Int_t plane = 0; + Int_t strip = 0; + Float_t mip = 0.; + Float_t hpl[43]; + memset(hpl,0,43*sizeof(Float_t)); + for ( Int_t i=0; iGetCaloLevel1()->istrip; i++ ){ + // + mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip); + // + nplane = 1 - view + 2 * plane; + if ( nplane > 37 ) nplane--; + if ( nplane < dgf ){ + hpl[nplane] += mip; + }; + // + }; // - nplane = 1 - view + 2 * plane; - if ( nplane > 37 ) nplane--; - if ( nplane < dgf ){ - hpl[nplane] += mip; + for (Int_t i=0; iGetAverageAt(i,erig)) * (hpl[j] - cf->GetAverageAt(j,erig)); + }; }; + } else { // - }; - // -// printf("3qui \n"); - // - for (Int_t i=0; iGetAverageAt(i,erig)) * (hpl[j] - cf->GetAverageAt(j,erig)); + // FULL CALORIMETER + // + // if ( rbi != 3 ) return; + printf(" matrix %i IEV %i \n",rbi,iev); + // fmatrix = cf->LoadFullMatrix(rbi); + // cf->LoadFullMatrix(rbi,fmatrix); + // fnmat = cf->LoadFullNMatrix(rbi); + // printf(" done \n"); + // printf(" start loop \n"); + // + CaloTrkVar *ct = L2->GetTrack(0)->GetCaloTrack(); + // + Int_t nplane = 0; + Int_t view = 0; + Int_t plane = 0; + Int_t strip = 0; + Float_t mip = 0.; + // + // Int_t mindgf = 48; + // Int_t dgf = 143; + // Int_t mindgf = 0; //tutto + // Int_t dgf = 191; //tutto + // Int_t mindgf = 94; + // Int_t dgf = 96; + // Int_t mindgf = 84; + // Int_t dgf = 106; + Int_t mindgf = 0; + // Int_t dgf = 43; + Int_t dgf = 11; + Int_t cs = 0; + Int_t cd = 0; + Int_t mstrip = 0; + // + // Float_t mipv[43][43]; + // memset(mipv,0,43*43*sizeof(Float_t)); + // Float_t mipv[43][31]; + // memset(mipv,0,43*31*sizeof(Float_t)); + Float_t mipv[43][11]; + memset(mipv,0,43*11*sizeof(Float_t)); + // + for ( Int_t i=0; iGetCaloLevel1()->istrip; i++ ){ + // + mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip); + // + nplane = 1 - view + 2 * plane; + if ( nplane > 37 ) nplane--; + // + cs = ct->tibar[plane][view] - 1; + // + cd = 95 - cs; + // + mstrip = cd + strip; + // + Int_t lstr = cf->ConvertStrip(mstrip); + mipv[nplane][lstr] += mip; + // }; + // + Float_t mip1 = 1.; + Int_t cs1; + Int_t cd1; + Float_t mip2 = 1.; + Int_t cs2; + Int_t cd2; + Int_t mi = -1; + Int_t mj = -1; + Int_t nn1 = 0; + Int_t pl1 = 0; + Int_t vi1 = 0; + Int_t nn2 = 0; + Int_t pl2 = 0; + Int_t vi2 = 0; + Int_t mstrip1min = 0; + Int_t mstrip1max = 0; + Int_t mstrip2min = 0; + Int_t mstrip2max = 0; + // + Int_t toto = 0; + // + // + // + Int_t therigb = rbi; + // + if ( erig < rig[0] ){ + therigb = 0; + }; + if ( erig >= rig[nbin-2] ){ + therigb = nbin - 3; + }; + // + for (Int_t nplane1=0; nplane1<43; nplane1++){ + if ( nplane1 >= 37 ) nn1 = nplane1 + 1; + vi1 = 1; + if ( nn1%2 ) vi1 = 0; + pl1 = (nn1 - 1 + vi1)/2; + // + cs1 = ct->tibar[pl1][vi1] - 1; // convertire nplane in pl1 e vi1 + // + cd1 = 95 - cs1; + // + Int_t at1 = TMath::Max(0,(cd1+0)); + Int_t at2 = TMath::Min(190,(cd1+95)); + mstrip1min = cf->ConvertStrip(at1); + mstrip1max = cf->ConvertStrip(at2) + 1; + // mstrip1min = cf->ConvertStrip(TMath::Max(mindgf,(cd1+0))); + // mstrip1max = cf->ConvertStrip(TMath::Min(dgf,(cd1+95))) + 1; + // + 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); + // + for (Int_t mstrip1=mstrip1min; mstrip1GetFullAverageAt(nplane1,mstrip1,erig,therigb); + // + // mi = (nplane1 * 191) + mstrip1; + // mi = (nplane1 * 43) + mstrip1; + // mi = (nplane1 * 31) + mstrip1; + mi = (nplane1 * 11) + mstrip1; + // + // if ( mstrip1 > mstrip1min ) break; + // if ( mstrip1 > dgf ) break; + // if ( mstrip1 >= mindgf && mstrip1 <= dgf && mstrip1 >= mstrip1min && mstrip1 <= mstrip1max ){ + // + // finmat[nplane1][mstrip1]++; + (*fnmat[rbi])[nplane1][mstrip1] += 1.; + // + if ( mip1 != 0. ){ + // + for (Int_t nplane2=0; nplane2<43; nplane2++){ + // + if ( nplane2 >= 37 ) nn2 = nplane2 + 1; + vi2 = 1; + if ( nn2%2 ) vi2 = 0; + pl1 = (nn2 - 1 + vi2)/2; + // + cs2 = ct->tibar[pl2][vi2] - 1; + // + cd2 = 95 - cs2; + // + // mstrip2min = cd2 + 0; + // mstrip2max = cd2 + 95; + Int_t t1 = TMath::Max(0,(cd2+0)); + Int_t t2 = TMath::Min(190,(cd2+95)); + mstrip2min = cf->ConvertStrip(t1); + mstrip2max = cf->ConvertStrip(t2) + 1; + // + if ( nplane1 == 0 && nplane2 == 0 && mstrip1==mstrip1min ) printf(" mstrip2min %i mstrip2max %i \n",mstrip2min,mstrip2max); + // + for (Int_t mstrip2=mstrip2min; mstrip2GetFullAverageAt(nplane2,mstrip2,erig,therigb); + // + if ( mip2 != 0. ){ + // + // mj = (nplane2 * 191) + mstrip2; + // mj = (nplane2 * 43) + mstrip2; +// mj = (nplane2 * 31) + mstrip2; +// Int_t sh = -15 + nplane1; +// if ( sh > 15 ) sh -= 31*nplane1; +// // +// mj = (nplane2 * 31) + mstrip2 + sh; +// // +// if ( mj < 0 ) mj += 1333; +// if ( mj >= 1333 ) mj -= 1333; +// +// +// +// +// Int_t sh = -5 + nplane1; +// if ( sh > 5 ) sh -= 11*nplane1; +// // +// mj = (nplane2 * 11) + mstrip2 + sh; +// // +// if ( mj < 0 ) mj += 473; +// if ( mj >= 473 ) mj -= 473; +// // + // + mj = (nplane2 * 11) + mstrip2; + // +// printf(" mi %i mj %i sh %i \n",mi,mj,sh); + // + // mj++; + // + // if ( mstrip2 > mstrip2min ) break; + // if ( mstrip2 > dgf ) break; + // if ( mstrip2 >= mindgf && mstrip2 <= dgf && mstrip2 >= mstrip2min && mstrip2 <= mstrip2max ){ + // if ( mstrip1 >= mstrip1min && mstrip1 <= mstrip1max && mstrip2 >= mstrip2min && mstrip2 <= mstrip2max){ + // (*fmatrix[rbi])[mi][mj] += (mipv[nplane1][mstrip1] - cf->GetFullAverageAt(nplane1,mstrip1,erig)) * (mipv[nplane2][mstrip2] - cf->GetFullAverageAt(nplane2,mstrip2,erig)); + // (*fnmat[rbi])[mi][mj] += 1.; + (*fmatrix[rbi])[mi][mj] += (mip1 * mip2); // giusto + // (*fmatrix)[mi][mj] += (mip1 * mip2) * 1000000.; + toto++; + // (*fmatrix)[mi][mj] += 1.; + // cf->GetFullAverageAt(nplane1,mstrip1,erig,rbi); + // cf->GetFullAverageAt(nplane2,mstrip2,erig,rbi); + // (*fnmat)[mi][mj] += 1.; + // }; + }; + }; + }; + }; + }; + }; + // + printf(" toto = %i \n",toto); + printf("\n done \n"); + // printf(" write matrix \n"); + // cf->WriteFullMatrix(fmatrix, rbi); + // cf->WriteFullNMatrix(fnmat, rbi); + // printf(" done \n"); + // printf(" unload matrix \n"); + // cf->UnLoadFullMatrix(rbi); + // cf->UnLoadFullNMatrix(rbi); + // printf(" done \n"); + // printf(" delete matrix \n"); + // delete fmatrix; + // delete fnmat; + // printf(" done \n"); }; -// printf("4qui \n"); - // gObjectTable->Print(); } //=============================================================== @@ -323,37 +855,308 @@ // printf("Finished, calculating average and inverting matrices\n"); // - for (Int_t i=0; i 0. ){ - (*matrix[i])[ii][j] /= (*nmat[i])[ii][j]; - } else { - (*matrix[i])[ii][j] = 0.; + if ( !FULL ){ + for (Int_t i=0; i 0. ){ + (*matrix[i])[ii][j] /= (*nmat[i])[ii][j]; + } else { + (*matrix[i])[ii][j] = 0.; + }; }; }; + // + cf->WriteLongMatrix(matrix[i],i); + // + if ( matrix[i]->Determinant() == 0. ){ + printf("\n"); + for (Int_t ii=0; ii<43; ii++){ + for (Int_t j=0; j<43; j++){ + printf(" %.f",(*matrix[i])[ii][j]); + }; + printf("\n"); + }; + printf("\n"); + printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i); + } else { + Double_t det = 0.; + TMatrixD invmatrix = (TMatrixD)(matrix[i]->Invert(&det)); + printf(" Bin %i determinant is %f \n",i,det); + cf->WriteInvertedLongMatrix((TMatrixD)invmatrix,i); + }; }; + } else { // - cf->WriteLongMatrix(matrix[i],i); + // FULL // - if ( matrix[i]->Determinant() == 0. ){ - printf("\n"); + for (Int_t i=0; iLoadFullMatrix(i); + // fnmat = cf->LoadFullNMatrix(i); + // + // for (Int_t ii=0; ii 0. ){ + // // (*fmatrix[i])[ii][j] /= (*fnmat[i])[ii][j]; + // // } else { + // // (*fmatrix[i])[ii][j] = 0.; + // // }; + // if ( (*fnmat)[ii][j] > 0. ){ + // (*fmatrix)[ii][j] /= (*fnmat)[ii][j]; + // } else { + // (*fmatrix)[ii][j] = 0.; + // }; + // }; + // }; + // +// TMatrixD *mymat3 = new TMatrixD(129,129); +// TMatrixD *mymat5 = new TMatrixD(215,215); +// TMatrixD *mymat7 = new TMatrixD(301,301); +// TMatrixD *mymat9 = new TMatrixD(387,387); +// TMatrixD *mymat11 = new TMatrixD(473,473); +// TMatrixD *mymat17 = new TMatrixD(731,731); + // TMatrixF *mymat = new TMatrixF(129,129); + // TMatrixF *mymat = new TMatrixF(989,989); + Int_t i1 = -1; + Int_t j1 = -1; + // int mi,mj; + Int_t nonzero = 0; + Int_t nonzero1 = 0; for (Int_t ii=0; ii<43; ii++){ - for (Int_t j=0; j<43; j++){ - printf(" %.f",(*matrix[i])[ii][j]); + // for (Int_t j=0; j<191; j++){ + // for (Int_t j=0; j<43; j++){ + // for (Int_t j=0; j<31; j++){ + for (Int_t j=0; j<11; j++){ + // if ( (*fnmat[i])[ii][j] > 0. ){ + // (*fmatrix[i])[ii][j] /= (*fnmat[i])[ii][j]; + // } else { + // (*fmatrix[i])[ii][j] = 0.; + // }; + // i1 = (ii * 191) + j; + // i1 = (ii * 43) + j; + // i1 = (ii * 31) + j; + i1 = (ii * 11) + j; + // j1 = -1; + for (Int_t iij=0; iij<43; iij++){ + // for (Int_t jj=0; jj<191; jj++){ + // for (Int_t jj=0; jj<43; jj++){ + // for (Int_t jj=0; jj<31; jj++){ + for (Int_t jj=0; jj<11; jj++){ + // + // j1 = (iij * 191) + jj; + // j1 = (iij * 43) + jj; +// Int_t sh = -15 + ii; +// if ( sh > 15 ) sh -= 31*ii; +// // +// j1 = (iij * 31) + jj + sh; +// // +// if ( j1 < 0 ) j1 += 1333; +// if ( j1 >= 1333 ) j1 -= 1333; + // + // + // +// Int_t sh = -5 + ii; +// if ( sh > 5 ) sh -= 11*ii; +// // +// j1 = (iij * 11) + jj + sh; +// // +// if ( j1 < 0 ) j1 += 473; +// if ( j1 >= 473 ) j1 -= 473; +// + j1 = (iij * 11) + jj; + // j1++; + // if ( finmat[ii][j] > 0 ){ + // (*fmatrix)[i1][j1] /= finmat[ii][j]; + if ( (*fnmat[i])[ii][j] == 0. || (*fmatrix[i])[i1][j1] == 0. || !((*fmatrix[i])[i1][j1] == (*fmatrix[i])[i1][j1]) ){ + (*fmatrix[i])[i1][j1] = 1.; + } else { + (*fmatrix[i])[i1][j1] /= (*fnmat[i])[ii][j]; + nonzero++; + if ( i1 == 0 ) nonzero1++; + }; + // +// if ( j>=7 && j <=23 && jj >=7 && jj<=23 ){ +// Int_t mi17 = (ii*3) + j -7; +// Int_t mj17 = (iij*3) + jj -7; +// (*mymat17)[mi17][mj17] = (*fmatrix[i])[i1][j1]; +// }; +// if ( j>=10 && j <=20 && jj >=10 && jj<=20 ){ +// Int_t mi11 = (ii*3) + j -10; +// Int_t mj11 = (iij*3) + jj -10; +// (*mymat11)[mi11][mj11] = (*fmatrix[i])[i1][j1]; +// }; +// if ( j>=11 && j <=19 && jj >=11 && jj<=19 ){ +// Int_t mi9 = (ii*3) + j -11; +// Int_t mj9 = (iij*3) + jj -11; +// (*mymat9)[mi9][mj9] = (*fmatrix[i])[i1][j1]; +// }; +// if ( j>=12 && j <=18 && jj >=12 && jj<=18 ){ +// Int_t mi7 = (ii*3) + j -12; +// Int_t mj7 = (iij*3) + jj -12; +// (*mymat7)[mi7][mj7] = (*fmatrix[i])[i1][j1]; +// }; +// if ( j>=13 && j <=17 && jj >=13 && jj<=17 ){ +// Int_t mi5 = (ii*3) + j -13; +// Int_t mj5 = (iij*3) + jj -13; +// (*mymat5)[mi5][mj5] = (*fmatrix[i])[i1][j1]; +// }; +// if ( j>=14 && j <=16 && jj >=14 && jj<=16 ){ +// Int_t mi3 = (ii*3) + j -14; +// Int_t mj3 = (iij*3) + jj -14; +// (*mymat3)[mi3][mj3] = (*fmatrix[i])[i1][j1]; +// }; + + + // if ( j>=94 && j <=96 && jj >=94 && jj<=96 ){ + // mi = (ii*3) + j -94; + // mj = (iij*3) + jj -94; + // (*mymat)[mi][mj] = (*fmatrix)[i1][j1]; + // }; + + + // if ( j>=84 && j <=106 && jj >=84 && jj<=106 ){ + // mi = (ii*3) + j -84; + // mj = (iij*3) + jj -84; + // (*mymat)[mi][mj] = (*fmatrix)[i1][j1]; + // }; + + }; + }; }; - printf("\n"); }; - printf("\n"); - printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i); - } else { - Double_t det = 0.; - TMatrixD invmatrix = (TMatrixD)(matrix[i]->Invert(&det)); - printf(" Bin %i determinant is %f \n",i,det); - cf->WriteInvertedLongMatrix((TMatrixD)invmatrix,i); + // + printf(" Matrix has %i non-zero elements \n",nonzero); + // printf(" Matrix has %i non-zero elements on the first row\n",nonzero1); + // + // Bool_t BAD = false; + // for (Int_t ii=0; ii<43; ii++){ + // for (Int_t j=0; j<191; j++){ + // // + // i1 = (ii * 191) + j; + // // + // for (Int_t iij=0; iij<43; iij++){ + // for (Int_t jj=0; jj<191; jj++){ + // // + // j1 = (iij * 191) + jj; + // // + // // printf(" ROW %i COLUMN %i VALUE %f \n",i1,j1,(*fmatrix)[i1][j1]); + // if ( (*fmatrix)[i1][j1] == 0. || !((*fmatrix)[i1][j1]==(*fmatrix)[i1][j1]) ){ + // printf(" ROW %i COLUMN %i VALUE %f \n",i1,j1,(*fmatrix)[i1][j1]); + // printf(" che schifo! \n"); + // BAD = true; + // }; + // // + // }; + // }; + // }; + // }; + // // + // if ( BAD ) printf(" questa matrice fa cagare \n"); + // + // + // cf->WriteFullMatrix(fmatrix, i); + // cf->WriteFullNMatrix(fnmat, i); + //cf->WriteFullMatrix(fmatrix[i],i);// OK + // cf->WriteFullNMatrix(fnmat[i], i); //OK + // + // TDecompSVD svd(*fmatrix[i]); + // Bool_t ok = svd.Decompose(); + // + Double_t zero = (Double_t)0.0; + // + if ( fmatrix[i]->Determinant() == zero ){ + //if ( fmatrix->Determinant() == 0. ){ + printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i); + cf->WriteFullMatrix(fmatrix[i],i);// OK + cf->WriteFullNMatrix(fnmat[i], i); //OK + } else { + // }; + // if ( i == 3 ){ + // if ( ok ){ + // Double_t tol = 1E-20; + // TDecompSVD svd((*fmatrix)[i],tol); + // svd.Decompose(); + // TMatrixD svdInv = svd.Invert(); + // svdInv.Print("svdInv"); + // cout << "condition: " << svd.Condition() << endl; + // cf->WriteInvertedFullMatrix((TMatrixD)svdInv,999); + + Double_t det = 0.; + TMatrixD invmatrix = (TMatrixD)(fmatrix[i]->Invert(&det)); + printf(" Bin %i determinant is %f \n",i,det); + cf->WriteInvertedFullMatrix((TMatrixD)invmatrix,i); + }; + +// if ( mymat3->Determinant() == 0. ){ +// printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i); +// } else { +// Double_t det = 0.; +// TMatrixD invmatrix = (TMatrixD)(mymat3->Invert(&det)); +// printf(" Mymat3 determinant is %f \n",det); +// cf->WriteInvertedFullMatrix((TMatrixD)invmatrix,1103); +// }; +// cf->WriteFullMatrix(mymat3, 103); +// if ( mymat5->Determinant() == 0. ){ +// printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i); +// } else { +// Double_t det = 0.; +// TMatrixD invmatrix = (TMatrixD)(mymat5->Invert(&det)); +// printf(" Mymat5 determinant is %f \n",det); +// cf->WriteInvertedFullMatrix((TMatrixD)invmatrix,1105); +// }; +// cf->WriteFullMatrix(mymat5, 105); +// if ( mymat7->Determinant() == 0. ){ +// printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i); +// } else { +// Double_t det = 0.; +// TMatrixD invmatrix = (TMatrixD)(mymat7->Invert(&det)); +// printf(" Mymat7 determinant is %f \n",det); +// cf->WriteInvertedFullMatrix((TMatrixD)invmatrix,1107); +// }; +// cf->WriteFullMatrix(mymat7, 107); +// if ( mymat9->Determinant() == 0. ){ +// printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i); +// } else { +// Double_t det = 0.; +// TMatrixD invmatrix = (TMatrixD)(mymat9->Invert(&det)); +// printf(" Mymat3 determinant is %f \n",det); +// cf->WriteInvertedFullMatrix((TMatrixD)invmatrix,1109); +// }; +// cf->WriteFullMatrix(mymat9, 109); +// if ( mymat11->Determinant() == 0. ){ +// printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i); +// } else { +// Double_t det = 0.; +// TMatrixD invmatrix = (TMatrixD)(mymat11->Invert(&det)); +// printf(" Mymat11 determinant is %f \n",det); +// cf->WriteInvertedFullMatrix((TMatrixD)invmatrix,1111); +// }; +// cf->WriteFullMatrix(mymat11, 111); +// if ( mymat17->Determinant() == 0. ){ +// printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i); +// } else { +// Double_t det = 0.; +// TMatrixD invmatrix = (TMatrixD)(mymat17->Invert(&det)); +// printf(" Mymat3 determinant is %f \n",det); +// cf->WriteInvertedFullMatrix((TMatrixD)invmatrix,1117); +// }; +// cf->WriteFullMatrix(mymat17, 117); + + + // + // cf->UnLoadFullMatrix(i); + // cf->UnLoadFullNMatrix(i); + // delete fmatrix; + // delete fnmat; + // }; }; //