TPC/AliTPCclusterer.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 $Log$
  18 Revision 1.10  2002/10/14 14:57:43  hristov
  19 Merging the VirtualMC branch to the main development branch (HEAD)
  20
  21 Revision 1.7.6.1  2002/06/10 15:26:12  hristov
  22 Merged with v3-08-02
  23
  24 Revision 1.9  2002/05/09 06:57:09  kowal2
  25 Protection against nonexisting input tree
  26
  27 Revision 1.8  2002/03/29 06:57:45  kowal2
  28 Restored backward compatibility to use the hits from Dec. 2000 production.
  29
  30 Revision 1.7  2001/10/21 19:04:55  hristov
  31 Several patches were done to adapt the barel reconstruction to the multi-event case. Some memory leaks were corrected. (Yu.Belikov)
  32
  33 Revision 1.6  2001/08/30 09:28:48  hristov
  34 TTree names are explicitly set via SetName(name) and then Write() is called
  35
  36 Revision 1.5  2001/07/20 14:32:44  kowal2
  37 Processing of many events possible now
  38
  39 Revision 1.4  2001/04/17 08:06:27  hristov
  40 Possibility to define the magnetic field in the reconstruction (Yu.Belikov)
  41
  42 Revision 1.3  2000/10/05 16:14:01  kowal2
  43 Forward declarations.
  44
  45 Revision 1.2  2000/06/30 12:07:50  kowal2
  46 Updated from the TPC-PreRelease branch
  47
  48 Revision 1.1.2.1  2000/06/25 08:53:55  kowal2
  49 Splitted from AliTPCtracking
  50
  51 */
  52
  53 //-------------------------------------------------------
  54 //          Implementation of the TPC clusterer
  55 //
  56 //   Origin: Jouri Belikov, CERN, Jouri.Belikov@cern.ch
  57 //-------------------------------------------------------
  58
  59 #include "AliTPCclusterer.h"
  60 #include "AliTPCcluster.h"
  61 #include <TObjArray.h>
  62 #include <TFile.h>
  63 #include "AliTPCClustersArray.h"
  64 #include "AliTPCClustersRow.h"
  65 #include "AliDigits.h"
  66 #include "AliSimDigits.h"
  67 #include "AliTPCParam.h"
  68 #include <Riostream.h>
  69 #include <TTree.h>
  70
  71 void AliTPCclusterer::FindPeaks(Int_t k,Int_t max,
  72 AliBin *b,Int_t *idx,UInt_t *msk,Int_t& n) {
  73   //find local maxima
  74   if (n<31)
  75   if (IsMaximum(k,max,b)) {
  76     idx[n]=k; msk[n]=(2<<n);
  77     n++;
  78   }
  79   b[k].SetMask(0);
  80   if (b[k-max].GetMask()&1) FindPeaks(k-max,max,b,idx,msk,n);
  81   if (b[k-1  ].GetMask()&1) FindPeaks(k-1  ,max,b,idx,msk,n);
  82   if (b[k+max].GetMask()&1) FindPeaks(k+max,max,b,idx,msk,n);
  83   if (b[k+1  ].GetMask()&1) FindPeaks(k+1  ,max,b,idx,msk,n);
  84 }
  85
  86 void AliTPCclusterer::MarkPeak(Int_t k, Int_t max, AliBin *bins, UInt_t m) {
  87   //mark this peak
  88   UShort_t q=bins[k].GetQ();
  89
  90   bins[k].SetMask(bins[k].GetMask()|m);
  91
  92   if (bins[k-max].GetQ() <= q)
  93      if ((bins[k-max].GetMask()&m) == 0) MarkPeak(k-max,max,bins,m);
  94   if (bins[k-1  ].GetQ() <= q)
  95      if ((bins[k-1  ].GetMask()&m) == 0) MarkPeak(k-1  ,max,bins,m);
  96   if (bins[k+max].GetQ() <= q)
  97      if ((bins[k+max].GetMask()&m) == 0) MarkPeak(k+max,max,bins,m);
  98   if (bins[k+1  ].GetQ() <= q)
  99      if ((bins[k+1  ].GetMask()&m) == 0) MarkPeak(k+1  ,max,bins,m);
 100 }
 101
 102 void AliTPCclusterer::MakeCluster(Int_t k,Int_t max,AliBin *bins,UInt_t m,
 103 AliTPCcluster &c) {
 104   //make cluster using digits of this peak
 105   Float_t q=(Float_t)bins[k].GetQ();
 106   Int_t i=k/max, j=k-i*max;
 107
 108   c.SetQ(c.GetQ()+q);
 109   c.SetY(c.GetY()+i*q);
 110   c.SetZ(c.GetZ()+j*q);
 111   c.SetSigmaY2(c.GetSigmaY2()+i*i*q);
 112   c.SetSigmaZ2(c.GetSigmaZ2()+j*j*q);
 113
 114   bins[k].SetMask(0xFFFFFFFE);
 115
 116   if (bins[k-max].GetMask() == m) MakeCluster(k-max,max,bins,m,c);
 117   if (bins[k-1  ].GetMask() == m) MakeCluster(k-1  ,max,bins,m,c);
 118   if (bins[k+max].GetMask() == m) MakeCluster(k+max,max,bins,m,c);
 119   if (bins[k+1  ].GetMask() == m) MakeCluster(k+1  ,max,bins,m,c);
 120 }
 121
 122 //_____________________________________________________________________________
 123 void AliTPCclusterer::Digits2Clusters(const AliTPCParam *par, TFile *of, Int_t eventn)
 124 {
 125   //-----------------------------------------------------------------
 126   // This is a simple cluster finder.
 127   //-----------------------------------------------------------------
 128   TDirectory *savedir=gDirectory;
 129
 130   if (!of->IsOpen()) {
 131      cerr<<"AliTPC::Digits2Clusters(): output file not open !\n";
 132      return;
 133   }
 134
 135   const Int_t kMAXZ=par->GetMaxTBin()+2;
 136
 137   char  dname[100];
 138   char   cname[100];
 139   if (eventn==-1) {
 140
 141     // for backward compatibility
 142
 143     sprintf(dname,"TreeD_75x40_100x60_150x60");
 144     sprintf(cname,"TreeC_TPC");
 145   }
 146   else {
 147     sprintf(dname,"TreeD_75x40_100x60_150x60_%d",eventn);
 148     sprintf(cname,"TreeC_TPC_%d",eventn);
 149   }
 150   TTree *t = (TTree *)gDirectory->Get(dname);
 151
 152   if (!t) {
 153     cerr<<"Input tree with "<<dname<<" not found"<<endl;
 154     return;
 155   }
 156
 157   AliSimDigits digarr, *dummy=&digarr;
 158   t->GetBranch("Segment")->SetAddress(&dummy);
 159   Stat_t nentries = t->GetEntries();
 160
 161   of->cd();
 162
 163   ((AliTPCParam*)par)->Write(par->GetTitle());
 164   AliTPCClustersArray carray;
 165   carray.Setup(par);
 166   carray.SetClusterType("AliTPCcluster");
 167   carray.MakeTree();
 168
 169
 170
 171   Int_t nclusters=0;
 172
 173   for (Int_t n=0; n<nentries; n++) {
 174     t->GetEvent(n);
 175     Int_t sec, row;
 176     if (!par->AdjustSectorRow(digarr.GetID(),sec,row)) {
 177        cerr<<"AliTPC warning: invalid segment ID ! "<<digarr.GetID()<<endl;
 178        continue;
 179     }
 180
 181     AliTPCClustersRow *clrow=carray.CreateRow(sec,row);
 182
 183     Float_t rx=par->GetPadRowRadii(sec,row);
 184
 185     Int_t npads, sign;
 186     {
 187        const Int_t kNIS=par->GetNInnerSector(), kNOS=par->GetNOuterSector();
 188        if (sec < kNIS) {
 189           npads = par->GetNPadsLow(row);
 190           sign = (sec < kNIS/2) ? 1 : -1;
 191        } else {
 192           npads = par->GetNPadsUp(row);
 193           sign = ((sec-kNIS) < kNOS/2) ? 1 : -1;
 194        }
 195     }
 196
 197     const Int_t kMAXBIN=kMAXZ*(npads+2);
 198     AliBin *bins=new AliBin[kMAXBIN];
 199     for (Int_t ii=0;ii<kMAXBIN;ii++) {
 200        bins[ii].SetQ(0); bins[ii].SetMask(0xFFFFFFFE);
 201     }
 202
 203     digarr.First();
 204     do {
 205        Short_t dig=digarr.CurrentDigit();
 206        if (dig<=par->GetZeroSup()) continue;
 207        Int_t j=digarr.CurrentRow()+1, i=digarr.CurrentColumn()+1;
 208        bins[i*kMAXZ+j].SetQ(dig);
 209        bins[i*kMAXZ+j].SetMask(1);
 210     } while (digarr.Next());
 211
 212     Int_t ncl=0;
 213     for (Int_t i=0; i<kMAXBIN; i++) {
 214       if ((bins[i].GetMask()&1) == 0) continue;
 215       Int_t idx[32]; UInt_t msk[32]; Int_t npeaks=0;
 216       FindPeaks(i, kMAXZ, bins, idx, msk, npeaks);
 217
 218       if (npeaks>30) continue;
 219
 220       Int_t k,l;
 221       for (k=0; k<npeaks-1; k++){//mark adjacent peaks
 222         if (idx[k] < 0) continue; //this peak is already removed
 223         for (l=k+1; l<npeaks; l++) {
 224            if (idx[l] < 0) continue; //this peak is already removed
 225            Int_t ki=idx[k]/kMAXZ, kj=idx[k] - ki*kMAXZ;
 226            Int_t li=idx[l]/kMAXZ, lj=idx[l] - li*kMAXZ;
 227            Int_t di=TMath::Abs(ki - li);
 228            Int_t dj=TMath::Abs(kj - lj);
 229            if (di>1 || dj>1) continue;
 230            if (bins[idx[k]].GetQ() > bins[idx[l]].GetQ()) {
 231               msk[l]=msk[k];
 232               idx[l]*=-1;
 233            } else {
 234               msk[k]=msk[l];
 235               idx[k]*=-1;
 236               break;
 237            }
 238         }
 239       }
 240
 241       for (k=0; k<npeaks; k++) {
 242         MarkPeak(TMath::Abs(idx[k]), kMAXZ, bins, msk[k]);
 243       }
 244
 245       for (k=0; k<npeaks; k++) {
 246          if (idx[k] < 0) continue; //removed peak
 247          AliTPCcluster c;
 248          MakeCluster(idx[k], kMAXZ, bins, msk[k], c);
 249          if (c.GetQ() < 5) continue; //noise cluster
 250          c.SetY(c.GetY()/c.GetQ());
 251          c.SetZ(c.GetZ()/c.GetQ());
 252
 253          Float_t s2 = c.GetSigmaY2()/c.GetQ() - c.GetY()*c.GetY();
 254          Float_t w=par->GetPadPitchWidth(sec);
 255          c.SetSigmaY2((s2 + 1./12.)*w*w);
 256          if (s2 != 0.) {
 257            c.SetSigmaY2(c.GetSigmaY2()*0.108);
 258            if (sec<par->GetNInnerSector()) c.SetSigmaY2(c.GetSigmaY2()*2.07);
 259          }
 260
 261          s2 = c.GetSigmaZ2()/c.GetQ() - c.GetZ()*c.GetZ();
 262          w=par->GetZWidth();
 263          c.SetSigmaZ2((s2 + 1./12.)*w*w);
 264          if (s2 != 0.) {
 265            c.SetSigmaZ2(c.GetSigmaZ2()*0.169);
 266            if (sec<par->GetNInnerSector()) c.SetSigmaZ2(c.GetSigmaZ2()*1.77);
 267          }
 268
 269          c.SetY((c.GetY() - 0.5 - 0.5*npads)*par->GetPadPitchWidth(sec));
 270          c.SetZ(par->GetZWidth()*(c.GetZ()-1));
 271          c.SetZ(c.GetZ() - 3.*par->GetZSigma()); // PASA delay
 272          c.SetZ(sign*(par->GetZLength() - c.GetZ()));
 273
 274          if (rx<230./250.*TMath::Abs(c.GetZ())) continue;
 275
 276          Int_t ki=idx[k]/kMAXZ, kj=idx[k] - ki*kMAXZ;
 277          c.SetLabel(digarr.GetTrackID(kj-1,ki-1,0),0);
 278          c.SetLabel(digarr.GetTrackID(kj-1,ki-1,1),1);
 279          c.SetLabel(digarr.GetTrackID(kj-1,ki-1,2),2);
 280
 281          c.SetQ(bins[idx[k]].GetQ());
 282
 283          if (ki==1 || ki==npads || kj==1 || kj==kMAXZ-2) {
 284            c.SetSigmaY2(c.GetSigmaY2()*25.);
 285            c.SetSigmaZ2(c.GetSigmaZ2()*4.);
 286          }
 287          clrow->InsertCluster(&c); ncl++;
 288       }
 289     }
 290     carray.StoreRow(sec,row);
 291     carray.ClearRow(sec,row);
 292
 293     //cerr<<"sector, row, compressed digits, clusters: "
 294     //<<sec<<' '<<row<<' '<<digarr.GetSize()<<' '<<ncl<<"                  \r";
 295
 296     nclusters+=ncl;
 297
 298     delete[] bins;
 299   }
 300
 301   cerr<<"Number of found clusters : "<<nclusters<<"                        \n";
 302
 303   carray.GetTree()->SetName(cname);
 304   carray.GetTree()->Write();
 305   savedir->cd();
 306
 307   delete t;  //Thanks to Mariana Bondila
 308 }
 309