Modifications needed for the combined cluster / track finder (Sacha)
authormartinez <martinez@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 28 Nov 2005 12:17:28 +0000 (12:17 +0000)
committermartinez <martinez@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 28 Nov 2005 12:17:28 +0000 (12:17 +0000)
15 files changed:
MUON/AliMUONClusterFinderAZ.cxx
MUON/AliMUONClusterFinderAZ.h
MUON/AliMUONClusterReconstructor.cxx
MUON/AliMUONClusterReconstructor.h
MUON/AliMUONData.cxx
MUON/AliMUONDetElement.cxx [new file with mode: 0644]
MUON/AliMUONDetElement.h [new file with mode: 0644]
MUON/AliMUONEventRecoCombi.cxx [new file with mode: 0644]
MUON/AliMUONEventRecoCombi.h [new file with mode: 0644]
MUON/AliMUONReconstructor.cxx
MUON/AliMUONTrackK.cxx
MUON/AliMUONTrackK.h
MUON/AliMUONTrackReconstructor.cxx
MUON/MUONrecLinkDef.h
MUON/libMUONrec.pkg

index 3ed43b5..0195add 100644 (file)
@@ -56,11 +56,13 @@ AliMUONClusterFinderAZ::AliMUONClusterFinderAZ(Bool_t draw)
   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;
 }
 
 //_____________________________________________________________________________
@@ -104,20 +106,24 @@ void AliMUONClusterFinderAZ::EventLoop(Int_t nev, Int_t ch)
   //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;
@@ -128,7 +134,7 @@ next:
     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
@@ -216,7 +222,7 @@ void AliMUONClusterFinderAZ::AddPad(Int_t cath, Int_t digit)
     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;
@@ -416,6 +422,7 @@ Bool_t AliMUONClusterFinderAZ::CheckPrecluster(Int_t *nShown)
        fXyq[2][i] = -2;
        fnPads[cath]--;
       }
+      if (fDraw) fDraw->UpdateCluster(npad);
     } // if (nFlags > 2)
 
     // Check correlations of cathode charges
@@ -432,39 +439,54 @@ Bool_t AliMUONClusterFinderAZ::CheckPrecluster(Int_t *nShown)
       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];
@@ -472,41 +494,39 @@ Bool_t AliMUONClusterFinderAZ::CheckPrecluster(Int_t *nShown)
            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;
@@ -711,35 +731,61 @@ void AliMUONClusterFinderAZ::AdjustPixel(Float_t wxmin, Float_t wymin)
 {
   // 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;
 }
 
 //_____________________________________________________________________________
@@ -805,9 +851,9 @@ Bool_t AliMUONClusterFinderAZ::MainLoop(Int_t iSimple)
     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] << " "; }
@@ -2048,21 +2094,22 @@ Int_t AliMUONClusterFinderAZ::Fit(Int_t nfit, Int_t *clustFit, TObjArray **clust
 
   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;
 }  
 
@@ -2573,12 +2620,16 @@ void AliMUONClusterFinderAZ::AddVirtualPad()
          //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);
@@ -2658,8 +2709,10 @@ void AliMUONClusterFinderAZ::PadsInXandY(Int_t &nInX, Int_t &nInY)
   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];
 }
 
 //_____________________________________________________________________________
index 7730d20..f90a5c7 100644 (file)
@@ -35,6 +35,11 @@ public:
   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);
@@ -47,19 +52,24 @@ protected:
 
   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
index 624b5b1..09005e0 100644 (file)
@@ -103,7 +103,7 @@ AliMUONClusterReconstructor::~AliMUONClusterReconstructor(void)
   return;
 }
 //____________________________________________________________________
