ITS/UPGRADE/AliITSUClusterizer.cxx

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