+ //Int_t nPads = nCols * nRows;
+ // This is what we are interested in: The center of gravity of the best candidates
+ Float_t cogyvals[kMaxPads]; memset(cogyvals, 0, sizeof(Float_t)*kMaxPads);
+ Float_t cogzvals[kMaxPads]; memset(cogzvals, 0, sizeof(Float_t)*kMaxPads);
+ Float_t *cogy[kMaxRows];
+ Float_t *cogz[kMaxRows];
+
+ // Lookup-Table storing coordinates according to the bins
+ Float_t yLengths[kMaxCols]; memset(yLengths, 0, kMaxCols*sizeof(Float_t));
+ Float_t zLengths[kMaxRows]; memset(zLengths, 0, kMaxRows*sizeof(Float_t));
+ for(Int_t icnt = 0; icnt < nCols; icnt++){
+ yLengths[icnt] = pp->GetColPos(nCols - 1 - icnt) + binlength/2;
+ }
+ for(Int_t icnt = 0; icnt < nRows; icnt++){
+ zLengths[icnt] = pp->GetRowPos(icnt) - pp->GetRowSize(icnt)/2;
+ }
+
+ // A bitfield is used to mask the pads as usable
+ Short_t mask[kMaxCols]; memset(mask, 0 ,sizeof(Short_t) * kMaxCols);//bool mvals[kMaxPads];
+ for(UChar_t icount = 0; icount < nRows; icount++){
+ cogy[icount] = &cogyvals[icount*kMaxCols];
+ cogz[icount] = &cogzvals[icount*kMaxCols];
+ }
+ // In this array the array position of the best candidates will be stored
+ Int_t cand[AliTRDtrackerV1::kMaxTracksStack];
+ Float_t sigcands[AliTRDtrackerV1::kMaxTracksStack];
+
+ // helper variables
+ Int_t indices[kMaxPads]; memset(indices, -1, sizeof(Int_t)*kMaxPads);
+ Int_t nCandidates = 0;
+ Float_t norm, cogv;
+ // histogram filled -> Select best bins
+ Int_t nPads = nCols * nRows;
+ // take out all the bins which have less than 3 entries (faster sorting)
+ Int_t content[kMaxPads], dictionary[kMaxPads], nCont = 0, padnumber = 0;
+ Int_t *iter = &hvals[0], *citer = &content[0], *diter = &dictionary[0]; // iterators for preselection
+ const Int_t threshold = 2;
+ hvals[nPads] = -1; // termination for iterator
+ do{
+ if(*iter > threshold){
+ *(citer++) = *iter;
+ *(diter++) = padnumber;
+ nCont++;
+ }
+ padnumber++;
+ }while(*(++iter) != -1);
+ TMath::Sort(nCont, content, indices);
+
+ Int_t col, row, lower, lower1, upper, upper1;
+ for(Int_t ib = 0; ib < nCont; ib++){
+ if(nCandidates >= AliTRDtrackerV1::kMaxTracksStack){
+ AliDebug(1, Form("Number of seed candidates %d exceeded maximum allowed per stack %d", nCandidates, AliTRDtrackerV1::kMaxTracksStack));
+ break;
+ }
+ // Positions
+ row = dictionary[indices[ib]]/nCols;
+ col = dictionary[indices[ib]]%nCols;
+ // here will be the threshold condition:
+ if((mask[col] & (1 << row)) != 0) continue; // Pad is masked: continue
+ // if(histogram[row][col] < TMath::Max(threshold, 1)){ // of course at least one cluster is needed
+ // break; // number of clusters below threshold: break;
+ // }
+ // passing: Mark the neighbors
+ lower = TMath::Max(col - 1, 0); upper = TMath::Min(col + 2, nCols);
+ lower1 = TMath::Max(row - 1, 0); upper1 = TMath::Min(row + 2, nCols);
+ for(Int_t ic = lower; ic < upper; ++ic)
+ for(Int_t ir = lower1; ir < upper1; ++ir){
+ if(ic == col && ir == row) continue;
+ mask[ic] |= (1 << ir);
+ }
+ // Storing the position in an array
+ // testing for neigboring
+ cogv = 0;
+ norm = 0;
+ lower = TMath::Max(col - 1, 0);
+ upper = TMath::Min(col + 2, nCols);
+ for(Int_t inb = lower; inb < upper; ++inb){
+ cogv += yLengths[inb] * histogram[row][inb];
+ norm += histogram[row][inb];
+ }
+ cogy[row][col] = cogv / norm;
+ cogv = 0; norm = 0;
+ lower = TMath::Max(row - 1, 0);
+ upper = TMath::Min(row + 2, nRows);
+ for(Int_t inb = lower; inb < upper; ++inb){
+ cogv += zLengths[inb] * histogram[inb][col];
+ norm += histogram[inb][col];
+ }
+ cogz[row][col] = Float_t(cogv) / norm;
+ // passed the filter
+ cand[nCandidates] = row*nCols + col; // store the position of a passig candidate into an Array
+ sigcands[nCandidates] = sigmas[row][col] / histogram[row][col]; // never be a floating point exeption
+ // Analysis output
+ nCandidates++;
+ }
+ if(!nCandidates) return kFALSE;
+
+ Float_t pos[3], sig[2];
+ Short_t signal[7]; memset(&signal[0], 0, 7*sizeof(Short_t));
+
+ new(fakeLayer) AliTRDchamberTimeBin(layer, stack, sector, z0, zl);
+ fakeLayer->SetReconstructor(rec);
+ fakeLayer->SetNRows(nRows);
+ fakeLayer->SetOwner(kFALSE);
+ if(nCandidates){
+ UInt_t fakeIndex = 0;
+ for(Int_t ican = 0; ican < nCandidates; ican++){
+ row = cand[ican] / nCols;
+ col = cand[ican] % nCols;
+ //temporary
+ Int_t n = 0; Double_t x = 0., y = 0., z = 0.;
+ for(int itb=0; itb<AliTRDseedV1::kNtb; itb++){
+ if(!(nClusters = fTB[itb].GetNClusters())) continue;
+ for(Int_t incl = 0; incl < nClusters; incl++){
+ c = fTB[itb].GetCluster(incl);
+ if(c->GetPadRow() != row) continue;
+ if(TMath::Abs(c->GetPadCol() - col) > 2) continue;
+ x += c->GetX();
+ y += c->GetY();
+ z += c->GetZ();
+ n++;
+ }
+ }
+ if(!n) continue;
+ pos[0] = x/n;
+ pos[1] = y/n;
+ pos[2] = z/n;
+ sig[0] = .02;
+ sig[1] = sigcands[ican];
+ fakeLayer->InsertCluster(new AliTRDcluster(fDetector, 0., pos, sig, NULL, 3, signal, col, row, 0, 0, 0., 0), fakeIndex++);
+ }
+ }
+ fakeLayer->BuildIndices();
+ //fakeLayer->Print();
+
+ if(rec->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kTracker) >= 3){
+ //TMatrixD hist(nRows, nCols);
+ //for(Int_t i = 0; i < nRows; i++)
+ // for(Int_t j = 0; j < nCols; j++)
+ // hist(i,j) = histogram[i][j];
+ TTreeSRedirector &cstreamer = *rec->GetDebugStream(AliTRDrecoParam::kTracker);
+ cstreamer << "GetSeedingLayer"
+ << "layer=" << layer
+ << "ymin=" << ymin
+ << "ymax=" << ymax
+ << "zmin=" << zmin
+ << "zmax=" << zmax
+ << "L.=" << fakeLayer
+ //<< "Histogram.=" << &hist
+ << "\n";
+ }
+
+ return kTRUE;
+}
+
+
+//_______________________________________________________
+void AliTRDtrackingChamber::Print(Option_t *opt) const
+{
+ // Print the chamber status
+ if(!GetNClusters()) return;
+ AliInfo(Form("fDetector = %d", fDetector));
+ AliInfo(Form("fX0 = %7.3f", fX0));
+ const AliTRDchamberTimeBin *itb = &fTB[0];
+ for(Int_t jtb=0; jtb<AliTRDseedV1::kNtb; jtb++, itb++) (*itb).Print(opt);
+}
+
+
+//_______________________________________________________
+void AliTRDtrackingChamber::Update()
+{
+// Steer purging of used and shared clusters