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