/** * \file TrkLevel1.cpp * \author Elena Vannuccini */ #include #include using namespace std; //...................................... // F77 routines //...................................... extern "C" { // int readetaparam_(); float cog_(int*,int*); float pfaeta_(int*,float*); float pfaeta2_(int*,float*); float pfaeta3_(int*,float*); float pfaeta4_(int*,float*); float pfaetal_(int*,float*); float digsat_(int*); int npfastrips_(int*,float*); } //-------------------------------------- // // //-------------------------------------- TrkCluster::TrkCluster(){ // cout << "TrkCluster::TrkCluster()"< cut*sigma ). * @param force Falg to force the PFA strip-inclusion pattern (nstrip>0) * If nstrip<=0 only the inclusion cut is used to determine the cluster size. */ Float_t TrkCluster::GetSignal(Int_t nstrip, Float_t cut, Bool_t force){ if(CLlength<=0)return 0; Float_t s = 0; //----------------------------------- // inlcude strips with s > cut*sigma //----------------------------------- if( nstrip<=0 ){ // for(Int_t is = 0; is < CLlength; is++){ // Float_t scut = cut*clsigma[is]; // if(clsignal[is] > scut) s += clsignal[is]; // }; for(Int_t is = indmax+1; is < CLlength; is++){ Float_t scut = cut*clsigma[is]; if(clsignal[is] > scut) s += clsignal[is]; else break; }; for(Int_t is = indmax; is >=0; is--){ Float_t scut = cut*clsigma[is]; if(clsignal[is] > scut) s += clsignal[is]; else break; }; return s; }; //--------------------------------------------------- // inlcude strips with s > cut*sigma, up to nstrip. // strips are included in order of decreasing signal //--------------------------------------------------- if( !force ){ Int_t il = indmax; Int_t ir = indmax; Int_t inc = 0; if( clsignal[indmax] < cut*clsigma[indmax] ) return 0; while ( inc < nstrip ){ Float_t sl = -100000; Float_t sr = -100000; if( il >= 0 ) sl = clsignal[il]; if( ir < CLlength ) sr = clsignal[ir]; if( sl == sr && inc == 0 ){ s += clsignal[il]; //cout << inc<<" - "<< clsignal[il]<<" "<= sr && sl > cut*clsigma[il] && inc !=0 ){ s += sl;//cout << inc<<" - "<< clsignal[il]<<" "< cut*clsigma[ir] ){ s += sr;//cout << inc<<" - " << clsignal[ir]<<" "< signal of the central strip Float_t sc = clsignal[indmax]; // signal of adjacent strips Float_t sl1 = -9999.; Float_t sl2 = -9999.; Float_t sr1 = -9999.; Float_t sr2 = -9999.; if(indmax-1>=0) sl1 = clsignal[indmax-1]; if(indmax-2>=0) sl2 = clsignal[indmax-2]; if(indmax+1sr1 && sl1+sc!=0 )s = (sl1+sc); if( sl1 clsigma[indmax+1] && sc+sr1!=0 )s = (sc+sr1); } }else if(nstrip==3){ s = (sl1+sc+sr1); }else if(nstrip==4){ if( sl2>sr2 && sl2+sl1+sc+sr1!=0 )s = (sl2+sl1+sc+sr1); if( sl2 clsigma[indmax+2] && sl1+sc+sr1+sr2!=0 )s = (sl1+sc+sr1+sr2); } }else if(nstrip==5){ s = (sl1+sc+sr1); if(sl2 != -9999.)s += sl2; if(sr2 != -9999.)s += sr2; }else{ cout << "Float_t TrkCluster::GetSignal("< cut*sigma ). * If nstrip<=0 only the inclusion cut is used to determine the cluster size. */ Float_t TrkCluster::GetSignalToNoise(Int_t nstrip, Float_t cut){ if(CLlength<=0)return 0; Float_t sn = 0; if( nstrip<=0 ){ for(Int_t is = indmax+1; is < CLlength; is++){ Float_t scut = cut*clsigma[is]; if(clsignal[is] > scut) sn += clsignal[is]/clsigma[is]; else break; }; for(Int_t is = indmax; is >=0; is--){ Float_t scut = cut*clsigma[is]; if(clsignal[is] > scut) sn += clsignal[is]/clsigma[is]; else break; }; return sn; }; Int_t il = indmax; Int_t ir = indmax; Int_t inc = 0; if( clsignal[indmax] < cut*clsigma[indmax] ) return 0; while ( inc < nstrip ){ Float_t sl = -100000; Float_t sr = -100000; if( il >= 0 ) sl = clsignal[il]; if( ir < CLlength ) sr = clsignal[ir]; if( sl == sr && inc == 0 ){ sn += clsignal[il]/clsigma[il]; il--; ir++; }else if ( sl >= sr && sl > cut*clsigma[il] && inc !=0 ){ sn += sl/clsigma[il]; il--; }else if ( sl < sr && sr > cut*clsigma[ir] ){ sn += sr/clsigma[ir]; ir++; }else break; inc++; } return sn; }; /** * Evaluate the cluster multiplicity. * @param cut Inclusion cut. */ Int_t TrkCluster::GetMultiplicity(Float_t cut){ if(CLlength<=0)return 0; Int_t m = 0; for(Int_t is = indmax+1; is < CLlength; is++){ Float_t scut = cut*clsigma[is]; if(clsignal[is] > scut) m++; else break; }; for(Int_t is = indmax; is >=0; is--){ Float_t scut = cut*clsigma[is]; if(clsignal[is] > scut) m++; else break; }; return m; }; /** * True if the cluster contains bad strips. * @param nbad Number of strips around the maximum. */ Bool_t TrkCluster::IsBad(Int_t nbad){ if(CLlength<=0)return 0; Int_t il,ir; il = indmax; ir = indmax; for(Int_t i=1; i clsignal[ir+1])il--; else ir++; } } Int_t isbad = 0; for(Int_t i=il; i<=ir; i++)isbad += clbad[i]; return ( isbad != nbad ); }; /** * True if the cluster contains saturated strips. * @param nbad Number of strips around the maximum. */ Bool_t TrkCluster::IsSaturated(Int_t nbad){ if(CLlength<=0)return 0; Int_t il,ir; il = indmax; ir = indmax; for(Int_t i=1; i clsignal[ir+1])il--; else ir++; } } Int_t isbad = 0; for(Int_t i=il; i<=ir; i++){ if( IsX() && cladc[i] > 2980 )isbad++; if( IsY() && cladc[i] < 80 )isbad++; } return ( isbad != 0 ); } //-------------------------------------- // // //-------------------------------------- void TrkCluster::Dump(){ cout << "----- Cluster" << endl; cout << "View "<nclstr1 = 1; l1->view[0] = view; l1->ladder[0] = GetLadder(); l1->maxs[0] = maxs; l1->mult[0] = GetMultiplicity(); l1->dedx[0] = GetSignal(); l1->indstart[0] = 1; l1->indmax[0] = indmax+1; l1->totCLlength = CLlength; for(Int_t i=0; iclsignal[i] = clsignal[i]; l1->clsigma[i] = clsigma[i]; l1->cladc[i] = cladc[i]; l1->clbad[i] = clbad[i]; }; // return l1; }; //-------------------------------------- // // //-------------------------------------- /** * Evaluates the Center-Of-Gravity (COG) of the cluster, in strips, relative to the strip with the maximum signal (TrkCluster::maxs). * @param ncog Number of strips to evaluate COG. * If ncog=0, the COG of the cluster is evaluated according to the cluster multiplicity (defined by the inclusion cut). * If ncog>0, the COG is evaluated using ncog strips, even if they have a negative signal (according to G.Landi) * * (NB TrkCluster::GetLevel1Struct() showld be called first, in order to fill the F77 level1 common with this single cluster) */ Float_t TrkCluster::GetCOG(Int_t ncog){ int ic = 1; // GetLevel1Struct(); //Elena: dangerous... return cog_(&ncog,&ic); }; /** * Evaluates the Center-Of-Gravity (COG) of the cluster, in strips, relative to the strip with the maximum signal (TrkCluster::maxs), * choosing the number of strips according to the angle, as implemented for the eta-algorythm . * @param angle Projected angle in degree. */ Float_t TrkCluster::GetCOG(Float_t angle){ Int_t neta = 0; // Float_t eta = GetETA(0,angle); // for(neta=2; neta<10; neta++) if( eta == GetETA(neta,angle) ) break; // if(eta != GetETA(neta,angle) )cout << "Attenzione!! pasticcio "< 10. && fabs(angle) <= 15. ){ neta = 3; }else{ neta = 4; }; }; return GetCOG(neta); }; //-------------------------------------- // // //-------------------------------------- /** * Evaluates the cluster position, in pitch units, relative to the strip * with the maximum signal (TrkCluster::maxs), by applying the non-linear * ETA-algorythm. * @param neta Number of strips to evaluate ETA. * @param angle Projected (effective) angle between particle track and detector plane. * @landi flag to apply Landi correction * Implemented values of neta are 2,3,4. If neta=0, ETA2, ETA3 and ETA4 are applied according to the angle. * (NB TrkCluster::GetLevel1Struct() showld be called first, in order to fill the F77 level1 common with this single cluster) */ Float_t TrkCluster::GetETA(Int_t neta, float angle, bool landi){ // cout << "GetETA(neta,angle) "<< neta << " "<< angle; // LoadPfaParam(); TrkParams::Load(4); if( !TrkParams::IsLoaded(4) ){ cout << "Float_t TrkCluster::GetETA(Int_t neta, float angle, bool landi) --- ERROR --- p.f.a. parameters not loaded"<0)mask=1; cout<nclstr(); i++) ((TrkCluster *)t[i])->Dump(); } /** * \brief Dump processing status */ void TrkLevel1::StatusDump(int view){ cout << "DSP n. "<= 12)return false; return !(good[view]&flagmask); }; //-------------------------------------- // // //-------------------------------------- /** * Fills a TrkLevel1 object with values from a struct cTrkLevel1 (to get data from F77 common). */ void TrkLevel1::SetFromLevel1Struct(cTrkLevel1 *l1, Bool_t full){ // cout << "void TrkLevel1::SetFromLevel1Struct(cTrkLevel1 *l1, Bool_t full)"<nclstr1; i++){ t_cl->Clear(); // if( full || (!full && l1->whichtrack[i]) ){ t_cl->view = l1->view[i]; t_cl->maxs = l1->maxs[i]; if( full || (!full && l1->whichtrack[i]) ){ t_cl->indmax = l1->indmax[i] - l1->indstart[i]; Int_t from = l1->indstart[i] -1; Int_t to = l1->totCLlength ; if(i != l1->nclstr1-1)to = l1->indstart[i+1] -1 ; t_cl->CLlength = to - from ; t_cl->clsignal = new Float_t[t_cl->CLlength]; t_cl->clsigma = new Float_t[t_cl->CLlength]; t_cl->cladc = new Int_t[t_cl->CLlength]; t_cl->clbad = new Bool_t[t_cl->CLlength]; Int_t index = 0; for(Int_t is = from; is < to; is++ ){ t_cl->clsignal[index] = (Float_t) l1->clsignal[is]; t_cl->clsigma[index] = (Float_t) l1->clsigma[is]; t_cl->cladc[index] = (Int_t) l1->cladc[is]; t_cl->clbad[index] = (Bool_t) l1->clbad[is]; index++; }; } new(t[i]) TrkCluster(*t_cl); // <<< store cluster }; delete t_cl; // ------------------------- // ****general variables**** // ------------------------- for(Int_t i=0; i<12 ; i++){ good[i] = l1->good[i]; for(Int_t j=0; j<24 ; j++){ cn[j][i] = l1->cnev[j][i]; // cnrms[j][i] = l1->cnrmsev[j][i]; cnn[j][i] = l1->cnnev[j][i]; }; }; } /** * Fills a struct cTrkLevel1 with values from a TrkLevel1 object (to put data into a F77 common). */ void TrkLevel1::GetLevel1Struct(cTrkLevel1* l1) { // cTrkLevel1* l1 = &level1event_ ; for(Int_t i=0; i<12 ; i++){ l1->good[i] = good[i]; for(Int_t j=0; j<24 ; j++){ l1->cnev[j][i] = cn[j][i] ; l1->cnnev[j][i] = cnn[j][i] ; l1->cnrmsev[j][i] = 0. ; }; l1->fshower[i] = 0; }; l1->nclstr1=0; l1->totCLlength=0; Int_t index=0; if(Cluster){ Int_t i=0; for(Int_t ii=0;iiGetEntries();ii++){ TrkCluster *clu = GetCluster(ii); // ---------------------------------------- // attenzione!! // se il cluster non e` salvato (view = 0) // DEVE essere escluso dal common F77 // ---------------------------------------- if(clu->view != 0 ){ l1->view[i] = clu->view; l1->ladder[i] = clu->GetLadder(); l1->maxs[i] = clu->maxs; l1->mult[i] = clu->GetMultiplicity(); l1->dedx[i] = clu->GetSignal(); l1->indstart[i] = index+1; l1->indmax[i] = l1->indstart[i] + clu->indmax; l1->totCLlength += clu->CLlength; for(Int_t iw=0; iw < clu->CLlength; iw++){ l1->clsignal[index] = clu->clsignal[iw]; l1->clsigma[index] = clu->clsigma[iw]; l1->cladc[index] = clu->cladc[iw]; l1->clbad[index] = clu->clbad[iw]; index++; } i++; } } l1->nclstr1 = i; } // return l1; } //-------------------------------------- // // //-------------------------------------- void TrkLevel1::Clear(){ for(Int_t i=0; i<12 ; i++){ good[i] = -1; for(Int_t j=0; j<24 ; j++){ cn[j][i] = 0; cnn[j][i] = 0; }; }; // if(Cluster)Cluster->Clear("C"); if(Cluster)Cluster->Delete(); } //-------------------------------------- // // //-------------------------------------- void TrkLevel1::Delete(){ // Clear(); if(Cluster)Cluster->Delete(); if(Cluster)delete Cluster; } //-------------------------------------- // // //-------------------------------------- TrkCluster *TrkLevel1::GetCluster(int is){ if(!Cluster)return 0; if(is >= nclstr()){ cout << "** TrkLevel1::GetCluster(int) ** Cluster "<< is << " does not exits! " << endl; cout << "( Stored clusters nclstr() = "<< this->nclstr()<<" )" << endl; return 0; } TClonesArray &t = *(Cluster); TrkCluster *cluster = (TrkCluster*)t[is]; return cluster; } // int TrkLevel1::GetPfaNbinsAngle(){ // TrkParams::Load(4); // if( !TrkParams::IsLoaded(4) ){ // cout << "int TrkLevel1::GetPfaNbinsAngle() --- ERROR --- p.f.a. parameters not loaded"<level2 processing. * The level2 output is stored in a common block, which can be retrieved * by mean of the method TrkLevel2::SetFromLevel2Struct(). * NB If the TrkLevel1 object is readout from a tree, and the * TrkLevel1::ProcessEvent(int pfa) is used to reprocess the event, attention * should be payed to the fact that single clusters (clusters not associated * with any track) might not be stored. Full reprocessing should be done starting * from level0 data. */ //int TrkLevel1::ProcessEvent(int pfa){ int TrkLevel1::ProcessEvent(){ // cout << "int TrkLevel1::ProcessEvent()" << endl; TrkParams::Load( ); if( !TrkParams::IsLoaded() )return 0; GetLevel1Struct(); // analysisflight_(&pfa); // TrkParams::SetPFA(pfa); analysisflight_(); return 1; } ClassImp(TrkLevel1); ClassImp(TrkCluster);