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