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