X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=ITS%2FAliITSClusterFinder.cxx;h=aa3f4cf7169c5595dbbcfd3558b45ab9d7d883e3;hb=08b801a6de165ef29fb388dc4fc9c0516b7b9486;hp=885160100fb03cef5ad318f75d85dc1b691335a7;hpb=bf3f2830e911bfce646abd1592f6dd2801785c47;p=u%2Fmrichter%2FAliRoot.git diff --git a/ITS/AliITSClusterFinder.cxx b/ITS/AliITSClusterFinder.cxx index 885160100fb..aa3f4cf7169 100644 --- a/ITS/AliITSClusterFinder.cxx +++ b/ITS/AliITSClusterFinder.cxx @@ -12,217 +12,204 @@ * about the suitability of this software for any purpose. It is * * provided "as is" without express or implied warranty. * **************************************************************************/ -// -// Base Class used to find -// the reconstructed points for ITS -// See also AliITSClusterFinderSPD, AliITSClusterFinderSDD, -// AliITSClusterFinderSDD -// +//////////////////////////////////////////////////////////////////////////// +// // +// Base Class used to find // +// the reconstructed points for ITS // +// See also AliITSClusterFinderSPD, AliITSClusterFinderSDD, // +// AliITSClusterFinderSDD AliITSClusterFinderV2 // +//////////////////////////////////////////////////////////////////////////// + +#include "AliRun.h" #include "AliITSClusterFinder.h" +#include "AliITSRecPoint.h" #include "AliITSdigit.h" -#include "AliRun.h" -#include "AliITS.h" +#include "AliITSDetTypeRec.h" +#include "AliITSMap.h" +#include "AliITSgeomTGeo.h" +#include +#include +#include "AliMC.h" +#include "AliLog.h" + +using std::endl; ClassImp(AliITSClusterFinder) +extern AliRun *gAlice; + //---------------------------------------------------------------------- -AliITSClusterFinder::AliITSClusterFinder(){ +AliITSClusterFinder::AliITSClusterFinder(): +TObject(), +fModule(0), +fDigits(0), +fNdigits(0), +fDetTypeRec(0), +fClusters(0), +fMap(0), +fNPeaks(-1), +fNModules(AliITSgeomTGeo::GetNModules()), +fEvent(0), +fZmin(0), +fZmax(0), +fXmin(0), +fXmax(0), +fNClusters(0), +fRawID2ClusID(0) +{ // default cluster finder - - fSegmentation = 0; - fResponse = 0; - fMap = 0; - fDigits = 0; - fNdigits = 0; - fNRawClusters = 0; - fNperMax = 0; - fDeclusterFlag= 0; - fClusterSize = 0; - fNPeaks = 0; + // Input: + // none. + // Output: + // none. + // Return: + // A default constructed AliITSCulsterFinder + for(Int_t i=0; i<2200; i++){ + fNdet[i]=0; + fNlayer[i]=0; + } } //---------------------------------------------------------------------- -AliITSClusterFinder::AliITSClusterFinder(AliITSsegmentation *seg, - AliITSresponse *response, - TClonesArray *digits){ - // cluster finder - fSegmentation=seg; - fResponse=response; - fMap = 0; - - fDigits=digits; - fNdigits = fDigits->GetEntriesFast(); - - fNRawClusters=0; - - SetNperMax(); - SetClusterSize(); - SetDeclusterFlag(); - - fNPeaks=-1; +AliITSClusterFinder::AliITSClusterFinder(AliITSDetTypeRec* dettyp): +TObject(), +fModule(0), +fDigits(0), +fNdigits(0), +fDetTypeRec(dettyp), +fClusters(0), +fMap(0), +fNPeaks(-1), +fNModules(AliITSgeomTGeo::GetNModules()), +fEvent(0), +fZmin(0), +fZmax(0), +fXmin(0), +fXmax(0), +fNClusters(0), +fRawID2ClusID(0) +{ + // default cluster finder + // Standard constructor for cluster finder + // Input: + // AliITSsegmentation *seg The segmentation class to be used + // AliITSresponse *res The response class to be used + // Output: + // none. + // Return: + // A Standard constructed AliITSCulsterFinder + for(Int_t i=0; i<2200; i++){ + fNdet[i]=0; + fNlayer[i]=0; + } } //---------------------------------------------------------------------- -AliITSClusterFinder::~AliITSClusterFinder(){ - // destructor cluster finder +AliITSClusterFinder::AliITSClusterFinder(AliITSDetTypeRec* dettyp, + TClonesArray *digits): +TObject(), +fModule(0), +fDigits(digits), +fNdigits(0), +fDetTypeRec(dettyp), +fClusters(0), +fMap(0), +fNPeaks(-1), +fNModules(AliITSgeomTGeo::GetNModules()), +fEvent(0), +fZmin(0), +fZmax(0), +fXmin(0), +fXmax(0), +fNClusters(0), +fRawID2ClusID(0) +{ + // default cluster finder + // Standard + cluster finder constructor + // Input: + // AliITSsegmentation *seg The segmentation class to be used + // AliITSresponse *res The response class to be used + // TClonesArray *digits Array of digits to be used + // Output: + // none. + // Return: + // A Standard constructed AliITSCulsterFinder - // Zero local pointers. Other classes own these pointers. - fSegmentation = 0; - fResponse = 0; - fMap = 0; - fDigits = 0; - fNdigits = 0; - fNRawClusters = 0; - fNperMax = 0; - fDeclusterFlag= 0; - fClusterSize = 0; - fNPeaks = 0; -} -//__________________________________________________________________________ -AliITSClusterFinder::AliITSClusterFinder(const AliITSClusterFinder &source){ - // Copy Constructor - if(&source == this) return; - this->fDigits = source.fDigits; - this->fNdigits = source.fNdigits; - this->fResponse = source.fResponse; - this->fSegmentation = source.fSegmentation; - this->fNRawClusters = source.fNRawClusters; - this->fMap = source.fMap; - this->fNperMax = source.fNperMax; - this->fDeclusterFlag = source.fDeclusterFlag; - this->fClusterSize = source.fClusterSize; - this->fNPeaks = source.fNPeaks; - return; + fNdigits = fDigits->GetEntriesFast(); + for(Int_t i=0; i<2200; i++){ + fNdet[i]=0; + fNlayer[i]=0; + } } + //______________________________________________________________________ -AliITSClusterFinder& AliITSClusterFinder::operator=(const AliITSClusterFinder &source) { - // Assignment operator - - if(&source == this) return *this; - this->fDigits = source.fDigits; - this->fNdigits = source.fNdigits; - this->fResponse = source.fResponse; - this->fSegmentation = source.fSegmentation; - this->fNRawClusters = source.fNRawClusters; - this->fMap = source.fMap; - this->fNperMax = source.fNperMax; - this->fDeclusterFlag = source.fDeclusterFlag; - this->fClusterSize = source.fClusterSize; - this->fNPeaks = source.fNPeaks; - return *this; +AliITSClusterFinder::AliITSClusterFinder(const AliITSClusterFinder &source) : + TObject(source), + fModule(source.fModule), + fDigits(), + fNdigits(source.fNdigits), + fDetTypeRec(), + fClusters(), + fMap(), + fNPeaks(source.fNPeaks), + fNModules(source.fNModules), + fEvent(source.fEvent), + fZmin(source.fZmin), + fZmax(source.fZmax), + fXmin(source.fXmin), + fXmax(source.fXmax), + fNClusters(source.fNClusters), + fRawID2ClusID(source.fRawID2ClusID) +{ + // Copy constructor + // Copies are not allowed. The method is protected to avoid misuse. + AliError("Copy constructor not allowed\n"); } -//---------------------------------------------------------------------- -void AliITSClusterFinder::AddCluster(Int_t branch, AliITSRawCluster *c){ - // Add a raw cluster copy to the list - AliITS *iTS=(AliITS*)gAlice->GetModule("ITS"); - iTS->AddCluster(branch,c); - fNRawClusters++; -} -//---------------------------------------------------------------------- -void AliITSClusterFinder::AddCluster(Int_t branch, AliITSRawCluster *c, - AliITSRecPoint &rp){ - // Add a raw cluster copy to the list - - AliITS *iTS=(AliITS*)gAlice->GetModule("ITS"); - iTS->AddCluster(branch,c); - fNRawClusters++; - iTS->AddRecPoint(rp); -} + //______________________________________________________________________ -void AliITSClusterFinder::FindRawClusters(Int_t module){ - // Default Cluster finder. +//AliITSClusterFinder& AliITSClusterFinder::operator=(const AliITSClusterFinder& /* source */){ + // Assignment operator + // Assignment is not allowed. The method is protected to avoid misuse. +// Fatal("= operator","Assignment operator not allowed\n"); +// return *this; +//} + +//---------------------------------------------------------------------- +AliITSClusterFinder::~AliITSClusterFinder(){ + // destructor cluster finder // Input: - // Int_t module Module number for which culster are to be found. + // none. // Output: // none. // Return: // none. - const Int_t kelms = 10; - Int_t ndigits = fDigits->GetEntriesFast(); - TObjArray *digs = new TObjArray(ndigits); - TObjArray *clusts = new TObjArray(ndigits); // max # cluster - TObjArray *clust0=0; // A spacific cluster of digits - TObjArray *clust1=0; // A spacific cluster of digits - AliITSdigit *dig=0; // locat pointer to a digit - Int_t i=0,nc=0,j[4],k,k2=0; - - // Copy all digits for this module into a local TObjArray. - for(i=0;iAddAt(new AliITSdigit(*((AliITSdigit*)(fDigits->At(i)))),i); - digs->Sort(); - // First digit is a cluster. - i = 0; - nc = 0; - clusts->AddAt(new TObjArray(kelms),nc); - clust0 = (TObjArray*)(clusts->At(nc)); - clust0->AddAtFree(digs->At(i)); // move owner ship from digs to clusts - nc++; - for(i=1;iAt(j[0])); - // Add to existing cluster. Find which cluster this digis - for(k=0;kAt(k))); - if(clust0->IndexOf(dig)>=0) break; - } // end for k - if(k>=nc){ - Fatal("FindRawClusters","Digit not found as expected"); - } // end if - if(j[1]>=0){ - dig = (AliITSdigit*)(digs->At(j[1])); - // Add to existing cluster. Find which cluster this digis - for(k2=0;k2At(k2))); - if(clust1->IndexOf(dig)>=0) break; - } // end for k2 - if(k2>=nc){ - Fatal("FindRawClusters","Digit not found as expected"); - } // end if - } // end if j[1]>=0 - // Found cluster with neighboring digits add this one to it. - if(clust0==clust1){ // same cluster - clust0->AddAtFree(digs->At(i)); - clust0 = 0; // finished with cluster. zero for safty - clust1 = 0; // finished wit hcluster. zero for safty - }else{ // two different clusters which need to be merged. - clust0->AddAtFree(digs->At(i)); // Add digit to this cluster. - for(k=0;kGetEntriesFast();k++){ - // move clust1 into clust0 - clust0->AddAtFree(clust1->At(k));//move digit to this cluster - clust1->AddAt(0,k); // zero this one - } // end for k - delete clust1; - clusts->AddAt(0,k2); // zero array of clusters element clust1 - clust0 = 0; // finished with cluster. zero for safty - clust1 = 0; // finished wit hcluster. zero for safty - } // end if clust0==clust1 - }else{// New cluster - clusts->AddAt(new TObjArray(kelms),nc); - clust0 = ((TObjArray*)(clusts->At(nc))); - clust0->AddAtFree(digs->At(i));// move owner ship from digs to clusts - clust0 = 0; // finished with cluster. zero for safty - nc++; - } // End if IsNeighbor - } // end for i - // There are now nc clusters in clusts. Each element of clust is an - // array of digits which are clustered together. - - // For each cluster call detector specific CreateRecPoints - for(i=0;iAt(i)),module); - - // clean up at the end. - for(i=0;iAt(i)); - // Digits deleted below, so zero this TObjArray - for(k=0;kGetEntriesFast();k++) clust0->AddAt(0,k); - delete clust0; // Delete this TObjArray - clusts->AddAt(0,i); // Contents deleted above, so zero it. - } // end for i - delete clusts; // Delete this TObjArray/ - // Delete the digits then the TObjArray which containted them. - for(i=0;iAt(i))); - delete digs; + + if(fMap) {delete fMap;} + // Zero local pointers. Other classes own these pointers. + fMap = 0; + fDigits = 0; + fNdigits = 0; + fNPeaks = 0; + fDetTypeRec = 0; + } +//__________________________________________________________________________ +void AliITSClusterFinder::InitGeometry(){ + // + // Initialisation of ITS geometry + // + Int_t mmax=AliITSgeomTGeo::GetNModules(); + for (Int_t m=0; mAt(i)))->GetCoord1(); iz = ((AliITSdigit*)(digs->At(i)))->GetCoord2(); for(j=0;jAt(j)))->GetCoord1(); - jz = ((AliITSdigit*)(digs->At(j)))->GetCoord2(); - if(jx+1==ix && jz ==iz){n[0] = j;nei[0] = kTRUE;} - if(jx ==ix && jz+1==iz){n[1] = j;nei[1] = kTRUE;} - if(jx+1==ix && jz+1==iz){n[2] = j;nei[2] = kTRUE;} - if(jx+1==ix && jz-1==iz){n[3] = j;nei[3] = kTRUE;} + jx = ((AliITSdigit*)(digs->At(j)))->GetCoord1(); + jz = ((AliITSdigit*)(digs->At(j)))->GetCoord2(); + if(jx+1==ix && jz ==iz){n[0] = j;nei[0] = kTRUE;} + if(jx ==ix && jz+1==iz){n[1] = j;nei[1] = kTRUE;} + if(jx+1==ix && jz+1==iz){n[2] = j;nei[2] = kTRUE;} + if(jx+1==ix && jz-1==iz){n[3] = j;nei[3] = kTRUE;} } // end for k if(nei[0]||nei[1]) return kTRUE; if(kdiagonal&&(nei[2]||nei[3])) return kTRUE; // no Neighbors found. return kFALSE; } + +//______________________________________________________________________ +void AliITSClusterFinder::Print(ostream *os) const{ + //Standard output format for this class + // Inputs: + // ostream *os Output stream + // Output: + // ostream *os Output stream + // Return: + // none. + + *os << fModule<<","; + *os << fNdigits<<","; + *os << fNPeaks<> fModule; + *is >> fNdigits; + *is >> fNPeaks; +} +//______________________________________________________________________ +ostream &operator<<(ostream &os,AliITSClusterFinder &source){ + // Standard output streaming function. + // Inputs: + // ostream *os Output stream + // AliITSClusterFinder &source Class to be printed + // Output: + // ostream *os Output stream + // Return: + // none. + + source.Print(&os); + return os; +} +//______________________________________________________________________ +istream &operator>>(istream &is,AliITSClusterFinder &source){ + // Standard output streaming function. + // Inputs: + // istream *is Input stream + // AliITSClusterFinder &source Class to be read in. + // Output: + // istream *is Input stream + // Return: + // none. + + source.Read(&is); + return is; +} + +//______________________________________________________________________ +void AliITSClusterFinder::CheckLabels2(Int_t lab[10]) +{ + //------------------------------------------------------------ + // Tries to find mother's labels + //------------------------------------------------------------ + AliRunLoader *rl = AliRunLoader::Instance(); + if(!rl) return; + TTree *trK=(TTree*)rl->TreeK(); + if (!trK) return; + // + int labS[10]; + Int_t nlabels = 0; + Int_t ntracks = gAlice->GetMCApp()->GetNtrack(); + for (Int_t i=0;i<10;i++) if (lab[i]>=0) labS[nlabels++] = lab[i]; + if (nlabels==0) return; + // + float mom[10]; + for (Int_t i=0;i=ntracks) continue; + TParticle *part=(TParticle*)gAlice->GetMCApp()->Particle(label); + mom[i] = part->P(); + if (part->P() < 0.02) { // reduce soft particles from the same cluster + Int_t m=part->GetFirstMother(); + if (m<0) continue; // primary + // + if (part->GetStatusCode()>0) continue; + // + // if the parent is within the same cluster, reassign the label to it + for (int j=0;j3) { // only 3 labels are stored in cluster, sort in decreasing momentum + int ind[10],labSS[10]; + TMath::Sort(nlabels,mom,ind); + for (int i=nlabels;i--;) labSS[i] = labS[i]; + for (int i=0;iTreeK(); + if(trK){ + if(label<0) return; // In case of no label just exit + + Int_t ntracks = gAlice->GetMCApp()->GetNtrack(); + if (label>ntracks) return; + for (Int_t i=0;i<10;i++){ + // if (label<0) break; + if (lab[i]==label) break; + if (lab[i]<0) { + lab[i]= label; + break; + } + } + } +} + + +//______________________________________________________________________ +void AliITSClusterFinder:: +FindCluster(Int_t k,Int_t maxz,AliBin *bins,Int_t &n,Int_t *idx) { + //------------------------------------------------------------ + // returns an array of indices of digits belonging to the cluster + // (needed when the segmentation is not regular) + //------------------------------------------------------------ + if (n<200) idx[n++]=bins[k].GetIndex(); + bins[k].Use(); + + if (bins[k-maxz].IsNotUsed()) FindCluster(k-maxz,maxz,bins,n,idx); + if (bins[k-1 ].IsNotUsed()) FindCluster(k-1 ,maxz,bins,n,idx); + if (bins[k+maxz].IsNotUsed()) FindCluster(k+maxz,maxz,bins,n,idx); + if (bins[k+1 ].IsNotUsed()) FindCluster(k+1 ,maxz,bins,n,idx); + /* + if (bins[k-maxz-1].IsNotUsed()) FindCluster(k-maxz-1,maxz,bins,n,idx); + if (bins[k-maxz+1].IsNotUsed()) FindCluster(k-maxz+1,maxz,bins,n,idx); + if (bins[k+maxz-1].IsNotUsed()) FindCluster(k+maxz-1,maxz,bins,n,idx); + if (bins[k+maxz+1].IsNotUsed()) FindCluster(k+maxz+1,maxz,bins,n,idx); + */ +} + +//______________________________________________________________________ +Bool_t AliITSClusterFinder::IsMaximum(Int_t k,Int_t max,const AliBin *bins) { + //------------------------------------------------------------ + //is this a local maximum ? + //------------------------------------------------------------ + UShort_t q=bins[k].GetQ(); + if (q==1023) return kFALSE; + if (bins[k-max].GetQ() > q) return kFALSE; + if (bins[k-1 ].GetQ() > q) return kFALSE; + if (bins[k+max].GetQ() > q) return kFALSE; + if (bins[k+1 ].GetQ() > q) return kFALSE; + if (bins[k-max-1].GetQ() > q) return kFALSE; + if (bins[k+max-1].GetQ() > q) return kFALSE; + if (bins[k+max+1].GetQ() > q) return kFALSE; + if (bins[k-max+1].GetQ() > q) return kFALSE; + return kTRUE; +} + +//______________________________________________________________________ +void AliITSClusterFinder:: +FindPeaks(Int_t k,Int_t max,AliBin *b,Int_t *idx,UInt_t *msk,Int_t& n) { + //------------------------------------------------------------ + //find local maxima + //------------------------------------------------------------ + if (n<31) + if (IsMaximum(k,max,b)) { + idx[n]=k; msk[n]=(2<fXmax) fXmax=i; + if(jfZmax) fZmax=j; + } + c.SetQ(c.GetQ()+q); + c.SetY(c.GetY()+i*q); + c.SetZ(c.GetZ()+j*q); + c.SetSigmaY2(c.GetSigmaY2()+i*i*q); + c.SetSigmaZ2(c.GetSigmaZ2()+j*j*q); + + bins[k].SetMask(0xFFFFFFFE); + if (fRawID2ClusID) { // RS: Register cluster id in raw words list + int rwid = bins[k].GetRawID(); + if (fRawID2ClusID->GetSize()<=rwid) fRawID2ClusID->Set( (rwid+10)<<1 ); + (*fRawID2ClusID)[rwid] = fNClusters+1; // RS: store clID+1 as a reference to the cluster + } + if (bins[k-max].GetMask() == m) MakeCluster(k-max,max,bins,m,c); + if (bins[k-1 ].GetMask() == m) MakeCluster(k-1 ,max,bins,m,c); + if (bins[k+max].GetMask() == m) MakeCluster(k+max,max,bins,m,c); + if (bins[k+1 ].GetMask() == m) MakeCluster(k+1 ,max,bins,m,c); +}