X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=ITS%2FAliITSclustererV2.cxx;h=2e9bb92dc78e8c2d0685f6d94e6016363174fa65;hb=f60271ca17d635bd65a3f2e05b535da268b7dfe7;hp=19eec7d015fe7486dfbc40ed49a1c8bfd93f6cec;hpb=7f4044f0ab63faa935ed50dc5c54824439c59227;p=u%2Fmrichter%2FAliRoot.git diff --git a/ITS/AliITSclustererV2.cxx b/ITS/AliITSclustererV2.cxx index 19eec7d015f..2e9bb92dc78 100644 --- a/ITS/AliITSclustererV2.cxx +++ b/ITS/AliITSclustererV2.cxx @@ -10,13 +10,19 @@ #include "AliITSclustererV2.h" #include "AliITSclusterV2.h" +#include "AliRawReader.h" +#include "AliITSRawStreamSPD.h" +#include "AliITSRawStreamSDD.h" +#include "AliITSRawStreamSSD.h" -#include #include #include #include #include "AliITSgeom.h" -#include "AliITSdigit.h" +#include "AliITSdigitSPD.h" +#include "AliITSdigitSDD.h" +#include "AliITSdigitSSD.h" +#include "AliMC.h" ClassImp(AliITSclustererV2) @@ -32,16 +38,22 @@ AliITSclustererV2::AliITSclustererV2(const AliITSgeom *geom) { fI=0; Int_t mmax=geom->GetIndexMax(); - if (mmax>2200) {cerr<<"Too many ITS subdetectors !\n"; exit(1);} + if (mmax>2200) { + Fatal("AliITSclustererV2","Too many ITS subdetectors !"); + } Int_t m; for (m=0; mGetModuleId(m,lay,lad,det); Float_t x,y,z; g->GetTrans(lay,lad,det,x,y,z); Double_t rot[9]; g->GetRotMatrix(lay,lad,det,rot); - fYshift[m] = x*rot[0] + y*rot[1]; + Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi(); + Double_t ca=TMath::Cos(alpha), sa=TMath::Sin(alpha); + 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(); //SPD geometry fLastSPD1=g->GetModuleIndex(2,1,1)-1; @@ -79,28 +91,19 @@ AliITSclustererV2::AliITSclustererV2(const AliITSgeom *geom) { fHlSSD=2.00; fTanP=0.0275; fTanN=0.0075; + } -void AliITSclustererV2::Digits2Clusters(const TFile *in, TFile *out) { + +Int_t AliITSclustererV2::Digits2Clusters(TTree *dTree, TTree *cTree) { //------------------------------------------------------------ // This function creates ITS clusters //------------------------------------------------------------ Int_t ncl=0; - TDirectory *savedir=gDirectory; - - if (!out->IsOpen()) { - cerr<<"AliITSclustererV2::Digits2Clusters(): output file not open !\n"; - return; - } - - Char_t name[100]; - sprintf(name,"TreeD%d",fEvent); - //TTree *dTree=(TTree*)((TFile*)in)->Get(name); - TTree *dTree=gAlice->TreeD(); if (!dTree) { - cerr<<"Input tree "<SetBranchAddress("ITSDigitsSSD",&digitsSSD); - Int_t mmax=(Int_t)dTree->GetEntries(); - - out->cd(); - - sprintf(name,"TreeC_ITS_%d",fEvent); - TTree cTree(name,"ITS clusters V2"); TClonesArray *clusters=new TClonesArray("AliITSclusterV2",1000); - cTree.Branch("Clusters",&clusters); + TBranch *branch=cTree->GetBranch("Clusters"); + if (!branch) cTree->Branch("Clusters",&clusters); + else branch->SetAddress(&clusters); + + Int_t mmax=(Int_t)dTree->GetEntries(); 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(); + cTree->Fill(); digitsSPD->Clear(); digitsSDD->Clear(); digitsSSD->Clear(); clusters->Clear(); } - cTree.Write(); + + //cTree->Write(); delete clusters; @@ -148,9 +151,70 @@ void AliITSclustererV2::Digits2Clusters(const TFile *in, TFile *out) { //delete dTree; - cerr<<"Number of found clusters : "<cd(); + return 0; +} + +void AliITSclustererV2::Digits2Clusters(AliRawReader* rawReader) { + //------------------------------------------------------------ + // This function creates ITS clusters from raw data + //------------------------------------------------------------ + AliRunLoader* runLoader = AliRunLoader::GetRunLoader(); + if (!runLoader) { + Error("Digits2Clusters", "no run loader found"); + return; + } + runLoader->LoadKinematics(); + AliLoader* itsLoader = runLoader->GetLoader("ITSLoader"); + if (!itsLoader) { + Error("Digits2Clusters", "no loader for ITS found"); + return; + } + if (!itsLoader->TreeR()) itsLoader->MakeTree("R"); + TTree* cTree = itsLoader->TreeR(); + + TClonesArray *array=new TClonesArray("AliITSclusterV2",1000); + cTree->Branch("Clusters",&array); + delete array; + + TClonesArray** clusters = new TClonesArray*[fNModules]; + for (Int_t iModule = 0; iModule < fNModules; iModule++) { + clusters[iModule] = NULL; + } + // one TClonesArray per module + + rawReader->Reset(); + AliITSRawStreamSPD inputSPD(rawReader); + FindClustersSPD(&inputSPD, clusters); + + rawReader->Reset(); + AliITSRawStreamSDD inputSDD(rawReader); + FindClustersSDD(&inputSDD, clusters); + + rawReader->Reset(); + AliITSRawStreamSSD inputSSD(rawReader); + FindClustersSSD(&inputSSD, clusters); + + // write all clusters to the tree + Int_t nClusters = 0; + for (Int_t iModule = 0; iModule < fNModules; iModule++) { + array = clusters[iModule]; + if (!array) { + Error("Digits2Clusters", "data for module %d missing!", iModule); + array = new TClonesArray("AliITSclusterV2"); + } + cTree->SetBranchAddress("Clusters", &array); + cTree->Fill(); + nClusters += array->GetEntriesFast(); + delete array; + } + itsLoader->WriteRecPoints("OVERWRITE"); + + delete[] clusters; + + Info("Digits2Clusters", "total number of found clusters in ITS: %d\n", + nClusters); } //**** Fast clusters ******************************* @@ -162,83 +226,194 @@ void AliITSclustererV2::Digits2Clusters(const TFile *in, TFile *out) { #include "AliITSsimulationFastPoints.h" #include "AliITSRecPoint.h" +/* static void CheckLabels(Int_t lab[3]) { //------------------------------------------------------------ // Tries to find mother's labels //------------------------------------------------------------ Int_t label=lab[0]; if (label>=0) { - TParticle *part=(TParticle*)gAlice->Particle(label); + TParticle *part=(TParticle*)gAlice->GetMCApp()->Particle(label); label=-3; while (part->P() < 0.005) { Int_t m=part->GetFirstMother(); - if (m<0) {cerr<<"Primary momentum: "<P()<P()); + break; + } + if (part->GetStatusCode()>0) { + Info("CheckLabels","Primary momentum: %f",part->P()); + break; + } label=m; - part=(TParticle*)gAlice->Particle(label); + part=(TParticle*)gAlice->GetMCApp()->Particle(label); } if (lab[1]<0) lab[1]=label; else if (lab[2]<0) lab[2]=label; else ;//cerr<<"CheckLabels : No empty labels !\n"; } } +*/ +static void CheckLabels(Int_t lab[3]) { + //------------------------------------------------------------ + // Tries to find mother's labels + //------------------------------------------------------------ -static void RecPoints2Clusters(const TClonesArray *points, -Double_t ys, Double_t zs, Int_t idx, Int_t sgn, TClonesArray *clusters) { + 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) { //------------------------------------------------------------ - // Conversion AliITSRecPoint -> AliITSclusterV2 + // Conversion AliITSRecPoint -> AliITSclusterV2 for the ITS + // subdetector indexed by idx //------------------------------------------------------------ TClonesArray &cl=*clusters; Int_t ncl=points->GetEntriesFast(); for (Int_t i=0; iUncheckedAt(i); Float_t lp[5]; - lp[0]=-p->GetX()-ys; lp[0]*=sgn; //SPD1 - lp[1]=p->GetZ()+zs; - lp[2]=p->GetSigmaX2()*0.5; - lp[3]=p->GetSigmaZ2()*0.5; - lp[4]=p->GetQ()/12.; //to be roughly consistent with SSD + lp[0]=-(-p->GetX()+fYshift[idx]); if (idx<=fLastSPD1) lp[0]*=-1; //SPD1 + lp[1]= -p->GetZ()+fZshift[idx]; + lp[2]=p->GetSigmaX2(); + 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]=idx; + 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); } } -void AliITSclustererV2::Hits2Clusters(const TFile *in, TFile *out) { +Int_t AliITSclustererV2::Hits2Clusters(TTree *hTree, TTree *cTree) { //------------------------------------------------------------ // This function creates ITS clusters //------------------------------------------------------------ - TDirectory *savedir=gDirectory; - - if (!out->IsOpen()) { - cerr<<"AliITSclustererV2::Hits2Clusters: output file not open !\n"; - return; - } - if (!gAlice) { - cerr<<"AliITSclustererV2::Hits2Clusters : gAlice==0 !\n"; - return; + Error("Hits2Clusters","gAlice==0 !"); + return 1; } AliITS *its = (AliITS*)gAlice->GetModule("ITS"); if (!its) { - cerr<<"AliITSclustererV2::Hits2Clusters : Can't find the ITS !\n"; - return; + Error("Hits2Clusters","Can't find the ITS !"); + return 2; } AliITSgeom *geom=its->GetITSgeom(); Int_t mmax=geom->GetIndexMax(); its->InitModules(-1,mmax); - its->FillModules(gAlice->TreeH(),0); - - out->cd(); + its->FillModules(hTree,0); - Char_t name[100]; - sprintf(name,"TreeC_ITS_%d",fEvent); - TTree cTree(name,"ITS clusters V2"); TClonesArray *clusters=new TClonesArray("AliITSclusterV2",1000); - cTree.Branch("Clusters",&clusters); + TBranch *branch=cTree->GetBranch("Clusters"); + if (!branch) cTree->Branch("Clusters",&clusters); + else branch->SetAddress(&clusters); static TClonesArray *points=its->RecPoints(); AliITSsimulationFastPoints sim; @@ -247,21 +422,21 @@ void AliITSclustererV2::Hits2Clusters(const TFile *in, TFile *out) { AliITSmodule *mod=its->GetModule(m); sim.CreateFastRecPoints(mod,m,gRandom); - Int_t s = (m<=fLastSPD1) ? -1 : 1; - RecPoints2Clusters(points, fYshift[m], fZshift[m], fNdet[m], s, clusters); + RecPoints2Clusters(points, m, clusters); its->ResetRecPoints(); ncl+=clusters->GetEntriesFast(); - cTree.Fill(); + cTree->Fill(); clusters->Clear(); } - cTree.Write(); - cerr<<"Number of found fast clusters : "<Write(); delete clusters; - savedir->cd(); + return 0; } //*********************************** @@ -294,7 +469,8 @@ FindClustersSPD(const TClonesArray *digits, TClonesArray *clusters) { //------------------------------------------------------------ // Actual SPD cluster finder //------------------------------------------------------------ - const Int_t kMAXBIN=(fNzSPD+2)*(fNySPD+2); + Int_t kNzBins = fNzSPD + 2; + const Int_t kMAXBIN=kNzBins*(fNySPD+2); Int_t ndigits=digits->GetEntriesFast(); AliBin *bins=new AliBin[kMAXBIN]; @@ -305,27 +481,28 @@ FindClustersSPD(const TClonesArray *digits, TClonesArray *clusters) { d=(AliITSdigitSPD*)digits->UncheckedAt(k); Int_t i=d->GetCoord2()+1; //y Int_t j=d->GetCoord1()+1; - bins[i*fNzSPD+j].SetIndex(k); - bins[i*fNzSPD+j].SetMask(1); + bins[i*kNzBins+j].SetIndex(k); + bins[i*kNzBins+j].SetMask(1); } Int_t n=0; TClonesArray &cl=*clusters; for (k=0; kUncheckedAt(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]); @@ -333,39 +510,191 @@ 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; + 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; +} + +void AliITSclustererV2::FindClustersSPD(AliITSRawStream* input, + TClonesArray** clusters) +{ + //------------------------------------------------------------ + // Actual SPD cluster finder for raw data + //------------------------------------------------------------ + + Int_t nClustersSPD = 0; + Int_t kNzBins = fNzSPD + 2; + Int_t kNyBins = fNySPD + 2; + Int_t kMaxBin = kNzBins * kNyBins; + AliBin *binsSPD = new AliBin[kMaxBin]; + AliBin *binsSPDInit = new AliBin[kMaxBin]; + AliBin *bins = NULL; + + // read raw data input stream + while (kTRUE) { + Bool_t next = input->Next(); + if (!next || input->IsNewModule()) { + Int_t iModule = input->GetPrevModuleID(); + + // when all data from a module was read, search for clusters + if (bins) { + clusters[iModule] = new TClonesArray("AliITSclusterV2"); + Int_t nClusters = 0; + + for (Int_t iBin = 0; iBin < kMaxBin; iBin++) { + if (bins[iBin].IsUsed()) continue; + Int_t nBins = 0; + Int_t idxBins[200]; + FindCluster(iBin, kNzBins, bins, nBins, idxBins); + if (nBins == 200) { + Error("FindClustersSPD", "SPD: Too big cluster !\n"); + continue; + } + + Int_t label[4]; + label[0] = -2; + label[1] = -2; + label[2] = -2; +// label[3] = iModule; + label[3] = fNdet[iModule]; + + Int_t ymin = (idxBins[0] / kNzBins) - 1; + Int_t ymax = ymin; + Int_t zmin = (idxBins[0] % kNzBins) - 1; + Int_t zmax = zmin; + for (Int_t idx = 0; idx < nBins; idx++) { + Int_t iy = (idxBins[idx] / kNzBins) - 1; + Int_t iz = (idxBins[idx] % kNzBins) - 1; + if (ymin > iy) ymin = iy; + if (ymax < iy) ymax = iy; + if (zmin > iz) zmin = iz; + if (zmax < iz) zmax = iz; + } + + 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; + bins = NULL; + } - 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); + if (!next) break; + bins = binsSPD; + memcpy(binsSPD,binsSPDInit,sizeof(AliBin)*kMaxBin); + } - //CheckLabels(lab); - new (cl[n]) AliITSclusterV2(lab,lp); n++; + // fill the current digit into the bins array + Int_t index = (input->GetCoord2()+1) * kNzBins + (input->GetCoord1()+1); + bins[index].SetIndex(index); + bins[index].SetMask(1); + bins[index].SetQ(1); } - delete bins; + delete [] binsSPDInit; + delete [] binsSPD; + + Info("FindClustersSPD", "found clusters in ITS SPD: %d", nClustersSPD); } + Bool_t AliITSclustererV2::IsMaximum(Int_t k,Int_t max,const AliBin *bins) { //------------------------------------------------------------ //is this a local maximum ? @@ -442,51 +771,28 @@ MakeCluster(Int_t k,Int_t max,AliBin *bins,UInt_t m,AliITSclusterV2 &c) { } void AliITSclustererV2:: -FindClustersSDD(const TClonesArray *digits, TClonesArray *clusters) { +FindClustersSDD(AliBin* bins[2], Int_t nMaxBin, Int_t nzBins, + const TClonesArray *digits, TClonesArray *clusters) { //------------------------------------------------------------ // Actual SDD cluster finder //------------------------------------------------------------ - const Int_t kMAXBIN=(fNzSDD+2)*(fNySDD+2); - - AliBin *bins[2]; - bins[0]=new AliBin[kMAXBIN]; - bins[1]=new AliBin[kMAXBIN]; - - AliITSdigitSDD *d=0; - Int_t i, ndigits=digits->GetEntriesFast(); - for (i=0; iUncheckedAt(i); - Int_t y=d->GetCoord2()+1; //y - Int_t z=d->GetCoord1()+1; //z - Int_t q=d->GetSignal(); - if (z <= fNzSDD) { - bins[0][y*fNzSDD+z].SetQ(q); - bins[0][y*fNzSDD+z].SetMask(1); - bins[0][y*fNzSDD+z].SetIndex(i); - } else { - z-=fNzSDD; - bins[1][y*fNzSDD+z].SetQ(q); - bins[1][y*fNzSDD+z].SetMask(1); - bins[1][y*fNzSDD+z].SetIndex(i); - } - } - Int_t ncl=0; TClonesArray &cl=*clusters; for (Int_t s=0; s<2; s++) - for (i=0; i30) continue; + if (npeaks==0) continue; Int_t k,l; for (k=0; k1 || dj>1) continue; @@ -502,32 +808,64 @@ FindClustersSDD(const TClonesArray *digits, TClonesArray *clusters) { } 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); - } + MakeCluster(idx[k], nzBins, bins[s], msk[k], c); + //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); @@ -536,7 +874,12 @@ FindClustersSDD(const TClonesArray *digits, TClonesArray *clusters) { 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 ? @@ -545,109 +888,609 @@ FindClustersSDD(const TClonesArray *digits, TClonesArray *clusters) { z=(z-0.5)*fZpitchSDD; z-=fHlSDD; - y=-y-fYshift[fI]; - z= z+fZshift[fI]; + y=-(-y+fYshift[fI]); + z= -z+fZshift[fI]; c.SetY(y); c.SetZ(z); - - c.SetQ(q/20.); //to be consistent with the SSD charges - - AliBin *b=&bins[s][idx[k]]; - 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]; + 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]; + //CheckLabels(lab); + CheckLabels2(milab); + c.SetLabel(milab[0],0); + c.SetLabel(milab[1],1); + c.SetLabel(milab[2],2); + c.SetLayer(fNlayer[fI]); } - } - if (l0<0) { - b=&bins[s][idx[k]-fNzSDD]; - if (b->GetQ()>0) { - d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); - l0=(d->GetTracks())[0]; - } - } - if (l0<0) { - b=&bins[s][idx[k]+fNzSDD]; - if (b->GetQ()>0) { - d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); - l0=(d->GetTracks())[0]; - } - } - - if (l0<0) { - b=&bins[s][idx[k]+fNzSDD+1]; - if (b->GetQ()>0) { - d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); - l0=(d->GetTracks())[0]; - } - } - if (l0<0) { - b=&bins[s][idx[k]+fNzSDD-1]; - if (b->GetQ()>0) { - d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); - l0=(d->GetTracks())[0]; - } - } - if (l0<0) { - b=&bins[s][idx[k]-fNzSDD+1]; - if (b->GetQ()>0) { - d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); - l0=(d->GetTracks())[0]; - } - } - if (l0<0) { - b=&bins[s][idx[k]-fNzSDD-1]; - if (b->GetQ()>0) { - d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); - l0=(d->GetTracks())[0]; - } - } - - { - 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); - } - + } new (cl[ncl]) AliITSclusterV2(c); ncl++; } } +} + +void AliITSclustererV2:: +FindClustersSDD(const TClonesArray *digits, TClonesArray *clusters) { + //------------------------------------------------------------ + // Actual SDD cluster finder + //------------------------------------------------------------ + Int_t kNzBins = fNzSDD + 2; + const Int_t kMAXBIN=kNzBins*(fNySDD+2); + + AliBin *bins[2]; + bins[0]=new AliBin[kMAXBIN]; + bins[1]=new AliBin[kMAXBIN]; + + AliITSdigitSDD *d=0; + Int_t i, ndigits=digits->GetEntriesFast(); + for (i=0; iUncheckedAt(i); + 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); + bins[0][y*kNzBins+z].SetIndex(i); + } else { + z-=fNzSDD; + bins[1][y*kNzBins+z].SetQ(q); + bins[1][y*kNzBins+z].SetMask(1); + bins[1][y*kNzBins+z].SetIndex(i); + } + } + + FindClustersSDD(bins, kMAXBIN, kNzBins, digits, clusters); delete[] bins[0]; delete[] bins[1]; + +} + +void AliITSclustererV2::FindClustersSDD(AliITSRawStream* input, + TClonesArray** clusters) +{ + //------------------------------------------------------------ + // Actual SDD cluster finder for raw data + //------------------------------------------------------------ + Int_t nClustersSDD = 0; + Int_t kNzBins = fNzSDD + 2; + Int_t kMaxBin = kNzBins * (fNySDD+2); + 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) { + Bool_t next = input->Next(); + if (!next || input->IsNewModule()) { + Int_t iModule = input->GetPrevModuleID(); + + // when all data from a module was read, search for clusters + if (bins[0]) { + clusters[iModule] = new TClonesArray("AliITSclusterV2"); + fI = iModule; + FindClustersSDD(bins, kMaxBin, kNzBins, NULL, clusters[iModule]); + Int_t nClusters = clusters[iModule]->GetEntriesFast(); + nClustersSDD += nClusters; + bins[0] = bins[1] = NULL; + } + + if (!next) break; + 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 + 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(const TClonesArray *digits, TClonesArray *clusters) { +FindClustersSSD(Ali1Dcluster* neg, Int_t nn, + Ali1Dcluster* pos, Int_t np, + TClonesArray *clusters) { //------------------------------------------------------------ // Actual SSD cluster finder //------------------------------------------------------------ - Int_t smax=digits->GetEntriesFast(); - if (smax==0) return; + TClonesArray &cl=*clusters; + // + 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;} + 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)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]); + } + */ + } + } + } + +// for (Int_t i=0; i<1000; i++) delete [] pairs[i]; +// delete [] pairs; + +} + + +void AliITSclustererV2:: +FindClustersSSD(const TClonesArray *alldigits, TClonesArray *clusters) { + //------------------------------------------------------------ + // Actual SSD cluster finder + //------------------------------------------------------------ + 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}; - - TClonesArray &cl=*clusters; - + AliITSdigitSSD *d=(AliITSdigitSSD*)digits->UncheckedAt(0); q += d->GetSignal(); y += d->GetCoord2()*d->GetSignal(); @@ -657,38 +1500,74 @@ FindClustersSSD(const TClonesArray *digits, TClonesArray *clusters) { Int_t flag=d->GetCoord1(); Int_t *n=&nn; Ali1Dcluster *c=neg; - Int_t nd=0; + 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) { - cerr<< - "AliITSclustererV2::FindClustersSSD: Too many 1D clusters !\n"; + /* + 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) { - cerr<<"AliITSclustererV2::FindClustersSSD: Too many 1D clusters !\n"; + Error("FindClustersSSD","Too many 1D clusters !"); return; } 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(); @@ -699,6 +1578,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); @@ -706,80 +1588,115 @@ 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) { - cerr<<"AliITSclustererV2::FindClustersSSD: Too many 1D clusters !\n"; + 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); } (*n)++; if (*n==MAX) { - cerr<<"AliITSclustererV2::FindClustersSSD: Too many 1D clusters !\n"; + Error("FindClustersSSD","Too many 1D clusters !"); return; } - Float_t tanp=fTanP, tann=fTanN; - if (fI>fLastSSD1) {tann=fTanP; tanp=fTanN;} - - Int_t idet=fNdet[fI]; - Int_t ncl=0; - for (Int_t i=0; i 4) { - lp[2]*=9; - lp[3]*=9; - } - lp[4]=qbest; //Q + FindClustersSSD(neg, nn, pos, np, clusters); +} - //CheckLabels(lab); - new (cl[ncl]) AliITSclusterV2(lab,lp); ncl++; +void AliITSclustererV2::FindClustersSSD(AliITSRawStream* input, + TClonesArray** clusters) +{ + //------------------------------------------------------------ + // Actual SSD cluster finder for raw data + //------------------------------------------------------------ + Int_t nClustersSSD = 0; + 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->GetModuleID() != prevModule)|| + (strip-prevStrip > 1) || (flag != prevFlag)) && + (nDigits > 0)) { + if (nClusters[prevFlag] == MAX) { + Error("FindClustersSSD", "Too many 1D clusters !"); + return; + } + Ali1Dcluster& cluster = clusters1D[prevFlag][nClusters[prevFlag]++]; + cluster.SetY(y/q); + cluster.SetQ(q); + cluster.SetNd(nDigits); + cluster.SetLabels(lab); + + //Split suspiciously big cluster + 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.25*nDigits); + cluster2.SetQ(0.5*q); + cluster2.SetNd(nDigits); + cluster2.SetLabels(lab); } + y = q = 0.; + nDigits = 0; } - /* - 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++; + + if (!next || (input->GetModuleID() != prevModule)) { + Int_t iModule = prevModule; + + // when all data from a module was read, search for clusters + if (prevFlag >= 0) { + clusters[iModule] = new TClonesArray("AliITSclusterV2"); + fI = iModule; + FindClustersSSD(&clusters1D[0][0], nClusters[0], + &clusters1D[1][0], nClusters[1], clusters[iModule]); + Int_t nClusters = clusters[iModule]->GetEntriesFast(); + nClustersSSD += nClusters; + } + + if (!next) break; + nClusters[0] = nClusters[1] = 0; + y = q = 0.; + nDigits = 0; } - */ + + // add digit to current cluster + q += input->GetSignal(); + y += strip * input->GetSignal(); + nDigits++; + prevStrip = strip; + prevFlag = flag; + prevModule = input->GetModuleID(); + } + + Info("FindClustersSSD", "found clusters in ITS SSD: %d", nClustersSSD); } #else //V1 @@ -809,8 +1726,7 @@ FindClustersSPD(const TClonesArray *digits, TClonesArray *clusters) { static AliITSClusterFinderSPD cf(seg, (TClonesArray*)digits, points); cf.FindRawClusters(fI); - Int_t s = (fI<=fLastSPD1) ? -1 : 1; - RecPoints2Clusters(points, fYshift[fI], fZshift[fI], fNdet[fI], s, clusters); + RecPoints2Clusters(points, fI, clusters); its->ResetRecPoints(); } @@ -830,7 +1746,12 @@ FindClustersSDD(const TClonesArray *digits, TClonesArray *clusters) { cf(seg,resp,(TClonesArray*)digits,its->ClustersAddress(1)); cf.FindRawClusters(fI); - RecPoints2Clusters(points, fYshift[fI], fZshift[fI], fNdet[fI], 1, 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.); + } + RecPoints2Clusters(points, fI, clusters); its->ResetClusters(1); its->ResetRecPoints(); @@ -848,7 +1769,7 @@ FindClustersSSD(const TClonesArray *digits, TClonesArray *clusters) { static AliITSClusterFinderSSD cf(seg,(TClonesArray*)digits); cf.FindRawClusters(fI); - RecPoints2Clusters(points, fYshift[fI], fZshift[fI], fNdet[fI], 1, clusters); + RecPoints2Clusters(points, fI, clusters); its->ResetRecPoints(); }