X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=ITS%2FAliITSClusterFinderSDD.cxx;h=a73af85e4ef085716e94a0376a6083b99c424df2;hb=c6e3abc54c7d9d133d3410c9dfb640e3fb756a1b;hp=72dc8e16874a4d3edcdf2ffe791cd48c744add2c;hpb=22bba6479156012ec06fd53272ab3005487342de;p=u%2Fmrichter%2FAliRoot.git diff --git a/ITS/AliITSClusterFinderSDD.cxx b/ITS/AliITSClusterFinderSDD.cxx index 72dc8e16874..a73af85e4ef 100644 --- a/ITS/AliITSClusterFinderSDD.cxx +++ b/ITS/AliITSClusterFinderSDD.cxx @@ -13,8 +13,16 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ +#include #include "AliITSClusterFinderSDD.h" +#include "AliITSMapA1.h" +#include "AliITS.h" +#include "AliITSdigit.h" +#include "AliITSRawCluster.h" +#include "AliITSRecPoint.h" +#include "AliITSsegmentation.h" +#include "AliITSresponse.h" #include "AliRun.h" @@ -26,18 +34,21 @@ AliITSClusterFinderSDD::AliITSClusterFinderSDD (AliITSsegmentation *seg, AliITSresponse *response, TClonesArray *digits, TClonesArray *recp) { // constructor + fSegmentation=seg; fResponse=response; fDigits=digits; fClusters=recp; fNclusters= fClusters->GetEntriesFast(); - printf("SDD: fNclusters %d\n",fNclusters); SetCutAmplitude(); SetDAnode(); SetDTime(); - SetMap(); SetMinPeak(); - SetNCells(); + SetMinNCells(); + SetMaxNCells(); + SetTimeCorr(); + fMap=new AliITSMapA1(fSegmentation,fDigits,fCutAmplitude); + } //_____________________________________________________________________________ @@ -49,15 +60,25 @@ AliITSClusterFinderSDD::AliITSClusterFinderSDD() fDigits=0; fClusters=0; fNclusters=0; + fMap=0; SetCutAmplitude(); SetDAnode(); SetDTime(); - SetMap(); SetMinPeak(); - SetNCells(); + SetMinNCells(); + SetMaxNCells(); + SetTimeCorr(); } +//_____________________________________________________________________________ +AliITSClusterFinderSDD::~AliITSClusterFinderSDD() +{ + // destructor + + if(fMap) delete fMap; + +} //__________________________________________________________________________ AliITSClusterFinderSDD::AliITSClusterFinderSDD(const AliITSClusterFinderSDD &source){ // Copy Constructor @@ -68,8 +89,10 @@ AliITSClusterFinderSDD::AliITSClusterFinderSDD(const AliITSClusterFinderSDD &sou 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; } @@ -84,41 +107,14 @@ AliITSClusterFinderSDD& 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::SetMap() -{ - // set map - if(!fMap) fMap=new AliITSMapA2(fSegmentation); - -} -//_____________________________________________________________________________ -void AliITSClusterFinderSDD::FillMap() -{ - // fCoord1 = anode # - // fCoord2 = time sample - - if (!fDigits) return; - - Int_t ndigits = fDigits->GetEntriesFast(); - //printf("FillMap: ndigits %d\n",ndigits); - if (!ndigits) return; - - AliITSdigitSDD *dig; - Int_t ndig; - for(ndig=0; ndigUncheckedAt(ndig); - Double_t signal=dig->fSignal; - //printf("FillMap: ndig fCoord1 fCoord2 signal %d %d %d %f\n",ndig,dig->fCoord1,dig->fCoord2,signal); - fMap->SetHit(dig->fCoord1,dig->fCoord2,signal); - } - -} //_____________________________________________________________________________ void AliITSClusterFinderSDD::Find1DClusters() @@ -128,7 +124,6 @@ void AliITSClusterFinderSDD::Find1DClusters() AliITS *iTS=(AliITS*)gAlice->GetModule("ITS"); // retrieve the parameters - Int_t i; Int_t fNofMaps = fSegmentation->Npz(); Int_t fMaxNofSamples = fSegmentation->Npx(); Int_t fNofAnodes = fNofMaps/2; @@ -139,142 +134,146 @@ void AliITSClusterFinderSDD::Find1DClusters() Float_t anodePitch = fSegmentation->Dpz(dummy); // map the signal - FillMap(); + fMap->SetThreshold(fCutAmplitude); + fMap->FillMap(); - // Piergiorgio's code - do not subtract baseline since we start - // from digits and do not duplicate arrays, i.e. use fMap instead - // of Float_t fadc[2*fNofAnodes][fMaxNofSamples]; - + Float_t noise; + Float_t baseline; + fResponse->GetNoiseParam(noise,baseline); - Int_t nofFoundClusters = 0; + Float_t maxadc = fResponse->MaxAdc(); + Float_t topValue = fResponse->MagicValue(); + Float_t norm = maxadc/topValue; - Float_t **dfadc = new Float_t*[fNofMaps]; - for(i=0;iGetSignal(idx,l); - fadc1=(Float_t)fMap->GetSignal(idx,l-1); - if(l>0) dfadc[k][l-1] = fadc2-fadc1; - } // samples - } // anodes + 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= 60) lthrmina = 2; - // if(it >= 100) lthrmina = 3; - Int_t lthrmint = 2; - //if(it >= 60) lthrmint = 3; - //if(it >= 100) lthrmint = 4; - - Int_t lthra = 0; - Int_t lthrt = 0; - for(m=0;m<10;m++) { - Int_t id = it+m; - if(id>=fMaxNofSamples) break; - fadc=fMap->GetSignal(idx,id); - if(fadc > fadcmax) { - fadcmax = fadc; - if(fadc > fCutAmplitude) { lthra++; lthrt++; } - imax = id; - } - if(dfadc[k][id] > dfadcmax) { - dfadcmax = dfadc[k][id]; - imaxd = id; - } - } - it = imaxd; - // skip if no signal over threshold - if(fMap->GetSignal(idx,imax) < fCutAmplitude) {it++; continue;} - - if(k>0) { - if(fMap->GetSignal(idx-1,imax) > fCutAmplitude) lthra++; - } - if(kGetSignal(idx+1,imax) > fCutAmplitude) lthra++; - - if(imax>0) { - if(fMap->GetSignal(idx,imax-1) > fCutAmplitude) lthrt++; - } - if(imaxGetSignal(idx,imax+1) > fCutAmplitude) lthrt++; + Float_t fadcmax = 0.; + Float_t dfadcmax = 0.; + Int_t lthrmina = 1; + // if(it >= 60) lthrmina = 2; + // if(it >= 100) lthrmina = 3; + Int_t lthrmint = 2; + //if(it >= 60) lthrmint = 3; + //if(it >= 100) lthrmint = 4; + + Int_t lthra = 1; + Int_t lthrt = 0; - // cluster charge - Int_t tstart = it-1; + for(m=0;m<20;m++) { + Int_t id = it+m; + if(id>=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; - Bool_t ilcl = 0; - if(lthrt >= lthrmint && lthra >= lthrmina) ilcl = 1; + if(fMap->TestHit(idx,imax) == kEmpty) {it++; continue;} + + // cluster charge + Int_t tstart = it-1; - if(ilcl) { - nofFoundClusters++; - Int_t tstop = tstart; - Float_t dfadcmin = 10000.; - Int_t ij; - for(ij=0; ij<10; ij++) { - if(dfadc[k][it+ij] < dfadcmin) { - tstop = it+ij+1; - 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; - for(its=tstart; its<=tstop; its++) { - fadc=fMap->GetSignal(idx,its); - clusterCharge += fadc; - if(fadc > clusterPeakAmplitude) clusterPeakAmplitude = fadc; - clusterTime += fadc*its; - clusterMult++; - if(its == tstop) { + Bool_t ilcl = 0; + if(lthrt >= lthrmint && lthra >= lthrmina) ilcl = 1; + //printf("ilcl %d\n",ilcl); + + Float_t b,n; + fResponse->GetNoiseParam(n,b); + + if(ilcl) { + nofFoundClusters++; + Int_t tstop = tstart; + Float_t dfadcmin = 10000.; + Int_t ij; + for(ij=0; ij<20; ij++) { + if(dfadc[k][it+ij] < dfadcmin) { + tstop = it+ij+1; + 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; + clusterMult++; + if(its == tstop) { + // charge from ADC back to nA + //clusterCharge /= norm; + if(clusterCharge <= 0.) printf("clusterCharge %f norm %f\n",clusterCharge,norm); clusterTime /= (clusterCharge/fTimeStep); // ns clusterCharge *= (fTimeStep/160.); // keV - if(clusterTime > 58.2) clusterTime -= 58.2; // ns - /* - else { - cout << "Warning: cluster with negative time " << clusterTime << ", peak ampl.: " << clusterPeakAmplitude << ", mult: " << clusterMult << ", charge: " << clusterCharge << endl; - cout << "Anode: " << k << ", tstart: " << tstart << ", tstop: " << tstop << ", Charge: " << clusterCharge << endl; - } - */ - } - } - // cout << "Anode: " << k << ", tstart: " << tstart << ", tstop: " << tstop << ", Charge: " << clusterCharge << endl; - - Float_t clusteranodePath = (0.06 + clusterAnode - fNofAnodes/2)*anodePitch; - Float_t clusterDriftPath = clusterTime*fDriftSpeed; - clusterDriftPath = fSddLength-clusterDriftPath; - - if(clusterCharge < 0.) break; - - //printf("wing clusterMult clusterAnode clusterTime %d %f %f %f \n",j+1,clusterMult, clusterAnode, clusterTime); - - AliITSRawClusterSDD *clust = new AliITSRawClusterSDD(j+1,clusterAnode,clusterTime,clusterCharge,clusterPeakAmplitude,0.,0.,clusterDriftPath,clusteranodePath,clusterMult); - //fClusters->Add(point); - iTS->AddCluster(1,clust); - it = tstop; + if(clusterTime > fTimeCorr) clusterTime -= fTimeCorr; // ns + } + } + // cout << "Anode: " << k << ", tstart: " << tstart << ", tstop: " << tstop << ", Charge: " << clusterCharge << endl; + + Float_t clusteranodePath = (clusterAnode - fNofAnodes/2)*anodePitch; + Float_t clusterDriftPath = clusterTime*fDriftSpeed; + clusterDriftPath = fSddLength-clusterDriftPath; + + if(clusterCharge <= 0.) break; + + AliITSRawClusterSDD *clust = new AliITSRawClusterSDD(j+1,clusterAnode,clusterTime,clusterCharge,clusterPeakAmplitude,peakpos,0.,0.,clusterDriftPath,clusteranodePath,clusterMult); + iTS->AddCluster(1,clust); + it = tstop; } // ilcl it++; @@ -283,14 +282,10 @@ void AliITSClusterFinderSDD::Find1DClusters() } // anodes } // detectors (2) - Int_t nofClusters = fClusters->GetEntriesFast(); - nofClusters -= fNclusters; - - //printf("SDD- Find1Dclust: fNclusters nofClusters %d %d \n",fNclusters, nofClusters); - fMap->ClearMap(); + //fMap->ClearMap(); - for(i=0;iGetEntriesFast(); nofClusters -= fNclusters; - //printf("SDD- GroupClusters: fNclusters nofClusters %d %d \n",fNclusters, nofClusters); - AliITSRawClusterSDD *clusterI; AliITSRawClusterSDD *clusterJ; @@ -328,8 +321,8 @@ void AliITSClusterFinderSDD::GroupClusters() if(clusterI->T() < fTimeStep*10) fDAnode = 1.2; Bool_t pair = clusterI->Brother(clusterJ,fDAnode,fDTime); if(!pair) continue; - // clusterI->Print(); - // clusterJ->Print(); + // clusterI->PrintInfo(); + // clusterJ->PrintInfo(); clusterI->Add(clusterJ); label[j] = 1; fClusters->RemoveAt(j); @@ -351,8 +344,6 @@ void AliITSClusterFinderSDD::SelectClusters() Int_t nofClusters = fClusters->GetEntriesFast(); nofClusters -= fNclusters; - //printf("SDD- SelectClusters: fNclusters nofClusters %d %d \n",fNclusters, nofClusters); - Int_t i; for(i=0; iAt(i); @@ -364,6 +355,7 @@ void AliITSClusterFinderSDD::SelectClusters() Float_t amp = clusterI->PeakAmpl(); if(amp < fMinPeak) rmflg = 1; if(wy < fMinNCells) rmflg = 1; + //if(wy > fMaxNCells) rmflg = 1; if(rmflg) fClusters->RemoveAt(i); } // I clusters fClusters->Compress(); @@ -376,7 +368,6 @@ void AliITSClusterFinderSDD::SelectClusters() void AliITSClusterFinderSDD::GetRecPoints() { // get rec points - //static Int_t counter=0; AliITS *iTS=(AliITS*)gAlice->GetModule("ITS"); @@ -384,31 +375,49 @@ void AliITSClusterFinderSDD::GetRecPoints() Int_t nofClusters = fClusters->GetEntriesFast(); nofClusters -= fNclusters; - //printf("SDD- GetRecPoints: fNclusters nofClusters %d %d \n",fNclusters, nofClusters); - -// const Float_t kdEdXtoQ = 2.778e+2; // KeV -> number of e-hole pairs in Si 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 maxt=fSegmentation->Npx(); + Int_t ndigits=fDigits->GetEntriesFast(); for(i=0; iAt(i); + AliITSRawClusterSDD *clusterI = (AliITSRawClusterSDD*)fClusters->At(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"); + } + AliITSRecPoint rnew; rnew.SetX(clusterI->X()); rnew.SetZ(clusterI->Z()); rnew.SetQ(clusterI->Q()); // in KeV - should be ADC - //rnew.SetdEdX((clusterI->Q())/kdEdXtoQ); rnew.SetdEdX(kconvGeV*clusterI->Q()); rnew.SetSigmaX2(kRMSx*kRMSx); rnew.SetSigmaZ2(kRMSz*kRMSz); - rnew.SetProbability(1.); + 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); - //counter++; } // I clusters - //printf("counter %d\n",counter); + fMap->ClearMap(); } //_____________________________________________________________________________ @@ -421,7 +430,3 @@ void AliITSClusterFinderSDD::FindRawClusters() SelectClusters(); GetRecPoints(); } - - - -