X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=ITS%2FAliITSClusterFinderSDD.cxx;h=41115ba683c95faa734a3e0a610b338a4a0708a7;hb=4178b5c7ed8261863bd588cdd933cc6a20452437;hp=f0085137fb96955294af8f1190fc963bf64c25c0;hpb=9e1a0ddb8feb3fb19ca04ee19673f48c265a8fe6;p=u%2Fmrichter%2FAliRoot.git diff --git a/ITS/AliITSClusterFinderSDD.cxx b/ITS/AliITSClusterFinderSDD.cxx index f0085137fb9..41115ba683c 100644 --- a/ITS/AliITSClusterFinderSDD.cxx +++ b/ITS/AliITSClusterFinderSDD.cxx @@ -12,1484 +12,1041 @@ * about the suitability of this software for any purpose. It is * * provided "as is" without express or implied warranty. * **************************************************************************/ +/* + $Id$ + */ +/////////////////////////////////////////////////////////////////////////// +// 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 "AliITSresponse.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, TClonesArray *digits, TClonesArray *recp) -{ - // constructor - - fSegmentation=seg; - fResponse=response; - fDigits=digits; - fClusters=recp; - fNclusters= fClusters->GetEntriesFast(); - SetCutAmplitude(); - SetDAnode(); - SetDTime(); - SetMinPeak(); - SetMinNCells(); - SetMaxNCells(); - SetTimeCorr(); - SetMinCharge(); - fMap=new AliITSMapA1(fSegmentation,fDigits,fCutAmplitude); - +//______________________________________________________________________ +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() -{ - // constructor - fSegmentation=0; - fResponse=0; - fDigits=0; - fClusters=0; - fNclusters=0; - fMap=0; - SetCutAmplitude(); +//______________________________________________________________________ +AliITSClusterFinderSDD::AliITSClusterFinderSDD(AliITSDetTypeRec* dettyp, + TClonesArray *digits, + 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 + + SetDigits(digits); + SetClusters(recp); + SetCutAmplitude(fDetTypeRec->GetITSgeom()->GetStartSDD()); SetDAnode(); SetDTime(); - SetMinPeak(); + SetMinPeak((Int_t)((AliITSCalibrationSDD*)GetResp(fDetTypeRec->GetITSgeom()->GetStartSDD()))->GetNoiseAfterElectronics(0)*5); SetMinNCells(); SetMaxNCells(); SetTimeCorr(); SetMinCharge(); - + SetMap(new AliITSMapA1(GetSeg(),Digits(),fCutAmplitude)); } - -//_____________________________________________________________________________ -AliITSClusterFinderSDD::~AliITSClusterFinderSDD() -{ - // destructor - - if(fMap) delete fMap; - -} -//__________________________________________________________________________ -AliITSClusterFinderSDD::AliITSClusterFinderSDD(const AliITSClusterFinderSDD &source){ - // Copy Constructor - if(&source == this) return; - this->fClusters = source.fClusters ; - this->fNclusters = source.fNclusters ; - this->fMap = source.fMap ; - this->fCutAmplitude = source.fCutAmplitude ; - this->fDAnode = source.fDAnode ; - this->fDTime = source.fDTime ; - this->fTimeCorr = source.fTimeCorr ; - this->fMinPeak = source.fMinPeak ; - this->fMinNCells = source.fMinNCells ; - this->fMaxNCells = source.fMaxNCells ; - return; -} - -//_________________________________________________________________________ -AliITSClusterFinderSDD& - AliITSClusterFinderSDD::operator=(const AliITSClusterFinderSDD &source) { - // Assignment operator - if(&source == this) return *this; - this->fClusters = source.fClusters ; - this->fNclusters = source.fNclusters ; - this->fMap = source.fMap ; - this->fCutAmplitude = source.fCutAmplitude ; - this->fDAnode = source.fDAnode ; - this->fDTime = source.fDTime ; - this->fTimeCorr = source.fTimeCorr ; - this->fMinPeak = source.fMinPeak ; - this->fMinNCells = source.fMinNCells ; - this->fMaxNCells = source.fMaxNCells ; - return *this; +//______________________________________________________________________ +void AliITSClusterFinderSDD::SetCutAmplitude(Int_t mod,Double_t nsigma){ + // set the signal threshold for cluster finder + Double_t baseline,noiseAfterEl; + + AliITSresponseSDD* res = (AliITSresponseSDD*)((AliITSCalibrationSDD*)GetResp(mod))->GetResponse(); + const char *option=res->ZeroSuppOption(); + Int_t nanodes = GetResp(mod)->Wings()*GetResp(mod)->Channels()*GetResp(mod)->Chips(); + fCutAmplitude.Set(nanodes); + for(Int_t ian=0;ianGetNoiseAfterElectronics(ian); + if((strstr(option,"1D")) || (strstr(option,"2D"))){ + fCutAmplitude[ian] = (Int_t)(nsigma*noiseAfterEl); + } + else{ + baseline=GetResp(mod)->GetBaseline(ian); + fCutAmplitude[ian] = (Int_t)((baseline + nsigma*noiseAfterEl)); + } + } } - - -//_____________________________________________________________________________ - -void AliITSClusterFinderSDD::Find1DClusters() -{ - // find 1D clusters - - AliITS *iTS=(AliITS*)gAlice->GetModule("ITS"); +//______________________________________________________________________ +void AliITSClusterFinderSDD::Find1DClusters(){ + // find 1D clusters - // 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->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++; - } - - if(dfadc[k][id] > dfadcmax) { - dfadcmax = dfadc[k][id]; - imaxd = id; - } - } - 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]; - } - } - - Float_t clusterCharge = 0.; - Float_t clusterAnode = k+0.5; - Float_t clusterTime = 0.; - Float_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); - } - clusterTime += fadc*its; - if(fadc > 0) clusterMult++; - if(its == tstop) { - clusterTime /= (clusterCharge/fTimeStep); // ns - if(clusterTime > fTimeCorr) clusterTime -= fTimeCorr; // ns - } - } - - Float_t clusteranodePath = (clusterAnode - fNofAnodes/2)*anodePitch; - Float_t clusterDriftPath = clusterTime*fDriftSpeed; - clusterDriftPath = fSddLength-clusterDriftPath; - - if(clusterCharge <= 0.) break; - AliITSRawClusterSDD clust(j+1,clusterAnode,clusterTime,clusterCharge,clusterPeakAmplitude,peakpos,0.,0.,clusterDriftPath,clusteranodePath,clusterMult,0,0,0,0,0,0,0); - iTS->AddCluster(1,&clust); - it = tstop; - } // ilcl - - it++; - - } // while (samples) - } // anodes - } // detectors (2) - - //fMap->ClearMap(); - - for(i=0;iNpz(); + Int_t fMaxNofSamples = GetSeg()->Npx(); + Int_t fNofAnodes = fNofMaps/2; + Int_t dummy = 0; + Double_t fTimeStep = GetSeg()->Dpx(dummy); + Double_t fSddLength = GetSeg()->Dx(); + Double_t anodePitch = GetSeg()->Dpz(dummy); + AliITSCalibrationSDD* cal = (AliITSCalibrationSDD*)GetResp(fModule); + AliITSresponseSDD* res = (AliITSresponseSDD*)((AliITSCalibrationSDD*)GetResp(fModule))->GetResponse(); + const char *option=res->ZeroSuppOption(); + + // map the signal + Map()->ClearMap(); + Map()->SetThresholdArr(fCutAmplitude); + Map()->FillMap2(); + Int_t nofFoundClusters = 0; + Int_t i; + Double_t **dfadc = new Double_t*[fNofAnodes]; + for(i=0;iGetSignal(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)Map()->GetSignal(idx,id); + if(fadc > fadcmax) { fadcmax = fadc; imax = id;} + if(fadc > (float)fCutAmplitude[idx])lthrt++; + if(dfadc[k][id] > dfadcmax) { + dfadcmax = dfadc[k][id]; + imaxd = id; + } // end if + } // end for m + it = imaxd; + if(Map()->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; + Double_t dfadcmin = 10000.; + Int_t ij; + for(ij=0; ij<20; ij++) { + if(tstart+ij > 255) { tstop = 255; break; } + fadc=(float)Map()->GetSignal(idx,tstart+ij); + if((dfadc[k][tstart+ij] < dfadcmin) && + (fadc > fCutAmplitude[idx])) { + tstop = tstart+ij+5; + if(tstop > 255) tstop = 255; + dfadcmin = dfadc[k][it+ij]; + } // end if + } // end for ij + + Double_t clusterCharge = 0.; + Double_t clusterAnode = k+0.5; + Double_t clusterTime = 0.; + Int_t clusterMult = 0; + Double_t clusterPeakAmplitude = 0.; + Int_t its,peakpos = -1; + + for(its=tstart; its<=tstop; its++) { + fadc=(float)Map()->GetSignal(idx,its); + if(!((strstr(option,"1D")) || (strstr(option,"2D")))){ + Double_t baseline=GetResp(fModule)->GetBaseline(idx); + 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=Map()->GetHitIndex(idx,its); + Int_t shift = (int)(fTimeCorr/fTimeStep); + if(its>shift && its<(fMaxNofSamples-shift)) + 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++; + if(its == tstop) { + clusterTime /= (clusterCharge/fTimeStep); // ns + if(clusterTime>fTimeCorr) clusterTime -=fTimeCorr; + //ns + } // end if + } // end for its + + Double_t clusteranodePath = (clusterAnode - fNofAnodes/2)* + anodePitch; + Double_t clusterDriftPath = (Double_t)cal->GetDriftPath(clusterTime,clusteranodePath); + 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 + fDetTypeRec->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->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; - } - - 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 ); - } - } - else - { - 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++; - } - On = kFALSE; - } - } - } // samples - } // anodes - } // wings - -// cout << "# Rawclusters " << nClu << endl; - return; - +//______________________________________________________________________ +void AliITSClusterFinderSDD::Find1DClustersE(){ + // find 1D clusters + // retrieve the parameters + Int_t fNofMaps = GetSeg()->Npz(); + Int_t fMaxNofSamples = GetSeg()->Npx(); + Int_t fNofAnodes = fNofMaps/2; + Int_t dummy=0; + Double_t fTimeStep = GetSeg()->Dpx( dummy ); + Double_t fSddLength = GetSeg()->Dx(); + Double_t anodePitch = GetSeg()->Dpz( dummy ); + AliITSCalibrationSDD* cal = (AliITSCalibrationSDD*)GetResp(fModule); + Map()->ClearMap(); + Map()->SetThresholdArr( fCutAmplitude ); + Map()->FillMap2(); + + AliITSresponseSDD* res = (AliITSresponseSDD*)cal->GetResponse(); + const char *option=res->ZeroSuppOption(); + + 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(!((strstr(option,"1D")) || (strstr(option,"2D")))){ + Double_t baseline=GetResp(fModule)->GetBaseline(idx); + 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 = Map()->GetHitIndex( idx, l+shift ); + else + peakpos = Map()->GetHitIndex( idx, l ); + if( peakpos < 0) peakpos = Map()->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 + Double_t anodePath =(anode-fNofAnodes/2)*anodePitch; + + Double_t driftPath = (Double_t)cal->GetDriftPath(time,anodePath); + driftPath = fSddLength-driftPath; + AliITSRawClusterSDD clust(j+1,anode,time,charge, + fmax, peakpos,0.,0., + driftPath,anodePath, + nTsteps,start,stop, + start, stop, 1, k, k ); + fDetTypeRec->AddCluster( 1, &clust ); + if(AliDebugLevel()>=5) clust.PrintInfo(); + nClu++; + } // end if nTsteps + on = kFALSE; + } // end if on==kTRUE + } // end if fadc>0 + } // samples + } // anodes + } // wings + AliDebug(3,Form("# Rawclusters %d",nClu)); + 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 ) -{ - Int_t npeak = 0; // # peaks +//_______________________________________________________________________ +Int_t AliITSClusterFinderSDD::SearchPeak(Double_t *spect,Int_t xdim,Int_t zdim, + Int_t *peakX, Int_t *peakZ, + 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 ] && - 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++; - } - } - } - - // 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]; - } - } - } - - // make average of peak groups - for( i=0; i 1 ) - { - peakX[i] /= nFlag; - peakZ[i] /= nFlag; - } - } - - delete [] Flag; - return( npeak ); -} - - -void AliITSClusterFinderSDD::PeakFunc( Int_t xdim, Int_t zdim, Float_t *par, Float_t *spe, Float_t *Integral) -{ - Int_t Electronics = fResponse->Electronics(); // 1 = PASCAL, 2 = OLA - Int_t param_peak = 5; - // par -> paramiters.. - // 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 npeak = (Int_t)par[0]; - - memset( spe, 0, sizeof( Float_t )*zdim*xdim ); - - Int_t k = 1; - for( Int_t i=0; i 0.) ? par[k] * x2*x2 * exp( -2*x2+2. - z2 ) : 0.0; // RCCR - Float_t signal = 0.; - if(Electronics == 1) - signal = (x2 > 0.) ? par[k] * x2 * exp( -x2+1. - z2 ) : 0.0; - else if(Electronics == 2) //OLA - signal = par[k] * exp( -x2 - z2 ); - else - cout << "Wrong Electronics" << endl; - - spe[x*zdim+z] += signal; - if( Integral != 0 ) Integral[i] += signal; - } - } - k += param_peak; - } - return; + // 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=0 ) -{ - Int_t param_peak = 5; -// par -> paramiters.. -// 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 npeak = (Int_t)par[0]; - - memset( spe, 0, sizeof( Float_t )*zdim*xdim ); - - Int_t k = 1; - for( Int_t i=0; i 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 = GetResp(fModule)->GetElectronics(); // 1 = PASCAL, 2 = OLA + const Int_t knParam = 5; + Int_t npeak = (Int_t)par[0]; + + 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; + signal = par[k] * exp( -x2 - z2 ); + } else { + Warning("PeakFunc","Wrong SDD Electronics = %d", + electronics); + // 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::chisq( 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 = chisq( xdim, zdim, spe, speFit ); - - p2 = p1+delta; - param[k] = p2; - - PeakFunc( xdim, zdim, param, speFit ); - chisq2 = chisq( 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; - } - - 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; - } - param[k] = p3; - - // Constrain paramiters - Int_t kpos = (k-1) % param_peak; - switch( kpos ) - { - case 0 : - if( param[k] <= 20 ) param[k] = fMinPeak; - case 1 : - if( fabs( param[k] - prm0[k] ) > 1.5 ) param[k] = prm0[k]; - case 2 : - if( fabs( param[k] - prm0[k] ) > 1. ) param[k] = prm0[k]; - case 3 : - if( param[k] < .5 ) param[k] = .5; - case 4 : - if( param[k] < .288 ) param[k] = .288; // 1/sqrt(12) = 0.288 - }; - - PeakFunc( xdim, zdim, param, speFit ); - chisq3 = chisq( 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) % param_peak; - switch( kpos ) - { - case 0 : - if( param[k] <= 20 ) param[k] = fMinPeak; - case 1 : - if( fabs( param[k] - prm0[k] ) > 1.5 ) param[k] = prm0[k]; - case 2 : - if( fabs( param[k] - prm0[k] ) > 1. ) param[k] = prm0[k]; - case 3 : - if( param[k] < .5 ) param[k] = .5; - case 4 : - if( param[k] < .288 ) param[k] = .288; // 1/sqrt(12) = 0.288 - }; - - PeakFunc( xdim, zdim, param, speFit ); - chisqt = chisq( 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; - } - - // EVALUATE FIT AND CHI-SQUARED FOR OPTIMIZED PARAMETERS - PeakFunc( xdim, zdim, param, speFit ); - *chisqr = chisq( xdim, zdim, spe, speFit ); - return; +//_______________________________________________________________________ +void AliITSClusterFinderSDD::Minim( Int_t xdim, Int_t zdim, Double_t *param, + Double_t *prm0,Double_t *steprm, + Double_t *chisqr,Double_t *spe, + Double_t *speFit ){ + // + Int_t k, nnn, mmm, i; + Double_t p1, delta, d1, chisq1, p2, chisq2, t, p3, chisq3, a, b, p0, chisqt; + const Int_t knParam = 5; + Int_t npeak = (Int_t)param[0]; + for( k=1; k<(npeak*knParam+1); k++ ) prm0[k] = param[k]; + for( k=1; k<(npeak*knParam+1); k++ ){ + p1 = param[k]; + delta = steprm[k]; + d1 = delta; + // ENSURE THAT STEP SIZE IS SENSIBLY LARGER THAN MACHINE ROUND OFF + 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 ); + 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( TMath::Abs( param[k] - prm0[k] ) > 1.5 ) param[k] = prm0[k]; + break; + case 2 : + 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] > 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 = (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 = (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; + switch( kpos ){ + case 0 : + if( param[k] <= 20 ) param[k] = fMinPeak; + break; + case 1 : + if( TMath::Abs( param[k] - prm0[k] ) > 1.5 ) param[k] = prm0[k]; + break; + case 2 : + 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] > 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 ) -{ - const Float_t chilmt = 0.01; // relative accuracy - const Int_t nel = 3; // for parabolic minimization - const Int_t nstop = 50; // Max. iteration number - const Int_t param_peak = 5; - - Int_t npeak = (Int_t)param[0]; - - // RETURN IF NUMBER OF DEGREES OF FREEDOM IS NOT POSITIVE - if( (xdim*zdim - npeak*param_peak) <= 0 ) return( -1 ); - Float_t deg_free = (xdim*zdim - npeak*param_peak)-1; - - Int_t n, k, iter_num = 0; - Float_t *prm0 = new Float_t[npeak*param_peak+1]; - Float_t *step = new Float_t[npeak*param_peak+1]; - Float_t *schi = new Float_t[npeak*param_peak+1]; - Float_t *sprm[3]; - sprm[0] = new Float_t[npeak*param_peak+1]; - sprm[1] = new Float_t[npeak*param_peak+1]; - sprm[2] = new Float_t[npeak*param_peak+1]; - - Float_t chi0, chi1, reldif, a, b, prmin, dp; - - Float_t *speFit = new Float_t[ xdim*zdim ]; - PeakFunc( xdim, zdim, param, speFit ); - chi0 = chisq( xdim, zdim, spe, speFit ); - chi1 = chi0; - - - for( k=1; k<(npeak*param_peak+1); k++) prm0[k] = param[k]; - - for( k=1 ; k<(npeak*param_peak+1); k+=param_peak ) - { - 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; - } - - Int_t out = 0; - do - { - iter_num++; - 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) chilmt ) - { - *chir = (chi1>0) ? (float) TMath::Sqrt (chi1/deg_free) :0; - *niter = iter_num; - out = 0; - break; - } - - if( (reldif < (float)(5*chilmt)) && (iter_num > nstop) ) - { - *chir = (chi1>0) ?(float) TMath::Sqrt (chi1/deg_free):0; - *niter = iter_num; - out = 0; - break; - } - - if( iter_num > 5*nstop ) - { - *chir = (chi1>0) ?(float) TMath::Sqrt (chi1/deg_free):0; - *niter = iter_num; - out = 1; - break; - } - - if( iter_num <= nel ) continue; - - n = iter_num - (iter_num/nel)*nel; // EXTRAPOLATION LIMIT COUNTER N - if( n > 3 || n == 0 ) continue; - schi[n-1] = chi1; - for( k=1; k<(npeak*param_peak+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*param_peak+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 - } - - } while( kTRUE ); - - delete [] prm0; - delete [] step; - delete [] schi; - delete [] sprm[0]; - delete [] sprm[1]; - delete [] sprm[2]; - delete [] speFit; - - return( out ); +//_________________________________________________________________________ +Int_t AliITSClusterFinderSDD::NoLinearFit( Int_t xdim, Int_t zdim, + Double_t *param, Double_t *spe, + Int_t *niter, Double_t *chir ){ + // fit method from Comput. Phys. Commun 46(1987) 149 + 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 ); + Double_t degFree = (xdim*zdim - npeak*knParam)-1; + Int_t n, k, iterNum = 0; + 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; + 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 ) ? ((Double_t) TMath::Abs( 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++ ){ + 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); + 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( TMath::Abs(prmin-tmp2) > TMath::Abs(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; - - 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 - - // fill Map of signals - fMap->FillMap(); - - 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; - } - Int_t xdim = tstop-tstart+3; - Int_t zdim = astop-astart+3; - 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 np = 150; - Int_t peakX1[np]; - Int_t peakZ1[np]; - Float_t peakAmp1[np]; - 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 paramiters 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 +//______________________________________________________________________ +void AliITSClusterFinderSDD::ResolveClusters(){ + // The function to resolve clusters if the clusters overlapping exists + Int_t i; + // get number of clusters for this module + Int_t nofClusters = NClusters(); + nofClusters -= fNclusters; + Int_t fNofMaps = GetSeg()->Npz(); + Int_t fNofAnodes = fNofMaps/2; + //Int_t fMaxNofSamples = GetSeg()->Npx(); + Int_t dummy=0; + Double_t fTimeStep = GetSeg()->Dpx( dummy ); + Double_t fSddLength = GetSeg()->Dx(); + Double_t anodePitch = GetSeg()->Dpz( dummy ); + Int_t electronics =GetResp(fModule)->GetElectronics(); // 1 = PASCAL, 2 = OLA + AliITSCalibrationSDD* cal = (AliITSCalibrationSDD*)GetResp(fModule); + AliITSresponseSDD* res = (AliITSresponseSDD*)cal->GetResponse(); + const char *option=res->ZeroSuppOption(); + - 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 ); - } - fClusters->RemoveAt( j ); - delete [] par; + for( Int_t j=0; jAstart(); + 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 ) { + 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++ ){ + Double_t fadc = Map()->GetSignal( ianode, itime ); + if(!((strstr(option,"1D")) || (strstr(option,"2D")))){ + Double_t baseline=GetResp(fModule)->GetBaseline(ianode); + if( fadc > baseline ) fadc -= (Double_t)baseline; + else fadc = 0.; } - else cout <<" --- Peak not found!!!! minpeak=" << fMinPeak<< - " cluster peak=" << clusterJ->PeakAmpl() << endl << endl; - - delete [] sp; - } // cluster loop - - fClusters->Compress(); - fMap->ClearMap(); + 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]; + Double_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(); + 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 = Map()->GetHitIndex(newAnode,newiTime+shift ); + // clusterI.SetPeakPos( peakpos ); + + clusterI.SetPeakAmpl( peakAmp1[i] ); + Double_t newAnodef = peakZ[i] - 0.5 + astart; + Double_t newiTimef = peakX[i] - 1 + tstart; + if( wing == 2 ) newAnodef -= fNofAnodes; + Double_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 + + 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 ); + Float_t dp = cal->GetDriftPath(newiTimef,anodePath); + Double_t driftPath = fSddLength - (Double_t)dp; + Double_t sign = ( wing == 1 ) ? -1. : 1.; + Double_t xcoord = driftPath*sign * 0.0001; + Double_t zcoord = anodePath * 0.0001; + CorrectPosition(zcoord,xcoord); + clusterI.SetX( xcoord ); + clusterI.SetZ( zcoord ); + clusterI.SetAnode( newAnodef ); + clusterI.SetTime( newiTimef ); + clusterI.SetAsigma( sigma[i]*anodePitch ); + clusterI.SetTsigma( tau[i]*fTimeStep ); + clusterI.SetQ( integral[i] ); + + fDetTypeRec->AddCluster( 1, &clusterI ); + } // end for i + Clusters()->RemoveAt( j ); + delete [] par; + } else { // something odd + 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 + Clusters()->Compress(); +// Map()->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::GroupClusters(){ + // group clusters + Int_t dummy=0; + Double_t fTimeStep = GetSeg()->Dpx(dummy); + // get number of clusters for this module + Int_t nofClusters = NClusters(); + nofClusters -= fNclusters; + AliITSRawClusterSDD *clusterI; + AliITSRawClusterSDD *clusterJ; + Int_t *label = new Int_t [nofClusters]; + Int_t i,j; + for(i=0; iT() < 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; + if(AliDebugLevel()>=4){ + clusterI->PrintInfo(); + clusterJ->PrintInfo(); + } // end if AliDebugLevel + clusterI->Add(clusterJ); + label[j] = 1; + Clusters()->RemoveAt(j); + j=i; // <- Ernesto + } // J clusters + label[i] = 1; + } // I clusters + Clusters()->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(); - } - 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::SelectClusters(){ + // get number of clusters for this module + Int_t nofClusters = NClusters(); + + nofClusters -= fNclusters; + Int_t i; + for(i=0; iAnodes() != 0.) { + wy = ((Double_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) Clusters()->RemoveAt(i); + } // I clusters + 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; - } - -// 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 ! - } - 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 - - 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"); - } + // get number of clusters for this module + Int_t nofClusters = NClusters(); + nofClusters -= fNclusters; + 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=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; iPeakPos(); + if(idx>ndigits) Error("SDD: GetRecPoints","idx ndigits",idx,ndigits); + // try peak neighbours - to be done + if(idx&&idx<= ndigits) dig =(AliITSdigitSDD*)GetDigit(idx); + if(!dig) { + // try cog + GetSeg()->GetPadIxz(clusterI->X(),clusterI->Z(),ix,iz); + dig = (AliITSdigitSDD*)Map()->GetHit(iz-1,ix-1); + // if null try neighbours + 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 + + Int_t lab[4] = {-3141593,-3141593,-3141593,ind}; + if (dig) { + lab[0] = dig->GetTrack(0); + lab[1] = dig->GetTrack(1); + lab[2] = dig->GetTrack(2); + } + Float_t hit[5] = {clusterI->X(),clusterI->Z(),kRMSx*kRMSx,kRMSz*kRMSz,clusterI->Q()}; + Int_t info[3] = {0,0,lyr}; - 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],clusterI->X(),clusterI->Z()); - iTS->AddRecPoint(rnew); - } // I clusters + AliITSRecPoint rnew(lab,hit,info,kTRUE); + rnew.SetdEdX(kconvGeV*clusterI->Q()); - fMap->ClearMap(); + fDetTypeRec->AddRecPoint(rnew); + } // I clusters +// Map()->ClearMap(); } - -//_____________________________________________________________________________ - -void AliITSClusterFinderSDD::FindRawClusters(Int_t mod) -{ - // find raw clusters +//______________________________________________________________________ +void AliITSClusterFinderSDD::FindRawClusters(Int_t mod){ + // find raw clusters + + SetModule(mod); + SetCutAmplitude(mod); + Int_t nanodes=GetSeg()->Npz(); + Int_t noise=0; + for(Int_t i=0;iGetNoiseAfterElectronics(i)); + } + SetMinPeak((noise/nanodes)*5); Find1DClustersE(); GroupClusters(); SelectClusters(); - ResolveClustersE(); + ResolveClusters(); GetRecPoints(); } -//_____________________________________________________________________________ +//_______________________________________________________________________ +void AliITSClusterFinderSDD::PrintStatus() const{ + // 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[0] << 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; +} + +//_________________________________________________________________________ +void AliITSClusterFinderSDD::CorrectPosition(Double_t &z, Double_t&y){ + //correction of coordinates using the maps stored in the DB + + AliITSCalibrationSDD* cal = (AliITSCalibrationSDD*)GetResp(fModule); + static const Int_t nbint = cal->GetMapTimeNBin(); + static const Int_t nbina = cal->Chips()*cal->Channels(); + Float_t stepa = (GetSeg()->Dpz(0))/10000.; //anode pitch in cm + Float_t stept = (GetSeg()->Dx()/cal->GetMapTimeNBin()/2.)/10.; + + Int_t bint = TMath::Abs((Int_t)(y/stept)); + if(y>=0) bint+=(Int_t)(nbint/2.); + if(bint>nbint) AliError("Wrong bin number!"); + + Int_t bina = TMath::Abs((Int_t)(z/stepa)); + if(z>=0) bina+=(Int_t)(nbina/2.); + if(bina>nbina) AliError("Wrong bin number!"); -void AliITSClusterFinderSDD::Print() -{ - // Print SDD cluster finder Parameters + Double_t devz = (Double_t)cal->GetMapACell(bina,bint)/10000.; + Double_t devx = (Double_t)cal->GetMapTCell(bina,bint)/10000.; + z+=devz; + y+=devx; - 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; }