-void AliMUONClusterReconstructor::Digits2Clusters()
+void AliMUONClusterReconstructor::Digits2Clusters(Int_t chBeg)
 {
 
     TClonesArray *dig1, *dig2, *digAll;
@@ -122,14 +122,18 @@ void AliMUONClusterReconstructor::Digits2Clusters()
     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;
@@ -193,5 +197,7 @@ void AliMUONClusterReconstructor::Digits2Clusters()
 void AliMUONClusterReconstructor::Trigger2Trigger() 
 {
 // copy trigger from TreeD to TreeR
+
+  fMUONData->SetTreeAddress("GLT");
   fMUONData->GetTriggerD();
 }
index 087c85f..9114e4d 100644 (file)
 /////////////////////////////////////
 
 #include <TObject.h>
-#include "AliMUONClusterFinderVS.h" //AZ
+#include "AliMUONClusterFinderVS.h" 
 
 class AliLoader;
 class AliMUON;
 class AliMUONRawCluster;
-//AZ class AliMUONClusterFinderVS;
 class AliMUONData;
 class AliRawReader;
 
@@ -32,17 +31,15 @@ class AliMUONClusterReconstructor : public TObject
 
  
   // 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:
@@ -54,7 +51,6 @@ class AliMUONClusterReconstructor : public TObject
 
   AliMUONData*            fMUONData;           //! Data container for MUON subsystem 
   AliMUONClusterFinderVS* fRecModel;           //! cluster recontruction model
-  //AliMUONClusterFinderAZ* fRecModel;           //! cluster recontruction model
 
   // alice loader
   AliLoader* fLoader;
index 9bd25ad..f1abcb7 100644 (file)
@@ -983,7 +983,7 @@ void AliMUONData::SetTreeAddress(Option_t* option)
     }
 
   }
