ITS/UPGRADE/AliITSUClusterizer.cxx

   1 #include <TTree.h>
   2 #include <TFile.h>
   3 #include <TObjArray.h>
   4 #include <TParticle.h>
   5 #include <TMath.h>
   6 #include "AliRun.h"
   7 #include "AliLog.h"
   8 #include "AliMC.h"
   9 #include <AliITSUSegmentationPix.h>
  10 #include "AliITSUClusterizer.h"
  11 #include "AliITSUGeomTGeo.h"
  12 #include "AliITSUSegmentationPix.h"
  13 #include "AliITSdigit.h"
  14 #include "AliITSURecoParam.h"
  15 #include "AliITSUAux.h"
  16 using namespace TMath;
  17 using namespace AliITSUAux;
  18
  19 ClassImp(AliITSUClusterizer)
  20
  21 //______________________________________________________________________________
  22 AliITSUClusterizer::AliITSUClusterizer(Int_t initNRow)
  23 :  fVolID(-1)
  24   ,fAllowDiagonalClusterization(kFALSE)
  25   ,fSegm(0)
  26   ,fRecoParam(0)
  27   ,fInputDigits(0)
  28   ,fInputDigitsReadIndex(0)
  29   ,fLayerID(0)
  30   ,fNLabels(0)
  31   ,fRawData(kFALSE)
  32   ,fLorAngCorrection(0)
  33   ,fOutputClusters(0)
  34   ,fDigitFreelist(0)
  35   ,fPartFreelist(0)
  36   ,fCandFreelist(0)
  37   ,fDigitFreelistBptrFirst(0)
  38   ,fDigitFreelistBptrLast(0)
  39   ,fPartFreelistBptr(0)
  40   ,fCandFreelistBptr(0)
  41 #ifdef _ClusterTopology_
  42   ,fTopology(0)
  43   ,fMinCol(0)
  44   ,fMinRow(0)
  45 #endif //_ClusterTopology_
  46 {
  47   SetUniqueID(0);
  48   // c-tor
  49   SetNRow(initNRow);
  50 //
  51 #ifdef _ClusterTopology_
  52   AliInfo("*********************************************************************");
  53   AliInfo("ATTENTION: YOU ARE RUNNING IN SPECIAL MODE OF STORING CLUSTER PATTERN");
  54   AliInfo("*********************************************************************");
  55 #endif //_ClusterTopology_
  56
  57 //
  58 }
  59
  60 //______________________________________________________________________________
  61 void AliITSUClusterizer::SetSegmentation(const AliITSUSegmentationPix *segm)
  62 {
  63   // attach segmentation, if needed, reinitialize array
  64   fSegm = segm;
  65   SetNRow(fSegm->GetNRow()); // reinitialize if needed
  66
  67 }
  68
  69 //______________________________________________________________________________
  70 void AliITSUClusterizer::SetNRow(Int_t nr)
  71 {
  72   // update buffers
  73   int nrOld = GetUniqueID();
  74   if (nrOld>=nr) return;
  75   SetUniqueID(nr);
  76   while (fDigitFreelistBptrFirst) {
  77     AliITSUClusterizerClusterDigit *next = fDigitFreelistBptrFirst[kDigitChunkSize-1].fNext;
  78     delete[] fDigitFreelistBptrFirst;
  79     fDigitFreelistBptrFirst=next;
  80   }
  81   delete[] fPartFreelistBptr;
  82   delete[] fCandFreelistBptr;
  83   //
  84   fPartFreelist=fPartFreelistBptr = new AliITSUClusterizerClusterPart[nr+1];
  85   fCandFreelist=fCandFreelistBptr = new AliITSUClusterizerClusterCand[nr+1];
  86   for (int i=0;i<nr;++i) {
  87     fPartFreelistBptr[i].fNextInRow = &fPartFreelistBptr[i+1];
  88     fCandFreelistBptr[i].fNext      = &fCandFreelistBptr[i+1];
  89   }
  90 }
  91
  92 //______________________________________________________________________________
  93 AliITSUClusterizer::~AliITSUClusterizer()
  94 {
  95   // d-tor
  96   while (fDigitFreelistBptrFirst) {
  97     AliITSUClusterizerClusterDigit *next = fDigitFreelistBptrFirst[kDigitChunkSize-1].fNext;
  98     delete[] fDigitFreelistBptrFirst;
  99     fDigitFreelistBptrFirst=next;
 100   }
 101   delete[] fPartFreelistBptr;
 102   delete[] fCandFreelistBptr;
 103   //
 104 }
 105
 106 //______________________________________________________________________________
 107 AliITSUClusterizer::AliITSUClusterizerClusterDigit* AliITSUClusterizer::AllocDigitFreelist()
 108 {
 109   // allocate aux space
 110   AliITSUClusterizerClusterDigit *tmp = new AliITSUClusterizerClusterDigit[kDigitChunkSize];
 111   for (int i=0;i<kDigitChunkSize-2;++i) tmp[i].fNext=&tmp[i+1];
 112   tmp[kDigitChunkSize-2].fNext=0;
 113   tmp[kDigitChunkSize-1].fNext=0;
 114   if (!fDigitFreelistBptrFirst) fDigitFreelistBptrFirst=tmp;
 115   else                          fDigitFreelistBptrLast[kDigitChunkSize-1].fNext=tmp;
 116   fDigitFreelistBptrLast=tmp;
 117   return tmp;
 118 }
 119
 120 //______________________________________________________________________________
 121 AliITSUClusterizer::AliITSUClusterizerClusterDigit* AliITSUClusterizer::NextDigit()
 122 {
 123   // get next digit
 124   if (fInputDigitsReadIndex<fInputDigits->GetEntriesFast()) {
 125     AliITSdigit *tmp=static_cast<AliITSdigit*>(fInputDigits->UncheckedAt(fInputDigitsReadIndex++));
 126     AliITSUClusterizerClusterDigit *digit=AllocDigit();
 127     digit->fDigit=tmp;
 128     return digit;
 129   }
 130   else
 131     return 0;
 132 }
 133
 134 //______________________________________________________________________________
 135 void AliITSUClusterizer::AttachPartToCand(AliITSUClusterizerClusterCand *cand,AliITSUClusterizerClusterPart *part)
 136 {
 137   // attach part
 138   part->fParent = cand;
 139   part->fPrevInCluster = 0;
 140   part->fNextInCluster = cand->fFirstPart;
 141   if (cand->fFirstPart) cand->fFirstPart->fPrevInCluster = part;
 142   cand->fFirstPart=part;
 143 }
 144
 145 //______________________________________________________________________________
 146 void AliITSUClusterizer::MergeCands(AliITSUClusterizerClusterCand *a,AliITSUClusterizerClusterCand *b)
 147 {
 148   // merge cluster parts
 149   AliITSUClusterizerClusterPart *ipart=b->fFirstPart;
 150   AliITSUClusterizerClusterPart *jpart;
 151   do {
 152     jpart=ipart;
 153     jpart->fParent=a;
 154   } while ((ipart=ipart->fNextInCluster));
 155   jpart->fNextInCluster=a->fFirstPart;
 156   jpart->fNextInCluster->fPrevInCluster=jpart;
 157   a->fFirstPart=b->fFirstPart;
 158   // merge digits
 159   b->fLastDigit->fNext=a->fFirstDigit;
 160   a->fFirstDigit=b->fFirstDigit;
 161   //  DeallocCand(b);
 162 }
 163
 164 //______________________________________________________________________________
 165 void AliITSUClusterizer::Transform(AliITSUClusterPix *cluster,AliITSUClusterizerClusterCand *cand)
 166 {
 167   // convert set of digits to cluster data in LOCAL frame
 168   const double k1to12 = 1./12;
 169   static int maxLbinDigit = AliITSdigit::GetNTracks();
 170   //
 171   fNLabels = 0;
 172   Int_t n=0;
 173   cand->fLastDigit->fNext=0;
 174   double x=0,z=0,xmn=1e9,xmx=-1e9,zmn=1e9,zmx=-1e9,px=0,pz=0;
 175   float  cx,cz;
 176   int charge=0;
 177   //
 178   for (AliITSUClusterizerClusterDigit *idigit=cand->fFirstDigit;idigit;idigit=idigit->fNext) {
 179     AliITSdigit* dig = idigit->fDigit;
 180     fSegm->DetToLocal(dig->GetCoord2(),dig->GetCoord1(),cx,cz); // center of pixel
 181     //
 182     // account for possible diod shift
 183     double ddx,ddz, dx=fSegm->Dpx(dig->GetCoord2()), dz=fSegm->Dpz(dig->GetCoord1());
 184     fSegm->GetDiodShift(dig->GetCoord2(),dig->GetCoord1(),ddx,ddz);
 185     //
 186     charge += dig->GetSignal();
 187     x += cx;
 188     z += cz;
 189     if (cx<xmn) xmn=cx;
 190     if (cx>xmx) xmx=cx;
 191     if (cz<zmn) zmn=cz;
 192     if (cz>zmx) zmx=cz;
 193     x += ddx*dx;
 194     z += ddz*dz;
 195     px += dx;
 196     pz += dz;
 197     //
 198     if (!fRawData) {
 199       for(Int_t dlab=0;dlab<maxLbinDigit;dlab++){
 200         Int_t digitlab = (dig->GetTracks())[dlab];
 201         if(digitlab<0) continue;
 202         AddLabel(digitlab);
 203       }
 204     }
 205     //
 206     ++n;
 207   }
 208   Int_t nx=1,nz=1;
 209   double dx = xmx-xmn, dz = zmx-zmn;
 210   if (n>1) {
 211     double fac=1./n;
 212     x  *= fac;  // mean coordinates
 213     z  *= fac;
 214     px *= fac;  // mean pitch
 215     pz *= fac;
 216     nx = 1+Nint(dx/px);
 217     nz = 1+Nint(dz/pz);
 218   }
 219   x -= fLorAngCorrection;  // LorentzAngle correction
 220   cluster->SetX(x);
 221   cluster->SetZ(z);
 222   cluster->SetY(0);
 223   cluster->SetSigmaZ2(nz>1 ? dz*dz*k1to12 : pz*pz*k1to12);
 224   cluster->SetSigmaY2(nx>1 ? dx*dx*k1to12 : px*px*k1to12);
 225   cluster->SetSigmaYZ(0);
 226   cluster->SetFrameLoc();
 227   if (n>AliITSUClusterPix::kMaskNPix) n = AliITSUClusterPix::kMaskNPix;
 228   cluster->SetNxNzN(nx>AliITSUClusterPix::kMaskNX ? AliITSUClusterPix::kMaskNX : nx,
 229                     nz>AliITSUClusterPix::kMaskNZ ? AliITSUClusterPix::kMaskNZ : nz,
 230                     n>AliITSUClusterPix::kMaskNPix ? AliITSUClusterPix::kMaskNPix : n);
 231   cluster->SetQ(charge); // note: this is MC info
 232   //  if (cluster->GetNPix()!=n || cluster->GetNx()!=nx || cluster->GetNz()!=nz) {
 233   //  printf("PROBLEM: %d %d %d -> %d %d %d\n",nx,nz,n,
 234   //       cluster->GetNx(),cluster->GetNz(),cluster->GetNPix());
 235   //  cluster->Dump();
 236   //  }
 237   //
 238   if (!fRawData) {
 239     CheckLabels();
 240     int nl = Min(kMaxLabInCluster,fNLabels);
 241     for (int i=nl;i--;) cluster->SetLabel(fCurrLabels[i],i);
 242   }
 243   //
 244   // Set Volume id
 245   cluster->SetVolumeId(fVolID);
 246   //    printf("mod %d: (%.4lf,%.4lf)cm\n",fVolID,x,z);
 247   //
 248 #ifdef _ClusterTopology_
 249   FillClusterTopology(cand,cluster);
 250 #endif //_ClusterTopology_
 251   //
 252 }
 253
 254 //______________________________________________________________________________
 255 void AliITSUClusterizer::CloseCand(AliITSUClusterizerClusterCand *cand)
 256 {
 257   // finish cluster
 258   AliITSUClusterPix *cluster = (AliITSUClusterPix*)NextCluster();
 259   Transform(cluster,cand);
 260   DeallocDigits(cand->fFirstDigit,cand->fLastDigit);
 261   DeallocCand(cand);
 262 }
 263
 264 //______________________________________________________________________________
 265 void AliITSUClusterizer::ClosePart(AliITSUClusterizerClusterPart *part)
 266 {
 267   // finish cluster part
 268   AliITSUClusterizerClusterCand *cand=part->fParent;
 269   DetachPartFromCand(cand,part);
 270   DeallocPart(part);
 271   if (!cand->fFirstPart) CloseCand(cand);
 272 }
 273
 274 //______________________________________________________________________________
 275 void AliITSUClusterizer::CloseRemainingParts(AliITSUClusterizerClusterPart *part)
 276 {
 277   // finish what is left
 278   while (part) {
 279     AliITSUClusterizerClusterPart *next=part->fNextInRow;
 280     ClosePart(part);
 281     part=next;
 282   }
 283 }
 284
 285 //______________________________________________________________________________
 286 void AliITSUClusterizer::Clusterize()
 287 {
 288   // main algo for MC clustererization
 289   SetRawData(kFALSE);
 290   //
 291   AliITSUClusterizerClusterDigit *iDigit=NextDigit();
 292   AliITSUClusterizerClusterPart *iPrevRowBegin=0;
 293   AliITSUClusterizerClusterPart *iNextRowBegin=0;
 294   AliITSUClusterizerClusterPart *iPrevRow=0;
 295   AliITSUClusterizerClusterPart *iNextRow=0;
 296   Int_t lastV=0;
 297   //
 298   while (iDigit) {
 299     if (iDigit->fDigit->GetCoord2()!=lastV) {
 300       // NEW ROW
 301       if (iNextRow) iNextRow->fNextInRow=0;
 302       if (iPrevRowBegin) CloseRemainingParts(iPrevRowBegin);
 303       if (iDigit->fDigit->GetCoord2()==lastV+1) {
 304         iPrevRowBegin=iNextRowBegin;
 305         iPrevRow     =iNextRowBegin;
 306       }
 307       else {
 308         // there was an empty row
 309         CloseRemainingParts(iNextRowBegin);
 310         iPrevRowBegin=0;
 311         iPrevRow     =0;
 312       }
 313       iNextRowBegin=0;
 314       iNextRow     =0;
 315       lastV=iDigit->fDigit->GetCoord2();
 316     }
 317     // skip cluster parts before this digit
 318     int limCol = iDigit->fDigit->GetCoord1()-fAllowDiagonalClusterization;
 319     while (iPrevRow && iPrevRow->fUEnd < limCol) {
 320       iPrevRow=iPrevRow->fNextInRow;
 321     }
 322     // find the longest continous line of digits [iDigit,pDigit]=[iDigit,jDigit)
 323     AliITSUClusterizerClusterCand *cand=AllocCand();
 324     AliITSUClusterizerClusterDigit *pDigit=iDigit;
 325     AliITSUClusterizerClusterDigit *jDigit=NextDigit();
 326     cand->fFirstPart=0;
 327     cand->fFirstDigit=cand->fLastDigit=iDigit; // NB: first diggit is attached differently
 328     iDigit->fNext=0;
 329     Int_t lastU =iDigit->fDigit->GetCoord1();
 330     Int_t lastU1=lastU+1;
 331     while (jDigit && jDigit->fDigit->GetCoord1()==lastU1 && jDigit->fDigit->GetCoord2()==lastV) {
 332       pDigit=jDigit;
 333       jDigit=NextDigit();
 334       AttachDigitToCand(cand,pDigit);
 335       ++lastU1;
 336     }
 337     if (!fAllowDiagonalClusterization) --lastU1;
 338     AliITSUClusterizerClusterPart *part=AllocPart();
 339     part->fUBegin=lastU ;
 340     part->fUEnd  =lastU1;
 341     AttachPartToCand(cand,part);
 342     // merge all cluster candidates of the previous line touching this one,
 343     // advance to the last one, but keep that one the next active one
 344     AliITSUClusterizerClusterPart *jPrevRow=iPrevRow;
 345     while (jPrevRow && jPrevRow->fUBegin<=lastU1) {
 346       if (jPrevRow->fParent!=cand) {
 347         MergeCands(jPrevRow->fParent,cand);
 348         DeallocCand(cand);
 349         cand=jPrevRow->fParent;
 350       }
 351       iPrevRow=jPrevRow;
 352       jPrevRow=jPrevRow->fNextInRow;
 353     }
 354     if (iNextRow)
 355       iNextRow->fNextInRow=part;
 356     else
 357       iNextRowBegin=part;
 358     iNextRow=part;
 359     iDigit=jDigit;
 360   }
 361   // remove remaining cluster parts
 362   CloseRemainingParts(iPrevRowBegin);
 363   if (iNextRow) iNextRow->fNextInRow=0;
 364   CloseRemainingParts(iNextRowBegin);
 365   return;
 366 }
 367
 368 //______________________________________________________________________________
 369 void AliITSUClusterizer::PrepareLorentzAngleCorrection(Double_t bz)
 370 {
 371   // calculate parameters for Lorentz Angle correction. Must be called
 372   // after setting segmentation and recoparams
 373   fLorAngCorrection = 0.5*fRecoParam->GetTanLorentzAngle(fLayerID)*bz/kNominalBz*fSegm->Dy();
 374 }
 375
 376 //______________________________________________________________________
 377 void AliITSUClusterizer::CheckLabels()
 378 {
 379   // Tries to find mother's labels
 380   //
 381   if (fNLabels<1) return;
 382   AliRunLoader *rl = AliRunLoader::Instance();
 383   if(!rl) return;
 384   TTree *trK=(TTree*)rl->TreeK();
 385   if (!trK) return;
 386   //
 387   static int   labS[kMaxLabels];
 388   static float kine[kMaxLabels];
 389   Int_t nlabels = fNLabels;
 390   Int_t ntracks = gAlice->GetMCApp()->GetNtrack();
 391   for (Int_t i=fNLabels;i--;) labS[i] = fCurrLabels[i];
 392   //
 393   for (Int_t i=0;i<nlabels;i++) {
 394     Int_t label = labS[i];
 395     if (label>=ntracks) continue;
 396     TParticle *part=(TParticle*)gAlice->GetMCApp()->Particle(label);
 397     kine[i] = part->Energy() - part->GetCalcMass(); // kinetic energy
 398     if (kine[i] < 0.001) {    // reduce soft particles from the same cluster
 399       Int_t m=part->GetFirstMother();
 400       if (m<0) continue; // primary
 401       //
 402       if (part->GetStatusCode()>0) continue;
 403       //
 404       // if the parent is within the same cluster, assign parent's label
 405       for (int j=0;j<nlabels;j++) if (labS[j]==m) { labS[i] = m; break;}
 406     }
 407   }
 408   //
 409   if (nlabels>kMaxLabInCluster) { // only 3 labels are stored in cluster, sort in decreasing momentum
 410     static int ind[kMaxLabels],labSS[kMaxLabels];
 411     TMath::Sort(nlabels,kine,ind);
 412     for (int i=nlabels;i--;) labSS[i] = labS[i];
 413     for (int i=nlabels;i--;) labS[i] = labSS[ind[i]];
 414   }
 415   //
 416   //compress labels -- if multi-times the same
 417   for (Int_t i=0;i<nlabels;i++) fCurrLabels[i]=-2;
 418   fNLabels = 0;
 419   int j=0;
 420   for (int i=0;i<nlabels;i++) {
 421     for (j=fNLabels;j--;) if (labS[i]==fCurrLabels[j]) break; // the label already there
 422     if (j<0) fCurrLabels[fNLabels++] = labS[i];
 423   }
 424   //
 425 }
 426
 427 //------------------------------------------------------------------------
 428
 429 #ifdef _ClusterTopology_
 430 //
 431 //______________________________________________________________________
 432 void AliITSUClusterizer::FillClusterTopology(const AliITSUClusterizerClusterCand *cand, AliITSUClusterPix* cl) const
 433 {
 434   // fill cluster topology bit pattern
 435   //
 436   //
 437   UShort_t mnCol=0xffff,mxCol=0,mnRow=0xffff,mxRow=0;
 438   for (AliITSUClusterizerClusterDigit *idigit=cand->fFirstDigit;idigit;idigit=idigit->fNext) {
 439     AliITSdigit* dig = idigit->fDigit;
 440     if (mnCol>dig->GetCoord1()) mnCol = (UShort_t)dig->GetCoord1();
 441     if (mxCol<dig->GetCoord1()) mxCol = (UShort_t)dig->GetCoord1();
 442     if (mnRow>dig->GetCoord2()) mnRow = (UShort_t)dig->GetCoord2();
 443     if (mxRow<dig->GetCoord2()) mxRow = (UShort_t)dig->GetCoord2();
 444   }
 445   UShort_t nColSpan = UShort_t(mxCol-mnCol+1);
 446   UShort_t nRowSpan = UShort_t(mxRow-mnRow+1);
 447   UShort_t nColSpanW = nColSpan;
 448   UShort_t nRowSpanW = nRowSpan;
 449   if (nColSpan*nRowSpan>AliITSUClusterPix::kMaxPatternBits) { // need to store partial info
 450     // will crtail largest dimension
 451     if (nColSpan>nRowSpan) {
 452       if ( (nColSpanW=AliITSUClusterPix::kMaxPatternBits/nRowSpan)==0 ) {
 453         nColSpanW = 1;
 454         nRowSpanW = AliITSUClusterPix::kMaxPatternBits;
 455       }
 456     }
 457     else {
 458       if ( (nRowSpanW=AliITSUClusterPix::kMaxPatternBits/nColSpan)==0 ) {
 459         nRowSpanW = 1;
 460         nColSpanW = AliITSUClusterPix::kMaxPatternBits;
 461       }
 462     }
 463   }
 464   cl->ResetPattern();
 465   //
 466   cl->SetPatternRowSpan(nRowSpanW,nRowSpanW<nRowSpan);
 467   cl->SetPatternColSpan(nColSpanW,nColSpanW<nColSpan);
 468   cl->SetPatternMinRow(mnRow);
 469   cl->SetPatternMinCol(mnCol);
 470   //
 471   for (AliITSUClusterizerClusterDigit *idigit=cand->fFirstDigit;idigit;idigit=idigit->fNext) {
 472     AliITSdigit* dig = idigit->fDigit;
 473     UInt_t ir = dig->GetCoord2()-mnRow;
 474     UInt_t ic = dig->GetCoord1()-mnCol;
 475     if (ir<nRowSpanW && ic<nColSpanW) cl->SetPixel(ir,ic);
 476   }
 477   //
 478 }
 479
 480 #endif //_ClusterTopology_