| 7f4044f0 |
1 | //------------------------------------------------------------------------- |
| 2 | // Implementation of the ITS clusterer V2 class |
| 3 | // |
| 4 | // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch |
| 5 | //------------------------------------------------------------------------- |
| 6 | //uncomment the line below for running with V1 cluster finder classes |
| 7 | //#define V1 |
| 8 | |
| 88cb7938 |
9 | //PH 19/05/2003 This class hast to be adapted to NewIO |
| 10 | |
| ff5af022 |
11 | #include <stdlib.h> |
| 12 | |
| 7f4044f0 |
13 | #include "AliRun.h" |
| 14 | |
| 15 | #include "AliITSclustererV2.h" |
| 16 | #include "AliITSclusterV2.h" |
| c391f9d9 |
17 | #include "AliITSRawStreamSPD.h" |
| 18 | #include "AliITSRawStreamSDD.h" |
| 19 | #include "AliITSRawStreamSSD.h" |
| 7f4044f0 |
20 | |
| 21 | #include <Riostream.h> |
| 22 | #include <TFile.h> |
| 23 | #include <TTree.h> |
| 24 | #include <TClonesArray.h> |
| 25 | #include "AliITSgeom.h" |
| 26 | #include "AliITSdigit.h" |
| 27 | |
| 28 | ClassImp(AliITSclustererV2) |
| 29 | |
| 30 | extern AliRun *gAlice; |
| 31 | |
| 32 | AliITSclustererV2::AliITSclustererV2(const AliITSgeom *geom) { |
| 33 | //------------------------------------------------------------ |
| 34 | // Standard constructor |
| 35 | //------------------------------------------------------------ |
| 36 | AliITSgeom *g=(AliITSgeom*)geom; |
| 37 | |
| 38 | fEvent=0; |
| 39 | fI=0; |
| 40 | |
| 41 | Int_t mmax=geom->GetIndexMax(); |
| 42 | if (mmax>2200) {cerr<<"Too many ITS subdetectors !\n"; exit(1);} |
| 43 | Int_t m; |
| 44 | for (m=0; m<mmax; m++) { |
| 45 | Int_t lay,lad,det; g->GetModuleId(m,lay,lad,det); |
| 46 | Float_t x,y,z; g->GetTrans(lay,lad,det,x,y,z); |
| 47 | Double_t rot[9]; g->GetRotMatrix(lay,lad,det,rot); |
| 48 | fYshift[m] = x*rot[0] + y*rot[1]; |
| 49 | fZshift[m] = (Double_t)z; |
| 50 | fNdet[m] = (lad-1)*g->GetNdetectors(lay) + (det-1); |
| 51 | } |
| c391f9d9 |
52 | fNModules = g->GetIndexMax(); |
| 7f4044f0 |
53 | |
| 54 | //SPD geometry |
| 55 | fLastSPD1=g->GetModuleIndex(2,1,1)-1; |
| 56 | fNySPD=256; fNzSPD=160; |
| 57 | fYpitchSPD=0.0050; |
| 58 | fZ1pitchSPD=0.0425; fZ2pitchSPD=0.0625; |
| 59 | fHwSPD=0.64; fHlSPD=3.48; |
| 60 | fYSPD[0]=0.5*fYpitchSPD; |
| 61 | for (m=1; m<fNySPD; m++) fYSPD[m]=fYSPD[m-1]+fYpitchSPD; |
| 62 | fZSPD[0]=fZ1pitchSPD; |
| 63 | for (m=1; m<fNzSPD; m++) { |
| 64 | Double_t dz=fZ1pitchSPD; |
| 65 | if (m==31 || m==32 || m==63 || m==64 || m==95 || m==96 || |
| 66 | m==127 || m==128) dz=fZ2pitchSPD; |
| 67 | fZSPD[m]=fZSPD[m-1]+dz; |
| 68 | } |
| 69 | for (m=0; m<fNzSPD; m++) { |
| 70 | Double_t dz=0.5*fZ1pitchSPD; |
| 71 | if (m==31 || m==32 || m==63 || m==64 || m==95 || m==96 || |
| 72 | m==127 || m==128) dz=0.5*fZ2pitchSPD; |
| 73 | fZSPD[m]-=dz; |
| 74 | } |
| 75 | |
| 76 | //SDD geometry |
| 77 | fNySDD=256; fNzSDD=256; |
| 78 | fYpitchSDD=0.01825; |
| 79 | fZpitchSDD=0.02940; |
| 80 | fHwSDD=3.5085; fHlSDD=3.7632; |
| 81 | fYoffSDD=0.0425; |
| 82 | |
| 83 | //SSD geometry |
| 84 | fLastSSD1=g->GetModuleIndex(6,1,1)-1; |
| 85 | fYpitchSSD=0.0095; |
| 86 | fHwSSD=3.65; |
| 87 | fHlSSD=2.00; |
| 88 | fTanP=0.0275; |
| 89 | fTanN=0.0075; |
| 90 | } |
| 91 | |
| 92 | void AliITSclustererV2::Digits2Clusters(const TFile *in, TFile *out) { |
| 93 | //------------------------------------------------------------ |
| 94 | // This function creates ITS clusters |
| 95 | //------------------------------------------------------------ |
| 96 | Int_t ncl=0; |
| 97 | TDirectory *savedir=gDirectory; |
| 98 | |
| 99 | if (!out->IsOpen()) { |
| 100 | cerr<<"AliITSclustererV2::Digits2Clusters(): output file not open !\n"; |
| 101 | return; |
| 102 | } |
| 103 | |
| 104 | Char_t name[100]; |
| 105 | sprintf(name,"TreeD%d",fEvent); |
| 106 | |
| 107 | //TTree *dTree=(TTree*)((TFile*)in)->Get(name); |
| 108 | TTree *dTree=gAlice->TreeD(); |
| 109 | if (!dTree) { |
| 110 | cerr<<"Input tree "<<name<<" not found !\n"; |
| 111 | return; |
| 112 | } |
| 113 | |
| 114 | TClonesArray *digitsSPD=new TClonesArray("AliITSdigitSPD",3000); |
| 115 | dTree->SetBranchAddress("ITSDigitsSPD",&digitsSPD); |
| 116 | TClonesArray *digitsSDD=new TClonesArray("AliITSdigitSDD",3000); |
| 117 | dTree->SetBranchAddress("ITSDigitsSDD",&digitsSDD); |
| 118 | TClonesArray *digitsSSD=new TClonesArray("AliITSdigitSSD",3000); |
| 119 | dTree->SetBranchAddress("ITSDigitsSSD",&digitsSSD); |
| 120 | |
| 121 | Int_t mmax=(Int_t)dTree->GetEntries(); |
| 122 | |
| 123 | out->cd(); |
| 124 | |
| 125 | sprintf(name,"TreeC_ITS_%d",fEvent); |
| 126 | TTree cTree(name,"ITS clusters V2"); |
| 127 | TClonesArray *clusters=new TClonesArray("AliITSclusterV2",1000); |
| 128 | cTree.Branch("Clusters",&clusters); |
| 129 | |
| 130 | for (fI=0; fI<mmax; fI++) { |
| 131 | dTree->GetEvent(fI); |
| 132 | |
| 133 | if (digitsSPD->GetEntriesFast()!=0) |
| 134 | FindClustersSPD(digitsSPD,clusters); |
| 135 | else if(digitsSDD->GetEntriesFast()!=0) |
| 136 | FindClustersSDD(digitsSDD,clusters); |
| 137 | else if(digitsSSD->GetEntriesFast()!=0) |
| 138 | FindClustersSSD(digitsSSD,clusters); |
| 139 | |
| 140 | ncl+=clusters->GetEntriesFast(); |
| 141 | |
| 142 | cTree.Fill(); |
| 143 | |
| 144 | digitsSPD->Clear(); |
| 145 | digitsSDD->Clear(); |
| 146 | digitsSSD->Clear(); |
| 147 | clusters->Clear(); |
| 148 | } |
| 149 | cTree.Write(); |
| 150 | |
| 151 | delete clusters; |
| 152 | |
| 153 | delete digitsSPD; |
| 154 | delete digitsSDD; |
| 155 | delete digitsSSD; |
| 156 | |
| 157 | //delete dTree; |
| 158 | |
| 159 | cerr<<"Number of found clusters : "<<ncl<<endl; |
| 160 | |
| 161 | savedir->cd(); |
| 162 | } |
| 163 | |
| c391f9d9 |
164 | void AliITSclustererV2::Digits2Clusters(TFile *out) { |
| 165 | //------------------------------------------------------------ |
| 166 | // This function creates ITS clusters from raw data |
| 167 | //------------------------------------------------------------ |
| 168 | TDirectory *savedir=gDirectory; |
| 169 | if (!out->IsOpen()) { |
| 170 | cerr<<"AliITSclustererV2::Digits2Clusters(): output file not open !\n"; |
| 171 | return; |
| 172 | } |
| 173 | out->cd(); |
| 174 | |
| 175 | Char_t name[100]; |
| 176 | sprintf(name,"TreeC_ITS_%d",fEvent); |
| 177 | TTree cTree(name,"ITS clusters V2"); |
| 178 | TClonesArray *array=new TClonesArray("AliITSclusterV2",1000); |
| 179 | cTree.Branch("Clusters",&array); |
| 180 | delete array; |
| 181 | |
| f1e0c1c5 |
182 | TClonesArray** clusters = new TClonesArray*[fNModules]; |
| c391f9d9 |
183 | // one TClonesArray per module |
| 184 | |
| 185 | AliITSRawStreamSPD inputSPD; |
| 186 | FindClustersSPD(&inputSPD, clusters); |
| 187 | AliITSRawStreamSDD inputSDD; |
| 188 | FindClustersSDD(&inputSDD, clusters); |
| 189 | AliITSRawStreamSSD inputSSD; |
| 190 | FindClustersSSD(&inputSSD, clusters); |
| 191 | |
| 192 | // write all clusters to the tree |
| 193 | Int_t nClusters = 0; |
| 194 | for (Int_t iModule = 0; iModule < fNModules; iModule++) { |
| 195 | TClonesArray* array = clusters[iModule]; |
| 196 | if (!array) { |
| 197 | Error("Digits2Clusters", "data for module %d missing!", iModule); |
| 198 | array = new TClonesArray("AliITSclusterV2"); |
| 199 | } |
| 200 | cTree.SetBranchAddress("Clusters", &array); |
| 201 | cTree.Fill(); |
| 202 | nClusters += array->GetEntriesFast(); |
| 203 | delete array; |
| 204 | } |
| 205 | cTree.Write(); |
| 206 | |
| 207 | savedir->cd(); |
| 208 | |
| 209 | Info("Digits2Clusters", "total number of found clusters in ITS: %d\n", |
| 210 | nClusters); |
| 211 | } |
| 212 | |
| 7f4044f0 |
213 | //**** Fast clusters ******************************* |
| 214 | #include "TParticle.h" |
| 215 | |
| 216 | #include "AliITS.h" |
| 217 | #include "AliITSmodule.h" |
| 218 | #include "AliITSRecPoint.h" |
| 219 | #include "AliITSsimulationFastPoints.h" |
| 220 | #include "AliITSRecPoint.h" |
| 221 | |
| 222 | static void CheckLabels(Int_t lab[3]) { |
| 223 | //------------------------------------------------------------ |
| 224 | // Tries to find mother's labels |
| 225 | //------------------------------------------------------------ |
| 226 | Int_t label=lab[0]; |
| 227 | if (label>=0) { |
| 228 | TParticle *part=(TParticle*)gAlice->Particle(label); |
| 229 | label=-3; |
| 230 | while (part->P() < 0.005) { |
| 231 | Int_t m=part->GetFirstMother(); |
| 232 | if (m<0) {cerr<<"Primary momentum: "<<part->P()<<endl; break;} |
| 1cca57bf |
233 | if (part->GetStatusCode()>0) { |
| 234 | cerr<<"Primary momentum: "<<part->P()<<endl; break; |
| 235 | } |
| 7f4044f0 |
236 | label=m; |
| 237 | part=(TParticle*)gAlice->Particle(label); |
| 238 | } |
| 239 | if (lab[1]<0) lab[1]=label; |
| 240 | else if (lab[2]<0) lab[2]=label; |
| 241 | else ;//cerr<<"CheckLabels : No empty labels !\n"; |
| 242 | } |
| 243 | } |
| 244 | |
| 1cca57bf |
245 | void AliITSclustererV2::RecPoints2Clusters |
| 246 | (const TClonesArray *points, Int_t idx, TClonesArray *clusters) { |
| 7f4044f0 |
247 | //------------------------------------------------------------ |
| 1cca57bf |
248 | // Conversion AliITSRecPoint -> AliITSclusterV2 for the ITS |
| 249 | // subdetector indexed by idx |
| 7f4044f0 |
250 | //------------------------------------------------------------ |
| 251 | TClonesArray &cl=*clusters; |
| 252 | Int_t ncl=points->GetEntriesFast(); |
| 253 | for (Int_t i=0; i<ncl; i++) { |
| 254 | AliITSRecPoint *p = (AliITSRecPoint *)points->UncheckedAt(i); |
| 255 | Float_t lp[5]; |
| 1cca57bf |
256 | lp[0]=-p->GetX()-fYshift[idx]; if (idx<=fLastSPD1) lp[0]*=-1; //SPD1 |
| 257 | lp[1]=p->GetZ()+fZshift[idx]; |
| 258 | lp[2]=p->GetSigmaX2(); |
| 259 | lp[3]=p->GetSigmaZ2(); |
| 260 | lp[4]=p->GetQ(); |
| 7f4044f0 |
261 | Int_t lab[4]; |
| 262 | lab[0]=p->GetLabel(0); lab[1]=p->GetLabel(1); lab[2]=p->GetLabel(2); |
| 1cca57bf |
263 | lab[3]=fNdet[idx]; |
| 7f4044f0 |
264 | CheckLabels(lab); |
| 265 | new (cl[i]) AliITSclusterV2(lab,lp); |
| 266 | } |
| 267 | } |
| 268 | |
| 269 | void AliITSclustererV2::Hits2Clusters(const TFile *in, TFile *out) { |
| 270 | //------------------------------------------------------------ |
| 271 | // This function creates ITS clusters |
| 272 | //------------------------------------------------------------ |
| 273 | TDirectory *savedir=gDirectory; |
| 274 | |
| 275 | if (!out->IsOpen()) { |
| 276 | cerr<<"AliITSclustererV2::Hits2Clusters: output file not open !\n"; |
| 277 | return; |
| 278 | } |
| 279 | |
| 280 | if (!gAlice) { |
| 281 | cerr<<"AliITSclustererV2::Hits2Clusters : gAlice==0 !\n"; |
| 282 | return; |
| 283 | } |
| 284 | |
| 285 | AliITS *its = (AliITS*)gAlice->GetModule("ITS"); |
| 286 | if (!its) { |
| 287 | cerr<<"AliITSclustererV2::Hits2Clusters : Can't find the ITS !\n"; |
| 288 | return; |
| 289 | } |
| 290 | AliITSgeom *geom=its->GetITSgeom(); |
| 291 | Int_t mmax=geom->GetIndexMax(); |
| 292 | |
| 293 | its->InitModules(-1,mmax); |
| 294 | its->FillModules(gAlice->TreeH(),0); |
| 295 | |
| 296 | out->cd(); |
| 297 | |
| 298 | Char_t name[100]; |
| 299 | sprintf(name,"TreeC_ITS_%d",fEvent); |
| 300 | TTree cTree(name,"ITS clusters V2"); |
| 301 | TClonesArray *clusters=new TClonesArray("AliITSclusterV2",1000); |
| 302 | cTree.Branch("Clusters",&clusters); |
| 303 | |
| 304 | static TClonesArray *points=its->RecPoints(); |
| 305 | AliITSsimulationFastPoints sim; |
| 306 | Int_t ncl=0; |
| 307 | for (Int_t m=0; m<mmax; m++) { |
| 308 | AliITSmodule *mod=its->GetModule(m); |
| 309 | sim.CreateFastRecPoints(mod,m,gRandom); |
| 310 | |
| 1cca57bf |
311 | RecPoints2Clusters(points, m, clusters); |
| 7f4044f0 |
312 | its->ResetRecPoints(); |
| 313 | |
| 314 | ncl+=clusters->GetEntriesFast(); |
| 315 | cTree.Fill(); |
| 316 | clusters->Clear(); |
| 317 | } |
| 318 | cTree.Write(); |
| 319 | |
| 320 | cerr<<"Number of found fast clusters : "<<ncl<<endl; |
| 321 | |
| 322 | delete clusters; |
| 323 | |
| 324 | savedir->cd(); |
| 325 | } |
| 326 | |
| 327 | //*********************************** |
| 328 | |
| 329 | #ifndef V1 |
| 330 | |
| 331 | void AliITSclustererV2:: |
| 332 | FindCluster(Int_t k,Int_t maxz,AliBin *bins,Int_t &n,Int_t *idx) { |
| 333 | //------------------------------------------------------------ |
| 334 | // returns an array of indices of digits belonging to the cluster |
| 335 | // (needed when the segmentation is not regular) |
| 336 | //------------------------------------------------------------ |
| 337 | if (n<200) idx[n++]=bins[k].GetIndex(); |
| 338 | bins[k].Use(); |
| 339 | |
| 340 | if (bins[k-maxz].IsNotUsed()) FindCluster(k-maxz,maxz,bins,n,idx); |
| 341 | if (bins[k-1 ].IsNotUsed()) FindCluster(k-1 ,maxz,bins,n,idx); |
| 342 | if (bins[k+maxz].IsNotUsed()) FindCluster(k+maxz,maxz,bins,n,idx); |
| 343 | if (bins[k+1 ].IsNotUsed()) FindCluster(k+1 ,maxz,bins,n,idx); |
| 344 | /* |
| 345 | if (bins[k-maxz-1].IsNotUsed()) FindCluster(k-maxz-1,maxz,bins,n,idx); |
| 346 | if (bins[k-maxz+1].IsNotUsed()) FindCluster(k-maxz+1,maxz,bins,n,idx); |
| 347 | if (bins[k+maxz-1].IsNotUsed()) FindCluster(k+maxz-1,maxz,bins,n,idx); |
| 348 | if (bins[k+maxz+1].IsNotUsed()) FindCluster(k+maxz+1,maxz,bins,n,idx); |
| 349 | */ |
| 350 | } |
| 351 | |
| 352 | void AliITSclustererV2:: |
| 353 | FindClustersSPD(const TClonesArray *digits, TClonesArray *clusters) { |
| 354 | //------------------------------------------------------------ |
| 355 | // Actual SPD cluster finder |
| 356 | //------------------------------------------------------------ |
| c391f9d9 |
357 | Int_t kNzBins = fNzSPD + 2; |
| 358 | const Int_t kMAXBIN=kNzBins*(fNySPD+2); |
| 7f4044f0 |
359 | |
| 360 | Int_t ndigits=digits->GetEntriesFast(); |
| 361 | AliBin *bins=new AliBin[kMAXBIN]; |
| 362 | |
| 363 | Int_t k; |
| 364 | AliITSdigitSPD *d=0; |
| 365 | for (k=0; k<ndigits; k++) { |
| 366 | d=(AliITSdigitSPD*)digits->UncheckedAt(k); |
| 367 | Int_t i=d->GetCoord2()+1; //y |
| 368 | Int_t j=d->GetCoord1()+1; |
| c391f9d9 |
369 | bins[i*kNzBins+j].SetIndex(k); |
| 370 | bins[i*kNzBins+j].SetMask(1); |
| 7f4044f0 |
371 | } |
| 372 | |
| 373 | Int_t n=0; TClonesArray &cl=*clusters; |
| 374 | for (k=0; k<kMAXBIN; k++) { |
| 375 | if (!bins[k].IsNotUsed()) continue; |
| 376 | Int_t ni=0, idx[200]; |
| c391f9d9 |
377 | FindCluster(k,kNzBins,bins,ni,idx); |
| 7f4044f0 |
378 | if (ni==200) {cerr<<"SPD: Too big cluster !\n"; continue;} |
| 379 | |
| 380 | Int_t lab[4]; |
| 381 | lab[0]=-2; |
| 382 | lab[1]=-2; |
| 383 | lab[2]=-2; |
| 384 | lab[3]=fNdet[fI]; |
| 385 | |
| 386 | d=(AliITSdigitSPD*)digits->UncheckedAt(idx[0]); |
| 387 | Int_t ymin=d->GetCoord2(),ymax=ymin; |
| 388 | Int_t zmin=d->GetCoord1(),zmax=zmin; |
| 389 | Float_t y=0.,z=0.,q=0.; |
| 390 | for (Int_t l=0; l<ni; l++) { |
| 391 | d=(AliITSdigitSPD*)digits->UncheckedAt(idx[l]); |
| 392 | |
| 393 | if (ymin > d->GetCoord2()) ymin=d->GetCoord2(); |
| 394 | if (ymax < d->GetCoord2()) ymax=d->GetCoord2(); |
| 395 | if (zmin > d->GetCoord1()) zmin=d->GetCoord1(); |
| 396 | if (zmax < d->GetCoord1()) zmax=d->GetCoord1(); |
| 397 | |
| 398 | Int_t lab0=(d->GetTracks())[0]; |
| 399 | if (lab0>=0) { |
| 400 | if (lab[0]<0) { |
| 401 | lab[0]=lab0; |
| 402 | } else if (lab[1]<0) { |
| 403 | if (lab0!=lab[0]) lab[1]=lab0; |
| 404 | } else if (lab[2]<0) { |
| 405 | if (lab0!=lab[0]) |
| 406 | if (lab0!=lab[1]) lab[2]=lab0; |
| 407 | } |
| 408 | } |
| 409 | Float_t qq=d->GetSignal(); |
| 410 | y+=qq*fYSPD[d->GetCoord2()]; z+=qq*fZSPD[d->GetCoord1()]; q+=qq; |
| 411 | } |
| 412 | y/=q; z/=q; |
| 413 | y-=fHwSPD; z-=fHlSPD; |
| 414 | |
| 415 | Float_t lp[5]; |
| 416 | lp[0]=-y-fYshift[fI]; if (fI<=fLastSPD1) lp[0]=-lp[0]; |
| 417 | lp[1]= z+fZshift[fI]; |
| 418 | lp[2]= fYpitchSPD*fYpitchSPD/12.; |
| 419 | lp[3]= fZ1pitchSPD*fZ1pitchSPD/12.; |
| 420 | //lp[4]= q; |
| 421 | lp[4]= (zmax-zmin+1)*100 + (ymax-ymin+1); |
| 422 | |
| 423 | //CheckLabels(lab); |
| c391f9d9 |
424 | d=(AliITSdigitSPD*)digits->UncheckedAt(idx[0]); |
| 7f4044f0 |
425 | new (cl[n]) AliITSclusterV2(lab,lp); n++; |
| 426 | } |
| 427 | |
| 428 | delete bins; |
| 429 | } |
| 430 | |
| c391f9d9 |
431 | void AliITSclustererV2::FindClustersSPD(AliITSRawStream* input, |
| 432 | TClonesArray** clusters) |
| 433 | { |
| 434 | //------------------------------------------------------------ |
| 435 | // Actual SPD cluster finder for raw data |
| 436 | //------------------------------------------------------------ |
| 437 | |
| 438 | Int_t nClustersSPD = 0; |
| 439 | Int_t kNzBins = fNzSPD + 2; |
| 440 | Int_t kNyBins = fNySPD + 2; |
| 441 | Int_t kMaxBin = kNzBins * kNyBins; |
| 442 | AliBin* bins = NULL; |
| 443 | |
| 444 | // read raw data input stream |
| 445 | while (kTRUE) { |
| 446 | Bool_t next = input->Next(); |
| 447 | if (!next || input->IsNewModule()) { |
| 448 | Int_t iModule = input->GetPrevModuleID(); |
| 449 | |
| 450 | // when all data from a module was read, search for clusters |
| 451 | if (bins) { |
| 452 | clusters[iModule] = new TClonesArray("AliITSclusterV2"); |
| 453 | Int_t nClusters = 0; |
| 454 | |
| 455 | for (Int_t iBin = 0; iBin < kMaxBin; iBin++) { |
| 456 | if (bins[iBin].IsUsed()) continue; |
| 457 | Int_t nBins = 0; |
| 458 | Int_t idxBins[200]; |
| 459 | FindCluster(iBin, kNzBins, bins, nBins, idxBins); |
| 460 | if (nBins == 200) { |
| 461 | Error("FindClustersSPD", "SPD: Too big cluster !\n"); |
| 462 | continue; |
| 463 | } |
| 464 | |
| 465 | Int_t label[4]; |
| 466 | label[0] = -2; |
| 467 | label[1] = -2; |
| 468 | label[2] = -2; |
| 469 | // label[3] = iModule; |
| 470 | label[3] = fNdet[iModule]; |
| 471 | |
| 472 | Int_t ymin = (idxBins[0] / kNzBins) - 1; |
| 473 | Int_t ymax = ymin; |
| 474 | Int_t zmin = (idxBins[0] % kNzBins) - 1; |
| 475 | Int_t zmax = zmin; |
| 476 | Float_t y = 0.; |
| 477 | Float_t z = 0.; |
| 478 | Float_t q = 0.; |
| 479 | for (Int_t idx = 0; idx < nBins; idx++) { |
| 480 | Int_t iy = (idxBins[idx] / kNzBins) - 1; |
| 481 | Int_t iz = (idxBins[idx] % kNzBins) - 1; |
| 482 | if (ymin > iy) ymin = iy; |
| 483 | if (ymax < iy) ymax = iy; |
| 484 | if (zmin > iz) zmin = iz; |
| 485 | if (zmax < iz) zmax = iz; |
| 486 | |
| 487 | Float_t qBin = bins[idxBins[idx]].GetQ(); |
| 488 | y += qBin * fYSPD[iy]; |
| 489 | z += qBin * fZSPD[iz]; |
| 490 | q += qBin; |
| 491 | } |
| 492 | y /= q; |
| 493 | z /= q; |
| 494 | y -= fHwSPD; |
| 495 | z -= fHlSPD; |
| 496 | |
| 497 | Float_t hit[5]; // y, z, sigma(y)^2, sigma(z)^2, charge |
| 498 | hit[0] = -y-fYshift[iModule]; |
| 499 | if (iModule <= fLastSPD1) hit[0] = -hit[0]; |
| 500 | hit[1] = z+fZshift[iModule]; |
| 501 | hit[2] = fYpitchSPD*fYpitchSPD/12.; |
| 502 | hit[3] = fZ1pitchSPD*fZ1pitchSPD/12.; |
| 503 | // hit[4] = q; |
| 504 | hit[4] = (zmax-zmin+1)*100 + (ymax-ymin+1); |
| 505 | |
| 506 | //CheckLabels(label); |
| 507 | new (clusters[iModule]->AddrAt(nClusters)) |
| 508 | AliITSclusterV2(label, hit); |
| 509 | nClusters++; |
| 510 | } |
| 511 | |
| 512 | nClustersSPD += nClusters; |
| 513 | delete bins; |
| 514 | } |
| 515 | |
| 516 | if (!next) break; |
| 517 | bins = new AliBin[kMaxBin]; |
| 518 | } |
| 519 | |
| 520 | // fill the current digit into the bins array |
| 521 | Int_t index = (input->GetCoord2()+1) * kNzBins + (input->GetCoord1()+1); |
| 522 | bins[index].SetIndex(index); |
| 523 | bins[index].SetMask(1); |
| 524 | bins[index].SetQ(1); |
| 525 | } |
| 526 | |
| 527 | Info("FindClustersSPD", "found clusters in ITS SPD: %d", nClustersSPD); |
| 528 | } |
| 529 | |
| 530 | |
| 7f4044f0 |
531 | Bool_t AliITSclustererV2::IsMaximum(Int_t k,Int_t max,const AliBin *bins) { |
| 532 | //------------------------------------------------------------ |
| 533 | //is this a local maximum ? |
| 534 | //------------------------------------------------------------ |
| 535 | UShort_t q=bins[k].GetQ(); |
| 536 | if (q==1023) return kFALSE; |
| 537 | if (bins[k-max].GetQ() > q) return kFALSE; |
| 538 | if (bins[k-1 ].GetQ() > q) return kFALSE; |
| 539 | if (bins[k+max].GetQ() > q) return kFALSE; |
| 540 | if (bins[k+1 ].GetQ() > q) return kFALSE; |
| 541 | if (bins[k-max-1].GetQ() > q) return kFALSE; |
| 542 | if (bins[k+max-1].GetQ() > q) return kFALSE; |
| 543 | if (bins[k+max+1].GetQ() > q) return kFALSE; |
| 544 | if (bins[k-max+1].GetQ() > q) return kFALSE; |
| 545 | return kTRUE; |
| 546 | } |
| 547 | |
| 548 | void AliITSclustererV2:: |
| 549 | FindPeaks(Int_t k,Int_t max,AliBin *b,Int_t *idx,UInt_t *msk,Int_t& n) { |
| 550 | //------------------------------------------------------------ |
| 551 | //find local maxima |
| 552 | //------------------------------------------------------------ |
| 553 | if (n<31) |
| 554 | if (IsMaximum(k,max,b)) { |
| 555 | idx[n]=k; msk[n]=(2<<n); |
| 556 | n++; |
| 557 | } |
| 558 | b[k].SetMask(0); |
| 559 | if (b[k-max].GetMask()&1) FindPeaks(k-max,max,b,idx,msk,n); |
| 560 | if (b[k-1 ].GetMask()&1) FindPeaks(k-1 ,max,b,idx,msk,n); |
| 561 | if (b[k+max].GetMask()&1) FindPeaks(k+max,max,b,idx,msk,n); |
| 562 | if (b[k+1 ].GetMask()&1) FindPeaks(k+1 ,max,b,idx,msk,n); |
| 563 | } |
| 564 | |
| 565 | void AliITSclustererV2:: |
| 566 | MarkPeak(Int_t k, Int_t max, AliBin *bins, UInt_t m) { |
| 567 | //------------------------------------------------------------ |
| 568 | //mark this peak |
| 569 | //------------------------------------------------------------ |
| 570 | UShort_t q=bins[k].GetQ(); |
| 571 | |
| 572 | bins[k].SetMask(bins[k].GetMask()|m); |
| 573 | |
| 574 | if (bins[k-max].GetQ() <= q) |
| 575 | if ((bins[k-max].GetMask()&m) == 0) MarkPeak(k-max,max,bins,m); |
| 576 | if (bins[k-1 ].GetQ() <= q) |
| 577 | if ((bins[k-1 ].GetMask()&m) == 0) MarkPeak(k-1 ,max,bins,m); |
| 578 | if (bins[k+max].GetQ() <= q) |
| 579 | if ((bins[k+max].GetMask()&m) == 0) MarkPeak(k+max,max,bins,m); |
| 580 | if (bins[k+1 ].GetQ() <= q) |
| 581 | if ((bins[k+1 ].GetMask()&m) == 0) MarkPeak(k+1 ,max,bins,m); |
| 582 | } |
| 583 | |
| 584 | void AliITSclustererV2:: |
| 585 | MakeCluster(Int_t k,Int_t max,AliBin *bins,UInt_t m,AliITSclusterV2 &c) { |
| 586 | //------------------------------------------------------------ |
| 587 | //make cluster using digits of this peak |
| 588 | //------------------------------------------------------------ |
| 589 | Float_t q=(Float_t)bins[k].GetQ(); |
| 590 | Int_t i=k/max, j=k-i*max; |
| 591 | |
| 592 | c.SetQ(c.GetQ()+q); |
| 593 | c.SetY(c.GetY()+i*q); |
| 594 | c.SetZ(c.GetZ()+j*q); |
| 595 | c.SetSigmaY2(c.GetSigmaY2()+i*i*q); |
| 596 | c.SetSigmaZ2(c.GetSigmaZ2()+j*j*q); |
| 597 | |
| 598 | bins[k].SetMask(0xFFFFFFFE); |
| 599 | |
| 600 | if (bins[k-max].GetMask() == m) MakeCluster(k-max,max,bins,m,c); |
| 601 | if (bins[k-1 ].GetMask() == m) MakeCluster(k-1 ,max,bins,m,c); |
| 602 | if (bins[k+max].GetMask() == m) MakeCluster(k+max,max,bins,m,c); |
| 603 | if (bins[k+1 ].GetMask() == m) MakeCluster(k+1 ,max,bins,m,c); |
| 604 | } |
| 605 | |
| 606 | void AliITSclustererV2:: |
| c391f9d9 |
607 | FindClustersSDD(AliBin* bins[2], Int_t nMaxBin, Int_t nzBins, |
| 608 | const TClonesArray *digits, TClonesArray *clusters) { |
| 7f4044f0 |
609 | //------------------------------------------------------------ |
| 610 | // Actual SDD cluster finder |
| 611 | //------------------------------------------------------------ |
| 7f4044f0 |
612 | Int_t ncl=0; TClonesArray &cl=*clusters; |
| 613 | for (Int_t s=0; s<2; s++) |
| c391f9d9 |
614 | for (Int_t i=0; i<nMaxBin; i++) { |
| 7f4044f0 |
615 | if (bins[s][i].IsUsed()) continue; |
| 616 | Int_t idx[32]; UInt_t msk[32]; Int_t npeaks=0; |
| c391f9d9 |
617 | FindPeaks(i, nzBins, bins[s], idx, msk, npeaks); |
| 7f4044f0 |
618 | |
| 619 | if (npeaks>30) continue; |
| 620 | |
| 621 | Int_t k,l; |
| 622 | for (k=0; k<npeaks-1; k++){//mark adjacent peaks |
| 623 | if (idx[k] < 0) continue; //this peak is already removed |
| 624 | for (l=k+1; l<npeaks; l++) { |
| 625 | if (idx[l] < 0) continue; //this peak is already removed |
| c391f9d9 |
626 | Int_t ki=idx[k]/nzBins, kj=idx[k] - ki*nzBins; |
| 627 | Int_t li=idx[l]/nzBins, lj=idx[l] - li*nzBins; |
| 7f4044f0 |
628 | Int_t di=TMath::Abs(ki - li); |
| 629 | Int_t dj=TMath::Abs(kj - lj); |
| 630 | if (di>1 || dj>1) continue; |
| 631 | if (bins[s][idx[k]].GetQ() > bins[s][idx[l]].GetQ()) { |
| 632 | msk[l]=msk[k]; |
| 633 | idx[l]*=-1; |
| 634 | } else { |
| 635 | msk[k]=msk[l]; |
| 636 | idx[k]*=-1; |
| 637 | break; |
| 638 | } |
| 639 | } |
| 640 | } |
| 641 | |
| 642 | for (k=0; k<npeaks; k++) { |
| c391f9d9 |
643 | MarkPeak(TMath::Abs(idx[k]), nzBins, bins[s], msk[k]); |
| 7f4044f0 |
644 | } |
| 645 | |
| 646 | for (k=0; k<npeaks; k++) { |
| 647 | if (idx[k] < 0) continue; //removed peak |
| 648 | AliITSclusterV2 c; |
| c391f9d9 |
649 | MakeCluster(idx[k], nzBins, bins[s], msk[k], c); |
| 7f4044f0 |
650 | |
| 651 | //if (c.GetQ() < 200) continue; //noise cluster |
| 652 | |
| 653 | /* |
| 654 | Float_t s2 = c.GetSigmaY2()/c.GetQ() - c.GetY()*c.GetY(); |
| 655 | Float_t w=par->GetPadPitchWidth(sec); |
| 656 | c.SetSigmaY2((s2 + 1./12.)*w*w); |
| 657 | if (s2 != 0.) { |
| 658 | c.SetSigmaY2(c.GetSigmaY2()*0.108); |
| 659 | if (sec<par->GetNInnerSector()) c.SetSigmaY2(c.GetSigmaY2()*2.07); |
| 660 | } |
| 661 | |
| 662 | s2 = c.GetSigmaZ2()/c.GetQ() - c.GetZ()*c.GetZ(); |
| 663 | w=par->GetZWidth(); |
| 664 | c.SetSigmaZ2((s2 + 1./12.)*w*w); |
| 665 | if (s2 != 0.) { |
| 666 | c.SetSigmaZ2(c.GetSigmaZ2()*0.169); |
| 667 | if (sec<par->GetNInnerSector()) c.SetSigmaZ2(c.GetSigmaZ2()*1.77); |
| 668 | } |
| 669 | */ |
| 670 | |
| 671 | c.SetSigmaY2(0.0030*0.0030); |
| 672 | c.SetSigmaZ2(0.0020*0.0020); |
| 673 | c.SetDetectorIndex(fNdet[fI]); |
| 674 | |
| 675 | Float_t y=c.GetY(),z=c.GetZ(), q=c.GetQ(); |
| 676 | y/=q; z/=q; |
| 677 | |
| 678 | y=(y-0.5)*fYpitchSDD; |
| 679 | y-=fHwSDD; |
| 680 | y-=fYoffSDD; //delay ? |
| 681 | if (s) y=-y; |
| 682 | |
| 683 | z=(z-0.5)*fZpitchSDD; |
| 684 | z-=fHlSDD; |
| 685 | |
| 686 | y=-y-fYshift[fI]; |
| 687 | z= z+fZshift[fI]; |
| 688 | c.SetY(y); |
| 689 | c.SetZ(z); |
| 690 | |
| 691 | c.SetQ(q/20.); //to be consistent with the SSD charges |
| 692 | |
| c391f9d9 |
693 | if (digits) { |
| 694 | AliBin *b=&bins[s][idx[k]]; |
| 695 | AliITSdigitSDD* d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); |
| 696 | Int_t l0=(d->GetTracks())[0]; |
| 697 | if (l0<0) { |
| 698 | b=&bins[s][idx[k]-1]; |
| 699 | if (b->GetQ()>0) { |
| 700 | d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); |
| 701 | l0=(d->GetTracks())[0]; |
| 702 | } |
| 703 | } |
| 704 | if (l0<0) { |
| 705 | b=&bins[s][idx[k]+1]; |
| 706 | if (b->GetQ()>0) { |
| 707 | d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); |
| 708 | l0=(d->GetTracks())[0]; |
| 709 | } |
| 710 | } |
| 711 | if (l0<0) { |
| 712 | b=&bins[s][idx[k]-nzBins]; |
| 713 | if (b->GetQ()>0) { |
| 714 | d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); |
| 715 | l0=(d->GetTracks())[0]; |
| 716 | } |
| 717 | } |
| 718 | if (l0<0) { |
| 719 | b=&bins[s][idx[k]+nzBins]; |
| 720 | if (b->GetQ()>0) { |
| 721 | d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); |
| 722 | l0=(d->GetTracks())[0]; |
| 723 | } |
| 7f4044f0 |
724 | } |
| 7f4044f0 |
725 | |
| c391f9d9 |
726 | if (l0<0) { |
| 727 | b=&bins[s][idx[k]+nzBins+1]; |
| 728 | if (b->GetQ()>0) { |
| 729 | d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); |
| 730 | l0=(d->GetTracks())[0]; |
| 731 | } |
| 732 | } |
| 733 | if (l0<0) { |
| 734 | b=&bins[s][idx[k]+nzBins-1]; |
| 735 | if (b->GetQ()>0) { |
| 736 | d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); |
| 737 | l0=(d->GetTracks())[0]; |
| 738 | } |
| 739 | } |
| 740 | if (l0<0) { |
| 741 | b=&bins[s][idx[k]-nzBins+1]; |
| 742 | if (b->GetQ()>0) { |
| 743 | d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); |
| 744 | l0=(d->GetTracks())[0]; |
| 745 | } |
| 746 | } |
| 747 | if (l0<0) { |
| 748 | b=&bins[s][idx[k]-nzBins-1]; |
| 749 | if (b->GetQ()>0) { |
| 750 | d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); |
| 751 | l0=(d->GetTracks())[0]; |
| 752 | } |
| 753 | } |
| 7f4044f0 |
754 | |
| c391f9d9 |
755 | { |
| 756 | Int_t lab[3]; |
| 757 | lab[0]=(d->GetTracks())[0]; |
| 758 | lab[1]=(d->GetTracks())[1]; |
| 759 | lab[2]=(d->GetTracks())[2]; |
| 760 | //CheckLabels(lab); |
| 761 | c.SetLabel(lab[0],0); |
| 762 | c.SetLabel(lab[1],1); |
| 763 | c.SetLabel(lab[2],2); |
| 764 | } |
| 765 | } |
| 7f4044f0 |
766 | |
| 767 | new (cl[ncl]) AliITSclusterV2(c); ncl++; |
| 768 | } |
| 769 | } |
| c391f9d9 |
770 | } |
| 771 | |
| 772 | void AliITSclustererV2:: |
| 773 | FindClustersSDD(const TClonesArray *digits, TClonesArray *clusters) { |
| 774 | //------------------------------------------------------------ |
| 775 | // Actual SDD cluster finder |
| 776 | //------------------------------------------------------------ |
| 777 | Int_t kNzBins = fNzSDD + 2; |
| 778 | const Int_t kMAXBIN=kNzBins*(fNySDD+2); |
| 779 | |
| 780 | AliBin *bins[2]; |
| 781 | bins[0]=new AliBin[kMAXBIN]; |
| 782 | bins[1]=new AliBin[kMAXBIN]; |
| 783 | |
| 784 | AliITSdigitSDD *d=0; |
| 785 | Int_t i, ndigits=digits->GetEntriesFast(); |
| 786 | for (i=0; i<ndigits; i++) { |
| 787 | d=(AliITSdigitSDD*)digits->UncheckedAt(i); |
| 788 | Int_t y=d->GetCoord2()+1; //y |
| 789 | Int_t z=d->GetCoord1()+1; //z |
| 790 | Int_t q=d->GetSignal(); |
| 791 | if (z <= fNzSDD) { |
| 792 | bins[0][y*kNzBins+z].SetQ(q); |
| 793 | bins[0][y*kNzBins+z].SetMask(1); |
| 794 | bins[0][y*kNzBins+z].SetIndex(i); |
| 795 | } else { |
| 796 | z-=fNzSDD; |
| 797 | bins[1][y*kNzBins+z].SetQ(q); |
| 798 | bins[1][y*kNzBins+z].SetMask(1); |
| 799 | bins[1][y*kNzBins+z].SetIndex(i); |
| 800 | } |
| 801 | } |
| 802 | |
| 803 | FindClustersSDD(bins, kMAXBIN, kNzBins, digits, clusters); |
| 7f4044f0 |
804 | |
| 805 | delete[] bins[0]; |
| 806 | delete[] bins[1]; |
| 807 | } |
| 808 | |
| c391f9d9 |
809 | void AliITSclustererV2::FindClustersSDD(AliITSRawStream* input, |
| 810 | TClonesArray** clusters) |
| 811 | { |
| 812 | //------------------------------------------------------------ |
| 813 | // Actual SDD cluster finder for raw data |
| 814 | //------------------------------------------------------------ |
| 815 | Int_t nClustersSDD = 0; |
| 816 | Int_t kNzBins = fNzSDD + 2; |
| 817 | Int_t kMaxBin = kNzBins * (fNySDD+2); |
| 818 | AliBin* bins[2] = {NULL, NULL}; |
| 819 | |
| 820 | // read raw data input stream |
| 821 | while (kTRUE) { |
| 822 | Bool_t next = input->Next(); |
| 823 | if (!next || input->IsNewModule()) { |
| 824 | Int_t iModule = input->GetPrevModuleID(); |
| 825 | |
| 826 | // when all data from a module was read, search for clusters |
| 827 | if (bins[0]) { |
| 828 | clusters[iModule] = new TClonesArray("AliITSclusterV2"); |
| 829 | fI = iModule; |
| 830 | FindClustersSDD(bins, kMaxBin, kNzBins, NULL, clusters[iModule]); |
| 831 | Int_t nClusters = clusters[iModule]->GetEntriesFast(); |
| 832 | nClustersSDD += nClusters; |
| 833 | delete[] bins[0]; |
| 834 | delete[] bins[1]; |
| 835 | } |
| 836 | |
| 837 | if (!next) break; |
| 838 | bins[0] = new AliBin[kMaxBin]; |
| 839 | bins[1] = new AliBin[kMaxBin]; |
| 840 | } |
| 841 | |
| 842 | // fill the current digit into the bins array |
| 843 | Int_t iz = input->GetCoord1()+1; |
| 844 | Int_t side = ((iz <= fNzSDD) ? 0 : 1); |
| 845 | iz -= side*fNzSDD; |
| 846 | Int_t index = (input->GetCoord2()+1) * kNzBins + iz; |
| 847 | bins[side][index].SetQ(input->GetSignal()); |
| 848 | bins[side][index].SetMask(1); |
| 849 | bins[side][index].SetIndex(index); |
| 850 | } |
| 851 | |
| 852 | Info("FindClustersSDD", "found clusters in ITS SDD: %d", nClustersSDD); |
| 853 | } |
| 854 | |
| 855 | void AliITSclustererV2:: |
| 856 | FindClustersSSD(Ali1Dcluster* neg, Int_t nn, |
| 857 | Ali1Dcluster* pos, Int_t np, |
| 858 | TClonesArray *clusters) { |
| 859 | //------------------------------------------------------------ |
| 860 | // Actual SSD cluster finder |
| 861 | //------------------------------------------------------------ |
| 862 | TClonesArray &cl=*clusters; |
| 863 | |
| 864 | Int_t lab[4]={-2,-2,-2,-2}; |
| 865 | Float_t tanp=fTanP, tann=fTanN; |
| 866 | if (fI>fLastSSD1) {tann=fTanP; tanp=fTanN;} |
| 867 | |
| 868 | Int_t idet=fNdet[fI]; |
| 869 | Int_t ncl=0; |
| 870 | for (Int_t i=0; i<np; i++) { |
| 871 | //Float_t dq_min=1.e+33; |
| 872 | Float_t ybest=1000,zbest=1000,qbest=0; |
| 873 | Float_t yp=pos[i].GetY()*fYpitchSSD; |
| 874 | for (Int_t j=0; j<nn; j++) { |
| 875 | //if (pos[i].fTracks[0] != neg[j].fTracks[0]) continue; |
| 876 | Float_t yn=neg[j].GetY()*fYpitchSSD; |
| 877 | Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp); |
| 878 | Float_t yt=yn + tann*zt; |
| 879 | zt-=fHlSSD; yt-=fHwSSD; |
| 880 | if (TMath::Abs(yt)<fHwSSD+0.01) |
| 881 | if (TMath::Abs(zt)<fHlSSD+0.01) { |
| 882 | //if (TMath::Abs(pos[i].GetQ()-neg[j].GetQ())<dq_min) { |
| 883 | //dq_min=TMath::Abs(pos[i].GetQ()-neg[j].GetQ()); |
| 884 | ybest=yt; zbest=zt; |
| 885 | qbest=0.5*(pos[i].GetQ()+neg[j].GetQ()); |
| 886 | |
| 887 | lab[0]=pos[i].GetLabel(0); |
| 888 | lab[1]=pos[i].GetLabel(1); |
| 889 | lab[2]=neg[i].GetLabel(0); |
| 890 | lab[3]=(((i<<10) + j)<<10) + idet; // pos|neg|det |
| 891 | Float_t lp[5]; |
| 892 | lp[0]=-ybest-fYshift[fI]; |
| 893 | lp[1]= zbest+fZshift[fI]; |
| 894 | lp[2]=0.0025*0.0025; //SigmaY2 |
| 895 | lp[3]=0.110*0.110; //SigmaZ2 |
| 896 | if (pos[i].GetNd()+neg[j].GetNd() > 4) { |
| 897 | lp[2]*=9; |
| 898 | lp[3]*=9; |
| 899 | } |
| 900 | lp[4]=qbest; //Q |
| 901 | |
| 902 | //CheckLabels(lab); |
| 903 | new (cl[ncl]) AliITSclusterV2(lab,lp); ncl++; |
| 904 | } |
| 905 | } |
| 906 | /* |
| 907 | if (ybest<100) { |
| 908 | lab[3]=idet; |
| 909 | Float_t lp[5]; |
| 910 | lp[0]=-ybest-fYshift[fI]; |
| 911 | lp[1]= zbest+fZshift[fI]; |
| 912 | lp[2]=0.002*0.002; //SigmaY2 |
| 913 | lp[3]=0.080*0.080; //SigmaZ2 |
| 914 | lp[4]=qbest; //Q |
| 915 | // |
| 916 | new (cl[ncl]) AliITSclusterV2(lab,lp); ncl++; |
| 917 | } |
| 918 | */ |
| 919 | } |
| 920 | } |
| 921 | |
| 7f4044f0 |
922 | void AliITSclustererV2:: |
| 923 | FindClustersSSD(const TClonesArray *digits, TClonesArray *clusters) { |
| 924 | //------------------------------------------------------------ |
| 925 | // Actual SSD cluster finder |
| 926 | //------------------------------------------------------------ |
| 927 | Int_t smax=digits->GetEntriesFast(); |
| 928 | if (smax==0) return; |
| 929 | |
| 930 | const Int_t MAX=1000; |
| 931 | Int_t np=0, nn=0; |
| 932 | Ali1Dcluster pos[MAX], neg[MAX]; |
| 933 | Float_t y=0., q=0., qmax=0.; |
| 934 | Int_t lab[4]={-2,-2,-2,-2}; |
| 935 | |
| 7f4044f0 |
936 | AliITSdigitSSD *d=(AliITSdigitSSD*)digits->UncheckedAt(0); |
| 937 | q += d->GetSignal(); |
| 938 | y += d->GetCoord2()*d->GetSignal(); |
| 939 | qmax=d->GetSignal(); |
| 940 | lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2); |
| 941 | Int_t curr=d->GetCoord2(); |
| 942 | Int_t flag=d->GetCoord1(); |
| 943 | Int_t *n=&nn; |
| 944 | Ali1Dcluster *c=neg; |
| c391f9d9 |
945 | Int_t nd=1; |
| 7f4044f0 |
946 | for (Int_t s=1; s<smax; s++) { |
| 947 | d=(AliITSdigitSSD*)digits->UncheckedAt(s); |
| 948 | Int_t strip=d->GetCoord2(); |
| 949 | if ((strip-curr) > 1 || flag!=d->GetCoord1()) { |
| 950 | c[*n].SetY(y/q); |
| 951 | c[*n].SetQ(q); |
| 952 | c[*n].SetNd(nd); |
| 953 | c[*n].SetLabels(lab); |
| 954 | //Split suspiciously big cluster |
| 955 | if (nd>3) { |
| 956 | c[*n].SetY(y/q-0.5*nd); |
| 957 | c[*n].SetQ(0.5*q); |
| 958 | (*n)++; |
| 959 | if (*n==MAX) { |
| 960 | cerr<< |
| 961 | "AliITSclustererV2::FindClustersSSD: Too many 1D clusters !\n"; |
| 962 | return; |
| 963 | } |
| 964 | c[*n].SetY(y/q+0.5*nd); |
| 965 | c[*n].SetQ(0.5*q); |
| 966 | c[*n].SetNd(nd); |
| 967 | c[*n].SetLabels(lab); |
| 968 | } |
| 969 | (*n)++; |
| 970 | if (*n==MAX) { |
| 971 | cerr<<"AliITSclustererV2::FindClustersSSD: Too many 1D clusters !\n"; |
| 972 | return; |
| 973 | } |
| 974 | y=q=qmax=0.; |
| 975 | nd=0; |
| 976 | lab[0]=lab[1]=lab[2]=-2; |
| 977 | if (flag!=d->GetCoord1()) { n=&np; c=pos; } |
| 978 | } |
| 979 | flag=d->GetCoord1(); |
| 980 | q += d->GetSignal(); |
| 981 | y += d->GetCoord2()*d->GetSignal(); |
| 982 | nd++; |
| 983 | if (d->GetSignal()>qmax) { |
| 984 | qmax=d->GetSignal(); |
| 985 | lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2); |
| 986 | } |
| 987 | curr=strip; |
| 988 | } |
| 989 | c[*n].SetY(y/q); |
| 990 | c[*n].SetQ(q); |
| 991 | c[*n].SetNd(nd); |
| 992 | c[*n].SetLabels(lab); |
| 993 | //Split suspiciously big cluster |
| 994 | if (nd>3) { |
| 995 | c[*n].SetY(y/q-0.5*nd); |
| 996 | c[*n].SetQ(0.5*q); |
| 997 | (*n)++; |
| 998 | if (*n==MAX) { |
| 999 | cerr<<"AliITSclustererV2::FindClustersSSD: Too many 1D clusters !\n"; |
| 1000 | return; |
| 1001 | } |
| 1002 | c[*n].SetY(y/q+0.5*nd); |
| 1003 | c[*n].SetQ(0.5*q); |
| 1004 | c[*n].SetNd(nd); |
| 1005 | c[*n].SetLabels(lab); |
| 1006 | } |
| 1007 | (*n)++; |
| 1008 | if (*n==MAX) { |
| 1009 | cerr<<"AliITSclustererV2::FindClustersSSD: Too many 1D clusters !\n"; |
| 1010 | return; |
| 1011 | } |
| 1012 | |
| c391f9d9 |
1013 | FindClustersSSD(neg, nn, pos, np, clusters); |
| 1014 | } |
| 7f4044f0 |
1015 | |
| c391f9d9 |
1016 | void AliITSclustererV2::FindClustersSSD(AliITSRawStream* input, |
| 1017 | TClonesArray** clusters) |
| 1018 | { |
| 1019 | //------------------------------------------------------------ |
| 1020 | // Actual SSD cluster finder for raw data |
| 1021 | //------------------------------------------------------------ |
| 1022 | Int_t nClustersSSD = 0; |
| 1023 | const Int_t MAX = 1000; |
| 1024 | Ali1Dcluster clusters1D[2][MAX]; |
| 1025 | Int_t nClusters[2] = {0, 0}; |
| 1026 | Float_t q = 0.; |
| 1027 | Float_t y = 0.; |
| 1028 | Int_t nDigits = 0; |
| 1029 | Int_t prevStrip = -1; |
| 1030 | Int_t prevFlag = -1; |
| 1031 | |
| 1032 | // read raw data input stream |
| 1033 | while (kTRUE) { |
| 1034 | Bool_t next = input->Next(); |
| 1035 | |
| 1036 | // check if a new cluster starts |
| 1037 | Int_t strip = input->GetCoord2(); |
| 1038 | Int_t flag = input->GetCoord1(); |
| 1039 | if ((!next || input->IsNewModule() || |
| 1040 | (strip-prevStrip > 1) || (flag != prevFlag)) && |
| 1041 | (nDigits > 0)) { |
| 1042 | if (nClusters[prevFlag] == MAX) { |
| 1043 | Error("FindClustersSSD", "Too many 1D clusters !"); |
| 1044 | return; |
| 1045 | } |
| 1046 | Ali1Dcluster& cluster = clusters1D[prevFlag][nClusters[prevFlag]++]; |
| 1047 | cluster.SetY(y/q); |
| 1048 | cluster.SetQ(q); |
| 1049 | cluster.SetNd(nDigits); |
| 1050 | |
| 1051 | //Split suspiciously big cluster |
| 1052 | if (nDigits > 3) { |
| 1053 | cluster.SetY(y/q - 0.5*nDigits); |
| 1054 | cluster.SetQ(0.5*q); |
| 1055 | if (nClusters[prevFlag] == MAX) { |
| 1056 | Error("FindClustersSSD", "Too many 1D clusters !"); |
| 1057 | return; |
| 1058 | } |
| 1059 | Ali1Dcluster& cluster2 = clusters1D[prevFlag][nClusters[prevFlag]++]; |
| 1060 | cluster2.SetY(y/q + 0.5*nDigits); |
| 1061 | cluster2.SetQ(0.5*q); |
| 1062 | cluster2.SetNd(nDigits); |
| 7f4044f0 |
1063 | } |
| c391f9d9 |
1064 | y = q = 0.; |
| 1065 | nDigits = 0; |
| 7f4044f0 |
1066 | } |
| c391f9d9 |
1067 | |
| 1068 | if (!next || input->IsNewModule()) { |
| 1069 | Int_t iModule = input->GetPrevModuleID(); |
| 1070 | |
| 1071 | // when all data from a module was read, search for clusters |
| 1072 | if (prevFlag >= 0) { |
| 1073 | clusters[iModule] = new TClonesArray("AliITSclusterV2"); |
| 1074 | fI = iModule; |
| 1075 | FindClustersSSD(&clusters1D[0][0], nClusters[0], |
| 1076 | &clusters1D[1][0], nClusters[1], clusters[iModule]); |
| 1077 | Int_t nClusters = clusters[iModule]->GetEntriesFast(); |
| 1078 | nClustersSSD += nClusters; |
| 1079 | } |
| 1080 | |
| 1081 | if (!next) break; |
| 1082 | nClusters[0] = nClusters[1] = 0; |
| 1083 | y = q = 0.; |
| 1084 | nDigits = 0; |
| 7f4044f0 |
1085 | } |
| c391f9d9 |
1086 | |
| 1087 | // add digit to current cluster |
| 1088 | q += input->GetSignal(); |
| 1089 | y += strip * input->GetSignal(); |
| 1090 | nDigits++; |
| 1091 | prevStrip = strip; |
| 1092 | prevFlag = flag; |
| 7f4044f0 |
1093 | } |
| c391f9d9 |
1094 | |
| 1095 | Info("FindClustersSSD", "found clusters in ITS SSD: %d", nClustersSSD); |
| 7f4044f0 |
1096 | } |
| 1097 | |
| 1098 | #else //V1 |
| 1099 | |
| 1100 | #include "AliITSDetType.h" |
| 1101 | |
| 1102 | #include "AliITSsegmentationSPD.h" |
| 1103 | #include "AliITSClusterFinderSPD.h" |
| 1104 | |
| 1105 | #include "AliITSresponseSDD.h" |
| 1106 | #include "AliITSsegmentationSDD.h" |
| 1107 | #include "AliITSClusterFinderSDD.h" |
| 1108 | |
| 1109 | #include "AliITSsegmentationSSD.h" |
| 1110 | #include "AliITSClusterFinderSSD.h" |
| 1111 | |
| 1112 | |
| 1113 | void AliITSclustererV2:: |
| 1114 | FindClustersSPD(const TClonesArray *digits, TClonesArray *clusters) { |
| 1115 | //------------------------------------------------------------ |
| 1116 | // Actual SPD cluster finding based on AliITSClusterFinderSPD |
| 1117 | //------------------------------------------------------------ |
| 1118 | static AliITS *its=(AliITS*)gAlice->GetModule("ITS"); |
| 1119 | static TClonesArray *points=its->RecPoints(); |
| 1120 | static AliITSsegmentationSPD *seg= |
| 1121 | (AliITSsegmentationSPD *)its->DetType(0)->GetSegmentationModel(); |
| 1122 | static AliITSClusterFinderSPD cf(seg, (TClonesArray*)digits, points); |
| 1123 | |
| 1124 | cf.FindRawClusters(fI); |
| 1cca57bf |
1125 | RecPoints2Clusters(points, fI, clusters); |
| 7f4044f0 |
1126 | its->ResetRecPoints(); |
| 1127 | |
| 1128 | } |
| 1129 | |
| 1130 | void AliITSclustererV2:: |
| 1131 | FindClustersSDD(const TClonesArray *digits, TClonesArray *clusters) { |
| 1132 | //------------------------------------------------------------ |
| 1133 | // Actual SDD cluster finding based on AliITSClusterFinderSDD |
| 1134 | //------------------------------------------------------------ |
| 1135 | static AliITS *its=(AliITS*)gAlice->GetModule("ITS"); |
| 1136 | static TClonesArray *points=its->RecPoints(); |
| 1137 | static AliITSresponseSDD *resp= |
| 1138 | (AliITSresponseSDD *)its->DetType(1)->GetResponseModel(); |
| 1139 | static AliITSsegmentationSDD *seg= |
| 1140 | (AliITSsegmentationSDD *)its->DetType(1)->GetSegmentationModel(); |
| 1141 | static AliITSClusterFinderSDD |
| 1142 | cf(seg,resp,(TClonesArray*)digits,its->ClustersAddress(1)); |
| 1143 | |
| 1144 | cf.FindRawClusters(fI); |
| 1cca57bf |
1145 | Int_t nc=points->GetEntriesFast(); |
| 1146 | while (nc--) { //To be consistent with the SSD cluster charges |
| 1147 | AliITSRecPoint *p=(AliITSRecPoint*)points->UncheckedAt(nc); |
| 1148 | p->SetQ(p->GetQ/12.); |
| 1149 | } |
| 1150 | RecPoints2Clusters(points, fI, clusters); |
| 7f4044f0 |
1151 | its->ResetClusters(1); |
| 1152 | its->ResetRecPoints(); |
| 1153 | |
| 1154 | } |
| 1155 | |
| 1156 | void AliITSclustererV2:: |
| 1157 | FindClustersSSD(const TClonesArray *digits, TClonesArray *clusters) { |
| 1158 | //------------------------------------------------------------ |
| 1159 | // Actual SSD cluster finding based on AliITSClusterFinderSSD |
| 1160 | //------------------------------------------------------------ |
| 1161 | static AliITS *its=(AliITS*)gAlice->GetModule("ITS"); |
| 1162 | static TClonesArray *points=its->RecPoints(); |
| 1163 | static AliITSsegmentationSSD *seg= |
| 1164 | (AliITSsegmentationSSD *)its->DetType(2)->GetSegmentationModel(); |
| 1165 | static AliITSClusterFinderSSD cf(seg,(TClonesArray*)digits); |
| 1166 | |
| 1167 | cf.FindRawClusters(fI); |
| 1cca57bf |
1168 | RecPoints2Clusters(points, fI, clusters); |
| 7f4044f0 |
1169 | its->ResetRecPoints(); |
| 1170 | |
| 1171 | } |
| 1172 | |
| 1173 | #endif |
git://git.uio.no / u / mrichter / AliRoot.git / blame