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