/************************************************************************** * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * * * * Author: The ALICE Off-line Project. * * Contributors are mentioned in the code where appropriate. * * * * Permission to use, copy, modify and distribute this software and its * * documentation strictly for non-commercial purposes is hereby granted * * without fee, provided that the above copyright notice appears in all * * copies and that both the copyright notice and this permission notice * * appear in the supporting documentation. The authors make no claims * * about the suitability of this software for any purpose. It is * * provided "as is" without express or implied warranty. * **************************************************************************/ /* $Id$ $Log$ Revision 1.24 2002/04/24 22:02:31 nilsen New SDigits and Digits routines, and related changes, (including new noise values). */ #include #include #include #include #include "AliITSClusterFinderSDD.h" #include "AliITSMapA1.h" #include "AliITS.h" #include "AliITSdigit.h" #include "AliITSRawCluster.h" #include "AliITSRecPoint.h" #include "AliITSsegmentation.h" #include "AliITSresponseSDD.h" #include "AliRun.h" ClassImp(AliITSClusterFinderSDD) //______________________________________________________________________ AliITSClusterFinderSDD::AliITSClusterFinderSDD(AliITSsegmentation *seg, AliITSresponse *response, TClonesArray *digits, TClonesArray *recp){ // standard constructor fSegmentation = seg; fResponse = response; fDigits = digits; fClusters = recp; fNclusters = fClusters->GetEntriesFast(); SetCutAmplitude(); SetDAnode(); SetDTime(); SetMinPeak((Int_t)(((AliITSresponseSDD*)fResponse)->GetNoiseAfterElectronics()*5)); // SetMinPeak(); SetMinNCells(); SetMaxNCells(); SetTimeCorr(); SetMinCharge(); fMap = new AliITSMapA1(fSegmentation,fDigits,fCutAmplitude); } //______________________________________________________________________ AliITSClusterFinderSDD::AliITSClusterFinderSDD(){ // default constructor fSegmentation = 0; fResponse = 0; fDigits = 0; fClusters = 0; fNclusters = 0; fMap = 0; fCutAmplitude = 0; fDAnode = 0; fDTime = 0; fMinPeak = 0; fMinNCells = 0; fMaxNCells = 0; fTimeCorr = 0; fMinCharge = 0; /* SetDAnode(); SetDTime(); SetMinPeak((Int_t)(((AliITSresponseSDD*)fResponse)->GetNoiseAfterElectronics()*5)); SetMinNCells(); SetMaxNCells(); SetTimeCorr(); SetMinCharge(); */ } //____________________________________________________________________________ AliITSClusterFinderSDD::~AliITSClusterFinderSDD(){ // destructor if(fMap) delete fMap; } //______________________________________________________________________ void AliITSClusterFinderSDD::SetCutAmplitude(Float_t nsigma){ // set the signal threshold for cluster finder Float_t baseline,noise,noise_after_el; fResponse->GetNoiseParam(noise,baseline); noise_after_el = ((AliITSresponseSDD*)fResponse)->GetNoiseAfterElectronics(); fCutAmplitude = (Int_t)((baseline + nsigma*noise_after_el)); } //______________________________________________________________________ void AliITSClusterFinderSDD::Find1DClusters(){ // find 1D clusters static AliITS *iTS = (AliITS*)gAlice->GetModule("ITS"); // retrieve the parameters Int_t fNofMaps = fSegmentation->Npz(); Int_t fMaxNofSamples = fSegmentation->Npx(); Int_t fNofAnodes = fNofMaps/2; Int_t dummy = 0; Float_t fTimeStep = fSegmentation->Dpx(dummy); Float_t fSddLength = fSegmentation->Dx(); Float_t fDriftSpeed = fResponse->DriftSpeed(); Float_t anodePitch = fSegmentation->Dpz(dummy); // map the signal fMap->ClearMap(); fMap->SetThreshold(fCutAmplitude); fMap->FillMap(); Float_t noise; Float_t baseline; fResponse->GetNoiseParam(noise,baseline); Int_t nofFoundClusters = 0; Int_t i; Float_t **dfadc = new Float_t*[fNofAnodes]; for(i=0;iGetSignal(idx,l); if(l>0) fadc1=(Float_t)fMap->GetSignal(idx,l-1); if(l>0) dfadc[k][l-1] = fadc2-fadc1; } // samples } // anodes for(k=0;k=fMaxNofSamples) break; fadc=(float)fMap->GetSignal(idx,id); if(fadc > fadcmax) { fadcmax = fadc; imax = id;} if(fadc > (float)fCutAmplitude) { lthrt++; } // end if if(dfadc[k][id] > dfadcmax) { dfadcmax = dfadc[k][id]; imaxd = id; } // end if } // end for m it = imaxd; if(fMap->TestHit(idx,imax) == kEmpty) {it++; continue;} // cluster charge Int_t tstart = it-2; if(tstart < 0) tstart = 0; Bool_t ilcl = 0; if(lthrt >= lthrmint && lthra >= lthrmina) ilcl = 1; if(ilcl) { nofFoundClusters++; Int_t tstop = tstart; Float_t dfadcmin = 10000.; Int_t ij; for(ij=0; ij<20; ij++) { if(tstart+ij > 255) { tstop = 255; break; } fadc=(float)fMap->GetSignal(idx,tstart+ij); if((dfadc[k][tstart+ij] < dfadcmin) && (fadc > fCutAmplitude)) { tstop = tstart+ij+5; if(tstop > 255) tstop = 255; dfadcmin = dfadc[k][it+ij]; } // end if } // end for ij Float_t clusterCharge = 0.; Float_t clusterAnode = k+0.5; Float_t clusterTime = 0.; Int_t clusterMult = 0; Float_t clusterPeakAmplitude = 0.; Int_t its,peakpos = -1; Float_t n, baseline; fResponse->GetNoiseParam(n,baseline); for(its=tstart; its<=tstop; its++) { fadc=(float)fMap->GetSignal(idx,its); if(fadc>baseline) fadc -= baseline; else fadc = 0.; clusterCharge += fadc; // as a matter of fact we should take the peak // pos before FFT // to get the list of tracks !!! if(fadc > clusterPeakAmplitude) { clusterPeakAmplitude = fadc; //peakpos=fMap->GetHitIndex(idx,its); Int_t shift = (int)(fTimeCorr/fTimeStep); if(its>shift && its<(fMaxNofSamples-shift)) peakpos = fMap->GetHitIndex(idx,its+shift); else peakpos = fMap->GetHitIndex(idx,its); if(peakpos<0) peakpos =fMap->GetHitIndex(idx,its); } // end if clusterTime += fadc*its; if(fadc > 0) clusterMult++; if(its == tstop) { clusterTime /= (clusterCharge/fTimeStep); // ns if(clusterTime>fTimeCorr) clusterTime -=fTimeCorr; //ns } // end if } // end for its Float_t clusteranodePath = (clusterAnode - fNofAnodes/2)* anodePitch; Float_t clusterDriftPath = clusterTime*fDriftSpeed; clusterDriftPath = fSddLength-clusterDriftPath; if(clusterCharge <= 0.) break; AliITSRawClusterSDD clust(j+1,//i clusterAnode,clusterTime,//ff clusterCharge, //f clusterPeakAmplitude, //f peakpos, //i 0.,0.,clusterDriftPath,//fff clusteranodePath, //f clusterMult, //i 0,0,0,0,0,0,0);//7*i iTS->AddCluster(1,&clust); it = tstop; } // ilcl it++; } // while (samples) } // anodes } // detectors (2) for(i=0;iGetModule("ITS"); // retrieve the parameters Int_t fNofMaps = fSegmentation->Npz(); Int_t fMaxNofSamples = fSegmentation->Npx(); Int_t fNofAnodes = fNofMaps/2; Int_t dummy=0; Float_t fTimeStep = fSegmentation->Dpx( dummy ); Float_t fSddLength = fSegmentation->Dx(); Float_t fDriftSpeed = fResponse->DriftSpeed(); Float_t anodePitch = fSegmentation->Dpz( dummy ); Float_t n, baseline; fResponse->GetNoiseParam( n, baseline ); // map the signal fMap->ClearMap(); fMap->SetThreshold( fCutAmplitude ); fMap->FillMap(); Int_t nClu = 0; // cout << "Search cluster... "<< endl; for( Int_t j=0; j<2; j++ ){ for( Int_t k=0; kGetSignal( idx, l ); if( fadc > 0.0 ){ if( on == kFALSE && lGetSignal( idx, l+1 ); if( fadc1 < fadc ) continue; start = l; fmax = 0.; lmax = 0; time = 0.; charge = 0.; on = kTRUE; nTsteps = 0; } // end if on... nTsteps++ ; if( fadc > baseline ) fadc -= baseline; else fadc=0.; charge += fadc; time += fadc*l; if( fadc > fmax ){ fmax = fadc; lmax = l; Int_t shift = (Int_t)(fTimeCorr/fTimeStep + 0.5); if( l > shift && l < (fMaxNofSamples-shift) ) peakpos = fMap->GetHitIndex( idx, l+shift ); else peakpos = fMap->GetHitIndex( idx, l ); if( peakpos < 0) peakpos = fMap->GetHitIndex( idx, l ); } // end if fadc }else{ // end fadc>0 if( on == kTRUE ){ if( nTsteps > 2 ){ // min # of timesteps for a RawCluster // Found a RawCluster... Int_t stop = l-1; time /= (charge/fTimeStep); // ns // time = lmax*fTimeStep; // ns if( time > fTimeCorr ) time -= fTimeCorr; // ns Float_t anodePath = (anode - fNofAnodes/2)*anodePitch; Float_t driftPath = time*fDriftSpeed; driftPath = fSddLength-driftPath; AliITSRawClusterSDD clust(j+1,anode,time,charge, fmax, peakpos,0.,0., driftPath,anodePath, nTsteps,start,stop, start, stop, 1, k, k ); iTS->AddCluster( 1, &clust ); // clust.PrintInfo(); nClu++; } // end if nTsteps on = kFALSE; } // end if on==kTRUE } // end if fadc>0 } // samples } // anodes } // wings // cout << "# Rawclusters " << nClu << endl; return; } //_______________________________________________________________________ Int_t AliITSClusterFinderSDD::SearchPeak(Float_t *spect,Int_t xdim,Int_t zdim, Int_t *peakX, Int_t *peakZ, Float_t *peakAmp, Float_t minpeak ){ // search peaks on a 2D cluster Int_t npeak = 0; // # peaks Int_t i,j; // search peaks for( Int_t z=1; z= spect[(x+1)*zdim+z ] && sxz >= spect[(x-1)*zdim+z ] && sxz >= spect[x*zdim +z+1] && sxz >= spect[x*zdim +z-1] && sxz >= spect[(x+1)*zdim+z+1] && sxz >= spect[(x+1)*zdim+z-1] && sxz >= spect[(x-1)*zdim+z+1] && sxz >= spect[(x-1)*zdim+z-1] ){ // peak found peakX[npeak] = x; peakZ[npeak] = z; peakAmp[npeak] = sxz; npeak++; } // end if .... } // end for x } // end for z // search groups of peaks with same amplitude. Int_t *flag = new Int_t[npeak]; for( i=0; i 0 ) continue; if( peakAmp[i] == peakAmp[j] && TMath::Abs(peakX[i]-peakX[j])<=1 && TMath::Abs(peakZ[i]-peakZ[j])<=1 ){ if( flag[i] == 0) flag[i] = i+1; flag[j] = flag[i]; } // end if ... } // end for j } // end for i // make average of peak groups for( i=0; i 1 ){ peakX[i] /= nFlag; peakZ[i] /= nFlag; } // end fi nFlag } // end for i delete [] flag; return( npeak ); } //______________________________________________________________________ void AliITSClusterFinderSDD::PeakFunc( Int_t xdim, Int_t zdim, Float_t *par, Float_t *spe, Float_t *integral){ // function used to fit the clusters // par -> parameters.. // par[0] number of peaks. // for each peak i=1, ..., par[0] // par[i] = Ampl. // par[i+1] = xpos // par[i+2] = zpos // par[i+3] = tau // par[i+4] = sigma. Int_t electronics = fResponse->Electronics(); // 1 = PASCAL, 2 = OLA const Int_t knParam = 5; Int_t npeak = (Int_t)par[0]; memset( spe, 0, sizeof( Float_t )*zdim*xdim ); Int_t k = 1; for( Int_t i=0; i0.) ? par[k]*x2*exp(-x2+1.-z2) :0.0; // RCCR2 // signal =(x2>0.) ? par[k]*x2*x2*exp(-2*x2+2.-z2 ):0.0;//RCCR }else if( electronics == 2 ) { // OLA x2 = (x-par[k+1])*(x-par[k+1])/T2; signal = par[k] * exp( -x2 - z2 ); } else { cout << "Wrong SDD Electronics =" << electronics << endl; // exit( 1 ); } // end if electronicx spe[x*zdim+z] += signal; if( integral != 0 ) integral[i] += signal; } // end for x } // end for z k += knParam; } // end for i return; } //__________________________________________________________________________ Float_t AliITSClusterFinderSDD::ChiSqr( Int_t xdim, Int_t zdim, Float_t *spe, Float_t *speFit ){ // EVALUATES UNNORMALIZED CHI-SQUARED Float_t chi2 = 0.; for( Int_t z=0; z 1.0E-6 ) if ( fabs( delta/p1 ) < 1.0E-4 ) delta = p1/1000; else delta = (Float_t)1.0E-4; // EVALUATE CHI-SQUARED AT FIRST TWO SEARCH POINTS PeakFunc( xdim, zdim, param, speFit ); chisq1 = ChiSqr( xdim, zdim, spe, speFit ); p2 = p1+delta; param[k] = p2; PeakFunc( xdim, zdim, param, speFit ); chisq2 = ChiSqr( xdim, zdim, spe, speFit ); if( chisq1 < chisq2 ){ // REVERSE DIRECTION OF SEARCH IF CHI-SQUARED IS INCREASING delta = -delta; t = p1; p1 = p2; p2 = t; t = chisq1; chisq1 = chisq2; chisq2 = t; } // end if i = 1; nnn = 0; do { // INCREMENT param(K) UNTIL CHI-SQUARED STARTS TO INCREASE nnn++; p3 = p2 + delta; mmm = nnn - (nnn/5)*5; // multiplo de 5 if( mmm == 0 ){ d1 = delta; // INCREASE STEP SIZE IF STEPPING TOWARDS MINIMUM IS TOO SLOW delta *= 5; } // end if param[k] = p3; // Constrain paramiters Int_t kpos = (k-1) % knParam; switch( kpos ){ case 0 : if( param[k] <= 20 ) param[k] = fMinPeak; break; case 1 : if( fabs( param[k] - prm0[k] ) > 1.5 ) param[k] = prm0[k]; break; case 2 : if( fabs( param[k] - prm0[k] ) > 1. ) param[k] = prm0[k]; break; case 3 : if( param[k] < .5 ) param[k] = .5; break; case 4 : if( param[k] < .288 ) param[k] = .288; // 1/sqrt(12) = 0.288 if( param[k] > zdim*.5 ) param[k] = zdim*.5; break; }; // end switch PeakFunc( xdim, zdim, param, speFit ); chisq3 = ChiSqr( xdim, zdim, spe, speFit ); if( chisq3 < chisq2 && nnn < 50 ){ p1 = p2; p2 = p3; chisq1 = chisq2; chisq2 = chisq3; }else i=0; } while( i ); // FIND MINIMUM OF PARABOLA DEFINED BY LAST THREE POINTS a = chisq1*(p2-p3)+chisq2*(p3-p1)+chisq3*(p1-p2); b = chisq1*(p2*p2-p3*p3)+chisq2*(p3*p3-p1*p1)+chisq3*(p1*p1-p2*p2); if( a!=0 ) p0 = (Float_t)(0.5*b/a); else p0 = 10000; //--IN CASE OF NEARLY EQUAL CHI-SQUARED AND TOO SMALL STEP SIZE PREVENT // ERRONEOUS EVALUATION OF PARABOLA MINIMUM //---NEXT TWO LINES CAN BE OMITTED FOR HIGHER PRECISION MACHINES //dp = (Float_t) max (fabs(p3-p2), fabs(p2-p1)); //if( fabs( p2-p0 ) > dp ) p0 = p2; param[k] = p0; // Constrain paramiters Int_t kpos = (k-1) % knParam; switch( kpos ){ case 0 : if( param[k] <= 20 ) param[k] = fMinPeak; break; case 1 : if( fabs( param[k] - prm0[k] ) > 1.5 ) param[k] = prm0[k]; break; case 2 : if( fabs( param[k] - prm0[k] ) > 1. ) param[k] = prm0[k]; break; case 3 : if( param[k] < .5 ) param[k] = .5; break; case 4 : if( param[k] < .288 ) param[k] = .288; // 1/sqrt(12) = 0.288 if( param[k] > zdim*.5 ) param[k] = zdim*.5; break; }; // end switch PeakFunc( xdim, zdim, param, speFit ); chisqt = ChiSqr( xdim, zdim, spe, speFit ); // DO NOT ALLOW ERRONEOUS INTERPOLATION if( chisqt <= *chisqr ) *chisqr = chisqt; else param[k] = prm0[k]; // OPTIMIZE SEARCH STEP FOR EVENTUAL NEXT CALL OF MINIM steprm[k] = (param[k]-prm0[k])/5; if( steprm[k] >= d1 ) steprm[k] = d1/5; } // end for k // EVALUATE FIT AND CHI-SQUARED FOR OPTIMIZED PARAMETERS PeakFunc( xdim, zdim, param, speFit ); *chisqr = ChiSqr( xdim, zdim, spe, speFit ); return; } //_________________________________________________________________________ Int_t AliITSClusterFinderSDD::NoLinearFit( Int_t xdim, Int_t zdim, Float_t *param, Float_t *spe, Int_t *niter, Float_t *chir ){ // fit method from Comput. Phys. Commun 46(1987) 149 const Float_t kchilmt = 0.01; // relative accuracy const Int_t knel = 3; // for parabolic minimization const Int_t knstop = 50; // Max. iteration number const Int_t knParam = 5; Int_t npeak = (Int_t)param[0]; // RETURN IF NUMBER OF DEGREES OF FREEDOM IS NOT POSITIVE if( (xdim*zdim - npeak*knParam) <= 0 ) return( -1 ); Float_t degFree = (xdim*zdim - npeak*knParam)-1; Int_t n, k, iterNum = 0; Float_t *prm0 = new Float_t[npeak*knParam+1]; Float_t *step = new Float_t[npeak*knParam+1]; Float_t *schi = new Float_t[npeak*knParam+1]; Float_t *sprm[3]; sprm[0] = new Float_t[npeak*knParam+1]; sprm[1] = new Float_t[npeak*knParam+1]; sprm[2] = new Float_t[npeak*knParam+1]; Float_t chi0, chi1, reldif, a, b, prmin, dp; Float_t *speFit = new Float_t[ xdim*zdim ]; PeakFunc( xdim, zdim, param, speFit ); chi0 = ChiSqr( xdim, zdim, spe, speFit ); chi1 = chi0; for( k=1; k<(npeak*knParam+1); k++) prm0[k] = param[k]; for( k=1 ; k<(npeak*knParam+1); k+=knParam ){ step[k] = param[k] / 20.0 ; step[k+1] = param[k+1] / 50.0; step[k+2] = param[k+2] / 50.0; step[k+3] = param[k+3] / 20.0; step[k+4] = param[k+4] / 20.0; } // end for k Int_t out = 0; do{ iterNum++; chi0 = chi1; Minim( xdim, zdim, param, prm0, step, &chi1, spe, speFit ); reldif = ( chi1 > 0 ) ? ((Float_t) fabs( chi1-chi0)/chi1 ) : 0; // EXIT conditions if( reldif < (float) kchilmt ){ *chir = (chi1>0) ? (float) TMath::Sqrt (chi1/degFree) :0; *niter = iterNum; out = 0; break; } // end if if( (reldif < (float)(5*kchilmt)) && (iterNum > knstop) ){ *chir = (chi1>0) ?(float) TMath::Sqrt (chi1/degFree):0; *niter = iterNum; out = 0; break; } // end if if( iterNum > 5*knstop ){ *chir = (chi1>0) ?(float) TMath::Sqrt (chi1/degFree):0; *niter = iterNum; out = 1; break; } // end if if( iterNum <= knel ) continue; n = iterNum - (iterNum/knel)*knel; // EXTRAPOLATION LIMIT COUNTER N if( n > 3 || n == 0 ) continue; schi[n-1] = chi1; for( k=1; k<(npeak*knParam+1); k++ ) sprm[n-1][k] = param[k]; if( n != 3 ) continue; // -EVALUATE EXTRAPOLATED VALUE OF EACH PARAMETER BY FINDING MINIMUM OF // PARABOLA DEFINED BY LAST THREE CALLS OF MINIM for( k=1; k<(npeak*knParam+1); k++ ){ Float_t tmp0 = sprm[0][k]; Float_t tmp1 = sprm[1][k]; Float_t tmp2 = sprm[2][k]; a = schi[0]*(tmp1-tmp2) + schi[1]*(tmp2-tmp0); a += (schi[2]*(tmp0-tmp1)); b = schi[0]*(tmp1*tmp1-tmp2*tmp2); b += (schi[1]*(tmp2*tmp2-tmp0*tmp0)+(schi[2]* (tmp0*tmp0-tmp1*tmp1))); if ((double)a < 1.0E-6) prmin = 0; else prmin = (float) (0.5*b/a); dp = 5*(tmp2-tmp0); if( fabs(prmin-tmp2) > fabs(dp) ) prmin = tmp2+dp; param[k] = prmin; step[k] = dp/10; // OPTIMIZE SEARCH STEP } // end for k } while( kTRUE ); delete [] prm0; delete [] step; delete [] schi; delete [] sprm[0]; delete [] sprm[1]; delete [] sprm[2]; delete [] speFit; return( out ); } //______________________________________________________________________ void AliITSClusterFinderSDD::ResolveClustersE(){ // The function to resolve clusters if the clusters overlapping exists Int_t i; static AliITS *iTS = (AliITS*)gAlice->GetModule( "ITS" ); // get number of clusters for this module Int_t nofClusters = fClusters->GetEntriesFast(); nofClusters -= fNclusters; Int_t fNofMaps = fSegmentation->Npz(); Int_t fNofAnodes = fNofMaps/2; Int_t fMaxNofSamples = fSegmentation->Npx(); Int_t dummy=0; Double_t fTimeStep = fSegmentation->Dpx( dummy ); Double_t fSddLength = fSegmentation->Dx(); Double_t fDriftSpeed = fResponse->DriftSpeed(); Double_t anodePitch = fSegmentation->Dpz( dummy ); Float_t n, baseline; fResponse->GetNoiseParam( n, baseline ); Int_t electronics = fResponse->Electronics(); // 1 = PASCAL, 2 = OLA for( Int_t j=0; jAt(j); Int_t astart = clusterJ->Astart(); Int_t astop = clusterJ->Astop(); Int_t tstart = clusterJ->Tstartf(); Int_t tstop = clusterJ->Tstopf(); Int_t wing = (Int_t)clusterJ->W(); if( wing == 2 ){ astart += fNofAnodes; astop += fNofAnodes; } // end if Int_t xdim = tstop-tstart+3; Int_t zdim = astop-astart+3; if(xdim > 50 || zdim > 30) { cout << "Warning: xdim: " << xdim << ", zdim: " << zdim << endl; continue; } Float_t *sp = new Float_t[ xdim*zdim+1 ]; memset( sp, 0, sizeof(Float_t)*(xdim*zdim+1) ); // make a local map from cluster region for( Int_t ianode=astart; ianode<=astop; ianode++ ){ for( Int_t itime=tstart; itime<=tstop; itime++ ){ Float_t fadc = fMap->GetSignal( ianode, itime ); if( fadc > baseline ) fadc -= (Double_t)baseline; else fadc = 0.; Int_t index = (itime-tstart+1)*zdim+(ianode-astart+1); sp[index] = fadc; } // time loop } // anode loop // search peaks on cluster const Int_t kNp = 150; Int_t peakX1[kNp]; Int_t peakZ1[kNp]; Float_t peakAmp1[kNp]; Int_t npeak = SearchPeak(sp,xdim,zdim,peakX1,peakZ1,peakAmp1,fMinPeak); // if multiple peaks, split cluster if( npeak >= 1 ) { // cout << "npeak " << npeak << endl; // clusterJ->PrintInfo(); Float_t *par = new Float_t[npeak*5+1]; par[0] = (Float_t)npeak; // Initial parameters in cell dimentions Int_t k1 = 1; for( i=0; i charge for each peak PeakFunc( xdim, zdim, par, sp, integral ); k1 = 1; for( i=0; i shift && newiTime < (fMaxNofSamples-shift) ) shift = 0; Int_t peakpos = fMap->GetHitIndex( newAnode, newiTime+shift ); clusterI.SetPeakPos( peakpos ); clusterI.SetPeakAmpl( peakAmp1[i] ); Float_t newAnodef = peakZ[i] - 0.5 + astart; Float_t newiTimef = peakX[i] - 1 + tstart; if( wing == 2 ) newAnodef -= fNofAnodes; Float_t anodePath = (newAnodef - fNofAnodes/2)*anodePitch; newiTimef *= fTimeStep; if( newiTimef > fTimeCorr ) newiTimef -= fTimeCorr; if( electronics == 1 ){ // newiTimef *= 0.999438; // PASCAL // newiTimef += (6./fDriftSpeed - newiTimef/3000.); }else if( electronics == 2 ) newiTimef *= 0.99714; // OLA Float_t driftPath = fSddLength - newiTimef * fDriftSpeed; Float_t sign = ( wing == 1 ) ? -1. : 1.; clusterI.SetX( driftPath*sign * 0.0001 ); clusterI.SetZ( anodePath * 0.0001 ); clusterI.SetAnode( newAnodef ); clusterI.SetTime( newiTimef ); clusterI.SetAsigma( sigma[i]*anodePitch ); clusterI.SetTsigma( tau[i]*fTimeStep ); clusterI.SetQ( integral[i] ); // clusterI.PrintInfo(); iTS->AddCluster( 1, &clusterI ); } // end for i fClusters->RemoveAt( j ); delete [] par; } else { // something odd cout << " --- Peak not found!!!! minpeak=" << fMinPeak<< " cluster peak=" << clusterJ->PeakAmpl() << " module=" << fModule << endl; clusterJ->PrintInfo(); cout << " xdim=" << xdim-2 << " zdim=" << zdim-2 << endl << endl; } delete [] sp; } // cluster loop fClusters->Compress(); // fMap->ClearMap(); } //________________________________________________________________________ void AliITSClusterFinderSDD::GroupClusters(){ // group clusters Int_t dummy=0; Float_t fTimeStep = fSegmentation->Dpx(dummy); // get number of clusters for this module Int_t nofClusters = fClusters->GetEntriesFast(); nofClusters -= fNclusters; AliITSRawClusterSDD *clusterI; AliITSRawClusterSDD *clusterJ; Int_t *label = new Int_t [nofClusters]; Int_t i,j; for(i=0; iAt(i); clusterJ = (AliITSRawClusterSDD*) fClusters->At(j); // 1.3 good if(clusterI->T() < fTimeStep*60) fDAnode = 4.2; // TB 3.2 if(clusterI->T() < fTimeStep*10) fDAnode = 1.5; // TB 1. Bool_t pair = clusterI->Brother(clusterJ,fDAnode,fDTime); if(!pair) continue; // clusterI->PrintInfo(); // clusterJ->PrintInfo(); clusterI->Add(clusterJ); label[j] = 1; fClusters->RemoveAt(j); j=i; // <- Ernesto } // J clusters label[i] = 1; } // I clusters fClusters->Compress(); delete [] label; return; } //________________________________________________________________________ void AliITSClusterFinderSDD::SelectClusters(){ // get number of clusters for this module Int_t nofClusters = fClusters->GetEntriesFast(); nofClusters -= fNclusters; Int_t i; for(i=0; iAt(i); Int_t rmflg = 0; Float_t wy = 0.; if(clusterI->Anodes() != 0.) { wy = ((Float_t) clusterI->Samples())/clusterI->Anodes(); } // end if Int_t amp = (Int_t) clusterI->PeakAmpl(); Int_t cha = (Int_t) clusterI->Q(); if(amp < fMinPeak) rmflg = 1; if(cha < fMinCharge) rmflg = 1; if(wy < fMinNCells) rmflg = 1; //if(wy > fMaxNCells) rmflg = 1; if(rmflg) fClusters->RemoveAt(i); } // I clusters fClusters->Compress(); return; } //__________________________________________________________________________ void AliITSClusterFinderSDD::ResolveClusters(){ // The function to resolve clusters if the clusters overlapping exists /* AliITS *iTS=(AliITS*)gAlice->GetModule("ITS"); // get number of clusters for this module Int_t nofClusters = fClusters->GetEntriesFast(); nofClusters -= fNclusters; //cout<<"Resolve Cl: nofClusters, fNclusters ="<Npz(); Int_t fNofAnodes = fNofMaps/2; Int_t dummy=0; Double_t fTimeStep = fSegmentation->Dpx(dummy); Double_t fSddLength = fSegmentation->Dx(); Double_t fDriftSpeed = fResponse->DriftSpeed(); Double_t anodePitch = fSegmentation->Dpz(dummy); Float_t n, baseline; fResponse->GetNoiseParam(n,baseline); Float_t dzz_1A = anodePitch * anodePitch / 12; // fill Map of signals fMap->FillMap(); Int_t j,i,ii,ianode,anode,itime; Int_t wing,astart,astop,tstart,tstop,nanode; Double_t fadc,ClusterTime; Double_t q[400],x[400],z[400]; // digit charges and coordinates for(j=0; jAt(j); Int_t ndigits = 0; astart=clusterJ->Astart(); astop=clusterJ->Astop(); tstart=clusterJ->Tstartf(); tstop=clusterJ->Tstopf(); nanode=clusterJ->Anodes(); // <- Ernesto wing=(Int_t)clusterJ->W(); if(wing == 2) { astart += fNofAnodes; astop += fNofAnodes; } // end if // cout<<"astart,astop,tstart,tstop ="<GetSignal(ianode,itime); if(fadc>baseline) { fadc-=(Double_t)baseline; q[ndigits] = fadc*(fTimeStep/160); // KeV anode = ianode; if(wing == 2) anode -= fNofAnodes; z[ndigits] = (anode + 0.5 - fNofAnodes/2)*anodePitch; ClusterTime = itime*fTimeStep; if(ClusterTime > fTimeCorr) ClusterTime -= fTimeCorr;// ns x[ndigits] = fSddLength - ClusterTime*fDriftSpeed; if(wing == 1) x[ndigits] *= (-1); // cout<<"ianode,itime,fadc ="< 0.0001) continue; if (fabs(x2-x2_old) > 0.0001) continue; if (fabs(r-r_old)/5 > 0.001) continue; a1=r*qq*pitchx*pitchz/(2*TMath::ACos(-1.)*sigma2); Double_t a2=a1*(1-r)/r; qfit[0]=a1; xfit[0]=x1*cosa - z12*sina; zfit[0]=x1*sina + z12*cosa; qfit[1]=a2; xfit[1]=x2*cosa - z12*sina; zfit[1]=x2*sina + z12*cosa; nfhits=2; break; // Ok ! } // end for i if (i==33) cerr<<"No more iterations ! "<0. && elps< 0.3 && nfhits == 1) cout<<" small elps, nfh=1 =" <PrintInfo(); cout << " in: " << endl; for (i=0; iPeakAmpl(); Float_t peakpos = clusterJ->PeakPos(); Float_t clusteranodePath = (Anode - fNofAnodes/2)*anodePitch; Float_t clusterDriftPath = Time*fDriftSpeed; clusterDriftPath = fSddLength-clusterDriftPath; AliITSRawClusterSDD *clust = new AliITSRawClusterSDD(wing,Anode, Time,qfit[i], clusterPeakAmplitude,peakpos, 0.,0.,clusterDriftPath, clusteranodePath,clusterJ->Samples()/2 ,tstart,tstop,0,0,0,astart,astop); clust->PrintInfo(); iTS->AddCluster(1,clust); // cout<<"new cluster added: tstart,tstop,astart,astop,x,ncl =" // <RemoveAt(j); } // if nfhits = 2 } // cluster loop fClusters->Compress(); fMap->ClearMap(); */ return; } //______________________________________________________________________ void AliITSClusterFinderSDD::GetRecPoints(){ // get rec points static AliITS *iTS=(AliITS*)gAlice->GetModule("ITS"); // get number of clusters for this module Int_t nofClusters = fClusters->GetEntriesFast(); nofClusters -= fNclusters; const Float_t kconvGeV = 1.e-6; // GeV -> KeV const Float_t kconv = 1.0e-4; const Float_t kRMSx = 38.0*kconv; // microns->cm ITS TDR Table 1.3 const Float_t kRMSz = 28.0*kconv; // microns->cm ITS TDR Table 1.3 Int_t i; Int_t ix, iz, idx=-1; AliITSdigitSDD *dig=0; Int_t ndigits=fDigits->GetEntriesFast(); for(i=0; iAt(i); if(!clusterI) Error("SDD: GetRecPoints","i clusterI ",i,clusterI); if(clusterI) idx=clusterI->PeakPos(); if(idx>ndigits) Error("SDD: GetRecPoints","idx ndigits",idx,ndigits); // try peak neighbours - to be done if(idx&&idx<= ndigits) dig =(AliITSdigitSDD*)fDigits->UncheckedAt(idx); if(!dig) { // try cog fSegmentation->GetPadIxz(clusterI->X(),clusterI->Z(),ix,iz); dig = (AliITSdigitSDD*)fMap->GetHit(iz-1,ix-1); // if null try neighbours if (!dig) dig = (AliITSdigitSDD*)fMap->GetHit(iz-1,ix); if (!dig) dig = (AliITSdigitSDD*)fMap->GetHit(iz-1,ix+1); if (!dig) printf("SDD: cannot assign the track number!\n"); } // end if !dig AliITSRecPoint rnew; rnew.SetX(clusterI->X()); rnew.SetZ(clusterI->Z()); rnew.SetQ(clusterI->Q()); // in KeV - should be ADC rnew.SetdEdX(kconvGeV*clusterI->Q()); rnew.SetSigmaX2(kRMSx*kRMSx); rnew.SetSigmaZ2(kRMSz*kRMSz); if(dig) rnew.fTracks[0]=dig->fTracks[0]; if(dig) rnew.fTracks[1]=dig->fTracks[1]; if(dig) rnew.fTracks[2]=dig->fTracks[2]; //printf("SDD: i %d track1 track2 track3 %d %d %d x y %f %f\n", // i,rnew.fTracks[0],rnew.fTracks[1],rnew.fTracks[2],c // lusterI->X(),clusterI->Z()); iTS->AddRecPoint(rnew); } // I clusters // fMap->ClearMap(); } //______________________________________________________________________ void AliITSClusterFinderSDD::FindRawClusters(Int_t mod){ // find raw clusters fModule = mod; Find1DClustersE(); GroupClusters(); SelectClusters(); ResolveClustersE(); GetRecPoints(); } //_______________________________________________________________________ void AliITSClusterFinderSDD::Print(){ // Print SDD cluster finder Parameters cout << "**************************************************" << endl; cout << " Silicon Drift Detector Cluster Finder Parameters " << endl; cout << "**************************************************" << endl; cout << "Number of Clusters: " << fNclusters << endl; cout << "Anode Tolerance: " << fDAnode << endl; cout << "Time Tolerance: " << fDTime << endl; cout << "Time correction (electronics): " << fTimeCorr << endl; cout << "Cut Amplitude (threshold): " << fCutAmplitude << endl; cout << "Minimum Amplitude: " << fMinPeak << endl; cout << "Minimum Charge: " << fMinCharge << endl; cout << "Minimum number of cells/clusters: " << fMinNCells << endl; cout << "Maximum number of cells/clusters: " << fMaxNCells << endl; cout << "**************************************************" << endl; }