-  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) {
diff --git a/MUON/AliMUONDetElement.cxx b/MUON/AliMUONDetElement.cxx
new file mode 100644 (file)
index 0000000..7133a9a
--- /dev/null
@@ -0,0 +1,290 @@
+/**************************************************************************
+ * 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();
+}
+*/
diff --git a/MUON/AliMUONDetElement.h b/MUON/AliMUONDetElement.h
new file mode 100644 (file)
index 0000000..38a6d8e
--- /dev/null
@@ -0,0 +1,77 @@
+#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
diff --git a/MUON/AliMUONEventRecoCombi.cxx b/MUON/AliMUONEventRecoCombi.cxx
new file mode 100644 (file)
index 0000000..553d6c5
--- /dev/null
@@ -0,0 +1,224 @@
+/**************************************************************************
+ * 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;
+}
diff --git a/MUON/AliMUONEventRecoCombi.h b/MUON/AliMUONEventRecoCombi.h
new file mode 100644 (file)
index 0000000..985cb47
--- /dev/null
@@ -0,0 +1,50 @@
+#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
index 5fcc7b5..1d5f15f 100644 (file)
 #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"
@@ -64,12 +65,13 @@ void AliMUONReconstructor::Reconstruct(AliRunLoader* runLoader) const
   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);
   }
@@ -78,32 +80,43 @@ void AliMUONReconstructor::Reconstruct(AliRunLoader* runLoader) const
   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();
@@ -120,8 +133,18 @@ void AliMUONReconstructor::Reconstruct(AliRunLoader* runLoader) const
     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();
@@ -129,7 +152,7 @@ void AliMUONReconstructor::Reconstruct(AliRunLoader* runLoader) const
 
     dataEvent->ResetRawClusters();
     dataEvent->ResetTrigger();
-    dataEvent->ResetRecTracks(); 
+    dataEvent->ResetRecTracks();  
     dataEvent->ResetRecTriggerTracks();
 
   }
@@ -164,9 +187,9 @@ void AliMUONReconstructor::Reconstruct(AliRunLoader* runLoader, AliRawReader* ra
   loader->LoadDigits("RECREATE");
 
 
-  //   Loop over events 
+  //   Loop over events  
   Int_t iEvent = 0;
-           
+            
   while (rawReader->NextEvent()) {
     printf("Event %d\n",iEvent);
     runLoader->GetEvent(iEvent++);
@@ -178,29 +201,29 @@ void AliMUONReconstructor::Reconstruct(AliRunLoader* runLoader, AliRawReader* ra
     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");
@@ -220,8 +243,8 @@ void AliMUONReconstructor::Reconstruct(AliRunLoader* runLoader, AliRawReader* ra
     recoEvent->EventReconstruct();
     dataEvent->Fill("RT");
 
-    loader->WriteTracks("OVERWRITE"); 
+    loader->WriteTracks("OVERWRITE");  
+  
     //--------------------------- Resetting branches -----------------------
     dataCluster->ResetDigits();
     dataCluster->ResetRawClusters();
@@ -231,7 +254,7 @@ void AliMUONReconstructor::Reconstruct(AliRunLoader* runLoader, AliRawReader* ra
     dataEvent->ResetTrigger();
     dataEvent->ResetRecTracks();
     dataEvent->ResetRecTriggerTracks();
+  
   }
   loader->UnloadRecPoints();
   loader->UnloadTracks();
@@ -246,12 +269,12 @@ void AliMUONReconstructor::FillESD(AliRunLoader* runLoader, AliESD* esd) const
 {
   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;
@@ -267,7 +290,7 @@ void AliMUONReconstructor::FillESD(AliRunLoader* runLoader, AliESD* esd) const
   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)
@@ -304,21 +327,21 @@ void AliMUONReconstructor::FillESD(AliRunLoader* runLoader, AliESD* esd) const
   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++) {
 
@@ -352,21 +375,21 @@ void AliMUONReconstructor::FillESD(AliRunLoader* runLoader, AliESD* esd) const
     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;
index cb7a413..a4f3979 100644 (file)
 #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"
 
@@ -40,15 +42,14 @@ const Int_t AliMUONTrackK::fgkTriesMax = 10000;
 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");
 
   //__________________________________________________________________________
@@ -59,7 +60,6 @@ AliMUONTrackK::AliMUONTrackK()
   // Default constructor
 
   fgTrackReconstructor = NULL; // pointer to event reconstructor
-  fgMUON = NULL; // pointer to Muon module
   fgHitForRec = NULL; // pointer to points
   fgNOfPoints = 0; // number of points
 
@@ -81,16 +81,16 @@ AliMUONTrackK::AliMUONTrackK()
 
   //__________________________________________________________________________
 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;
@@ -109,7 +109,6 @@ AliMUONTrackK::AliMUONTrackK(AliMUONTrackReconstructor *TrackReconstructor, TClo
 
   //__________________________________________________________________________
 AliMUONTrackK::AliMUONTrackK(AliMUONSegment *segment)
-  //AZ: TObject()
   : AliMUONTrack(segment, segment, fgTrackReconstructor)
 {
   // Constructor from a segment
@@ -186,7 +185,7 @@ AliMUONTrackK::AliMUONTrackK(AliMUONSegment *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;
@@ -284,7 +283,6 @@ AliMUONTrackK::~AliMUONTrackK()
 
   //__________________________________________________________________________
 AliMUONTrackK::AliMUONTrackK (const AliMUONTrackK& source)
-  //AZ: TObject(source)
   : AliMUONTrack(source)
 {
 // Protected copy constructor
@@ -373,7 +371,6 @@ void AliMUONTrackK::EvalCovariance(Double_t dZ)
   (*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,
@@ -381,7 +378,6 @@ void AliMUONTrackK::EvalCovariance(Double_t dZ)
   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:");
@@ -395,9 +391,9 @@ Bool_t AliMUONTrackK::KalmanFilter(Int_t ichamBeg, Int_t ichamEnd, Bool_t Back,
   // 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;
@@ -433,7 +429,12 @@ Bool_t AliMUONTrackK::KalmanFilter(Int_t ichamBeg, Int_t ichamEnd, Bool_t Back,
        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();
@@ -445,12 +446,17 @@ Bool_t AliMUONTrackK::KalmanFilter(Int_t ichamBeg, Int_t ichamEnd, Bool_t Back,
     // 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
@@ -471,8 +477,25 @@ Bool_t AliMUONTrackK::KalmanFilter(Int_t ichamBeg, Int_t ichamEnd, Bool_t Back,
          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
@@ -579,6 +602,7 @@ Bool_t AliMUONTrackK::KalmanFilter(Int_t ichamBeg, Int_t ichamEnd, Bool_t Back,
       Back = kFALSE;
       fRecover =0;
       ichambOK = ((AliMUONHitForRec*)((*fTrackHitsPtr)[fNTrackHits-1]))->GetChamberNumber();
+      //? if (fgTrackReconstructor->GetTrackMethod() != 3) currIndx = fgHitForRec->IndexOf(fSkipHit);
       currIndx = fgHitForRec->IndexOf(fSkipHit);
     }
 
@@ -595,7 +619,6 @@ Bool_t AliMUONTrackK::KalmanFilter(Int_t ichamBeg, Int_t ichamEnd, Bool_t Back,
        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
@@ -603,7 +626,7 @@ Bool_t AliMUONTrackK::KalmanFilter(Int_t ichamBeg, Int_t ichamEnd, Bool_t Back,
       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 {
@@ -684,19 +707,7 @@ void AliMUONTrackK::ParPropagation(Double_t zEnd)
   //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)
 {
@@ -717,7 +728,6 @@ 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();
@@ -755,7 +765,6 @@ void AliMUONTrackK::WeightPropagation(Double_t zEnd, Bool_t smooth)
 
   // 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);
@@ -786,7 +795,7 @@ void AliMUONTrackK::WeightPropagation(Double_t zEnd, Bool_t smooth)
 }
 
   //__________________________________________________________________________
-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
@@ -795,15 +804,35 @@ Bool_t AliMUONTrackK::FindPoint(Int_t ichamb, Double_t zEnd, Int_t currIndx, Int
   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 {
@@ -818,12 +847,20 @@ Bool_t AliMUONTrackK::FindPoint(Int_t ichamb, Double_t zEnd, Int_t currIndx, Int
   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) {
@@ -849,6 +886,13 @@ Bool_t AliMUONTrackK::FindPoint(Int_t ichamb, Double_t zEnd, Int_t currIndx, Int
        }
        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());
 
@@ -893,6 +937,7 @@ Bool_t AliMUONTrackK::FindPoint(Int_t ichamb, Double_t zEnd, Int_t currIndx, Int
          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) {
@@ -927,7 +972,6 @@ Bool_t AliMUONTrackK::FindPoint(Int_t ichamb, Double_t zEnd, Int_t currIndx, Int
            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--) {
@@ -967,7 +1011,6 @@ Bool_t AliMUONTrackK::FindPoint(Int_t ichamb, Double_t zEnd, Int_t currIndx, Int
            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; 
            }
@@ -984,7 +1027,6 @@ Bool_t AliMUONTrackK::FindPoint(Int_t ichamb, Double_t zEnd, Int_t currIndx, Int
     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) {
@@ -1024,9 +1066,8 @@ void AliMUONTrackK::TryPoint(TMatrixD &point, const TMatrixD &pointWeight, TMatr
   // 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:");
   }
@@ -1054,8 +1095,6 @@ void AliMUONTrackK::MSThin(Int_t sign)
   // 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 {
@@ -1072,8 +1111,6 @@ void AliMUONTrackK::MSThin(Int_t sign)
 
   (*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;
@@ -1303,8 +1340,8 @@ void AliMUONTrackK::Branson(void)
   // 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:");
   }
@@ -1504,15 +1541,27 @@ vertex:
     // 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", "   ");
     }
@@ -1532,7 +1581,6 @@ vertex:
   /* 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 {
@@ -1607,7 +1655,6 @@ void AliMUONTrackK::MSLine(Double_t dZ, Double_t x0)
   // Get weight matrix
   *fWeight = *fCovariance;
   if (fWeight->Determinant() != 0) {
-    //  fWeight->Invert();
     Int_t ifail;
     mnvertLocal(&((*fWeight)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
   } else {
@@ -2046,7 +2093,6 @@ void AliMUONTrackK::Outlier()
   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;
index c500d12..7c79a3d 100644 (file)
@@ -11,7 +11,7 @@
 /// \brief Kalman track in MUON arm of ALICE
 
 #include <TObject.h>
-#include "AliMUONTrack.h" //AZ
+#include "AliMUONTrack.h" 
 
 class TArrayD;
 class TMatrixD;
@@ -19,10 +19,9 @@ class AliMUONTrackReconstructor;
 class TClonesArray;
 class TObjArray;
 class AliMUONSegment;
-class AliMUON;
 class AliMUONHitForRec;
+class AliMUONEventRecoCombi;
 
-//AZ class AliMUONTrackK : public TObject {
 class AliMUONTrackK : public AliMUONTrack {
 
  public:
@@ -30,7 +29,6 @@ class AliMUONTrackK : public AliMUONTrack {
   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
 
@@ -82,9 +80,10 @@ class AliMUONTrackK : public AliMUONTrack {
  
   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
@@ -123,7 +122,7 @@ class AliMUONTrackK : public AliMUONTrack {
   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);
index 0baedd7..3c6c932 100644 (file)
@@ -13,6 +13,8 @@
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
+/* $Id$ */
+
 ////////////////////////////////////
 //
 // MUON track reconstructor in ALICE (class renamed from AliMUONEventReconstructor)
@@ -34,7 +36,7 @@
 #include <Riostream.h>
 #include <TDirectory.h>
 #include <TFile.h>
-#include <TMatrixD.h> //AZ
+#include <TMatrixD.h>
 
 #include "AliMUONTrackReconstructor.h"
 #include "AliMUON.h"
@@ -52,7 +54,7 @@
 #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"
@@ -102,7 +104,6 @@ AliMUONTrackReconstructor::AliMUONTrackReconstructor(AliLoader* loader)
   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 ????
 
@@ -477,7 +478,7 @@ void AliMUONTrackReconstructor::MakeEventToBeReconstructed(void)
     // 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
   }
