X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=ITS%2FAliITSClusterFinderSDD.cxx;h=6430c48401dbf5250c004ecdd72cfcc3b610c72b;hb=8e50d897557aadaff83f0809747a23e28dfab8ca;hp=bb42a2e36da415f880d0c7ab0fdf50da4f1638d5;hpb=b9d0a01d7a0723a09071b0b56200d72f59a9c2b6;p=u%2Fmrichter%2FAliRoot.git diff --git a/ITS/AliITSClusterFinderSDD.cxx b/ITS/AliITSClusterFinderSDD.cxx index bb42a2e36da..6430c48401d 100644 --- a/ITS/AliITSClusterFinderSDD.cxx +++ b/ITS/AliITSClusterFinderSDD.cxx @@ -14,137 +14,109 @@ **************************************************************************/ /* $Id$ - $Log$ - Revision 1.23.4.2 2002/10/14 13:14:07 hristov - Updating VirtualMC to v3-09-02 - - Revision 1.26 2002/09/09 17:23:28 nilsen - Minor changes in support of changes to AliITSdigitS?D class'. - - Revision 1.25 2002/05/10 22:29:40 nilsen - Change my Massimo Masera in the default constructor to bring things into - compliance. - - Revision 1.24 2002/04/24 22:02:31 nilsen - New SDigits and Digits routines, and related changes, (including new - noise values). - */ +/////////////////////////////////////////////////////////////////////////// +// Cluster finder // +// for Silicon // +// Drift Detector // +////////////////////////////////////////////////////////////////////////// -#include -#include -#include -#include #include "AliITSClusterFinderSDD.h" #include "AliITSMapA1.h" -#include "AliITS.h" -#include "AliITSdigit.h" -#include "AliITSRawCluster.h" +#include "AliITSRawClusterSDD.h" #include "AliITSRecPoint.h" -#include "AliITSsegmentation.h" -#include "AliITSresponseSDD.h" -#include "AliRun.h" +#include "AliITSdigitSDD.h" +#include "AliITSDetTypeRec.h" +#include "AliITSCalibrationSDD.h" +#include "AliITSsegmentationSDD.h" +#include "AliLog.h" ClassImp(AliITSClusterFinderSDD) //______________________________________________________________________ -AliITSClusterFinderSDD::AliITSClusterFinderSDD(AliITSsegmentation *seg, - AliITSresponse *response, +AliITSClusterFinderSDD::AliITSClusterFinderSDD(): +AliITSClusterFinder(), +fNclusters(0), +fDAnode(0.0), +fDTime(0.0), +fTimeCorr(0.0), +fCutAmplitude(0), +fMinPeak(0), +fMinCharge(0), +fMinNCells(0), +fMaxNCells(0){ + // default constructor +} +//______________________________________________________________________ +AliITSClusterFinderSDD::AliITSClusterFinderSDD(AliITSDetTypeRec* dettyp, TClonesArray *digits, - TClonesArray *recp){ + TClonesArray *recp): +AliITSClusterFinder(dettyp), +fNclusters(0), +fDAnode(0.0), +fDTime(0.0), +fTimeCorr(0.0), +fCutAmplitude(0), +fMinPeak(0), +fMinCharge(0), +fMinNCells(0), +fMaxNCells(0){ // standard constructor - fSegmentation = seg; - fResponse = response; - fDigits = digits; - fClusters = recp; - fNclusters = fClusters->GetEntriesFast(); - SetCutAmplitude(); + SetDigits(digits); + SetClusters(recp); + SetCutAmplitude(fDetTypeRec->GetITSgeom()->GetStartSDD()); SetDAnode(); SetDTime(); - SetMinPeak((Int_t)(((AliITSresponseSDD*)fResponse)->GetNoiseAfterElectronics()*5)); + SetMinPeak((Int_t)(((AliITSCalibrationSDD*)GetResp(fDetTypeRec->GetITSgeom()->GetStartSDD()))-> + GetNoiseAfterElectronics()*5)); // SetMinPeak(); SetMinNCells(); SetMaxNCells(); SetTimeCorr(); SetMinCharge(); - fMap = new AliITSMapA1(fSegmentation,fDigits,fCutAmplitude); + SetMap(new AliITSMapA1(GetSeg(),Digits(),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){ +void AliITSClusterFinderSDD::SetCutAmplitude(Int_t mod,Double_t nsigma){ // set the signal threshold for cluster finder - Float_t baseline,noise,noise_after_el; + Double_t baseline,noise,noiseAfterEl; - fResponse->GetNoiseParam(noise,baseline); - noise_after_el = ((AliITSresponseSDD*)fResponse)->GetNoiseAfterElectronics(); - fCutAmplitude = (Int_t)((baseline + nsigma*noise_after_el)); + GetResp(mod)->GetNoiseParam(noise,baseline); + noiseAfterEl = ((AliITSCalibrationSDD*)GetResp(mod))->GetNoiseAfterElectronics(); + fCutAmplitude = (Int_t)((baseline + nsigma*noiseAfterEl)); } //______________________________________________________________________ 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 fNofMaps = GetSeg()->Npz(); + Int_t fMaxNofSamples = GetSeg()->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); + Double_t fTimeStep = GetSeg()->Dpx(dummy); + Double_t fSddLength = GetSeg()->Dx(); + Double_t fDriftSpeed = GetResp(fModule)->GetDriftSpeed(); + Double_t anodePitch = GetSeg()->Dpz(dummy); // map the signal - fMap->ClearMap(); - fMap->SetThreshold(fCutAmplitude); - fMap->FillMap(); + Map()->ClearMap(); + Map()->SetThreshold(fCutAmplitude); + Map()->FillMap(); - Float_t noise; - Float_t baseline; - fResponse->GetNoiseParam(noise,baseline); + Double_t noise; + Double_t baseline; + GetResp(fModule)->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); + fadc2=(Double_t)Map()->GetSignal(idx,l); + if(l>0) fadc1=(Double_t)Map()->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); + fadc=(float)Map()->GetSignal(idx,id); if(fadc > fadcmax) { fadcmax = fadc; imax = id;} - if(fadc > (float)fCutAmplitude) { - lthrt++; - } // end if + if(fadc > (float)fCutAmplitude)lthrt++; 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;} + if(Map()->TestHit(idx,imax) == kEmpty) {it++; continue;} // cluster charge Int_t tstart = it-2; if(tstart < 0) tstart = 0; @@ -197,11 +167,11 @@ void AliITSClusterFinderSDD::Find1DClusters(){ if(ilcl) { nofFoundClusters++; Int_t tstop = tstart; - Float_t dfadcmin = 10000.; + Double_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); + fadc=(float)Map()->GetSignal(idx,tstart+ij); if((dfadc[k][tstart+ij] < dfadcmin) && (fadc > fCutAmplitude)) { tstop = tstart+ij+5; @@ -210,16 +180,16 @@ void AliITSClusterFinderSDD::Find1DClusters(){ } // end if } // end for ij - Float_t clusterCharge = 0.; - Float_t clusterAnode = k+0.5; - Float_t clusterTime = 0.; + Double_t clusterCharge = 0.; + Double_t clusterAnode = k+0.5; + Double_t clusterTime = 0.; Int_t clusterMult = 0; - Float_t clusterPeakAmplitude = 0.; + Double_t clusterPeakAmplitude = 0.; Int_t its,peakpos = -1; - Float_t n, baseline; - fResponse->GetNoiseParam(n,baseline); + Double_t n, baseline; + GetResp(fModule)->GetNoiseParam(n,baseline); for(its=tstart; its<=tstop; its++) { - fadc=(float)fMap->GetSignal(idx,its); + fadc=(float)Map()->GetSignal(idx,its); if(fadc>baseline) fadc -= baseline; else fadc = 0.; clusterCharge += fadc; @@ -228,12 +198,12 @@ void AliITSClusterFinderSDD::Find1DClusters(){ // to get the list of tracks !!! if(fadc > clusterPeakAmplitude) { clusterPeakAmplitude = fadc; - //peakpos=fMap->GetHitIndex(idx,its); + //peakpos=Map()->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); + peakpos = Map()->GetHitIndex(idx,its+shift); + else peakpos = Map()->GetHitIndex(idx,its); + if(peakpos<0) peakpos =Map()->GetHitIndex(idx,its); } // end if clusterTime += fadc*its; if(fadc > 0) clusterMult++; @@ -244,9 +214,9 @@ void AliITSClusterFinderSDD::Find1DClusters(){ } // end if } // end for its - Float_t clusteranodePath = (clusterAnode - fNofAnodes/2)* - anodePitch; - Float_t clusterDriftPath = clusterTime*fDriftSpeed; + Double_t clusteranodePath = (clusterAnode - fNofAnodes/2)* + anodePitch; + Double_t clusterDriftPath = clusterTime*fDriftSpeed; clusterDriftPath = fSddLength-clusterDriftPath; if(clusterCharge <= 0.) break; AliITSRawClusterSDD clust(j+1,//i @@ -258,7 +228,7 @@ void AliITSClusterFinderSDD::Find1DClusters(){ clusteranodePath, //f clusterMult, //i 0,0,0,0,0,0,0);//7*i - iTS->AddCluster(1,&clust); + fDetTypeRec->AddCluster(1,&clust); it = tstop; } // ilcl it++; @@ -271,29 +241,25 @@ void AliITSClusterFinderSDD::Find1DClusters(){ return; } - - - //______________________________________________________________________ void AliITSClusterFinderSDD::Find1DClustersE(){ // 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 fNofMaps = GetSeg()->Npz(); + Int_t fMaxNofSamples = GetSeg()->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 ); + Double_t fTimeStep = GetSeg()->Dpx( dummy ); + Double_t fSddLength = GetSeg()->Dx(); + Double_t fDriftSpeed = GetResp(fModule)->GetDriftSpeed(); + Double_t anodePitch = GetSeg()->Dpz( dummy ); + Double_t n, baseline; + GetResp(fModule)->GetNoiseParam( n, baseline ); // map the signal - fMap->ClearMap(); - fMap->SetThreshold( fCutAmplitude ); - fMap->FillMap(); - + Map()->ClearMap(); + Map()->SetThreshold( fCutAmplitude ); + Map()->FillMap(); + Int_t nClu = 0; // cout << "Search cluster... "<< endl; for( Int_t j=0; j<2; j++ ){ @@ -302,18 +268,18 @@ void AliITSClusterFinderSDD::Find1DClustersE(){ Bool_t on = kFALSE; Int_t start = 0; Int_t nTsteps = 0; - Float_t fmax = 0.; + Double_t fmax = 0.; Int_t lmax = 0; - Float_t charge = 0.; - Float_t time = 0.; - Float_t anode = k+0.5; + Double_t charge = 0.; + Double_t time = 0.; + Double_t anode = k+0.5; Int_t peakpos = -1; for( Int_t l=0; lGetSignal( idx, l ); + Double_t fadc = (Double_t)Map()->GetSignal( idx, l ); if( fadc > 0.0 ){ if( on == kFALSE && lGetSignal( idx, l+1 ); + Double_t fadc1 = (Double_t)Map()->GetSignal( idx, l+1 ); if( fadc1 < fadc ) continue; start = l; fmax = 0.; @@ -333,10 +299,10 @@ void AliITSClusterFinderSDD::Find1DClustersE(){ lmax = l; Int_t shift = (Int_t)(fTimeCorr/fTimeStep + 0.5); if( l > shift && l < (fMaxNofSamples-shift) ) - peakpos = fMap->GetHitIndex( idx, l+shift ); + peakpos = Map()->GetHitIndex( idx, l+shift ); else - peakpos = fMap->GetHitIndex( idx, l ); - if( peakpos < 0) peakpos = fMap->GetHitIndex( idx, l ); + peakpos = Map()->GetHitIndex( idx, l ); + if( peakpos < 0) peakpos = Map()->GetHitIndex(idx,l); } // end if fadc }else{ // end fadc>0 if( on == kTRUE ){ @@ -347,16 +313,16 @@ void AliITSClusterFinderSDD::Find1DClustersE(){ 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; + Double_t anodePath =(anode-fNofAnodes/2)*anodePitch; + Double_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(); + fDetTypeRec->AddCluster( 1, &clust ); + if(AliDebugLevel()>=5) clust.PrintInfo(); nClu++; } // end if nTsteps on = kFALSE; @@ -365,22 +331,22 @@ void AliITSClusterFinderSDD::Find1DClustersE(){ } // samples } // anodes } // wings - // cout << "# Rawclusters " << nClu << endl; + AliDebug(3,Form("# Rawclusters %d",nClu)); return; } //_______________________________________________________________________ -Int_t AliITSClusterFinderSDD::SearchPeak(Float_t *spect,Int_t xdim,Int_t zdim, +Int_t AliITSClusterFinderSDD::SearchPeak(Double_t *spect,Int_t xdim,Int_t zdim, Int_t *peakX, Int_t *peakZ, - Float_t *peakAmp, Float_t minpeak ){ + Double_t *peakAmp, Double_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 ] && @@ -438,8 +404,8 @@ Int_t AliITSClusterFinderSDD::SearchPeak(Float_t *spect,Int_t xdim,Int_t zdim, return( npeak ); } //______________________________________________________________________ -void AliITSClusterFinderSDD::PeakFunc( Int_t xdim, Int_t zdim, Float_t *par, - Float_t *spe, Float_t *integral){ +void AliITSClusterFinderSDD::PeakFunc( Int_t xdim, Int_t zdim, Double_t *par, + Double_t *spe, Double_t *integral){ // function used to fit the clusters // par -> parameters.. // par[0] number of peaks. @@ -449,32 +415,33 @@ void AliITSClusterFinderSDD::PeakFunc( Int_t xdim, Int_t zdim, Float_t *par, // par[i+2] = zpos // par[i+3] = tau // par[i+4] = sigma. - Int_t electronics = fResponse->Electronics(); // 1 = PASCAL, 2 = OLA + Int_t electronics = GetResp(fModule)->GetElectronics(); // 1 = PASCAL, 2 = OLA const Int_t knParam = 5; Int_t npeak = (Int_t)par[0]; - memset( spe, 0, sizeof( Float_t )*zdim*xdim ); + memset( spe, 0, sizeof( Double_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; + x2 = (x-par[k+1])*(x-par[k+1])/t2; signal = par[k] * exp( -x2 - z2 ); } else { - cout << "Wrong SDD Electronics =" << electronics << endl; + Warning("PeakFunc","Wrong SDD Electronics = %d", + electronics); // exit( 1 ); } // end if electronicx spe[x*zdim+z] += signal; @@ -486,27 +453,27 @@ void AliITSClusterFinderSDD::PeakFunc( Int_t xdim, Int_t zdim, Float_t *par, return; } //__________________________________________________________________________ -Float_t AliITSClusterFinderSDD::ChiSqr( Int_t xdim, Int_t zdim, Float_t *spe, - Float_t *speFit ){ +Double_t AliITSClusterFinderSDD::ChiSqr( Int_t xdim, Int_t zdim, Double_t *spe, + Double_t *speFit ) const{ // EVALUATES UNNORMALIZED CHI-SQUARED - Float_t chi2 = 0.; + Double_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; + if( TMath::Abs( p1 ) > 1.0E-6 ) + if ( TMath::Abs( delta/p1 ) < 1.0E-4 ) delta = p1/1000; + else delta = (Double_t)1.0E-4; // EVALUATE CHI-SQUARED AT FIRST TWO SEARCH POINTS PeakFunc( xdim, zdim, param, speFit ); chisq1 = ChiSqr( xdim, zdim, spe, speFit ); @@ -553,16 +520,16 @@ void AliITSClusterFinderSDD::Minim( Int_t xdim, Int_t zdim, Float_t *param, if( param[k] <= 20 ) param[k] = fMinPeak; break; case 1 : - if( fabs( param[k] - prm0[k] ) > 1.5 ) param[k] = prm0[k]; + if( TMath::Abs( param[k] - prm0[k] ) > 1.5 ) param[k] = prm0[k]; break; case 2 : - if( fabs( param[k] - prm0[k] ) > 1. ) param[k] = prm0[k]; + if( TMath::Abs( 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] < .288 ) param[k] = .288;// 1/sqrt(12) = 0.288 if( param[k] > zdim*.5 ) param[k] = zdim*.5; break; }; // end switch @@ -578,13 +545,13 @@ void AliITSClusterFinderSDD::Minim( Int_t xdim, Int_t zdim, Float_t *param, // 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); + if( a!=0 ) p0 = (Double_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; + //dp = (Double_t) max (TMath::Abs(p3-p2), TMath::Abs(p2-p1)); + //if( TMath::Abs( p2-p0 ) > dp ) p0 = p2; param[k] = p0; // Constrain paramiters Int_t kpos = (k-1) % knParam; @@ -593,10 +560,10 @@ void AliITSClusterFinderSDD::Minim( Int_t xdim, Int_t zdim, Float_t *param, if( param[k] <= 20 ) param[k] = fMinPeak; break; case 1 : - if( fabs( param[k] - prm0[k] ) > 1.5 ) param[k] = prm0[k]; + if( TMath::Abs( param[k] - prm0[k] ) > 1.5 ) param[k] = prm0[k]; break; case 2 : - if( fabs( param[k] - prm0[k] ) > 1. ) param[k] = prm0[k]; + if( TMath::Abs( param[k] - prm0[k] ) > 1. ) param[k] = prm0[k]; break; case 3 : if( param[k] < .5 ) param[k] = .5; @@ -622,27 +589,27 @@ void AliITSClusterFinderSDD::Minim( Int_t xdim, Int_t zdim, Float_t *param, } //_________________________________________________________________________ Int_t AliITSClusterFinderSDD::NoLinearFit( Int_t xdim, Int_t zdim, - Float_t *param, Float_t *spe, - Int_t *niter, Float_t *chir ){ + Double_t *param, Double_t *spe, + Int_t *niter, Double_t *chir ){ // fit method from Comput. Phys. Commun 46(1987) 149 - const Float_t kchilmt = 0.01; // relative accuracy + const Double_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; + Double_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 ]; + Double_t *prm0 = new Double_t[npeak*knParam+1]; + Double_t *step = new Double_t[npeak*knParam+1]; + Double_t *schi = new Double_t[npeak*knParam+1]; + Double_t *sprm[3]; + sprm[0] = new Double_t[npeak*knParam+1]; + sprm[1] = new Double_t[npeak*knParam+1]; + sprm[2] = new Double_t[npeak*knParam+1]; + Double_t chi0, chi1, reldif, a, b, prmin, dp; + Double_t *speFit = new Double_t[ xdim*zdim ]; PeakFunc( xdim, zdim, param, speFit ); chi0 = ChiSqr( xdim, zdim, spe, speFit ); chi1 = chi0; @@ -659,7 +626,7 @@ Int_t AliITSClusterFinderSDD::NoLinearFit( Int_t xdim, Int_t zdim, iterNum++; chi0 = chi1; Minim( xdim, zdim, param, prm0, step, &chi1, spe, speFit ); - reldif = ( chi1 > 0 ) ? ((Float_t) fabs( chi1-chi0)/chi1 ) : 0; + reldif = ( chi1 > 0 ) ? ((Double_t) TMath::Abs( chi1-chi0)/chi1 ) : 0; // EXIT conditions if( reldif < (float) kchilmt ){ *chir = (chi1>0) ? (float) TMath::Sqrt (chi1/degFree) :0; @@ -688,9 +655,9 @@ Int_t AliITSClusterFinderSDD::NoLinearFit( Int_t xdim, Int_t zdim, // -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]; + Double_t tmp0 = sprm[0][k]; + Double_t tmp1 = sprm[1][k]; + Double_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); @@ -699,7 +666,7 @@ Int_t AliITSClusterFinderSDD::NoLinearFit( Int_t xdim, Int_t zdim, 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; + if( TMath::Abs(prmin-tmp2) > TMath::Abs(dp) ) prmin = tmp2+dp; param[k] = prmin; step[k] = dp/10; // OPTIMIZE SEARCH STEP } // end for k @@ -715,28 +682,27 @@ Int_t AliITSClusterFinderSDD::NoLinearFit( Int_t xdim, Int_t zdim, } //______________________________________________________________________ -void AliITSClusterFinderSDD::ResolveClustersE(){ +void AliITSClusterFinderSDD::ResolveClusters(){ // 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(); + Int_t nofClusters = NClusters(); nofClusters -= fNclusters; - Int_t fNofMaps = fSegmentation->Npz(); + Int_t fNofMaps = GetSeg()->Npz(); Int_t fNofAnodes = fNofMaps/2; - Int_t fMaxNofSamples = fSegmentation->Npx(); + //Int_t fMaxNofSamples = GetSeg()->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 + Double_t fTimeStep = GetSeg()->Dpx( dummy ); + Double_t fSddLength = GetSeg()->Dx(); + Double_t fDriftSpeed = GetResp(fModule)->GetDriftSpeed(); + Double_t anodePitch = GetSeg()->Dpz( dummy ); + Double_t n, baseline; + GetResp(fModule)->GetNoiseParam( n, baseline ); + Int_t electronics =GetResp(fModule)->GetElectronics(); // 1 = PASCAL, 2 = OLA for( Int_t j=0; jAt(j); + AliITSRawClusterSDD *clusterJ=(AliITSRawClusterSDD*) Cluster(j); Int_t astart = clusterJ->Astart(); Int_t astop = clusterJ->Astop(); Int_t tstart = clusterJ->Tstartf(); @@ -748,14 +714,17 @@ void AliITSClusterFinderSDD::ResolveClustersE(){ } // 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) ); + if( xdim > 50 || zdim > 30 ) { + Warning("ResolveClusters","xdim: %d , zdim: %d ",xdim,zdim); + continue; + } + Double_t *sp = new Double_t[ xdim*zdim+1 ]; + memset( sp, 0, sizeof(Double_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 ); + Double_t fadc = Map()->GetSignal( ianode, itime ); if( fadc > baseline ) fadc -= (Double_t)baseline; else fadc = 0.; Int_t index = (itime-tstart+1)*zdim+(ianode-astart+1); @@ -767,64 +736,64 @@ void AliITSClusterFinderSDD::ResolveClustersE(){ const Int_t kNp = 150; Int_t peakX1[kNp]; Int_t peakZ1[kNp]; - Float_t peakAmp1[kNp]; + Double_t peakAmp1[kNp]; Int_t npeak = SearchPeak(sp,xdim,zdim,peakX1,peakZ1,peakAmp1,fMinPeak); // if multiple peaks, split cluster - if( npeak >= 1 ) - { + if( npeak >= 1 ){ // cout << "npeak " << npeak << endl; // clusterJ->PrintInfo(); - Float_t *par = new Float_t[npeak*5+1]; - par[0] = (Float_t)npeak; + Double_t *par = new Double_t[npeak*5+1]; + par[0] = (Double_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 ); + + // Int_t newiTime = peakX1[i]-1 + tstart; + // Int_t shift = (Int_t)(fTimeCorr/fTimeStep + 0.5); + // if( newiTime > shift && newiTime < (fMaxNofSamples-shift) ) + // shift = 0; + // Int_t peakpos = Map()->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; + Double_t newAnodef = peakZ[i] - 0.5 + astart; + Double_t newiTimef = peakX[i] - 1 + tstart; if( wing == 2 ) newAnodef -= fNofAnodes; - Float_t anodePath = (newAnodef - fNofAnodes/2)*anodePitch; + Double_t anodePath = (newAnodef - fNofAnodes/2)*anodePitch; newiTimef *= fTimeStep; if( newiTimef > fTimeCorr ) newiTimef -= fTimeCorr; if( electronics == 1 ){ @@ -832,8 +801,27 @@ void AliITSClusterFinderSDD::ResolveClustersE(){ // 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.; + + Int_t timeBin = (Int_t)(newiTimef/fTimeStep+0.5); + Int_t peakpos = Map()->GetHitIndex( newAnode, timeBin ); + if( peakpos < 0 ) { + for( Int_t ii=0; ii<3; ii++ ) { + peakpos = Map()->GetHitIndex( newAnode, timeBin+ii ); + if( peakpos > 0 ) break; + peakpos = Map()->GetHitIndex( newAnode, timeBin-ii ); + if( peakpos > 0 ) break; + } + } + + if( peakpos < 0 ) { + //Warning("ResolveClusters", + // "Digit not found for cluster"); + //if(AliDebugLevel()>=3) clusterI.PrintInfo(); + continue; + } + clusterI.SetPeakPos( peakpos ); + Double_t driftPath = fSddLength - newiTimef * fDriftSpeed; + Double_t sign = ( wing == 1 ) ? -1. : 1.; clusterI.SetX( driftPath*sign * 0.0001 ); clusterI.SetZ( anodePath * 0.0001 ); clusterI.SetAnode( newAnodef ); @@ -841,32 +829,31 @@ void AliITSClusterFinderSDD::ResolveClustersE(){ clusterI.SetAsigma( sigma[i]*anodePitch ); clusterI.SetTsigma( tau[i]*fTimeStep ); clusterI.SetQ( integral[i] ); - // clusterI.PrintInfo(); - iTS->AddCluster( 1, &clusterI ); + + fDetTypeRec->AddCluster( 1, &clusterI ); } // end for i - fClusters->RemoveAt( j ); + Clusters()->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; + Warning( "ResolveClusters", + "--- Peak not found!!!! minpeak=%d ,cluster peak= %f" + " , module= %d", + fMinPeak, clusterJ->PeakAmpl(),GetModule()); + clusterJ->PrintInfo(); + Warning( "ResolveClusters"," xdim= %d zdim= %d", xdim-2, zdim-2 ); } delete [] sp; } // cluster loop - fClusters->Compress(); -// fMap->ClearMap(); + Clusters()->Compress(); +// Map()->ClearMap(); } - - //________________________________________________________________________ void AliITSClusterFinderSDD::GroupClusters(){ // group clusters Int_t dummy=0; - Float_t fTimeStep = fSegmentation->Dpx(dummy); + Double_t fTimeStep = GetSeg()->Dpx(dummy); // get number of clusters for this module - Int_t nofClusters = fClusters->GetEntriesFast(); + Int_t nofClusters = NClusters(); nofClusters -= fNclusters; AliITSRawClusterSDD *clusterI; AliITSRawClusterSDD *clusterJ; @@ -877,23 +864,25 @@ void AliITSClusterFinderSDD::GroupClusters(){ if(label[i] != 0) continue; for(j=i+1; jAt(i); - clusterJ = (AliITSRawClusterSDD*) fClusters->At(j); + clusterI = (AliITSRawClusterSDD*) Cluster(i); + clusterJ = (AliITSRawClusterSDD*) Cluster(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(); + if(AliDebugLevel()>=4){ + clusterI->PrintInfo(); + clusterJ->PrintInfo(); + } // end if AliDebugLevel clusterI->Add(clusterJ); label[j] = 1; - fClusters->RemoveAt(j); + Clusters()->RemoveAt(j); j=i; // <- Ernesto } // J clusters label[i] = 1; } // I clusters - fClusters->Compress(); + Clusters()->Compress(); delete [] label; return; @@ -901,16 +890,16 @@ void AliITSClusterFinderSDD::GroupClusters(){ //________________________________________________________________________ void AliITSClusterFinderSDD::SelectClusters(){ // get number of clusters for this module - Int_t nofClusters = fClusters->GetEntriesFast(); + Int_t nofClusters = NClusters(); nofClusters -= fNclusters; Int_t i; for(i=0; iAt(i); + AliITSRawClusterSDD *clusterI =(AliITSRawClusterSDD*) Cluster(i); Int_t rmflg = 0; - Float_t wy = 0.; + Double_t wy = 0.; if(clusterI->Anodes() != 0.) { - wy = ((Float_t) clusterI->Samples())/clusterI->Anodes(); + wy = ((Double_t) clusterI->Samples())/clusterI->Anodes(); } // end if Int_t amp = (Int_t) clusterI->PeakAmpl(); Int_t cha = (Int_t) clusterI->Q(); @@ -918,377 +907,79 @@ void AliITSClusterFinderSDD::SelectClusters(){ if(cha < fMinCharge) rmflg = 1; if(wy < fMinNCells) rmflg = 1; //if(wy > fMaxNCells) rmflg = 1; - if(rmflg) fClusters->RemoveAt(i); + if(rmflg) Clusters()->RemoveAt(i); } // I clusters - fClusters->Compress(); + Clusters()->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(); + Int_t nofClusters = NClusters(); 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,j; + const Double_t kconvGeV = 1.e-6; // GeV -> KeV + const Double_t kconv = 1.0e-4; + const Double_t kRMSx = 38.0*kconv; // microns->cm ITS TDR Table 1.3 + const Double_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(); + Int_t ndigits=NDigits(); + + Int_t lay,lad,det; + fDetTypeRec->GetITSgeom()->GetModuleId(fModule,lay,lad,det); + Int_t ind=(lad-1)*fDetTypeRec->GetITSgeom()->GetNdetectors(lay)+(det-1); + Int_t lyr=(lay-1); + + for(i=0; iAt(i); + AliITSRawClusterSDD *clusterI = (AliITSRawClusterSDD*)Cluster(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(idx&&idx<= ndigits) dig =(AliITSdigitSDD*)GetDigit(idx); if(!dig) { // try cog - fSegmentation->GetPadIxz(clusterI->X(),clusterI->Z(),ix,iz); - dig = (AliITSdigitSDD*)fMap->GetHit(iz-1,ix-1); + GetSeg()->GetPadIxz(clusterI->X(),clusterI->Z(),ix,iz); + dig = (AliITSdigitSDD*)Map()->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) dig = (AliITSdigitSDD*)Map()->GetHit(iz-1,ix); + if (!dig) dig = (AliITSdigitSDD*)Map()->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()); + AliITSRecPoint rnew(fDetTypeRec->GetITSgeom()); + rnew.SetXZ(fModule,clusterI->X(),clusterI->Z()); rnew.SetQ(clusterI->Q()); // in KeV - should be ADC rnew.SetdEdX(kconvGeV*clusterI->Q()); - rnew.SetSigmaX2(kRMSx*kRMSx); + rnew.SetSigmaDetLocX2(kRMSx*kRMSx); rnew.SetSigmaZ2(kRMSz*kRMSz); - if(dig){ - rnew.fTracks[0] = dig->fTracks[0]; - rnew.fTracks[1] = -3; - rnew.fTracks[2] = -3; - j=1; - while(rnew.fTracks[0]==dig->fTracks[j] && - jGetNTracks()) j++; - if(jGetNTracks()){ - rnew.fTracks[1] = dig->fTracks[j]; - while(rnew.fTracks[1]==dig->fTracks[j] && - jGetNTracks()) j++; - if(jGetNTracks()) rnew.fTracks[2] = dig->fTracks[j]; - } // end if - } // end if - //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); + + if(dig) rnew.SetLabel(dig->GetTrack(0),0); + if(dig) rnew.SetLabel(dig->GetTrack(1),1); + if(dig) rnew.SetLabel(dig->GetTrack(2),2); + rnew.SetDetectorIndex(ind); + rnew.SetLayer(lyr); + fDetTypeRec->AddRecPoint(rnew); } // I clusters -// fMap->ClearMap(); +// Map()->ClearMap(); } //______________________________________________________________________ void AliITSClusterFinderSDD::FindRawClusters(Int_t mod){ // find raw clusters - fModule = mod; - + SetModule(mod); Find1DClustersE(); GroupClusters(); SelectClusters(); - ResolveClustersE(); + ResolveClusters(); GetRecPoints(); } //_______________________________________________________________________ -void AliITSClusterFinderSDD::Print(){ +void AliITSClusterFinderSDD::PrintStatus() const{ // Print SDD cluster finder Parameters cout << "**************************************************" << endl;