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   UChar_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   cluster->SetNxNzN(nx,nz,n);
 228   cluster->SetQ(charge); // note: this is MC info
 229   //
 230   if (!fRawData) {
 231     CheckLabels();
 232     int nl = Min(kMaxLabInCluster,fNLabels);
 233     for (int i=nl;i--;) cluster->SetLabel(fCurrLabels[i],i);
 234   }
 235   //
 236   // Set Volume id
 237   cluster->SetVolumeId(fVolID);
 238   //    printf("mod %d: (%.4lf,%.4lf)cm\n",fVolID,x,z);
 239   //
 240 #ifdef _ClusterTopology_
 241   FillClusterTopology(cand,cluster);
 242 #endif //_ClusterTopology_
 243   //
 244 }
 245
 246 //______________________________________________________________________________
 247 void AliITSUClusterizer::CloseCand(AliITSUClusterizerClusterCand *cand)
 248 {
 249   // finish cluster
 250   AliITSUClusterPix *cluster = (AliITSUClusterPix*)NextCluster();
 251   Transform(cluster,cand);
 252   DeallocDigits(cand->fFirstDigit,cand->fLastDigit);
 253   DeallocCand(cand);
 254 }
 255
 256 //______________________________________________________________________________
 257 void AliITSUClusterizer::ClosePart(AliITSUClusterizerClusterPart *part)
 258 {
 259   // finish cluster part
 260   AliITSUClusterizerClusterCand *cand=part->fParent;
 261   DetachPartFromCand(cand,part);
 262   DeallocPart(part);
 263   if (!cand->fFirstPart) CloseCand(cand);
 264 }
 265
 266 //______________________________________________________________________________
 267 void AliITSUClusterizer::CloseRemainingParts(AliITSUClusterizerClusterPart *part)
 268 {
 269   // finish what is left
 270   while (part) {
 271     AliITSUClusterizerClusterPart *next=part->fNextInRow;
 272     ClosePart(part);
 273     part=next;
 274   }
 275 }
 276
 277 //______________________________________________________________________________
 278 void AliITSUClusterizer::Clusterize()
 279 {
 280   // main algo for MC clustererization
 281   SetRawData(kFALSE);
 282   //
 283   AliITSUClusterizerClusterDigit *iDigit=NextDigit();
 284   AliITSUClusterizerClusterPart *iPrevRowBegin=0;
 285   AliITSUClusterizerClusterPart *iNextRowBegin=0;
 286   AliITSUClusterizerClusterPart *iPrevRow=0;
 287   AliITSUClusterizerClusterPart *iNextRow=0;
 288   Int_t lastV=0;
 289   //
 290   while (iDigit) {
 291     if (iDigit->fDigit->GetCoord2()!=lastV) {
 292       // NEW ROW
 293       if (iNextRow) iNextRow->fNextInRow=0;
 294       if (iPrevRowBegin) CloseRemainingParts(iPrevRowBegin);
 295       if (iDigit->fDigit->GetCoord2()==lastV+1) {
 296         iPrevRowBegin=iNextRowBegin;
 297         iPrevRow     =iNextRowBegin;
 298       }
 299       else {
 300         // there was an empty row
 301         CloseRemainingParts(iNextRowBegin);
 302         iPrevRowBegin=0;
 303         iPrevRow     =0;
 304       }
 305       iNextRowBegin=0;
 306       iNextRow     =0;
 307       lastV=iDigit->fDigit->GetCoord2();
 308     }
 309     // skip cluster parts before this digit
 310     int limCol = iDigit->fDigit->GetCoord1()-fAllowDiagonalClusterization;
 311     while (iPrevRow && iPrevRow->fUEnd < limCol) {
 312       iPrevRow=iPrevRow->fNextInRow;
 313     }
 314     // find the longest continous line of digits [iDigit,pDigit]=[iDigit,jDigit)
 315     AliITSUClusterizerClusterCand *cand=AllocCand();
 316     AliITSUClusterizerClusterDigit *pDigit=iDigit;
 317     AliITSUClusterizerClusterDigit *jDigit=NextDigit();
 318     cand->fFirstPart=0;
 319     cand->fFirstDigit=cand->fLastDigit=iDigit; // NB: first diggit is attached differently
 320     iDigit->fNext=0;
 321     Int_t lastU =iDigit->fDigit->GetCoord1();
 322     Int_t lastU1=lastU+1;
 323     while (jDigit && jDigit->fDigit->GetCoord1()==lastU1 && jDigit->fDigit->GetCoord2()==lastV) {
 324       pDigit=jDigit;
 325       jDigit=NextDigit();
 326       AttachDigitToCand(cand,pDigit);
 327       ++lastU1;
 328     }
 329     if (!fAllowDiagonalClusterization) --lastU1;
 330     AliITSUClusterizerClusterPart *part=AllocPart();
 331     part->fUBegin=lastU ;
 332     part->fUEnd  =lastU1;
 333     AttachPartToCand(cand,part);
 334     // merge all cluster candidates of the previous line touching this one,
 335     // advance to the last one, but keep that one the next active one
 336     AliITSUClusterizerClusterPart *jPrevRow=iPrevRow;
 337     while (jPrevRow && jPrevRow->fUBegin<=lastU1) {
 338       if (jPrevRow->fParent!=cand) {
 339         MergeCands(jPrevRow->fParent,cand);
 340         DeallocCand(cand);
 341         cand=jPrevRow->fParent;
 342       }
 343       iPrevRow=jPrevRow;
 344       jPrevRow=jPrevRow->fNextInRow;
 345     }
 346     if (iNextRow)
 347       iNextRow->fNextInRow=part;
 348     else
 349       iNextRowBegin=part;
 350     iNextRow=part;
 351     iDigit=jDigit;
 352   }
 353   // remove remaining cluster parts
 354   CloseRemainingParts(iPrevRowBegin);
 355   if (iNextRow) iNextRow->fNextInRow=0;
 356   CloseRemainingParts(iNextRowBegin);
 357   return;
 358 }
 359
 360 //______________________________________________________________________________
 361 void AliITSUClusterizer::PrepareLorentzAngleCorrection(Double_t bz)
 362 {
 363   // calculate parameters for Lorentz Angle correction. Must be called
 364   // after setting segmentation and recoparams
 365   fLorAngCorrection = 0.5*fRecoParam->GetTanLorentzAngle(fLayerID)*bz/kNominalBz*fSegm->Dy();
 366 }
 367
 368 //______________________________________________________________________
 369 void AliITSUClusterizer::CheckLabels()
 370 {
 371   // Tries to find mother's labels
 372   //
 373   if (fNLabels<1) return;
 374   AliRunLoader *rl = AliRunLoader::Instance();
 375   if(!rl) return;
 376   TTree *trK=(TTree*)rl->TreeK();
 377   if (!trK) return;
 378   //
 379   static int   labS[kMaxLabels];
 380   static float kine[kMaxLabels];
 381   Int_t nlabels = fNLabels;
 382   Int_t ntracks = gAlice->GetMCApp()->GetNtrack();
 383   for (Int_t i=fNLabels;i--;) labS[i] = fCurrLabels[i];
 384   //
 385   for (Int_t i=0;i<nlabels;i++) {
 386     Int_t label = labS[i];
 387     if (label>=ntracks) continue;
 388     TParticle *part=(TParticle*)gAlice->GetMCApp()->Particle(label);
 389     kine[i] = part->Energy() - part->GetCalcMass(); // kinetic energy
 390     if (kine[i] < 0.02) {    // reduce soft particles from the same cluster
 391       Int_t m=part->GetFirstMother();
 392       if (m<0) continue; // primary
 393       //
 394       if (part->GetStatusCode()>0) continue;
 395       //
 396       // if the parent is within the same cluster, assign parent's label
 397       for (int j=0;j<nlabels;j++) if (labS[j]==m) { labS[i] = m; break;}
 398     }
 399   }
 400   //
 401   if (nlabels>kMaxLabInCluster) { // only 3 labels are stored in cluster, sort in decreasing momentum
 402     static int ind[kMaxLabels],labSS[kMaxLabels];
 403     TMath::Sort(nlabels,kine,ind);
 404     for (int i=nlabels;i--;) labSS[i] = labS[i];
 405     for (int i=nlabels;i--;) labS[i] = labSS[ind[i]];
 406   }
 407   //
 408   //compress labels -- if multi-times the same
 409   for (Int_t i=0;i<nlabels;i++) fCurrLabels[i]=-2;
 410   fNLabels = 0;
 411   int j=0;
 412   for (int i=0;i<nlabels;i++) {
 413     for (j=fNLabels;j--;) if (labS[i]==fCurrLabels[j]) break; // the label already there
 414     if (j<0) fCurrLabels[fNLabels++] = labS[i];
 415   }
 416   //
 417 }
 418
 419 //------------------------------------------------------------------------
 420
 421 #ifdef _ClusterTopology_
 422 //
 423 //______________________________________________________________________
 424 void AliITSUClusterizer::FillClusterTopology(const AliITSUClusterizerClusterCand *cand, AliITSUClusterPix* cl) const
 425 {
 426   // fill cluster topology bit pattern
 427   //
 428   //
 429   UShort_t mnCol=0xffff,mxCol=0,mnRow=0xffff,mxRow=0;
 430   for (AliITSUClusterizerClusterDigit *idigit=cand->fFirstDigit;idigit;idigit=idigit->fNext) {
 431     AliITSdigit* dig = idigit->fDigit;
 432     if (mnCol>dig->GetCoord1()) mnCol = (UShort_t)dig->GetCoord1();
 433     if (mxCol<dig->GetCoord1()) mxCol = (UShort_t)dig->GetCoord1();
 434     if (mnRow>dig->GetCoord2()) mnRow = (UShort_t)dig->GetCoord2();
 435     if (mxRow<dig->GetCoord2()) mxRow = (UShort_t)dig->GetCoord2();
 436   }
 437   UShort_t nColSpan = UShort_t(mxCol-mnCol+1);
 438   UShort_t nRowSpan = UShort_t(mxRow-mnRow+1);
 439   UShort_t nColSpanW = nColSpan;
 440   UShort_t nRowSpanW = nRowSpan;
 441   if (nColSpan*nRowSpan>AliITSUClusterPix::kMaxPatternBits) { // need to store partial info
 442     // will crtail largest dimension
 443     if (nColSpan>nRowSpan) {
 444       if ( (nColSpanW=AliITSUClusterPix::kMaxPatternBits/nRowSpan)==0 ) {
 445         nColSpanW = 1;
 446         nRowSpanW = AliITSUClusterPix::kMaxPatternBits;
 447       }
 448     }
 449     else {
 450       if ( (nRowSpanW=AliITSUClusterPix::kMaxPatternBits/nColSpan)==0 ) {
 451         nRowSpanW = 1;
 452         nColSpanW = AliITSUClusterPix::kMaxPatternBits;
 453       }
 454     }
 455   }
 456   cl->ResetPattern();
 457   //
 458   cl->SetPatternRowSpan(nRowSpanW,nRowSpanW<nRowSpan);
 459   cl->SetPatternColSpan(nColSpanW,nColSpanW<nColSpan);
 460   cl->SetPatternMinRow(mnRow);
 461   cl->SetPatternMinCol(mnCol);
 462   //
 463   for (AliITSUClusterizerClusterDigit *idigit=cand->fFirstDigit;idigit;idigit=idigit->fNext) {
 464     AliITSdigit* dig = idigit->fDigit;
 465     UInt_t ir = dig->GetCoord2()-mnRow;
 466     UInt_t ic = dig->GetCoord1()-mnCol;
 467     if (ir<nRowSpanW && ic<nColSpanW) cl->SetPixel(ir,ic);
 468   }
 469   //
 470 }
 471
 472 #endif //_ClusterTopology_