@@ -711,12 +712,15 @@ void AliMUONTrackReconstructor::AddHitsForRecFromRawClusters(TTree* TR)
   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++) {
@@ -741,6 +745,8 @@ void AliMUONTrackReconstructor::AddHitsForRecFromRawClusters(TTree* TR)
       //  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);
@@ -754,7 +760,7 @@ void AliMUONTrackReconstructor::AddHitsForRecFromRawClusters(TTree* TR)
        hitForRec->Dump();}
     } // end of cluster loop
   } // end of chamber loop
-  SortHitsForRecWithIncreasingChamber(); //AZ 
+  SortHitsForRecWithIncreasingChamber(); 
   return;
 }
 
@@ -766,9 +772,8 @@ void AliMUONTrackReconstructor::MakeSegments(void)
   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;
@@ -903,7 +908,7 @@ void AliMUONTrackReconstructor::MakeTracks(void)
   // 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
@@ -1188,6 +1193,7 @@ void AliMUONTrackReconstructor::FollowTracks(void)
   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");
@@ -1611,8 +1617,7 @@ void AliMUONTrackReconstructor::MakeTrackCandidatesK(void)
     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;
@@ -1631,8 +1636,6 @@ void AliMUONTrackReconstructor::FollowTracksK(void)
   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();
 
