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