| | 1 | //------------------------------------------------------------------------- |
| | 2 | // Implementation of the ITS clusterer V2 class |
| | 3 | // |
| | 4 | // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch |
| | 5 | //------------------------------------------------------------------------- |
| | 6 | |
| | 7 | |
| | 8 | #include "AliRun.h" |
| | 9 | |
| | 10 | #include "AliITSclustererV2.h" |
| | 11 | #include "AliITSclusterV2.h" |
| | 12 | #include "AliITSDetTypeRec.h" |
| | 13 | #include "AliRawReader.h" |
| | 14 | #include "AliITSRawStreamSPD.h" |
| | 15 | #include "AliITSRawStreamSDD.h" |
| | 16 | #include "AliITSRawStreamSSD.h" |
| | 17 | |
| | 18 | #include <TFile.h> |
| | 19 | #include <TTree.h> |
| | 20 | #include <TClonesArray.h> |
| | 21 | #include "AliITSgeomTGeo.h" |
| | 22 | #include "AliITSdigitSPD.h" |
| | 23 | #include "AliITSdigitSDD.h" |
| | 24 | #include "AliITSdigitSSD.h" |
| | 25 | #include "AliMC.h" |
| | 26 | |
| | 27 | ClassImp(AliITSclustererV2) |
| | 28 | |
| | 29 | extern AliRun *gAlice; |
| | 30 | |
| | 31 | AliITSclustererV2::AliITSclustererV2(): |
| | 32 | fNModules(0), |
| | 33 | fEvent(0), |
| | 34 | fI(0), |
| | 35 | fLastSPD1(0), |
| | 36 | fNySPD(0), |
| | 37 | fNzSPD(0), |
| | 38 | fYpitchSPD(0), |
| | 39 | fZ1pitchSPD(0), |
| | 40 | fZ2pitchSPD(0), |
| | 41 | fHwSPD(0), |
| | 42 | fHlSPD(0), |
| | 43 | fNySDD(0), |
| | 44 | fNzSDD(0), |
| | 45 | fYpitchSDD(0), |
| | 46 | fZpitchSDD(0), |
| | 47 | fHwSDD(0), |
| | 48 | fHlSDD(0), |
| | 49 | fYoffSDD(0), |
| | 50 | fLastSSD1(0), |
| | 51 | fYpitchSSD(0), |
| | 52 | fHwSSD(0), |
| | 53 | fHlSSD(0), |
| | 54 | fTanP(0), |
| | 55 | fTanN(0){ |
| | 56 | //default constructor |
| | 57 | } |
| | 58 | AliITSclustererV2::AliITSclustererV2(const Char_t *geom): |
| | 59 | fNModules(AliITSgeomTGeo::GetNModules()), |
| | 60 | fEvent(0), |
| | 61 | fI(0), |
| | 62 | fLastSPD1(AliITSgeomTGeo::GetModuleIndex(2,1,1)-1), |
| | 63 | fNySPD(256), |
| | 64 | fNzSPD(160), |
| | 65 | fYpitchSPD(0.0050), |
| | 66 | fZ1pitchSPD(0.0425), |
| | 67 | fZ2pitchSPD(0.0625), |
| | 68 | fHwSPD(0.64), |
| | 69 | fHlSPD(3.48), |
| | 70 | fNySDD(256), |
| | 71 | fNzSDD(256), |
| | 72 | fYpitchSDD(0.01825), |
| | 73 | fZpitchSDD(0.02940), |
| | 74 | fHwSDD(3.5085), |
| | 75 | fHlSDD(3.7632), |
| | 76 | fYoffSDD(0.0425), |
| | 77 | fLastSSD1(AliITSgeomTGeo::GetModuleIndex(6,1,1)-1), |
| | 78 | fYpitchSSD(0.0095), |
| | 79 | fHwSSD(3.65), |
| | 80 | fHlSSD(2.00), |
| | 81 | fTanP(0.0275), |
| | 82 | fTanN(0.0075) { |
| | 83 | //------------------------------------------------------------ |
| | 84 | // Standard constructor |
| | 85 | //------------------------------------------------------------ |
| | 86 | if (geom) { |
| | 87 | AliWarning("\"geom\" is actually a dummy argument !"); |
| | 88 | } |
| | 89 | Int_t m; |
| | 90 | for (m=0; m<fNModules; m++) { |
| | 91 | Int_t lay,lad,det; AliITSgeomTGeo::GetModuleId(m,lay,lad,det); |
| | 92 | const TGeoHMatrix *tm=AliITSgeomTGeo::GetTracking2LocalMatrix(m); |
| | 93 | fYshift[m] = (tm->Inverse()).GetTranslation()[1]; |
| | 94 | fZshift[m] = (tm->Inverse()).GetTranslation()[2]; |
| | 95 | fNdet[m] = (lad-1)*AliITSgeomTGeo::GetNDetectors(lay) + (det-1); |
| | 96 | fNlayer[m] = lay-1; |
| | 97 | } |
| | 98 | |
| | 99 | //SPD geometry |
| | 100 | fYSPD[0]=0.5*fYpitchSPD; |
| | 101 | for (m=1; m<fNySPD; m++) fYSPD[m]=fYSPD[m-1]+fYpitchSPD; |
| | 102 | fZSPD[0]=fZ1pitchSPD; |
| | 103 | for (m=1; m<fNzSPD; m++) { |
| | 104 | Double_t dz=fZ1pitchSPD; |
| | 105 | if (m==31 || m==32 || m==63 || m==64 || m==95 || m==96 || |
| | 106 | m==127 || m==128) dz=fZ2pitchSPD; |
| | 107 | fZSPD[m]=fZSPD[m-1]+dz; |
| | 108 | } |
| | 109 | for (m=0; m<fNzSPD; m++) { |
| | 110 | Double_t dz=0.5*fZ1pitchSPD; |
| | 111 | if (m==31 || m==32 || m==63 || m==64 || m==95 || m==96 || |
| | 112 | m==127 || m==128) dz=0.5*fZ2pitchSPD; |
| | 113 | fZSPD[m]-=dz; |
| | 114 | } |
| | 115 | |
| | 116 | } |
| | 117 | |
| | 118 | |
| | 119 | Int_t AliITSclustererV2::Digits2Clusters(TTree *dTree, TTree *cTree) { |
| | 120 | //------------------------------------------------------------ |
| | 121 | // This function creates ITS clusters |
| | 122 | //------------------------------------------------------------ |
| | 123 | Int_t ncl=0; |
| | 124 | |
| | 125 | if (!dTree) { |
| | 126 | Error("Digits2Clusters","Can't get the tree with digits !"); |
| | 127 | return 1; |
| | 128 | } |
| | 129 | |
| | 130 | TClonesArray *digitsSPD=new TClonesArray("AliITSdigitSPD",3000); |
| | 131 | dTree->SetBranchAddress("ITSDigitsSPD",&digitsSPD); |
| | 132 | TClonesArray *digitsSDD=new TClonesArray("AliITSdigitSDD",3000); |
| | 133 | dTree->SetBranchAddress("ITSDigitsSDD",&digitsSDD); |
| | 134 | TClonesArray *digitsSSD=new TClonesArray("AliITSdigitSSD",3000); |
| | 135 | dTree->SetBranchAddress("ITSDigitsSSD",&digitsSSD); |
| | 136 | |
| | 137 | TClonesArray *clusters=new TClonesArray("AliITSclusterV2",1000); |
| | 138 | TBranch *branch=cTree->GetBranch("Clusters"); |
| | 139 | if (!branch) cTree->Branch("Clusters",&clusters); |
| | 140 | else branch->SetAddress(&clusters); |
| | 141 | |
| | 142 | Int_t mmax=(Int_t)dTree->GetEntries(); |
| | 143 | if (mmax!=fNModules) { |
| | 144 | Error("Digits2Clusters","Number of entries != number of modules !"); |
| | 145 | return 1; |
| | 146 | } |
| | 147 | |
| | 148 | for (fI=0; fI<mmax; fI++) { |
| | 149 | dTree->GetEvent(fI); |
| | 150 | |
| | 151 | if (digitsSPD->GetEntriesFast()!=0) |
| | 152 | FindClustersSPD(digitsSPD,clusters); |
| | 153 | else |
| | 154 | if(digitsSDD->GetEntriesFast()!=0) |
| | 155 | FindClustersSDD(digitsSDD,clusters); |
| | 156 | else if(digitsSSD->GetEntriesFast()!=0) |
| | 157 | FindClustersSSD(digitsSSD,clusters); |
| | 158 | |
| | 159 | ncl+=clusters->GetEntriesFast(); |
| | 160 | |
| | 161 | cTree->Fill(); |
| | 162 | |
| | 163 | digitsSPD->Clear(); |
| | 164 | digitsSDD->Clear(); |
| | 165 | digitsSSD->Clear(); |
| | 166 | clusters->Clear(); |
| | 167 | } |
| | 168 | |
| | 169 | //cTree->Write(); |
| | 170 | |
| | 171 | delete clusters; |
| | 172 | |
| | 173 | delete digitsSPD; |
| | 174 | delete digitsSDD; |
| | 175 | delete digitsSSD; |
| | 176 | |
| | 177 | //delete dTree; |
| | 178 | |
| | 179 | Info("Digits2Clusters","Number of found clusters : %d",ncl); |
| | 180 | |
| | 181 | return 0; |
| | 182 | } |
| | 183 | |
| | 184 | void AliITSclustererV2::Digits2Clusters(AliRawReader* rawReader) { |
| | 185 | //------------------------------------------------------------ |
| | 186 | // This function creates ITS clusters from raw data |
| | 187 | //------------------------------------------------------------ |
| | 188 | AliRunLoader* runLoader = AliRunLoader::GetRunLoader(); |
| | 189 | if (!runLoader) { |
| | 190 | Error("Digits2Clusters", "no run loader found"); |
| | 191 | return; |
| | 192 | } |
| | 193 | runLoader->LoadKinematics(); |
| | 194 | AliLoader* itsLoader = runLoader->GetLoader("ITSLoader"); |
| | 195 | if (!itsLoader) { |
| | 196 | Error("Digits2Clusters", "no loader for ITS found"); |
| | 197 | return; |
| | 198 | } |
| | 199 | if (!itsLoader->TreeR()) itsLoader->MakeTree("R"); |
| | 200 | TTree* cTree = itsLoader->TreeR(); |
| | 201 | |
| | 202 | TClonesArray *array=new TClonesArray("AliITSclusterV2",1000); |
| | 203 | cTree->Branch("Clusters",&array); |
| | 204 | delete array; |
| | 205 | |
| | 206 | TClonesArray** clusters = new TClonesArray*[fNModules]; |
| | 207 | for (Int_t iModule = 0; iModule < fNModules; iModule++) { |
| | 208 | clusters[iModule] = NULL; |
| | 209 | } |
| | 210 | // one TClonesArray per module |
| | 211 | |
| | 212 | rawReader->Reset(); |
| | 213 | AliITSRawStreamSPD inputSPD(rawReader); |
| | 214 | FindClustersSPD(&inputSPD, clusters); |
| | 215 | |
| | 216 | rawReader->Reset(); |
| | 217 | AliITSRawStreamSDD inputSDD(rawReader); |
| | 218 | FindClustersSDD(&inputSDD, clusters); |
| | 219 | |
| | 220 | rawReader->Reset(); |
| | 221 | AliITSRawStreamSSD inputSSD(rawReader); |
| | 222 | FindClustersSSD(&inputSSD, clusters); |
| | 223 | |
| | 224 | // write all clusters to the tree |
| | 225 | Int_t nClusters = 0; |
| | 226 | TClonesArray *emptyArray=new TClonesArray("AliITSclusterV2"); |
| | 227 | for (Int_t iModule = 0; iModule < fNModules; iModule++) { |
| | 228 | array = clusters[iModule]; |
| | 229 | if (!array) { |
| | 230 | Error("Digits2Clusters", "data for module %d missing!", iModule); |
| | 231 | array = emptyArray; |
| | 232 | } |
| | 233 | cTree->SetBranchAddress("Clusters", &array); |
| | 234 | cTree->Fill(); |
| | 235 | nClusters += array->GetEntriesFast(); |
| | 236 | } |
| | 237 | delete emptyArray; |
| | 238 | |
| | 239 | itsLoader->WriteRecPoints("OVERWRITE"); |
| | 240 | |
| | 241 | delete[] clusters; |
| | 242 | |
| | 243 | Info("Digits2Clusters", "total number of found clusters in ITS: %d\n", |
| | 244 | nClusters); |
| | 245 | } |
| | 246 | |
| | 247 | //**** Fast clusters ******************************* |
| | 248 | #include "TParticle.h" |
| | 249 | |
| | 250 | //#include "AliITS.h" |
| | 251 | #include "AliITSmodule.h" |
| | 252 | #include "AliITSRecPoint.h" |
| | 253 | #include "AliITSsimulationFastPoints.h" |
| | 254 | #include "AliITSRecPoint.h" |
| | 255 | |
| | 256 | |
| | 257 | static void CheckLabels(Int_t lab[3]) { |
| | 258 | //------------------------------------------------------------ |
| | 259 | // Tries to find mother's labels |
| | 260 | //------------------------------------------------------------ |
| | 261 | AliRunLoader *rl = AliRunLoader::GetRunLoader(); |
| | 262 | TTree *trK=(TTree*)rl->TreeK(); |
| | 263 | |
| | 264 | if(trK){ |
| | 265 | Int_t label=lab[0]; |
| | 266 | if (label>=0) { |
| | 267 | TParticle *part=(TParticle*)gAlice->GetMCApp()->Particle(label); |
| | 268 | label=-3; |
| | 269 | while (part->P() < 0.005) { |
| | 270 | Int_t m=part->GetFirstMother(); |
| | 271 | if (m<0) { |
| | 272 | Info("CheckLabels","Primary momentum: %f",part->P()); |
| | 273 | break; |
| | 274 | } |
| | 275 | if (part->GetStatusCode()>0) { |
| | 276 | Info("CheckLabels","Primary momentum: %f",part->P()); |
| | 277 | break; |
| | 278 | } |
| | 279 | label=m; |
| | 280 | part=(TParticle*)gAlice->GetMCApp()->Particle(label); |
| | 281 | } |
| | 282 | if (lab[1]<0) lab[1]=label; |
| | 283 | else if (lab[2]<0) lab[2]=label; |
| | 284 | else ;//cerr<<"CheckLabels : No empty labels !\n"; |
| | 285 | } |
| | 286 | } |
| | 287 | } |
| | 288 | |
| | 289 | /* |
| | 290 | static void CheckLabels(Int_t lab[3]) { |
| | 291 | //------------------------------------------------------------ |
| | 292 | // Tries to find mother's labels |
| | 293 | //------------------------------------------------------------ |
| | 294 | |
| | 295 | if(lab[0]<0 && lab[1]<0 && lab[2]<0) return; // In case of no labels just exit |
| | 296 | |
| | 297 | Int_t ntracks = gAlice->GetMCApp()->GetNtrack(); |
| | 298 | for (Int_t i=0;i<3;i++){ |
| | 299 | Int_t label = lab[i]; |
| | 300 | if (label>=0 && label<ntracks) { |
| | 301 | TParticle *part=(TParticle*)gAlice->GetMCApp()->Particle(label); |
| | 302 | if (part->P() < 0.005) { |
| | 303 | Int_t m=part->GetFirstMother(); |
| | 304 | if (m<0) { |
| | 305 | continue; |
| | 306 | } |
| | 307 | if (part->GetStatusCode()>0) { |
| | 308 | continue; |
| | 309 | } |
| | 310 | lab[i]=m; |
| | 311 | } |
| | 312 | } |
| | 313 | } |
| | 314 | |
| | 315 | } |
| | 316 | */ |
| | 317 | static void CheckLabels2(Int_t lab[10]) { |
| | 318 | //------------------------------------------------------------ |
| | 319 | // Tries to find mother's labels |
| | 320 | //------------------------------------------------------------ |
| | 321 | AliRunLoader *rl = AliRunLoader::GetRunLoader(); |
| | 322 | TTree *trK=(TTree*)rl->TreeK(); |
| | 323 | |
| | 324 | if(trK){ |
| | 325 | Int_t nlabels =0; |
| | 326 | for (Int_t i=0;i<10;i++) if (lab[i]>=0) nlabels++; |
| | 327 | if(nlabels == 0) return; // In case of no labels just exit |
| | 328 | |
| | 329 | Int_t ntracks = gAlice->GetMCApp()->GetNtrack(); |
| | 330 | |
| | 331 | for (Int_t i=0;i<10;i++){ |
| | 332 | Int_t label = lab[i]; |
| | 333 | if (label>=0 && label<ntracks) { |
| | 334 | TParticle *part=(TParticle*)gAlice->GetMCApp()->Particle(label); |
| | 335 | if (part->P() < 0.02) { |
| | 336 | Int_t m=part->GetFirstMother(); |
| | 337 | if (m<0) { |
| | 338 | continue; |
| | 339 | } |
| | 340 | if (part->GetStatusCode()>0) { |
| | 341 | continue; |
| | 342 | } |
| | 343 | lab[i]=m; |
| | 344 | } |
| | 345 | else |
| | 346 | if (part->P() < 0.12 && nlabels>3) { |
| | 347 | lab[i]=-2; |
| | 348 | nlabels--; |
| | 349 | } |
| | 350 | } |
| | 351 | else{ |
| | 352 | if ( (label>ntracks||label <0) && nlabels>3) { |
| | 353 | lab[i]=-2; |
| | 354 | nlabels--; |
| | 355 | } |
| | 356 | } |
| | 357 | } |
| | 358 | if (nlabels>3){ |
| | 359 | for (Int_t i=0;i<10;i++){ |
| | 360 | if (nlabels>3){ |
| | 361 | Int_t label = lab[i]; |
| | 362 | if (label>=0 && label<ntracks) { |
| | 363 | TParticle *part=(TParticle*)gAlice->GetMCApp()->Particle(label); |
| | 364 | if (part->P() < 0.1) { |
| | 365 | lab[i]=-2; |
| | 366 | nlabels--; |
| | 367 | } |
| | 368 | } |
| | 369 | } |
| | 370 | } |
| | 371 | } |
| | 372 | |
| | 373 | //compress labels -- if multi-times the same |
| | 374 | Int_t lab2[10]; |
| | 375 | for (Int_t i=0;i<10;i++) lab2[i]=-2; |
| | 376 | for (Int_t i=0;i<10 ;i++){ |
| | 377 | if (lab[i]<0) continue; |
| | 378 | for (Int_t j=0;j<10 &&lab2[j]!=lab[i];j++){ |
| | 379 | if (lab2[j]<0) { |
| | 380 | lab2[j]= lab[i]; |
| | 381 | break; |
| | 382 | } |
| | 383 | } |
| | 384 | } |
| | 385 | for (Int_t j=0;j<10;j++) lab[j]=lab2[j]; |
| | 386 | } |
| | 387 | } |
| | 388 | |
| | 389 | static void AddLabel(Int_t lab[10], Int_t label) { |
| | 390 | // add label of the particle in the kine tree which originates this cluster |
| | 391 | // (used for reconstruction of MC data only, for comparison purposes) |
| | 392 | AliRunLoader *rl = AliRunLoader::GetRunLoader(); |
| | 393 | TTree *trK=(TTree*)rl->TreeK(); |
| | 394 | |
| | 395 | if(trK){ |
| | 396 | if(label<0) return; // In case of no label just exit |
| | 397 | |
| | 398 | Int_t ntracks = gAlice->GetMCApp()->GetNtrack(); |
| | 399 | if (label>ntracks) return; |
| | 400 | for (Int_t i=0;i<10;i++){ |
| | 401 | // if (label<0) break; |
| | 402 | if (lab[i]==label) break; |
| | 403 | if (lab[i]<0) { |
| | 404 | lab[i]= label; |
| | 405 | break; |
| | 406 | } |
| | 407 | } |
| | 408 | } |
| | 409 | } |
| | 410 | |
| | 411 | void AliITSclustererV2::RecPoints2Clusters |
| | 412 | (const TClonesArray *points, Int_t idx, TClonesArray *clusters) { |
| | 413 | //------------------------------------------------------------ |
| | 414 | // Conversion AliITSRecPoint -> AliITSclusterV2 for the ITS |
| | 415 | // subdetector indexed by idx |
| | 416 | //------------------------------------------------------------ |
| | 417 | TClonesArray &cl=*clusters; |
| | 418 | Int_t ncl=points->GetEntriesFast(); |
| | 419 | for (Int_t i=0; i<ncl; i++) { |
| | 420 | AliITSRecPoint *p = (AliITSRecPoint *)points->UncheckedAt(i); |
| | 421 | Float_t lp[5]; |
| | 422 | lp[0]=-(-p->GetDetLocalX()+fYshift[idx]); if (idx<=fLastSPD1) lp[0]*=-1; //SPD1 |
| | 423 | lp[1]= -p->GetZ()+fZshift[idx]; |
| | 424 | lp[2]=p->GetSigmaDetLocX2(); |
| | 425 | lp[3]=p->GetSigmaZ2(); |
| | 426 | lp[4]=p->GetQ()*36./23333.; //electrons -> ADC |
| | 427 | Int_t lab[4]; |
| | 428 | lab[0]=p->GetLabel(0); lab[1]=p->GetLabel(1); lab[2]=p->GetLabel(2); |
| | 429 | lab[3]=fNdet[idx]; |
| | 430 | CheckLabels(lab); |
| | 431 | Int_t dummy[3]={0,0,0}; |
| | 432 | new (cl[i]) AliITSclusterV2(lab,lp, dummy); |
| | 433 | } |
| | 434 | } |
| | 435 | |
| | 436 | //*********************************** |
| | 437 | |
| | 438 | |
| | 439 | void AliITSclustererV2:: |
| | 440 | FindCluster(Int_t k,Int_t maxz,AliBin *bins,Int_t &n,Int_t *idx) { |
| | 441 | //------------------------------------------------------------ |
| | 442 | // returns an array of indices of digits belonging to the cluster |
| | 443 | // (needed when the segmentation is not regular) |
| | 444 | //------------------------------------------------------------ |
| | 445 | if (n<200) idx[n++]=bins[k].GetIndex(); |
| | 446 | bins[k].Use(); |
| | 447 | |
| | 448 | if (bins[k-maxz].IsNotUsed()) FindCluster(k-maxz,maxz,bins,n,idx); |
| | 449 | if (bins[k-1 ].IsNotUsed()) FindCluster(k-1 ,maxz,bins,n,idx); |
| | 450 | if (bins[k+maxz].IsNotUsed()) FindCluster(k+maxz,maxz,bins,n,idx); |
| | 451 | if (bins[k+1 ].IsNotUsed()) FindCluster(k+1 ,maxz,bins,n,idx); |
| | 452 | /* |
| | 453 | if (bins[k-maxz-1].IsNotUsed()) FindCluster(k-maxz-1,maxz,bins,n,idx); |
| | 454 | if (bins[k-maxz+1].IsNotUsed()) FindCluster(k-maxz+1,maxz,bins,n,idx); |
| | 455 | if (bins[k+maxz-1].IsNotUsed()) FindCluster(k+maxz-1,maxz,bins,n,idx); |
| | 456 | if (bins[k+maxz+1].IsNotUsed()) FindCluster(k+maxz+1,maxz,bins,n,idx); |
| | 457 | */ |
| | 458 | } |
| | 459 | |
| | 460 | void AliITSclustererV2:: |
| | 461 | FindClustersSPD(const TClonesArray *digits, TClonesArray *clusters) { |
| | 462 | //------------------------------------------------------------ |
| | 463 | // Actual SPD cluster finder |
| | 464 | //------------------------------------------------------------ |
| | 465 | Int_t kNzBins = fNzSPD + 2; |
| | 466 | const Int_t kMAXBIN=kNzBins*(fNySPD+2); |
| | 467 | |
| | 468 | Int_t ndigits=digits->GetEntriesFast(); |
| | 469 | AliBin *bins=new AliBin[kMAXBIN]; |
| | 470 | |
| | 471 | Int_t k; |
| | 472 | AliITSdigitSPD *d=0; |
| | 473 | for (k=0; k<ndigits; k++) { |
| | 474 | d=(AliITSdigitSPD*)digits->UncheckedAt(k); |
| | 475 | Int_t i=d->GetCoord2()+1; //y |
| | 476 | Int_t j=d->GetCoord1()+1; |
| | 477 | bins[i*kNzBins+j].SetIndex(k); |
| | 478 | bins[i*kNzBins+j].SetMask(1); |
| | 479 | } |
| | 480 | |
| | 481 | Int_t n=0; TClonesArray &cl=*clusters; |
| | 482 | for (k=0; k<kMAXBIN; k++) { |
| | 483 | if (!bins[k].IsNotUsed()) continue; |
| | 484 | Int_t ni=0, idx[200]; |
| | 485 | FindCluster(k,kNzBins,bins,ni,idx); |
| | 486 | if (ni==200) { |
| | 487 | Info("FindClustersSPD","Too big cluster !"); |
| | 488 | continue; |
| | 489 | } |
| | 490 | Int_t milab[10]; |
| | 491 | for (Int_t ilab=0;ilab<10;ilab++){ |
| | 492 | milab[ilab]=-2; |
| | 493 | } |
| | 494 | |
| | 495 | d=(AliITSdigitSPD*)digits->UncheckedAt(idx[0]); |
| | 496 | Int_t ymin=d->GetCoord2(),ymax=ymin; |
| | 497 | Int_t zmin=d->GetCoord1(),zmax=zmin; |
| | 498 | |
| | 499 | for (Int_t l=0; l<ni; l++) { |
| | 500 | d=(AliITSdigitSPD*)digits->UncheckedAt(idx[l]); |
| | 501 | |
| | 502 | if (ymin > d->GetCoord2()) ymin=d->GetCoord2(); |
| | 503 | if (ymax < d->GetCoord2()) ymax=d->GetCoord2(); |
| | 504 | if (zmin > d->GetCoord1()) zmin=d->GetCoord1(); |
| | 505 | if (zmax < d->GetCoord1()) zmax=d->GetCoord1(); |
| | 506 | // MI addition - find all labels in cluster |
| | 507 | for (Int_t dlab=0;dlab<10;dlab++){ |
| | 508 | Int_t digitlab = (d->GetTracks())[dlab]; |
| | 509 | if (digitlab<0) continue; |
| | 510 | AddLabel(milab,digitlab); |
| | 511 | } |
| | 512 | if (milab[9]>0) CheckLabels2(milab); |
| | 513 | } |
| | 514 | CheckLabels2(milab); |
| | 515 | // |
| | 516 | //Int_t idy = (fNlayer[fI]==0)? 2:3; |
| | 517 | //for (Int_t iz=zmin; iz<=zmax;iz+=2) |
| | 518 | //Int_t idy = (ymax-ymin)/4.; // max 2 clusters |
| | 519 | Int_t idy = 0; // max 2 clusters |
| | 520 | if (fNlayer[fI]==0 &&idy<3) idy=3; |
| | 521 | if (fNlayer[fI]==1 &&idy<4) idy=4; |
| | 522 | Int_t idz =3; |
| | 523 | for (Int_t iz=zmin; iz<=zmax;iz+=idz) |
| | 524 | for (Int_t iy=ymin; iy<=ymax;iy+=idy){ |
| | 525 | // |
| | 526 | Int_t nodigits =0; |
| | 527 | Float_t y=0.,z=0.,q=0.; |
| | 528 | for (Int_t l=0; l<ni; l++) { |
| | 529 | d=(AliITSdigitSPD*)digits->UncheckedAt(idx[l]); |
| | 530 | if (zmax-zmin>=idz || ymax-ymin>=idy){ |
| | 531 | if (TMath::Abs( d->GetCoord2()-iy)>0.75*idy) continue; |
| | 532 | if (TMath::Abs( d->GetCoord1()-iz)>0.75*idz) continue; |
| | 533 | } |
| | 534 | nodigits++; |
| | 535 | Float_t qq=d->GetSignal(); |
| | 536 | y+=qq*fYSPD[d->GetCoord2()]; z+=qq*fZSPD[d->GetCoord1()]; q+=qq; |
| | 537 | |
| | 538 | } |
| | 539 | if (nodigits==0) continue; |
| | 540 | y/=q; z/=q; |
| | 541 | y-=fHwSPD; z-=fHlSPD; |
| | 542 | |
| | 543 | Float_t lp[5]; |
| | 544 | lp[0]=-(-y+fYshift[fI]); if (fI<=fLastSPD1) lp[0]=-lp[0]; |
| | 545 | lp[1]= -z+fZshift[fI]; |
| | 546 | // Float_t factor=TMath::Max(double(ni-3.),1.5); |
| | 547 | Float_t factory=TMath::Max(ymax-ymin,1); |
| | 548 | Float_t factorz=TMath::Max(zmax-zmin,1); |
| | 549 | factory*= factory; |
| | 550 | factorz*= factorz; |
| | 551 | //lp[2]= (fYpitchSPD*fYpitchSPD/12.)*factory; |
| | 552 | //lp[3]= (fZ1pitchSPD*fZ1pitchSPD/12.)*factorz; |
| | 553 | lp[2]= (fYpitchSPD*fYpitchSPD/12.); |
| | 554 | lp[3]= (fZ1pitchSPD*fZ1pitchSPD/12.); |
| | 555 | //lp[4]= q; |
| | 556 | lp[4]= (zmax-zmin+1)*100 + (ymax-ymin+1); |
| | 557 | |
| | 558 | milab[3]=fNdet[fI]; |
| | 559 | d=(AliITSdigitSPD*)digits->UncheckedAt(idx[0]); |
| | 560 | Int_t info[3] = {ymax-ymin+1,zmax-zmin+1,fNlayer[fI]}; |
| | 561 | new (cl[n]) AliITSclusterV2(milab,lp,info); n++; |
| | 562 | } |
| | 563 | } |
| | 564 | |
| | 565 | delete [] bins; |
| | 566 | } |
| | 567 | |
| | 568 | void AliITSclustererV2::FindClustersSPD(AliITSRawStream* input, |
| | 569 | TClonesArray** clusters) |
| | 570 | { |
| | 571 | //------------------------------------------------------------ |
| | 572 | // Actual SPD cluster finder for raw data |
| | 573 | //------------------------------------------------------------ |
| | 574 | |
| | 575 | Int_t nClustersSPD = 0; |
| | 576 | Int_t kNzBins = fNzSPD + 2; |
| | 577 | Int_t kNyBins = fNySPD + 2; |
| | 578 | Int_t kMaxBin = kNzBins * kNyBins; |
| | 579 | AliBin *binsSPD = new AliBin[kMaxBin]; |
| | 580 | AliBin *binsSPDInit = new AliBin[kMaxBin]; |
| | 581 | AliBin *bins = NULL; |
| | 582 | |
| | 583 | // read raw data input stream |
| | 584 | while (kTRUE) { |
| | 585 | Bool_t next = input->Next(); |
| | 586 | if (!next || input->IsNewModule()) { |
| | 587 | Int_t iModule = input->GetPrevModuleID(); |
| | 588 | |
| | 589 | // when all data from a module was read, search for clusters |
| | 590 | if (bins) { |
| | 591 | clusters[iModule] = new TClonesArray("AliITSclusterV2"); |
| | 592 | Int_t nClusters = 0; |
| | 593 | |
| | 594 | for (Int_t iBin = 0; iBin < kMaxBin; iBin++) { |
| | 595 | if (bins[iBin].IsUsed()) continue; |
| | 596 | Int_t nBins = 0; |
| | 597 | Int_t idxBins[200]; |
| | 598 | FindCluster(iBin, kNzBins, bins, nBins, idxBins); |
| | 599 | if (nBins == 200) { |
| | 600 | Error("FindClustersSPD", "SPD: Too big cluster !\n"); |
| | 601 | continue; |
| | 602 | } |
| | 603 | |
| | 604 | Int_t label[4]; |
| | 605 | label[0] = -2; |
| | 606 | label[1] = -2; |
| | 607 | label[2] = -2; |
| | 608 | // label[3] = iModule; |
| | 609 | label[3] = fNdet[iModule]; |
| | 610 | |
| | 611 | Int_t ymin = (idxBins[0] / kNzBins) - 1; |
| | 612 | Int_t ymax = ymin; |
| | 613 | Int_t zmin = (idxBins[0] % kNzBins) - 1; |
| | 614 | Int_t zmax = zmin; |
| | 615 | for (Int_t idx = 0; idx < nBins; idx++) { |
| | 616 | Int_t iy = (idxBins[idx] / kNzBins) - 1; |
| | 617 | Int_t iz = (idxBins[idx] % kNzBins) - 1; |
| | 618 | if (ymin > iy) ymin = iy; |
| | 619 | if (ymax < iy) ymax = iy; |
| | 620 | if (zmin > iz) zmin = iz; |
| | 621 | if (zmax < iz) zmax = iz; |
| | 622 | } |
| | 623 | |
| | 624 | Int_t idy = 0; // max 2 clusters |
| | 625 | if ((iModule <= fLastSPD1) &&idy<3) idy=3; |
| | 626 | if ((iModule > fLastSPD1) &&idy<4) idy=4; |
| | 627 | Int_t idz =3; |
| | 628 | for (Int_t iiz=zmin; iiz<=zmax;iiz+=idz) |
| | 629 | for (Int_t iiy=ymin; iiy<=ymax;iiy+=idy){ |
| | 630 | // |
| | 631 | Int_t ndigits =0; |
| | 632 | Float_t y=0.,z=0.,q=0.; |
| | 633 | for (Int_t idx = 0; idx < nBins; idx++) { |
| | 634 | Int_t iy = (idxBins[idx] / kNzBins) - 1; |
| | 635 | Int_t iz = (idxBins[idx] % kNzBins) - 1; |
| | 636 | if (zmax-zmin>=idz || ymax-ymin>=idy){ |
| | 637 | if (TMath::Abs(iy-iiy)>0.75*idy) continue; |
| | 638 | if (TMath::Abs(iz-iiz)>0.75*idz) continue; |
| | 639 | } |
| | 640 | ndigits++; |
| | 641 | Float_t qBin = bins[idxBins[idx]].GetQ(); |
| | 642 | y += qBin * fYSPD[iy]; |
| | 643 | z += qBin * fZSPD[iz]; |
| | 644 | q += qBin; |
| | 645 | } |
| | 646 | if (ndigits==0) continue; |
| | 647 | y /= q; |
| | 648 | z /= q; |
| | 649 | y -= fHwSPD; |
| | 650 | z -= fHlSPD; |
| | 651 | |
| | 652 | Float_t hit[5]; // y, z, sigma(y)^2, sigma(z)^2, charge |
| | 653 | hit[0] = -(-y+fYshift[iModule]); |
| | 654 | if (iModule <= fLastSPD1) hit[0] = -hit[0]; |
| | 655 | hit[1] = -z+fZshift[iModule]; |
| | 656 | hit[2] = fYpitchSPD*fYpitchSPD/12.; |
| | 657 | hit[3] = fZ1pitchSPD*fZ1pitchSPD/12.; |
| | 658 | // hit[4] = q; |
| | 659 | hit[4] = (zmax-zmin+1)*100 + (ymax-ymin+1); |
| | 660 | // CheckLabels(label); |
| | 661 | Int_t info[3]={ymax-ymin+1,zmax-zmin+1,fNlayer[iModule]}; |
| | 662 | new (clusters[iModule]->AddrAt(nClusters)) |
| | 663 | AliITSclusterV2(label, hit,info); |
| | 664 | nClusters++; |
| | 665 | } |
| | 666 | } |
| | 667 | |
| | 668 | nClustersSPD += nClusters; |
| | 669 | bins = NULL; |
| | 670 | } |
| | 671 | |
| | 672 | if (!next) break; |
| | 673 | bins = binsSPD; |
| | 674 | memcpy(binsSPD,binsSPDInit,sizeof(AliBin)*kMaxBin); |
| | 675 | } |
| | 676 | |
| | 677 | // fill the current digit into the bins array |
| | 678 | Int_t index = (input->GetCoord2()+1) * kNzBins + (input->GetCoord1()+1); |
| | 679 | bins[index].SetIndex(index); |
| | 680 | bins[index].SetMask(1); |
| | 681 | bins[index].SetQ(1); |
| | 682 | } |
| | 683 | |
| | 684 | delete [] binsSPDInit; |
| | 685 | delete [] binsSPD; |
| | 686 | |
| | 687 | Info("FindClustersSPD", "found clusters in ITS SPD: %d", nClustersSPD); |
| | 688 | } |
| | 689 | |
| | 690 | |
| | 691 | Bool_t AliITSclustererV2::IsMaximum(Int_t k,Int_t max,const AliBin *bins) { |
| | 692 | //------------------------------------------------------------ |
| | 693 | //is this a local maximum ? |
| | 694 | //------------------------------------------------------------ |
| | 695 | UShort_t q=bins[k].GetQ(); |
| | 696 | if (q==1023) return kFALSE; |
| | 697 | if (bins[k-max].GetQ() > q) return kFALSE; |
| | 698 | if (bins[k-1 ].GetQ() > q) return kFALSE; |
| | 699 | if (bins[k+max].GetQ() > q) return kFALSE; |
| | 700 | if (bins[k+1 ].GetQ() > q) return kFALSE; |
| | 701 | if (bins[k-max-1].GetQ() > q) return kFALSE; |
| | 702 | if (bins[k+max-1].GetQ() > q) return kFALSE; |
| | 703 | if (bins[k+max+1].GetQ() > q) return kFALSE; |
| | 704 | if (bins[k-max+1].GetQ() > q) return kFALSE; |
| | 705 | return kTRUE; |
| | 706 | } |
| | 707 | |
| | 708 | void AliITSclustererV2:: |
| | 709 | FindPeaks(Int_t k,Int_t max,AliBin *b,Int_t *idx,UInt_t *msk,Int_t& n) { |
| | 710 | //------------------------------------------------------------ |
| | 711 | //find local maxima |
| | 712 | //------------------------------------------------------------ |
| | 713 | if (n<31) |
| | 714 | if (IsMaximum(k,max,b)) { |
| | 715 | idx[n]=k; msk[n]=(2<<n); |
| | 716 | n++; |
| | 717 | } |
| | 718 | b[k].SetMask(0); |
| | 719 | if (b[k-max].GetMask()&1) FindPeaks(k-max,max,b,idx,msk,n); |
| | 720 | if (b[k-1 ].GetMask()&1) FindPeaks(k-1 ,max,b,idx,msk,n); |
| | 721 | if (b[k+max].GetMask()&1) FindPeaks(k+max,max,b,idx,msk,n); |
| | 722 | if (b[k+1 ].GetMask()&1) FindPeaks(k+1 ,max,b,idx,msk,n); |
| | 723 | } |
| | 724 | |
| | 725 | void AliITSclustererV2:: |
| | 726 | MarkPeak(Int_t k, Int_t max, AliBin *bins, UInt_t m) { |
| | 727 | //------------------------------------------------------------ |
| | 728 | //mark this peak |
| | 729 | //------------------------------------------------------------ |
| | 730 | UShort_t q=bins[k].GetQ(); |
| | 731 | |
| | 732 | bins[k].SetMask(bins[k].GetMask()|m); |
| | 733 | |
| | 734 | if (bins[k-max].GetQ() <= q) |
| | 735 | if ((bins[k-max].GetMask()&m) == 0) MarkPeak(k-max,max,bins,m); |
| | 736 | if (bins[k-1 ].GetQ() <= q) |
| | 737 | if ((bins[k-1 ].GetMask()&m) == 0) MarkPeak(k-1 ,max,bins,m); |
| | 738 | if (bins[k+max].GetQ() <= q) |
| | 739 | if ((bins[k+max].GetMask()&m) == 0) MarkPeak(k+max,max,bins,m); |
| | 740 | if (bins[k+1 ].GetQ() <= q) |
| | 741 | if ((bins[k+1 ].GetMask()&m) == 0) MarkPeak(k+1 ,max,bins,m); |
| | 742 | } |
| | 743 | |
| | 744 | void AliITSclustererV2:: |
| | 745 | MakeCluster(Int_t k,Int_t max,AliBin *bins,UInt_t m,AliITSclusterV2 &c) { |
| | 746 | //------------------------------------------------------------ |
| | 747 | //make cluster using digits of this peak |
| | 748 | //------------------------------------------------------------ |
| | 749 | Float_t q=(Float_t)bins[k].GetQ(); |
| | 750 | Int_t i=k/max, j=k-i*max; |
| | 751 | |
| | 752 | c.SetQ(c.GetQ()+q); |
| | 753 | c.SetY(c.GetY()+i*q); |
| | 754 | c.SetZ(c.GetZ()+j*q); |
| | 755 | c.SetSigmaY2(c.GetSigmaY2()+i*i*q); |
| | 756 | c.SetSigmaZ2(c.GetSigmaZ2()+j*j*q); |
| | 757 | |
| | 758 | bins[k].SetMask(0xFFFFFFFE); |
| | 759 | |
| | 760 | if (bins[k-max].GetMask() == m) MakeCluster(k-max,max,bins,m,c); |
| | 761 | if (bins[k-1 ].GetMask() == m) MakeCluster(k-1 ,max,bins,m,c); |
| | 762 | if (bins[k+max].GetMask() == m) MakeCluster(k+max,max,bins,m,c); |
| | 763 | if (bins[k+1 ].GetMask() == m) MakeCluster(k+1 ,max,bins,m,c); |
| | 764 | } |
| | 765 | |
| | 766 | void AliITSclustererV2:: |
| | 767 | FindClustersSDD(AliBin* bins[2], Int_t nMaxBin, Int_t nzBins, |
| | 768 | const TClonesArray *digits, TClonesArray *clusters) { |
| | 769 | //------------------------------------------------------------ |
| | 770 | // Actual SDD cluster finder |
| | 771 | //------------------------------------------------------------ |
| | 772 | Int_t ncl=0; TClonesArray &cl=*clusters; |
| | 773 | for (Int_t s=0; s<2; s++) |
| | 774 | for (Int_t i=0; i<nMaxBin; i++) { |
| | 775 | if (bins[s][i].IsUsed()) continue; |
| | 776 | Int_t idx[32]; UInt_t msk[32]; Int_t npeaks=0; |
| | 777 | FindPeaks(i, nzBins, bins[s], idx, msk, npeaks); |
| | 778 | |
| | 779 | if (npeaks>30) continue; |
| | 780 | if (npeaks==0) continue; |
| | 781 | |
| | 782 | Int_t k,l; |
| | 783 | for (k=0; k<npeaks-1; k++){//mark adjacent peaks |
| | 784 | if (idx[k] < 0) continue; //this peak is already removed |
| | 785 | for (l=k+1; l<npeaks; l++) { |
| | 786 | if (idx[l] < 0) continue; //this peak is already removed |
| | 787 | Int_t ki=idx[k]/nzBins, kj=idx[k] - ki*nzBins; |
| | 788 | Int_t li=idx[l]/nzBins, lj=idx[l] - li*nzBins; |
| | 789 | Int_t di=TMath::Abs(ki - li); |
| | 790 | Int_t dj=TMath::Abs(kj - lj); |
| | 791 | if (di>1 || dj>1) continue; |
| | 792 | if (bins[s][idx[k]].GetQ() > bins[s][idx[l]].GetQ()) { |
| | 793 | msk[l]=msk[k]; |
| | 794 | idx[l]*=-1; |
| | 795 | } else { |
| | 796 | msk[k]=msk[l]; |
| | 797 | idx[k]*=-1; |
| | 798 | break; |
| | 799 | } |
| | 800 | } |
| | 801 | } |
| | 802 | |
| | 803 | for (k=0; k<npeaks; k++) { |
| | 804 | MarkPeak(TMath::Abs(idx[k]), nzBins, bins[s], msk[k]); |
| | 805 | } |
| | 806 | |
| | 807 | for (k=0; k<npeaks; k++) { |
| | 808 | if (idx[k] < 0) continue; //removed peak |
| | 809 | AliITSclusterV2 c; |
| | 810 | MakeCluster(idx[k], nzBins, bins[s], msk[k], c); |
| | 811 | //mi change |
| | 812 | Int_t milab[10]; |
| | 813 | for (Int_t ilab=0;ilab<10;ilab++){ |
| | 814 | milab[ilab]=-2; |
| | 815 | } |
| | 816 | Int_t maxi=0,mini=0,maxj=0,minj=0; |
| | 817 | //AliBin *bmax=&bins[s][idx[k]]; |
| | 818 | //Float_t max = TMath::Max(TMath::Abs(bmax->GetQ())/5.,3.); |
| | 819 | Float_t max=3; |
| | 820 | for (Int_t di=-2; di<=2;di++) |
| | 821 | for (Int_t dj=-3;dj<=3;dj++){ |
| | 822 | Int_t index = idx[k]+di+dj*nzBins; |
| | 823 | if (index<0) continue; |
| | 824 | if (index>=nMaxBin) continue; |
| | 825 | AliBin *b=&bins[s][index]; |
| | 826 | if (TMath::Abs(b->GetQ())>max){ |
| | 827 | if (di>maxi) maxi=di; |
| | 828 | if (di<mini) mini=di; |
| | 829 | if (dj>maxj) maxj=dj; |
| | 830 | if (dj<minj) minj=dj; |
| | 831 | // |
| | 832 | if(digits) { |
| | 833 | if (TMath::Abs(di)<2&&TMath::Abs(dj)<2){ |
| | 834 | AliITSdigitSDD* d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); |
| | 835 | for (Int_t itrack=0;itrack<10;itrack++){ |
| | 836 | Int_t track = (d->GetTracks())[itrack]; |
| | 837 | if (track>=0) { |
| | 838 | AddLabel(milab, track); |
| | 839 | } |
| | 840 | } |
| | 841 | } |
| | 842 | } |
| | 843 | } |
| | 844 | } |
| | 845 | |
| | 846 | /* |
| | 847 | Float_t s2 = c.GetSigmaY2()/c.GetQ() - c.GetY()*c.GetY(); |
| | 848 | Float_t w=par->GetPadPitchWidth(sec); |
| | 849 | c.SetSigmaY2(s2); |
| | 850 | if (s2 != 0.) { |
| | 851 | c.SetSigmaY2(c.GetSigmaY2()*0.108); |
| | 852 | if (sec<par->GetNInnerSector()) c.SetSigmaY2(c.GetSigmaY2()*2.07); |
| | 853 | } |
| | 854 | s2 = c.GetSigmaZ2()/c.GetQ() - c.GetZ()*c.GetZ(); |
| | 855 | w=par->GetZWidth(); |
| | 856 | c.SetSigmaZ2(s2); |
| | 857 | |
| | 858 | if (s2 != 0.) { |
| | 859 | c.SetSigmaZ2(c.GetSigmaZ2()*0.169); |
| | 860 | if (sec<par->GetNInnerSector()) c.SetSigmaZ2(c.GetSigmaZ2()*1.77); |
| | 861 | } |
| | 862 | */ |
| | 863 | |
| | 864 | Float_t y=c.GetY(),z=c.GetZ(), q=c.GetQ(); |
| | 865 | y/=q; z/=q; |
| | 866 | // |
| | 867 | //Float_t s2 = c.GetSigmaY2()/c.GetQ() - y*y; |
| | 868 | // c.SetSigmaY2(s2); |
| | 869 | //s2 = c.GetSigmaZ2()/c.GetQ() - z*z; |
| | 870 | //c.SetSigmaZ2(s2); |
| | 871 | // |
| | 872 | y=(y-0.5)*fYpitchSDD; |
| | 873 | y-=fHwSDD; |
| | 874 | y-=fYoffSDD; //delay ? |
| | 875 | if (s) y=-y; |
| | 876 | |
| | 877 | z=(z-0.5)*fZpitchSDD; |
| | 878 | z-=fHlSDD; |
| | 879 | |
| | 880 | y=-(-y+fYshift[fI]); |
| | 881 | z= -z+fZshift[fI]; |
| | 882 | |
| | 883 | q/=12.7; //this WAS consistent with SSD. To be reassessed |
| | 884 | // 23-MAR-2007 |
| | 885 | Float_t hit[5] = {y, z, 0.0030*0.0030, 0.0020*0.0020, q}; |
| | 886 | Int_t info[3] = {maxj-minj+1, maxi-mini+1, fNlayer[fI]}; |
| | 887 | |
| | 888 | if (c.GetQ() < 20.) continue; //noise cluster |
| | 889 | |
| | 890 | if (digits) { |
| | 891 | // AliBin *b=&bins[s][idx[k]]; |
| | 892 | // AliITSdigitSDD* d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); |
| | 893 | { |
| | 894 | //Int_t lab[3]; |
| | 895 | //lab[0]=(d->GetTracks())[0]; |
| | 896 | //lab[1]=(d->GetTracks())[1]; |
| | 897 | //lab[2]=(d->GetTracks())[2]; |
| | 898 | //CheckLabels(lab); |
| | 899 | CheckLabels2(milab); |
| | 900 | } |
| | 901 | } |
| | 902 | milab[3]=fNdet[fI]; |
| | 903 | |
| | 904 | AliITSclusterV2 *cc = new (cl[ncl]) AliITSclusterV2(milab,hit,info); |
| | 905 | cc->SetType(npeaks); |
| | 906 | ncl++; |
| | 907 | } |
| | 908 | } |
| | 909 | } |
| | 910 | |
| | 911 | void AliITSclustererV2:: |
| | 912 | FindClustersSDD(const TClonesArray *digits, TClonesArray *clusters) { |
| | 913 | //------------------------------------------------------------ |
| | 914 | // Actual SDD cluster finder |
| | 915 | //------------------------------------------------------------ |
| | 916 | Int_t kNzBins = fNzSDD + 2; |
| | 917 | const Int_t kMAXBIN=kNzBins*(fNySDD+2); |
| | 918 | |
| | 919 | AliBin *bins[2]; |
| | 920 | bins[0]=new AliBin[kMAXBIN]; |
| | 921 | bins[1]=new AliBin[kMAXBIN]; |
| | 922 | |
| | 923 | AliITSdigitSDD *d=0; |
| | 924 | Int_t i, ndigits=digits->GetEntriesFast(); |
| | 925 | for (i=0; i<ndigits; i++) { |
| | 926 | d=(AliITSdigitSDD*)digits->UncheckedAt(i); |
| | 927 | Int_t y=d->GetCoord2()+1; //y |
| | 928 | Int_t z=d->GetCoord1()+1; //z |
| | 929 | Int_t q=d->GetSignal(); |
| | 930 | if (q<3) continue; |
| | 931 | |
| | 932 | if (z <= fNzSDD) { |
| | 933 | bins[0][y*kNzBins+z].SetQ(q); |
| | 934 | bins[0][y*kNzBins+z].SetMask(1); |
| | 935 | bins[0][y*kNzBins+z].SetIndex(i); |
| | 936 | } else { |
| | 937 | z-=fNzSDD; |
| | 938 | bins[1][y*kNzBins+z].SetQ(q); |
| | 939 | bins[1][y*kNzBins+z].SetMask(1); |
| | 940 | bins[1][y*kNzBins+z].SetIndex(i); |
| | 941 | } |
| | 942 | } |
| | 943 | |
| | 944 | FindClustersSDD(bins, kMAXBIN, kNzBins, digits, clusters); |
| | 945 | |
| | 946 | delete[] bins[0]; |
| | 947 | delete[] bins[1]; |
| | 948 | |
| | 949 | } |
| | 950 | |
| | 951 | void AliITSclustererV2::FindClustersSDD(AliITSRawStream* input, |
| | 952 | TClonesArray** clusters) |
| | 953 | { |
| | 954 | //------------------------------------------------------------ |
| | 955 | // Actual SDD cluster finder for raw data |
| | 956 | //------------------------------------------------------------ |
| | 957 | Int_t nClustersSDD = 0; |
| | 958 | Int_t kNzBins = fNzSDD + 2; |
| | 959 | Int_t kMaxBin = kNzBins * (fNySDD+2); |
| | 960 | AliBin *binsSDDInit = new AliBin[kMaxBin]; |
| | 961 | AliBin *binsSDD1 = new AliBin[kMaxBin]; |
| | 962 | AliBin *binsSDD2 = new AliBin[kMaxBin]; |
| | 963 | AliBin *bins[2] = {NULL, NULL}; |
| | 964 | |
| | 965 | // read raw data input stream |
| | 966 | while (kTRUE) { |
| | 967 | Bool_t next = input->Next(); |
| | 968 | if (!next || input->IsNewModule()) { |
| | 969 | Int_t iModule = input->GetPrevModuleID(); |
| | 970 | |
| | 971 | // when all data from a module was read, search for clusters |
| | 972 | if (bins[0]) { |
| | 973 | clusters[iModule] = new TClonesArray("AliITSclusterV2"); |
| | 974 | fI = iModule; |
| | 975 | FindClustersSDD(bins, kMaxBin, kNzBins, NULL, clusters[iModule]); |
| | 976 | Int_t nClusters = clusters[iModule]->GetEntriesFast(); |
| | 977 | nClustersSDD += nClusters; |
| | 978 | bins[0] = bins[1] = NULL; |
| | 979 | } |
| | 980 | |
| | 981 | if (!next) break; |
| | 982 | bins[0]=binsSDD1; |
| | 983 | bins[1]=binsSDD2; |
| | 984 | memcpy(binsSDD1,binsSDDInit,sizeof(AliBin)*kMaxBin); |
| | 985 | memcpy(binsSDD2,binsSDDInit,sizeof(AliBin)*kMaxBin); |
| | 986 | } |
| | 987 | |
| | 988 | // fill the current digit into the bins array |
| | 989 | if(input->GetSignal()>=3) { |
| | 990 | Int_t iz = input->GetCoord1()+1; |
| | 991 | Int_t side = ((iz <= fNzSDD) ? 0 : 1); |
| | 992 | iz -= side*fNzSDD; |
| | 993 | Int_t index = (input->GetCoord2()+1) * kNzBins + iz; |
| | 994 | bins[side][index].SetQ(input->GetSignal()); |
| | 995 | bins[side][index].SetMask(1); |
| | 996 | bins[side][index].SetIndex(index); |
| | 997 | } |
| | 998 | } |
| | 999 | |
| | 1000 | delete[] binsSDD1; |
| | 1001 | delete[] binsSDD2; |
| | 1002 | delete[] binsSDDInit; |
| | 1003 | |
| | 1004 | Info("FindClustersSDD", "found clusters in ITS SDD: %d", nClustersSDD); |
| | 1005 | } |
| | 1006 | |
| | 1007 | |
| | 1008 | |
| | 1009 | void AliITSclustererV2:: |
| | 1010 | FindClustersSSD(Ali1Dcluster* neg, Int_t nn, |
| | 1011 | Ali1Dcluster* pos, Int_t np, |
| | 1012 | TClonesArray *clusters) { |
| | 1013 | //------------------------------------------------------------ |
| | 1014 | // Actual SSD cluster finder |
| | 1015 | //------------------------------------------------------------ |
| | 1016 | TClonesArray &cl=*clusters; |
| | 1017 | // |
| | 1018 | Float_t tanp=fTanP, tann=fTanN; |
| | 1019 | if (fI>fLastSSD1) {tann=fTanP; tanp=fTanN;} |
| | 1020 | Int_t idet=fNdet[fI]; |
| | 1021 | Int_t ncl=0; |
| | 1022 | // |
| | 1023 | Int_t negativepair[30000]; |
| | 1024 | Int_t cnegative[3000]; |
| | 1025 | Int_t cused1[3000]; |
| | 1026 | Int_t positivepair[30000]; |
| | 1027 | Int_t cpositive[3000]; |
| | 1028 | Int_t cused2[3000]; |
| | 1029 | for (Int_t i=0;i<3000;i++) {cnegative[i]=0; cused1[i]=0;} |
| | 1030 | for (Int_t i=0;i<3000;i++) {cpositive[i]=0; cused2[i]=0;} |
| | 1031 | for (Int_t i=0;i<30000;i++){negativepair[i]=positivepair[i]=0;} |
| | 1032 | static Short_t pairs[1000][1000]; |
| | 1033 | memset(pairs,0,sizeof(Short_t)*1000000); |
| | 1034 | // Short_t ** pairs = new Short_t*[1000]; |
| | 1035 | // for (Int_t i=0; i<1000; i++) { |
| | 1036 | // pairs[i] = new Short_t[1000]; |
| | 1037 | // memset(pairs[i],0,sizeof(Short_t)*1000); |
| | 1038 | // } |
| | 1039 | // |
| | 1040 | // find available pairs |
| | 1041 | // |
| | 1042 | for (Int_t i=0; i<np; i++) { |
| | 1043 | Float_t yp=pos[i].GetY()*fYpitchSSD; |
| | 1044 | if (pos[i].GetQ()<3) continue; |
| | 1045 | for (Int_t j=0; j<nn; j++) { |
| | 1046 | if (neg[j].GetQ()<3) continue; |
| | 1047 | Float_t yn=neg[j].GetY()*fYpitchSSD; |
| | 1048 | Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp); |
| | 1049 | Float_t yt=yn + tann*zt; |
| | 1050 | zt-=fHlSSD; yt-=fHwSSD; |
| | 1051 | if (TMath::Abs(yt)<fHwSSD+0.01) |
| | 1052 | if (TMath::Abs(zt)<fHlSSD+0.01*(neg[j].GetNd()+pos[i].GetNd())) { |
| | 1053 | negativepair[i*10+cnegative[i]] =j; //index |
| | 1054 | positivepair[j*10+cpositive[j]] =i; |
| | 1055 | cnegative[i]++; //counters |
| | 1056 | cpositive[j]++; |
| | 1057 | pairs[i][j]=100; |
| | 1058 | } |
| | 1059 | } |
| | 1060 | } |
| | 1061 | // |
| | 1062 | for (Int_t i=0; i<np; i++) { |
| | 1063 | Float_t yp=pos[i].GetY()*fYpitchSSD; |
| | 1064 | if (pos[i].GetQ()<3) continue; |
| | 1065 | for (Int_t j=0; j<nn; j++) { |
| | 1066 | if (neg[j].GetQ()<3) continue; |
| | 1067 | if (cpositive[j]&&cnegative[i]) continue; |
| | 1068 | Float_t yn=neg[j].GetY()*fYpitchSSD; |
| | 1069 | Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp); |
| | 1070 | Float_t yt=yn + tann*zt; |
| | 1071 | zt-=fHlSSD; yt-=fHwSSD; |
| | 1072 | if (TMath::Abs(yt)<fHwSSD+0.1) |
| | 1073 | if (TMath::Abs(zt)<fHlSSD+0.15) { |
| | 1074 | if (cnegative[i]==0) pos[i].SetNd(100); // not available pair |
| | 1075 | if (cpositive[j]==0) neg[j].SetNd(100); // not available pair |
| | 1076 | negativepair[i*10+cnegative[i]] =j; //index |
| | 1077 | positivepair[j*10+cpositive[j]] =i; |
| | 1078 | cnegative[i]++; //counters |
| | 1079 | cpositive[j]++; |
| | 1080 | pairs[i][j]=100; |
| | 1081 | } |
| | 1082 | } |
| | 1083 | } |
| | 1084 | // |
| | 1085 | Float_t lp[5]; |
| | 1086 | Int_t milab[10]; |
| | 1087 | Double_t ratio; |
| | 1088 | |
| | 1089 | // |
| | 1090 | // sign gold tracks |
| | 1091 | // |
| | 1092 | for (Int_t ip=0;ip<np;ip++){ |
| | 1093 | Float_t ybest=1000,zbest=1000,qbest=0; |
| | 1094 | // |
| | 1095 | // select gold clusters |
| | 1096 | if ( (cnegative[ip]==1) && cpositive[negativepair[10*ip]]==1){ |
| | 1097 | Float_t yp=pos[ip].GetY()*fYpitchSSD; |
| | 1098 | Int_t j = negativepair[10*ip]; |
| | 1099 | ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ()); |
| | 1100 | // |
| | 1101 | Float_t yn=neg[j].GetY()*fYpitchSSD; |
| | 1102 | Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp); |
| | 1103 | Float_t yt=yn + tann*zt; |
| | 1104 | zt-=fHlSSD; yt-=fHwSSD; |
| | 1105 | ybest=yt; zbest=zt; |
| | 1106 | qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ()); |
| | 1107 | lp[0]=-(-ybest+fYshift[fI]); |
| | 1108 | lp[1]= -zbest+fZshift[fI]; |
| | 1109 | lp[2]=0.0025*0.0025; //SigmaY2 |
| | 1110 | lp[3]=0.110*0.110; //SigmaZ2 |
| | 1111 | |
| | 1112 | lp[4]=qbest; //Q |
| | 1113 | for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2; |
| | 1114 | for (Int_t ilab=0;ilab<3;ilab++){ |
| | 1115 | milab[ilab] = pos[ip].GetLabel(ilab); |
| | 1116 | milab[ilab+3] = neg[j].GetLabel(ilab); |
| | 1117 | } |
| | 1118 | // |
| | 1119 | CheckLabels2(milab); |
| | 1120 | milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det |
| | 1121 | Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fI]}; |
| | 1122 | AliITSclusterV2 * cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info); |
| | 1123 | ncl++; |
| | 1124 | cl2->SetChargeRatio(ratio); |
| | 1125 | cl2->SetType(1); |
| | 1126 | pairs[ip][j]=1; |
| | 1127 | if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster |
| | 1128 | cl2->SetType(2); |
| | 1129 | pairs[ip][j]=2; |
| | 1130 | } |
| | 1131 | cused1[ip]++; |
| | 1132 | cused2[j]++; |
| | 1133 | } |
| | 1134 | } |
| | 1135 | |
| | 1136 | for (Int_t ip=0;ip<np;ip++){ |
| | 1137 | Float_t ybest=1000,zbest=1000,qbest=0; |
| | 1138 | // |
| | 1139 | // |
| | 1140 | // select "silber" cluster |
| | 1141 | if ( cnegative[ip]==1 && cpositive[negativepair[10*ip]]==2){ |
| | 1142 | Int_t in = negativepair[10*ip]; |
| | 1143 | Int_t ip2 = positivepair[10*in]; |
| | 1144 | if (ip2==ip) ip2 = positivepair[10*in+1]; |
| | 1145 | Float_t pcharge = pos[ip].GetQ()+pos[ip2].GetQ(); |
| | 1146 | if (TMath::Abs(pcharge-neg[in].GetQ())<10){ |
| | 1147 | // |
| | 1148 | // add first pair |
| | 1149 | if (pairs[ip][in]==100){ // |
| | 1150 | Float_t yp=pos[ip].GetY()*fYpitchSSD; |
| | 1151 | Float_t yn=neg[in].GetY()*fYpitchSSD; |
| | 1152 | Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp); |
| | 1153 | Float_t yt=yn + tann*zt; |
| | 1154 | zt-=fHlSSD; yt-=fHwSSD; |
| | 1155 | ybest =yt; zbest=zt; |
| | 1156 | qbest =pos[ip].GetQ(); |
| | 1157 | lp[0]=-(-ybest+fYshift[fI]); |
| | 1158 | lp[1]= -zbest+fZshift[fI]; |
| | 1159 | lp[2]=0.0025*0.0025; //SigmaY2 |
| | 1160 | lp[3]=0.110*0.110; //SigmaZ2 |
| | 1161 | |
| | 1162 | lp[4]=qbest; //Q |
| | 1163 | for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2; |
| | 1164 | for (Int_t ilab=0;ilab<3;ilab++){ |
| | 1165 | milab[ilab] = pos[ip].GetLabel(ilab); |
| | 1166 | milab[ilab+3] = neg[in].GetLabel(ilab); |
| | 1167 | } |
| | 1168 | // |
| | 1169 | CheckLabels2(milab); |
| | 1170 | ratio = (pos[ip].GetQ()-neg[in].GetQ())/(pos[ip].GetQ()+neg[in].GetQ()); |
| | 1171 | milab[3]=(((ip<<10) + in)<<10) + idet; // pos|neg|det |
| | 1172 | Int_t info[3] = {pos[ip].GetNd(),neg[in].GetNd(),fNlayer[fI]}; |
| | 1173 | AliITSclusterV2 * cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info); |
| | 1174 | ncl++; |
| | 1175 | cl2->SetChargeRatio(ratio); |
| | 1176 | cl2->SetType(5); |
| | 1177 | pairs[ip][in] = 5; |
| | 1178 | if ((pos[ip].GetNd()+neg[in].GetNd())>6){ //multi cluster |
| | 1179 | cl2->SetType(6); |
| | 1180 | pairs[ip][in] = 6; |
| | 1181 | } |
| | 1182 | } |
| | 1183 | // |
| | 1184 | // add second pair |
| | 1185 | |
| | 1186 | // if (!(cused1[ip2] || cused2[in])){ // |
| | 1187 | if (pairs[ip2][in]==100){ |
| | 1188 | Float_t yp=pos[ip2].GetY()*fYpitchSSD; |
| | 1189 | Float_t yn=neg[in].GetY()*fYpitchSSD; |
| | 1190 | Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp); |
| | 1191 | Float_t yt=yn + tann*zt; |
| | 1192 | zt-=fHlSSD; yt-=fHwSSD; |
| | 1193 | ybest =yt; zbest=zt; |
| | 1194 | qbest =pos[ip2].GetQ(); |
| | 1195 | lp[0]=-(-ybest+fYshift[fI]); |
| | 1196 | lp[1]= -zbest+fZshift[fI]; |
| | 1197 | lp[2]=0.0025*0.0025; //SigmaY2 |
| | 1198 | lp[3]=0.110*0.110; //SigmaZ2 |
| | 1199 | |
| | 1200 | lp[4]=qbest; //Q |
| | 1201 | for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2; |
| | 1202 | for (Int_t ilab=0;ilab<3;ilab++){ |
| | 1203 | milab[ilab] = pos[ip2].GetLabel(ilab); |
| | 1204 | milab[ilab+3] = neg[in].GetLabel(ilab); |
| | 1205 | } |
| | 1206 | // |
| | 1207 | CheckLabels2(milab); |
| | 1208 | ratio = (pos[ip2].GetQ()-neg[in].GetQ())/(pos[ip2].GetQ()+neg[in].GetQ()); |
| | 1209 | milab[3]=(((ip2<<10) + in)<<10) + idet; // pos|neg|det |
| | 1210 | Int_t info[3] = {pos[ip2].GetNd(),neg[in].GetNd(),fNlayer[fI]}; |
| | 1211 | AliITSclusterV2 *cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info); |
| | 1212 | ncl++; |
| | 1213 | cl2->SetChargeRatio(ratio); |
| | 1214 | cl2->SetType(5); |
| | 1215 | pairs[ip2][in] =5; |
| | 1216 | if ((pos[ip2].GetNd()+neg[in].GetNd())>6){ //multi cluster |
| | 1217 | cl2->SetType(6); |
| | 1218 | pairs[ip2][in] =6; |
| | 1219 | } |
| | 1220 | } |
| | 1221 | cused1[ip]++; |
| | 1222 | cused1[ip2]++; |
| | 1223 | cused2[in]++; |
| | 1224 | } |
| | 1225 | } |
| | 1226 | } |
| | 1227 | |
| | 1228 | // |
| | 1229 | for (Int_t jn=0;jn<nn;jn++){ |
| | 1230 | if (cused2[jn]) continue; |
| | 1231 | Float_t ybest=1000,zbest=1000,qbest=0; |
| | 1232 | // select "silber" cluster |
| | 1233 | if ( cpositive[jn]==1 && cnegative[positivepair[10*jn]]==2){ |
| | 1234 | Int_t ip = positivepair[10*jn]; |
| | 1235 | Int_t jn2 = negativepair[10*ip]; |
| | 1236 | if (jn2==jn) jn2 = negativepair[10*ip+1]; |
| | 1237 | Float_t pcharge = neg[jn].GetQ()+neg[jn2].GetQ(); |
| | 1238 | // |
| | 1239 | if (TMath::Abs(pcharge-pos[ip].GetQ())<10){ |
| | 1240 | // |
| | 1241 | // add first pair |
| | 1242 | // if (!(cused1[ip]||cused2[jn])){ |
| | 1243 | if (pairs[ip][jn]==100){ |
| | 1244 | Float_t yn=neg[jn].GetY()*fYpitchSSD; |
| | 1245 | Float_t yp=pos[ip].GetY()*fYpitchSSD; |
| | 1246 | Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp); |
| | 1247 | Float_t yt=yn + tann*zt; |
| | 1248 | zt-=fHlSSD; yt-=fHwSSD; |
| | 1249 | ybest =yt; zbest=zt; |
| | 1250 | qbest =neg[jn].GetQ(); |
| | 1251 | lp[0]=-(-ybest+fYshift[fI]); |
| | 1252 | lp[1]= -zbest+fZshift[fI]; |
| | 1253 | lp[2]=0.0025*0.0025; //SigmaY2 |
| | 1254 | lp[3]=0.110*0.110; //SigmaZ2 |
| | 1255 | |
| | 1256 | lp[4]=qbest; //Q |
| | 1257 | for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2; |
| | 1258 | for (Int_t ilab=0;ilab<3;ilab++){ |
| | 1259 | milab[ilab] = pos[ip].GetLabel(ilab); |
| | 1260 | milab[ilab+3] = neg[jn].GetLabel(ilab); |
| | 1261 | } |
| | 1262 | // |
| | 1263 | CheckLabels2(milab); |
| | 1264 | ratio = (pos[ip].GetQ()-neg[jn].GetQ())/(pos[ip].GetQ()+neg[jn].GetQ()); |
| | 1265 | milab[3]=(((ip<<10) + jn)<<10) + idet; // pos|neg|det |
| | 1266 | Int_t info[3] = {pos[ip].GetNd(),neg[jn].GetNd(),fNlayer[fI]}; |
| | 1267 | AliITSclusterV2 * cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info); |
| | 1268 | ncl++; |
| | 1269 | cl2->SetChargeRatio(ratio); |
| | 1270 | cl2->SetType(7); |
| | 1271 | pairs[ip][jn] =7; |
| | 1272 | if ((pos[ip].GetNd()+neg[jn].GetNd())>6){ //multi cluster |
| | 1273 | cl2->SetType(8); |
| | 1274 | pairs[ip][jn]=8; |
| | 1275 | } |
| | 1276 | } |
| | 1277 | // |
| | 1278 | // add second pair |
| | 1279 | // if (!(cused1[ip]||cused2[jn2])){ |
| | 1280 | if (pairs[ip][jn2]==100){ |
| | 1281 | Float_t yn=neg[jn2].GetY()*fYpitchSSD; |
| | 1282 | Double_t yp=pos[ip].GetY()*fYpitchSSD; |
| | 1283 | Double_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp); |
| | 1284 | Double_t yt=yn + tann*zt; |
| | 1285 | zt-=fHlSSD; yt-=fHwSSD; |
| | 1286 | ybest =yt; zbest=zt; |
| | 1287 | qbest =neg[jn2].GetQ(); |
| | 1288 | lp[0]=-(-ybest+fYshift[fI]); |
| | 1289 | lp[1]= -zbest+fZshift[fI]; |
| | 1290 | lp[2]=0.0025*0.0025; //SigmaY2 |
| | 1291 | lp[3]=0.110*0.110; //SigmaZ2 |
| | 1292 | |
| | 1293 | lp[4]=qbest; //Q |
| | 1294 | for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2; |
| | 1295 | for (Int_t ilab=0;ilab<3;ilab++){ |
| | 1296 | milab[ilab] = pos[ip].GetLabel(ilab); |
| | 1297 | milab[ilab+3] = neg[jn2].GetLabel(ilab); |
| | 1298 | } |
| | 1299 | // |
| | 1300 | CheckLabels2(milab); |
| | 1301 | ratio = (pos[ip].GetQ()-neg[jn2].GetQ())/(pos[ip].GetQ()+neg[jn2].GetQ()); |
| | 1302 | milab[3]=(((ip<<10) + jn2)<<10) + idet; // pos|neg|det |
| | 1303 | Int_t info[3] = {pos[ip].GetNd(),neg[jn2].GetNd(),fNlayer[fI]}; |
| | 1304 | AliITSclusterV2* cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info); |
| | 1305 | ncl++; |
| | 1306 | cl2->SetChargeRatio(ratio); |
| | 1307 | pairs[ip][jn2]=7; |
| | 1308 | cl2->SetType(7); |
| | 1309 | if ((pos[ip].GetNd()+neg[jn2].GetNd())>6){ //multi cluster |
| | 1310 | cl2->SetType(8); |
| | 1311 | pairs[ip][jn2]=8; |
| | 1312 | } |
| | 1313 | } |
| | 1314 | cused1[ip]++; |
| | 1315 | cused2[jn]++; |
| | 1316 | cused2[jn2]++; |
| | 1317 | } |
| | 1318 | } |
| | 1319 | } |
| | 1320 | |
| | 1321 | for (Int_t ip=0;ip<np;ip++){ |
| | 1322 | Float_t ybest=1000,zbest=1000,qbest=0; |
| | 1323 | // |
| | 1324 | // 2x2 clusters |
| | 1325 | // |
| | 1326 | if ( (cnegative[ip]<5) && cpositive[negativepair[10*ip]]<5){ |
| | 1327 | Float_t minchargediff =4.; |
| | 1328 | Int_t j=-1; |
| | 1329 | for (Int_t di=0;di<cnegative[ip];di++){ |
| | 1330 | Int_t jc = negativepair[ip*10+di]; |
| | 1331 | Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ(); |
| | 1332 | if (TMath::Abs(chargedif)<minchargediff){ |
| | 1333 | j =jc; |
| | 1334 | minchargediff = TMath::Abs(chargedif); |
| | 1335 | } |
| | 1336 | } |
| | 1337 | if (j<0) continue; // not proper cluster |
| | 1338 | Int_t count =0; |
| | 1339 | for (Int_t di=0;di<cnegative[ip];di++){ |
| | 1340 | Int_t jc = negativepair[ip*10+di]; |
| | 1341 | Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ(); |
| | 1342 | if (TMath::Abs(chargedif)<minchargediff+3.) count++; |
| | 1343 | } |
| | 1344 | if (count>1) continue; // more than one "proper" cluster for positive |
| | 1345 | // |
| | 1346 | count =0; |
| | 1347 | for (Int_t dj=0;dj<cpositive[j];dj++){ |
| | 1348 | Int_t ic = positivepair[j*10+dj]; |
| | 1349 | Float_t chargedif = pos[ic].GetQ()-neg[j].GetQ(); |
| | 1350 | if (TMath::Abs(chargedif)<minchargediff+3.) count++; |
| | 1351 | } |
| | 1352 | if (count>1) continue; // more than one "proper" cluster for negative |
| | 1353 | |
| | 1354 | Int_t jp = 0; |
| | 1355 | |
| | 1356 | count =0; |
| | 1357 | for (Int_t dj=0;dj<cnegative[jp];dj++){ |
| | 1358 | Int_t ic = positivepair[jp*10+dj]; |
| | 1359 | Float_t chargedif = pos[ic].GetQ()-neg[jp].GetQ(); |
| | 1360 | if (TMath::Abs(chargedif)<minchargediff+4.) count++; |
| | 1361 | } |
| | 1362 | if (count>1) continue; |
| | 1363 | if (pairs[ip][j]<100) continue; |
| | 1364 | // |
| | 1365 | //almost gold clusters |
| | 1366 | Float_t yp=pos[ip].GetY()*fYpitchSSD; |
| | 1367 | Float_t yn=neg[j].GetY()*fYpitchSSD; |
| | 1368 | Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp); |
| | 1369 | Float_t yt=yn + tann*zt; |
| | 1370 | zt-=fHlSSD; yt-=fHwSSD; |
| | 1371 | ybest=yt; zbest=zt; |
| | 1372 | qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ()); |
| | 1373 | lp[0]=-(-ybest+fYshift[fI]); |
| | 1374 | lp[1]= -zbest+fZshift[fI]; |
| | 1375 | lp[2]=0.0025*0.0025; //SigmaY2 |
| | 1376 | lp[3]=0.110*0.110; //SigmaZ2 |
| | 1377 | lp[4]=qbest; //Q |
| | 1378 | for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2; |
| | 1379 | for (Int_t ilab=0;ilab<3;ilab++){ |
| | 1380 | milab[ilab] = pos[ip].GetLabel(ilab); |
| | 1381 | milab[ilab+3] = neg[j].GetLabel(ilab); |
| | 1382 | } |
| | 1383 | // |
| | 1384 | CheckLabels2(milab); |
| | 1385 | ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ()); |
| | 1386 | milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det |
| | 1387 | Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fI]}; |
| | 1388 | AliITSclusterV2 * cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info); |
| | 1389 | ncl++; |
| | 1390 | cl2->SetChargeRatio(ratio); |
| | 1391 | cl2->SetType(10); |
| | 1392 | pairs[ip][j]=10; |
| | 1393 | if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster |
| | 1394 | cl2->SetType(11); |
| | 1395 | pairs[ip][j]=11; |
| | 1396 | } |
| | 1397 | cused1[ip]++; |
| | 1398 | cused2[j]++; |
| | 1399 | } |
| | 1400 | |
| | 1401 | } |
| | 1402 | |
| | 1403 | // |
| | 1404 | for (Int_t i=0; i<np; i++) { |
| | 1405 | Float_t ybest=1000,zbest=1000,qbest=0; |
| | 1406 | Float_t yp=pos[i].GetY()*fYpitchSSD; |
| | 1407 | if (pos[i].GetQ()<3) continue; |
| | 1408 | for (Int_t j=0; j<nn; j++) { |
| | 1409 | // for (Int_t di = 0;di<cpositive[i];di++){ |
| | 1410 | // Int_t j = negativepair[10*i+di]; |
| | 1411 | if (neg[j].GetQ()<3) continue; |
| | 1412 | if (cused2[j]||cused1[i]) continue; |
| | 1413 | if (pairs[i][j]>0 &&pairs[i][j]<100) continue; |
| | 1414 | ratio = (pos[i].GetQ()-neg[j].GetQ())/(pos[i].GetQ()+neg[j].GetQ()); |
| | 1415 | Float_t yn=neg[j].GetY()*fYpitchSSD; |
| | 1416 | Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp); |
| | 1417 | Float_t yt=yn + tann*zt; |
| | 1418 | zt-=fHlSSD; yt-=fHwSSD; |
| | 1419 | if (TMath::Abs(yt)<fHwSSD+0.01) |
| | 1420 | if (TMath::Abs(zt)<fHlSSD+0.01*(neg[j].GetNd()+pos[i].GetNd())) { |
| | 1421 | ybest=yt; zbest=zt; |
| | 1422 | qbest=0.5*(pos[i].GetQ()+neg[j].GetQ()); |
| | 1423 | lp[0]=-(-ybest+fYshift[fI]); |
| | 1424 | lp[1]= -zbest+fZshift[fI]; |
| | 1425 | lp[2]=0.0025*0.0025; //SigmaY2 |
| | 1426 | lp[3]=0.110*0.110; //SigmaZ2 |
| | 1427 | |
| | 1428 | lp[4]=qbest; //Q |
| | 1429 | for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2; |
| | 1430 | for (Int_t ilab=0;ilab<3;ilab++){ |
| | 1431 | milab[ilab] = pos[i].GetLabel(ilab); |
| | 1432 | milab[ilab+3] = neg[j].GetLabel(ilab); |
| | 1433 | } |
| | 1434 | // |
| | 1435 | CheckLabels2(milab); |
| | 1436 | milab[3]=(((i<<10) + j)<<10) + idet; // pos|neg|det |
| | 1437 | Int_t info[3] = {pos[i].GetNd(),neg[j].GetNd(),fNlayer[fI]}; |
| | 1438 | AliITSclusterV2 * cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info); |
| | 1439 | ncl++; |
| | 1440 | cl2->SetChargeRatio(ratio); |
| | 1441 | cl2->SetType(100+cpositive[j]+cnegative[i]); |
| | 1442 | //cl2->SetType(0); |
| | 1443 | /* |
| | 1444 | if (pairs[i][j]<100){ |
| | 1445 | printf("problem:- %d\n", pairs[i][j]); |
| | 1446 | } |
| | 1447 | if (cnegative[i]<2&&cpositive[j]<2){ |
| | 1448 | printf("problem:- %d\n", pairs[i][j]); |
| | 1449 | } |
| | 1450 | */ |
| | 1451 | } |
| | 1452 | } |
| | 1453 | } |
| | 1454 | |
| | 1455 | // for (Int_t i=0; i<1000; i++) delete [] pairs[i]; |
| | 1456 | // delete [] pairs; |
| | 1457 | |
| | 1458 | } |
| | 1459 | |
| | 1460 | |
| | 1461 | void AliITSclustererV2:: |
| | 1462 | FindClustersSSD(const TClonesArray *alldigits, TClonesArray *clusters) { |
| | 1463 | //------------------------------------------------------------ |
| | 1464 | // Actual SSD cluster finder |
| | 1465 | //------------------------------------------------------------ |
| | 1466 | Int_t smaxall=alldigits->GetEntriesFast(); |
| | 1467 | if (smaxall==0) return; |
| | 1468 | TObjArray *digits = new TObjArray; |
| | 1469 | for (Int_t i=0;i<smaxall; i++){ |
| | 1470 | AliITSdigitSSD *d=(AliITSdigitSSD*)alldigits->UncheckedAt(i); |
| | 1471 | if (d->GetSignal()<3) continue; |
| | 1472 | digits->AddLast(d); |
| | 1473 | } |
| | 1474 | Int_t smax = digits->GetEntriesFast(); |
| | 1475 | if (smax==0) return; |
| | 1476 | |
| | 1477 | const Int_t MAX=1000; |
| | 1478 | Int_t np=0, nn=0; |
| | 1479 | Ali1Dcluster pos[MAX], neg[MAX]; |
| | 1480 | Float_t y=0., q=0., qmax=0.; |
| | 1481 | Int_t lab[4]={-2,-2,-2,-2}; |
| | 1482 | |
| | 1483 | AliITSdigitSSD *d=(AliITSdigitSSD*)digits->UncheckedAt(0); |
| | 1484 | q += d->GetSignal(); |
| | 1485 | y += d->GetCoord2()*d->GetSignal(); |
| | 1486 | qmax=d->GetSignal(); |
| | 1487 | lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2); |
| | 1488 | Int_t curr=d->GetCoord2(); |
| | 1489 | Int_t flag=d->GetCoord1(); |
| | 1490 | Int_t *n=&nn; |
| | 1491 | Ali1Dcluster *c=neg; |
| | 1492 | Int_t nd=1; |
| | 1493 | Int_t milab[10]; |
| | 1494 | for (Int_t ilab=0;ilab<10;ilab++){ |
| | 1495 | milab[ilab]=-2; |
| | 1496 | } |
| | 1497 | milab[0]=d->GetTrack(0); milab[1]=d->GetTrack(1); milab[2]=d->GetTrack(2); |
| | 1498 | |
| | 1499 | for (Int_t s=1; s<smax; s++) { |
| | 1500 | d=(AliITSdigitSSD*)digits->UncheckedAt(s); |
| | 1501 | Int_t strip=d->GetCoord2(); |
| | 1502 | if ((strip-curr) > 1 || flag!=d->GetCoord1()) { |
| | 1503 | c[*n].SetY(y/q); |
| | 1504 | c[*n].SetQ(q); |
| | 1505 | c[*n].SetNd(nd); |
| | 1506 | CheckLabels2(milab); |
| | 1507 | c[*n].SetLabels(milab); |
| | 1508 | //Split suspiciously big cluster |
| | 1509 | /* |
| | 1510 | if (nd>10&&nd<16){ |
| | 1511 | c[*n].SetY(y/q-0.3*nd); |
| | 1512 | c[*n].SetQ(0.5*q); |
| | 1513 | (*n)++; |
| | 1514 | if (*n==MAX) { |
| | 1515 | Error("FindClustersSSD","Too many 1D clusters !"); |
| | 1516 | return; |
| | 1517 | } |
| | 1518 | c[*n].SetY(y/q-0.0*nd); |
| | 1519 | c[*n].SetQ(0.5*q); |
| | 1520 | c[*n].SetNd(nd); |
| | 1521 | (*n)++; |
| | 1522 | if (*n==MAX) { |
| | 1523 | Error("FindClustersSSD","Too many 1D clusters !"); |
| | 1524 | return; |
| | 1525 | } |
| | 1526 | // |
| | 1527 | c[*n].SetY(y/q+0.3*nd); |
| | 1528 | c[*n].SetQ(0.5*q); |
| | 1529 | c[*n].SetNd(nd); |
| | 1530 | c[*n].SetLabels(milab); |
| | 1531 | } |
| | 1532 | else{ |
| | 1533 | */ |
| | 1534 | if (nd>4&&nd<25) { |
| | 1535 | c[*n].SetY(y/q-0.25*nd); |
| | 1536 | c[*n].SetQ(0.5*q); |
| | 1537 | (*n)++; |
| | 1538 | if (*n==MAX) { |
| | 1539 | Error("FindClustersSSD","Too many 1D clusters !"); |
| | 1540 | return; |
| | 1541 | } |
| | 1542 | c[*n].SetY(y/q+0.25*nd); |
| | 1543 | c[*n].SetQ(0.5*q); |
| | 1544 | c[*n].SetNd(nd); |
| | 1545 | c[*n].SetLabels(milab); |
| | 1546 | } |
| | 1547 | (*n)++; |
| | 1548 | if (*n==MAX) { |
| | 1549 | Error("FindClustersSSD","Too many 1D clusters !"); |
| | 1550 | return; |
| | 1551 | } |
| | 1552 | y=q=qmax=0.; |
| | 1553 | nd=0; |
| | 1554 | lab[0]=lab[1]=lab[2]=-2; |
| | 1555 | // |
| | 1556 | for (Int_t ilab=0;ilab<10;ilab++){ |
| | 1557 | milab[ilab]=-2; |
| | 1558 | } |
| | 1559 | // |
| | 1560 | if (flag!=d->GetCoord1()) { n=&np; c=pos; } |
| | 1561 | } |
| | 1562 | flag=d->GetCoord1(); |
| | 1563 | q += d->GetSignal(); |
| | 1564 | y += d->GetCoord2()*d->GetSignal(); |
| | 1565 | nd++; |
| | 1566 | if (d->GetSignal()>qmax) { |
| | 1567 | qmax=d->GetSignal(); |
| | 1568 | lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2); |
| | 1569 | } |
| | 1570 | for (Int_t ilab=0;ilab<10;ilab++) { |
| | 1571 | if (d->GetTrack(ilab)>=0) AddLabel(milab, (d->GetTrack(ilab))); |
| | 1572 | } |
| | 1573 | curr=strip; |
| | 1574 | } |
| | 1575 | c[*n].SetY(y/q); |
| | 1576 | c[*n].SetQ(q); |
| | 1577 | c[*n].SetNd(nd); |
| | 1578 | c[*n].SetLabels(lab); |
| | 1579 | //Split suspiciously big cluster |
| | 1580 | if (nd>4 && nd<25) { |
| | 1581 | c[*n].SetY(y/q-0.25*nd); |
| | 1582 | c[*n].SetQ(0.5*q); |
| | 1583 | (*n)++; |
| | 1584 | if (*n==MAX) { |
| | 1585 | Error("FindClustersSSD","Too many 1D clusters !"); |
| | 1586 | return; |
| | 1587 | } |
| | 1588 | c[*n].SetY(y/q+0.25*nd); |
| | 1589 | c[*n].SetQ(0.5*q); |
| | 1590 | c[*n].SetNd(nd); |
| | 1591 | c[*n].SetLabels(lab); |
| | 1592 | } |
| | 1593 | (*n)++; |
| | 1594 | if (*n==MAX) { |
| | 1595 | Error("FindClustersSSD","Too many 1D clusters !"); |
| | 1596 | return; |
| | 1597 | } |
| | 1598 | |
| | 1599 | FindClustersSSD(neg, nn, pos, np, clusters); |
| | 1600 | } |
| | 1601 | |
| | 1602 | void AliITSclustererV2::FindClustersSSD(AliITSRawStream* input, |
| | 1603 | TClonesArray** clusters) |
| | 1604 | { |
| | 1605 | //------------------------------------------------------------ |
| | 1606 | // Actual SSD cluster finder for raw data |
| | 1607 | //------------------------------------------------------------ |
| | 1608 | Int_t nClustersSSD = 0; |
| | 1609 | const Int_t MAX = 1000; |
| | 1610 | Ali1Dcluster clusters1D[2][MAX]; |
| | 1611 | Int_t nClusters[2] = {0, 0}; |
| | 1612 | Int_t lab[3]={-2,-2,-2}; |
| | 1613 | Float_t q = 0.; |
| | 1614 | Float_t y = 0.; |
| | 1615 | Int_t nDigits = 0; |
| | 1616 | Int_t prevStrip = -1; |
| | 1617 | Int_t prevFlag = -1; |
| | 1618 | Int_t prevModule = -1; |
| | 1619 | |
| | 1620 | // read raw data input stream |
| | 1621 | while (kTRUE) { |
| | 1622 | Bool_t next = input->Next(); |
| | 1623 | |
| | 1624 | if(input->GetSignal()<3 && next) continue; |
| | 1625 | // check if a new cluster starts |
| | 1626 | Int_t strip = input->GetCoord2(); |
| | 1627 | Int_t flag = input->GetCoord1(); |
| | 1628 | if ((!next || (input->GetModuleID() != prevModule)|| |
| | 1629 | (strip-prevStrip > 1) || (flag != prevFlag)) && |
| | 1630 | (nDigits > 0)) { |
| | 1631 | if (nClusters[prevFlag] == MAX) { |
| | 1632 | Error("FindClustersSSD", "Too many 1D clusters !"); |
| | 1633 | return; |
| | 1634 | } |
| | 1635 | Ali1Dcluster& cluster = clusters1D[prevFlag][nClusters[prevFlag]++]; |
| | 1636 | cluster.SetY(y/q); |
| | 1637 | cluster.SetQ(q); |
| | 1638 | cluster.SetNd(nDigits); |
| | 1639 | cluster.SetLabels(lab); |
| | 1640 | |
| | 1641 | //Split suspiciously big cluster |
| | 1642 | if (nDigits > 4&&nDigits < 25) { |
| | 1643 | cluster.SetY(y/q - 0.25*nDigits); |
| | 1644 | cluster.SetQ(0.5*q); |
| | 1645 | if (nClusters[prevFlag] == MAX) { |
| | 1646 | Error("FindClustersSSD", "Too many 1D clusters !"); |
| | 1647 | return; |
| | 1648 | } |
| | 1649 | Ali1Dcluster& cluster2 = clusters1D[prevFlag][nClusters[prevFlag]++]; |
| | 1650 | cluster2.SetY(y/q + 0.25*nDigits); |
| | 1651 | cluster2.SetQ(0.5*q); |
| | 1652 | cluster2.SetNd(nDigits); |
| | 1653 | cluster2.SetLabels(lab); |
| | 1654 | } |
| | 1655 | y = q = 0.; |
| | 1656 | nDigits = 0; |
| | 1657 | } |
| | 1658 | |
| | 1659 | if (!next || (input->GetModuleID() != prevModule)) { |
| | 1660 | Int_t iModule = prevModule; |
| | 1661 | |
| | 1662 | // when all data from a module was read, search for clusters |
| | 1663 | if (prevFlag >= 0) { |
| | 1664 | clusters[iModule] = new TClonesArray("AliITSclusterV2"); |
| | 1665 | fI = iModule; |
| | 1666 | FindClustersSSD(&clusters1D[0][0], nClusters[0], |
| | 1667 | &clusters1D[1][0], nClusters[1], clusters[iModule]); |
| | 1668 | Int_t noClusters = clusters[iModule]->GetEntriesFast(); |
| | 1669 | nClustersSSD += noClusters; |
| | 1670 | } |
| | 1671 | |
| | 1672 | if (!next) break; |
| | 1673 | nClusters[0] = nClusters[1] = 0; |
| | 1674 | y = q = 0.; |
| | 1675 | nDigits = 0; |
| | 1676 | } |
| | 1677 | |
| | 1678 | // add digit to current cluster |
| | 1679 | q += input->GetSignal(); |
| | 1680 | y += strip * input->GetSignal(); |
| | 1681 | nDigits++; |
| | 1682 | prevStrip = strip; |
| | 1683 | prevFlag = flag; |
| | 1684 | prevModule = input->GetModuleID(); |
| | 1685 | |
| | 1686 | } |
| | 1687 | |
| | 1688 | Info("FindClustersSSD", "found clusters in ITS SSD: %d", nClustersSSD); |
| | 1689 | } |
| | 1690 | |
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental