--- calo/flight/CaloFranzini/src/Calib.cpp 2007/12/13 17:08:09 1.1 +++ calo/flight/CaloFranzini/src/Calib.cpp 2008/01/11 15:27:13 1.4 @@ -6,6 +6,7 @@ #include #include #include +#include #include #include // @@ -14,17 +15,38 @@ // 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[18]; -Float_t qqplane[18][43]; -Int_t nnqplane[18][43]; - -TArrayF *qplane[18]; -TArrayI *nqplane[18]; -TMatrixD *matrix[18]; -TMatrixD *nmat[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]; + +TArrayF *qplane[17]; +TArrayI *nqplane[17]; +TMatrixD *matrix[17]; +TMatrixD *nmat[17]; + +TMatrixD *fqplane; +TMatrixD *fnqplane; +//TMatrixD *fqplane[17]; +//TMatrixD *fnqplane[17]; +//TMatrixF *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 ){ @@ -40,7 +62,9 @@ // 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 && @@ -66,20 +90,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() && @@ -87,6 +110,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; @@ -104,6 +177,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 @@ -115,11 +197,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->GetCaloLevel1()->qtotpl(0) > 7. ) 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 + }; + // + 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; } @@ -134,7 +267,19 @@ void CreateHistos( PamLevel2* event , TString file){ cf = new CaloFranzini(event); - cf->CreateMatrixFile(file.Data()); + // + if ( MATRIX ){ + cf->UpdateMatrixFile(file.Data()); + cf->LoadBin(); + if ( !FULL ){ + cf->LoadLong(); + } else { + cf->LoadFull(); + }; + } else { + cf->CreateMatrixFile(file.Data()); + }; + // // nbin = 18; rig[0] = 0.1; @@ -156,105 +301,236 @@ rig[16] = 200.; rig[17] = 400.; // -// memset(qplane, 0, 44*18*sizeof(Float_t)); - memset(rmean, 0, 18*sizeof(Float_t)); - memset(qqplane, 0, 43*18*sizeof(Float_t)); - memset(nnqplane, 0, 43*18*sizeof(Float_t)); - // - for (Int_t i=0; i < 18 ; 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=0; i < 18 ; i++){ -// for (Int_t j = 0; j < 43;i++){ -// (*qplane[i])[j] = 0.; -// qplane[i]->Reset(); -// nqplane[i]->Reset(); -// for (Int_t k = 0; k < 43;k++){ -// (*matrix[i])[j][k] = 0.; -// (*nmat[i])[j][k] = 0.; -// }; -// }; -// }; + memset(rmean, 0, 17*sizeof(Float_t)); + memset(ntot, 0, 17*sizeof(Int_t)); + // memset(finmat, 0, 43*191*sizeof(Int_t)); + // + for (Int_t i=0; i < 17 ; i++){ + //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); + fnmat[i] = new TMatrixF(43,191); +// 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);// + // + cf->WriteFullMean(fqplane, i); + cf->WriteFullNMean(fnqplane, i); + delete fqplane; + delete fnqplane; + // + }; + }; + }; // } //=============================================================== void FindAverage( PamLevel2* L2, int iev ){ // - // printf("SELEZIONATO \n"); 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[i] && erig < rig[i+1] ){ rbi = i; break; }; }; // - // printf("SELEZIONATO rbi %i erig %f rig[0] %f rig[nbin] %f \n",rbi, erig, rig[0],rig[nbin-1]); - // if ( erig < rig[0] ) return; if ( erig >= rig[nbin-1] ) return; // rmean[rbi] += erig; + ntot[rbi]++; // - Int_t dgf = 43; - // -// for (Int_t i=0; i < 22; i++){ -// if ( L2->GetTrack(0)->GetCaloTrack()->tibar[i][1] < 0 ){ -// dgf = 2 * i; - // break; - // }; - // if ( L2->GetTrack(0)->GetCaloTrack()->tibar[i][0] < 0 ){ - // dgf = 1 + 2 * i; - // break; - // }; - // }; - // -// if ( dgf < 43 && dgf > 37 ) dgf--; - // - 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 + // + fqplane = cf->LoadFullAverage(rbi); + fnqplane = cf->LoadFullNAverage(rbi); + 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 cs = 0; + Int_t cd = 0; + Int_t mstrip = 0; + // + for (Int_t j=0; j<2; j++){ + for (Int_t i=0; i<21; i++){ + nplane = 1 - j + 2*i; + if ( nplane > 37 ) nplane--; + // + cs = ct->tibar[i][j] - 1; + // + cd = 95 - cs; + // + 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.; + }; + }; + }; // - mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip); // - nplane = 1 - view + 2 * plane; - if ( nplane > 37 ) nplane--; - // printf(" mip %f view %i plane %i => nplane %i \n",mip,view,plane,nplane); - if ( nplane < dgf ){ - // printf(" PRIMA: qpl %f nqpl %i \n",(*qplane[rbi])[nplane],(*nqplane[rbi])[nplane]); - (*qplane[rbi])[nplane] += mip; - qqplane[rbi][nplane] += mip; - // printf(" DOPO: qpl %f nqpl %i \n",(*qplane[rbi])[nplane],(*nqplane[rbi])[nplane]); + 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; + // + // (*fqplane[rbi])[nplane][mstrip] += mip; + (*fqplane)[nplane][mstrip] += mip; + // }; // + cf->WriteFullMean(fqplane, rbi); + cf->WriteFullNMean(fnqplane, rbi); + cf->UnLoadFullAverage(rbi); + cf->UnLoadFullNAverage(rbi); + delete fqplane; + delete fnqplane; + // }; } +void CalculateAverage(){ + // + 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++){ +// 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)[j][k] > 0 ){ + if ( (*fqplane)[j][k] == 0. ) (*fqplane)[j][k] = 0.7; + (*fqplane)[j][k] /= (Float_t)(*fnqplane)[j][k]; + } else { + (*fqplane)[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]); + }; + }; + cf->WriteFullMean(fqplane, i); + cf->WriteFullNMean(fnqplane, 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"); + // + // +} + //=============================================================== void FindMatrix( PamLevel2* L2, int iev ){ // - // printf("SELEZIONATO \n"); 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; - // -// for (Int_t i=0; i < 22; i++){ -// if ( L2->GetTrack(0)->GetCaloTrack()->tibar[i][1] < 0 ){ -// dgf = 2 * i; - // break; - // }; - // if ( L2->GetTrack(0)->GetCaloTrack()->tibar[i][0] < 0 ){ - // dgf = 1 + 2 * i; - // break; - // }; -// }; - // -// if ( dgf < 43 && dgf > 37 ) dgf--; - // - for (Int_t i=0; iGetCaloLevel1()->istrip; i++ ){ // - mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip); + // Fill the estrip matrix // - nplane = 1 - view + 2 * plane; - if ( nplane > 37 ) nplane--; - if ( nplane < dgf ){ - hpl[nplane] += mip; + 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; + }; + // }; // - }; - // - for (Int_t i=0; iGetAverageAt(i,erig)) * (hpl[j] - cf->GetAverageAt(j,erig)); + }; }; - }; -} - -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.; + } else { + // + // FULL CALORIMETER + // + printf(" Retrieve 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; + Int_t dgf = 191; + // Int_t mindgf = 85; + // Int_t dgf = 115; + Int_t cs = 0; + Int_t cd = 0; + Int_t mstrip = 0; + // + Float_t mipv[43][191]; + memset(mipv,0,43*191*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; + // + mipv[nplane][mstrip] = 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; + // + 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; + // + mstrip1min = TMath::Max(mindgf,(cd1+0)); + mstrip1max = TMath::Min(dgf,(cd1+95)) + 1; + // + if ( nplane1 == 0 || nplane1 == 42 ) printf(" pl %i mstrip1min %i mstrip1max %i \n",nplane1,mstrip1min,mstrip1max); + // + for (Int_t mstrip1=mstrip1min; mstrip1GetFullAverageAt(nplane1,mstrip1,erig,rbi); + // + mi = (nplane1 * 191) + 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; + mstrip2min = TMath::Max(mindgf,(cd2+0)); + mstrip2max = TMath::Min(dgf,(cd2+95)) + 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,rbi); + // + if ( mip2 != 0. ){ + // + mj = (nplane2 * 191) + mstrip2; + // 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)[mi][mj] += (mip1 * mip2); // giusto + toto++; + // (*fmatrix)[mi][mj] += 1.; + // cf->GetFullAverageAt(nplane1,mstrip1,erig,rbi); + // cf->GetFullAverageAt(nplane2,mstrip2,erig,rbi); + // (*fnmat)[mi][mj] += 1.; + // }; + }; + }; + }; + }; }; - printf(" BIN %i plane %i average energy %f qplane %f nqplane %i \n",i,j,rmean[i],(*qplane[i])[j],(*nqplane[i])[j]); -// if ( nnqplane[i][j] > 0 ){ -// qqplane[i][j] /= (Float_t)nnqplane[i][j]; -// } else { -// qqplane[i][j] = 0.; -// }; -// printf(" BBIN %i plane %i average energy %f qplane %f nqplane %i \n",i,j,rmean[i],qqplane[i][j],nnqplane[i][j]); }; - }; - - cf->WriteNumBin(nbin); - - TArrayF *rigbin = new TArrayF(18, rig); - cf->WriteRigBin(rigbin); - - TArrayF *rmeanbin = new TArrayF(18, rmean); - TFile *file = cf->GetFile(); - file->cd(); - // rigbin->Write("binrig"); - file->WriteObject(&(*rmeanbin), "binrigmean"); - -// cf->WriteRigBin(rigbin); - - for (Int_t i=0; iWriteLongMean(qplane[i], i); - + // + 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"); }; } + //=============================================================== // Save histograms // @@ -370,70 +755,166 @@ // //=============================================================== void SaveHistos(){ - - printf("Finished, calculating average and inverting matrices\n"); - - - -// cf->WriteNumBin(nbin); - -// TArrayF *rigbin = new TArrayF(18, rig); -// cf->WriteRigBin(rigbin); - - for (Int_t i=0; iWriteLongMean(qplane[i], i); - + // + if ( MATRIX ){ + // + printf("Finished, calculating average and inverting matrices\n"); // - // determine the average matrix - // - for (Int_t ii=0; ii<43; ii++){ - for (Int_t j=0; j<43; j++){ - // if ( i == 1 ) printf(" num %f \n",(*nmat[i])[ii][j]); - if ( (*nmat[i])[ii][j] > 0. ){ - // if ( i == 1 ) printf(" value %f \n",(*matrix[i])[ii][j]); - (*matrix[i])[ii][j] /= (*nmat[i])[ii][j]; + 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 { - (*matrix[i])[ii][j] = 0.; + 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); }; }; - }; - -// if ( i == 2 ){ -// printf("\n"); -// for (Int_t ii=0; ii<43; ii++){ -// for (Int_t j=0; j<43; j++){ -// printf(" %.f",(*matrix[i])[ii][j]); + } else { + // + // FULL + // + // 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.; +// }; +// }; // }; -// printf("\n"); -// }; -// printf("\n"); -// }; - - - 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"); +// + Int_t i1 = -1; + Int_t j1 = -1; + Int_t nonzero = 0; + Int_t nonzero1 = 0; + for (Int_t ii=0; ii<43; ii++){ + for (Int_t j=0; j<191; 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; + // j1 = -1; + for (Int_t iij=0; iij<43; iij++){ + for (Int_t jj=0; jj<191; jj++){ + // + j1 = (iij * 191) + jj; + // j1++; + // if ( finmat[ii][j] > 0 ){ + // (*fmatrix)[i1][j1] /= finmat[ii][j]; + if ( (*fnmat[i])[ii][j] == 0. || (*fmatrix)[i1][j1] == 0. || !((*fmatrix)[i1][j1] == (*fmatrix)[i1][j1]) ){ + (*fmatrix)[i1][j1] = 0.; + } else { + (*fmatrix)[i1][j1] /= (*fnmat[i])[ii][j]; + nonzero++; + if ( i1 == 0 ) nonzero1++; + }; + }; + }; + }; + }; + // + 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],i); + cf->WriteFullMatrix(fmatrix, i); + // cf->WriteFullNMatrix(fnmat, i); + cf->WriteFullNMatrix(fnmat[i], i); + // + // if ( fmatrix[i]->Determinant() == 0. ){ + if ( fmatrix->Determinant() == 0. ){ + printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i); + } else { + Double_t det = 0.; + // TMatrixF invmatrix = (TMatrixF)(fmatrix[i]->Invert(&det)); + // printf(" Bin %i determinant is %f \n",i,det); + // cf->WriteInvertedFullMatrix((TMatrixF)invmatrix,i); + TMatrixF invmatrix = (TMatrixF)(fmatrix->Invert(&det)); + printf(" Bin %i determinant is %f \n",i,det); + cf->WriteInvertedFullMatrix((TMatrixF)invmatrix,i); + }; + // + cf->UnLoadFullMatrix(i); + // cf->UnLoadFullNMatrix(i); + delete fmatrix; + // delete fnmat; + // }; - 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(" done, closing file and exiting\n"); + // }; - - printf(" done, closing file and exiting\n"); - + // cf->CloseMatrixFile(); - + // cf->Delete(); + // }