ITS/AliITSClusterFinder.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 //                                                                        //
  17 // Base Class used to find                                                //
  18 // the reconstructed points for ITS                                       //
  19 // See also AliITSClusterFinderSPD, AliITSClusterFinderSDD,               //
  20 // AliITSClusterFinderSDD  AliITSClusterFinderV2                          //
  21 ////////////////////////////////////////////////////////////////////////////
  22
  23 #include "AliITSClusterFinder.h"
  24 #include "AliITSRecPoint.h"
  25 #include "AliITSdigit.h"
  26 #include "AliITSDetTypeRec.h"
  27 #include "AliITSMap.h"
  28
  29 ClassImp(AliITSClusterFinder)
  30
  31 //----------------------------------------------------------------------
  32 AliITSClusterFinder::AliITSClusterFinder():
  33 TObject(),
  34 fDebug(0),
  35 fModule(0),
  36 fDigits(0),
  37 fNdigits(0),
  38 fDetTypeRec(0),
  39 fClusters(0),
  40 fNRawClusters(0),
  41 fMap(0),
  42 fNperMax(0),
  43 fDeclusterFlag(0),
  44 fClusterSize(0),
  45 fNPeaks(-1){
  46     // default cluster finder
  47     // Input:
  48     //   none.
  49     // Output:
  50     //   none.
  51     // Return:
  52     //   A default constructed AliITSCulsterFinder
  53 }
  54 //----------------------------------------------------------------------
  55 AliITSClusterFinder::AliITSClusterFinder(AliITSDetTypeRec* dettyp):
  56 TObject(),
  57 fDebug(0),
  58 fModule(0),
  59 fDigits(0),
  60 fNdigits(0),
  61 fClusters(0),
  62 fNRawClusters(0),
  63 fMap(0),
  64 fNperMax(0),
  65 fDeclusterFlag(0),
  66 fClusterSize(0),
  67 fNPeaks(-1){
  68     // Standard constructor for cluster finder
  69     // Input:
  70     //   AliITSsegmentation *seg  The segmentation class to be used
  71     //   AliITSresponse     *res  The response class to be used
  72     // Output:
  73     //   none.
  74     // Return:
  75     //   A Standard constructed AliITSCulsterFinder
  76
  77     SetNperMax();
  78     SetClusterSize();
  79     SetDeclusterFlag();
  80     fDetTypeRec = dettyp;
  81 }
  82 //----------------------------------------------------------------------
  83 AliITSClusterFinder::AliITSClusterFinder(AliITSDetTypeRec* dettyp,
  84                                          TClonesArray *digits):
  85 TObject(),
  86 fDebug(0),
  87 fModule(0),
  88 fDigits(digits),
  89 fNdigits(0),
  90 fDetTypeRec(dettyp),
  91 fClusters(0),
  92 fNRawClusters(0),
  93 fMap(0),
  94 fNperMax(0),
  95 fDeclusterFlag(0),
  96 fClusterSize(0),
  97 fNPeaks(-1){
  98     // Standard + cluster finder constructor
  99     // Input:
 100     //   AliITSsegmentation *seg  The segmentation class to be used
 101     //   AliITSresponse     *res  The response class to be used
 102     //   TClonesArray    *digits  Array of digits to be used
 103     // Output:
 104     //   none.
 105     // Return:
 106     //   A Standard constructed AliITSCulsterFinder
 107
 108     fNdigits = fDigits->GetEntriesFast();
 109     SetNperMax();
 110     SetClusterSize();
 111     SetDeclusterFlag();
 112 }
 113 //______________________________________________________________________
 114 AliITSClusterFinder::AliITSClusterFinder(const AliITSClusterFinder &source) : TObject(source) {
 115   // Copy constructor
 116   // Copies are not allowed. The method is protected to avoid misuse.
 117   Fatal("AliITSClusterFinder","Copy constructor not allowed\n");
 118 }
 119
 120 //______________________________________________________________________
 121 AliITSClusterFinder& AliITSClusterFinder::operator=(const AliITSClusterFinder& /* source */){
 122   // Assignment operator
 123   // Assignment is not allowed. The method is protected to avoid misuse.
 124   Fatal("= operator","Assignment operator not allowed\n");
 125   return *this;
 126 }
 127
 128 //----------------------------------------------------------------------
 129 AliITSClusterFinder::~AliITSClusterFinder(){
 130     // destructor cluster finder
 131     // Input:
 132     //   none.
 133     // Output:
 134     //   none.
 135     // Return:
 136     //   none.
 137
 138     if(fMap) {delete fMap;}
 139     // Zero local pointers. Other classes own these pointers.
 140     fMap          = 0;
 141     fDigits       = 0;
 142     fNdigits      = 0;
 143     fNRawClusters = 0;
 144     fNperMax      = 0;
 145     fDeclusterFlag= 0;
 146     fClusterSize  = 0;
 147     fNPeaks       = 0;
 148     fDetTypeRec   = 0;
 149
 150 }
 151 //__________________________________________________________________________
 152 void AliITSClusterFinder::InitGeometry(){
 153
 154
 155  //Initialisation of ITS geometry
 156   if(!fDetTypeRec->GetITSgeom()) {
 157     Error("InitGeometry","ITS geom is null!");
 158     return;
 159   }
 160   Int_t mmax=fDetTypeRec->GetITSgeom()->GetIndexMax();
 161   if (mmax>2200) {
 162     Fatal("AliITSClusterFinder","Too many ITS subdetectors !");
 163   }
 164   Int_t m;
 165   for (m=0; m<mmax; m++) {
 166     Int_t lay,lad,det; fDetTypeRec->GetITSgeom()->GetModuleId(m,lay,lad,det);
 167     Float_t x,y,z;     fDetTypeRec->GetITSgeom()->GetTrans(lay,lad,det,x,y,z);
 168     Double_t rot[9];   fDetTypeRec->GetITSgeom()->GetRotMatrix(lay,lad,det,rot);
 169     Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
 170     Double_t ca=TMath::Cos(alpha), sa=TMath::Sin(alpha);
 171     fYshift[m] = x*ca + y*sa;
 172     fZshift[m] = (Double_t)z;
 173     fNdet[m] = (lad-1)*fDetTypeRec->GetITSgeom()->GetNdetectors(lay) + (det-1);
 174     fNlayer[m] = lay-1;
 175   }
 176 }
 177
 178
 179
 180 //----------------------------------------------------------------------
 181 void AliITSClusterFinder::AddCluster(Int_t branch, AliITSRawCluster *c){
 182     // Add a raw cluster copy to the list
 183     // Input:
 184     //   Int_t       branch  The branch to which the cluster is to be added to
 185     //   AliITSRawCluster *c The cluster to be added to the array of clusters
 186     // Output:
 187     //   none.
 188     // Return:
 189     //   none.
 190
 191    if(!fDetTypeRec) {
 192     Error("AddCluster","fDetTypeRec is null!");
 193     return;
 194   }
 195   fDetTypeRec->AddCluster(branch,c);
 196   fNRawClusters++;
 197 }
 198 //----------------------------------------------------------------------
 199 void AliITSClusterFinder::AddCluster(Int_t branch, AliITSRawCluster *c,
 200                                      AliITSRecPoint &rp){
 201     // Add a raw cluster copy to the list and the RecPoint
 202     // Input:
 203     //   Int_t       branch  The branch to which the cluster is to be added to
 204     //   AliITSRawCluster *c The cluster to be added to the array of clusters
 205     //   AliITSRecPoint  &rp The RecPoint to be added to the array of RecPoints
 206     // Output:
 207     //   none.
 208     // Return:
 209     //   none.
 210   if(!fDetTypeRec) {
 211     Error("AddCluster","fDetTypeRec is null!");
 212     return;
 213   }
 214
 215   fDetTypeRec->AddCluster(branch,c);
 216   fNRawClusters++;
 217   fDetTypeRec->AddRecPoint(rp);
 218
 219 }
 220 /*
 221 //______________________________________________________________________
 222 void AliITSClusterFinder::CheckLabels(Int_t lab[3]) {
 223   //------------------------------------------------------------
 224   // Tries to find mother's labels
 225   //------------------------------------------------------------
 226
 227   if(lab[0]<0 && lab[1]<0 && lab[2]<0) return; // In case of no labels just exit
 228   // Check if simulation
 229   AliMC* mc = gAlice->GetMCApp();
 230   if(!mc)return;
 231
 232   Int_t ntracks = mc->GetNtrack();
 233   for (Int_t i=0;i<3;i++){
 234     Int_t label = lab[i];
 235     if (label>=0 && label<ntracks) {
 236       TParticle *part=(TParticle*)mc->Particle(label);
 237       if (part->P() < 0.005) {
 238         Int_t m=part->GetFirstMother();
 239         if (m<0) {
 240           continue;
 241         }
 242         if (part->GetStatusCode()>0) {
 243           continue;
 244         }
 245         lab[i]=m;
 246       }
 247     }
 248   }
 249
 250 }
 251 */
 252 //______________________________________________________________________
 253 void AliITSClusterFinder::FindRawClusters(Int_t module){
 254     // Default Cluster finder.
 255     // Input:
 256     //   Int_t module   Module number for which culster are to be found.
 257     // Output:
 258     //   none.
 259     // Return:
 260     //   none.
 261     const Int_t kelms = 10;
 262     Int_t ndigits = fDigits->GetEntriesFast();
 263     TObjArray *digs = new TObjArray(ndigits);
 264     TObjArray *clusts = new TObjArray(ndigits); // max # cluster
 265     TObjArray *clust0=0; // A spacific cluster of digits
 266     TObjArray *clust1=0; // A spacific cluster of digits
 267     AliITSdigit *dig=0; // locat pointer to a digit
 268     Int_t i=0,nc=0,j[4],k,k2=0;
 269
 270     // Copy all digits for this module into a local TObjArray.
 271     for(i=0;i<ndigits;i++) digs->AddAt(new AliITSdigit(*(GetDigit(i))),i);
 272     digs->Sort();
 273     // First digit is a cluster.
 274     i  = 0;
 275     nc = 0;
 276     clusts->AddAt(new TObjArray(kelms),nc);
 277     clust0 = (TObjArray*)(clusts->At(nc));
 278     clust0->AddAtFree(digs->At(i)); // move owner ship from digs to clusts
 279     nc++;
 280     for(i=1;i<ndigits;i++){
 281         if(IsNeighbor(digs,i,j)){
 282             dig = (AliITSdigit*)(digs->At(j[0]));
 283             // Add to existing cluster. Find which cluster this digis
 284             for(k=0;k<nc;k++){
 285                 clust0 = ((TObjArray*)(clusts->At(k)));
 286                 if(clust0->IndexOf(dig)>=0) break;
 287             } // end for k
 288             if(k>=nc){
 289                 Fatal("FindRawClusters","Digit not found as expected");
 290             } // end if
 291             if(j[1]>=0){
 292                 dig = (AliITSdigit*)(digs->At(j[1]));
 293                 // Add to existing cluster. Find which cluster this digis
 294                 for(k2=0;k2<nc;k2++){
 295                     clust1 = ((TObjArray*)(clusts->At(k2)));
 296                     if(clust1->IndexOf(dig)>=0) break;
 297                 } // end for k2
 298                 if(k2>=nc){
 299                     Fatal("FindRawClusters","Digit not found as expected");
 300                 } // end if
 301             } // end if j[1]>=0
 302             // Found cluster with neighboring digits add this one to it.
 303             if(clust0==clust1){ // same cluster
 304                 clust0->AddAtFree(digs->At(i));
 305                 clust0 = 0; // finished with cluster. zero for safty
 306                 clust1 = 0; // finished wit hcluster. zero for safty
 307             }else{ // two different clusters which need to be merged.
 308                 clust0->AddAtFree(digs->At(i)); // Add digit to this cluster.
 309                 for(k=0;k<clust1->GetEntriesFast();k++){
 310                     // move clust1 into clust0
 311                     //move digit to this cluster
 312                     clust0->AddAtFree(clust1->At(k));
 313                     clust1->AddAt(0,k); // zero this one
 314                 } // end for k
 315                 delete clust1;
 316                 clusts->AddAt(0,k2); // zero array of clusters element clust1
 317                 clust0 = 0; // finished with cluster. zero for safty
 318                 clust1 = 0; // finished wit hcluster. zero for safty
 319             } // end if clust0==clust1
 320         }else{// New cluster
 321             clusts->AddAt(new TObjArray(kelms),nc);
 322             clust0 = ((TObjArray*)(clusts->At(nc)));
 323             // move owner ship from digs to clusts
 324             clust0->AddAtFree(digs->At(i));
 325             clust0 = 0; // finished with cluster. zero for safty
 326             nc++;
 327         } // End if IsNeighbor
 328     } // end for i
 329     // There are now nc clusters in clusts. Each element of clust is an
 330     // array of digits which are clustered together.
 331
 332     // For each cluster call detector specific CreateRecPoints
 333     for(i=0;i<nc;i++) CreateRecPoints((TObjArray*)(clusts->At(i)),module);
 334
 335     // clean up at the end.
 336     for(i=0;i<nc;i++){
 337         clust0 =(TObjArray*)(clusts->At(i));
 338         // Digits deleted below, so zero this TObjArray
 339         for(k=0;k<clust0->GetEntriesFast();k++) clust0->AddAt(0,k);
 340         delete clust0; // Delete this TObjArray
 341         clusts->AddAt(0,i); // Contents deleted above, so zero it.
 342     } // end for i
 343     delete clusts; // Delete this TObjArray/
 344     // Delete the digits then the TObjArray which containted them.
 345     for(i=0;i<ndigits;i++) delete ((AliITSdigit*)(digs->At(i)));
 346     delete digs;
 347 }
 348 //______________________________________________________________________
 349 Bool_t AliITSClusterFinder::IsNeighbor(TObjArray *digs,Int_t i,Int_t n[])const{
 350     // Locagical function which checks to see if digit i has a neighbor.
 351     // If so, then it returns kTRUE and its neighbor index j.
 352     // This routine checks if the digits are side by side or one before the
 353     // other. Requires that the array of digits be in proper order.
 354     // Returns kTRUE in the following cases.
 355     //                 ji   0j   if kdiagonal  j0    0i
 356     //                 00   0i   if kdiagonal  0i    j0
 357     // Inputs:
 358     //    TObjArray *digs   Array to search for neighbors in
 359     //    Int_t      i      Index of digit for which we are searching for
 360     //                      a neighbor of.
 361     // Output:
 362     //    Int_t      j[4]   Index of one or more of the digits which is a
 363     //                      neighbor of digit a index i.
 364     // Return:
 365     //    Bool_t            kTRUE if a neighbor was found kFALSE otherwise.
 366     Int_t ix,jx,iz,jz,j;
 367     const Bool_t kdiagonal=kFALSE;
 368     Bool_t nei[4];
 369
 370     // No neighbors found if array empty.
 371     if(digs->GetEntriesFast()<=0) return kFALSE;
 372     // can not be a digit with first element or elements out or range
 373     if(i<=0 || i>=digs->GetEntriesFast()) return kFALSE;
 374
 375     for(j=0;j<4;j++){n[j] = -1;nei[j]=kFALSE;}
 376     ix = ((AliITSdigit*)(digs->At(i)))->GetCoord1();
 377     iz = ((AliITSdigit*)(digs->At(i)))->GetCoord2();
 378     for(j=0;j<i;j++){
 379         jx = ((AliITSdigit*)(digs->At(j)))->GetCoord1();
 380         jz = ((AliITSdigit*)(digs->At(j)))->GetCoord2();
 381         if(jx+1==ix && jz  ==iz){n[0] = j;nei[0] = kTRUE;}
 382         if(jx  ==ix && jz+1==iz){n[1] = j;nei[1] = kTRUE;}
 383         if(jx+1==ix && jz+1==iz){n[2] = j;nei[2] = kTRUE;}
 384         if(jx+1==ix && jz-1==iz){n[3] = j;nei[3] = kTRUE;}
 385     } // end for k
 386     if(nei[0]||nei[1]) return kTRUE;
 387     if(kdiagonal&&(nei[2]||nei[3])) return kTRUE;
 388     // no Neighbors found.
 389     return kFALSE;
 390 }
 391
 392 //______________________________________________________________________
 393 void AliITSClusterFinder::Print(ostream *os) const{
 394     //Standard output format for this class
 395     // Inputs:
 396     //    ostream *os   Output stream
 397     // Output:
 398     //    ostream *os   Output stream
 399     // Return:
 400     //    none.
 401
 402     *os << fDebug<<",";
 403     *os << fModule<<",";
 404     *os << fNdigits<<",";
 405     *os << fNRawClusters<<",";
 406     *os << fNperMax<<",";
 407     *os << fDeclusterFlag<<",";
 408     *os << fClusterSize<<",";
 409     *os << fNPeaks<<endl;
 410 }
 411 /*
 412 //______________________________________________________________________
 413 void AliITSClusterFinder::RecPoints2Clusters
 414 (const TClonesArray *points, Int_t idx, TClonesArray *clusters) {
 415   //------------------------------------------------------------
 416   // Conversion AliITSRecPoint -> AliITSclusterV2 for the ITS
 417   // subdetector indexed by idx
 418   //------------------------------------------------------------
 419   if(!fDetTypeRec->GetITSgeom()) {
 420     Error("RecPoints2Clusters","ITS geom is null!");
 421     return;
 422   }
 423   Int_t lastSPD1=fDetTypeRec->GetITSgeom()->GetModuleIndex(2,1,1)-1;
 424   TClonesArray &cl=*clusters;
 425   Int_t ncl=points->GetEntriesFast();
 426   for (Int_t i=0; i<ncl; i++) {
 427     AliITSRecPoint *p = (AliITSRecPoint *)points->UncheckedAt(i);
 428     Float_t lp[5];
 429     lp[0]=-(-p->GetX()+fYshift[idx]); if (idx<=lastSPD1) lp[0]*=-1; //SPD1
 430     lp[1]=  -p->GetZ()+fZshift[idx];
 431     lp[2]=p->GetSigmaX2();
 432     lp[3]=p->GetSigmaZ2();
 433     lp[4]=p->GetQ()*36./23333.;  //electrons -> ADC
 434     Int_t lab[4];
 435     lab[0]=p->GetLabel(0); lab[1]=p->GetLabel(1); lab[2]=p->GetLabel(2);
 436     lab[3]=fNdet[idx];
 437     CheckLabels(lab);
 438     Int_t dummy[3]={0,0,0};
 439     new (cl[i]) AliITSclusterV2(lab,lp, dummy);
 440   }
 441 }
 442 */
 443 //______________________________________________________________________
 444 void AliITSClusterFinder::Read(istream *is)  {
 445     //Standard input for this class
 446     // Inputs:
 447     //    istream *is   Input stream
 448     // Output:
 449     //    istream *is   Input stream
 450     // Return:
 451     //    none.
 452
 453     *is >> fDebug;
 454     *is >> fModule;
 455     *is >> fNdigits;
 456     *is >> fNRawClusters;
 457     *is >> fNperMax;
 458     *is >> fDeclusterFlag;
 459     *is >> fClusterSize;
 460     *is >> fNPeaks;
 461 }
 462 //______________________________________________________________________
 463 ostream &operator<<(ostream &os,AliITSClusterFinder &source){
 464     // Standard output streaming function.
 465     // Inputs:
 466     //    ostream             *os     Output stream
 467     //    AliITSClusterFinder &source Class to be printed
 468     // Output:
 469     //    ostream             *os     Output stream
 470     // Return:
 471     //    none.
 472
 473     source.Print(&os);
 474     return os;
 475 }
 476 //______________________________________________________________________
 477 istream &operator>>(istream &is,AliITSClusterFinder &source){
 478     // Standard output streaming function.
 479     // Inputs:
 480     //    istream              *is      Input stream
 481     //     AliITSClusterFinder &source  Class to be read in.
 482     // Output:
 483     //    istream              *is      Input stream
 484     // Return:
 485     //    none.
 486
 487     source.Read(&is);
 488     return is;
 489 }
 490 /*
 491 void AliITSClusterFinder::RecPoints2Clusters
 492 (const TClonesArray *points, Int_t idx, TClonesArray *clusters) {
 493   //------------------------------------------------------------
 494   // Conversion AliITSRecPoint -> AliITSclusterV2 for the ITS
 495   // subdetector indexed by idx
 496   //------------------------------------------------------------
 497   TClonesArray &cl=*clusters;
 498   Int_t ncl=points->GetEntriesFast();
 499   for (Int_t i=0; i<ncl; i++) {
 500     AliITSRecPoint *p = (AliITSRecPoint *)points->UncheckedAt(i);
 501     Float_t lp[5];
 502     lp[0]=-(-p->GetX()+fYshift[idx]); if (idx<=fLastSPD1) lp[0]*=-1; //SPD1
 503     lp[1]=  -p->GetZ()+fZshift[idx];
 504     lp[2]=p->GetSigmaX2();
 505     lp[3]=p->GetSigmaZ2();
 506     lp[4]=p->GetQ()*36./23333.;  //electrons -> ADC
 507     Int_t lab[4];
 508     lab[0]=p->GetLabel(0); lab[1]=p->GetLabel(1); lab[2]=p->GetLabel(2);
 509     lab[3]=fNdet[idx];
 510     CheckLabels(lab);
 511     Int_t dummy[3]={0,0,0};
 512     new (cl[i]) AliITSclusterV2(lab,lp, dummy);
 513   }
 514 }
 515 */