@@ -1653,8 +1656,8 @@ void AliMUONTrackReconstructor::FollowTracksK(void)
        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)
     }
@@ -1689,8 +1692,7 @@ void AliMUONTrackReconstructor::FollowTracksK(void)
     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);
@@ -1699,7 +1701,6 @@ void AliMUONTrackReconstructor::FollowTracksK(void)
         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;
@@ -1718,7 +1719,6 @@ void AliMUONTrackReconstructor::FollowTracksK(void)
     }
 
     if (ok) {
-      //AZ trackK->SetTrackDir(-1);
       trackK->SetTrackDir(1);
       trackK->SetBPFlag(kFALSE);
       ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kFALSE,zDipole1,zDipole2);
@@ -1847,11 +1847,13 @@ void AliMUONTrackReconstructor::SetTrackMethod(Int_t iTrackMethod)
 {
   // 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;
 
 }
index ea99ff4..1eb8907 100644 (file)
@@ -21,6 +21,8 @@
 #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+;
index 40d9d6a..a438270 100644 (file)
@@ -23,7 +23,9 @@ SRCS:= AliMUONClusterReconstructor.cxx \
        AliMUONSubEventTracker.cxx \
        AliMUONRawData.cxx \
        AliMUONRawStream.cxx \
-       AliMUONClusterDrawAZ.cxx
+       AliMUONClusterDrawAZ.cxx \
+       AliMUONDetElement.cxx \
+       AliMUONEventRecoCombi.cxx
  
 
 HDRS:= $(SRCS:.cxx=.h)