]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/AliITSClusterFinderSSD.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
corrected test for presence of rotation in ProcessLogicalVolume()
[u/mrichter/AliRoot.git] / ITS / AliITSClusterFinderSSD.cxx
CommitLineData
b0f5e3fc 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/*
17
18Adding rekonstruction facilities
19Piotr Krzysztof Skowronski
20December 1999.
766ef0c8 21
22//Piotr Krzysztof Skowronski
23//Warsaw University of Technology
24//skowron@if.pw.edu.pl
25
b0f5e3fc 26*/
27
28/*
2915 -18 V 2000
30Eroor counting routines
31Automatic combination routines improved (traps)
32
33*/
34
3c9a27a7 35#include <iostream.h>
36
37
b0f5e3fc 38#include "AliRun.h"
e8189707 39#include "AliITS.h"
3c9a27a7 40#include "AliITSdigit.h"
41#include "AliITSRawCluster.h"
42#include "AliITSRecPoint.h"
e8189707 43#include "AliITSMapA1.h"
b0f5e3fc 44#include "AliITSClusterFinderSSD.h"
e8189707 45#include "AliITSclusterSSD.h"
46#include "AliITSpackageSSD.h"
3c9a27a7 47#include "AliITSsegmentation.h"
b669392e 48
114da3a4 49const Bool_t AliITSClusterFinderSSD::fgkSIDEP=kTRUE;
50const Bool_t AliITSClusterFinderSSD::fgkSIDEN=kFALSE;
51
3c9a27a7 52const Int_t debug=0;
b0f5e3fc 53
54ClassImp(AliITSClusterFinderSSD)
55
56//____________________________________________________________________
57//
58// Constructor
59//____________________________________________________________________
60//
61
62
63AliITSClusterFinderSSD::AliITSClusterFinderSSD(AliITSsegmentation *seg, TClonesArray *digits, TClonesArray *recp)
64{
766ef0c8 65//Standard constructor
3c9a27a7 66
b0f5e3fc 67 fSegmentation=seg;
68 fDigits=digits;
69 fRecPoints=recp;
70
e8189707 71 fMap = new AliITSMapA1(fSegmentation,fDigits);
72
b0f5e3fc 73 fITS=(AliITS*)gAlice->GetModule("ITS");
74
75 fClusterP = new TClonesArray ("AliITSclusterSSD",200);
76 fNClusterP =0;
77
78 fClusterN= new TClonesArray ("AliITSclusterSSD",200);
79 fNClusterN =0;
80
81 fPackages = new TClonesArray ("AliITSpackageSSD",200); //packages
82 fNPackages = 0;
83
84
85 fDigitsIndexP = new TArrayI(300);
86 fNDigitsP = 0;
87
88 fDigitsIndexN = new TArrayI(300);
89 fNDigitsN = 0;
90
91 SetAlpha1(1000);
92 SetAlpha2(1000);
93 SetAlpha3(1000);
94
95
96 fPitch = fSegmentation->Dpx(0);
97 Float_t StereoP,StereoN;
98 fSegmentation->Angles(StereoP,StereoN);
99 fTanP=TMath::Tan(StereoP);
100 fTanN=TMath::Tan(StereoN);
e8189707 101
be33dccb 102 fPNsignalRatio=7./8.; // warning: hard-wired number
e8189707 103
b0f5e3fc 104}
105
106//-------------------------------------------------------
107AliITSClusterFinderSSD::~AliITSClusterFinderSSD() {
766ef0c8 108// Default destructor
3c9a27a7 109
b0f5e3fc 110 delete fClusterP;
111 delete fClusterN;
112 delete fPackages;
113 delete fDigitsIndexP;
114 delete fDigitsIndexN;
e8189707 115
116 delete fMap;
117
b0f5e3fc 118}
119
120//-------------------------------------------------------
121void AliITSClusterFinderSSD::InitReconstruction()
122{
766ef0c8 123// initialization of the cluster finder
b0f5e3fc 124
125 register Int_t i; //iterator
126
e8189707 127 for (i=0;i<fNClusterP;i++)
b0f5e3fc 128 {
129 fClusterP->RemoveAt(i);
130 }
131 fNClusterP =0;
e8189707 132 for (i=0;i<fNClusterN;i++)
b0f5e3fc 133 {
134 fClusterN->RemoveAt(i);
135 }
136 fNClusterN=0;
137
e8189707 138 for (i=0;i<fNPackages;i++)
b0f5e3fc 139 {
140 fPackages->RemoveAt(i);
141 }
142
143 fNPackages = 0;
144 fNDigitsP=0;
145 fNDigitsN=0;
146
147 Float_t StereoP,StereoN;
148 fSegmentation->Angles(StereoP,StereoN);
149
150 CalcStepFactor (StereoP,StereoN);
151
152 if (debug) cout<<"fSFF = "<<fSFF<<" fSFB = "<<fSFB<<"\n";
153}
154
155
156//---------------------------------------------
157void AliITSClusterFinderSSD::FindRawClusters()
158{
b0f5e3fc 159// This function findes out all clusters belonging to one module
160// 1. Zeroes all space after previous module reconstruction
161// 2. Finds all neighbouring digits
162// 3. If necesery, resolves for each group of neighbouring digits
163// how many clusters creates it.
164// 4. Creates packages
165// 5. Creates clusters
166
167 InitReconstruction(); //ad. 1
e8189707 168 fMap->FillMap();
b0f5e3fc 169 FillDigitsIndex();
170 SortDigits();
171 FindNeighbouringDigits(); //ad. 2
172 SeparateOverlappedClusters(); //ad. 3
173 ClustersToPackages(); //ad. 4
174 ConsumeClusters();
175 PackagesToPoints(); //ad. 5
176 ReconstructNotConsumedClusters();
e8189707 177
178 fMap->ClearMap();
b0f5e3fc 179}
180
181
182//-------------------------------------------------
183void AliITSClusterFinderSSD::FindNeighbouringDigits()
184{
766ef0c8 185//If there are any digits on this side, create 1st Cluster,
186// add to it this digit, and increment number of clusters
b0f5e3fc 187
188
b0f5e3fc 189 register Int_t i;
190
b0f5e3fc 191
192 if ((fNDigitsP>0 ) && (fNDigitsN > 0 )) {
193
194 Int_t currentstripNo;
195 Int_t *dbuffer = new Int_t [300]; //buffer for strip numbers
196 Int_t dnumber; //curent number of digits in buffer
197 TArrayI &lDigitsIndexP = *fDigitsIndexP;
198 TArrayI &lDigitsIndexN = *fDigitsIndexN;
199 TObjArray &lDigits=*fDigits;
200 TClonesArray &lClusterP = *fClusterP;
201 TClonesArray &lClusterN = *fClusterN;
202
203 //process P side
204 dnumber = 1;
205 dbuffer[0]=lDigitsIndexP[0];
206 //If next digit is a neighbour of previous, adds to last cluster this digit
e8189707 207 for (i=1; i<fNDigitsP; i++) {
b0f5e3fc 208 //reads new digit
209 currentstripNo = ((AliITSdigitSSD*)lDigits[lDigitsIndexP[i]])->
210 GetStripNumber();
211 //check if it is a neighbour of a previous one
212 if ( (((AliITSdigitSSD*)lDigits[lDigitsIndexP[i-1]])->GetStripNumber())
213 == (currentstripNo-1) ) dbuffer[dnumber++]=lDigitsIndexP[i];
214 else {
215 //create a new one side cluster
e8189707 216 new(lClusterP[fNClusterP++]) AliITSclusterSSD(dnumber,dbuffer,fDigits,fgkSIDEP);
b0f5e3fc 217 dbuffer[0]=lDigitsIndexP[i];
218 dnumber = 1;
219 }
220 } // end loop over fNDigitsP
e8189707 221 new(lClusterP[fNClusterP++]) AliITSclusterSSD(dnumber,dbuffer,fDigits,fgkSIDEP);
b0f5e3fc 222
223
224 //process N side
225 //for comments, see above
226 dnumber = 1;
227 dbuffer[0]=lDigitsIndexN[0];
228 //If next digit is a neighbour of previous, adds to last cluster this digit
e8189707 229 for (i=1; i<fNDigitsN; i++) {
b0f5e3fc 230 currentstripNo = ((AliITSdigitSSD*)(lDigits[lDigitsIndexN[i]]))->
231 GetStripNumber();
232 if ( (((AliITSdigitSSD*)lDigits[lDigitsIndexN[i-1]])->GetStripNumber())
233 == (currentstripNo-1) ) dbuffer[dnumber++]=lDigitsIndexN[i];
234 else {
e8189707 235 new(lClusterN[fNClusterN++]) AliITSclusterSSD(dnumber,dbuffer,fDigits,fgkSIDEN);
b0f5e3fc 236 dbuffer[0]=lDigitsIndexN[i];
237 dnumber = 1;
238 }
239 } // end loop over fNDigitsN
e8189707 240 new(lClusterN[fNClusterN++]) AliITSclusterSSD(dnumber,dbuffer,fDigits,fgkSIDEN);
b0f5e3fc 241 delete [] dbuffer;
242
243 } // end condition on NDigits
244
245 if (debug) cout<<"\n Found clusters: fNClusterP = "<<fNClusterP<<" fNClusterN ="<<fNClusterN<<"\n";
246
247}
248/**********************************************************************/
249
250
251void AliITSClusterFinderSSD::SeparateOverlappedClusters()
252{
766ef0c8 253// overlapped clusters separation
b0f5e3fc 254
e8189707 255 register Int_t i; //iterator
b0f5e3fc 256
257 Float_t factor=0.75; // How many percent must be lower signal
258 // on the middle one digit
259 // from its neighbours
260 Int_t signal0; //signal on the strip before the current one
261 Int_t signal1; //signal on the current one signal
262 Int_t signal2; //signal on the strip after the current one
263 TArrayI *splitlist; // List of splits
264 Int_t numerofsplits=0; // number of splits
265 Int_t initPsize = fNClusterP; //initial size of the arrays
266 Int_t initNsize = fNClusterN; //we have to keep it because it will grow
267 // in this function and it doasn't make
268 // sense to pass through it again
269
270 splitlist = new TArrayI(300);
271
e8189707 272 for (i=0;i<initPsize;i++)
b0f5e3fc 273 {
274 if (( ((AliITSclusterSSD*)(*fClusterP)[i])->GetNumOfDigits())==1) continue;
275 if (( ((AliITSclusterSSD*)(*fClusterP)[i])->GetNumOfDigits())==2) continue;
276 Int_t nj=(((AliITSclusterSSD*)(*fClusterP)[i])->GetNumOfDigits()-1);
e8189707 277 for (Int_t j=1; j<nj; j++)
b0f5e3fc 278 {
279 signal1=((AliITSclusterSSD*)(*fClusterP)[i])->GetDigitSignal(j);
280 signal0=((AliITSclusterSSD*)(*fClusterP)[i])->GetDigitSignal(j-1);
281 signal2=((AliITSclusterSSD*)(*fClusterP)[i])->GetDigitSignal(j+1);
282 //if signal is less then factor*signal of its neighbours
283 if ( (signal1<(factor*signal0)) && (signal1<(factor*signal2)) )
284 {
285 (*splitlist)[numerofsplits++]=j;
286 }
287 } // end loop over number of digits
288 //split this cluster if necessary
3c9a27a7 289 if(numerofsplits>0) SplitCluster(splitlist,numerofsplits,i,fgkSIDEP);
290 numerofsplits=0;
b0f5e3fc 291
292 //in signed places (splitlist)
293 } // end loop over clusters on Pside
294
e8189707 295 for (i=0;i<initNsize;i++) {
b0f5e3fc 296 if (( ((AliITSclusterSSD*)(*fClusterN)[i])->GetNumOfDigits())==1) continue;
297 if (( ((AliITSclusterSSD*)(*fClusterN)[i])->GetNumOfDigits())==2) continue;
298 Int_t nj=(((AliITSclusterSSD*)(*fClusterN)[i])->GetNumOfDigits()-1);
e8189707 299 for (Int_t j=1; j<nj; j++)
b0f5e3fc 300 {
301 signal1=((AliITSclusterSSD*)(*fClusterN)[i])->GetDigitSignal(j);
302 signal0=((AliITSclusterSSD*)(*fClusterN)[i])->GetDigitSignal(j-1);
303 signal2=((AliITSclusterSSD*)(*fClusterN)[i])->GetDigitSignal(j+1);
304 //if signal is less then factor*signal of its neighbours
305 if ( (signal1<(factor*signal0)) && (signal1<(factor*signal2)) )
306 (*splitlist)[numerofsplits++]=j;
307 } // end loop over number of digits
308 //split this cluster into more clusters
e8189707 309 if(numerofsplits>0) SplitCluster(splitlist,numerofsplits,i,fgkSIDEN);
b0f5e3fc 310 numerofsplits=0; //in signed places (splitlist)
311 } // end loop over clusters on Nside
312
313 delete splitlist;
314}
315
316//-------------------------------------------------------
317void AliITSClusterFinderSSD::SplitCluster(TArrayI *list, Int_t nsplits, Int_t index, Bool_t side)
318{
b0f5e3fc 319 //This function splits one side cluster into more clusters
320 //number of splits is defined by "nsplits"
321 //Place of splits are defined in the TArray "list"
322
323 // For further optimisation: Replace this function by two
324 // specialised ones (each for one side)
325 // save one "if"
326
327 //For comlete comments see AliITSclusterSSD::SplitCluster
328
329
330 register Int_t i; //iterator
331
332 AliITSclusterSSD* curentcluster;
333 Int_t *tmpdigits = new Int_t[100];
334 Int_t NN;
335 // side true means P side
336 if (side) {
337 curentcluster =((AliITSclusterSSD*)((*fClusterP)[index])) ;
e8189707 338 for (i = nsplits; i>0 ;i--) {
b0f5e3fc 339 NN=curentcluster->SplitCluster((*list)[(i-1)],tmpdigits);
340 new ((*fClusterP)[fNClusterP]) AliITSclusterSSD(NN,tmpdigits,fDigits,side);
341 ( (AliITSclusterSSD*)((*fClusterP)[fNClusterP]) )->
342 SetLeftNeighbour(kTRUE);
343 //if left cluster had neighbour on the right before split
344 //new should have it too
345 if ( curentcluster->GetRightNeighbour() )
346 ( (AliITSclusterSSD*)((*fClusterP)[fNClusterP]) )->
347 SetRightNeighbour(kTRUE);
348 else curentcluster->SetRightNeighbour(kTRUE);
349 fNClusterP++;
350 } // end loop over nplits
351 } else {
352 curentcluster =((AliITSclusterSSD*)((*fClusterN)[index]));
e8189707 353 for (i = nsplits; i>0 ;i--) {
b0f5e3fc 354 NN=curentcluster->SplitCluster((*list)[(i-1)],tmpdigits);
355 new ((*fClusterN)[fNClusterN]) AliITSclusterSSD(NN,tmpdigits,fDigits,side);
356 ((AliITSclusterSSD*)((*fClusterN)[fNClusterN]))->
357 SetRightNeighbour(kTRUE);
358 if (curentcluster->GetRightNeighbour())
359 ( (AliITSclusterSSD*)( (*fClusterN)[fNClusterN]) )->
360 SetRightNeighbour(kTRUE);
361 else curentcluster->SetRightNeighbour(kTRUE);
362 fNClusterN++;
363 } // end loop over nplits
364 } // end if side
365 delete []tmpdigits;
366
367}
368
369
370//-------------------------------------------------
371Int_t AliITSClusterFinderSSD::SortDigitsP(Int_t start, Int_t end)
372{
766ef0c8 373// sort digits on the P side
b0f5e3fc 374
b0f5e3fc 375
376 Int_t right;
377 Int_t left;
378 if (start != (end - 1) )
379 {
380 left=this->SortDigitsP(start,(start+end)/2);
381 right=this->SortDigitsP((start+end)/2,end);
382 return (left || right);
383 }
384 else
385 {
386 left = ((AliITSdigitSSD*)((*fDigits)[(*fDigitsIndexP)[start]]))->GetStripNumber();
387 right= ((AliITSdigitSSD*)((*fDigits)[(*fDigitsIndexP)[end]]))->GetStripNumber();
388 if( left > right )
389 {
390 Int_t tmp = (*fDigitsIndexP)[start];
391 (*fDigitsIndexP)[start]=(*fDigitsIndexP)[end];
392 (*fDigitsIndexP)[end]=tmp;
393 return 1;
394 }
395 else return 0;
396 }
397}
398
399
400//--------------------------------------------------
401
402Int_t AliITSClusterFinderSSD::SortDigitsN(Int_t start, Int_t end)
403{
766ef0c8 404// sort digits on the N side
b0f5e3fc 405
406 Int_t right;
407 Int_t left;
408 if (start != (end - 1))
409 {
410 left=this->SortDigitsN(start,(start+end)/2);
411 right=this->SortDigitsN((start+end)/2,end);
412 return (left || right);
413 }
414 else
415 {
416 left =((AliITSdigitSSD*)((*fDigits)[(*fDigitsIndexN)[start]]))->GetStripNumber();
417 right=((AliITSdigitSSD*)((*fDigits)[(*fDigitsIndexN)[end]]))->GetStripNumber();
418 if ( left > right )
419 {
420 Int_t tmp = (*fDigitsIndexN)[start];
421 (*fDigitsIndexN)[start]=(*fDigitsIndexN)[end];
422 (*fDigitsIndexN)[end]=tmp;
423 return 1;
424 }else return 0;
425 }
426}
427
428
429//------------------------------------------------
430void AliITSClusterFinderSSD::FillDigitsIndex()
431{
b0f5e3fc 432 //Fill the indexes of the clusters belonging to a given ITS module
b0f5e3fc 433
434 Int_t PNs=0, NNs=0;
435 Int_t tmp,bit,k;
436 Int_t N;
437 Int_t i;
438
439 N = fDigits->GetEntriesFast();
440
441 Int_t* PSidx = new Int_t [N*sizeof(Int_t)];
442 Int_t* NSidx = new Int_t [N*sizeof(Int_t)];
443 if (fDigitsIndexP==NULL) fDigitsIndexP = new TArrayI(N);
444 if (fDigitsIndexN==NULL) fDigitsIndexN = new TArrayI(N);
445
446 AliITSdigitSSD *dig;
447
e8189707 448 for ( i = 0 ; i< N; i++ ) {
b0f5e3fc 449 dig=(AliITSdigitSSD*)fDigits->UncheckedAt(i);
450 if(dig->IsSideP()) {
451 bit=1;
452 tmp=dig->GetStripNumber();
453 // I find this totally unnecessary - it's just a
454 // CPU consuming double check
e8189707 455 for( k=0;k<PNs;k++)
b0f5e3fc 456 {
457 if (tmp==PSidx[k])
458 {
459 if (debug) cout<<"Such a digit exists \n";
460 bit=0;
461 }
462 }
463 // end comment
464 if(bit) {
465 fDigitsIndexP->AddAt(i,fNDigitsP++);
466 PSidx[PNs++]=tmp;
467 }
468 } else {
469 bit=1;
470 tmp=dig->GetStripNumber();
471 // same as above
e8189707 472 for( k=0;k<NNs;k++)
b0f5e3fc 473 {
474 if (tmp==NSidx[k])
475 {
476 if (debug) cout<<"Such a digit exists \n";
477 bit=0;
478 }
479 }
480 // end comment
481 if (bit) {
482 fDigitsIndexN->AddAt(i,fNDigitsN++);
483 NSidx[NNs++] =tmp;
484 }
485 }
486 }
487
3c9a27a7 488 delete [] PSidx;
489 delete [] NSidx;
b0f5e3fc 490
491 if (debug) cout<<"Digits : P = "<<fNDigitsP<<" N = "<<fNDigitsN<<endl;
492
493}
494
495
496//-------------------------------------------
497
498void AliITSClusterFinderSSD::SortDigits()
499{
766ef0c8 500// sort digits
b0f5e3fc 501
502 Int_t i;
503 if(fNDigitsP>1)
e8189707 504 for (i=0;i<fNDigitsP-1;i++)
b0f5e3fc 505 if (SortDigitsP(0,(fNDigitsP-1-i))==0) break;
506
507 if(fNDigitsN>1)
e8189707 508 for (i=0;i<fNDigitsN-1;i++)
b0f5e3fc 509 if(SortDigitsN(0,(fNDigitsN-1-i))==0) break;
510}
511
512
513
514//----------------------------------------------
515void AliITSClusterFinderSSD::FillClIndexArrays(Int_t* arrayP, Int_t *arrayN)
516{
766ef0c8 517// fill cluster index array
b0f5e3fc 518
519 register Int_t i;
e8189707 520 for (i=0; i<fNClusterP;i++)
b0f5e3fc 521 {
522 arrayP[i]=i;
523 }
e8189707 524 for (i=0; i<fNClusterN;i++)
b0f5e3fc 525 {
526 arrayN[i]=i;
527 }
528}
529
530
531//------------------------------------------------------
532void AliITSClusterFinderSSD::SortClusters(Int_t* arrayP, Int_t *arrayN)
533{
766ef0c8 534// sort clusters
b0f5e3fc 535
536 Int_t i;
537 if(fNClusterP>1)
e8189707 538 for (i=0;i<fNClusterP-1;i++)
b0f5e3fc 539 if (SortClustersP(0,(fNClusterP-1),arrayP)==0) break;
540
541
542 if(fNClusterN>1)
e8189707 543 for (i=0;i<fNClusterN-1;i++)
b0f5e3fc 544 if (SortClustersN(0,(fNClusterN-1),arrayN)==0) break;
545
546}
547
548
549
550//---------------------------------------------------
551Int_t AliITSClusterFinderSSD::SortClustersP(Int_t start, Int_t end, Int_t *array)
552{
766ef0c8 553//Sort P side clusters
b0f5e3fc 554
555
556 Int_t right;
557 Int_t left;
558 if (start != (end - 1) ) {
559 left=this->SortClustersP(start,(start+end)/2,array);
560 right=this->SortClustersP((start+end)/2,end,array);
561 return (left || right);
562 } else {
563 left =((AliITSclusterSSD*)((*fClusterP)[array[start]]))->
564 GetDigitStripNo(0);
565 right=((AliITSclusterSSD*)((*fClusterP)[array[ end ]]))->
566 GetDigitStripNo(0);
567 if(left>right) {
568 Int_t tmp = array[start];
569 array[start]=array[end];
570 array[end]=tmp;
571 return 1;
572 } else return 0;
573 }
574
575
576}
577
578
579
580//-------------------------------------------------------
581Int_t AliITSClusterFinderSSD::SortClustersN(Int_t start, Int_t end, Int_t *array)
582{
766ef0c8 583//Sort N side clusters
b0f5e3fc 584
585
b0f5e3fc 586 Int_t right;
587 Int_t left;
588
589 if (start != (end - 1) ) {
590 left=this->SortClustersN(start,(start+end)/2,array);
591 right=this->SortClustersN((start+end)/2,end,array);
592 return (left || right);
593 } else {
594 left =((AliITSclusterSSD*)((*fClusterN)[array[start]]))->
595 GetDigitStripNo(0);
596 right=((AliITSclusterSSD*)((*fClusterN)[array[ end ]]))->
597 GetDigitStripNo(0);
598 if( left > right) {
599 Int_t tmp = array[start];
600 array[start]=array[end];
601 array[end]=tmp;
602 return 1;
603 } else return 0;
604 }
605
606}
607
608
609
610//-------------------------------------------------------
611void AliITSClusterFinderSSD::ClustersToPackages()
612{
766ef0c8 613// fill packages
b0f5e3fc 614
615 Int_t *oneSclP = new Int_t[fNClusterP]; //I want to have sorted 1 S clusters
616 Int_t *oneSclN = new Int_t[fNClusterN]; //I can not sort it in TClonesArray
617 //so, I create table of indexes and
618 //sort it
619 //I do not use TArrayI on purpose
620 // MB: well, that's not true that one
621 //cannot sort objects in TClonesArray
622
e8189707 623 register Int_t i; //iterator
b0f5e3fc 624 AliITSpackageSSD *currentpkg;
625 AliITSclusterSSD *currentP;
626 AliITSclusterSSD *currentN;
627 AliITSclusterSSD *fcurrN;
628
629 Int_t lastNclIndex=0; //index of last N sidecluster
630 Int_t Nclpack=0; //number of P clusters in current package
631 Int_t Pclpack=0; //number of N clusters in current package
632 Int_t firstStripP;
633 Int_t lastStripP;
634 Int_t firstStripN;
635 Int_t lastStripN;
636 Int_t tmplastNclIndex;
637
638//Fills in One Side Clusters Index Array
639 FillClIndexArrays(oneSclP,oneSclN);
640//Sorts filled Arrays
641 SortClusters(oneSclP,oneSclN);
642
643
644 fNPackages=1;
645 new ((*fPackages)[0]) AliITSpackageSSD(fClusterP,fClusterN);
646 currentpkg = (AliITSpackageSSD*)((*fPackages)[0]);
647
648 //main loop on sorted P side clusters
e8189707 649 Int_t k;
650 for (i=0;i<fNClusterP;i++) {
b0f5e3fc 651 //Take new P side cluster
652 currentP = GetPSideCluster(oneSclP[i]);
653 currentN = GetNSideCluster(oneSclN[lastNclIndex]);
654 // take a new P cluster or not ?
655 if(IsCrossing(currentP,currentN)) {
656 // don't take a new P cluster
657 Pclpack++;
658 currentP->AddCross(oneSclN[lastNclIndex]);
659 currentN->AddCross(oneSclP[i]);
660 currentpkg->AddPSideCluster(oneSclP[i]);
661 if (Nclpack==0) {
662 currentpkg->AddNSideCluster(oneSclN[lastNclIndex]);
663 Nclpack++;
664 }
665 //check if previous N side clusters crosses with it too
666 for(k=1;k<fSFB+1;k++) {
667 //check if we are not out of array
668 if ((lastNclIndex-k)>-1) {
669 fcurrN = GetNSideCluster( oneSclN[lastNclIndex-k] );
670 if( IsCrossing(currentP,fcurrN) ) {
671 currentP->AddCross(oneSclN[lastNclIndex-k]);
672 fcurrN->AddCross(oneSclP[i]);
673 } else break; //There is no sense to check rest of clusters
674 } else break;
675 }
676 tmplastNclIndex=lastNclIndex;
677 //Check if next N side clusters crosses with it too
e8189707 678 for (Int_t k=1;k<fSFF+1;k++) {
b0f5e3fc 679 //check if we are not out of array
680 if ((tmplastNclIndex+k)<fNClusterN) {
681 fcurrN = GetNSideCluster( oneSclN[tmplastNclIndex+k] );
682 if(IsCrossing(currentP,fcurrN) ) {
683 lastNclIndex++;
684 fcurrN->AddCross(oneSclP[i]);
685 currentP->AddCross(oneSclN[tmplastNclIndex+k]);
686 currentpkg->AddNSideCluster(oneSclN[lastNclIndex]);
687 Nclpack++;
688 }else break;
689 } else break;
690 }
691
692 // end of package
693 //if( IsCrossing )
694 } else {
695 lastStripP = currentP->GetLastDigitStripNo();
696 lastStripN = currentN->GetLastDigitStripNo();
697
698 if((lastStripN-fSFF) < lastStripP ) {
699 //Take new PCluster
700 if((Nclpack>0)&& (Pclpack>0)) {
701 new ((*fPackages)[fNPackages]) AliITSpackageSSD(fClusterP,fClusterN);
702 currentpkg = (AliITSpackageSSD*)((*fPackages)[fNPackages++]);
703 }
704
705 //so we have to take another cluster on side N and check it
706 //we have to check until taken cluster's last strip will
707 //be > last of this cluster(P)
708 //so it means that there is no corresponding cluster on side N
709 //to this cluster (P)
710 //so we have to continue main loop with new "lastNclIndex"
711
712 Nclpack=0;
713 Pclpack=0;
714 //We are not sure that next N cluter will cross with this P cluster
715 //There might one, or more, clusters that do not have any
716 //corresponding cluster on the other side
717 for(;;) {
718 lastNclIndex++;
719 //Check if we are not out of array
720 if (lastNclIndex<fNClusterN) {
721 currentN = GetNSideCluster(oneSclN[lastNclIndex]);
722 if ( IsCrossing(currentP, currentN) ){
723 Nclpack++;
724 Pclpack++;
725 currentP->AddCross(oneSclN[lastNclIndex]);
726 currentN->AddCross(oneSclP[i]);
727 currentpkg->AddNSideCluster(oneSclN[lastNclIndex]);
728 currentpkg->AddPSideCluster(oneSclP[i]);
729 //Check, if next N side clusters crosses with it too
730 tmplastNclIndex=lastNclIndex;
e8189707 731 for (Int_t k=1;k<fSFF+1;k++) {
b0f5e3fc 732 //Check if we are not out of array
733 if ((tmplastNclIndex+k)<fNClusterN) {
734 fcurrN = GetNSideCluster(oneSclN[tmplastNclIndex+k]);
735 if( IsCrossing(currentP, fcurrN) ) {
736 Nclpack++;
737 lastNclIndex++;
738 currentP->AddCross(oneSclN[tmplastNclIndex+k]);
739 fcurrN->AddCross(oneSclP[i]);
740 currentpkg->AddNSideCluster(oneSclN[tmplastNclIndex+k]);
741 }else break;
742 } else break; //we are out of array
743 } // end loop
744 break;
745 } else {
746 firstStripP = currentP->GetFirstDigitStripNo();
747 firstStripN = currentN->GetFirstDigitStripNo();
748 if((firstStripN+fSFB)>=firstStripP) break;
749 else continue;
750 }
751 } else goto EndOfFunction;
752 } // end for(;;)
753 } else //EndOfPackage
754 {
755 continue;
756 } //else EndOfPackage
757
758 }//else IsCrossing
759
760 }//main for
761
762 EndOfFunction:
763 if ((Nclpack<1) || (Pclpack<1)) fNPackages--;
764
765 delete oneSclP;
766 delete oneSclN;
767
768}
769
770
771
772//-----------------------------------------------
773void AliITSClusterFinderSSD::PackagesToPoints()
774{
766ef0c8 775// find recpoints from packages
b0f5e3fc 776
777 register Int_t i;
778 AliITSpackageSSD *currentpkg;
779 Int_t NumNcl;
780 Int_t NumPcl;
781 Int_t clusterIndex;
782 Bool_t clusterSide;
783
e8189707 784 for (i=0;i<fNPackages;i++) {
b0f5e3fc 785 //get pointer to next package
786 currentpkg = (AliITSpackageSSD*)((*fPackages)[i]);
787 NumNcl = currentpkg->GetNumOfClustersN();
788 NumPcl = currentpkg->GetNumOfClustersP();
789 if(debug) cout<<"New Package\nNumPcl ="<<NumPcl<<" NumNcl ="<<NumNcl<<"\n";
790
791 while(((NumPcl>=2)&&(NumNcl>2))||((NumPcl>2)&&(NumNcl>=2))) {
792 //This is case of "big" pacakage
793 //conditions see 3 lines above
794 //if in the big package exists cluster with one cross
795 //we can reconstruct this point without any geometrical ambiguities
796 if(currentpkg->GetClusterWithOneCross(clusterIndex, clusterSide) )
797 {
798 ResolveClusterWithOneCross(currentpkg, clusterIndex, clusterSide);
799 } else if (clusterIndex == -2) {
800 NumPcl = 0;
801 NumNcl = 0;
802 break;
803 } else {
804 if ( (NumNcl==NumPcl) && (NumPcl<10)) {
805 //if there is no cluster with one cross
806 //we can resolve whole package at once
807 //by finding best combination
808 //we can make combinations when NumNcl==NumPcl
809 if (ResolvePackageBestCombin(currentpkg)) {
810 //sometimes creating combinations fail,
811 //but it happens very rarely
812 NumPcl = 0;
813 NumNcl = 0;
814 break;
815 } else ResolveOneBestMatchingPoint(currentpkg);
816 } else {
817 ResolveOneBestMatchingPoint(currentpkg);
818 }
819 }
820 NumNcl = currentpkg->GetNumOfClustersN();
821 NumPcl = currentpkg->GetNumOfClustersP();
822 } // end while
823 if ((NumPcl==1)&&(NumNcl==1)) {
824 ResolveSimplePackage(currentpkg);
825 continue;
826 }
827 if (NumPcl==1) {
828 ResolvePackageWithOnePSideCluster(currentpkg);
829 continue;
830 }
831 if (NumNcl==1) {
832 ResolvePackageWithOneNSideCluster(currentpkg);
833 continue;
834 }
835 if ((NumPcl==2)&&(NumNcl==2)) {
836 ResolveTwoForTwoPackage(currentpkg);
837 continue;
838 }
839 if ((NumPcl==0)&&(NumNcl>0)) { }
840 if ((NumNcl==0)&&(NumPcl>0)) { }
841
842 } // end loop over fNPackages
843
844
845}
846
847
848//----------------------------------------------------------
849
850void AliITSClusterFinderSSD::
851ResolveClusterWithOneCross(AliITSpackageSSD *currentpkg, Int_t clusterIndex, Bool_t clSide)
852{
766ef0c8 853// resolve clusters with one cross
854/*
855
856ie:
857 \ /
858 \ /
859 \/
860 /\
861 / \
862 / \
863 / \
864 / \
865 / \
866*/
867
868
b0f5e3fc 869
e8189707 870 if (clSide == fgkSIDEP) ResolvePClusterWithOneCross(currentpkg,clusterIndex);
b0f5e3fc 871 else ResolveNClusterWithOneCross(currentpkg,clusterIndex);
872
873}
874
875
876//---------------------------------------------------------
877void AliITSClusterFinderSSD::
878ResolvePClusterWithOneCross(AliITSpackageSSD *pkg, Int_t clusterIndex)
879{
766ef0c8 880// resolve clusters with one cross on P side
b0f5e3fc 881/*
882There is cluster (side P) which crosses with only one cluster on the other
883side (N)
884
885ie:
886 \ / \/ /
887 \ / /\ /
888 \/ / \/
889 /\ / /\
890 / \/ / \ .....
891 / /\ / \
892 / / \/ \
893 / / /\ \
894 / / / \ \
895*/
896
897
898 AliITSclusterSSD * clusterP;
899 AliITSclusterSSD * clusterN;
900
901 Float_t posClusterP; //Cluster P position in strip coordinates
902 Float_t posClusterN; //Cluster N position in strip coordinates
903
904 Float_t posErrorClusterP;
905 Float_t posErrorClusterN;
906
907 Float_t sigClusterP;
908 Float_t sigClusterN;
909
910 Float_t sigErrorClusterP;
911 Float_t sigErrorClusterN;
912
913 Float_t ps;
914 Float_t ns;
915
916 Float_t Chicomb;
917 Int_t clusterIdx;
918
919 if (debug) cout<<"ResolvePClusterWithOneCross\n";
920
921 clusterP=pkg->GetPSideCluster(clusterIndex);
922 posClusterP = GetClusterZ(clusterP);
923 posErrorClusterP = clusterP->GetPositionError();
924 sigClusterP = ps = clusterP->GetTotalSignal();
925 sigErrorClusterP = clusterP->GetTotalSignalError();
926 //carefully, it returns index in TClonesArray
927
928 clusterIdx = clusterP->GetCross(0);
929 clusterN=this->GetNSideCluster(clusterIdx);
930 ns = clusterN->GetTotalSignal();
931 posClusterN = GetClusterZ(clusterN);
932 posErrorClusterN = clusterN->GetPositionError();
e8189707 933 pkg->DelCluster(clusterIndex,fgkSIDEP);
be33dccb 934 sigClusterN = ps/fPNsignalRatio;
b0f5e3fc 935 // there is no sonse to check how signal ratio is far from perfect
936 // matching line if the if below it is true
937 if (ns < sigClusterN) {
938 sigClusterN=ns;
be33dccb 939 if (debug) cout<<"n1 < p1/fPNsignalRatio";
b0f5e3fc 940 if (debug) cout<<"Attempting to del cluster N "<<clusterIdx<<" ... \n";
e8189707 941 pkg->DelClusterOI(clusterIdx,fgkSIDEN);
b0f5e3fc 942 } else {
943 //Let's see how signal ratio is far from perfect matching line
944 Chicomb = DistToPML(ps,ns);
945 if (debug) cout<<"Chic "<<Chicomb<<"\n";
766ef0c8 946 if (Chicomb > fAlpha2) {
b0f5e3fc 947 //it is near, so we can risk throwing this cluster away too
948 if (debug) cout<<"Attempting to del cluster N "<<clusterIdx<<"...\n";
e8189707 949 pkg->DelClusterOI(clusterIdx,fgkSIDEN);
b0f5e3fc 950 } else {
951 clusterN->CutTotalSignal(sigClusterN);
952 if (debug) cout <<"Signal cut |||||||||||||\n";
953 }
954 }
955 sigErrorClusterN= clusterN->GetTotalSignalError();
956 CreateNewRecPoint(posClusterP,posErrorClusterP,posClusterN,posErrorClusterN,
957 sigClusterP+sigClusterN, sigErrorClusterN+sigErrorClusterP,
958 clusterP, clusterN, 0.75);
959
960}
961
962
963//------------------------------------------------------
964void AliITSClusterFinderSSD::
965ResolveNClusterWithOneCross(AliITSpackageSSD *pkg, Int_t clusterIndex)
966{
766ef0c8 967// resolve clusters with one cross on N side
968
b0f5e3fc 969
970 AliITSclusterSSD * clusterP;
971 AliITSclusterSSD * clusterN;
972
973 Float_t posClusterP; //Cluster P position in strip coordinates
974 Float_t posClusterN; //Cluster N position in strip coordinates
975
976 Float_t posErrorClusterP;
977 Float_t posErrorClusterN;
978
979 Float_t sigClusterP;
980 Float_t sigClusterN;
981
982 Float_t sigErrorClusterP;
983 Float_t sigErrorClusterN;
984
985 Float_t ps;
986 Float_t ns;
987
988 Float_t Chicomb;
989 Int_t clusterIdx;
990
991 if (debug) cout<<"ResolveNClusterWithOneCross\n";
992
993 clusterN=pkg->GetNSideCluster(clusterIndex);
994 posClusterN = GetClusterZ(clusterN);
995 posErrorClusterN = clusterN->GetPositionError();
996 sigClusterN = ns = clusterN->GetTotalSignal();
997 sigErrorClusterN = clusterN->GetTotalSignalError();
998 //carefully, it returns index in TClonesArray
999 clusterIdx = clusterN->GetCross(0);
1000 clusterP=this->GetPSideCluster(clusterIdx);
1001 ps = clusterP->GetTotalSignal();
1002 posClusterP = GetClusterZ(clusterP);
1003 posErrorClusterP = clusterP->GetPositionError();
e8189707 1004 pkg->DelCluster(clusterIndex,fgkSIDEN);
be33dccb 1005 sigClusterP=ns*fPNsignalRatio;
b0f5e3fc 1006 // there is no sonse to check how signal ratio is far from perfect
1007 // matching line if the if below it is true
1008 if (ps < sigClusterP) {
1009 sigClusterP = ps;
be33dccb 1010 if (debug) cout<<"ps < ns*fPNsignalRatio";
b0f5e3fc 1011 if (debug) cout<<"Attempting to del cluster P "<<clusterIdx<<" ... \n";
e8189707 1012 pkg->DelClusterOI(clusterIdx,fgkSIDEP);
b0f5e3fc 1013 } else {
1014 //Let's see how signal ratio is far from perfect matching line
1015 Chicomb = DistToPML(ps,ns);
1016 if (debug) cout<<"Chic "<<Chicomb<<"\n";
766ef0c8 1017 if (Chicomb > fAlpha2) {
b0f5e3fc 1018 //it is near, so we can risk frowing this cluster away too
1019 if (debug) cout<<"Attempting to del cluster P "<<clusterIdx<<"...\n";
e8189707 1020 pkg->DelClusterOI(clusterIdx,fgkSIDEP);
b0f5e3fc 1021 } else {
1022 clusterN->CutTotalSignal(sigClusterP);
1023 if (debug) cout <<"Signal cut ------------\n";
1024 }
1025 }
1026 sigErrorClusterP= clusterP->GetTotalSignalError();
1027 CreateNewRecPoint( posClusterP,posErrorClusterP,posClusterN,posErrorClusterN,
1028 sigClusterP+sigClusterN,sigErrorClusterN+sigErrorClusterP,
1029 clusterP, clusterN, 0.75);
1030
1031
1032}
1033
1034
1035
1036//-------------------------------------------------------
1037
1038Bool_t AliITSClusterFinderSSD::
1039ResolvePackageBestCombin(AliITSpackageSSD *pkg)
1040{
766ef0c8 1041//find best combination
1042
b0f5e3fc 1043 if (debug) cout<<"NumNcl==NumPcl ("<<pkg->GetNumOfClustersN()
1044 <<"=="<<pkg->GetNumOfClustersP()<<"); Generating combinations ... \n";
1045
1046
1047 AliITSclusterSSD * clusterP;
1048 AliITSclusterSSD * clusterN;
1049
1050 Int_t Ncombin; //Number of combinations
1051 Int_t itera; //iterator
1052 Int_t sizet=1; //size of array to allocate
e8189707 1053
b0f5e3fc 1054 Int_t NP = pkg->GetNumOfClustersP();
e8189707 1055 for (itera =2; itera <= NP ;itera ++) {
b0f5e3fc 1056 sizet=sizet*itera;
1057 if (sizet > 10000) {
1058 sizet=10000;
1059 break;
1060 }
1061 }
1062
e8189707 1063 Int_t** combin = new Int_t*[sizet]; //2D array to keep combinations in
b0f5e3fc 1064
e8189707 1065 for (itera =0; itera <sizet;itera++) {
b0f5e3fc 1066 combin[itera] = new Int_t[NP+1];
1067 }
1068
1069 pkg->GetAllCombinations(combin,Ncombin,sizet);
1070
1071 if (Ncombin==0) {
1072 if (debug) cout<<"No combin Error";
1073 return kFALSE;
1074 }
1075//calculate which combination fits the best to perfect matching line
1076 Int_t bc = GetBestComb(combin,Ncombin,NP,pkg);
1077 if (debug) cout<<" bc "<<bc <<"\n";
1078
e8189707 1079 for (itera =0; itera < NP; itera ++) {
b0f5e3fc 1080 clusterP = pkg->GetPSideCluster(itera);
1081 //carefully here
1082 //becase AliITSclusterSSD::GetCross returns index in
1083 //AliITSmoduleSSD.fClusterP, which is put to "combin"
1084 clusterN = GetNSideCluster(combin[bc][itera]);
1085 CreateNewRecPoint(clusterP, clusterN, 0.75);
1086 }
1087
e8189707 1088 for (itera =0; itera <sizet;itera++) {
b0f5e3fc 1089 delete [](combin[itera]);
1090 }
1091 delete [] combin;
1092 return kTRUE;
1093
1094}
1095
1096
1097//----------------------------------------------------
1098void AliITSClusterFinderSSD::
1099ResolveOneBestMatchingPoint(AliITSpackageSSD *pkg)
1100{
766ef0c8 1101// find best matching point
b0f5e3fc 1102
1103 Int_t ni, pi;
1104
1105 Int_t prvP, prvN;
1106 Int_t nextP, nextN=0;
1107
1108
1109 AliITSclusterSSD * clusterP;
1110 AliITSclusterSSD * clusterN;
1111
1112 Bool_t split = kFALSE;
1113
1114 if (debug) cout<<"ResolveOneBestMatchingPoint \n";
1115
1116 GetBestMatchingPoint(pi, ni, pkg);
1117
1118 clusterP = GetPSideCluster(pi);
1119 clusterN = GetNSideCluster(ni);
1120
1121 CreateNewRecPoint(clusterP, clusterN, 0.75);
1122
1123 if ((nextP=pkg->GetNextPIdx(pi))!=-1)
1124 if ((nextN=pkg->GetNextNIdx(ni))!=-1)
1125 if ((prvP=pkg->GetPrvPIdx(pi))!=-1)
1126 if ((prvN=pkg->GetPrvNIdx(ni))!=-1)
1127 if( !(GetPSideCluster(prvP)->IsCrossingWith(nextN)))
1128 if( !(GetPSideCluster(nextP)->IsCrossingWith(prvN)))
1129 {
1130 split = kTRUE;
1131 }
1132
1133
e8189707 1134 pkg->DelClusterOI(pi, fgkSIDEP);
1135 pkg->DelClusterOI(ni, fgkSIDEN);
3c9a27a7 1136
b0f5e3fc 1137 if (split) {
1138 if (debug) cout<<"spltting package ...\n";
1139 new ((*fPackages)[fNPackages]) AliITSpackageSSD(fClusterP,fClusterN);
1140 pkg->
1141 SplitPackage(nextP,nextN,(AliITSpackageSSD*)((*fPackages)[fNPackages++]));
1142 }
1143
1144}
1145
1146
1147//--------------------------------------------------
1148void AliITSClusterFinderSSD::ResolveSimplePackage(AliITSpackageSSD *pkg)
1149{
766ef0c8 1150// resolve simple package
b0f5e3fc 1151
1152 AliITSclusterSSD * clusterP;
1153 AliITSclusterSSD * clusterN;
1154
1155 clusterP = pkg->GetPSideCluster(0);
1156
1157 clusterN = pkg->GetNSideCluster(0);
1158
1159 CreateNewRecPoint(clusterP, clusterN, 1.0);
1160
1161
1162}
1163
1164
1165//--------------------------------------------------
1166void AliITSClusterFinderSSD:: ResolvePackageWithOnePSideCluster(AliITSpackageSSD *pkg) {
766ef0c8 1167// resolve P side clusters from packages
b0f5e3fc 1168
1169/*
1170 \ \ \ /
1171 \ \ \/
1172 \ \ /\
1173 \ \/ \
1174 \ /\ \
1175 \/ \ \
1176 /\ \ \
1177 / \ \ \
1178
1179*/
1180
1181
1182/************************/
1183/**************************/
1184/*XP SP itg jest tylko jeden nie musi byc tablica i w peetli nie trzeba po niej iterowcac*/
1185/***************************/
1186
1187
1188 Int_t k;
1189 AliITSclusterSSD * clusterP;
1190 AliITSclusterSSD * clusterN;
1191
1192 Int_t NN=pkg->GetNumOfClustersN();
1193 Float_t sumsig = 0;
1194
1195 Float_t XP;
1196 Float_t XPerr;
1197
1198 Float_t * XN = new Float_t[NN];
1199 Float_t * XNerr = new Float_t[NN];
1200
1201 Float_t * SP = new Float_t[NN];
1202 Float_t * SN = new Float_t[NN];
1203
1204 Float_t * SPerr = new Float_t[NN];
1205 Float_t * SNerr = new Float_t[NN];
1206
1207 Float_t p1;
1208 Float_t p1err;
1209
1210 clusterP = pkg->GetPSideCluster(0);
1211 p1 = clusterP->GetTotalSignal();
1212
1213 XP = GetClusterZ(clusterP);
1214 XPerr = clusterP->GetPositionError();
1215 p1err = clusterP->GetTotalSignalError();
1216
e8189707 1217 for (k=0;k<NN;k++) {
b0f5e3fc 1218 SN[k] = pkg->GetNSideCluster(k)->GetTotalSignal();
1219 SNerr[k] = pkg->GetNSideCluster(k)->GetTotalSignalError();
1220 sumsig += SN[k];
1221 }
e8189707 1222 for (k=0;k<NN;k++) {
b0f5e3fc 1223 clusterN = pkg->GetNSideCluster(k);
1224 SP[k]= p1*SN[k]/sumsig;
1225 SPerr[k] = p1err*SN[k]/sumsig;
1226 XN[k]=GetClusterZ(clusterN);
1227 XNerr[k]= clusterN->GetPositionError();
1228 CreateNewRecPoint(XP,XPerr, XN[k],XNerr[k], SP[k]+SN[k], SPerr[k]+SNerr[k], clusterP, clusterN, 1.0);
1229
1230 }
1231
3c9a27a7 1232 delete [] XN;
1233 delete [] XNerr;
1234 delete [] SP;
1235 delete [] SN;
1236 delete [] SPerr;
1237 delete [] SNerr;
1238
1239
b0f5e3fc 1240}
1241
1242
1243//---------------------------------------------------------
1244void AliITSClusterFinderSSD::ResolvePackageWithOneNSideCluster(AliITSpackageSSD *pkg) {
766ef0c8 1245// resolve N side clusters from packages
b0f5e3fc 1246
1247/*
1248 \ / / /
1249 \ / / /
1250 \/ / /
1251 /\ / /
1252 / \/ /
1253 / /\ /
1254 / / \/
1255 / / /\
1256 / / / \
1257*/
1258
1259
1260 Int_t k;
1261 AliITSclusterSSD * clusterP;
1262 AliITSclusterSSD * clusterN;
1263
1264 Int_t NP=pkg->GetNumOfClustersP();
1265 Float_t sumsig = 0;
1266
1267 Float_t * XP = new Float_t[NP];
1268 Float_t * XPerr = new Float_t[NP];
1269
1270 Float_t XN;
1271 Float_t XNerr;
1272
1273 Float_t * SP = new Float_t[NP];
1274 Float_t * SN = new Float_t[NP];
1275
1276 Float_t * SPerr = new Float_t[NP];
1277 Float_t * SNerr = new Float_t[NP];
1278
1279 Float_t n1;
1280 Float_t n1err;
1281
1282 clusterN = pkg->GetNSideCluster(0);
1283 n1 = clusterN->GetTotalSignal();
1284
1285 XN = GetClusterZ(clusterN);
1286 XNerr = clusterN->GetPositionError();
1287
1288 n1err=clusterN->GetTotalSignalError();
1289
e8189707 1290 for (k=0;k<NP;k++) {
b0f5e3fc 1291 SP[k] = pkg->GetPSideCluster(k)->GetTotalSignal();
1292 sumsig += SP[k];
1293 SPerr[k] = pkg->GetPSideCluster(k)->GetTotalSignalError();
1294 }
1295
e8189707 1296 for (k=0;k<NP;k++) {
b0f5e3fc 1297 clusterP = pkg->GetPSideCluster(k);
1298 SN[k]= n1*SP[k]/sumsig;
1299 XP[k]=GetClusterZ(clusterP);
1300 XPerr[k]= clusterP->GetPositionError();
1301 SNerr[k] = n1err*SP[k]/sumsig;
1302 CreateNewRecPoint(XP[k],XPerr[k], XN,XNerr, SP[k]+SN[k], SPerr[k]+SNerr[k],clusterP, clusterN, 1.0);
1303 }
1304
3c9a27a7 1305 delete [] XP;
1306 delete [] XPerr;
1307 delete [] SP;
1308 delete [] SN;
1309 delete [] SPerr;
1310 delete [] SNerr;
b0f5e3fc 1311
1312}
1313
1314/**********************************************************************/
1315/********* 2 X 2 *********************************************/
1316/**********************************************************************/
1317
1318void AliITSClusterFinderSSD::
1319ResolveTwoForTwoPackage(AliITSpackageSSD *pkg)
1320{
766ef0c8 1321// resolve 2x2 packages
b0f5e3fc 1322
1323 AliITSclusterSSD *clusterP1 = pkg->GetPSideCluster(0);
1324 AliITSclusterSSD *clusterP2 = pkg->GetPSideCluster(1);
1325 AliITSclusterSSD *clusterN1 = pkg->GetNSideCluster(0);
1326 AliITSclusterSSD *clusterN2 = pkg->GetNSideCluster(1);
1327
1328 AliITSclusterSSD *tmp;
1329
1330 Float_t p1sig;
1331 Float_t p2sig;
1332 Float_t n1sig;
1333 Float_t n2sig;
1334
1335 Float_t p1sigErr;
1336 Float_t p2sigErr;
1337 Float_t n1sigErr;
1338 Float_t n2sigErr;
1339
1340 Float_t ZposP1;
1341 Float_t ZposP2;
1342 Float_t ZposN1;
1343 Float_t ZposN2;
1344
1345 Float_t ZposErrP1;
1346 Float_t ZposErrP2;
1347 Float_t ZposErrN1;
1348 Float_t ZposErrN2;
1349
1350 Float_t D12;
1351
1352 Float_t Chicomb1;
1353 Float_t Chicomb2;
1354
1355 if(clusterP1->GetDigitStripNo(0) > clusterP2->GetDigitStripNo(0)) {
1356 if (debug) cout<<"P strips flip\n";
1357 tmp = clusterP1;
1358 clusterP1 = clusterP2;
1359 clusterP2 = tmp;
1360 }
1361 if(clusterN1->GetDigitStripNo(0) > clusterN2->GetDigitStripNo(0)) {
1362 if (debug) cout<<"N strips flip\n";
1363 tmp = clusterN1;
1364 clusterN1 = clusterN2;
1365 clusterN2 = tmp;
1366 }
1367 p1sig = clusterP1->GetTotalSignal();
1368 p2sig = clusterP2->GetTotalSignal();
1369 n1sig = clusterN1->GetTotalSignal();
1370 n2sig = clusterN2->GetTotalSignal();
1371
1372
1373 p1sigErr = clusterP1->GetTotalSignalError();
1374 n1sigErr = clusterN1->GetTotalSignalError();
1375 p2sigErr = clusterP2->GetTotalSignalError();
1376 n2sigErr = clusterN2->GetTotalSignalError();
1377
1378 ZposP1 = clusterP1->GetPosition();
1379 ZposN1 = clusterN1->GetPosition();
1380 ZposP2 = clusterP2->GetPosition();
1381 ZposN2 = clusterN2->GetPosition();
1382
1383 ZposErrP1 = clusterP1->GetPositionError();
1384 ZposErrN1 = clusterN1->GetPositionError();
1385 ZposErrP2 = clusterP2->GetPositionError();
1386 ZposErrN2 = clusterN2->GetPositionError();
1387
1388 //in this may be two types:
1389 // 1.When each cluster crosses with 2 clusters on the other side
1390 // gives 2 ghosts and 2 real points
1391 //
1392 // 2.When two clusters (one an each side) crosses with only one on
1393 // the other side and two crosses (one on the each side) with
1394 // two on the other gives 2 or 3 points,
1395
1396 if (debug) cout<<"2 for 2 ambiguity ...";
1397 /***************************/
1398 /*First sort of combination*/
1399 /***************************/
1400
1401 if((clusterP1->GetCrossNo()==2) && (clusterN1->GetCrossNo()==2)) {
1402 if (debug) cout<<"All clusters has 2 crosses\n";
1403 D12 = TMath::Sqrt((Float_t)((p1sig -p2sig)*(p1sig -p2sig) + (n1sig -n2sig)*(n1sig -n2sig)));
1404
1405 if(debug) cout<<"Distance between points in P(N) plane D12 = "<<D12<<"\n";
1406
1407 /*********************************************/
1408 /*resolving ambiguities: */
1409 /*we count Chicomb's and we take combination */
1410 /*giving lower Chicomb as a real points */
1411 /*Keep only better combinantion */
1412 /*********************************************/
1413
766ef0c8 1414 if (D12 > (fAlpha3*17.8768)) {
b0f5e3fc 1415 if (debug) cout<<"decided to take only one pair \n";
1416 Chicomb1 = DistToPML(p1sig,n1sig) + DistToPML(p2sig,n2sig);
1417 Chicomb2 = DistToPML(p2sig,n1sig) + DistToPML(p1sig,n2sig);
1418 if (debug) {
1419 cout<<" 00 11 combination : "<<Chicomb1<<" 01 10 combination : "<<Chicomb2<<" \n";
1420 cout<<"p1 = "<<p1sig<<" n1 = "<<n1sig<<" p2 = "<<p2sig<<" n2 = "<<n2sig<<"\n";
1421 }
1422 if (Chicomb1 < Chicomb2) {
1423 if (debug) cout<<"00 11";
1424 CreateNewRecPoint(ZposP1,ZposErrP1, ZposN1,ZposErrN1, p1sig+n1sig, p1sigErr+n1sigErr, clusterP1, clusterN1, 0.75);
1425 CreateNewRecPoint(ZposP2,ZposErrP2, ZposN2,ZposErrN2, p2sig+n2sig, p2sigErr+n2sigErr, clusterP2, clusterN2, 0.75);
1426
1427
1428 //second cominantion
1429 } else {
1430 if (debug) cout<<"01 10";
1431 CreateNewRecPoint(ZposP1,0, ZposN2,0, p1sig+n2sig, p1sigErr+n2sigErr, clusterP1, clusterN2, 0.75);
1432 CreateNewRecPoint(ZposP2,0, ZposN1,0, p2sig+n1sig, p2sigErr+n1sigErr, clusterP2, clusterN1, 0.75);
1433 } //end second combinantion
766ef0c8 1434 //if (D12 > fAlpha3*17.8768)
b0f5e3fc 1435 //keep all combinations
1436 } else {
1437 if (debug) cout<<"We decide to take all points\n";
1438 CreateNewRecPoint(ZposP1,ZposErrP1, ZposN1,ZposErrN1, p1sig+n1sig, p1sigErr+n1sigErr, clusterP1, clusterN1, 0.5);
1439 CreateNewRecPoint(ZposP2,ZposErrP2, ZposN2,ZposErrN2, p2sig+n2sig, p2sigErr+n2sigErr, clusterP2, clusterN2, 0.5);
1440 CreateNewRecPoint(ZposP1,ZposErrP1, ZposN2,ZposErrN2, p1sig+n2sig, p1sigErr+n2sigErr, clusterP1, clusterN2, 0.5);
1441 CreateNewRecPoint(ZposP2,ZposErrP2, ZposN1,ZposErrN1, p2sig+n1sig, p2sigErr+n1sigErr, clusterP2, clusterN1, 0.5);
1442 }
1443 }
1444 // ad.2 Second type of combination
1445 else {
1446 Chicomb1 = DistToPML(p1sig,n1sig) + DistToPML(p2sig,n2sig);
1447 if (debug) cout<<"\nhere can be reconstructed 3 points: chicomb = "<<Chicomb1<<"\n";
1448
766ef0c8 1449 if (Chicomb1<fAlpha1) {
b0f5e3fc 1450 if (debug) cout<<"\nWe decided to take 3rd point";
1451 if (clusterP1->GetCrossNo()==1) {
1452 if (debug) cout<<"... P1 has one cross\n";
be33dccb 1453 n1sig = p1sig/fPNsignalRatio;
1454 p2sig = n2sig*fPNsignalRatio;
b0f5e3fc 1455
1456 clusterN1->CutTotalSignal(n1sig);
1457 clusterP2->CutTotalSignal(p2sig);
1458
1459 CreateNewRecPoint(ZposP2,ZposErrP2, ZposN1,ZposErrN1, p2sig+n1sig, p2sigErr+n1sigErr, clusterP2, clusterN1, 0.5);
1460
1461 n1sig = clusterN1->GetTotalSignal();
1462 p2sig = clusterP2->GetTotalSignal();
1463
1464 } else {
1465 if (debug) cout<<"... N1 has one cross\n";
1466
be33dccb 1467 n2sig=p2sig/fPNsignalRatio;
1468 p1sig=n1sig*fPNsignalRatio;
b0f5e3fc 1469
1470 clusterN2->CutTotalSignal(n2sig);
1471 clusterP1->CutTotalSignal(p1sig);
1472
1473 CreateNewRecPoint(ZposP1,ZposErrP1, ZposN2,ZposErrN2, p1sig+n2sig, p1sigErr+n2sigErr, clusterP1, clusterN2, 0.5);
1474
1475 n2sig=clusterN2->GetTotalSignal();
1476 p1sig=clusterP1->GetTotalSignal();
1477 }
1478 } else {
1479 if (debug) cout<<"\nWe decided NOT to take 3rd point\n";
1480 }
1481
1482 CreateNewRecPoint(ZposP1,ZposErrP1, ZposN1,ZposErrN1, p1sig+n1sig, p1sigErr+n1sigErr, clusterP1, clusterN1, 1.0);
1483 CreateNewRecPoint(ZposP2,ZposErrP2, ZposN2,ZposErrN2, p2sig+n2sig, p2sigErr+n2sigErr, clusterP2, clusterN2, 1.0);
1484
1485 }
1486
1487}
1488
1489
1490
1491//------------------------------------------------------
1492Bool_t AliITSClusterFinderSSD::
1493CreateNewRecPoint(Float_t P, Float_t dP, Float_t N, Float_t dN,
1494 Float_t Sig,Float_t dSig,
1495 AliITSclusterSSD *clusterP, AliITSclusterSSD *clusterN,
1496 Stat_t prob)
1497{
766ef0c8 1498// create the recpoints
b0f5e3fc 1499
1500 const Float_t kdEdXtoQ = 2.778e+8;
1501 const Float_t kconv = 1.0e-4;
1502 const Float_t kRMSx = 20.0*kconv; // microns->cm ITS TDR Table 1.3
1503 const Float_t kRMSz = 830.0*kconv; // microns->cm ITS TDR Table 1.3
1504
1505
1506 Float_t p=P;
1507 Float_t n=N;
e8189707 1508 Int_t stripP, stripN;
1509 Int_t sigP, sigN;
1510 AliITSdigitSSD *digP, *digN;
b0f5e3fc 1511 if (GetCrossing(P,N)) {
1512 GetCrossingError(dP,dN);
1513 AliITSRawClusterSSD cnew;
1514 Int_t nstripsP=clusterP->GetNumOfDigits();
1515 Int_t nstripsN=clusterN->GetNumOfDigits();
1516 cnew.fMultiplicity=nstripsP;
1517 cnew.fMultiplicityN=nstripsN;
1518 /*
1519 if (nstripsP > 100) {
1520 printf("multiplicity > 100 - increase the dimension of the arrays %d\n",nstripsP);
1521 nstripsP=100;
1522 }
e8189707 1523 for(int i=0;i<nstripsP;i++) {
b0f5e3fc 1524 // check if 'clusterP->GetDigitStripNo(i)' returns the digit index
1525 cnew.fIndexMap[i] = clusterP->GetDigitStripNo(i);
1526 }
1527 if (nstripsN > 100) {
1528 printf("multiplicity > 100 - increase the dimension of the arrays %d\n",nstripsN);
1529 nstripsN=100;
1530 }
e8189707 1531 for(int i=0;i<nstripsN;i++) {
b0f5e3fc 1532 // check if 'clusterN->GetDigitStripNo(i)' returns the digit index
1533 cnew.fIndexMapN[i] = clusterN->GetDigitStripNo(i);
1534 }
1535 */
1536 cnew.fQErr=dSig;
1537 //cnew.fProbability=(float)prob;
1538 fITS->AddCluster(2,&cnew);
e8189707 1539 fSegmentation->GetPadIxz(P,N,stripP,stripN);
1540 digP = (AliITSdigitSSD*)fMap->GetHit(1,stripP);
1541 digN = (AliITSdigitSSD*)fMap->GetHit(0,stripN);
3c9a27a7 1542 printf("SSD: digP digN %p %p\n",digP,digN);
1543 if(digP) sigP = digP->fSignal;
1544 else sigP=0;
1545 if(digN) sigN = digN->fSignal;
1546 else sigN=0;
1547 if (!digP && !digN) {
1548 Error("CreateNewRecPoint","cannot find the digit!");
1549 return kFALSE;
1550 }
b0f5e3fc 1551 // add the rec point info
1552 AliITSRecPoint rnew;
1553 rnew.SetX(P*kconv);
1554 rnew.SetZ(N*kconv);
1555 rnew.SetQ(Sig);
1556 rnew.SetdEdX(Sig/kdEdXtoQ);
1557 rnew.SetSigmaX2( kRMSx* kRMSx);
1558 rnew.SetSigmaZ2( kRMSz* kRMSz);
3c9a27a7 1559 //rnew.SetProbability((float)prob);
e8189707 1560 if(sigP > sigN) {
1561 rnew.fTracks[0]=digP->fTracks[0];
1562 rnew.fTracks[1]=digP->fTracks[1];
1563 rnew.fTracks[2]=digP->fTracks[2];
3c9a27a7 1564 //printf("sigP>sigN: %d %d %d\n",digP->fTracks[0],digP->fTracks[1],digP->fTracks[2]);
e8189707 1565 } else {
1566 rnew.fTracks[0]=digN->fTracks[0];
1567 rnew.fTracks[1]=digN->fTracks[1];
1568 rnew.fTracks[2]=digN->fTracks[2];
3c9a27a7 1569 //printf("sigN>sigP: %d %d %d\n",digN->fTracks[0],digN->fTracks[1],digN->fTracks[2]);
e8189707 1570 }
3c9a27a7 1571 //printf("SSD: track1 track2 track3 X Z %d %d %d %f %f\n",rnew.fTracks[0],rnew.fTracks[1],rnew.fTracks[2],P*kconv,N*kconv);
b0f5e3fc 1572 fITS->AddRecPoint(rnew);
1573 /*
1574 // it was
1575 fPointsM->AddLast( (TObject*)
1576 new AliITSRecPointSSD(P,dP,0,0,N,dN,Sig,dSig,fLayer,fLad,fDet,clusterP,clusterN,prob) );
1577 */
1578
1579 if(debug) cout<<"\n"<<": ImpactPoint Created: X = "<<P<<" Z = "<<N<<"; P = "<<p<<" N = "<<n<<"\n";
1580 return kTRUE;
1581 } else {
1582 if (debug) {
1583 cout<<"\n"<<": ATENTION : stupid ImpactPoit Point X = "<<P<<" Z = "<<N<<"; P = "<<p<<" N = "<<n<<"\n";
1584 }
1585 return kFALSE;
1586 }
1587}
1588
1589
1590
1591/**********************************************************************/
1592Bool_t AliITSClusterFinderSSD::
1593CreateNewRecPoint(AliITSclusterSSD *clusterP, AliITSclusterSSD *clusterN, Stat_t prob)
1594{
766ef0c8 1595// create recpoints
1596
b0f5e3fc 1597 Float_t posClusterP; //Cluster P position in strip coordinates
1598 Float_t posClusterN; //Cluster N position in strip coordinates
1599
1600 Float_t posErrorClusterP;
1601 Float_t posErrorClusterN;
1602
1603 Float_t sigClusterP;
1604 Float_t sigClusterN;
1605
1606 Float_t sigErrorClusterP;
1607 Float_t sigErrorClusterN;
1608
1609 posClusterP = clusterP->GetPosition();
1610 posErrorClusterP = clusterP->GetPositionError();
1611 sigClusterP = clusterP->GetTotalSignal();
1612 sigErrorClusterP = clusterP->GetTotalSignalError();
1613
1614 posClusterN = clusterN->GetPosition();
1615 posErrorClusterN = clusterN->GetPositionError();
1616 sigClusterN = clusterN->GetTotalSignal();
1617 sigErrorClusterN = clusterN->GetTotalSignalError();
1618
1619 return CreateNewRecPoint( posClusterP,posErrorClusterP, posClusterN,posErrorClusterN,
1620 sigClusterP+sigClusterN, sigErrorClusterN+sigErrorClusterP,
1621 clusterP, clusterN, prob);
1622
1623}
1624
1625//--------------------------------------------------
1626Bool_t AliITSClusterFinderSSD::IsCrossing(AliITSclusterSSD* p, AliITSclusterSSD* n)
1627{
766ef0c8 1628// check for crossings
1629
b0f5e3fc 1630 Float_t x = p->GetPosition();
1631 Float_t y = n->GetPosition();
1632 return GetCrossing(x,y);
1633}
1634
1635
1636//----------------------------------------------
1637Bool_t AliITSClusterFinderSSD::Strip2Local( Float_t stripP, Float_t stripN, Float_t &Z,Float_t &X)
1638{
766ef0c8 1639// convert between strip numbers and local coordinates
b0f5e3fc 1640
1641/*
1642 Z = (stripN-stripP)*fFactorOne;
1643 X = (stripN + stripP - fStripZeroOffset)*fFactorTwo;
1644*/
1645
1646 Float_t P =stripP;
1647 Float_t N =stripN;
1648
1649 GetCrossing(P,N);
1650 X=P;
1651 Z=N;
1652 if (debug) cout<<"P="<<stripP<<" N="<<stripN<<" X = "<<X<<" Z = "<<Z<<"\n";
e8189707 1653 if ((Z<2.1)&&(Z>-2.1)&&(X<3.65)&&(X>-3.65)) return kTRUE;
1654 else return kFALSE;
b0f5e3fc 1655
1656}
1657
1658
1659//-----------------------------------------------------------
1660Float_t AliITSClusterFinderSSD::GetClusterZ(AliITSclusterSSD* clust)
1661{
766ef0c8 1662// get Z coordinate of the cluster
b0f5e3fc 1663
1664 return clust->GetPosition();
1665
1666}
1667
1668//-------------------------------------------------------------
1669Int_t AliITSClusterFinderSSD::GetBestComb
1670(Int_t** comb,Int_t Ncomb, Int_t Ncl, AliITSpackageSSD * pkg)
1671{
b0f5e3fc 1672//returns index of best combination in "comb"
1673//comb : sets of combinations
1674// in given combination on place "n" is index of
1675// Nside cluster coresponding to "n" P side cluster
1676//
1677//Ncomb : number of combinations == number of rows in "comb"
1678//
1679//Ncl : number of clusters in each row == number of columns in "comb"
1680//
1681//pkg : package
1682
1683
1684 Float_t currbval=-1; //current best value,
1685 //starting value is set to number negative,
1686 //to be changed in first loop turn automatically
1687
1688 Int_t curridx=0; //index of combination giving currently the best chi^2
1689 Float_t chi;
1690 Float_t ps, ns; //signal of P cluster and N cluster
1691
e8189707 1692 for (Int_t i=0;i<Ncomb;i++)
b0f5e3fc 1693 {
1694 chi=0;
1695
e8189707 1696 for (Int_t j=0;j<Ncl;j++)
b0f5e3fc 1697 {
1698 ps = pkg->GetPSideCluster(j)->GetTotalSignal(); //carrefully here, different functions
1699 ns = GetNSideCluster(comb[i][j])->GetTotalSignal();
1700 chi+=DistToPML(ps,ns);
1701 }
1702
1703 if (currbval==-1) currbval = chi;
1704 if (chi<currbval)
1705 {
1706 curridx = i;
1707 currbval =chi;
1708 }
1709 }
1710
1711
1712 return curridx;
1713
1714}
1715
1716//--------------------------------------------------
1717void AliITSClusterFinderSSD::GetBestMatchingPoint
1718(Int_t & ip, Int_t & in, AliITSpackageSSD* pkg )
1719{
766ef0c8 1720// find the best matching point
b0f5e3fc 1721
1722 if (debug) pkg->PrintClusters();
1723
1724 Float_t ps, ns; //signals on side p & n
1725
1726 Int_t nc; // number of crosses in given p side cluster
1727 Int_t n,p;
1728 AliITSclusterSSD *curPcl, *curNcl; //current p side cluster
1729 Float_t bestchi, chi;
1730 ip=p=pkg->GetPSideClusterIdx(0);
1731 in=n=pkg->GetNSideClusterIdx(0);
1732
1733 bestchi=DistToPML( pkg->GetPSideCluster(0)->GetTotalSignal(),
1734 pkg->GetNSideCluster(0)->GetTotalSignal() );
e8189707 1735 for (Int_t i = 0; i< pkg->GetNumOfClustersP(); i++)
b0f5e3fc 1736 {
1737
1738 p = pkg->GetPSideClusterIdx(i);
1739 curPcl= GetPSideCluster(p);
1740 nc=curPcl->GetCrossNo();
1741 ps=curPcl->GetTotalSignal();
1742
e8189707 1743 for (Int_t j = 0; j< nc; j++)
b0f5e3fc 1744 {
1745 n=curPcl->GetCross(j);
1746 curNcl= GetNSideCluster(n);
1747 ns=curNcl->GetTotalSignal();
1748 chi = DistToPML(ps, ns);
1749
1750 if (chi <= bestchi)
1751 {
1752 bestchi = chi;
1753 in = n;
1754 ip = p;
1755 }
1756 }
1757 }
1758
1759}
1760
1761
1762//------------------------------------------------------
1763void AliITSClusterFinderSSD::CalcStepFactor(Float_t Psteo, Float_t Nsteo )
1764{
766ef0c8 1765// calculate the step factor for matching clusters
b0f5e3fc 1766
1767
1768 // 95 is the pitch, 4000 - dimension along z ?
1769 /*
1770 fSFF = ( (Int_t) (Psteo*40000/95 ) );// +1;
1771 fSFB = ( (Int_t) (Nsteo*40000/95 ) );// +1;
1772 */
1773
1774
1775 Float_t dz=fSegmentation->Dz();
1776
1777 fSFF = ( (Int_t) (Psteo*dz/fPitch ) );// +1;
1778 fSFB = ( (Int_t) (Nsteo*dz/fPitch ) );// +1;
1779
1780}
1781
1782
1783//-----------------------------------------------------------
1784AliITSclusterSSD* AliITSClusterFinderSSD::GetPSideCluster(Int_t idx)
1785{
766ef0c8 1786// get P side clusters
b0f5e3fc 1787
1788
b0f5e3fc 1789
1790
1791 if((idx<0)||(idx>=fNClusterP))
1792 {
1793 printf("AliITSClusterFinderSSD::GetPSideCluster : index out of range\n");
1794 return 0;
1795 }
1796 else
1797 {
1798 return (AliITSclusterSSD*)((*fClusterP)[idx]);
1799 }
1800}
1801
1802//-------------------------------------------------------
1803AliITSclusterSSD* AliITSClusterFinderSSD::GetNSideCluster(Int_t idx)
1804{
766ef0c8 1805// get N side clusters
b0f5e3fc 1806
1807 if((idx<0)||(idx>=fNClusterN))
1808 {
1809 printf("AliITSClusterFinderSSD::GetNSideCluster : index out of range\n");
1810 return 0;
1811 }
1812 else
1813 {
1814 return (AliITSclusterSSD*)((*fClusterN)[idx]);
1815 }
1816}
1817
1818//--------------------------------------------------------
1819AliITSclusterSSD* AliITSClusterFinderSSD::GetCluster(Int_t idx, Bool_t side)
1820{
766ef0c8 1821// Get cluster
b0f5e3fc 1822
1823 return (side) ? GetPSideCluster(idx) : GetNSideCluster(idx);
1824}
1825
1826
1827//--------------------------------------------------------
1828void AliITSClusterFinderSSD::ConsumeClusters()
1829{
766ef0c8 1830// remove clusters from the list of unmatched clusters
e8189707 1831
1832 for ( Int_t i=0;i<fNPackages;i++)
b0f5e3fc 1833 {
1834 ((AliITSpackageSSD*)((*fPackages)[i]))->ConsumeClusters();
1835 }
1836
1837}
1838
1839
1840//--------------------------------------------------------
1841void AliITSClusterFinderSSD::ReconstructNotConsumedClusters()
1842{
766ef0c8 1843// reconstruct remaining non-crossing clusters
1844
b0f5e3fc 1845 Int_t i;
1846 AliITSclusterSSD *cluster;
1847 Float_t pos;
1848 Float_t sig;
1849 Float_t sigerr;
1850 Float_t x1,x2,z1,z2;
1851
1852 Float_t Dx = fSegmentation->Dx();
1853 Float_t Dz = fSegmentation->Dz();
1854
1855 const Float_t kdEdXtoQ = 2.778e+8; // GeV -> number of e-hole pairs
1856 const Float_t kconv = 1.0e-4;
1857 const Float_t kRMSx = 20.0*kconv; // microns->cm ITS TDR Table 1.3
1858 const Float_t kRMSz = 830.0*kconv; // microns->cm ITS TDR Table 1.3
1859
e8189707 1860 Int_t stripP,stripN;
1861 AliITSdigitSSD *dig;
1862 for (i=0;i<fNClusterP;i++)
b0f5e3fc 1863 {
3c9a27a7 1864 //printf("P side cluster: cluster number %d\n",i);
b0f5e3fc 1865 cluster = GetPSideCluster(i);
1866 if (!cluster->IsConsumed())
1867 {
1868 if( (pos = cluster->GetPosition()) < fSFB)
1869 {
1870 sig = cluster->GetTotalSignal();
1871
be33dccb 1872 sig += sig/fPNsignalRatio;
b0f5e3fc 1873
1874 sigerr = cluster->GetTotalSignalError();
1875 x1 = -Dx/2 + pos *fPitch;
1876 z1 = -Dz/2;
1877
1878 x2 = pos;
1879 z2 = 1;
1880 GetCrossing (x2,z2);
3c9a27a7 1881 fSegmentation->GetPadIxz((x1+x2)/2,(z1+z2)/2,stripP,stripN);
1882 dig = (AliITSdigitSSD*)fMap->GetHit(1,stripP);
1883 if (!dig) printf("SSD: check the digit! dig %p\n",dig);
1884 if(!dig) continue;
b0f5e3fc 1885 AliITSRawClusterSSD cnew;
1886 Int_t nstripsP=cluster->GetNumOfDigits();
1887 cnew.fMultiplicity=nstripsP;
1888 cnew.fMultiplicityN=0;
1889 /*
1890 if (nstripsP > 100) {
1891 printf("multiplicity > 100 - increase the dimension of the arrays %d\n",nstripsP);
1892 nstripsP=100;
1893 }
e8189707 1894 for(int k=0;k<nstripsP;k++) {
b0f5e3fc 1895 // check if 'clusterP->GetDigitStripNo(i)' returns
1896 // the digit index and not smth else
1897 cnew.fIndexMap[k] = cluster->GetDigitStripNo(k);
1898 }
1899 */
1900 cnew.fQErr=sigerr;
1901 //cnew.fProbability=0.75;
1902 fITS->AddCluster(2,&cnew);
1903 // add the rec point info
1904 AliITSRecPoint rnew;
1905 rnew.SetX(kconv*(x1+x2)/2);
1906 rnew.SetZ(kconv*(z1+z2)/2);
1907 rnew.SetQ(sig);
1908 rnew.SetdEdX(sig/kdEdXtoQ);
1909 rnew.SetSigmaX2( kRMSx* kRMSx);
1910 rnew.SetSigmaZ2( kRMSz* kRMSz);
3c9a27a7 1911 //rnew.SetProbability(0.75);
e8189707 1912 rnew.fTracks[0]=dig->fTracks[0];
1913 rnew.fTracks[1]=dig->fTracks[1];
1914 rnew.fTracks[2]=dig->fTracks[2];
3c9a27a7 1915 //printf("digP: %d %d %d\n",dig->fTracks[0],dig->fTracks[1],dig->fTracks[2]);
1916 //printf("SSD - P : track1 track2 track3 X Z %d %d %d %f %f\n",rnew.fTracks[0],rnew.fTracks[1],rnew.fTracks[2],kconv*(x1+x2)/2,kconv*(z1+z2)/2);
b0f5e3fc 1917 fITS->AddRecPoint(rnew);
1918 /*
1919 fPointsM->AddLast( (TObject*)
1920 new AliITSRecPointSSD((x1+x2)/2 ,0,0,0,(z1+z2)/2,0,sig,sigerr,fLayer,fLad,fDet,cluster, 0.75) );
1921 */
1922
1923 if(debug) cout<<"\n"<<": SINGLE SIDE ImpactPoint Created: X = "
1924 <<(x1+x2)/2<<" Z = "<<(z1+z2)/2<<"; Pos = "<<pos;
1925 }
1926 }
1927 }
1928
e8189707 1929 for (i=0;i<fNClusterN;i++)
b0f5e3fc 1930 {
3c9a27a7 1931 //printf("N side cluster: cluster number %d\n",i);
b0f5e3fc 1932 cluster = GetNSideCluster(i);
1933 if (!cluster->IsConsumed())
1934 {
1935 if( (pos = cluster->GetPosition()) < fSFF)
1936 {
1937 sig = cluster->GetTotalSignal();
1938
be33dccb 1939 sig += sig*fPNsignalRatio;
b0f5e3fc 1940
1941 sigerr = cluster->GetTotalSignalError();
1942
1943 x1 = -Dx/2 + pos *fPitch;
1944 z1 = Dz/2;
1945
1946 x2 = 1;
1947 z2 = pos;
1948
1949 GetCrossing (x2,z2);
3c9a27a7 1950 fSegmentation->GetPadIxz((x1+x2)/2,(z1+z2)/2,stripP,stripN);
1951 dig = (AliITSdigitSSD*)fMap->GetHit(0,stripN);
1952 if (!dig) printf("SSD: check the digit! dig %p\n",dig);
1953 if(!dig) continue;
b0f5e3fc 1954 AliITSRawClusterSSD cnew;
1955 Int_t nstripsN=cluster->GetNumOfDigits();
1956 cnew.fMultiplicity=0;
1957 cnew.fMultiplicityN=nstripsN;
1958 /*
1959 if (nstripsN > 100) {
1960 printf("multiplicity > 100 - increase the dimension of the arrays %d\n",nstripsN);
1961 nstripsN=100;
1962 }
e8189707 1963 for(int k=0;k<nstripsN;k++) {
b0f5e3fc 1964 // check if 'clusterP->GetDigitStripNo(i)' returns
1965 // the digit index and not smth else
1966 cnew.fIndexMapN[k] = cluster->GetDigitStripNo(k);
1967 }
1968 */
1969 cnew.fQErr=sigerr;
1970 //cnew.fProbability=0.75;
1971 fITS->AddCluster(2,&cnew);
b0f5e3fc 1972 AliITSRecPoint rnew;
1973 rnew.SetX(kconv*(x1+x2)/2);
1974 rnew.SetZ(kconv*(z1+z2)/2);
1975 rnew.SetQ(sig);
1976 rnew.SetdEdX(sig/kdEdXtoQ);
1977 rnew.SetSigmaX2( kRMSx* kRMSx);
1978 rnew.SetSigmaZ2( kRMSz* kRMSz);
3c9a27a7 1979 //rnew.SetProbability(0.75);
e8189707 1980 rnew.fTracks[0]=dig->fTracks[0];
1981 rnew.fTracks[1]=dig->fTracks[1];
1982 rnew.fTracks[2]=dig->fTracks[2];
3c9a27a7 1983 //printf("digN: %d %d %d\n",dig->fTracks[0],dig->fTracks[1],dig->fTracks[2]);
1984 //printf("SSD -N : track1 track2 track3 X Z %d %d %d %f %f\n",rnew.fTracks[0],rnew.fTracks[1],rnew.fTracks[2],kconv*(x1+x2)/2,kconv*(z1+z2)/2);
b0f5e3fc 1985 fITS->AddRecPoint(rnew);
1986 /*
1987 fPointsM->AddLast( (TObject*)
1988 new AliITSRecPointSSD((x1+x2)/2 ,0,0,0,(z1+z2)/2,0,sig,sigerr,fLayer,fLad,fDet,cluster, 0.75) );
1989 */
1990
1991 if(debug) cout<<"\n"<<": SINGLE SIDE ImpactPoint Created: X = "
1992 <<(x1+x2)/2<<" Z = "<<(z1+z2)/2<<"; Pos = "<<pos;
1993
1994 }
1995 }
1996 }
1997
1998
1999}
2000
2001//_______________________________________________________________________
2002
2003Bool_t AliITSClusterFinderSSD::GetCrossing (Float_t &P, Float_t &N)
2004{
766ef0c8 2005// get crossing
b0f5e3fc 2006
2007 Float_t Dx = fSegmentation->Dx();
2008 Float_t Dz = fSegmentation->Dz();
2009
2010 //P==X
2011 //N==Z
2012
2013 Float_t dx = 0.1;
2014 P *= fPitch;
2015 N *= fPitch;
2016
2017 //P = (kZ * fTan + N + P)/2.0; // x coordinate
2018 //N = kZ - (P-N)/fTan; // z coordinate
2019
2020 if(fTanP + fTanN == 0) return kFALSE;
2021
2022 N = (N - P + fTanN * Dz) / (fTanP + fTanN); // X coordinate
2023 P = P + fTanP * N; // Y coordinate
2024
2025 P -= Dx/2;
2026 N -= Dz/2;
2027
2028 // N = -(N - P + kZ/2*(fTanP + fTanN))/(fTanP + fTanN);
2029 // P = (fTanP*(N-kX/2-kZ*fTanN/2) + fTanN*(P-kX/2-kZ*fTanP/2) )/(fTanP + fTanN);
2030
2031 if ( ( N < -(Dz/2+dx) ) || ( N > (Dz/2+dx) ) ) return kFALSE;
2032 if ( ( P < -(Dx/2+dx) ) || ( P > (Dx/2+dx) ) ) return kFALSE;
2033
2034 return kTRUE;
2035}
2036
2037
2038//_________________________________________________________________________
2039
2040void AliITSClusterFinderSSD::GetCrossingError(Float_t& dP, Float_t& dN)
2041{
766ef0c8 2042// get crossing error
2043
b0f5e3fc 2044 Float_t dz, dx;
2045
2046 dz = TMath::Abs(( dP + dN )*fPitch/(fTanP + fTanN) );
2047 dx = fPitch*(TMath::Abs(dP*(1 - fTanP/(fTanP + fTanN))) +
2048 TMath::Abs(dN *fTanP/(fTanP + fTanN) ));
2049
2050 dN = dz;
2051 dP = dx;
2052}