]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - ITS/AliITSClusterFinderSSD.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Fix bug in FillDigitsIndex
[u/mrichter/AliRoot.git] / ITS / AliITSClusterFinderSSD.cxx
... / ...
CommitLineData
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.
21
22//Piotr Krzysztof Skowronski
23//Warsaw University of Technology
24//skowron@if.pw.edu.pl
25
26*/
27
28/*
2915 -18 V 2000
30Eroor counting routines
31Automatic combination routines improved (traps)
32
33*/
34
35#include <iostream.h>
36#include <TArrayI.h>
37
38
39#include "AliRun.h"
40#include "AliITS.h"
41#include "AliITSdigit.h"
42#include "AliITSRawCluster.h"
43#include "AliITSRecPoint.h"
44#include "AliITSMapA1.h"
45#include "AliITSClusterFinderSSD.h"
46#include "AliITSclusterSSD.h"
47#include "AliITSpackageSSD.h"
48#include "AliITSsegmentation.h"
49
50const Bool_t AliITSClusterFinderSSD::fgkSIDEP=kTRUE;
51const Bool_t AliITSClusterFinderSSD::fgkSIDEN=kFALSE;
52
53const Int_t debug=0;
54
55ClassImp(AliITSClusterFinderSSD)
56
57//____________________________________________________________________
58//
59// Constructor
60//____________________________________________________________________
61//
62
63
64AliITSClusterFinderSSD::AliITSClusterFinderSSD(AliITSsegmentation *seg, TClonesArray *digits)
65{
66//Standard constructor
67
68 fSegmentation=seg;
69 fDigits=digits;
70
71 fMap = new AliITSMapA1(fSegmentation,fDigits);
72
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);
101
102 fPNsignalRatio=7./8.; // warning: hard-wired number
103
104}
105
106//-------------------------------------------------------
107AliITSClusterFinderSSD::~AliITSClusterFinderSSD() {
108// Default destructor
109
110 delete fClusterP;
111 delete fClusterN;
112 delete fPackages;
113 delete fDigitsIndexP;
114 delete fDigitsIndexN;
115
116 delete fMap;
117
118}
119
120//-------------------------------------------------------
121void AliITSClusterFinderSSD::InitReconstruction()
122{
123// initialization of the cluster finder
124
125 register Int_t i; //iterator
126
127 for (i=0;i<fNClusterP;i++)
128 {
129 fClusterP->RemoveAt(i);
130 }
131 fNClusterP =0;
132 for (i=0;i<fNClusterN;i++)
133 {
134 fClusterN->RemoveAt(i);
135 }
136 fNClusterN=0;
137
138 for (i=0;i<fNPackages;i++)
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{
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
168 fMap->FillMap();
169 FillDigitsIndex();
170 SortDigits();
171 FindNeighbouringDigits(); //ad. 2
172 SeparateOverlappedClusters(); //ad. 3
173 ClustersToPackages(); //ad. 4
174 ConsumeClusters();
175 PackagesToPoints(); //ad. 5
176 ReconstructNotConsumedClusters();
177
178 fMap->ClearMap();
179}
180
181
182//-------------------------------------------------
183void AliITSClusterFinderSSD::FindNeighbouringDigits()
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
189 register Int_t i;
190
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 = *(Digits());
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
207 for (i=1; i<fNDigitsP; i++) {
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
216 new(lClusterP[fNClusterP++]) AliITSclusterSSD(dnumber,dbuffer,Digits(),fgkSIDEP);
217 dbuffer[0]=lDigitsIndexP[i];
218 dnumber = 1;
219 }
220 } // end loop over fNDigitsP
221 new(lClusterP[fNClusterP++]) AliITSclusterSSD(dnumber,dbuffer,Digits(),fgkSIDEP);
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
229 for (i=1; i<fNDigitsN; i++) {
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 {
235 new(lClusterN[fNClusterN++]) AliITSclusterSSD(dnumber,dbuffer,Digits(),fgkSIDEN);
236 dbuffer[0]=lDigitsIndexN[i];
237 dnumber = 1;
238 }
239 } // end loop over fNDigitsN
240 new(lClusterN[fNClusterN++]) AliITSclusterSSD(dnumber,dbuffer,Digits(),fgkSIDEN);
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{
253// overlapped clusters separation
254
255 register Int_t i; //iterator
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
272 for (i=0;i<initPsize;i++)
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);
277 for (Int_t j=1; j<nj; j++)
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
289 if(numerofsplits>0) SplitCluster(splitlist,numerofsplits,i,fgkSIDEP);
290 numerofsplits=0;
291
292 //in signed places (splitlist)
293 } // end loop over clusters on Pside
294
295 for (i=0;i<initNsize;i++) {
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);
299 for (Int_t j=1; j<nj; j++)
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
309 if(numerofsplits>0) SplitCluster(splitlist,numerofsplits,i,fgkSIDEN);
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{
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])) ;
338 for (i = nsplits; i>0 ;i--) {
339 NN=curentcluster->SplitCluster((*list)[(i-1)],tmpdigits);
340 new ((*fClusterP)[fNClusterP]) AliITSclusterSSD(NN,tmpdigits,Digits(),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]));
353 for (i = nsplits; i>0 ;i--) {
354 NN=curentcluster->SplitCluster((*list)[(i-1)],tmpdigits);
355 new ((*fClusterN)[fNClusterN]) AliITSclusterSSD(NN,tmpdigits,Digits(),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{
373// sort digits on the P side
374
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*)((*(Digits()))[(*fDigitsIndexP)[start]]))->GetStripNumber();
387 right= ((AliITSdigitSSD*)((*(Digits()))[(*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{
404// sort digits on the N side
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*)((*(Digits()))[(*fDigitsIndexN)[start]]))->GetStripNumber();
417 right=((AliITSdigitSSD*)((*(Digits()))[(*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{
432 //Fill the indexes of the clusters belonging to a given ITS module
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
448 for ( i = 0 ; i< N; i++ ) {
449 dig = (AliITSdigitSSD*)GetDigit(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
455 for( k=0;k<PNs;k++)
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
472 for( k=0;k<NNs;k++)
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
488 delete [] PSidx;
489 delete [] NSidx;
490
491 if (debug) cout<<"Digits : P = "<<fNDigitsP<<" N = "<<fNDigitsN<<endl;
492
493}
494
495
496//-------------------------------------------
497
498void AliITSClusterFinderSSD::SortDigits()
499{
500// sort digits
501
502 Int_t i;
503 if(fNDigitsP>1)
504 for (i=0;i<fNDigitsP-1;i++)
505 if (SortDigitsP(0,(fNDigitsP-1-i))==0) break;
506
507 if(fNDigitsN>1)
508 for (i=0;i<fNDigitsN-1;i++)
509 if(SortDigitsN(0,(fNDigitsN-1-i))==0) break;
510}
511
512
513
514//----------------------------------------------
515void AliITSClusterFinderSSD::FillClIndexArrays(Int_t* arrayP, Int_t *arrayN)
516{
517// fill cluster index array
518
519 register Int_t i;
520 for (i=0; i<fNClusterP;i++)
521 {
522 arrayP[i]=i;
523 }
524 for (i=0; i<fNClusterN;i++)
525 {
526 arrayN[i]=i;
527 }
528}
529
530
531//------------------------------------------------------
532void AliITSClusterFinderSSD::SortClusters(Int_t* arrayP, Int_t *arrayN)
533{
534// sort clusters
535
536 Int_t i;
537 if(fNClusterP>1)
538 for (i=0;i<fNClusterP-1;i++)
539 if (SortClustersP(0,(fNClusterP-1),arrayP)==0) break;
540
541
542 if(fNClusterN>1)
543 for (i=0;i<fNClusterN-1;i++)
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{
553//Sort P side clusters
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{
583//Sort N side clusters
584
585
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{
613// fill packages
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
623 register Int_t i; //iterator
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
649 Int_t k;
650 for (i=0;i<fNClusterP;i++) {
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
678 for (Int_t k=1;k<fSFF+1;k++) {
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;
731 for (Int_t k=1;k<fSFF+1;k++) {
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{
775// find recpoints from packages
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
784 for (i=0;i<fNPackages;i++) {
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{
853// resolve clusters with one cross
854/*
855
856ie:
857 \ /
858 \ /
859 \/
860 /\
861 / \
862 / \
863 / \
864 / \
865 / \
866*/
867
868
869
870 if (clSide == fgkSIDEP) ResolvePClusterWithOneCross(currentpkg,clusterIndex);
871 else ResolveNClusterWithOneCross(currentpkg,clusterIndex);
872
873}
874
875
876//---------------------------------------------------------
877void AliITSClusterFinderSSD::
878ResolvePClusterWithOneCross(AliITSpackageSSD *pkg, Int_t clusterIndex)
879{
880// resolve clusters with one cross on P side
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();
933 pkg->DelCluster(clusterIndex,fgkSIDEP);
934 sigClusterN = ps/fPNsignalRatio;
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;
939 if (debug) cout<<"n1 < p1/fPNsignalRatio";
940 if (debug) cout<<"Attempting to del cluster N "<<clusterIdx<<" ... \n";
941 pkg->DelClusterOI(clusterIdx,fgkSIDEN);
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";
946 if (Chicomb > fAlpha2) {
947 //it is near, so we can risk throwing this cluster away too
948 if (debug) cout<<"Attempting to del cluster N "<<clusterIdx<<"...\n";
949 pkg->DelClusterOI(clusterIdx,fgkSIDEN);
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{
967// resolve clusters with one cross on N side
968
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();
1004 pkg->DelCluster(clusterIndex,fgkSIDEN);
1005 sigClusterP=ns*fPNsignalRatio;
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;
1010 if (debug) cout<<"ps < ns*fPNsignalRatio";
1011 if (debug) cout<<"Attempting to del cluster P "<<clusterIdx<<" ... \n";
1012 pkg->DelClusterOI(clusterIdx,fgkSIDEP);
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";
1017 if (Chicomb > fAlpha2) {
1018 //it is near, so we can risk frowing this cluster away too
1019 if (debug) cout<<"Attempting to del cluster P "<<clusterIdx<<"...\n";
1020 pkg->DelClusterOI(clusterIdx,fgkSIDEP);
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{
1041//find best combination
1042
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
1053
1054 Int_t NP = pkg->GetNumOfClustersP();
1055 for (itera =2; itera <= NP ;itera ++) {
1056 sizet=sizet*itera;
1057 if (sizet > 10000) {
1058 sizet=10000;
1059 break;
1060 }
1061 }
1062
1063 Int_t** combin = new Int_t*[sizet]; //2D array to keep combinations in
1064
1065 for (itera =0; itera <sizet;itera++) {
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
1079 for (itera =0; itera < NP; itera ++) {
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
1088 for (itera =0; itera <sizet;itera++) {
1089 delete [](combin[itera]);
1090 }
1091 delete [] combin;
1092 return kTRUE;
1093
1094}
1095
1096
1097//----------------------------------------------------
1098void AliITSClusterFinderSSD::
1099ResolveOneBestMatchingPoint(AliITSpackageSSD *pkg)
1100{
1101// find best matching point
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
1134 pkg->DelClusterOI(pi, fgkSIDEP);
1135 pkg->DelClusterOI(ni, fgkSIDEN);
1136
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{
1150// resolve simple package
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) {
1167// resolve P side clusters from packages
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
1217 for (k=0;k<NN;k++) {
1218 SN[k] = pkg->GetNSideCluster(k)->GetTotalSignal();
1219 SNerr[k] = pkg->GetNSideCluster(k)->GetTotalSignalError();
1220 sumsig += SN[k];
1221 }
1222 for (k=0;k<NN;k++) {
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
1232 delete [] XN;
1233 delete [] XNerr;
1234 delete [] SP;
1235 delete [] SN;
1236 delete [] SPerr;
1237 delete [] SNerr;
1238
1239
1240}
1241
1242
1243//---------------------------------------------------------
1244void AliITSClusterFinderSSD::ResolvePackageWithOneNSideCluster(AliITSpackageSSD *pkg) {
1245// resolve N side clusters from packages
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
1290 for (k=0;k<NP;k++) {
1291 SP[k] = pkg->GetPSideCluster(k)->GetTotalSignal();
1292 sumsig += SP[k];
1293 SPerr[k] = pkg->GetPSideCluster(k)->GetTotalSignalError();
1294 }
1295
1296 for (k=0;k<NP;k++) {
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
1305 delete [] XP;
1306 delete [] XPerr;
1307 delete [] SP;
1308 delete [] SN;
1309 delete [] SPerr;
1310 delete [] SNerr;
1311
1312}
1313
1314/**********************************************************************/
1315/********* 2 X 2 *********************************************/
1316/**********************************************************************/
1317
1318void AliITSClusterFinderSSD::
1319ResolveTwoForTwoPackage(AliITSpackageSSD *pkg)
1320{
1321// resolve 2x2 packages
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
1414 if (D12 > (fAlpha3*17.8768)) {
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
1434 //if (D12 > fAlpha3*17.8768)
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
1449 if (Chicomb1<fAlpha1) {
1450 if (debug) cout<<"\nWe decided to take 3rd point";
1451 if (clusterP1->GetCrossNo()==1) {
1452 if (debug) cout<<"... P1 has one cross\n";
1453 n1sig = p1sig/fPNsignalRatio;
1454 p2sig = n2sig*fPNsignalRatio;
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
1467 n2sig=p2sig/fPNsignalRatio;
1468 p1sig=n1sig*fPNsignalRatio;
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{
1498// create the recpoints
1499
1500
1501 const Float_t kdEdXtoQ = 2.778e+8;
1502 const Float_t kconv = 1.0e-4;
1503 const Float_t kRMSx = 20.0*kconv; // microns->cm ITS TDR Table 1.3
1504 const Float_t kRMSz = 830.0*kconv; // microns->cm ITS TDR Table 1.3
1505
1506
1507 Float_t p=P;
1508 Float_t n=N;
1509 Int_t stripP, stripN;
1510 Int_t sigP, sigN;
1511 AliITSdigitSSD *digP, *digN;
1512 if (GetCrossing(P,N)) {
1513 GetCrossingError(dP,dN);
1514 AliITSRawClusterSSD cnew;
1515 Int_t nstripsP=clusterP->GetNumOfDigits();
1516 Int_t nstripsN=clusterN->GetNumOfDigits();
1517 cnew.fMultiplicity=nstripsP;
1518 cnew.fMultiplicityN=nstripsN;
1519 /*
1520 if (nstripsP > 100) {
1521 printf("multiplicity > 100 - increase the dimension of the arrays %d\n",nstripsP);
1522 nstripsP=100;
1523 }
1524 for(int i=0;i<nstripsP;i++) {
1525 // check if 'clusterP->GetDigitStripNo(i)' returns the digit index
1526 cnew.fIndexMap[i] = clusterP->GetDigitStripNo(i);
1527 }
1528 if (nstripsN > 100) {
1529 printf("multiplicity > 100 - increase the dimension of the arrays %d\n",nstripsN);
1530 nstripsN=100;
1531 }
1532 for(int i=0;i<nstripsN;i++) {
1533 // check if 'clusterN->GetDigitStripNo(i)' returns the digit index
1534 cnew.fIndexMapN[i] = clusterN->GetDigitStripNo(i);
1535 }
1536 */
1537 cnew.fQErr=dSig;
1538 //cnew.fProbability=(float)prob;
1539 fITS->AddCluster(2,&cnew);
1540 fSegmentation->GetPadIxz(P,N,stripP,stripN);
1541 digP = (AliITSdigitSSD*)fMap->GetHit(1,stripP);
1542 digN = (AliITSdigitSSD*)fMap->GetHit(0,stripN);
1543 //printf("SSD: digP digN %p %p\n",digP,digN);
1544 if(digP) sigP = digP->fSignal;
1545 else sigP=0;
1546 if(digN) sigN = digN->fSignal;
1547 else sigN=0;
1548 if (!digP && !digN) {
1549 Error("CreateNewRecPoint","cannot find the digit!");
1550 return kFALSE;
1551 }
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);
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 //printf("sigP>sigN: %d %d %d\n",digP->fTracks[0],digP->fTracks[1],digP->fTracks[2]);
1566 } else {
1567 rnew.fTracks[0]=digN->fTracks[0];
1568 rnew.fTracks[1]=digN->fTracks[1];
1569 rnew.fTracks[2]=digN->fTracks[2];
1570 //printf("sigN>sigP: %d %d %d\n",digN->fTracks[0],digN->fTracks[1],digN->fTracks[2]);
1571 }
1572 //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);
1573 fITS->AddRecPoint(rnew);
1574 /*
1575 // it was
1576 fPointsM->AddLast( (TObject*)
1577 new AliITSRecPointSSD(P,dP,0,0,N,dN,Sig,dSig,fLayer,fLad,fDet,clusterP,clusterN,prob) );
1578 */
1579
1580 if(debug) cout<<"\n"<<": ImpactPoint Created: X = "<<P<<" Z = "<<N<<"; P = "<<p<<" N = "<<n<<"\n";
1581 return kTRUE;
1582 } else {
1583 if (debug) {
1584 cout<<"\n"<<": ATENTION : stupid ImpactPoit Point X = "<<P<<" Z = "<<N<<"; P = "<<p<<" N = "<<n<<"\n";
1585 }
1586 return kFALSE;
1587 }
1588}
1589
1590
1591
1592/**********************************************************************/
1593Bool_t AliITSClusterFinderSSD::
1594CreateNewRecPoint(AliITSclusterSSD *clusterP, AliITSclusterSSD *clusterN, Stat_t prob)
1595{
1596// create recpoints
1597
1598 Float_t posClusterP; //Cluster P position in strip coordinates
1599 Float_t posClusterN; //Cluster N position in strip coordinates
1600
1601 Float_t posErrorClusterP;
1602 Float_t posErrorClusterN;
1603
1604 Float_t sigClusterP;
1605 Float_t sigClusterN;
1606
1607 Float_t sigErrorClusterP;
1608 Float_t sigErrorClusterN;
1609
1610 posClusterP = clusterP->GetPosition();
1611 posErrorClusterP = clusterP->GetPositionError();
1612 sigClusterP = clusterP->GetTotalSignal();
1613 sigErrorClusterP = clusterP->GetTotalSignalError();
1614
1615 posClusterN = clusterN->GetPosition();
1616 posErrorClusterN = clusterN->GetPositionError();
1617 sigClusterN = clusterN->GetTotalSignal();
1618 sigErrorClusterN = clusterN->GetTotalSignalError();
1619
1620 return CreateNewRecPoint( posClusterP,posErrorClusterP, posClusterN,posErrorClusterN,
1621 sigClusterP+sigClusterN, sigErrorClusterN+sigErrorClusterP,
1622 clusterP, clusterN, prob);
1623
1624}
1625
1626//--------------------------------------------------
1627Bool_t AliITSClusterFinderSSD::IsCrossing(AliITSclusterSSD* p, AliITSclusterSSD* n)
1628{
1629// check for crossings
1630
1631 Float_t x = p->GetPosition();
1632 Float_t y = n->GetPosition();
1633 return GetCrossing(x,y);
1634}
1635
1636
1637//----------------------------------------------
1638Bool_t AliITSClusterFinderSSD::Strip2Local( Float_t stripP, Float_t stripN, Float_t &Z,Float_t &X)
1639{
1640// convert between strip numbers and local coordinates
1641
1642/*
1643 Z = (stripN-stripP)*fFactorOne;
1644 X = (stripN + stripP - fStripZeroOffset)*fFactorTwo;
1645*/
1646
1647 Float_t P =stripP;
1648 Float_t N =stripN;
1649
1650 GetCrossing(P,N);
1651 X=P;
1652 Z=N;
1653 if (debug) cout<<"P="<<stripP<<" N="<<stripN<<" X = "<<X<<" Z = "<<Z<<"\n";
1654 if ((Z<2.1)&&(Z>-2.1)&&(X<3.65)&&(X>-3.65)) return kTRUE;
1655 else return kFALSE;
1656
1657}
1658
1659
1660//-----------------------------------------------------------
1661Float_t AliITSClusterFinderSSD::GetClusterZ(AliITSclusterSSD* clust)
1662{
1663// get Z coordinate of the cluster
1664
1665 return clust->GetPosition();
1666
1667}
1668
1669//-------------------------------------------------------------
1670Int_t AliITSClusterFinderSSD::GetBestComb
1671(Int_t** comb,Int_t Ncomb, Int_t Ncl, AliITSpackageSSD * pkg)
1672{
1673//returns index of best combination in "comb"
1674//comb : sets of combinations
1675// in given combination on place "n" is index of
1676// Nside cluster coresponding to "n" P side cluster
1677//
1678//Ncomb : number of combinations == number of rows in "comb"
1679//
1680//Ncl : number of clusters in each row == number of columns in "comb"
1681//
1682//pkg : package
1683
1684
1685 Float_t currbval=-1; //current best value,
1686 //starting value is set to number negative,
1687 //to be changed in first loop turn automatically
1688
1689 Int_t curridx=0; //index of combination giving currently the best chi^2
1690 Float_t chi;
1691 Float_t ps, ns; //signal of P cluster and N cluster
1692
1693 for (Int_t i=0;i<Ncomb;i++)
1694 {
1695 chi=0;
1696
1697 for (Int_t j=0;j<Ncl;j++)
1698 {
1699 ps = pkg->GetPSideCluster(j)->GetTotalSignal(); //carrefully here, different functions
1700 ns = GetNSideCluster(comb[i][j])->GetTotalSignal();
1701 chi+=DistToPML(ps,ns);
1702 }
1703
1704 if (currbval==-1) currbval = chi;
1705 if (chi<currbval)
1706 {
1707 curridx = i;
1708 currbval =chi;
1709 }
1710 }
1711
1712
1713 return curridx;
1714
1715}
1716
1717//--------------------------------------------------
1718void AliITSClusterFinderSSD::GetBestMatchingPoint
1719(Int_t & ip, Int_t & in, AliITSpackageSSD* pkg )
1720{
1721// find the best matching point
1722
1723 if (debug) pkg->PrintClusters();
1724
1725 Float_t ps, ns; //signals on side p & n
1726
1727 Int_t nc; // number of crosses in given p side cluster
1728 Int_t n,p;
1729 AliITSclusterSSD *curPcl, *curNcl; //current p side cluster
1730 Float_t bestchi, chi;
1731 ip=p=pkg->GetPSideClusterIdx(0);
1732 in=n=pkg->GetNSideClusterIdx(0);
1733
1734 bestchi=DistToPML( pkg->GetPSideCluster(0)->GetTotalSignal(),
1735 pkg->GetNSideCluster(0)->GetTotalSignal() );
1736 for (Int_t i = 0; i< pkg->GetNumOfClustersP(); i++)
1737 {
1738
1739 p = pkg->GetPSideClusterIdx(i);
1740 curPcl= GetPSideCluster(p);
1741 nc=curPcl->GetCrossNo();
1742 ps=curPcl->GetTotalSignal();
1743
1744 for (Int_t j = 0; j< nc; j++)
1745 {
1746 n=curPcl->GetCross(j);
1747 curNcl= GetNSideCluster(n);
1748 ns=curNcl->GetTotalSignal();
1749 chi = DistToPML(ps, ns);
1750
1751 if (chi <= bestchi)
1752 {
1753 bestchi = chi;
1754 in = n;
1755 ip = p;
1756 }
1757 }
1758 }
1759
1760}
1761
1762
1763//------------------------------------------------------
1764void AliITSClusterFinderSSD::CalcStepFactor(Float_t Psteo, Float_t Nsteo )
1765{
1766// calculate the step factor for matching clusters
1767
1768
1769 // 95 is the pitch, 4000 - dimension along z ?
1770 /*
1771 fSFF = ( (Int_t) (Psteo*40000/95 ) );// +1;
1772 fSFB = ( (Int_t) (Nsteo*40000/95 ) );// +1;
1773 */
1774
1775
1776 Float_t dz=fSegmentation->Dz();
1777
1778 fSFF = ( (Int_t) (Psteo*dz/fPitch ) );// +1;
1779 fSFB = ( (Int_t) (Nsteo*dz/fPitch ) );// +1;
1780
1781}
1782
1783
1784//-----------------------------------------------------------
1785AliITSclusterSSD* AliITSClusterFinderSSD::GetPSideCluster(Int_t idx)
1786{
1787// get P side clusters
1788
1789
1790
1791
1792 if((idx<0)||(idx>=fNClusterP))
1793 {
1794 printf("AliITSClusterFinderSSD::GetPSideCluster : index out of range\n");
1795 return 0;
1796 }
1797 else
1798 {
1799 return (AliITSclusterSSD*)((*fClusterP)[idx]);
1800 }
1801}
1802
1803//-------------------------------------------------------
1804AliITSclusterSSD* AliITSClusterFinderSSD::GetNSideCluster(Int_t idx)
1805{
1806// get N side clusters
1807
1808 if((idx<0)||(idx>=fNClusterN))
1809 {
1810 printf("AliITSClusterFinderSSD::GetNSideCluster : index out of range\n");
1811 return 0;
1812 }
1813 else
1814 {
1815 return (AliITSclusterSSD*)((*fClusterN)[idx]);
1816 }
1817}
1818
1819//--------------------------------------------------------
1820AliITSclusterSSD* AliITSClusterFinderSSD::GetCluster(Int_t idx, Bool_t side)
1821{
1822// Get cluster
1823
1824 return (side) ? GetPSideCluster(idx) : GetNSideCluster(idx);
1825}
1826
1827
1828//--------------------------------------------------------
1829void AliITSClusterFinderSSD::ConsumeClusters()
1830{
1831// remove clusters from the list of unmatched clusters
1832
1833 for ( Int_t i=0;i<fNPackages;i++)
1834 {
1835 ((AliITSpackageSSD*)((*fPackages)[i]))->ConsumeClusters();
1836 }
1837
1838}
1839
1840
1841//--------------------------------------------------------
1842void AliITSClusterFinderSSD::ReconstructNotConsumedClusters()
1843{
1844// reconstruct remaining non-crossing clusters
1845
1846 //printf("ReconstructNotConsumedClusters !!!\n");
1847
1848 Int_t i;
1849 AliITSclusterSSD *cluster;
1850 Float_t pos;
1851 Float_t sig;
1852 Float_t sigerr;
1853 Float_t x1,x2,z1,z2;
1854
1855 Float_t Dx = fSegmentation->Dx();
1856 Float_t Dz = fSegmentation->Dz();
1857
1858 const Float_t kdEdXtoQ = 2.778e+8; // GeV -> number of e-hole pairs
1859 const Float_t kconv = 1.0e-4;
1860 const Float_t kRMSx = 20.0*kconv; // microns->cm ITS TDR Table 1.3
1861 const Float_t kRMSz = 830.0*kconv; // microns->cm ITS TDR Table 1.3
1862
1863 Int_t stripP,stripN;
1864 AliITSdigitSSD *dig;
1865 for (i=0;i<fNClusterP;i++)
1866 {
1867 //printf("P side cluster: cluster number %d\n",i);
1868 cluster = GetPSideCluster(i);
1869 if (!cluster->IsConsumed())
1870 {
1871 if( (pos = cluster->GetPosition()) < fSFB)
1872 {
1873 sig = cluster->GetTotalSignal();
1874
1875 sig += sig/fPNsignalRatio;
1876
1877 sigerr = cluster->GetTotalSignalError();
1878 x1 = -Dx/2 + pos *fPitch;
1879 z1 = -Dz/2;
1880
1881 x2 = pos;
1882 z2 = 1;
1883 GetCrossing (x2,z2);
1884 fSegmentation->GetPadIxz((x1+x2)/2,(z1+z2)/2,stripP,stripN);
1885 dig = (AliITSdigitSSD*)fMap->GetHit(1,stripP);
1886 //printf("dig %p\n",dig);
1887 if (!dig) printf("SSD: check the digit! dig %p\n",dig);
1888 if(!dig) continue;
1889 AliITSRawClusterSSD cnew;
1890 Int_t nstripsP=cluster->GetNumOfDigits();
1891 cnew.fMultiplicity=nstripsP;
1892 cnew.fMultiplicityN=0;
1893 /*
1894 if (nstripsP > 100) {
1895 printf("multiplicity > 100 - increase the dimension of the arrays %d\n",nstripsP);
1896 nstripsP=100;
1897 }
1898 for(int k=0;k<nstripsP;k++) {
1899 // check if 'clusterP->GetDigitStripNo(i)' returns
1900 // the digit index and not smth else
1901 cnew.fIndexMap[k] = cluster->GetDigitStripNo(k);
1902 }
1903 */
1904 cnew.fQErr=sigerr;
1905 //cnew.fProbability=0.75;
1906 //printf("AddCluster: %p nstripsP %d\n",&cnew,nstripsP);
1907 fITS->AddCluster(2,&cnew);
1908 // add the rec point info
1909 AliITSRecPoint rnew;
1910 rnew.SetX(kconv*(x1+x2)/2);
1911 rnew.SetZ(kconv*(z1+z2)/2);
1912 rnew.SetQ(sig);
1913 rnew.SetdEdX(sig/kdEdXtoQ);
1914 rnew.SetSigmaX2( kRMSx* kRMSx);
1915 rnew.SetSigmaZ2( kRMSz* kRMSz);
1916 //rnew.SetProbability(0.75);
1917 rnew.fTracks[0]=dig->fTracks[0];
1918 rnew.fTracks[1]=dig->fTracks[1];
1919 rnew.fTracks[2]=dig->fTracks[2];
1920 //printf("digP: %d %d %d\n",dig->fTracks[0],dig->fTracks[1],dig->fTracks[2]);
1921 //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);
1922 fITS->AddRecPoint(rnew);
1923 /*
1924 fPointsM->AddLast( (TObject*)
1925 new AliITSRecPointSSD((x1+x2)/2 ,0,0,0,(z1+z2)/2,0,sig,sigerr,fLayer,fLad,fDet,cluster, 0.75) );
1926 */
1927
1928 if(debug) cout<<"\n"<<": SINGLE SIDE ImpactPoint Created: X = "
1929 <<(x1+x2)/2<<" Z = "<<(z1+z2)/2<<"; Pos = "<<pos;
1930 }
1931 }
1932 }
1933
1934 for (i=0;i<fNClusterN;i++)
1935 {
1936 //printf("N side cluster: cluster number %d\n",i);
1937 cluster = GetNSideCluster(i);
1938 if (!cluster->IsConsumed())
1939 {
1940 if( (pos = cluster->GetPosition()) < fSFF)
1941 {
1942 sig = cluster->GetTotalSignal();
1943
1944 sig += sig*fPNsignalRatio;
1945
1946 sigerr = cluster->GetTotalSignalError();
1947
1948 x1 = -Dx/2 + pos *fPitch;
1949 z1 = Dz/2;
1950
1951 x2 = 1;
1952 z2 = pos;
1953
1954 GetCrossing (x2,z2);
1955 fSegmentation->GetPadIxz((x1+x2)/2,(z1+z2)/2,stripP,stripN);
1956 dig = (AliITSdigitSSD*)fMap->GetHit(0,stripN);
1957 if (!dig) printf("SSD: check the digit! dig %p\n",dig);
1958 if(!dig) continue;
1959 AliITSRawClusterSSD cnew;
1960 Int_t nstripsN=cluster->GetNumOfDigits();
1961 cnew.fMultiplicity=0;
1962 cnew.fMultiplicityN=nstripsN;
1963 /*
1964 if (nstripsN > 100) {
1965 printf("multiplicity > 100 - increase the dimension of the arrays %d\n",nstripsN);
1966 nstripsN=100;
1967 }
1968 for(int k=0;k<nstripsN;k++) {
1969 // check if 'clusterP->GetDigitStripNo(i)' returns
1970 // the digit index and not smth else
1971 cnew.fIndexMapN[k] = cluster->GetDigitStripNo(k);
1972 }
1973 */
1974 cnew.fQErr=sigerr;
1975 //cnew.fProbability=0.75;
1976 fITS->AddCluster(2,&cnew);
1977 AliITSRecPoint rnew;
1978 rnew.SetX(kconv*(x1+x2)/2);
1979 rnew.SetZ(kconv*(z1+z2)/2);
1980 rnew.SetQ(sig);
1981 rnew.SetdEdX(sig/kdEdXtoQ);
1982 rnew.SetSigmaX2( kRMSx* kRMSx);
1983 rnew.SetSigmaZ2( kRMSz* kRMSz);
1984 //rnew.SetProbability(0.75);
1985 rnew.fTracks[0]=dig->fTracks[0];
1986 rnew.fTracks[1]=dig->fTracks[1];
1987 rnew.fTracks[2]=dig->fTracks[2];
1988 //printf("digN: %d %d %d\n",dig->fTracks[0],dig->fTracks[1],dig->fTracks[2]);
1989 //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);
1990 fITS->AddRecPoint(rnew);
1991 /*
1992 fPointsM->AddLast( (TObject*)
1993 new AliITSRecPointSSD((x1+x2)/2 ,0,0,0,(z1+z2)/2,0,sig,sigerr,fLayer,fLad,fDet,cluster, 0.75) );
1994 */
1995
1996 if(debug) cout<<"\n"<<": SINGLE SIDE ImpactPoint Created: X = "
1997 <<(x1+x2)/2<<" Z = "<<(z1+z2)/2<<"; Pos = "<<pos;
1998
1999 }
2000 }
2001 }
2002
2003
2004}
2005
2006//_______________________________________________________________________
2007
2008Bool_t AliITSClusterFinderSSD::GetCrossing (Float_t &P, Float_t &N)
2009{
2010// get crossing
2011
2012 Float_t Dx = fSegmentation->Dx();
2013 Float_t Dz = fSegmentation->Dz();
2014
2015 //P==X
2016 //N==Z
2017
2018 Float_t dx = 0.1;
2019 P *= fPitch;
2020 N *= fPitch;
2021
2022 //P = (kZ * fTan + N + P)/2.0; // x coordinate
2023 //N = kZ - (P-N)/fTan; // z coordinate
2024
2025 if(fTanP + fTanN == 0) return kFALSE;
2026
2027 N = (N - P + fTanN * Dz) / (fTanP + fTanN); // X coordinate
2028 P = P + fTanP * N; // Y coordinate
2029
2030 P -= Dx/2;
2031 N -= Dz/2;
2032
2033 // N = -(N - P + kZ/2*(fTanP + fTanN))/(fTanP + fTanN);
2034 // P = (fTanP*(N-kX/2-kZ*fTanN/2) + fTanN*(P-kX/2-kZ*fTanP/2) )/(fTanP + fTanN);
2035
2036 if ( ( N < -(Dz/2+dx) ) || ( N > (Dz/2+dx) ) ) return kFALSE;
2037 if ( ( P < -(Dx/2+dx) ) || ( P > (Dx/2+dx) ) ) return kFALSE;
2038
2039 return kTRUE;
2040}
2041
2042
2043//_________________________________________________________________________
2044
2045void AliITSClusterFinderSSD::GetCrossingError(Float_t& dP, Float_t& dN)
2046{
2047// get crossing error
2048
2049 Float_t dz, dx;
2050
2051 dz = TMath::Abs(( dP + dN )*fPitch/(fTanP + fTanN) );
2052 dx = fPitch*(TMath::Abs(dP*(1 - fTanP/(fTanP + fTanN))) +
2053 TMath::Abs(dN *fTanP/(fTanP + fTanN) ));
2054
2055 dN = dz;
2056 dP = dx;
2057}