if (!fgMinuit) fgMinuit = new TMinuit(8);
fPixArray = new TObjArray(20);
fDebug = 0; //0;
+ fReco = 1;
if (draw) {
fDebug = 1;
+ fReco = 0;
fDraw = new AliMUONClusterDrawAZ(this);
}
- AliWarning("*** Running AZ cluster finder ***");
+ cout << " *** Running AZ cluster finder *** " << endl;
}
//_____________________________________________________________________________
//AZ fResponse = AliMUONClusterInput::Instance()->Response();
Int_t ndigits[2] = {9,9}, nShown[2] = {0};
- for (Int_t i = 0; i < 2; i++) {
- for (Int_t j = 0; j < fgkDim; j++) { fUsed[i][j] = kFALSE; }
+ if (fReco != 2) { // skip initialization for the combined cluster / track
+ fCathBeg = fPadBeg[0] = fPadBeg[1] = 0;
+ for (Int_t i = 0; i < 2; i++) {
+ for (Int_t j = 0; j < fgkDim; j++) { fUsed[i][j] = kFALSE; }
+ }
}
next:
+ if (fReco == 2 && (nShown[0] || nShown[1])) return; // only one precluster for the combined finder
if (ndigits[0] == nShown[0] && ndigits[1] == nShown[1]) return;
Float_t xpad, ypad, zpad, zpad0;
Bool_t first = kTRUE;
if (fDebug) cout << " *** Event # " << nev << " chamber: " << ch << endl;
fnPads[0] = fnPads[1] = 0;
- for (Int_t i = 0; i < fgkDim; i++) { fPadIJ[1][i] = 0; }
+ for (Int_t i = 0; i < fgkDim; i++) fPadIJ[1][i] = 0;
- for (Int_t iii = 0; iii < 2; iii++) {
+ for (Int_t iii = fCathBeg; iii < 2; iii++) {
Int_t cath = TMath::Odd(iii);
ndigits[cath] = AliMUONClusterInput::Instance()->NDigits(cath);
if (!ndigits[0] && !ndigits[1]) return;
Int_t digit;
Bool_t eEOC = kTRUE; // end-of-cluster
- for (digit = 0; digit < ndigits[cath]; digit++) {
+ for (digit = fPadBeg[cath]; digit < ndigits[cath]; digit++) {
mdig = AliMUONClusterInput::Instance()->Digit(cath,digit);
if (first) {
// Find first unused pad
fUsed[cath][digit] = kTRUE;
return;
}
- Int_t isec = fSegmentation[cath]->Sector(fInput->DetElemId(),mdig->PadX(), mdig->PadY());
+ Int_t isec = fSegmentation[cath]->Sector(fInput->DetElemId(), mdig->PadX(), mdig->PadY());
Int_t nPads = fnPads[0] + fnPads[1];
fXyq[0][nPads] = xpad;
fXyq[1][nPads] = ypad;
fXyq[2][i] = -2;
fnPads[cath]--;
}
+ if (fDraw) fDraw->UpdateCluster(npad);
} // if (nFlags > 2)
// Check correlations of cathode charges
if ((over[0] || over[1]) && TMath::Abs(sum[0]-sum[1])/(sum[0]+sum[1])*2 > 1) { // 3 times difference
if (fDebug) cout << " Release " << endl;
// Big difference
- cath = sum[0]>sum[1] ? 0 : 1;
- Int_t imax = 0;
- Double_t cmax=-1;
+ cath = sum[0] > sum[1] ? 0 : 1;
+ Int_t imax = 0, imin = 0;
+ Double_t cmax = -1, cmin = 9999, dxMin = 0, dyMin = 0;
Double_t *dist = new Double_t[npad];
- for (Int_t i=0; i<npad; i++) {
- if (fPadIJ[0][i] != cath) continue;
+ for (Int_t i = 0; i < npad; i++) {
+ if (fPadIJ[0][i] != cath || fXyq[2][i] < 0) continue;
+ if (fXyq[2][i] < cmin) {
+ cmin = fXyq[2][i];
+ imin = i;
+ }
if (fXyq[2][i] < cmax) continue;
cmax = fXyq[2][i];
imax = i;
}
// Arrange pads according to their distance to the max,
// normalized to the pad size
- for (Int_t i=0; i<npad; i++) {
+ for (Int_t i = 0; i < npad; i++) {
dist[i] = 0;
- if (fPadIJ[0][i] != cath) continue;
+ if (fPadIJ[0][i] != cath || fXyq[2][i] < 0) continue;
if (i == imax) continue;
- if (fXyq[2][i] < 0) continue;
- dist[i] = (fXyq[0][i]-fXyq[0][imax])*(fXyq[0][i]-fXyq[0][imax])/
- fXyq[3][imax]/fXyq[3][imax]/4;
- dist[i] += (fXyq[1][i]-fXyq[1][imax])*(fXyq[1][i]-fXyq[1][imax])/
- fXyq[4][imax]/fXyq[4][imax]/4;
- dist[i] = TMath::Sqrt (dist[i]);
+ Double_t dx = (fXyq[0][i] - fXyq[0][imax]) / fXyq[3][imax] / 2;
+ Double_t dy = (fXyq[1][i] - fXyq[1][imax]) / fXyq[4][imax] / 2;
+ dist[i] = TMath::Sqrt (dx * dx + dy * dy);
+ if (i == imin) {
+ cmin = dist[i]; // distance to the pad with minimum charge
+ dxMin = dx;
+ dyMin = dy;
+ }
}
TMath::Sort(npad, dist, flags, kFALSE); // in increasing order
Int_t indx;
Double_t xmax = -1;
- for (Int_t i=0; i<npad; i++) {
+ for (Int_t i = 0; i < npad; i++) {
indx = flags[i];
- if (fPadIJ[0][indx] != cath) continue;
- if (fXyq[2][indx] < 0) continue;
- if (fXyq[2][indx] <= cmax || TMath::Abs(dist[indx]-xmax)<1.e-3) {
+ if (fPadIJ[0][indx] != cath || fXyq[2][indx] < 0) continue;
+ if (dist[indx] > cmin) {
+ // Farther than the minimum pad
+ Double_t dx = (fXyq[0][indx] - fXyq[0][imax]) / fXyq[3][imax] / 2;
+ Double_t dy = (fXyq[1][indx] - fXyq[1][imax]) / fXyq[4][imax] / 2;
+ dx *= dxMin;
+ dy *= dyMin;
+ if (dx >= 0 && dy >= 0) continue;
+ if (TMath::Abs(dx) > TMath::Abs(dy) && dx >= 0) continue;
+ if (TMath::Abs(dy) > TMath::Abs(dx) && dy >= 0) continue;
+ }
+ if (fXyq[2][indx] <= cmax || TMath::Abs(dist[indx]-xmax) < 1.e-3) {
// Release pads
- if (TMath::Abs(dist[indx]-xmax)<1.e-3)
+ if (TMath::Abs(dist[indx]-xmax) < 1.e-3)
cmax = TMath::Max((Double_t)(fXyq[2][indx]),cmax);
else cmax = fXyq[2][indx];
xmax = dist[indx];
fUsed[cath][digit] = kFALSE;
fXyq[2][indx] = -2;
fnPads[cath]--;
- // xmax = dist[i]; // Bug?
- } else {
- // Check pad overlaps once more
- for (Int_t i=0; i<npad; i++) flags[i] = 0;
- for (Int_t i=0; i<npad; i++) {
- if (fXyq[2][i] < 0) continue;
- if (fPadIJ[0][i] != i1) continue;
- xy1[0] = fXyq[0][i] - fXyq[3][i];
- xy1[1] = fXyq[0][i] + fXyq[3][i];
- xy1[2] = fXyq[1][i] - fXyq[4][i];
- xy1[3] = fXyq[1][i] + fXyq[4][i];
- for (Int_t j=0; j<npad; j++) {
- if (fXyq[2][j] < 0) continue;
- if (fPadIJ[0][j] != i2) continue;
- if (!Overlap(xy1, j, xy12, 0)) continue;
- flags[i] = flags[j] = 1; // mark overlapped pads
- } // for (Int_t j=0;
- } // for (Int_t i=0;
- nFlags=0;
- for (Int_t i=0; i<npad; i++) {
- if (fXyq[2][i] < 0 || flags[i]) continue;
- nFlags ++;
- }
- if (nFlags == fnPads[0] + fnPads[1]) {
- // No overlap
- for (Int_t i=0; i<npad; i++) {
- if (fXyq[2][i] < 0 || fPadIJ[0][i] != cath) continue;
- fXyq[2][i] = -2;
- fnPads[cath]--;
- }
- }
- break;
+ }
+ } // for (Int_t i = 0; i < npad;
+
+ // Check pad overlaps once more
+ for (Int_t j = 0; j < npad; j++) flags[j] = 0;
+ for (Int_t k = 0; k < npad; k++) {
+ if (fXyq[2][k] < 0 || fPadIJ[0][k] != i1) continue;
+ xy1[0] = fXyq[0][k] - fXyq[3][k];
+ xy1[1] = fXyq[0][k] + fXyq[3][k];
+ xy1[2] = fXyq[1][k] - fXyq[4][k];
+ xy1[3] = fXyq[1][k] + fXyq[4][k];
+ for (Int_t j = 0; j < npad; j++) {
+ if (fXyq[2][j] < 0) continue;
+ if (fPadIJ[0][j] != i2) continue;
+ if (!Overlap(xy1, j, xy12, 0)) continue;
+ flags[k] = flags[j] = 1; // mark overlapped pads
+ } // for (Int_t j = 0;
+ } // for (Int_t k = 0;
+ nFlags = 0;
+ for (Int_t j = 0; j < npad; j++) {
+ if (fXyq[2][j] < 0 || flags[j]) continue;
+ nFlags ++;
+ }
+ if (nFlags == fnPads[0] + fnPads[1]) {
+ // No overlap
+ for (Int_t j = 0; j < npad; j++) {
+ if (fXyq[2][j] < 0 || fPadIJ[0][j] != cath) continue;
+ fXyq[2][j] = -2;
+ fnPads[cath]--;
}
- }
+ }
delete [] dist; dist = 0;
+ if (fDraw) fDraw->UpdateCluster(npad);
} // TMath::Abs(sum[0]-sum[1])...
} // if (fnPads[0] && fnPads[1])
delete [] flags; flags = 0;
{
// Check if some pixels have large size (adjust if necessary)
- Int_t nx, ny;
- Int_t nPix = fPixArray->GetEntriesFast();
- AliMUONPixel *pixPtr, *pixPtr1, pix;
+ Int_t n1[2], n2[2], iOK = 1, nPix = fPixArray->GetEntriesFast();
+ AliMUONPixel *pixPtr, pix;
+ Double_t xy0[2] = {9999, 9999}, wxy[2], dist[2];
// Check if large pixel size
- for (Int_t i=0; i<nPix; i++) {
+ for (Int_t i = 0; i < nPix; i++) {
pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(i);
if (pixPtr->Charge() < 1) continue; // discarded pixel
- if (pixPtr->Size(0)-wxmin > 1.e-4 || pixPtr->Size(1)-wymin > 1.e-4) {
- if (fDebug) cout << " Different " << pixPtr->Size(0) << " " << wxmin << " " << pixPtr->Size(1) << " " << wymin << endl;
- pix = *pixPtr;
- nx = TMath::Nint (pix.Size(0)/wxmin);
- ny = TMath::Nint (pix.Size(1)/wymin);
- pix.Shift(0, -pix.Size(0)-wxmin);
- pix.Shift(1, -pix.Size(1)-wymin);
- pix.SetSize(0, wxmin);
- pix.SetSize(1, wymin);
- for (Int_t ii=0; ii<nx; ii++) {
- pix.Shift(0, wxmin*2);
- for (Int_t jj=0; jj<ny; jj++) {
- pix.Shift(1, wymin*2);
- pixPtr1 = new AliMUONPixel(pix);
- fPixArray->Add((TObject*)pixPtr1);
- }
+ if (pixPtr->Size(0) - wxmin < 1.e-4) {
+ if (xy0[0] > 9998) xy0[0] = pixPtr->Coord(0); // position of a "normal" pixel
+ if (pixPtr->Size(1) - wymin < 1.e-4) {
+ if (xy0[1] > 9998) xy0[1] = pixPtr->Coord(1); // position of a "normal" pixel
+ continue;
+ } else iOK = 0; // large pixel
+ } else {
+ iOK = 0; // large pixel
+ if (xy0[1] > 9998 && pixPtr->Size(1) - wymin < 1.e-4) xy0[1] = pixPtr->Coord(1); // "normal" pixel
+ }
+ if (xy0[0] < 9998 && xy0[1] < 9998) break;
+ }
+ if (iOK) return;
+
+ wxy[0] = wxmin;
+ wxy[1] = wymin;
+ //cout << xy0[0] << " " << xy0[1] << endl;
+ for (Int_t i = 0; i < nPix; i++) {
+ pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(i);
+ if (pixPtr->Charge() < 1) continue; // discarded pixel
+ n1[0] = n1[1] = 999;
+ n2[0] = n2[1] = 1;
+ for (Int_t j = 0; j < 2; j++) {
+ if (pixPtr->Size(j) - wxy[j] < 1.e-4) continue;
+ dist[j] = (pixPtr->Coord(j) - xy0[j]) / wxy[j] / 2; // normalized distance to "normal" pixel
+ n2[j] = TMath::Nint (pixPtr->Size(j) / wxy[j]);
+ n1[j] = n2[j] == 1 ? TMath::Nint(dist[j]) : (Int_t)dist[j];
+ }
+ if (n1[0] > 998 && n1[1] > 998) continue;
+ if (fDebug) cout << " Different " << pixPtr->Size(0) << " " << wxy[0] << " "
+ << pixPtr->Size(1) << " " << wxy[1] <<endl;
+
+ if (n2[0] > 2 || n2[1] > 2) { cout << n2[0] << " " << n2[1] << endl; AliFatal("Too large pixel."); }
+ //cout << n1[0] << " " << n2[0] << " " << n1[1] << " " << n2[1] << endl;
+ pix = *pixPtr;
+ pix.SetSize(0, wxy[0]); pix.SetSize(1, wxy[1]);
+ //pixPtr->Print();
+ for (Int_t ii = 0; ii < n2[0]; ii++) {
+ if (n1[0] < 999) pix.SetCoord(0, xy0[0] + (n1[0] + TMath::Sign(1.,dist[0]) * ii) * 2 * wxy[0]);
+ for (Int_t jj = 0; jj < n2[1]; jj++) {
+ if (n1[1] < 999) pix.SetCoord(1, xy0[1] + (n1[1] + TMath::Sign(1.,dist[1]) * jj) * 2 * wxy[1]);
+ fPixArray->Add(new AliMUONPixel(pix));
+ //pix.Print();
}
- pixPtr->SetCharge(0);
}
- } // for (Int_t i=0; i<nPix;
- return;
+ pixPtr->SetCharge(0);
+ } // for (Int_t i = 0; i < nPix;
}
//_____________________________________________________________________________
Double_t xylim[4] = {999, 999, 999, 999};
for (Int_t ipix=0; ipix<nPix; ipix++) {
pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(ipix);
+ //cout << ipix+1; pixPtr->Print();
for (Int_t i=0; i<4; i++)
xylim[i] = TMath::Min (xylim[i], (i%2 ? -1 : 1)*pixPtr->Coord(i/2));
- //cout << ipix+1; pixPtr->Print();
}
for (Int_t i=0; i<4; i++) {
xylim[i] -= pixPtr->Size(i/2); if (fDebug) cout << (i%2 ? -1 : 1)*xylim[i] << " "; }
Int_t indx;
fnPads[1] -= nVirtual;
- Double_t coef = 0;
- //for (Int_t j=0; j<nfit; j++) {
- for (Int_t j=nfit-1; j>=0; j--) {
- indx = j<2 ? j*2 : j*2+1;
- if (nfit == 1) coef = 1;
- else coef = j==nfit-1 ? parOk[indx+2] : 1-coef;
- coef = TMath::Max (coef, 0.);
- if (nfit == 3 && j < 2) coef = j==1 ? coef*parOk[indx+2] : coef - parOk[7];
- coef = TMath::Max (coef, 0.);
- AddRawCluster (parOk[indx], parOk[indx+1], coef*fQtot, errOk[indx], nfit0+10*nfit, tracks,
- //sigCand[maxSeed[j]][0], sigCand[maxSeed[j]][1]);
- //sigCand[0][0], sigCand[0][1], dist[j]);
- sigCand[0][0], sigCand[0][1], dist[TMath::LocMin(nfit,dist)]);
- }
- if (fDraw) fDraw->FillMuon(nfit, parOk, errOk);
+ if (!fDraw) {
+ Double_t coef = 0;
+ //for (Int_t j=0; j<nfit; j++) {
+ for (Int_t j=nfit-1; j>=0; j--) {
+ indx = j<2 ? j*2 : j*2+1;
+ if (nfit == 1) coef = 1;
+ else coef = j==nfit-1 ? parOk[indx+2] : 1-coef;
+ coef = TMath::Max (coef, 0.);
+ if (nfit == 3 && j < 2) coef = j==1 ? coef*parOk[indx+2] : coef - parOk[7];
+ coef = TMath::Max (coef, 0.);
+ AddRawCluster (parOk[indx], parOk[indx+1], coef*fQtot, errOk[indx], nfit0+10*nfit, tracks,
+ //sigCand[maxSeed[j]][0], sigCand[maxSeed[j]][1]);
+ //sigCand[0][0], sigCand[0][1], dist[j]);
+ sigCand[0][0], sigCand[0][1], dist[TMath::LocMin(nfit,dist)]);
+ }
+ } else fDraw->FillMuon(nfit, parOk, errOk);
return nfit;
}
//if (iPad && TMath::Abs(fXyq[1][npads]-fXyq[1][iAddY]) < fXyq[4][iAddY]) continue;
//if (TMath::Abs(fXyq[0][npads]) < 1 && TMath::Abs(fXyq[1][npads]) < 1) continue; // strange case (something with pad mapping)
if (maxpad[0][0] < 0 || mirror && maxpad[1][0] >= 0) {
- if (!iPad) fXyq[2][npads] = TMath::Min (sigmax[1]/20, 5.);
- else fXyq[2][npads] = TMath::Min (aamax[1]/20, 5.);
+ //if (!iPad) fXyq[2][npads] = TMath::Min (sigmax[1]/20, 5.);
+ //else fXyq[2][npads] = TMath::Min (aamax[1]/20, 5.);
+ if (!iPad) fXyq[2][npads] = TMath::Min (sigmax[1]/10, (double)fResponse->ZeroSuppression());
+ else fXyq[2][npads] = TMath::Min (aamax[1]/10, (double)fResponse->ZeroSuppression());
}
else {
- if (!iPad) fXyq[2][npads] = TMath::Min (sigmax[0]/20, 5.);
- else fXyq[2][npads] = TMath::Min (aamax[0]/20, 5.);
+ //if (!iPad) fXyq[2][npads] = TMath::Min (sigmax[0]/20, 5.);
+ //else fXyq[2][npads] = TMath::Min (aamax[0]/20, 5.);
+ if (!iPad) fXyq[2][npads] = TMath::Min (sigmax[0]/10, (double)fResponse->ZeroSuppression());
+ else fXyq[2][npads] = TMath::Min (aamax[0]/10, (double)fResponse->ZeroSuppression());
}
fXyq[2][npads] = TMath::Max (fXyq[2][npads], (float)1);
//isec = fSegmentation[cath]->Sector(fInput->DetElemId(),ix, iy);
if (fnPads[1]) { delete [] xPad1; delete [] yPad1; delete [] nPad1; }
//nInY = TMath::Max (npady[0], npady[1]);
//nInX = TMath::Max (npadx[0], npadx[1]);
- nInY = wMinY[0] < wMinY[1] ? npady[0] : npady[1];
- nInX = wMinX[0] < wMinX[1] ? npadx[0] : npadx[1];
+ if (TMath::Abs (wMinY[0] - wMinY[1]) < 1.e-3) nInY = TMath::Max (npady[0], npady[1]);
+ else nInY = wMinY[0] < wMinY[1] ? npady[0] : npady[1];
+ if (TMath::Abs (wMinX[0] - wMinX[1]) < 1.e-3) nInX = TMath::Max (npadx[0], npadx[1]);
+ else nInX = wMinX[0] < wMinX[1] ? npadx[0] : npadx[1];
}
//_____________________________________________________________________________
Int_t GetIJ(Int_t indx, Int_t iPad) const { return fPadIJ[indx][iPad]; }
Float_t GetXyq(Int_t indx, Int_t iPad) const { return fXyq[indx][iPad]; }
Float_t GetZpad() { return fZpad; }
+ Bool_t GetUsed(Int_t cath, Int_t dig) const { return fUsed[cath][dig]; }
+ void SetUsed(Int_t cath, Int_t dig) { fUsed[cath][dig] = kTRUE; } // mark used digits
+ void SetUnused(Int_t cath, Int_t dig) { fUsed[cath][dig] = kFALSE; } // unmark digits
+ void SetReco(Int_t iReco) { fReco = iReco; } // set reco flag
+ void SetStart(Int_t iCath, Int_t iPad) { fCathBeg = iCath; fPadBeg[0] = fPadBeg[1] = 0; fPadBeg[fCathBeg] = iPad; } // start
protected:
AliMUONClusterFinderAZ(const AliMUONClusterFinderAZ& rhs);
static AliMUONClusterFinderAZ* fgClusterFinder; // the ClusterFinderAZ instance
- Int_t fnPads[2]; // ! number of pads in the cluster on 2 cathodes
- Float_t fXyq[7][fgkDim]; // ! pad information
- Int_t fPadIJ[2][fgkDim]; // ! pad information
+ Int_t fnPads[2]; // ! number of pads in the cluster on 2 cathodes
+ Float_t fXyq[7][fgkDim]; // ! pad information
+ Int_t fPadIJ[2][fgkDim]; // ! pad information
AliMUONGeometrySegmentation *fSegmentation[2]; // ! new segmentation
- AliMUONResponse *fResponse;// ! response
- Float_t fZpad; // ! z-coordinate of the hit
- Int_t fNpar; // ! number of fit parameters
- Double_t fQtot; // ! total cluster charge
- Int_t fReco; // ! =1 if run reco with writing of reconstructed clusters
+ AliMUONResponse *fResponse; // ! response
+ Float_t fZpad; // ! z-coordinate of the hit
+ Int_t fNpar; // ! number of fit parameters
+ Double_t fQtot; // ! total cluster charge
+ Int_t fReco; // ! !=0 if run reco with writing of reconstructed clusters
+ Int_t fCathBeg; // ! starting cathode (for combined cluster / track reco)
+ Int_t fPadBeg[2]; // ! starting pads (for combined cluster / track reco)
static TMinuit* fgMinuit; // ! Fitter
Bool_t fUsed[2][fgkDim]; // ! flags for used pads
+ //TH2F* fHist[4]; // ! histograms
AliMUONClusterDrawAZ *fDraw; // ! drawing object
+ //Int_t fnMu; // ! number of muons passing thru the selected area
+ //Double_t fxyMu[2][7]; // ! muon information
TObjArray* fPixArray; // ! collection of pixels
Int_t fnCoupled; // ! number of coupled clusters in precluster
Int_t fDebug; // ! debug level
return;
}
//____________________________________________________________________
-void AliMUONClusterReconstructor::Digits2Clusters()
+void AliMUONClusterReconstructor::Digits2Clusters(Int_t chBeg)
{
TClonesArray *dig1, *dig2, *digAll;
Int_t n2;
Int_t n1;
- for (Int_t ich = 0; ich < AliMUONConstants::NTrackingCh(); ich++) {
+ fMUONData->ResetDigits(); //AZ
+ fMUONData->GetDigits(); //AZ
+
+ for (Int_t ich = chBeg; ich < AliMUONConstants::NTrackingCh(); ich++) {
id.Reset();
n1 = 0;
n2 = 0;
//cathode 0 & 1
- fMUONData->ResetDigits();
- fMUONData->GetDigits();
+ //fMUONData->ResetDigits();
+ //fMUONData->GetDigits();
+
muonDigits = fMUONData->Digits(ich);
ndig = muonDigits->GetEntriesFast();
TClonesArray &lDigit = *digAll;
void AliMUONClusterReconstructor::Trigger2Trigger()
{
// copy trigger from TreeD to TreeR
+
+ fMUONData->SetTreeAddress("GLT");
fMUONData->GetTriggerD();
}
/////////////////////////////////////
#include <TObject.h>
-#include "AliMUONClusterFinderVS.h" //AZ
+#include "AliMUONClusterFinderVS.h"
class AliLoader;
class AliMUON;
class AliMUONRawCluster;
-//AZ class AliMUONClusterFinderVS;
class AliMUONData;
class AliRawReader;
// Cluster Finding & Trigger
- virtual void Digits2Clusters();
+ virtual void Digits2Clusters(Int_t chBeg = 0);
virtual void Trigger2Trigger() ;
// pointer to data container
AliMUONData* GetMUONData() {return fMUONData;}
// Reco Model
AliMUONClusterFinderVS* GetRecoModel() {return fRecModel;}
- // AliMUONClusterFinderAZ* GetRecoModel() {return fRecModel;}
//AZ void SetRecoModel(AliMUONClusterFinderVS* rec) {fRecModel = rec;}
void SetRecoModel(AliMUONClusterFinderVS* rec) {if (fRecModel) delete fRecModel; fRecModel = rec;} //AZ
- // void SetRecoModel(AliMUONClusterFinderAZ* rec) {fRecModel = rec;}
protected:
AliMUONData* fMUONData; //! Data container for MUON subsystem
AliMUONClusterFinderVS* fRecModel; //! cluster recontruction model
- //AliMUONClusterFinderAZ* fRecModel; //! cluster recontruction model
// alice loader
AliLoader* fLoader;
}
}
- if ( TreeR() && fRawClusters && cRC) {
+ if ( TreeR() && fRawClusters && cRC && !strstr(cRC,"RCC")) {
for (int i=0; i<AliMUONConstants::NTrackingCh(); i++) {
sprintf(branchname,"%sRawClusters%d",GetName(),i+1);
if (fRawClusters) {
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+/* $Id$ */
+
+#include <TObjArray.h>
+#include <TClonesArray.h>
+#include "AliMUONDetElement.h"
+#include "AliMUON.h"
+#include "AliMUONSegmentation.h"
+#include "AliMUONDigit.h"
+#include "AliMUONHitMapA1.h"
+#include "AliMUONData.h"
+#include "AliMUONRawCluster.h"
+#include "AliMUONHitForRec.h"
+#include "AliMUONClusterInput.h"
+#include "AliMUONClusterFinderAZ.h"
+#include "AliRun.h"
+#include "AliLog.h"
+
+ClassImp(AliMUONDetElement) // Class implementation in ROOT context
+ FILE *lun = 0x0; //fopen("hitmap.dat","w");
+
+//_________________________________________________________________________
+AliMUONDetElement::AliMUONDetElement()
+ : TObject()
+{
+// Default constructor
+ for (Int_t i = 0; i < 2; i++) {
+ fHitMap[i] = NULL;
+ fDigits[i] = NULL;
+ fSeg[i] = NULL;
+ }
+ fRawClus = fHitsForRec = NULL;
+ fRecModel = NULL;
+}
+
+//_________________________________________________________________________
+AliMUONDetElement::AliMUONDetElement(Int_t idDE, AliMUONDigit *dig, AliMUONClusterFinderAZ *recModel)
+ : TObject()
+{
+ // Constructor
+ fidDE = idDE;
+ fChamber = fidDE / 100 - 1;
+ fDigits[0] = new TObjArray(10);
+ fDigits[1] = new TObjArray(10);
+ fRawClus = new TObjArray(10);
+ fHitsForRec = new TClonesArray("AliMUONHitForRec",10);
+ fNHitsForRec = 0;
+ fRecModel = recModel;
+ AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON");
+ AliMUONSegmentation *pSegmentation = pMUON->GetSegmentation();
+ fSeg[0] = pSegmentation->GetModuleSegmentation(fChamber, 0);
+ fSeg[1] = pSegmentation->GetModuleSegmentation(fChamber, 1);
+ Float_t x, y, z;
+ fSeg[dig->Cathode()]->GetPadC(fidDE, dig->PadX(), dig->PadY(), x, y, z);
+ fZ = z;
+ AddDigit(dig);
+}
+
+//_________________________________________________________________________
+AliMUONDetElement::~AliMUONDetElement()
+{
+ // Destructor
+ for (Int_t i = 0; i < 2; i++) {
+ delete fHitMap[i]; fHitMap[i] = NULL;
+ delete fDigits[i]; fDigits[i] = NULL;
+ }
+ if (fRawClus) { fRawClus->Delete(); delete fRawClus; fRawClus = 0; }
+ //if (fRawClus) { delete fRawClus; fRawClus = 0; }
+ delete fHitsForRec; fHitsForRec = 0;
+}
+
+//_________________________________________________________________________
+AliMUONDetElement::AliMUONDetElement (const AliMUONDetElement& rhs)
+ : TObject(rhs)
+{
+// Copy constructor
+
+ AliFatal("Not implemented.");
+}
+
+//________________________________________________________________________
+AliMUONDetElement & AliMUONDetElement::operator = (const AliMUONDetElement& rhs)
+{
+// Assignement operator
+
+ if (this == &rhs) return *this;
+ AliFatal( "Not implemented.");
+ return *this;
+}
+
+//_________________________________________________________________________
+Int_t AliMUONDetElement::Compare(const TObject* detElem) const
+{
+ // "Compare" function to sort in Z (towards interaction point)
+ // Returns -1 (0, +1) if charge of current pixel
+ // is greater than (equal to, less than) charge of pixel
+ if (fZ > ((AliMUONDetElement*)detElem)->Z()) return(+1);
+ else if (fZ == ((AliMUONDetElement*)detElem)->Z()) return( 0);
+ else return(-1);
+}
+
+//_________________________________________________________________________
+void AliMUONDetElement::Fill(AliMUONData */*data*/)
+{
+ // Fill hit maps
+ fLeft[0] = fDigits[0]->GetEntriesFast();
+ fLeft[1] = fDigits[1]->GetEntriesFast();
+
+ fHitMap[0] = new AliMUONHitMapA1(fidDE, fSeg[0], fDigits[0]);
+ fHitMap[1] = new AliMUONHitMapA1(fidDE, fSeg[1], fDigits[1]);
+ fHitMap[0]->FillHits();
+ fHitMap[1]->FillHits();
+
+ // The part below is just for debugging (fill rec. points already found)
+ /*
+ fLeft[0] = fLeft[1] = 0;
+ TClonesArray *rawClus = data->RawClusters(fChamber);
+ cout << rawClus << " " << rawClus->GetEntriesFast() << endl;
+ for (Int_t i = 0; i < rawClus->GetEntriesFast(); i++) {
+ AliMUONRawCluster *recP = (AliMUONRawCluster*) rawClus->UncheckedAt(i);
+ cout << fChamber << " " << recP->GetZ(0) << " " << recP->GetZ(1) << " " << fZ << endl;
+ if (TMath::Abs(recP->GetZ(0)-fZ) > 0.5) continue;
+ if (!Inside(recP->GetX(0), recP->GetY(0), recP->GetZ(0))) continue;
+ AddHitForRec(recP); // add hit for rec.
+ rawClus->RemoveAt(i); // remove
+ }
+ cout << fHitsForRec->GetEntriesFast() << endl;
+ rawClus->Compress();
+ */
+}
+
+//_________________________________________________________________________
+void AliMUONDetElement::AddDigit(AliMUONDigit *dig)
+{
+ // Add digit
+
+ fDigits[dig->Cathode()]->Add(dig);
+}
+
+//_________________________________________________________________________
+Bool_t AliMUONDetElement::Inside(Double_t x, Double_t y, Double_t z) const
+{
+ // Check if point is inside detection element
+
+ Int_t ix, iy;
+ for (Int_t i = 0; i < 2; i++) {
+ if (!fSeg[i]) continue;
+ fSeg[i]->GetPadI(fidDE, x, y, z, ix, iy);
+ //cout << x << " " << y << " " << z << " " << fChamber << " " << ix << " " << iy << " " << fSeg[i]->Npx(fidDE) << " " << fSeg[i]->Npy(fidDE) /*<< " " << fSeg[i]->GetPadI(fidDE, x, y, z, ix, iy)*/ << endl;
+ if (ix > 0 && iy > 0 && ix <= fSeg[i]->Npx(fidDE) && iy <= fSeg[i]->Npy(fidDE)) return kTRUE;
+ }
+ // Check for edge effect (extrapolated track "right outside" det. elem. boundaries (+- 1cm in X and Y)
+ for (Int_t i = 0; i < 2; i++) {
+ if (!fSeg[i]) continue;
+ for (Int_t idx = -1; idx < 2; idx++) {
+ Double_t x1 = x + 1. * idx;
+ for (Int_t idy = -1; idy < 2; idy++) {
+ if (idx == 0 && idy == 0) continue;
+ Double_t y1 = y + 1. * idy;
+ fSeg[i]->GetPadI(fidDE, x1, y1, z, ix, iy);
+ //cout << x1 << " " << y1 << " " << z << " " << fChamber << " " << ix << " " << iy << " " << fSeg[i]->Npx(fidDE) << " " << fSeg[i]->Npy(fidDE) /*<< " " << fSeg[i]->GetPadI(fidDE, x, y, z, ix, iy)*/ << endl;
+ if (ix > 0 && iy > 0 && ix <= fSeg[i]->Npx(fidDE) && iy <= fSeg[i]->Npy(fidDE)) return kTRUE;
+ }
+ }
+ }
+ return kFALSE;
+}
+
+//_________________________________________________________________________
+void AliMUONDetElement::ClusterReco(Double_t xTrack, Double_t yTrack)
+{
+ // Run cluster reconstruction around point (x,y)
+
+ if (fLeft[0] == 0 && fLeft[1] == 0) return; // all digits have been used
+ Float_t dx, dy;
+ dx = dy = 5; // 5 cm for now
+ AliMUONClusterInput::Instance()->SetDigits(fChamber, fidDE,
+ (TClonesArray*)fDigits[0], (TClonesArray*)fDigits[1]);
+
+ // Mark used pads
+ for (Int_t cath = 0; cath < 2; cath++) {
+ if (fDigits[cath]->GetEntriesFast() == 0) continue; // empty cathode
+ for (Int_t i = 0; i < fDigits[cath]->GetEntriesFast(); i++) {
+ if (fLeft[cath] == 0) { fRecModel->SetUsed(cath,i); continue; }
+ AliMUONDigit *dig = (AliMUONDigit*) fDigits[cath]->UncheckedAt(i);
+ //cout << i << " " << dig->PadX() << " " << dig->PadY() << " " << fHitMap[cath]->TestHit(dig->PadX(), dig->PadY()) << endl;
+ if (fHitMap[cath]->TestHit(dig->PadX(), dig->PadY()) == kUsed) fRecModel->SetUsed(cath,i);
+ else fRecModel->SetUnused(cath,i);
+ }
+ }
+
+ fRecModel->ResetRawClusters();
+
+ for (Int_t cath = 0; cath < 2; cath++) {
+ if (fDigits[cath]->GetEntriesFast() == 0) continue; // empty cathode
+ // Loop over pads
+ for (fSeg[cath]->FirstPad(fidDE, xTrack, yTrack, fZ, dx, dy);
+ fSeg[cath]->MorePads(fidDE);
+ fSeg[cath]->NextPad(fidDE)) {
+ if (fLeft[cath] == 0) break;
+ //cout << cath << " " << fSeg[cath]->Ix() << " " << fSeg[cath]->Iy() << " " << fSeg[cath]->DetElemId() << " " << fHitMap[cath]->TestHit(fSeg[cath]->Ix(), fSeg[cath]->Iy()) << endl;
+ if (fHitMap[cath]->TestHit(fSeg[cath]->Ix(), fSeg[cath]->Iy()) == kEmpty ||
+ fHitMap[cath]->TestHit(fSeg[cath]->Ix(), fSeg[cath]->Iy()) == kUsed) continue;
+
+ // Set starting pad
+ for (Int_t j = 0; j < fDigits[cath]->GetEntriesFast(); j++) {
+ AliMUONDigit *dig = (AliMUONDigit*) fDigits[cath]->UncheckedAt(j);
+ if (dig->PadX() != fSeg[cath]->Ix() || dig->PadY() != fSeg[cath]->Iy()) continue;
+ //cout << fidDE << " " << j << " " << fSeg[cath]->Ix() << " " << fSeg[cath]->Iy() << endl;
+ fRecModel->SetStart(cath, j);
+ break;
+ }
+
+ fRecModel->FindRawClusters();
+ Int_t nClusEnd = fRecModel->GetRawClusters()->GetEntriesFast();
+ //cout << " ***nclus: " << nClusEnd << endl;
+ for (Int_t i = 0; i < nClusEnd; i++) {
+ AliMUONRawCluster *clus = (AliMUONRawCluster*) fRecModel->GetRawClusters()->UncheckedAt(i);
+ AddHitForRec(clus); // add hit for rec.
+ //cout << clus->GetX(0) << " " << clus->GetY(0) << endl;
+ }
+ // Mark used pads
+ for (Int_t cath1 = 0; cath1 < 2; cath1++) {
+ for (Int_t j = 0; j < fDigits[cath1]->GetEntriesFast(); j++) {
+ if (fLeft[cath1] == 0) break;
+ AliMUONDigit *dig = (AliMUONDigit*) fDigits[cath1]->UncheckedAt(j);
+ Float_t x, y, z;
+ fSeg[cath1]->GetPadC(fidDE,dig->PadX(),dig->PadY(),x,y,z);
+ //cout << "clus " << cath1 << " " << fLeft[cath1] << " " << dig->PadX() << " " << dig->PadY() << " " << x << " " << y << " " << z << " " << fRecModel->GetUsed(cath1,j) << endl;
+ if (!fRecModel->GetUsed(cath1,j)) continue;
+ if (fHitMap[cath1]->TestHit(dig->PadX(), dig->PadY()) == kUsed) continue;
+ fHitMap[cath1]->FlagHit(dig->PadX(), dig->PadY());
+ if (lun) fprintf(lun," %d %d %d %d \n", cath1, fidDE, dig->PadX(), dig->PadY());
+ fLeft[cath1]--;
+ }
+ }
+ } // for (fSeg[cath]->FirstPad(...
+ } // for (Int_t cath = 0;
+}
+
+//_________________________________________________________________________
+void AliMUONDetElement::AddHitForRec(AliMUONRawCluster *clus)
+{
+ // Make HitForRec from raw cluster (rec. point)
+
+ fRawClus->Add(new AliMUONRawCluster(*clus));
+ AliMUONHitForRec *hitForRec =
+ new ((*fHitsForRec)[fNHitsForRec++]) AliMUONHitForRec(clus);
+
+ // more information into HitForRec
+ // resolution: info should be already in raw cluster and taken from it ????
+ hitForRec->SetBendingReso2(-1); //fBendingResolution * fBendingResolution);
+ hitForRec->SetNonBendingReso2(-1); //fNonBendingResolution * fNonBendingResolution);
+ // original raw cluster
+ hitForRec->SetChamberNumber(fChamber);
+ hitForRec->SetZ(clus->GetZ(0));
+ //hitForRec->SetHitNumber(-(fIndex+1)*100000-fNHitsForRec+1);
+ hitForRec->SetHitNumber(-(fIndex+1)*100000-fRawClus->GetEntriesFast()+1);
+ //delete clus; // for now
+}
+
+/*
+//_________________________________________________________________________
+Int_t AliMUONDetElement::GetMapElem(AliMUONDigit *digit)
+{
+ Int_t cath = digit->Cathode();
+ return 0;
+
+}
+
+//_________________________________________________________________________
+void AliMUONDetElement::SetMapElem(const AliMUONDigit *digit, Int_t flag)
+{
+ Int_t cath = digit->Cathode();
+}
+*/
--- /dev/null
+#ifndef ALIMUONDETELEMENT_H
+#define ALIMUONDETELEMENT_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+/* $Id$ */
+
+/// \ingroup rec
+/// \class AliMUONDetElement
+/// \brief Detection element object containing information for combined cluster / track finder in MUON arm
+
+#include <TObject.h>
+class TObjArray;
+class TClonesArray;
+class AliMUONDigit;
+class AliMUONHitMapA1;
+class AliMUONGeometrySegmentation;
+class AliMUONData;
+class AliMUONRawCluster;
+class AliMUONClusterFinderAZ;
+
+class AliMUONDetElement : public TObject
+{
+ public:
+
+ AliMUONDetElement();
+ AliMUONDetElement(Int_t idDE, AliMUONDigit *dig, AliMUONClusterFinderAZ *recModel); // constructor
+ AliMUONDetElement(const AliMUONDetElement & rhs); // copy constructor
+ AliMUONDetElement& operator = (const AliMUONDetElement& rhs); // assignment operator
+ virtual ~AliMUONDetElement(); // Destructor
+
+ Int_t IdDE(void) const { return fidDE; } // det. elem. ID
+ Int_t Chamber(void) const { return fChamber; } // chamber No
+ Double_t Z(void) const { return fZ; } // Z-coordinate
+ Int_t Left(Int_t cath) const { return fLeft[cath]; } // number of unused digits
+ //Int_t GetMapElem(AliMUONDigit *digit) const; // get map element
+ TObjArray *Digits(Int_t cath) const { return fDigits[cath]; } // array of digits
+ TObjArray *RawClusters() const { return fRawClus; } // array of raw clusters
+ TClonesArray *HitsForRec() const { return fHitsForRec; } // hits for rec.
+ Int_t NHitsForRec() const { return fNHitsForRec; } // No. of hits for rec.
+ Bool_t Inside(Double_t x, Double_t y, Double_t z) const; // check if point inside DE
+
+ void SetID(Int_t idDE) { fidDE = idDE; } // set det. elem. ID
+ void SetIndex(Int_t index) { fIndex = index; } // set position index
+ void SetZ(Double_t z) { fZ = z; } // set Z-coord.
+ void SetLeft(Int_t cath, Int_t left) { fLeft[cath] = left; } // set No. of digits
+ //void SetMapElem(const AliMUONDigit *digit, Int_t flag); // set map element
+ void AddDigit(AliMUONDigit *dig); // add digit
+ void Fill(AliMUONData *data); // fill hit maps
+ void ClusterReco(Double_t x, Double_t y); // run cluster reco around (x,y)
+ void AddHitForRec(AliMUONRawCluster *clus); // make HitForRec
+ // What is necessary for sorting TObjArray's
+ Bool_t IsSortable() const { return kTRUE; }
+ Int_t Compare(const TObject* detElem) const; // "Compare" function for sorting
+
+ protected:
+
+ private:
+
+ Int_t fidDE; // det. elem. ID
+ Int_t fIndex; // det. elem. position index in container
+ Int_t fChamber; // chamber No
+ Double_t fZ; // det. elem. Z-coordinate
+ Int_t fLeft[2]; // numbers of digits not used for clustering
+ Int_t fNHitsForRec; // number of hits for rec.
+ AliMUONGeometrySegmentation* fSeg[2]; // segmentation
+ AliMUONHitMapA1 *fHitMap[2]; // map of digits
+ TObjArray *fDigits[2]; // container of digits from this det. elem.
+ TObjArray *fRawClus; // raw clusters
+ TClonesArray *fHitsForRec; // HitForRec's
+ AliMUONClusterFinderAZ *fRecModel; // cluster finder
+
+ // Functions
+
+ ClassDef(AliMUONDetElement,0) // detection element object
+ };
+#endif
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+/* $Id$ */
+
+#include <Riostream.h>
+#include <TClonesArray.h>
+#include <TArrayS.h>
+#include <TArrayD.h>
+#include "AliMUONEventRecoCombi.h"
+#include "AliMUONData.h"
+#include "AliMUONDetElement.h"
+#include "AliMUONDigit.h"
+#include "AliMUONHitForRec.h"
+#include "AliMUONRawCluster.h"
+#include "AliMUONTrackK.h"
+#include "AliMUONTrackReconstructor.h"
+#include "AliMUONConstants.h"
+#include "AliLoader.h"
+
+ClassImp(AliMUONEventRecoCombi)
+
+AliMUONEventRecoCombi* AliMUONEventRecoCombi::fgRecoCombi = 0;
+
+//_________________________________________________________________________
+AliMUONEventRecoCombi::AliMUONEventRecoCombi() : TObject()
+{
+ // Ctor
+
+ fDetElems = new TClonesArray("AliMUONDetElement", 20);
+ fZ = new TArrayD(20);
+ fNZ = 0;
+ fDEvsZ = NULL;
+}
+
+//_________________________________________________________________________
+AliMUONEventRecoCombi* AliMUONEventRecoCombi::Instance()
+{
+// return pointer to the singleton instance
+
+ if (fgRecoCombi == 0) {
+ fgRecoCombi = new AliMUONEventRecoCombi();
+ }
+ //fDetElems = new TClonesArray("AliMUONDetElement", 20);
+ //fZ = new TArrayD(20);
+ //fNZ = 0;
+ return fgRecoCombi;
+}
+
+//_________________________________________________________________________
+AliMUONEventRecoCombi::~AliMUONEventRecoCombi()
+{
+ // Destructor
+ delete fDetElems;
+ delete fZ;
+ delete [] fDEvsZ;
+}
+
+//_________________________________________________________________________
+void AliMUONEventRecoCombi::FillEvent(AliMUONData *dataCluster, AliMUONData *dataEvent, AliMUONClusterFinderAZ *recModel)
+{
+ // Fill event information
+
+ // Clear previous event
+ fDetElems->Delete();
+ for (Int_t i = 0; i < fNZ; i++) delete [] fDEvsZ[i];
+ delete [] fDEvsZ; fDEvsZ = NULL;
+ fNZ = -1;
+
+ Int_t nDetElem = 0;
+ for (Int_t ich = 0; ich < 6; ich++) {
+ // loop over chambers 0-5
+ TClonesArray *digs = dataCluster->Digits(ich);
+ //cout << ich << " " << digs << " " << digs->GetEntriesFast() << endl;
+ Int_t idDE = -1;
+ for (Int_t i = 0; i < digs->GetEntriesFast(); i++) {
+ AliMUONDigit *dig = (AliMUONDigit*) digs->UncheckedAt(i);
+ if (dig->DetElemId() != idDE) {
+ idDE = dig->DetElemId();
+ new ((*fDetElems)[nDetElem++]) AliMUONDetElement(idDE, dig, recModel);
+ }
+ else ((AliMUONDetElement*)fDetElems->UncheckedAt(nDetElem-1))->AddDigit(dig);
+ }
+ }
+
+ // Sort according to Z
+ fDetElems->Sort();
+ //cout << nDetElem << endl;
+ // Fill det. elems. position index in the container
+ for (Int_t i = 0; i < nDetElem; i++)
+ ((AliMUONDetElement*)fDetElems->UncheckedAt(i))->SetIndex(i);
+
+ // Find groups of det. elements with the same Z
+ Double_t z0 = -99999;
+ TArrayS *nPerZ = new TArrayS(20);
+ for (Int_t i = 0; i < nDetElem; i++) {
+ AliMUONDetElement *detElem = (AliMUONDetElement*) fDetElems->UncheckedAt(i);
+ detElem->Fill(dataCluster);
+ //cout << i << " " << detElem->Z() << endl;
+ if (detElem->Z() - z0 < 0.5) {
+ // the same Z
+ (*nPerZ)[fNZ]++;
+ } else {
+ if (fZ->GetSize() <= fNZ) fZ->Set(fZ->GetSize()+10);
+ if (nPerZ->GetSize() <= fNZ) nPerZ->Set(nPerZ->GetSize()+10);
+ (*fZ)[++fNZ] = detElem->Z();
+ z0 = detElem->Z();
+ (*nPerZ)[fNZ]++;
+ }
+ }
+ fNZ++;
+ /*
+ cout << fNZ << endl;
+ for (Int_t i = 0; i < 7; i++) {
+ cout << i << " " << dataCluster->RawClusters(i)->GetEntriesFast() << endl;
+ }
+ */
+
+ // Build list of DE locations vs Z
+ fDEvsZ = new Int_t* [fNZ];
+ Int_t iPos = 0;
+ for (Int_t i = 0; i < fNZ; i++) {
+ Int_t *idPerZ = new Int_t[(*nPerZ)[i]+1];
+ for (Int_t j = 1; j < (*nPerZ)[i]+1; j++) idPerZ[j] = iPos++;
+ idPerZ[0] = (*nPerZ)[i]; // number of DE's as first element of array
+ fDEvsZ[i] = idPerZ;
+ //cout << (*nPerZ)[i] << " ";
+ }
+ //cout << endl;
+ delete nPerZ;
+
+ // Fill rec. point container for stations 4 and 5
+ //cout << dataEvent->TreeR() << endl;
+ //dataEvent->MakeBranch("RC");
+ dataEvent->SetTreeAddress("RCC");
+ for (Int_t ch = 6; ch < 10; ch++) {
+ TClonesArray *raw = dataCluster->RawClusters(ch);
+ //cout << raw->GetEntriesFast() << " " << dataEvent->RawClusters(ch) << endl;
+ for (Int_t i = 0; i < raw->GetEntriesFast(); i++) {
+ AliMUONRawCluster *clus = (AliMUONRawCluster*) raw->UncheckedAt(i);
+ dataEvent->AddRawCluster(ch, *clus);
+ }
+ }
+ //dataCluster->SetTreeAddress("RC");
+}
+
+//_________________________________________________________________________
+void AliMUONEventRecoCombi::FillRecP(AliMUONData *dataCluster, AliMUONTrackReconstructor *recoTrack)
+{
+ // Fill rec. points used for tracking from det. elems
+
+ TClonesArray *tracks = recoTrack->GetRecTracksPtr();
+ for (Int_t i = 0; i < recoTrack->GetNRecTracks(); i++) {
+ AliMUONTrackK *track = (AliMUONTrackK*) tracks->UncheckedAt(i);
+ TObjArray *hits = track->GetHitOnTrack();
+ for (Int_t j = 0; j < track->GetNTrackHits(); j++) {
+ AliMUONHitForRec *hit = (AliMUONHitForRec*) hits->UncheckedAt(j);
+ if (hit->GetHitNumber() >= 0) continue;
+ // Combined cluster / track finder
+ Int_t index = -hit->GetHitNumber() / 100000;
+ Int_t iPos = -hit->GetHitNumber() - index * 100000;
+ AliMUONRawCluster *clus = (AliMUONRawCluster*) DetElem(index-1)->RawClusters()->UncheckedAt(iPos);
+ //cout << j << " " << iPos << " " << clus << " " << index << " " << DetElem(index-1)->Chamber() << endl;
+ dataCluster->AddRawCluster(DetElem(index-1)->Chamber(), *clus);
+ }
+ }
+ /*
+ for (Int_t ch = 0; ch < 10; ch++) {
+ TClonesArray *raw = dataCluster->RawClusters(ch);
+ cout << ch << " " << raw->GetEntriesFast() << endl;
+ }
+ */
+ // Reset raw cluster tree
+ /*
+ char branchname[30];
+ TBranch * branch = 0x0;
+ if ( dataCluster->TreeR()) {
+ if ( dataCluster->IsTriggerBranchesInTree() ) {
+ // Branch per branch filling
+ for (int i=0; i<AliMUONConstants::NTrackingCh(); i++) {
+ sprintf(branchname,"%sRawClusters%d",dataCluster->GetName(),i+1);
+ branch = dataCluster->TreeR()->GetBranch(branchname);
+ //branch->Fill();
+ //branch->Reset();
+ //branch->Clear();
+ branch->Delete();
+ }
+ }
+ //else TreeR()->Fill();
+ else dataCluster->TreeR()->Reset();
+ }
+ */
+}
+
+//_________________________________________________________________________
+Int_t AliMUONEventRecoCombi::IZfromHit(AliMUONHitForRec *hit) const
+{
+ // Get Iz of det. elem. from the hit
+
+ Int_t index = -hit->GetHitNumber() / 100000 - 1, iz0 = -1;
+ for (Int_t iz = 0; iz < fNZ; iz++) {
+ Int_t *pDEatZ = DEatZ(iz);
+ Int_t nDetElem = pDEatZ[-1];
+ for (Int_t j = 0; j < nDetElem; j++) {
+ if (pDEatZ[j] != index) continue;
+ iz0 = iz;
+ break;
+ }
+ if (iz0 >= 0) break;
+ }
+ return iz0;
+}
--- /dev/null
+#ifndef ALIMUONEVENTRECOCOMBI_H
+#define ALIMUONEVENTRECOCOMBI_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+/* $Id$ */
+
+/// \ingroup rec
+/// \class AliMUONEventRecoCombi
+/// \brief Combined cluster / track finder in MUON arm of ALICE
+
+#include <TObject.h>
+#include <TArrayD.h>
+class TClonesArray;
+class AliMUONData;
+class AliMUONDetElement;
+class AliMUONTrackReconstructor;
+class AliMUONClusterFinderAZ;
+class AliMUONHitForRec;
+class AliLoader;
+
+class AliMUONEventRecoCombi : public TObject
+{
+ public:
+ virtual ~AliMUONEventRecoCombi();
+ static AliMUONEventRecoCombi* Instance();
+ void FillEvent(AliMUONData *dataCluster, AliMUONData *dataEvent, AliMUONClusterFinderAZ *recModel); // fill event info
+ void FillRecP(AliMUONData *dataCluster, AliMUONTrackReconstructor *recoTrack); // fill used rec. points from det. elems
+
+ Int_t Nz() const { return fNZ; } // number of DE different Z-positions
+ Double_t Z(Int_t iz) const { return (*fZ)[iz]; } // Z of DE
+ Int_t *DEatZ(Int_t iz) const { return fDEvsZ[iz]+1; } // list of DE's at Z
+ AliMUONDetElement *DetElem(Int_t iPos) const { return (AliMUONDetElement*) fDetElems->UncheckedAt(iPos); }
+ Int_t IZfromHit(AliMUONHitForRec *hit) const; // IZ from Hit
+
+ protected:
+ AliMUONEventRecoCombi();
+ //AliMUONEventRecoCombi(const AliMUONEventRecoCombi& rhs);
+ //AliMUONEventRecoCombi & operator = (const AliMUONEventRecoCombi& rhs);
+
+ private:
+ static AliMUONEventRecoCombi* fgRecoCombi; // ! singleton instance
+ TClonesArray *fDetElems; // ! array of Det. Elem. objects
+ TArrayD *fZ; // ! array of det. elem. Z-coordinates
+ Int_t fNZ; // ! number of different Z's
+ Int_t **fDEvsZ; // ! list of DE's vs Z-coordinates
+
+ ClassDef(AliMUONEventRecoCombi, 0) // Combined cluster/track finder steering class
+ };
+#endif
#include "AliGenEventHeader.h"
#include "AliESD.h"
#include "AliMUONReconstructor.h"
-
+
#include "AliMUONData.h"
#include "AliMUONTrackReconstructor.h"
#include "AliMUONClusterReconstructor.h"
#include "AliMUONClusterFinderVS.h"
#include "AliMUONClusterFinderAZ.h"
+#include "AliMUONEventRecoCombi.h"
#include "AliMUONTrack.h"
#include "AliMUONTrackParam.h"
#include "AliMUONTriggerTrack.h"
AliMUONTrackReconstructor* recoEvent = new AliMUONTrackReconstructor(loader);
AliMUONData* dataEvent = recoEvent->GetMUONData();
if (strstr(GetOption(),"Kalman")) recoEvent->SetTrackMethod(2); // Kalman
+ else if (strstr(GetOption(),"Combi")) recoEvent->SetTrackMethod(3); // Combined cluster / track
else recoEvent->SetTrackMethod(1); // original
AliMUONClusterReconstructor* recoCluster = new AliMUONClusterReconstructor(loader);
AliMUONData* dataCluster = recoCluster->GetMUONData();
AliMUONClusterFinderVS *recModel = recoCluster->GetRecoModel();
- if (strstr(GetOption(),"AZ")) {
+ if (strstr(GetOption(),"AZ") || strstr(GetOption(),"Combi")) {
recModel = (AliMUONClusterFinderVS*) new AliMUONClusterFinderAZ();
recoCluster->SetRecoModel(recModel);
}
loader->LoadDigits("READ");
loader->LoadRecPoints("RECREATE");
loader->LoadTracks("RECREATE");
-
- // Loop over events
+
+ Int_t chBeg = recoEvent->GetTrackMethod() == 3 ? 6 : 0;
+ // Loop over events
for(Int_t ievent = 0; ievent < nEvents; ievent++) {
printf("Event %d\n",ievent);
runLoader->GetEvent(ievent);
//----------------------- digit2cluster & Trigger2Trigger -------------------
if (!loader->TreeR()) loader->MakeRecPointsContainer();
-
+
// tracking branch
- dataCluster->MakeBranch("RC");
- dataCluster->SetTreeAddress("RC");
- dataCluster->SetTreeAddress("D,");
+ if (recoEvent->GetTrackMethod() != 3) {
+ dataCluster->MakeBranch("RC");
+ dataCluster->SetTreeAddress("D,RC");
+ } else {
+ dataCluster->SetTreeAddress("D");
+ dataCluster->SetTreeAddress("RCC");
+ }
// Important for avoiding a memory leak when reading digits ( to be investigated more in detail)
- // In any case the reading of GLT is need for the Trigger2Tigger method below
+ // In any case the reading of GLT is needed for the Trigger2Tigger method below
dataCluster->SetTreeAddress("GLT");
- recoCluster->Digits2Clusters();
- dataCluster->Fill("RC");
+
+ recoCluster->Digits2Clusters(chBeg);
+ if (recoEvent->GetTrackMethod() == 3) {
+ // Combined cluster / track finder
+ AliMUONEventRecoCombi::Instance()->FillEvent(dataCluster, dataEvent, (AliMUONClusterFinderAZ*)recModel);
+ ((AliMUONClusterFinderAZ*) recModel)->SetReco(2);
+ }
+ else dataCluster->Fill("RC");
// trigger branch
dataCluster->MakeBranch("TC");
dataCluster->SetTreeAddress("TC");
- recoCluster->Trigger2Trigger();
+ recoCluster->Trigger2Trigger();
dataCluster->Fill("TC");
- loader->WriteRecPoints("OVERWRITE");
+ //AZ loader->WriteRecPoints("OVERWRITE");
//---------------------------- Track & TriggerTrack ---------------------
if (!loader->TreeT()) loader->MakeTracksContainer();
recoEvent->EventReconstruct();
dataEvent->Fill("RT");
- loader->WriteTracks("OVERWRITE");
-
+ loader->WriteTracks("OVERWRITE");
+
+ if (recoEvent->GetTrackMethod() == 3) {
+ // Combined cluster / track
+ ((AliMUONClusterFinderAZ*) recModel)->SetReco(1);
+ dataCluster->MakeBranch("RC");
+ dataCluster->SetTreeAddress("RC");
+ AliMUONEventRecoCombi::Instance()->FillRecP(dataCluster, recoEvent);
+ dataCluster->Fill("RC");
+ }
+ loader->WriteRecPoints("OVERWRITE");
+
//--------------------------- Resetting branches -----------------------
dataCluster->ResetDigits();
dataCluster->ResetRawClusters();
dataEvent->ResetRawClusters();
dataEvent->ResetTrigger();
- dataEvent->ResetRecTracks();
+ dataEvent->ResetRecTracks();
dataEvent->ResetRecTriggerTracks();
}
loader->LoadDigits("RECREATE");
- // Loop over events
+ // Loop over events
Int_t iEvent = 0;
-
+
while (rawReader->NextEvent()) {
printf("Event %d\n",iEvent);
runLoader->GetEvent(iEvent++);
dataCluster->MakeBranch("D");
dataCluster->SetTreeAddress("D");
rawData->ReadTrackerDDL(rawReader);
- dataCluster->Fill("D");
+ dataCluster->Fill("D");
// trigger branch
dataCluster->MakeBranch("GLT");
dataCluster->SetTreeAddress("GLT");
rawData->ReadTriggerDDL(rawReader);
- dataCluster->Fill("GLT");
+ dataCluster->Fill("GLT");
loader->WriteDigits("OVERWRITE");
//----------------------- digit2cluster & Trigger2Trigger -------------------
if (!loader->TreeR()) loader->MakeRecPointsContainer();
-
+
// tracking branch
dataCluster->MakeBranch("RC");
dataCluster->SetTreeAddress("RC");
- recoCluster->Digits2Clusters();
- dataCluster->Fill("RC");
+ recoCluster->Digits2Clusters();
+ dataCluster->Fill("RC");
// trigger branch
dataCluster->MakeBranch("TC");
dataCluster->SetTreeAddress("TC");
- recoCluster->Trigger2Trigger();
+ recoCluster->Trigger2Trigger();
dataCluster->Fill("TC");
loader->WriteRecPoints("OVERWRITE");
recoEvent->EventReconstruct();
dataEvent->Fill("RT");
- loader->WriteTracks("OVERWRITE");
-
+ loader->WriteTracks("OVERWRITE");
+
//--------------------------- Resetting branches -----------------------
dataCluster->ResetDigits();
dataCluster->ResetRawClusters();
dataEvent->ResetTrigger();
dataEvent->ResetRecTracks();
dataEvent->ResetRecTriggerTracks();
-
+
}
loader->UnloadRecPoints();
loader->UnloadTracks();
{
TClonesArray* recTracksArray = 0;
TClonesArray* recTrigTracksArray = 0;
-
+
AliLoader* loader = runLoader->GetLoader("MUONLoader");
loader->LoadTracks("READ");
AliMUONData* muonData = new AliMUONData(loader,"MUON","MUON");
- // declaration
+ // declaration
Int_t iEvent;// nPart;
Int_t nTrackHits;// nPrimary;
Double_t fitFmin;
AliMUONTriggerTrack* recTriggerTrack = 0;
// TParticle* particle = new TParticle();
// AliGenEventHeader* header = 0;
- iEvent = runLoader->GetEventNumber();
+ iEvent = runLoader->GetEventNumber();
runLoader->GetEvent(iEvent);
// vertex calculation (maybe it exists already somewhere else)
recTrigTracksArray = muonData->RecTriggerTracks();
// ready global trigger pattern from first track
- if (recTrigTracksArray)
+ if (recTrigTracksArray)
recTriggerTrack = (AliMUONTriggerTrack*) recTrigTracksArray->First();
if (recTriggerTrack) trigPat = recTriggerTrack->GetGTPattern();
//printf(">>> Event %d Number of Recconstructed tracks %d \n",iEvent, nrectracks);
-
+
// -------------------- tracks-------------
muonData->SetTreeAddress("RT");
muonData->GetRecTracks();
recTracksArray = muonData->RecTracks();
-
+
Int_t nRecTracks = 0;
if (recTracksArray)
nRecTracks = (Int_t) recTracksArray->GetEntriesFast(); //
-
+
// loop over tracks
for (Int_t iRecTracks = 0; iRecTracks < nRecTracks; iRecTracks++) {
theESDTrack->SetMatchTrigger(matchTrigger);
theESDTrack->SetChi2MatchTrigger(chi2MatchTrigger);
- // storing ESD MUON Track into ESD Event
- if (nRecTracks != 0)
+ // storing ESD MUON Track into ESD Event
+ if (nRecTracks != 0)
esd->AddMuonTrack(theESDTrack);
} // end loop tracks
// add global trigger pattern
- if (nRecTracks != 0)
+ if (nRecTracks != 0)
esd->SetTrigger(trigPat);
// reset muondata
muonData->ResetRecTracks();
muonData->ResetRecTriggerTracks();
- //} // end loop on event
- loader->UnloadTracks();
+ //} // end loop on event
+ loader->UnloadTracks();
// if (!header)
// runLoader->UnloadKinematics();
delete theESDTrack;
#include <TMatrixD.h>
#include "AliMUONTrackK.h"
-#include "AliCallf77.h"
#include "AliMUON.h"
+
#include "AliMUONTrackReconstructor.h"
#include "AliMagF.h"
#include "AliMUONSegment.h"
#include "AliMUONHitForRec.h"
#include "AliMUONRawCluster.h"
#include "AliMUONTrackParam.h"
+#include "AliMUONEventRecoCombi.h"
+#include "AliMUONDetElement.h"
#include "AliRun.h"
#include "AliLog.h"
const Double_t AliMUONTrackK::fgkEpsilon = 0.002;
void mnvertLocal(Double_t* a, Int_t l, Int_t m, Int_t n, Int_t& ifail); // from AliMUONTrack
-//void mnvertLocalK(Double_t* a, Int_t l, Int_t m, Int_t n, Int_t& ifail);
ClassImp(AliMUONTrackK) // Class implementation in ROOT context
-Int_t AliMUONTrackK::fgDebug = -1;
+ Int_t AliMUONTrackK::fgDebug = -1; //-1;
Int_t AliMUONTrackK::fgNOfPoints = 0;
-AliMUON* AliMUONTrackK::fgMUON = NULL;
AliMUONTrackReconstructor* AliMUONTrackK::fgTrackReconstructor = NULL;
TClonesArray* AliMUONTrackK::fgHitForRec = NULL;
+AliMUONEventRecoCombi *AliMUONTrackK::fgCombi = NULL;
//FILE *lun1 = fopen("window.dat","w");
//__________________________________________________________________________
// Default constructor
fgTrackReconstructor = NULL; // pointer to event reconstructor
- fgMUON = NULL; // pointer to Muon module
fgHitForRec = NULL; // pointer to points
fgNOfPoints = 0; // number of points
//__________________________________________________________________________
AliMUONTrackK::AliMUONTrackK(AliMUONTrackReconstructor *TrackReconstructor, TClonesArray *hitForRec)
- //AZ: TObject()
: AliMUONTrack()
{
// Constructor
if (!TrackReconstructor) return;
fgTrackReconstructor = TrackReconstructor; // pointer to event reconstructor
- fgMUON = (AliMUON*) gAlice->GetModule("MUON"); // pointer to Muon module
fgHitForRec = hitForRec; // pointer to points
fgNOfPoints = fgHitForRec->GetEntriesFast(); // number of points
+ fgCombi = NULL;
+ if (fgTrackReconstructor->GetTrackMethod() == 3) fgCombi = AliMUONEventRecoCombi::Instance();
fStartSegment = NULL;
fTrackHitsPtr = NULL;
//__________________________________________________________________________
AliMUONTrackK::AliMUONTrackK(AliMUONSegment *segment)
- //AZ: TObject()
: AliMUONTrack(segment, segment, fgTrackReconstructor)
{
// Constructor from a segment
fParSmooth = NULL;
if (fgDebug < 0 ) return;
- cout << fgTrackReconstructor->GetBendingMomentumFromImpactParam(segment->GetBendingImpact()) << " " << 1/(*fTrackPar)(4,0) << " ";
+ cout << fgTrackReconstructor->GetNRecTracks()-1 << " " << fgTrackReconstructor->GetBendingMomentumFromImpactParam(segment->GetBendingImpact()) << " " << 1/(*fTrackPar)(4,0) << " ";
if (fgTrackReconstructor->GetRecTrackRefHits()) {
// from track ref. hits
cout << ((AliMUONHitForRec*)((*fTrackHitsPtr)[0]))->GetTTRTrack() << "<-->" << ((AliMUONHitForRec*)((*fTrackHitsPtr)[1]))->GetTTRTrack() << " @ " << fStartSegment->GetHitForRec1()->GetChamberNumber() << endl;
//__________________________________________________________________________
AliMUONTrackK::AliMUONTrackK (const AliMUONTrackK& source)
- //AZ: TObject(source)
: AliMUONTrack(source)
{
// Protected copy constructor
(*fWeight)(1,3) = dBdX*(*fWeight)(1,1); // <xb>
(*fWeight)(3,1) = (*fWeight)(1,3);
- //(*fWeight)(4,4) = ((*fTrackPar)(4,0)*0.2)*((*fTrackPar)(4,0)*0.2); // error 20%
(*fWeight)(4,4) = ((*fTrackPar)(4,0)*0.5)*((*fTrackPar)(4,0)*0.5); // error 50%
// check whether the Invert method returns flag if matrix cannot be inverted,
if (fWeight->Determinant() != 0) {
// fWeight->Invert();
Int_t ifail;
- //mnvertLocalK(&((*fWeight)(0,0)), fgkSize,fgkSize,fgkSize,ifailWeight);
mnvertLocal(&((*fWeight)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
} else {
AliWarning(" Determinant fWeight=0:");
// Follows track through detector stations
Bool_t miss, success;
Int_t ichamb, iFB, iMin, iMax, dChamb, ichambOK;
- Int_t ihit, firstIndx, lastIndx, currIndx;
+ Int_t ihit, firstIndx = -1, lastIndx = -1, currIndx = -1, iz0 = -1, iz = -1;
Double_t zEnd, dChi2;
- AliMUONHitForRec *hitAdd, *firstHit, *lastHit, *hit;
+ AliMUONHitForRec *hitAdd, *firstHit = NULL, *lastHit = NULL, *hit;
AliMUONRawCluster *clus;
TClonesArray *rawclusters;
hit = 0; clus = 0; rawclusters = 0;
fRecover = 0;
ichambOK = ((AliMUONHitForRec*)((*fTrackHitsPtr)[fNTrackHits-1]))->GetChamberNumber();
//if (ichambOK >= 8 && ((AliMUONHitForRec*)((*fTrackHitsPtr)[0]))->GetChamberNumber() == 6) ichambOK = 6;
- currIndx = fgHitForRec->IndexOf(fSkipHit);
+ if (fgTrackReconstructor->GetTrackMethod() == 3 &&
+ fSkipHit->GetHitNumber() < 0) {
+ iz0 = fgCombi->IZfromHit(fSkipHit);
+ currIndx = -1;
+ }
+ else currIndx = fgHitForRec->IndexOf(fSkipHit);
} else {
// outlier
fTrackHitsPtr->Compress();
// Get indices of the 1'st and last hits on the track candidate
firstHit = (AliMUONHitForRec*) fTrackHitsPtr->First();
lastHit = (AliMUONHitForRec*) fTrackHitsPtr->Last();
- firstIndx = fgHitForRec->IndexOf(firstHit);
- lastIndx = fgHitForRec->IndexOf(lastHit);
- currIndx = TMath::Abs (TMath::Max(firstIndx*iFB,lastIndx*iFB));
+ if (fgTrackReconstructor->GetTrackMethod() == 3 &&
+ lastHit->GetHitNumber() < 0) iz0 = fgCombi->IZfromHit(lastHit);
+ else {
+ firstIndx = fgHitForRec->IndexOf(firstHit);
+ lastIndx = fgHitForRec->IndexOf(lastHit);
+ currIndx = TMath::Abs (TMath::Max(firstIndx*iFB,lastIndx*iFB));
+ }
}
- while (ichamb>=iMin && ichamb<=iMax) {
+ if (iz0 < 0) iz0 = iFB;
+ while (ichamb >= iMin && ichamb <= iMax) {
// Find the closest hit in Z, not belonging to the current plane
if (Back) {
// backpropagation
break;
}
}
+
+ // Combined cluster / track finder
+ if (zEnd > 999 && iFB < 0 && fgTrackReconstructor->GetTrackMethod() == 3) {
+ currIndx = -2;
+ AliMUONHitForRec hitTmp;
+ for (iz = iz0 - iFB; iz < fgCombi->Nz(); iz++) {
+ if (TMath::Abs(fgCombi->Z(iz)-fPosition) < 0.5) continue;
+ Int_t *pDEatZ = fgCombi->DEatZ(iz);
+ //cout << iz << " " << fgCombi->Z(iz) << endl;
+ zEnd = fgCombi->Z(iz);
+ iz0 = iz;
+ AliMUONDetElement *detElem = fgCombi->DetElem(pDEatZ[0]);
+ hitAdd = &hitTmp;
+ hitAdd->SetChamberNumber(detElem->Chamber());
+ //hitAdd = (AliMUONHitForRec*) detElem->HitsForRec()->First();
+ if (hitAdd) break;
+ }
+ }
}
- //AZ(Z->-Z) if (zEnd<-999 && ichamb==ichamEnd) endOfProp = 1; // end-of-propagation
if (zEnd>999 && ichamb==ichamEnd) endOfProp = 1; // end-of-propagation
else {
// Check if there is a chamber without hits
Back = kFALSE;
fRecover =0;
ichambOK = ((AliMUONHitForRec*)((*fTrackHitsPtr)[fNTrackHits-1]))->GetChamberNumber();
+ //? if (fgTrackReconstructor->GetTrackMethod() != 3) currIndx = fgHitForRec->IndexOf(fSkipHit);
currIndx = fgHitForRec->IndexOf(fSkipHit);
}
TryPoint(point,pointWeight,trackParTmp,dChi2);
*fTrackPar = trackParTmp;
*fWeight += pointWeight;
- //AZ(z->-z) if (fTrackDir < 0) AddMatrices (this, dChi2, hitAdd);
if (fTrackDir > 0) AddMatrices (this, dChi2, hitAdd);
fChi2 += dChi2; // Chi2
} else *fTrackPar = *fTrackParNew; // adjust end point to the last hit
currIndx ++;
} else {
// forward propagator
- if (miss || !FindPoint(ichamb,zEnd,currIndx,iFB,hitAdd)) {
+ if (miss || !FindPoint(ichamb, zEnd, currIndx, iFB, hitAdd, iz)) {
// missing point
*fTrackPar = *fTrackParNew;
} else {
//fTrackParNew->Print();
return;
}
-/*
- //__________________________________________________________________________
-void AliMUONTrackK::WeightPropagation(void)
-{
- // Propagation of the weight matrix
- // W = DtWD, where D is Jacobian
- // !!! not implemented TMatrixD weight1(*fJacob,TMatrixD::kAtBA,*fWeight); // DtWD
- TMatrixD weight1(*fWeight,TMatrixD::kMult,*fJacob); // WD
- *fWeight = TMatrixD(*fJacob,TMatrixD::kTransposeMult,weight1); // DtWD
- return;
-}
-*/
//__________________________________________________________________________
void AliMUONTrackK::WeightPropagation(Double_t zEnd, Bool_t smooth)
{
// check whether the Invert method returns flag if matrix cannot be inverted,
// and do not calculate the Determinant in that case !!!!
if (fCovariance->Determinant() != 0) {
- // fCovariance->Invert();
Int_t ifail;
mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
//fCovariance->Print();
// Save for smoother
if (!smooth) return; // do not use smoother
- //AZ(z->-z) if (fTrackDir > 0) return; // only for propagation towards int. point
if (fTrackDir < 0) return; // only for propagation towards int. point
TMatrixD *tmp = new TMatrixD(*fTrackParNew); // extrapolated parameters
fParExtrap->Add(tmp);
}
//__________________________________________________________________________
-Bool_t AliMUONTrackK::FindPoint(Int_t ichamb, Double_t zEnd, Int_t currIndx, Int_t iFB, AliMUONHitForRec *&hitAdd)
+Bool_t AliMUONTrackK::FindPoint(Int_t ichamb, Double_t zEnd, Int_t currIndx, Int_t iFB, AliMUONHitForRec *&hitAdd, Int_t iz)
{
// Picks up point within a window for the chamber No ichamb
// Split the track if there are more than 1 hit
TClonesArray *trackPtr;
AliMUONHitForRec *hit, *hitLoop;
AliMUONTrackK *trackK;
+ AliMUONDetElement *detElem = NULL;
Bool_t ok = kFALSE;
+
+ if (fgTrackReconstructor->GetTrackMethod() == 3 && iz >= 0) {
+ // Combined cluster / track finder
+ // Check if extrapolated track passes thru det. elems. at iz
+ Int_t *pDEatZ = fgCombi->DEatZ(iz);
+ Int_t nDetElem = pDEatZ[-1];
+ //cout << fgCombi->Z(iz) << " " << nDetElem << endl;
+ for (Int_t i = 0; i < nDetElem; i++) {
+ detElem = fgCombi->DetElem(pDEatZ[i]);
+ if (detElem->Inside((*fTrackParNew)(1,0), (*fTrackParNew)(0,0), fPosition)) {
+ detElem->ClusterReco((*fTrackParNew)(1,0), (*fTrackParNew)(0,0));
+ hitAdd = (AliMUONHitForRec*) detElem->HitsForRec()->First();
+ ok = kTRUE;
+ break;
+ }
+ }
+ if (!ok) return ok; // outside det. elem.
+ ok = kFALSE;
+ }
+
//sigmaB = fgTrackReconstructor->GetBendingResolution(); // bending resolution
//sigmaNonB = fgTrackReconstructor->GetNonBendingResolution(); // non-bending resolution
*fCovariance = *fWeight;
// check whether the Invert method returns flag if matrix cannot be inverted,
// and do not calculate the Determinant in that case !!!!
if (fCovariance->Determinant() != 0) {
- // fCovariance->Invert();
Int_t ifail;
mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
} else {
TMatrixD point(fgkSize,1);
TMatrixD trackPar = point;
TMatrixD trackParTmp = point;
- Int_t nHitsOK = 0;
+ Int_t nHitsOK = 0, ihitB = currIndx, ihitE = fgNOfPoints, iDhit = iFB;
Double_t zLast;
zLast = ((AliMUONHitForRec*)fTrackHitsPtr->Last())->GetZ();
+ if (fgTrackReconstructor->GetTrackMethod() == 3 && detElem) {
+ ihitB = 0;
+ ihitE = detElem->NHitsForRec();
+ iDhit = 1;
+ }
+
- for (ihit=currIndx; ihit>=0 && ihit<fgNOfPoints; ihit+=iFB) {
- hit = (AliMUONHitForRec*) ((*fgHitForRec)[ihit]);
+ for (ihit = ihitB; ihit >= 0 && ihit < ihitE; ihit+=iDhit) {
+ if (fgTrackReconstructor->GetTrackMethod() != 3 || !detElem)
+ hit = (AliMUONHitForRec*) ((*fgHitForRec)[ihit]);
+ else hit = (AliMUONHitForRec*) detElem->HitsForRec()->UncheckedAt(ihit);
if (hit->GetChamberNumber() == ichamb) {
//if (TMath::Abs(hit->GetZ()-zEnd) < 0.1) {
if (TMath::Abs(hit->GetZ()-zEnd) < 0.5) {
}
y = hit->GetBendingCoor();
x = hit->GetNonBendingCoor();
+ if (hit->GetBendingReso2() < 0) {
+ // Combined cluster / track finder
+ hit->SetBendingReso2(fgTrackReconstructor->GetBendingResolution()*
+ fgTrackReconstructor->GetBendingResolution());
+ hit->SetNonBendingReso2(fgTrackReconstructor->GetNonBendingResolution()*
+ fgTrackReconstructor->GetNonBendingResolution());
+ }
windowB = fgkNSigma*TMath::Sqrt((*fCovariance)(0,0)+hit->GetBendingReso2());
windowNonB = fgkNSigma*TMath::Sqrt((*fCovariance)(1,1)+hit->GetNonBendingReso2());
pointWeight(1,1) = 1/hit->GetNonBendingReso2();
TryPoint(point,pointWeight,trackPar,dChi2);
if (TMath::Abs(1./(trackPar)(4,0)) < fgTrackReconstructor->GetMinBendingMomentum()) continue; // p < p_min - next hit
+ // if (TMath::Sign(1.,(trackPar)(4,0)*(*fTrackPar)(4,0)) < 0) continue; // change of sign
ok = kTRUE;
nHitsOK++;
//if (nHitsOK > -1) {
cout << (*(trackK->fCovariance))(0,0) << " " << (*(trackK->fCovariance))(1,1) << " " << (*fCovariance)(0,0) << " " << (*fCovariance)(1,1) << endl;
*/
// Add filtered matrices for the smoother
- //AZ(z->-z) if (fTrackDir < 0) {
if (fTrackDir > 0) {
if (nHitsOK > 2) { // check for position adjustment
for (Int_t i=trackK->fNSteps-1; i>=0; i--) {
hit->SetNTrackHits(hit->GetNTrackHits()+1); // mark hit as being on track
if (ichamb == 9) {
// the last chamber
- //AZ(z->-z) trackK->fTrackDir = -1;
trackK->fTrackDir = 1;
trackK->fBPFlag = kTRUE;
}
fChi2 += dChi2Tmp; // Chi2
fPosition = savePosition;
// Add filtered matrices for the smoother
- //AZ(z->-z) if (fTrackDir < 0) {
if (fTrackDir > 0) {
for (Int_t i=fNSteps-1; i>=0; i--) {
if (TMath::Abs(fPosition-(*fSteps)[i]) > 0.1) {
// check whether the Invert method returns flag if matrix cannot be inverted,
// and do not calculate the Determinant in that case !!!!
if (wu.Determinant() != 0) {
- // wu.Invert();
- Int_t ifail;
- mnvertLocal(&((wu)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
+ Int_t ifail;
+ mnvertLocal(&((wu)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
} else {
AliWarning(" Determinant wu=0:");
}
// check whether the Invert method returns flag if matrix cannot be inverted,
// and do not calculate the Determinant in that case !!!!
if (fWeight->Determinant() != 0) {
- //fWeight->Invert(); // covariance
-
Int_t ifail;
mnvertLocal(&((*fWeight)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
} else {
(*fWeight)(2,2) += sign*theta0/cosBeta*theta0/cosBeta; // alpha
(*fWeight)(3,3) += sign*theta0*theta0; // beta
- //fWeight->Invert(); // weight
-
Int_t ifail;
mnvertLocal(&((*fWeight)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
return;
// Get covariance matrix
*fCovariance = *fWeight;
if (fCovariance->Determinant() != 0) {
- Int_t ifail;
- mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
+ Int_t ifail;
+ mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
} else {
AliWarning(" Determinant fCovariance=0:");
}
// from raw clusters
for (Int_t i1=0; i1<fNTrackHits; i1++) {
hit = (AliMUONHitForRec*) ((*fTrackHitsPtr)[i1]);
- rawclusters = fgTrackReconstructor->GetMUONData()->RawClusters(hit->GetChamberNumber());
- clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->GetHitNumber());
+ if (hit->GetHitNumber() < 0) { // combined cluster / track finder
+ Int_t index = -hit->GetHitNumber() / 100000;
+ Int_t iPos = -hit->GetHitNumber() - index * 100000;
+ clus = (AliMUONRawCluster*) fgCombi->DetElem(index-1)->RawClusters()->UncheckedAt(iPos);
+ } else {
+ rawclusters = fgTrackReconstructor->GetMUONData()->RawClusters(hit->GetChamberNumber());
+ clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->GetHitNumber());
+ }
printf ("%4d", clus->GetTrack(1));
}
cout << endl;
for (Int_t i1=0; i1<fNTrackHits; i1++) {
hit = (AliMUONHitForRec*) ((*fTrackHitsPtr)[i1]);
- rawclusters = fgTrackReconstructor->GetMUONData()->RawClusters(hit->GetChamberNumber());
- clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->GetHitNumber());
+ if (hit->GetHitNumber() < 0) { // combined cluster / track finder
+ Int_t index = -hit->GetHitNumber() / 100000;
+ Int_t iPos = -hit->GetHitNumber() - index * 100000;
+ clus = (AliMUONRawCluster*) fgCombi->DetElem(index-1)->RawClusters()->UncheckedAt(iPos);
+ } else {
+ rawclusters = fgTrackReconstructor->GetMUONData()->RawClusters(hit->GetChamberNumber());
+ clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->GetHitNumber());
+ }
if (clus->GetTrack(2) != -1) printf ("%4d", clus->GetTrack(2));
else printf ("%4s", " ");
}
/* Not needed - covariance matrix is not interesting to anybody
*fCovariance = *fWeight;
if (fCovariance->Determinant() != 0) {
- // fCovariance->Invert();
Int_t ifail;
mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
} else {
// Get weight matrix
*fWeight = *fCovariance;
if (fWeight->Determinant() != 0) {
- // fWeight->Invert();
Int_t ifail;
mnvertLocal(&((*fWeight)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
} else {
trackK = new ((*trackPtr)[nRecTracks]) AliMUONTrackK(NULL, NULL);
*trackK = *this;
fgTrackReconstructor->SetNRecTracks(nRecTracks+1);
- //AZ(z->-z) trackK->fTrackDir = -1;
trackK->fTrackDir = 1;
trackK->fRecover = 2;
trackK->fSkipHit = hit;
/// \brief Kalman track in MUON arm of ALICE
#include <TObject.h>
-#include "AliMUONTrack.h" //AZ
+#include "AliMUONTrack.h"
class TArrayD;
class TMatrixD;
class TClonesArray;
class TObjArray;
class AliMUONSegment;
-class AliMUON;
class AliMUONHitForRec;
+class AliMUONEventRecoCombi;
-//AZ class AliMUONTrackK : public TObject {
class AliMUONTrackK : public AliMUONTrack {
public:
AliMUONTrackK(); // Default constructor
virtual ~AliMUONTrackK(); // Destructor
- //AliMUONTrackK(const AliMUONTrackReconstructor *TrackReconstructor, const AliMUONHitForRec *hitForRec); // Constructor
AliMUONTrackK(AliMUONTrackReconstructor *TrackReconstructor, TClonesArray *hitForRec); // Constructor
AliMUONTrackK(AliMUONSegment *segment); // Constructor from a segment
static Int_t fgDebug; // debug level
static Int_t fgNOfPoints; // number of points in event
- static AliMUON *fgMUON; // pointer to MUON module
+ //static AliMUON *fgMUON; // pointer to MUON module
static AliMUONTrackReconstructor *fgTrackReconstructor; // pointer to event reconstructor
static TClonesArray *fgHitForRec; // pointer to hits
+ static AliMUONEventRecoCombi *fgCombi; // pointer to combined cluster/track finder
AliMUONSegment *fStartSegment; // seed segment
Double_t fPosition; // Z-coordinate of track
void WeightPropagation(Double_t zEnd, Bool_t smooth);
void MSThin(Int_t sign);
void MSLine(Double_t dZ, Double_t X0);
- Bool_t FindPoint(Int_t ichamb, Double_t zEnd, Int_t currIndx, Int_t iFB, AliMUONHitForRec *&hitAdd);
+ Bool_t FindPoint(Int_t ichamb, Double_t zEnd, Int_t currIndx, Int_t iFB, AliMUONHitForRec *&hitAdd, Int_t iz);
void TryPoint(TMatrixD &point, const TMatrixD &pointWeight, TMatrixD &trackParTmp, Double_t &dChi2);
void SetGeantParam(Double_t *VGeant3, Int_t iFB);
void GetFromGeantParam(Double_t *VGeant3, Int_t iFB);
* provided "as is" without express or implied warranty. *
**************************************************************************/
+/* $Id$ */
+
////////////////////////////////////
//
// MUON track reconstructor in ALICE (class renamed from AliMUONEventReconstructor)
#include <Riostream.h>
#include <TDirectory.h>
#include <TFile.h>
-#include <TMatrixD.h> //AZ
+#include <TMatrixD.h>
#include "AliMUONTrackReconstructor.h"
#include "AliMUON.h"
#include "AliRun.h" // for gAlice
#include "AliRunLoader.h"
#include "AliLoader.h"
-#include "AliMUONTrackK.h" //AZ
+#include "AliMUONTrackK.h"
#include "AliMC.h"
#include "AliLog.h"
#include "AliTrackReference.h"
fNRecTracks = 0; // really needed or GetEntriesFast sufficient ????
// Memory allocation for the TClonesArray of hits on reconstructed tracks
// Is 100 the right size ????
- //AZ fRecTrackHitsPtr = new TClonesArray("AliMUONTrack", 100);
fRecTrackHitsPtr = new TClonesArray("AliMUONTrackHit", 100);
fNRecTrackHits = 0; // really needed or GetEntriesFast sufficient ????
// TreeR assumed to be be "prepared" in calling function
// by "MUON->GetTreeR(nev)" ????
TTree *treeR = fLoader->TreeR();
- fMUONData->SetTreeAddress("RC");
+ //AZ? fMUONData->SetTreeAddress("RC");
AddHitsForRecFromRawClusters(treeR);
// No sorting: it is done automatically in the previous function
}
TClonesArray *rawclusters;
AliDebug(1,"Enter AddHitsForRecFromRawClusters");
- nTRentries = Int_t(TR->GetEntries());
- if (nTRentries != 1) {
- AliError(Form("nTRentries = %d not equal to 1 ",nTRentries));
- exit(0);
+ if (fTrackMethod != 3) { //AZ
+ fMUONData->SetTreeAddress("RC"); //AZ
+ nTRentries = Int_t(TR->GetEntries());
+ if (nTRentries != 1) {
+ AliError(Form("nTRentries = %d not equal to 1 ",nTRentries));
+ exit(0);
+ }
+ fLoader->TreeR()->GetEvent(0); // only one entry
}
- fLoader->TreeR()->GetEvent(0); // only one entry
// Loop over tracking chambers
for (Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) {
// resolution: info should be already in raw cluster and taken from it ????
hitForRec->SetBendingReso2(fBendingResolution * fBendingResolution);
hitForRec->SetNonBendingReso2(fNonBendingResolution * fNonBendingResolution);
+ //hitForRec->SetBendingReso2(clus->GetErrY() * clus->GetErrY());
+ //hitForRec->SetNonBendingReso2(clus->GetErrX() * clus->GetErrX());
// original raw cluster
hitForRec->SetChamberNumber(ch);
hitForRec->SetHitNumber(iclus);
hitForRec->Dump();}
} // end of cluster loop
} // end of chamber loop
- SortHitsForRecWithIncreasingChamber(); //AZ
+ SortHitsForRecWithIncreasingChamber();
return;
}
AliDebug(1,"Enter MakeSegments");
ResetSegments();
// Loop over stations
- Int_t nb = (fTrackMethod == 2) ? 3 : 0; //AZ
- //AZ for (Int_t st = 0; st < fgkMaxMuonTrackingStations; st++)
- for (Int_t st = nb; st < AliMUONConstants::NTrackingCh()/2; st++) //AZ
+ Int_t nb = (fTrackMethod != 1) ? 3 : 0; //AZ
+ for (Int_t st = nb; st < AliMUONConstants::NTrackingCh()/2; st++)
MakeSegmentsPerStation(st);
if (AliLog::GetGlobalDebugLevel() > 1) {
cout << "end of MakeSegments" << endl;
// The order may be important for the following Reset's
ResetTracks();
ResetTrackHits();
- if (fTrackMethod == 2) { //AZ - Kalman filter
+ if (fTrackMethod != 1) { //AZ - Kalman filter
MakeTrackCandidatesK();
if (fRecTracksPtr->GetEntriesFast() == 0) return;
// Follow tracks in stations(1..) 3, 2 and 1
Int_t ch = -1, chInStation, chBestHit = -1, iHit, iSegment, station, trackIndex;
Double_t bestChi2, chi2, dZ1, dZ2, dZ3, maxSigma2Distance, mcsFactor;
Double_t bendingMomentum, chi2Norm = 0.;
+
// local maxSigma2Distance, for easy increase in testing
maxSigma2Distance = fMaxSigma2Distance;
AliDebug(2,"Enter FollowTracks");
for (iseg=0; iseg<fNSegments[istat]; iseg++) {
// Transform segments to tracks and evaluate covariance matrix
segment = (AliMUONSegment*) ((*fSegmentsPtr[istat])[iseg]);
- trackK = new ((*fRecTracksPtr)[fNRecTracks]) AliMUONTrackK(segment);
- fNRecTracks++;
+ trackK = new ((*fRecTracksPtr)[fNRecTracks++]) AliMUONTrackK(segment);
} // for (iseg=0;...)
} // for (istat=4;...)
return;
TClonesArray *rawclusters;
clus = 0; rawclusters = 0;
- //AZ(z->-z) zDipole1 = GetSimpleBPosition() - GetSimpleBLength()/2;
- //AZ(z->-z) zDipole2 = zDipole1 + GetSimpleBLength();
zDipole1 = GetSimpleBPosition() + GetSimpleBLength()/2;
zDipole2 = zDipole1 - GetSimpleBLength();
rawclusters = fMUONData->RawClusters(hit->GetChamberNumber());
clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->
GetHitNumber());
- printf("%3d", clus->GetTrack(1)-1);
- if (clus->GetTrack(2) != 0) printf("%3d \n", clus->GetTrack(2)-1);
+ printf(" %d", clus->GetTrack(1));
+ if (clus->GetTrack(2) != -1) printf(" %d \n", clus->GetTrack(2));
else printf("\n");
} // if (fRecTrackRefHits)
}
if (hit) ichamBeg = hit->GetChamberNumber();
ichamEnd = 0;
// Check propagation direction
- if (hit == NULL) ichamBeg = ichamEnd;
- //AZ(z->-z) else if (trackK->GetTrackDir() > 0) {
+ if (!hit) { ichamBeg = ichamEnd; AliFatal(" ??? "); }
else if (trackK->GetTrackDir() < 0) {
ichamEnd = 9; // forward propagation
ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kFALSE,zDipole1,zDipole2);
ichamEnd = 6; // backward propagation
// Change weight matrix and zero fChi2 for backpropagation
trackK->StartBack();
- //AZ trackK->SetTrackDir(-1);
trackK->SetTrackDir(1);
ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kTRUE,zDipole1,zDipole2);
ichamBeg = ichamEnd;
}
if (ok) {
- //AZ trackK->SetTrackDir(-1);
trackK->SetTrackDir(1);
trackK->SetBPFlag(kFALSE);
ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kFALSE,zDipole1,zDipole2);
{
// Set track method and recreate track container if necessary
- fTrackMethod = iTrackMethod == 2 ? 2 : 1;
- if (fTrackMethod == 2) {
+ fTrackMethod = TMath::Min (iTrackMethod, 3);
+ fTrackMethod = TMath::Max (fTrackMethod, 1);
+ if (fTrackMethod != 1) {
if (fRecTracksPtr) delete fRecTracksPtr;
fRecTracksPtr = new TClonesArray("AliMUONTrackK", 10);
- cout << " *** Tracking with the Kalman filter *** " << endl;
+ if (fTrackMethod == 2) cout << " *** Tracking with the Kalman filter *** " << endl;
+ else cout << " *** Combined cluster / track finder ***" << endl;
} else cout << " *** Traditional tracking *** " << endl;
}
#pragma link C++ class AliMUONTrackHit+;
#pragma link C++ class AliMUONSegment+;
#pragma link C++ class AliMUONClusterDrawAZ+;
+#pragma link C++ class AliMUONDetElement+;
+#pragma link C++ class AliMUONEventRecoCombi+;
// raw data
#pragma link C++ class AliMUONDDLTrigger+;
AliMUONSubEventTracker.cxx \
AliMUONRawData.cxx \
AliMUONRawStream.cxx \
- AliMUONClusterDrawAZ.cxx
+ AliMUONClusterDrawAZ.cxx \
+ AliMUONDetElement.cxx \
+ AliMUONEventRecoCombi.cxx
HDRS:= $(SRCS:.cxx=.h)
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff