X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=ITS%2FAliITSclustererV2.cxx;h=2cafb669215d0cf0cb1936e08686d691713e38ed;hb=ecc64c3f29d0ba10f97cb96ee82bd71df6ecfe60;hp=43d477589743b1d6873cf2eabdf11dc2c9d3145f;hpb=b6087704af91ee48c31cce3d9aa82824c67fba3a;p=u%2Fmrichter%2FAliRoot.git diff --git a/ITS/AliITSclustererV2.cxx b/ITS/AliITSclustererV2.cxx index 43d47758974..2cafb669215 100644 --- a/ITS/AliITSclustererV2.cxx +++ b/ITS/AliITSclustererV2.cxx @@ -10,6 +10,7 @@ #include "AliITSclustererV2.h" #include "AliITSclusterV2.h" +#include "AliITSDetTypeRec.h" #include "AliRawReader.h" #include "AliITSRawStreamSPD.h" #include "AliITSRawStreamSDD.h" @@ -19,7 +20,9 @@ #include #include #include "AliITSgeom.h" -#include "AliITSdigit.h" +#include "AliITSdigitSPD.h" +#include "AliITSdigitSDD.h" +#include "AliITSdigitSSD.h" #include "AliMC.h" ClassImp(AliITSclustererV2) @@ -49,6 +52,7 @@ AliITSclustererV2::AliITSclustererV2(const AliITSgeom *geom) { fYshift[m] = x*ca + y*sa; fZshift[m] = (Double_t)z; fNdet[m] = (lad-1)*g->GetNdetectors(lay) + (det-1); + fNlayer[m] = lay-1; } fNModules = g->GetIndexMax(); @@ -88,8 +92,10 @@ AliITSclustererV2::AliITSclustererV2(const AliITSgeom *geom) { fHlSSD=2.00; fTanP=0.0275; fTanN=0.0075; + } + Int_t AliITSclustererV2::Digits2Clusters(TTree *dTree, TTree *cTree) { //------------------------------------------------------------ // This function creates ITS clusters @@ -114,17 +120,22 @@ Int_t AliITSclustererV2::Digits2Clusters(TTree *dTree, TTree *cTree) { else branch->SetAddress(&clusters); Int_t mmax=(Int_t)dTree->GetEntries(); + if (mmax!=fNModules) { + Error("Digits2Clusters","Number of entries != number of modules !"); + return 1; + } for (fI=0; fIGetEvent(fI); if (digitsSPD->GetEntriesFast()!=0) - FindClustersSPD(digitsSPD,clusters); - else if(digitsSDD->GetEntriesFast()!=0) - FindClustersSDD(digitsSDD,clusters); - else if(digitsSSD->GetEntriesFast()!=0) - FindClustersSSD(digitsSSD,clusters); - + FindClustersSPD(digitsSPD,clusters); + else + if(digitsSDD->GetEntriesFast()!=0) + FindClustersSDD(digitsSDD,clusters); + else if(digitsSSD->GetEntriesFast()!=0) + FindClustersSSD(digitsSSD,clusters); + ncl+=clusters->GetEntriesFast(); cTree->Fill(); @@ -159,6 +170,7 @@ void AliITSclustererV2::Digits2Clusters(AliRawReader* rawReader) { Error("Digits2Clusters", "no run loader found"); return; } + runLoader->LoadKinematics(); AliLoader* itsLoader = runLoader->GetLoader("ITSLoader"); if (!itsLoader) { Error("Digits2Clusters", "no loader for ITS found"); @@ -213,12 +225,13 @@ void AliITSclustererV2::Digits2Clusters(AliRawReader* rawReader) { //**** Fast clusters ******************************* #include "TParticle.h" -#include "AliITS.h" +//#include "AliITS.h" #include "AliITSmodule.h" #include "AliITSRecPoint.h" #include "AliITSsimulationFastPoints.h" #include "AliITSRecPoint.h" + static void CheckLabels(Int_t lab[3]) { //------------------------------------------------------------ // Tries to find mother's labels @@ -246,6 +259,118 @@ static void CheckLabels(Int_t lab[3]) { } } +/* +static void CheckLabels(Int_t lab[3]) { + //------------------------------------------------------------ + // Tries to find mother's labels + //------------------------------------------------------------ + + if(lab[0]<0 && lab[1]<0 && lab[2]<0) return; // In case of no labels just exit + + Int_t ntracks = gAlice->GetMCApp()->GetNtrack(); + for (Int_t i=0;i<3;i++){ + Int_t label = lab[i]; + if (label>=0 && labelGetMCApp()->Particle(label); + if (part->P() < 0.005) { + Int_t m=part->GetFirstMother(); + if (m<0) { + continue; + } + if (part->GetStatusCode()>0) { + continue; + } + lab[i]=m; + } + } + } + +} +*/ +static void CheckLabels2(Int_t lab[10]) { + //------------------------------------------------------------ + // Tries to find mother's labels + //------------------------------------------------------------ + Int_t nlabels =0; + for (Int_t i=0;i<10;i++) if (lab[i]>=0) nlabels++; + if(nlabels == 0) return; // In case of no labels just exit + + Int_t ntracks = gAlice->GetMCApp()->GetNtrack(); + + for (Int_t i=0;i<10;i++){ + Int_t label = lab[i]; + if (label>=0 && labelGetMCApp()->Particle(label); + if (part->P() < 0.02) { + Int_t m=part->GetFirstMother(); + if (m<0) { + continue; + } + if (part->GetStatusCode()>0) { + continue; + } + lab[i]=m; + } + else + if (part->P() < 0.12 && nlabels>3) { + lab[i]=-2; + nlabels--; + } + } + else{ + if ( (label>ntracks||label <0) && nlabels>3) { + lab[i]=-2; + nlabels--; + } + } + } + if (nlabels>3){ + for (Int_t i=0;i<10;i++){ + if (nlabels>3){ + Int_t label = lab[i]; + if (label>=0 && labelGetMCApp()->Particle(label); + if (part->P() < 0.1) { + lab[i]=-2; + nlabels--; + } + } + } + } + } + + //compress labels -- if multi-times the same + Int_t lab2[10]; + for (Int_t i=0;i<10;i++) lab2[i]=-2; + for (Int_t i=0;i<10 ;i++){ + if (lab[i]<0) continue; + for (Int_t j=0;j<10 &&lab2[j]!=lab[i];j++){ + if (lab2[j]<0) { + lab2[j]= lab[i]; + break; + } + } + } + for (Int_t j=0;j<10;j++) lab[j]=lab2[j]; + +} + +static void AddLabel(Int_t lab[10], Int_t label) { + + 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 AliITSclustererV2::RecPoints2Clusters (const TClonesArray *points, Int_t idx, TClonesArray *clusters) { //------------------------------------------------------------ @@ -257,68 +382,20 @@ void AliITSclustererV2::RecPoints2Clusters for (Int_t i=0; iUncheckedAt(i); Float_t lp[5]; - lp[0]=-(-p->GetX()+fYshift[idx]); if (idx<=fLastSPD1) lp[0]*=-1; //SPD1 + lp[0]=-(-p->GetDetLocalX()+fYshift[idx]); if (idx<=fLastSPD1) lp[0]*=-1; //SPD1 lp[1]= -p->GetZ()+fZshift[idx]; - lp[2]=p->GetSigmaX2(); + lp[2]=p->GetSigmaDetLocX2(); lp[3]=p->GetSigmaZ2(); lp[4]=p->GetQ()*36./23333.; //electrons -> ADC Int_t lab[4]; lab[0]=p->GetLabel(0); lab[1]=p->GetLabel(1); lab[2]=p->GetLabel(2); lab[3]=fNdet[idx]; CheckLabels(lab); - new (cl[i]) AliITSclusterV2(lab,lp); + Int_t dummy[3]={0,0,0}; + new (cl[i]) AliITSclusterV2(lab,lp, dummy); } } -Int_t AliITSclustererV2::Hits2Clusters(TTree *hTree, TTree *cTree) { - //------------------------------------------------------------ - // This function creates ITS clusters - //------------------------------------------------------------ - if (!gAlice) { - Error("Hits2Clusters","gAlice==0 !"); - return 1; - } - - AliITS *its = (AliITS*)gAlice->GetModule("ITS"); - if (!its) { - Error("Hits2Clusters","Can't find the ITS !"); - return 2; - } - AliITSgeom *geom=its->GetITSgeom(); - Int_t mmax=geom->GetIndexMax(); - - its->InitModules(-1,mmax); - its->FillModules(hTree,0); - - TClonesArray *clusters=new TClonesArray("AliITSclusterV2",1000); - TBranch *branch=cTree->GetBranch("Clusters"); - if (!branch) cTree->Branch("Clusters",&clusters); - else branch->SetAddress(&clusters); - - static TClonesArray *points=its->RecPoints(); - AliITSsimulationFastPoints sim; - Int_t ncl=0; - for (Int_t m=0; mGetModule(m); - sim.CreateFastRecPoints(mod,m,gRandom); - - RecPoints2Clusters(points, m, clusters); - its->ResetRecPoints(); - - ncl+=clusters->GetEntriesFast(); - cTree->Fill(); - clusters->Clear(); - } - - Info("Hits2Clusters","Number of found fast clusters : %d",ncl); - - //cTree->Write(); - - delete clusters; - - return 0; -} - //*********************************** #ifndef V1 @@ -374,16 +451,15 @@ FindClustersSPD(const TClonesArray *digits, TClonesArray *clusters) { Info("FindClustersSPD","Too big cluster !"); continue; } - Int_t lab[4]; - lab[0]=-2; - lab[1]=-2; - lab[2]=-2; - lab[3]=fNdet[fI]; + Int_t milab[10]; + for (Int_t ilab=0;ilab<10;ilab++){ + milab[ilab]=-2; + } d=(AliITSdigitSPD*)digits->UncheckedAt(idx[0]); Int_t ymin=d->GetCoord2(),ymax=ymin; Int_t zmin=d->GetCoord1(),zmax=zmin; - Float_t y=0.,z=0.,q=0.; + for (Int_t l=0; lUncheckedAt(idx[l]); @@ -391,37 +467,65 @@ FindClustersSPD(const TClonesArray *digits, TClonesArray *clusters) { if (ymax < d->GetCoord2()) ymax=d->GetCoord2(); if (zmin > d->GetCoord1()) zmin=d->GetCoord1(); if (zmax < d->GetCoord1()) zmax=d->GetCoord1(); - - Int_t lab0=(d->GetTracks())[0]; - if (lab0>=0) { - if (lab[0]<0) { - lab[0]=lab0; - } else if (lab[1]<0) { - if (lab0!=lab[0]) lab[1]=lab0; - } else if (lab[2]<0) { - if (lab0!=lab[0]) - if (lab0!=lab[1]) lab[2]=lab0; - } - } - Float_t qq=d->GetSignal(); - y+=qq*fYSPD[d->GetCoord2()]; z+=qq*fZSPD[d->GetCoord1()]; q+=qq; + // MI addition - find all labels in cluster + for (Int_t dlab=0;dlab<10;dlab++){ + Int_t digitlab = (d->GetTracks())[dlab]; + if (digitlab<0) continue; + AddLabel(milab,digitlab); + } + if (milab[9]>0) CheckLabels2(milab); } - y/=q; z/=q; - y-=fHwSPD; z-=fHlSPD; - - Float_t lp[5]; - lp[0]=-(-y+fYshift[fI]); if (fI<=fLastSPD1) lp[0]=-lp[0]; - lp[1]= -z+fZshift[fI]; - lp[2]= fYpitchSPD*fYpitchSPD/12.; - lp[3]= fZ1pitchSPD*fZ1pitchSPD/12.; - //lp[4]= q; - lp[4]= (zmax-zmin+1)*100 + (ymax-ymin+1); - - //CheckLabels(lab); - d=(AliITSdigitSPD*)digits->UncheckedAt(idx[0]); - new (cl[n]) AliITSclusterV2(lab,lp); n++; + CheckLabels2(milab); + // + //Int_t idy = (fNlayer[fI]==0)? 2:3; + //for (Int_t iz=zmin; iz<=zmax;iz+=2) + //Int_t idy = (ymax-ymin)/4.; // max 2 clusters + Int_t idy = 0; // max 2 clusters + if (fNlayer[fI]==0 &&idy<3) idy=3; + if (fNlayer[fI]==1 &&idy<4) idy=4; + Int_t idz =3; + for (Int_t iz=zmin; iz<=zmax;iz+=idz) + for (Int_t iy=ymin; iy<=ymax;iy+=idy){ + // + Int_t ndigits =0; + Float_t y=0.,z=0.,q=0.; + for (Int_t l=0; lUncheckedAt(idx[l]); + if (zmax-zmin>=idz || ymax-ymin>=idy){ + if (TMath::Abs( d->GetCoord2()-iy)>0.75*idy) continue; + if (TMath::Abs( d->GetCoord1()-iz)>0.75*idz) continue; + } + ndigits++; + Float_t qq=d->GetSignal(); + y+=qq*fYSPD[d->GetCoord2()]; z+=qq*fZSPD[d->GetCoord1()]; q+=qq; + + } + if (ndigits==0) continue; + y/=q; z/=q; + y-=fHwSPD; z-=fHlSPD; + + Float_t lp[5]; + lp[0]=-(-y+fYshift[fI]); if (fI<=fLastSPD1) lp[0]=-lp[0]; + lp[1]= -z+fZshift[fI]; + // Float_t factor=TMath::Max(double(ni-3.),1.5); + Float_t factory=TMath::Max(ymax-ymin,1); + Float_t factorz=TMath::Max(zmax-zmin,1); + factory*= factory; + factorz*= factorz; + //lp[2]= (fYpitchSPD*fYpitchSPD/12.)*factory; + //lp[3]= (fZ1pitchSPD*fZ1pitchSPD/12.)*factorz; + lp[2]= (fYpitchSPD*fYpitchSPD/12.); + lp[3]= (fZ1pitchSPD*fZ1pitchSPD/12.); + //lp[4]= q; + lp[4]= (zmax-zmin+1)*100 + (ymax-ymin+1); + + milab[3]=fNdet[fI]; + d=(AliITSdigitSPD*)digits->UncheckedAt(idx[0]); + Int_t info[3] = {ymax-ymin+1,zmax-zmin+1,fNlayer[fI]}; + new (cl[n]) AliITSclusterV2(milab,lp,info); n++; + } } - + delete [] bins; } @@ -436,7 +540,9 @@ void AliITSclustererV2::FindClustersSPD(AliITSRawStream* input, Int_t kNzBins = fNzSPD + 2; Int_t kNyBins = fNySPD + 2; Int_t kMaxBin = kNzBins * kNyBins; - AliBin* bins = NULL; + AliBin *binsSPD = new AliBin[kMaxBin]; + AliBin *binsSPDInit = new AliBin[kMaxBin]; + AliBin *bins = NULL; // read raw data input stream while (kTRUE) { @@ -470,9 +576,6 @@ void AliITSclustererV2::FindClustersSPD(AliITSRawStream* input, Int_t ymax = ymin; Int_t zmin = (idxBins[0] % kNzBins) - 1; Int_t zmax = zmin; - Float_t y = 0.; - Float_t z = 0.; - Float_t q = 0.; for (Int_t idx = 0; idx < nBins; idx++) { Int_t iy = (idxBins[idx] / kNzBins) - 1; Int_t iz = (idxBins[idx] % kNzBins) - 1; @@ -480,38 +583,59 @@ void AliITSclustererV2::FindClustersSPD(AliITSRawStream* input, if (ymax < iy) ymax = iy; if (zmin > iz) zmin = iz; if (zmax < iz) zmax = iz; - - Float_t qBin = bins[idxBins[idx]].GetQ(); - y += qBin * fYSPD[iy]; - z += qBin * fZSPD[iz]; - q += qBin; } - y /= q; - z /= q; - y -= fHwSPD; - z -= fHlSPD; - - Float_t hit[5]; // y, z, sigma(y)^2, sigma(z)^2, charge - hit[0] = -y-fYshift[iModule]; - if (iModule <= fLastSPD1) hit[0] = -hit[0]; - hit[1] = z+fZshift[iModule]; - hit[2] = fYpitchSPD*fYpitchSPD/12.; - hit[3] = fZ1pitchSPD*fZ1pitchSPD/12.; -// hit[4] = q; - hit[4] = (zmax-zmin+1)*100 + (ymax-ymin+1); - - //CheckLabels(label); - new (clusters[iModule]->AddrAt(nClusters)) - AliITSclusterV2(label, hit); - nClusters++; + + Int_t idy = 0; // max 2 clusters + if ((iModule <= fLastSPD1) &&idy<3) idy=3; + if ((iModule > fLastSPD1) &&idy<4) idy=4; + Int_t idz =3; + for (Int_t iiz=zmin; iiz<=zmax;iiz+=idz) + for (Int_t iiy=ymin; iiy<=ymax;iiy+=idy){ + // + Int_t ndigits =0; + Float_t y=0.,z=0.,q=0.; + for (Int_t idx = 0; idx < nBins; idx++) { + Int_t iy = (idxBins[idx] / kNzBins) - 1; + Int_t iz = (idxBins[idx] % kNzBins) - 1; + if (zmax-zmin>=idz || ymax-ymin>=idy){ + if (TMath::Abs(iy-iiy)>0.75*idy) continue; + if (TMath::Abs(iz-iiz)>0.75*idz) continue; + } + ndigits++; + Float_t qBin = bins[idxBins[idx]].GetQ(); + y += qBin * fYSPD[iy]; + z += qBin * fZSPD[iz]; + q += qBin; + } + if (ndigits==0) continue; + y /= q; + z /= q; + y -= fHwSPD; + z -= fHlSPD; + + Float_t hit[5]; // y, z, sigma(y)^2, sigma(z)^2, charge + hit[0] = -(-y+fYshift[iModule]); + if (iModule <= fLastSPD1) hit[0] = -hit[0]; + hit[1] = -z+fZshift[iModule]; + hit[2] = fYpitchSPD*fYpitchSPD/12.; + hit[3] = fZ1pitchSPD*fZ1pitchSPD/12.; + // hit[4] = q; + hit[4] = (zmax-zmin+1)*100 + (ymax-ymin+1); + // CheckLabels(label); + Int_t info[3]={ymax-ymin+1,zmax-zmin+1,fNlayer[iModule]}; + new (clusters[iModule]->AddrAt(nClusters)) + AliITSclusterV2(label, hit,info); + nClusters++; + } } nClustersSPD += nClusters; - delete bins; + bins = NULL; } if (!next) break; - bins = new AliBin[kMaxBin]; + bins = binsSPD; + memcpy(binsSPD,binsSPDInit,sizeof(AliBin)*kMaxBin); } // fill the current digit into the bins array @@ -521,6 +645,9 @@ void AliITSclustererV2::FindClustersSPD(AliITSRawStream* input, bins[index].SetQ(1); } + delete [] binsSPDInit; + delete [] binsSPD; + Info("FindClustersSPD", "found clusters in ITS SPD: %d", nClustersSPD); } @@ -614,6 +741,7 @@ FindClustersSDD(AliBin* bins[2], Int_t nMaxBin, Int_t nzBins, FindPeaks(i, nzBins, bins[s], idx, msk, npeaks); if (npeaks>30) continue; + if (npeaks==0) continue; Int_t k,l; for (k=0; kGetPadPitchWidth(sec); - c.SetSigmaY2((s2 + 1./12.)*w*w); - if (s2 != 0.) { - c.SetSigmaY2(c.GetSigmaY2()*0.108); - if (secGetNInnerSector()) c.SetSigmaY2(c.GetSigmaY2()*2.07); - } - - s2 = c.GetSigmaZ2()/c.GetQ() - c.GetZ()*c.GetZ(); - w=par->GetZWidth(); - c.SetSigmaZ2((s2 + 1./12.)*w*w); - if (s2 != 0.) { - c.SetSigmaZ2(c.GetSigmaZ2()*0.169); - if (secGetNInnerSector()) c.SetSigmaZ2(c.GetSigmaZ2()*1.77); - } + //mi change + Int_t milab[10]; + for (Int_t ilab=0;ilab<10;ilab++){ + milab[ilab]=-2; + } + Int_t maxi=0,mini=0,maxj=0,minj=0; + //AliBin *bmax=&bins[s][idx[k]]; + //Float_t max = TMath::Max(TMath::Abs(bmax->GetQ())/5.,3.); + Float_t max=3; + for (Int_t di=-2; di<=2;di++) + for (Int_t dj=-3;dj<=3;dj++){ + Int_t index = idx[k]+di+dj*nzBins; + if (index<0) continue; + if (index>=nMaxBin) continue; + AliBin *b=&bins[s][index]; + if (TMath::Abs(b->GetQ())>max){ + if (di>maxi) maxi=di; + if (dimaxj) maxj=dj; + if (djUncheckedAt(b->GetIndex()); + for (Int_t itrack=0;itrack<10;itrack++){ + Int_t track = (d->GetTracks())[itrack]; + if (track>=0) { + AddLabel(milab, track); + } + } + } + } + } + } + + /* + Float_t s2 = c.GetSigmaY2()/c.GetQ() - c.GetY()*c.GetY(); + Float_t w=par->GetPadPitchWidth(sec); + c.SetSigmaY2(s2); + if (s2 != 0.) { + c.SetSigmaY2(c.GetSigmaY2()*0.108); + if (secGetNInnerSector()) c.SetSigmaY2(c.GetSigmaY2()*2.07); + } + s2 = c.GetSigmaZ2()/c.GetQ() - c.GetZ()*c.GetZ(); + w=par->GetZWidth(); + c.SetSigmaZ2(s2); + + if (s2 != 0.) { + c.SetSigmaZ2(c.GetSigmaZ2()*0.169); + if (secGetNInnerSector()) c.SetSigmaZ2(c.GetSigmaZ2()*1.77); + } */ c.SetSigmaY2(0.0030*0.0030); @@ -671,7 +831,12 @@ FindClustersSDD(AliBin* bins[2], Int_t nMaxBin, Int_t nzBins, Float_t y=c.GetY(),z=c.GetZ(), q=c.GetQ(); y/=q; z/=q; - + // + //Float_t s2 = c.GetSigmaY2()/c.GetQ() - y*y; + // c.SetSigmaY2(s2); + //s2 = c.GetSigmaZ2()/c.GetQ() - z*z; + //c.SetSigmaZ2(s2); + // y=(y-0.5)*fYpitchSDD; y-=fHwSDD; y-=fYoffSDD; //delay ? @@ -684,83 +849,29 @@ FindClustersSDD(AliBin* bins[2], Int_t nMaxBin, Int_t nzBins, z= -z+fZshift[fI]; c.SetY(y); c.SetZ(z); - - c.SetQ(q/20.); //to be consistent with the SSD charges - - if (digits) { - AliBin *b=&bins[s][idx[k]]; - AliITSdigitSDD* d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); - Int_t l0=(d->GetTracks())[0]; - if (l0<0) { - b=&bins[s][idx[k]-1]; - if (b->GetQ()>0) { - d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); - l0=(d->GetTracks())[0]; - } - } - if (l0<0) { - b=&bins[s][idx[k]+1]; - if (b->GetQ()>0) { - d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); - l0=(d->GetTracks())[0]; - } - } - if (l0<0) { - b=&bins[s][idx[k]-nzBins]; - if (b->GetQ()>0) { - d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); - l0=(d->GetTracks())[0]; - } - } - if (l0<0) { - b=&bins[s][idx[k]+nzBins]; - if (b->GetQ()>0) { - d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); - l0=(d->GetTracks())[0]; - } - } - - if (l0<0) { - b=&bins[s][idx[k]+nzBins+1]; - if (b->GetQ()>0) { - d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); - l0=(d->GetTracks())[0]; - } - } - if (l0<0) { - b=&bins[s][idx[k]+nzBins-1]; - if (b->GetQ()>0) { - d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); - l0=(d->GetTracks())[0]; - } - } - if (l0<0) { - b=&bins[s][idx[k]-nzBins+1]; - if (b->GetQ()>0) { - d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); - l0=(d->GetTracks())[0]; - } - } - if (l0<0) { - b=&bins[s][idx[k]-nzBins-1]; - if (b->GetQ()>0) { - d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); - l0=(d->GetTracks())[0]; - } - } - + c.SetNy(maxj-minj+1); + c.SetNz(maxi-mini+1); + c.SetType(npeaks); + c.SetQ(q/12.7); //to be consistent with the SSD charges + + if (c.GetQ() < 20.) continue; //noise cluster + + if (digits) { + // AliBin *b=&bins[s][idx[k]]; + // AliITSdigitSDD* d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); { - Int_t lab[3]; - lab[0]=(d->GetTracks())[0]; - lab[1]=(d->GetTracks())[1]; - lab[2]=(d->GetTracks())[2]; + //Int_t lab[3]; + //lab[0]=(d->GetTracks())[0]; + //lab[1]=(d->GetTracks())[1]; + //lab[2]=(d->GetTracks())[2]; //CheckLabels(lab); - c.SetLabel(lab[0],0); - c.SetLabel(lab[1],1); - c.SetLabel(lab[2],2); + CheckLabels2(milab); + c.SetLabel(milab[0],0); + c.SetLabel(milab[1],1); + c.SetLabel(milab[2],2); + c.SetLayer(fNlayer[fI]); } } - new (cl[ncl]) AliITSclusterV2(c); ncl++; } } @@ -785,6 +896,8 @@ FindClustersSDD(const TClonesArray *digits, TClonesArray *clusters) { Int_t y=d->GetCoord2()+1; //y Int_t z=d->GetCoord1()+1; //z Int_t q=d->GetSignal(); + if (q<3) continue; + if (z <= fNzSDD) { bins[0][y*kNzBins+z].SetQ(q); bins[0][y*kNzBins+z].SetMask(1); @@ -801,6 +914,7 @@ FindClustersSDD(const TClonesArray *digits, TClonesArray *clusters) { delete[] bins[0]; delete[] bins[1]; + } void AliITSclustererV2::FindClustersSDD(AliITSRawStream* input, @@ -812,7 +926,10 @@ void AliITSclustererV2::FindClustersSDD(AliITSRawStream* input, Int_t nClustersSDD = 0; Int_t kNzBins = fNzSDD + 2; Int_t kMaxBin = kNzBins * (fNySDD+2); - AliBin* bins[2] = {NULL, NULL}; + AliBin *binsSDDInit = new AliBin[kMaxBin]; + AliBin *binsSDD1 = new AliBin[kMaxBin]; + AliBin *binsSDD2 = new AliBin[kMaxBin]; + AliBin *bins[2] = {NULL, NULL}; // read raw data input stream while (kTRUE) { @@ -827,28 +944,37 @@ void AliITSclustererV2::FindClustersSDD(AliITSRawStream* input, FindClustersSDD(bins, kMaxBin, kNzBins, NULL, clusters[iModule]); Int_t nClusters = clusters[iModule]->GetEntriesFast(); nClustersSDD += nClusters; - delete[] bins[0]; - delete[] bins[1]; + bins[0] = bins[1] = NULL; } if (!next) break; - bins[0] = new AliBin[kMaxBin]; - bins[1] = new AliBin[kMaxBin]; + bins[0]=binsSDD1; + bins[1]=binsSDD2; + memcpy(binsSDD1,binsSDDInit,sizeof(AliBin)*kMaxBin); + memcpy(binsSDD2,binsSDDInit,sizeof(AliBin)*kMaxBin); } // fill the current digit into the bins array - Int_t iz = input->GetCoord1()+1; - Int_t side = ((iz <= fNzSDD) ? 0 : 1); - iz -= side*fNzSDD; - Int_t index = (input->GetCoord2()+1) * kNzBins + iz; - bins[side][index].SetQ(input->GetSignal()); - bins[side][index].SetMask(1); - bins[side][index].SetIndex(index); + if(input->GetSignal()>=3) { + Int_t iz = input->GetCoord1()+1; + Int_t side = ((iz <= fNzSDD) ? 0 : 1); + iz -= side*fNzSDD; + Int_t index = (input->GetCoord2()+1) * kNzBins + iz; + bins[side][index].SetQ(input->GetSignal()); + bins[side][index].SetMask(1); + bins[side][index].SetIndex(index); + } } + delete[] binsSDD1; + delete[] binsSDD2; + delete[] binsSDDInit; + Info("FindClustersSDD", "found clusters in ITS SDD: %d", nClustersSDD); } + + void AliITSclustererV2:: FindClustersSSD(Ali1Dcluster* neg, Int_t nn, Ali1Dcluster* pos, Int_t np, @@ -857,79 +983,472 @@ FindClustersSSD(Ali1Dcluster* neg, Int_t nn, // Actual SSD cluster finder //------------------------------------------------------------ TClonesArray &cl=*clusters; - - Int_t lab[4]={-2,-2,-2,-2}; + // Float_t tanp=fTanP, tann=fTanN; if (fI>fLastSSD1) {tann=fTanP; tanp=fTanN;} - Int_t idet=fNdet[fI]; Int_t ncl=0; + // + Int_t negativepair[30000]; + Int_t cnegative[3000]; + Int_t cused1[3000]; + Int_t positivepair[30000]; + Int_t cpositive[3000]; + Int_t cused2[3000]; + for (Int_t i=0;i<3000;i++) {cnegative[i]=0; cused1[i]=0;} + for (Int_t i=0;i<3000;i++) {cpositive[i]=0; cused2[i]=0;} + for (Int_t i=0;i<30000;i++){negativepair[i]=positivepair[i]=0;} + static Short_t pairs[1000][1000]; + memset(pairs,0,sizeof(Short_t)*1000000); +// Short_t ** pairs = new Short_t*[1000]; +// for (Int_t i=0; i<1000; i++) { +// pairs[i] = new Short_t[1000]; +// memset(pairs[i],0,sizeof(Short_t)*1000); +// } + // + // find available pairs + // + for (Int_t i=0; iSetChargeRatio(ratio); + cl2->SetType(1); + pairs[ip][j]=1; + if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster + cl2->SetType(2); + pairs[ip][j]=2; + } + cused1[ip]++; + cused2[j]++; + } + } + + for (Int_t ip=0;ipSetChargeRatio(ratio); + cl2->SetType(5); + pairs[ip][in] = 5; + if ((pos[ip].GetNd()+neg[in].GetNd())>6){ //multi cluster + cl2->SetType(6); + pairs[ip][in] = 6; + } + } + // + // add second pair + + // if (!(cused1[ip2] || cused2[in])){ // + if (pairs[ip2][in]==100){ + Float_t yp=pos[ip2].GetY()*fYpitchSSD; + Float_t yn=neg[in].GetY()*fYpitchSSD; + Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp); + Float_t yt=yn + tann*zt; + zt-=fHlSSD; yt-=fHwSSD; + ybest =yt; zbest=zt; + qbest =pos[ip2].GetQ(); + lp[0]=-(-ybest+fYshift[fI]); + lp[1]= -zbest+fZshift[fI]; + lp[2]=0.0025*0.0025; //SigmaY2 + lp[3]=0.110*0.110; //SigmaZ2 + + lp[4]=qbest; //Q + for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2; + for (Int_t ilab=0;ilab<3;ilab++){ + milab[ilab] = pos[ip2].GetLabel(ilab); + milab[ilab+3] = neg[in].GetLabel(ilab); + } + // + CheckLabels2(milab); + ratio = (pos[ip2].GetQ()-neg[in].GetQ())/(pos[ip2].GetQ()+neg[in].GetQ()); + milab[3]=(((ip2<<10) + in)<<10) + idet; // pos|neg|det + Int_t info[3] = {pos[ip2].GetNd(),neg[in].GetNd(),fNlayer[fI]}; + AliITSclusterV2 *cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info); + ncl++; + cl2->SetChargeRatio(ratio); + cl2->SetType(5); + pairs[ip2][in] =5; + if ((pos[ip2].GetNd()+neg[in].GetNd())>6){ //multi cluster + cl2->SetType(6); + pairs[ip2][in] =6; + } + } + cused1[ip]++; + cused1[ip2]++; + cused2[in]++; + } + } + } + + // + for (Int_t jn=0;jnSetChargeRatio(ratio); + cl2->SetType(7); + pairs[ip][jn] =7; + if ((pos[ip].GetNd()+neg[jn].GetNd())>6){ //multi cluster + cl2->SetType(8); + pairs[ip][jn]=8; + } + } + // + // add second pair + // if (!(cused1[ip]||cused2[jn2])){ + if (pairs[ip][jn2]==100){ + Float_t yn=neg[jn2].GetY()*fYpitchSSD; + Double_t yp=pos[ip].GetY()*fYpitchSSD; + Double_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp); + Double_t yt=yn + tann*zt; + zt-=fHlSSD; yt-=fHwSSD; + ybest =yt; zbest=zt; + qbest =neg[jn2].GetQ(); + lp[0]=-(-ybest+fYshift[fI]); + lp[1]= -zbest+fZshift[fI]; + lp[2]=0.0025*0.0025; //SigmaY2 + lp[3]=0.110*0.110; //SigmaZ2 + + lp[4]=qbest; //Q + for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2; + for (Int_t ilab=0;ilab<3;ilab++){ + milab[ilab] = pos[ip].GetLabel(ilab); + milab[ilab+3] = neg[jn2].GetLabel(ilab); + } + // + CheckLabels2(milab); + ratio = (pos[ip].GetQ()-neg[jn2].GetQ())/(pos[ip].GetQ()+neg[jn2].GetQ()); + milab[3]=(((ip<<10) + jn2)<<10) + idet; // pos|neg|det + Int_t info[3] = {pos[ip].GetNd(),neg[jn2].GetNd(),fNlayer[fI]}; + AliITSclusterV2* cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info); + ncl++; + cl2->SetChargeRatio(ratio); + pairs[ip][jn2]=7; + cl2->SetType(7); + if ((pos[ip].GetNd()+neg[jn2].GetNd())>6){ //multi cluster + cl2->SetType(8); + pairs[ip][jn2]=8; + } + } + cused1[ip]++; + cused2[jn]++; + cused2[jn2]++; + } + } + } + + for (Int_t ip=0;ip1) continue; // more than one "proper" cluster for positive + // + count =0; + for (Int_t dj=0;dj1) continue; // more than one "proper" cluster for negative + + Int_t jp = 0; + + count =0; + for (Int_t dj=0;dj1) continue; + if (pairs[ip][j]<100) continue; + // + //almost gold clusters + Float_t yp=pos[ip].GetY()*fYpitchSSD; + Float_t yn=neg[j].GetY()*fYpitchSSD; + Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp); + Float_t yt=yn + tann*zt; + zt-=fHlSSD; yt-=fHwSSD; + ybest=yt; zbest=zt; + qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ()); + lp[0]=-(-ybest+fYshift[fI]); + lp[1]= -zbest+fZshift[fI]; + lp[2]=0.0025*0.0025; //SigmaY2 + lp[3]=0.110*0.110; //SigmaZ2 + lp[4]=qbest; //Q + for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2; + for (Int_t ilab=0;ilab<3;ilab++){ + milab[ilab] = pos[ip].GetLabel(ilab); + milab[ilab+3] = neg[j].GetLabel(ilab); + } + // + CheckLabels2(milab); + ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ()); + milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det + Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fI]}; + AliITSclusterV2 * cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info); + ncl++; + cl2->SetChargeRatio(ratio); + cl2->SetType(10); + pairs[ip][j]=10; + if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster + cl2->SetType(11); + pairs[ip][j]=11; + } + cused1[ip]++; + cused2[j]++; + } + + } + + // for (Int_t i=0; i0 &&pairs[i][j]<100) continue; + ratio = (pos[i].GetQ()-neg[j].GetQ())/(pos[i].GetQ()+neg[j].GetQ()); Float_t yn=neg[j].GetY()*fYpitchSSD; Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp); Float_t yt=yn + tann*zt; zt-=fHlSSD; yt-=fHwSSD; if (TMath::Abs(yt) 4) { - lp[2]*=9; - lp[3]*=9; - } - lp[4]=qbest; //Q - //CheckLabels(lab); - new (cl[ncl]) AliITSclusterV2(lab,lp); ncl++; + lp[4]=qbest; //Q + for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2; + for (Int_t ilab=0;ilab<3;ilab++){ + milab[ilab] = pos[i].GetLabel(ilab); + milab[ilab+3] = neg[j].GetLabel(ilab); + } + // + CheckLabels2(milab); + milab[3]=(((i<<10) + j)<<10) + idet; // pos|neg|det + Int_t info[3] = {pos[i].GetNd(),neg[j].GetNd(),fNlayer[fI]}; + AliITSclusterV2 * cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info); + ncl++; + cl2->SetChargeRatio(ratio); + cl2->SetType(100+cpositive[j]+cnegative[i]); + //cl2->SetType(0); + /* + if (pairs[i][j]<100){ + printf("problem:- %d\n", pairs[i][j]); + } + if (cnegative[i]<2&&cpositive[j]<2){ + printf("problem:- %d\n", pairs[i][j]); + } + */ } } - /* - if (ybest<100) { - lab[3]=idet; - Float_t lp[5]; - lp[0]=-ybest-fYshift[fI]; - lp[1]= zbest+fZshift[fI]; - lp[2]=0.002*0.002; //SigmaY2 - lp[3]=0.080*0.080; //SigmaZ2 - lp[4]=qbest; //Q - // - new (cl[ncl]) AliITSclusterV2(lab,lp); ncl++; - } - */ } + +// for (Int_t i=0; i<1000; i++) delete [] pairs[i]; +// delete [] pairs; + } + void AliITSclustererV2:: -FindClustersSSD(const TClonesArray *digits, TClonesArray *clusters) { +FindClustersSSD(const TClonesArray *alldigits, TClonesArray *clusters) { //------------------------------------------------------------ // Actual SSD cluster finder //------------------------------------------------------------ - Int_t smax=digits->GetEntriesFast(); + Int_t smaxall=alldigits->GetEntriesFast(); + if (smaxall==0) return; + TObjArray *digits = new TObjArray; + for (Int_t i=0;iUncheckedAt(i); + if (d->GetSignal()<3) continue; + digits->AddLast(d); + } + Int_t smax = digits->GetEntriesFast(); if (smax==0) return; - + const Int_t MAX=1000; Int_t np=0, nn=0; Ali1Dcluster pos[MAX], neg[MAX]; Float_t y=0., q=0., qmax=0.; Int_t lab[4]={-2,-2,-2,-2}; - + AliITSdigitSSD *d=(AliITSdigitSSD*)digits->UncheckedAt(0); q += d->GetSignal(); y += d->GetCoord2()*d->GetSignal(); @@ -940,28 +1459,60 @@ FindClustersSSD(const TClonesArray *digits, TClonesArray *clusters) { Int_t *n=&nn; Ali1Dcluster *c=neg; Int_t nd=1; + Int_t milab[10]; + for (Int_t ilab=0;ilab<10;ilab++){ + milab[ilab]=-2; + } + milab[0]=d->GetTrack(0); milab[1]=d->GetTrack(1); milab[2]=d->GetTrack(2); + for (Int_t s=1; sUncheckedAt(s); + d=(AliITSdigitSSD*)digits->UncheckedAt(s); Int_t strip=d->GetCoord2(); if ((strip-curr) > 1 || flag!=d->GetCoord1()) { c[*n].SetY(y/q); c[*n].SetQ(q); c[*n].SetNd(nd); - c[*n].SetLabels(lab); + CheckLabels2(milab); + c[*n].SetLabels(milab); //Split suspiciously big cluster - if (nd>3) { - c[*n].SetY(y/q-0.5*nd); - c[*n].SetQ(0.5*q); - (*n)++; - if (*n==MAX) { - Error("FindClustersSSD","Too many 1D clusters !"); + /* + if (nd>10&&nd<16){ + c[*n].SetY(y/q-0.3*nd); + c[*n].SetQ(0.5*q); + (*n)++; + if (*n==MAX) { + Error("FindClustersSSD","Too many 1D clusters !"); return; - } - c[*n].SetY(y/q+0.5*nd); - c[*n].SetQ(0.5*q); - c[*n].SetNd(nd); - c[*n].SetLabels(lab); - } + } + c[*n].SetY(y/q-0.0*nd); + c[*n].SetQ(0.5*q); + c[*n].SetNd(nd); + (*n)++; + if (*n==MAX) { + Error("FindClustersSSD","Too many 1D clusters !"); + return; + } + // + c[*n].SetY(y/q+0.3*nd); + c[*n].SetQ(0.5*q); + c[*n].SetNd(nd); + c[*n].SetLabels(milab); + } + else{ + */ + if (nd>4&&nd<25) { + c[*n].SetY(y/q-0.25*nd); + c[*n].SetQ(0.5*q); + (*n)++; + if (*n==MAX) { + Error("FindClustersSSD","Too many 1D clusters !"); + return; + } + c[*n].SetY(y/q+0.25*nd); + c[*n].SetQ(0.5*q); + c[*n].SetNd(nd); + c[*n].SetLabels(milab); + } (*n)++; if (*n==MAX) { Error("FindClustersSSD","Too many 1D clusters !"); @@ -970,6 +1521,11 @@ FindClustersSSD(const TClonesArray *digits, TClonesArray *clusters) { y=q=qmax=0.; nd=0; lab[0]=lab[1]=lab[2]=-2; + // + for (Int_t ilab=0;ilab<10;ilab++){ + milab[ilab]=-2; + } + // if (flag!=d->GetCoord1()) { n=&np; c=pos; } } flag=d->GetCoord1(); @@ -980,6 +1536,9 @@ FindClustersSSD(const TClonesArray *digits, TClonesArray *clusters) { qmax=d->GetSignal(); lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2); } + for (Int_t ilab=0;ilab<10;ilab++) { + if (d->GetTrack(ilab)>=0) AddLabel(milab, (d->GetTrack(ilab))); + } curr=strip; } c[*n].SetY(y/q); @@ -987,15 +1546,15 @@ FindClustersSSD(const TClonesArray *digits, TClonesArray *clusters) { c[*n].SetNd(nd); c[*n].SetLabels(lab); //Split suspiciously big cluster - if (nd>3) { - c[*n].SetY(y/q-0.5*nd); + if (nd>4 && nd<25) { + c[*n].SetY(y/q-0.25*nd); c[*n].SetQ(0.5*q); (*n)++; if (*n==MAX) { Error("FindClustersSSD","Too many 1D clusters !"); return; } - c[*n].SetY(y/q+0.5*nd); + c[*n].SetY(y/q+0.25*nd); c[*n].SetQ(0.5*q); c[*n].SetNd(nd); c[*n].SetLabels(lab); @@ -1019,20 +1578,23 @@ void AliITSclustererV2::FindClustersSSD(AliITSRawStream* input, const Int_t MAX = 1000; Ali1Dcluster clusters1D[2][MAX]; Int_t nClusters[2] = {0, 0}; + Int_t lab[3]={-2,-2,-2}; Float_t q = 0.; Float_t y = 0.; Int_t nDigits = 0; Int_t prevStrip = -1; Int_t prevFlag = -1; + Int_t prevModule = -1; // read raw data input stream while (kTRUE) { Bool_t next = input->Next(); + if(input->GetSignal()<3 && next) continue; // check if a new cluster starts Int_t strip = input->GetCoord2(); Int_t flag = input->GetCoord1(); - if ((!next || input->IsNewModule() || + if ((!next || (input->GetModuleID() != prevModule)|| (strip-prevStrip > 1) || (flag != prevFlag)) && (nDigits > 0)) { if (nClusters[prevFlag] == MAX) { @@ -1043,26 +1605,28 @@ void AliITSclustererV2::FindClustersSSD(AliITSRawStream* input, cluster.SetY(y/q); cluster.SetQ(q); cluster.SetNd(nDigits); + cluster.SetLabels(lab); //Split suspiciously big cluster - if (nDigits > 3) { - cluster.SetY(y/q - 0.5*nDigits); + if (nDigits > 4&&nDigits < 25) { + cluster.SetY(y/q - 0.25*nDigits); cluster.SetQ(0.5*q); if (nClusters[prevFlag] == MAX) { Error("FindClustersSSD", "Too many 1D clusters !"); return; } Ali1Dcluster& cluster2 = clusters1D[prevFlag][nClusters[prevFlag]++]; - cluster2.SetY(y/q + 0.5*nDigits); + cluster2.SetY(y/q + 0.25*nDigits); cluster2.SetQ(0.5*q); cluster2.SetNd(nDigits); + cluster2.SetLabels(lab); } y = q = 0.; nDigits = 0; } - if (!next || input->IsNewModule()) { - Int_t iModule = input->GetPrevModuleID(); + if (!next || (input->GetModuleID() != prevModule)) { + Int_t iModule = prevModule; // when all data from a module was read, search for clusters if (prevFlag >= 0) { @@ -1086,6 +1650,8 @@ void AliITSclustererV2::FindClustersSSD(AliITSRawStream* input, nDigits++; prevStrip = strip; prevFlag = flag; + prevModule = input->GetModuleID(); + } Info("FindClustersSSD", "found clusters in ITS SSD: %d", nClustersSSD); @@ -1094,7 +1660,7 @@ void AliITSclustererV2::FindClustersSSD(AliITSRawStream* input, #else //V1 #include "AliITSDetType.h" - +#include "AliITS.h" #include "AliITSsegmentationSPD.h" #include "AliITSClusterFinderSPD.h" @@ -1141,7 +1707,7 @@ FindClustersSDD(const TClonesArray *digits, TClonesArray *clusters) { Int_t nc=points->GetEntriesFast(); while (nc--) { //To be consistent with the SSD cluster charges AliITSRecPoint *p=(AliITSRecPoint*)points->UncheckedAt(nc); - p->SetQ(p->GetQ/12.); + p->SetQ(p->GetQ()/12.); } RecPoints2Clusters(points, fI, clusters); its->ResetClusters(1);