MUON/AliMUONDetElement.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 // --------------------------
  19 // Class AliMUONDetElement
  20 // --------------------------
  21 // Detection element object containing information for combined
  22 // cluster / track finder in the MUON arm
  23 // Author: Alexander Zinchenko, JINR Dubna
  24
  25 #include "AliMUONDetElement.h"
  26 #include "AliMUON.h"
  27 #include "AliMUONSegmentation.h"
  28 #include "AliMUONDigit.h"
  29 #include "AliMUONHitMapA1.h"
  30 #include "AliMUONData.h"
  31 #include "AliMUONRawCluster.h"
  32 #include "AliMUONHitForRec.h"
  33 #include "AliMUONClusterInput.h"
  34 #include "AliMUONClusterFinderAZ.h"
  35 #include "AliMUONGeometryModuleTransformer.h"
  36 #include "AliMUONVGeometryDESegmentation.h"
  37 #include "AliMpVSegmentation.h"
  38
  39 #include "AliRun.h"
  40 #include "AliLog.h"
  41
  42 #include <TObjArray.h>
  43 #include <TClonesArray.h>
  44 #include <TVector2.h>
  45
  46 /// \cond CLASSIMP
  47 ClassImp(AliMUONDetElement) // Class implementation in ROOT context
  48 /// \endcond
  49   FILE *lun = 0x0; //fopen("hitmap.dat","w");
  50
  51 //_________________________________________________________________________
  52 AliMUONDetElement::AliMUONDetElement()
  53   : TObject(),
  54     fidDE(0),
  55     fIndex(0),
  56     fChamber(0),
  57     fZ(0.),
  58     fNHitsForRec(0),
  59     fRawClus(0x0),
  60     fHitsForRec(0x0),
  61     fRecModel(0x0)
  62 {
  63   /// Default constructor
  64   for (Int_t i = 0; i < 2; i++) {
  65     fHitMap[i] = NULL;
  66     fDigits[i] = NULL;
  67     fSeg[i] = NULL;
  68   }
  69 }
  70
  71 //_________________________________________________________________________
  72 AliMUONDetElement::AliMUONDetElement(Int_t idDE, AliMUONDigit *dig, AliMUONClusterFinderAZ *recModel)
  73   : TObject(),
  74     fidDE(idDE),
  75     fIndex(0),
  76     fChamber(idDE / 100 - 1),
  77     fZ(0.),
  78     fNHitsForRec(0),
  79     fRawClus(new TObjArray(10)),
  80     fHitsForRec(new TObjArray(10)),
  81     fRecModel(recModel)
  82 {
  83   /// Constructor
  84   fRawClus->SetOwner(kTRUE);
  85   fHitsForRec->SetOwner(kTRUE);
  86   fDigits[0] = new TObjArray(20);
  87   fDigits[1] = new TObjArray(20);
  88   fDigits[0]->SetOwner(kTRUE);
  89   fDigits[1]->SetOwner(kTRUE);
  90   AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON");
  91   AliMUONSegmentation *pSegmentation = pMUON->GetSegmentation();
  92   fSeg[0] = pSegmentation->GetModuleSegmentationByDEId(idDE, 0);
  93   fSeg[1] = pSegmentation->GetModuleSegmentationByDEId(idDE, 1);
  94   /*
  95   Float_t x, y, z;
  96   fSeg[dig->Cathode()]->GetPadC(fidDE, dig->PadX(), dig->PadY(), x, y, z);
  97   fZ = z;
  98   */
  99   AddDigit(dig);
 100 }
 101
 102 //_________________________________________________________________________
 103 AliMUONDetElement::~AliMUONDetElement()
 104 {
 105   /// Destructor
 106
 107   for (Int_t i = 0; i < 2; i++) {
 108     delete fHitMap[i];
 109     delete fDigits[i];
 110   }
 111   delete fRawClus;
 112   for (Int_t i = 0; i < fNHitsForRec; i++) {
 113     Int_t ntracks = ((AliMUONHitForRec*)fHitsForRec->UncheckedAt(i))->GetNTrackHits();
 114     //cout << ntracks << endl;
 115     if (ntracks) fHitsForRec->RemoveAt(i);
 116   }
 117   delete fHitsForRec;
 118   //fHitsForRec->SetOwner(kFALSE);
 119   //fHitsForRec->Clear(); cout << " here " << fHitsForRec->GetEntriesFast() << endl; delete fHitsForRec; fHitsForRec = 0;
 120 }
 121
 122 //_________________________________________________________________________
 123 Int_t AliMUONDetElement::Compare(const TObject* detElem) const
 124 {
 125   /// "Compare" function to sort in Z (towards interaction point)
 126   /// Returns -1 (0, +1) if charge of current pixel
 127   /// is greater than (equal to, less than) charge of pixel
 128   if (fZ > ((AliMUONDetElement*)detElem)->Z()) return(+1);
 129   else if (fZ == ((AliMUONDetElement*)detElem)->Z()) return( 0);
 130   else return(-1);
 131 }
 132
 133 //_________________________________________________________________________
 134 void AliMUONDetElement::Fill(AliMUONData */*data*/)
 135 {
 136   /// Fill hit maps
 137   fLeft[0] = fDigits[0]->GetEntriesFast();
 138   fLeft[1] = fDigits[1]->GetEntriesFast();
 139
 140   Int_t  npx0  = fSeg[0]->Npx(fidDE)+1;
 141   Int_t  npy0  = fSeg[0]->Npy(fidDE)+1;
 142   fHitMap[0] = new AliMUONHitMapA1(npx0, npy0, fDigits[0]);
 143   fHitMap[0]->FillHits();
 144
 145   Int_t  npx1  = fSeg[1]->Npx(fidDE)+1;
 146   Int_t  npy1  = fSeg[1]->Npy(fidDE)+1;
 147   fHitMap[1] = new AliMUONHitMapA1(npx1, npy1, fDigits[1]);
 148   fHitMap[1]->FillHits();
 149
 150   // The part below is just for debugging (fill rec. points already found)
 151   /*
 152   fLeft[0] = fLeft[1] = 0;
 153   TClonesArray *rawClus = data->RawClusters(fChamber);
 154   cout << rawClus << " " << rawClus->GetEntriesFast() << endl;
 155   for (Int_t i = 0; i < rawClus->GetEntriesFast(); i++) {
 156     AliMUONRawCluster *recP = (AliMUONRawCluster*) rawClus->UncheckedAt(i);
 157     cout << fChamber << " " << recP->GetZ(0) << " " << recP->GetZ(1) << " " << fZ << endl;
 158     if (TMath::Abs(recP->GetZ(0)-fZ) > 0.5) continue;
 159     if (!Inside(recP->GetX(0), recP->GetY(0), recP->GetZ(0))) continue;
 160     AddHitForRec(recP); // add hit for rec.
 161     rawClus->RemoveAt(i); // remove
 162   }
 163   cout << fHitsForRec->GetEntriesFast() << endl;
 164   rawClus->Compress();
 165   */
 166 }
 167
 168 //_________________________________________________________________________
 169 void AliMUONDetElement::AddDigit(AliMUONDigit *dig)
 170 {
 171   /// Add digit
 172
 173   //fDigits[dig->Cathode()]->Add(dig);
 174   fDigits[dig->Cathode()]->Add(new AliMUONDigit(*dig));
 175   Float_t x, y, z;
 176   fSeg[dig->Cathode()]->GetPadC(fidDE, dig->PadX(), dig->PadY(), x, y, z);
 177   fZ += z;
 178   /*
 179   Int_t ndig = fDigits[dig->Cathode()]->GetEntriesFast();
 180   new((*fDigits[dig->Cathode()])[ndig]) AliMUONDigit(*dig);
 181   */
 182 }
 183
 184 //_________________________________________________________________________
 185 Bool_t AliMUONDetElement::Inside(Double_t x, Double_t y, Double_t z, Double_t dx, Double_t dy) const
 186 {
 187   /// Check if point is inside detection element
 188
 189   dx = TMath::Max (dx, 1.);
 190   dy = TMath::Max (dy, 1.);
 191   Int_t ix, iy;
 192   ix = iy = 0;
 193   Double_t xl, yl, zl;
 194   for (Int_t i = 0; i < 2; i++) {
 195     if (!fSeg[i]) continue;
 196     if (i == 0 || fSeg[0] == 0x0) {
 197       fSeg[i]->GetTransformer()->Global2Local(fidDE, x, y, z, xl, yl, zl);
 198       //cout << xl << " " << yl << " " << zl << endl;
 199       if (TMath::Abs (zl) > 0.5) return kFALSE;
 200       //if (fChamber < 4 && (xl < 0 || yl < 0)) return kFALSE;
 201     }
 202     fSeg[i]->GetPadI(fidDE, x, y, z, ix, iy);
 203     //cout << x << " " << y << " " << z << " " << fChamber << " " << ix << " " << iy << " " << fSeg[i]->Npx(fidDE) << " " << fSeg[i]->Npy(fidDE) /*<< " " << fSeg[i]->GetPadI(fidDE, x, y, z, ix, iy)*/ << " " << fSeg[i]->HasPad(fidDE, (Float_t)x, (Float_t)y, (Float_t)z) << endl;
 204     //if (ix > 0 && iy > 0 && ix <= fSeg[i]->Npx(fidDE) && iy <= fSeg[i]->Npy(fidDE)) return kTRUE;
 205     if (fSeg[i]->HasPad(fidDE, (Float_t)x, (Float_t)y, (Float_t)z)) return kTRUE;
 206   }
 207   //Double_t xl, yl, zl;
 208   //fSeg[0]->GetTransformer()->Global2Local(fidDE, x, y, z, xl, yl, zl);
 209   //cout << xl << " " << yl << " " << zl << endl;
 210
 211   Int_t cath = fSeg[0] ? 0 : 1;
 212   const AliMpVSegmentation *mpSeg = fSeg[cath]->GetDESegmentation(fidDE)->GetMpSegmentation();
 213   TVector2 dim = mpSeg->Dimensions();
 214   if (fChamber < 4) dim *= 2.; // quadrants
 215   //cout << fChamber << " " << dim.X() << " " << dim.Y() << endl;
 216   if (TMath::Abs (xl) - dim.X() > dx) return kFALSE;
 217   if (TMath::Abs (yl) - dim.Y() > dy) return kFALSE;
 218
 219   Int_t flagX = 0, flagY = 0;
 220   if (TMath::Abs (xl) - dim.X() < -2) flagX = 1;
 221   if (TMath::Abs (yl) - dim.Y() < -2) flagY = 1;
 222   if (flagX && flagY) flagX = flagY = 0; // both coordinates "inside" det. elem. limits
 223
 224   // Check for edge effect (extrapolated track "right outside" det. elem. boundaries
 225   // (+- dx * cath in X and Y)
 226   dx /= 2;
 227   dy /= 2;
 228   for (Int_t i = 0; i < 2; i++) {
 229     if (!fSeg[i]) continue;
 230     for (Int_t idx = -1; idx < 2; idx++) {
 231       if (flagX && idx) continue;
 232       Double_t x1 = x + dx * idx * (i + 1);
 233       for (Int_t idy = -1; idy < 2; idy++) {
 234         if (idx == 0 && idy == 0) continue;
 235         if (flagY && idy) continue;
 236         Double_t y1 = y + dy * idy * (i + 1);
 237         fSeg[i]->GetPadI(fidDE, x1, y1, z, ix, iy);
 238         //cout << x1 << " " << y1 << " " << z << " " << fChamber << " " << ix << " " << iy << " " << fSeg[i]->Npx(fidDE) << " " << fSeg[i]->Npy(fidDE) /*<< " " << fSeg[i]->GetPadI(fidDE, x, y, z, ix, iy)*/ << " " << fSeg[i]->HasPad(fidDE, (Float_t)x1, (Float_t)y1, (Float_t)z) << endl;
 239         //if (ix > 0 && iy > 0 && ix <= fSeg[i]->Npx(fidDE) && iy <= fSeg[i]->Npy(fidDE)) return kTRUE;
 240         if (fSeg[i]->HasPad(fidDE, (Float_t)x1, (Float_t)y1, (Float_t)z)) return kTRUE;
 241       }
 242     }
 243   }
 244   return kFALSE;
 245 }
 246
 247 //_________________________________________________________________________
 248 void AliMUONDetElement::ClusterReco(Double_t xTrack, Double_t yTrack, Double_t dx, Double_t dy)
 249 {
 250   /// Run cluster reconstruction around point (x,y) inside window (dx,dy)
 251
 252   if (fLeft[0] == 0 && fLeft[1] == 0) return; // all digits have been used
 253   AliMUONClusterInput::Instance()->SetDigits(fChamber, fidDE,
 254              (TClonesArray*)fDigits[0], (TClonesArray*)fDigits[1]);
 255
 256   // Mark used pads
 257   for (Int_t cath = 0; cath < 2; cath++) {
 258     Int_t ndig = fDigits[cath]->GetEntriesFast();
 259     if (ndig == 0) continue; // empty cathode
 260     //cout << fidDE << " " << cath << " " << ndig << endl;
 261     for (Int_t i = 0; i < ndig; i++) {
 262       if (fLeft[cath] == 0) { fRecModel->SetUsed(cath,i); continue; }
 263       AliMUONDigit *dig = (AliMUONDigit*) fDigits[cath]->UncheckedAt(i);
 264       //cout << i << " " << dig << " " /*<< dig->PadX() << " " << dig->PadY() << " " << fHitMap[cath]->TestHit(dig->PadX(), dig->PadY()) */<< endl;
 265       if (fHitMap[cath]->TestHit(dig->PadX(), dig->PadY()) == kUsed) fRecModel->SetUsed(cath,i);
 266       else fRecModel->SetUnused(cath,i);
 267     }
 268   }
 269
 270   fRecModel->ResetRawClusters();
 271
 272   for (Int_t cath = 0; cath < 2; cath++) {
 273     Int_t ndig = fDigits[cath]->GetEntriesFast();
 274     if (ndig == 0) continue; // empty cathode
 275     // Loop over pads
 276     for (fSeg[cath]->FirstPad(fidDE, xTrack, yTrack, fZ, dx, dy);
 277          fSeg[cath]->MorePads(fidDE);
 278          fSeg[cath]->NextPad(fidDE)) {
 279       if (fLeft[cath] == 0) break;
 280       //cout << cath << " " << fSeg[cath]->Ix() << " " << fSeg[cath]->Iy() << " " << fSeg[cath]->DetElemId() << " " << fHitMap[cath]->TestHit(fSeg[cath]->Ix(), fSeg[cath]->Iy()) << endl;
 281       if (fHitMap[cath]->TestHit(fSeg[cath]->Ix(), fSeg[cath]->Iy()) == kEmpty ||
 282           fHitMap[cath]->TestHit(fSeg[cath]->Ix(), fSeg[cath]->Iy()) == kUsed) continue;
 283
 284       Int_t iOK = 0;
 285       // Set starting pad
 286       for (Int_t j = 0; j < ndig; j++) {
 287         AliMUONDigit *dig = (AliMUONDigit*) fDigits[cath]->UncheckedAt(j);
 288         if (dig->PadX() != fSeg[cath]->Ix() || dig->PadY() != fSeg[cath]->Iy()) continue;
 289         //cout << fidDE << " " << j << " " << fDigits[cath]->GetEntriesFast() << " " << fSeg[cath]->Ix() << " " << fSeg[cath]->Iy() << endl;
 290         fRecModel->SetStart(cath, j);
 291         iOK = 1;
 292         break;
 293       }
 294       if (!iOK) continue;
 295
 296       fRecModel->FindRawClusters();
 297       Int_t nClusEnd = fRecModel->GetRawClusters()->GetEntriesFast();
 298       //cout << " ***nclus: " << nClusEnd << endl;
 299       for (Int_t i = 0; i < nClusEnd; i++) {
 300         AliMUONRawCluster *clus = (AliMUONRawCluster*) fRecModel->GetRawClusters()->UncheckedAt(i);
 301         AddHitForRec(clus); // add hit for rec.
 302         //cout << clus->GetX(0) << " " << clus->GetY(0) << endl;
 303       }
 304       // Mark used pads
 305       for (Int_t cath1 = 0; cath1 < 2; cath1++) {
 306         Int_t ndig1 = fDigits[cath1]->GetEntriesFast();
 307         for (Int_t j = 0; j < ndig1; j++) {
 308           if (fLeft[cath1] == 0) break;
 309           AliMUONDigit *dig = (AliMUONDigit*) fDigits[cath1]->UncheckedAt(j);
 310           Float_t x, y, z;
 311           fSeg[cath1]->GetPadC(fidDE,dig->PadX(),dig->PadY(),x,y,z);
 312           //cout << "clus " << cath1 << " " << fLeft[cath1] << " " << dig->PadX() << " " << dig->PadY() << " " << x << " " << y << " " << z << " " << fRecModel->GetUsed(cath1,j) << endl;
 313           if (!fRecModel->GetUsed(cath1,j)) continue;
 314           if (fHitMap[cath1]->TestHit(dig->PadX(), dig->PadY()) == kUsed) continue;
 315           fHitMap[cath1]->FlagHit(dig->PadX(), dig->PadY());
 316           if (lun) fprintf(lun," %d %d %d %d \n", cath1, fidDE, dig->PadX(), dig->PadY());
 317           fLeft[cath1]--;
 318         }
 319       }
 320     } // for (fSeg[cath]->FirstPad(...
 321   } // for (Int_t cath = 0;
 322 }
 323
 324 //_________________________________________________________________________
 325 void AliMUONDetElement::AddHitForRec(AliMUONRawCluster *clus)
 326 {
 327   /// Make HitForRec from raw cluster (rec. point)
 328
 329   fRawClus->Add(new AliMUONRawCluster(*clus));
 330   //AliMUONHitForRec *hitForRec =
 331   //new ((*fHitsForRec)[fNHitsForRec++]) AliMUONHitForRec(clus);
 332   AliMUONHitForRec *hitForRec = new AliMUONHitForRec(clus);
 333   fHitsForRec->Add(hitForRec);
 334   fNHitsForRec++;
 335
 336   // more information into HitForRec
 337   //  resolution: info should be already in raw cluster and taken from it ????
 338   //hitForRec->SetBendingReso2(-1); //fBendingResolution * fBendingResolution);
 339   //hitForRec->SetNonBendingReso2(-1); //fNonBendingResolution * fNonBendingResolution);
 340   hitForRec->SetBendingReso2(clus->GetErrY() * clus->GetErrY());
 341   hitForRec->SetNonBendingReso2(clus->GetErrX() * clus->GetErrX());
 342   //  original raw cluster
 343   hitForRec->SetChamberNumber(fChamber);
 344   hitForRec->SetZ(clus->GetZ(0));
 345   //hitForRec->SetHitNumber(-(fIndex+1)*100000-fNHitsForRec+1);
 346   hitForRec->SetHitNumber(-(fIndex+1)*100000-fRawClus->GetEntriesFast()+1);
 347   //delete clus; // for now
 348 }
 349
 350 /*
 351 //_________________________________________________________________________
 352 Int_t AliMUONDetElement::GetMapElem(AliMUONDigit *digit)
 353 {
 354   Int_t cath = digit->Cathode();
 355   return 0;
 356
 357 }
 358
 359 //_________________________________________________________________________
 360 void AliMUONDetElement::SetMapElem(const AliMUONDigit *digit, Int_t flag)
 361 {
 362   Int_t cath = digit->Cathode();
 363 }
 364 */