Line data Source code
1 : /**************************************************************************
2 : * Copyright(c) 2007-2009, 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 : /* $Id$ */
17 :
18 : ////////////////////////////////////////////////////////////////////////////
19 : // Implementation of the ITS clusterer V2 class //
20 : // //
21 : // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch //
22 : // Last revision: 13-05-09 Enrico Fragiacomo //
23 : // enrico.fragiacomo@ts.infn.it //
24 : // //
25 : ///////////////////////////////////////////////////////////////////////////
26 :
27 : #include "AliITSClusterFinderV2SSD.h"
28 :
29 : #include <Riostream.h>
30 : #include <TGeoGlobalMagField.h>
31 :
32 : #include "AliLog.h"
33 : #include "AliMagF.h"
34 : #include "AliITSRecPoint.h"
35 : #include "AliITSRecPointContainer.h"
36 : #include "AliITSgeomTGeo.h"
37 : #include "AliRawReader.h"
38 : #include "AliITSRawStreamSSD.h"
39 : #include <TClonesArray.h>
40 : #include <TCollection.h>
41 : #include "AliITSdigitSSD.h"
42 : #include "AliITSReconstructor.h"
43 : #include "AliITSCalibrationSSD.h"
44 : #include "AliITSsegmentationSSD.h"
45 :
46 : Short_t *AliITSClusterFinderV2SSD::fgPairs = 0x0;
47 : Int_t AliITSClusterFinderV2SSD::fgPairsSize = 0;
48 : const Float_t AliITSClusterFinderV2SSD::fgkThreshold = 5.;
49 :
50 : const Float_t AliITSClusterFinderV2SSD::fgkCosmic2008StripShifts[16][9] =
51 : {{-0.35,-0.35,-0.35,-0.35,-0.35,-0.35,-0.35,-0.35,-0.35}, // DDL 512
52 : {-0.35,-0.35,-0.35,-0.35,-0.35,-0.35,-0.35,-0.35,-0.35}, // DDL 513
53 : {-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15}, // DDL 514
54 : {-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15}, // DDL 515
55 : { 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00}, // DDL 516
56 : { 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00}, // DDL 517
57 : {-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15}, // DDL 518
58 : {-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15}, // DDL 519
59 : {-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.25,-0.15}, // DDL 520
60 : {-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15}, // DDL 521
61 : {-0.10,-0.10,-0.10,-0.40,-0.40,-0.40,-0.10,-0.10,-0.45}, // DDL 522
62 : {-0.10,-0.10,-0.10,-0.35,-0.35,-0.35,-0.10,-0.35,-0.50}, // DDL 523
63 : { 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00}, // DDL 524
64 : { 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00}, // DDL 525
65 : { 0.35, 0.35, 0.35, 0.35, 0.35, 0.35, 0.35, 0.35, 0.35}, // DDL 526
66 : { 0.45, 0.45, 0.45, 0.45, 0.45, 0.45, 0.45, 0.45, 0.45}}; // DDL 527
67 :
68 118 : ClassImp(AliITSClusterFinderV2SSD)
69 :
70 :
71 14 : AliITSClusterFinderV2SSD::AliITSClusterFinderV2SSD(AliITSDetTypeRec* dettyp):AliITSClusterFinder(dettyp),fLastSSD1(AliITSgeomTGeo::GetModuleIndex(6,1,1)-1), fLorentzShiftP(0), fLorentzShiftN(0)
72 10 : {
73 : //Default constructor
74 : static AliITSRecoParam *repa = NULL;
75 2 : if(!repa){
76 4 : repa = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
77 2 : if(!repa){
78 0 : repa = AliITSRecoParam::GetHighFluxParam();
79 0 : AliWarning("Using default AliITSRecoParam class");
80 : }
81 : }
82 :
83 2 : if (repa->GetCorrectLorentzAngleSSD()) {
84 8 : AliMagF* field = dynamic_cast<AliMagF*>(TGeoGlobalMagField::Instance()->GetField());
85 2 : if (field == 0) {
86 0 : AliError("Cannot get magnetic field from TGeoGlobalMagField");
87 : }
88 : else {
89 2 : Float_t bField = field->SolenoidField();
90 : // NB: spatial shift has opposite sign for lay 5 and 6, but strip numbering also changes direction, so no sign-change
91 : // Shift due to ExB on drift N-side, units: strip width
92 2 : fLorentzShiftP = -repa->GetTanLorentzAngleElectronsSSD() * 150.e-4/95.e-4 * bField / 5.0;
93 : // Shift due to ExB on drift P-side, units: strip width
94 2 : fLorentzShiftN = -repa->GetTanLorentzAngleHolesSSD() * 150.e-4/95.e-4 * bField / 5.0;
95 10 : AliDebug(1,Form("bField %f Lorentz Shift P-side %f N-side %f",bField,fLorentzShiftN,fLorentzShiftP));
96 : }
97 2 : }
98 4 : }
99 :
100 : //______________________________________________________________________
101 0 : AliITSClusterFinderV2SSD::AliITSClusterFinderV2SSD(const AliITSClusterFinderV2SSD &cf) : AliITSClusterFinder(cf), fLastSSD1(cf.fLastSSD1), fLorentzShiftP(cf.fLorentzShiftP), fLorentzShiftN(cf.fLorentzShiftN)
102 0 : {
103 : // Copy constructor
104 0 : }
105 :
106 : //______________________________________________________________________
107 : AliITSClusterFinderV2SSD& AliITSClusterFinderV2SSD::operator=(const AliITSClusterFinderV2SSD& cf ){
108 : // Assignment operator
109 :
110 0 : this->~AliITSClusterFinderV2SSD();
111 0 : new(this) AliITSClusterFinderV2SSD(cf);
112 0 : return *this;
113 0 : }
114 :
115 :
116 : void AliITSClusterFinderV2SSD::FindRawClusters(Int_t mod){
117 :
118 : //Find clusters V2
119 11048 : SetModule(mod);
120 5524 : FindClustersSSD(fDigits);
121 :
122 5524 : }
123 :
124 : void AliITSClusterFinderV2SSD::FindClustersSSD(TClonesArray *alldigits) {
125 : //------------------------------------------------------------
126 : // Actual SSD cluster finder
127 : //------------------------------------------------------------
128 11048 : Int_t smaxall=alldigits->GetEntriesFast();
129 5524 : if (smaxall==0) return;
130 :
131 :
132 : //---------------------------------------
133 : // load recoparam and calibration
134 : //
135 : static AliITSRecoParam *repa = NULL;
136 5524 : if(!repa){
137 1 : repa = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
138 1 : if(!repa){
139 0 : repa = AliITSRecoParam::GetHighFluxParam();
140 0 : AliWarning("Using default AliITSRecoParam class");
141 0 : }
142 : }
143 :
144 5524 : AliITSCalibrationSSD* cal = (AliITSCalibrationSSD*)GetResp(fModule);
145 : Float_t gain=0;
146 : Float_t noise=0;
147 : //---------------------------------------
148 :
149 :
150 : //------------------------------------
151 : // fill the digits array with zero-suppression condition
152 : // Signal is converted in KeV
153 : //
154 5524 : TObjArray digits;
155 35978 : for (Int_t i=0;i<smaxall; i++){
156 12465 : AliITSdigitSSD *d=(AliITSdigitSSD*)alldigits->UncheckedAt(i);
157 :
158 23841 : if(d->IsSideP()) noise = cal->GetNoiseP(d->GetStripNumber());
159 6777 : else noise = cal->GetNoiseN(d->GetStripNumber());
160 26620 : if (d->GetSignal()<3.*noise) continue;
161 :
162 20403 : if(d->IsSideP()) gain = cal->GetGainP(d->GetStripNumber());
163 5961 : else gain = cal->GetGainN(d->GetStripNumber());
164 :
165 21550 : Float_t q=gain*d->GetSignal(); //
166 10775 : q=cal->ADCToKeV(q); // converts the charge in KeV from ADC units
167 10775 : d->SetSignal(Int_t(q));
168 :
169 10775 : digits.AddLast(d);
170 10775 : }
171 5524 : Int_t smax = digits.GetEntriesFast();
172 5845 : if (smax==0) return;
173 : //------------------------------------
174 :
175 :
176 : const Int_t kMax=1000;
177 5203 : Int_t np=0, nn=0;
178 5203 : Ali1Dcluster pos[kMax], neg[kMax];
179 : Float_t y=0., q=0., qmax=0.;
180 5203 : Int_t lab[4]={-2,-2,-2,-2};
181 : Bool_t flag5 = 0;
182 :
183 : /*
184 : cout<<"-----------------------------"<<endl;
185 : cout<<"this is module "<<fModule;
186 : cout<<endl;
187 : cout<<endl;
188 :
189 : Int_t layer = 4;
190 : if (fModule>fLastSSD1)
191 : layer = 5;
192 : */
193 : //--------------------------------------------------------
194 : // start 1D-clustering from the first digit in the digits array
195 : //
196 5203 : AliITSdigitSSD *d=(AliITSdigitSSD*)digits.UncheckedAt(0);
197 10406 : q += d->GetSignal();
198 15609 : y += d->GetCoord2()*d->GetSignal();
199 10406 : qmax=d->GetSignal();
200 20812 : lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2);
201 :
202 10406 : if(d->IsSideP()) {
203 8532 : noise = cal->GetNoiseP(d->GetStripNumber());
204 3329 : gain = cal->GetGainP(d->GetStripNumber());
205 3329 : }
206 : else {
207 1874 : noise = cal->GetNoiseN(d->GetStripNumber());
208 1874 : gain = cal->GetGainN(d->GetStripNumber());
209 : }
210 5203 : noise*=gain;
211 5203 : noise=cal->ADCToKeV(noise); // converts noise in KeV from ADC units
212 :
213 5311 : if(qmax>fgkThreshold*noise) flag5=1; // seed for the cluster
214 :
215 : /*
216 : cout<<d->GetSignal()<<" "<<noise<<" "<<flag5<<" "<<
217 : d->GetCoord1()<<" "<<d->GetCoord2()<<endl;
218 : */
219 :
220 5203 : Int_t curr=d->GetCoord2();
221 5203 : Int_t flag=d->GetCoord1();
222 :
223 : // Note: the first side which will be processed is supposed to be the
224 : // P-side which is neg
225 : Int_t *n=&nn;
226 5203 : Ali1Dcluster *c=neg;
227 7077 : if(flag) {n=&np; c=pos;} // in case we have only Nstrips (P was bad!)
228 :
229 : Int_t nd=1;
230 5203 : Int_t milab[10];
231 114466 : for (Int_t ilab=0;ilab<10;ilab++){
232 52030 : milab[ilab]=-2;
233 : }
234 20812 : milab[0]=d->GetTrack(0); milab[1]=d->GetTrack(1); milab[2]=d->GetTrack(2);
235 :
236 :
237 : //----------------------------------------------------------
238 : // search for neighboring digits
239 : //
240 26753 : for (Int_t s=1; s<smax; s++) {
241 5572 : d=(AliITSdigitSSD*)digits.UncheckedAt(s);
242 5572 : Int_t strip=d->GetCoord2();
243 :
244 : // if digits is not a neighbour or side did not change
245 : // and at least one of the previous digits met the seed condition
246 : // then creates a new 1D cluster
247 8566 : if ( ( ((strip-curr) > 1) || (flag!=d->GetCoord1()) ) ) {
248 :
249 5379 : if(flag5) {
250 : //cout<<"here1"<<endl;
251 : Float_t dLorentz = 0;
252 220 : if (!flag) { // P-side is neg clust
253 155 : dLorentz = fLorentzShiftN;
254 155 : }
255 : else { // N-side is p clust
256 65 : dLorentz = fLorentzShiftP;
257 : }
258 220 : c[*n].SetY(y/q+dLorentz);
259 220 : c[*n].SetQ(q);
260 220 : c[*n].SetNd(nd);
261 220 : CheckLabels2(milab);
262 220 : c[*n].SetLabels(milab);
263 :
264 220 : if(repa->GetUseUnfoldingInClusterFinderSSD()==kTRUE) {
265 : // Note: fUseUnfoldingInClusterFinderSSD=kFALSE by default in RecoParam
266 :
267 : //Split suspiciously big cluster
268 0 : if (nd>4&&nd<25) {
269 0 : c[*n].SetY(y/q-0.25*nd+dLorentz);
270 0 : c[*n].SetQ(0.5*q);
271 0 : (*n)++;
272 0 : if (*n==kMax) {
273 0 : Error("FindClustersSSD","Too many 1D clusters !");
274 0 : return;
275 : }
276 0 : c[*n].SetY(y/q+0.25*nd+dLorentz);
277 0 : c[*n].SetQ(0.5*q);
278 0 : c[*n].SetNd(nd);
279 0 : c[*n].SetLabels(milab);
280 0 : }
281 :
282 : } // unfolding is on
283 :
284 220 : (*n)++;
285 220 : if (*n==kMax) {
286 0 : Error("FindClustersSSD","Too many 1D clusters !");
287 0 : return;
288 : }
289 :
290 220 : } // flag5 set
291 :
292 : // reset everything
293 : y=q=qmax=0.;
294 : nd=0;
295 : flag5=0;
296 5379 : lab[0]=lab[1]=lab[2]=-2;
297 118338 : for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
298 :
299 : // if side changed from P to N, switch to pos 1D clusters
300 : // (if for some reason the side changed from N to P then do the opposite)
301 10758 : if (flag!=d->GetCoord1())
302 3954 : { if(!flag) {n=&np; c=pos;} else {n=&nn; c=neg;} }
303 :
304 : } // end create new 1D cluster from previous neighboring digits
305 :
306 : // continues adding digits to the previous cluster
307 : // or start a new one
308 5572 : flag=d->GetCoord1();
309 11144 : q += d->GetSignal();
310 16716 : y += d->GetCoord2()*d->GetSignal();
311 5572 : nd++;
312 :
313 11144 : if(d->IsSideP()) {
314 7057 : noise = cal->GetNoiseP(d->GetStripNumber());
315 1485 : gain = cal->GetGainP(d->GetStripNumber());
316 1485 : }
317 : else {
318 4087 : noise = cal->GetNoiseN(d->GetStripNumber());
319 4087 : gain = cal->GetGainN(d->GetStripNumber());
320 : }
321 5572 : noise*=gain;
322 5572 : noise=cal->ADCToKeV(noise); // converts the charge in KeV from ADC units
323 :
324 11491 : if(d->GetSignal()>fgkThreshold*noise) flag5=1;
325 :
326 : /*
327 : cout<<d->GetSignal()<<" "<<noise<<" "<<flag5<<" "<<
328 : d->GetCoord1()<<" "<<d->GetCoord2()<<endl;
329 : */
330 :
331 11144 : if (d->GetSignal()>qmax) {
332 10920 : qmax=d->GetSignal();
333 21840 : lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2);
334 5460 : }
335 122584 : for (Int_t ilab=0;ilab<10;ilab++) {
336 112374 : if (d->GetTrack(ilab)>=0) AddLabel(milab, (d->GetTrack(ilab)));
337 : }
338 : curr=strip;
339 :
340 :
341 5572 : } // loop over digits, no more digits in the digits array
342 :
343 :
344 : // add the last 1D cluster
345 5203 : if(flag5) {
346 :
347 : // cout<<"here2"<<endl;
348 : Float_t dLorentz = 0;
349 95 : if (!flag) { // P-side is neg clust
350 4 : dLorentz = fLorentzShiftN;
351 4 : }
352 : else { // N-side is p clust
353 91 : dLorentz = fLorentzShiftP;
354 : }
355 :
356 95 : c[*n].SetY(y/q + dLorentz);
357 95 : c[*n].SetQ(q);
358 95 : c[*n].SetNd(nd);
359 95 : c[*n].SetLabels(lab);
360 :
361 95 : if(repa->GetUseUnfoldingInClusterFinderSSD()==kTRUE) {
362 :
363 : //Split suspiciously big cluster
364 0 : if (nd>4 && nd<25) {
365 0 : c[*n].SetY(y/q-0.25*nd + dLorentz);
366 0 : c[*n].SetQ(0.5*q);
367 0 : (*n)++;
368 0 : if (*n==kMax) {
369 0 : Error("FindClustersSSD","Too many 1D clusters !");
370 0 : return;
371 : }
372 0 : c[*n].SetY(y/q+0.25*nd + dLorentz);
373 0 : c[*n].SetQ(0.5*q);
374 0 : c[*n].SetNd(nd);
375 0 : c[*n].SetLabels(lab);
376 0 : }
377 : } // unfolding is on
378 :
379 95 : (*n)++;
380 95 : if (*n==kMax) {
381 0 : Error("FindClustersSSD","Too many 1D clusters !");
382 0 : return;
383 : }
384 :
385 95 : } // if flag5 last 1D cluster added
386 :
387 :
388 : //------------------------------------------------------
389 : // call FindClustersSSD to pair neg and pos 1D clusters
390 : // and create recpoints from the crosses
391 : // Note1: neg are Pside and pos are Nside!!
392 : // Note2: if there are no Pside digits nn=0 (bad strips!!) (same for Nside)
393 : //
394 : // cout<<nn<<" Pside and "<<np<<" Nside clusters"<<endl;
395 :
396 5203 : AliITSRecPointContainer* rpc = AliITSRecPointContainer::Instance();
397 5203 : if (nn*np > 0) {
398 137 : TClonesArray* clusters = rpc->UncheckedGetClusters(fModule);
399 137 : clusters->Clear();
400 137 : FindClustersSSD(neg, nn, pos, np, clusters);
401 137 : TIter itr(clusters);
402 : AliITSRecPoint *irp;
403 849 : while ((irp = (AliITSRecPoint*)itr.Next())) fDetTypeRec->AddRecPoint(*irp);
404 137 : }
405 : //-----------------------------------------------------
406 21454 : }
407 :
408 :
409 : void AliITSClusterFinderV2SSD::RawdataToClusters(AliRawReader* rawReader){
410 :
411 : //------------------------------------------------------------
412 : // This function creates ITS clusters from raw data
413 : //------------------------------------------------------------
414 8 : fNClusters = 0;
415 4 : rawReader->Reset();
416 4 : AliITSRawStreamSSD inputSSD(rawReader);
417 4 : FindClustersSSD(&inputSSD);
418 :
419 4 : }
420 :
421 :
422 : void AliITSClusterFinderV2SSD::FindClustersSSD(AliITSRawStreamSSD* input)
423 : {
424 : //------------------------------------------------------------
425 : // Actual SSD cluster finder for raw data
426 : //------------------------------------------------------------
427 :
428 8 : AliITSRecPointContainer* rpc = AliITSRecPointContainer::Instance();
429 : static AliITSRecoParam *repa = NULL;
430 4 : if(!repa){
431 1 : repa = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
432 1 : if(!repa){
433 0 : repa = AliITSRecoParam::GetHighFluxParam();
434 0 : AliWarning("Using default AliITSRecoParam class");
435 0 : }
436 : }
437 4 : if (fRawID2ClusID) { // RS: reset references from 1D clusters to rawID's
438 0 : fRawIDRef[0].Reset();
439 0 : fRawIDRef[1].Reset();
440 0 : }
441 : Int_t nClustersSSD = 0;
442 : const Int_t kNADC = 12;
443 : const Int_t kMaxADCClusters = 1000;
444 :
445 4 : Int_t strips[kNADC][2][kMaxADCClusters][3]; // [ADC],[side],[istrip], [0]=istrip [1]=signal [2]=rawID (for embedding, RS)
446 4 : Int_t nStrips[kNADC][2];
447 :
448 104 : for( int i=0; i<kNADC; i++ ){
449 48 : nStrips[i][0] = 0;
450 48 : nStrips[i][1] = 0;
451 : }
452 :
453 : Int_t ddl = -1;
454 : Int_t ad = -1;
455 :
456 : //*
457 : //* Loop over modules DDL+AD
458 : //*
459 : int countRW = 0; //RS
460 4 : if (fRawID2ClusID) fRawID2ClusID->Reset(); //RS if array was provided, we shall store the rawID -> ClusterID
461 :
462 : while (kTRUE) {
463 :
464 12469 : bool next = input->Next();
465 :
466 : //*
467 : //* Continue if corrupted input
468 : //*
469 :
470 12473 : if( (!next)&&(input->flag) ){
471 0 : AliWarning(Form("ITSClustersFinderSSD: Corrupted data: warning from RawReader"));
472 0 : continue;
473 : }
474 :
475 12469 : Int_t newDDL = input->GetDDL();
476 12469 : Int_t newAD = input->GetAD();
477 :
478 12469 : if( next ){
479 12465 : if( newDDL<0 || newDDL>15 ){
480 0 : AliWarning(Form("ITSClustersFinderSSD: Corrupted data: wrong DDL number (%d)",newDDL));
481 0 : continue;
482 : }
483 :
484 12465 : if( newAD<1 || newAD>9 ){
485 0 : AliWarning(Form("ITSClustersFinderSSD: Corrupted data: wrong AD number (%d)",newAD));
486 0 : continue;
487 : }
488 : }
489 :
490 49804 : bool newModule = ( !next || ddl!= newDDL || ad!=newAD );
491 :
492 12469 : if( newModule && ddl>=0 && ad>=0 ){
493 :
494 : //*
495 : //* Reconstruct the previous block of 12 modules --- actual clusterfinder
496 : //*
497 : //cout<<endl;
498 14690 : for( int adc = 0; adc<kNADC; adc++ ){
499 :
500 : //* 1D clusterfinder
501 :
502 6780 : Ali1Dcluster clusters1D[2][kMaxADCClusters]; // per ADC, per side
503 6780 : Int_t nClusters1D[2] = {0,0};
504 : //int nstat[2] = {0,0};
505 6780 : fModule = AliITSRawStreamSSD::GetModuleNumber(ddl, (ad - 1) * 12 + adc );
506 :
507 6780 : if( fModule<0 ){
508 : // AliWarning(Form("ITSClustersFinderSSD: Corrupted data: module (ddl %d ad %d adc %d) not found in the map",ddl,ad,adc));
509 : //CM channels are always present even everything is suppressed
510 117 : continue;
511 : }
512 :
513 : /* Int_t layer = 4;
514 : if (fModule>fLastSSD1)
515 : layer = 5;
516 : */
517 :
518 6663 : AliITSCalibrationSSD* cal = (AliITSCalibrationSSD*)fDetTypeRec->GetCalibrationModel(fModule);
519 6663 : if( !cal ){
520 0 : AliWarning(Form("ITSClustersFinderSSD: No calibration found for module (ddl %d ad %d adc %d)",ddl,ad,adc));
521 0 : continue;
522 : }
523 :
524 : Float_t dStrip = 0;
525 :
526 6663 : if( repa->GetUseCosmicRunShiftsSSD()) { // Special condition for 2007/2008 cosmic data
527 0 : dStrip = fgkCosmic2008StripShifts[ddl][ad-1];
528 0 : if (TMath::Abs(dStrip) > 1.5){
529 0 : AliWarning(Form("Indexing error in Cosmic calibration: ddl = %d, dStrip %f\n",ddl,dStrip));
530 : dStrip = 0;
531 0 : }
532 : }
533 :
534 39978 : for( int side=0; side<=1; side++ ){
535 :
536 13326 : Int_t lab[3]={-2,-2,-2};
537 : Float_t q = 0.;
538 : Float_t y = 0.;
539 : Int_t nDigits = 0;
540 : Int_t ostrip = -2;
541 : Bool_t snFlag = 0;
542 :
543 : Float_t dLorentz = 0;
544 13326 : if (side==0) { // P-side is neg clust
545 6663 : dLorentz = fLorentzShiftN;
546 6663 : }
547 : else { // N-side is pos clust
548 6663 : dLorentz = fLorentzShiftP;
549 : }
550 :
551 13326 : Int_t n = nStrips[adc][side];
552 78234 : for( int istr = 0; istr<n+1; istr++ ){
553 :
554 : bool stripOK = 1;
555 : Int_t strip=0, rwID = 0;
556 : Float_t signal=0.0, noise=0.0, gain=0.0;
557 :
558 25791 : if( istr<n ){
559 12465 : strip = strips[adc][side][istr][0];
560 12465 : signal = strips[adc][side][istr][1];
561 12465 : rwID = strips[adc][side][istr][2]; // RS
562 : //cout<<"strip "<<adc<<" / "<<side<<": "<<strip<<endl;
563 :
564 12465 : if( cal ){
565 37395 : noise = side ?cal->GetNoiseN(strip) :cal->GetNoiseP(strip);
566 37395 : gain = side ?cal->GetGainN(strip) :cal->GetGainP(strip);
567 37349 : stripOK = ( noise>=1. && signal>=3.0*noise
568 : //&& !cal->IsPChannelBad(strip)
569 : );
570 12465 : }
571 : } else stripOK = 0; // end of data
572 :
573 51780 : bool newCluster = ( (abs(strip-ostrip)!=1) || !stripOK );
574 :
575 25791 : if( newCluster ){
576 :
577 : //* Store the previous cluster
578 :
579 36135 : if( nDigits>0 && q>0 && snFlag ){
580 :
581 317 : if (nClusters1D[side] >= kMaxADCClusters-1 ) {
582 0 : AliWarning("HLT ClustersFinderSSD: Too many 1D clusters !");
583 0 : }else {
584 :
585 317 : Ali1Dcluster &cluster = clusters1D[side][nClusters1D[side]++];
586 317 : cluster.SetY( y / q + dStrip + dLorentz);
587 317 : cluster.SetQ(q);
588 317 : cluster.SetNd(nDigits);
589 317 : cluster.SetLabels(lab);
590 : //cout<<"cluster 1D side "<<side<<": y= "<<y<<" q= "<<q<<" d="<<dStrip<<" Y="<<cluster.GetY()<<endl;
591 : //Split suspiciously big cluster
592 :
593 317 : if( repa->GetUseUnfoldingInClusterFinderSSD()
594 317 : && nDigits > 4 && nDigits < 25
595 : ){
596 0 : cluster.SetY(y/q + dStrip - 0.25*nDigits + dLorentz);
597 0 : cluster.SetQ(0.5*q);
598 0 : Ali1Dcluster& cluster2 = clusters1D[side][nClusters1D[side]++];
599 0 : cluster2.SetY(y/q + dStrip + 0.25*nDigits + dLorentz);
600 0 : cluster2.SetQ(0.5*q);
601 0 : cluster2.SetNd(nDigits);
602 0 : cluster2.SetLabels(lab);
603 0 : } // unfolding is on
604 : }
605 : }
606 : y = q = 0.;
607 : nDigits = 0;
608 : snFlag = 0;
609 :
610 25598 : } //* End store the previous cluster
611 :
612 25791 : if( stripOK ){ // add new signal to the cluster
613 : // signal = (Int_t) (signal * gain); // signal is corrected for gain
614 11191 : if( signal>fgkThreshold*noise) snFlag = 1;
615 10730 : signal = signal * gain; // signal is corrected for gain
616 : // if( cal ) signal = (Int_t) cal->ADCToKeV( signal ); // signal is converted in KeV
617 21460 : if( cal ) signal = cal->ADCToKeV( signal ); // signal is converted in KeV
618 10730 : q += signal; // add digit to current cluster
619 10730 : y += strip * signal;
620 10730 : nDigits++;
621 : //nstat[side]++;
622 : ostrip = strip;
623 10730 : if (fRawID2ClusID) fRawIDRef[side].AddReference(nClusters1D[side],rwID);
624 :
625 : }
626 : } //* end loop over strips
627 :
628 13326 : } //* end loop over ADC sides
629 :
630 :
631 : //* 2D clusterfinder
632 6949 : if( nClusters1D[0] && nClusters1D[1] && fModule>=0 ){
633 137 : TClonesArray* clusters = rpc->UncheckedGetClusters(fModule);
634 137 : FindClustersSSD( clusters1D[0], nClusters1D[0], clusters1D[1], nClusters1D[1], clusters);
635 137 : Int_t nClustersn = clusters->GetEntriesFast();
636 137 : nClustersSSD += nClustersn;
637 137 : }
638 :
639 : //cout<<"SG: "<<ddl<<" "<<ad<<" "<<adc<<": strips "<<nstat[0]<<"+"<<nstat[1]<<", clusters 1D= "<<nClusters1D[0]<<" + "<<nClusters1D[1]<<", 2D= "<<clusters.size()<<endl;
640 :
641 13443 : }//* end loop over adc
642 :
643 565 : }//* end of reconstruction of previous block of 12 modules
644 :
645 12469 : if( newModule ){
646 :
647 : //*
648 : //* Clean up arrays and set new module
649 : //*
650 :
651 14794 : for( int i=0; i<kNADC; i++ ){
652 6828 : nStrips[i][0] = 0;
653 6828 : nStrips[i][1] = 0;
654 : }
655 : ddl = newDDL;
656 : ad = newAD;
657 569 : }
658 :
659 :
660 : //*
661 : //* Exit main loop when there is no more input
662 : //*
663 :
664 12473 : if( !next ) break;
665 :
666 : //*
667 : //* Fill the current strip information
668 : //*
669 :
670 12465 : Int_t adc = input->GetADC();
671 24930 : if( adc<0 || adc>=kNADC+2 || (adc>5&&adc<8) ){
672 0 : AliWarning(Form("HLT ClustersFinderSSD: Corrupted data: wrong adc number (%d)", adc));
673 0 : continue;
674 : }
675 :
676 18808 : if( adc>7 ) adc-= 2; // shift ADC numbers 8-13 to 6-11
677 :
678 12465 : Bool_t side = input->GetSideFlag();
679 12465 : Int_t strip = input->GetStrip();
680 12465 : Int_t signal = input->GetSignal();
681 :
682 :
683 : //cout<<"SSD: "<<ddl<<" "<<ad<<" "<<adc<<" "<<side<<" "<<strip<<" : "<<signal<<endl;
684 :
685 12465 : if( strip>767 ){
686 0 : AliWarning(Form("HLT ClustersFinderSSD: Corrupted data: wrong strip number (ddl %d ad %d adc %d side %d, strip %d",
687 : ddl, ad, adc, side,strip) );
688 0 : continue;
689 : }
690 12465 : if (strip < 0) continue;
691 :
692 12465 : int &n = nStrips[adc][side];
693 12465 : if( n >0 ){
694 4538 : Int_t oldStrip = strips[adc][side][n-1][0];
695 :
696 4538 : if( strip==oldStrip ){
697 0 : AliWarning(Form("HLT ClustersFinderSSD: Corrupted data: duplicated signal: ddl %d ad %d adc %d, side %d, strip %d",
698 : ddl, ad, adc, side, strip ));
699 0 : continue;
700 : }
701 4538 : }
702 12465 : strips[adc][side][n][0] = strip;
703 12465 : strips[adc][side][n][1] = signal;
704 12465 : strips[adc][side][n][2] = countRW;
705 12465 : n++;
706 :
707 : //cout<<"SSD: "<<input->GetDDL()<<" "<<input->GetAD()<<" "
708 : //<<input->GetADC()<<" "<<input->GetSideFlag()<<" "<<((int)input->GetStrip())<<" "<<strip<<" : "<<input->GetSignal()<<endl;
709 : //
710 12465 : countRW++; //RS
711 12465 : } //* End main loop over the input
712 :
713 12 : AliDebug(1,Form("found clusters in ITS SSD: %d", nClustersSSD));
714 4 : }
715 :
716 :
717 : void AliITSClusterFinderV2SSD::
718 : FindClustersSSD(const Ali1Dcluster* neg, Int_t nn,
719 : const Ali1Dcluster* pos, Int_t np,
720 : TClonesArray *clusters) {
721 : //------------------------------------------------------------
722 : // Actual SSD cluster finder
723 : //------------------------------------------------------------
724 :
725 548 : const TGeoHMatrix *mT2L=AliITSgeomTGeo::GetTracking2LocalMatrix(fModule);
726 :
727 : //---------------------------------------
728 : // load recoparam
729 : //
730 : static AliITSRecoParam *repa = NULL;
731 274 : if(!repa){
732 2 : repa = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
733 2 : if(!repa){
734 0 : repa = AliITSRecoParam::GetHighFluxParam();
735 0 : AliWarning("Using default AliITSRecoParam class");
736 0 : }
737 : }
738 :
739 : // TClonesArray &cl=*clusters;
740 :
741 274 : AliITSsegmentationSSD *seg = static_cast<AliITSsegmentationSSD*>(fDetTypeRec->GetSegmentationModel(2));
742 548 : if (fModule>fLastSSD1)
743 412 : seg->SetLayer(6);
744 : else
745 136 : seg->SetLayer(5);
746 :
747 274 : Float_t hwSSD = seg->Dx()*1e-4/2;
748 274 : Float_t hlSSD = seg->Dz()*1e-4/2;
749 :
750 274 : Int_t idet=fNdet[fModule];
751 : Int_t ncl=0;
752 :
753 : //
754 274 : Int_t *cnegative = new Int_t[np];
755 274 : Int_t *cused1 = new Int_t[np];
756 274 : Int_t *negativepair = new Int_t[10*np];
757 274 : Int_t *cpositive = new Int_t[nn];
758 274 : Int_t *cused2 = new Int_t[nn];
759 274 : Int_t *positivepair = new Int_t[10*nn];
760 1144 : for (Int_t i=0;i<np;i++) {cnegative[i]=0; cused1[i]=0;}
761 1144 : for (Int_t i=0;i<nn;i++) {cpositive[i]=0; cused2[i]=0;}
762 6508 : for (Int_t i=0;i<10*np;i++) {negativepair[i]=0;}
763 6508 : for (Int_t i=0;i<10*nn;i++) {positivepair[i]=0;}
764 :
765 274 : if ((np*nn) > fgPairsSize) {
766 :
767 10 : delete [] fgPairs;
768 6 : fgPairsSize = 2*np*nn;
769 6 : fgPairs = new Short_t[fgPairsSize];
770 6 : }
771 274 : memset(fgPairs,0,sizeof(Short_t)*np*nn);
772 :
773 : //
774 : // find available pairs
775 : //
776 : Int_t ncross = 0;
777 1144 : for (Int_t i=0; i<np; i++) {
778 298 : Float_t yp=pos[i].GetY();
779 596 : if ( (pos[i].GetQ()>0) && (pos[i].GetQ()<3) ) continue;
780 1304 : for (Int_t j=0; j<nn; j++) {
781 708 : if ( (neg[j].GetQ()>0) && (neg[j].GetQ()<3) ) continue;
782 354 : Float_t yn=neg[j].GetY();
783 :
784 354 : Float_t xt, zt;
785 354 : seg->GetPadCxz(yn, yp, xt, zt);
786 : //cout<<yn<<" "<<yp<<" "<<xt<<" "<<zt<<endl;
787 :
788 354 : if (TMath::Abs(xt)<hwSSD)
789 354 : if (TMath::Abs(zt)<hlSSD) {
790 298 : Int_t in = i*10+cnegative[i];
791 298 : Int_t ip = j*10+cpositive[j];
792 596 : if ((in < 10*np) && (ip < 10*nn)) {
793 298 : negativepair[in] =j; //index
794 298 : positivepair[ip] =i;
795 298 : cnegative[i]++; //counters
796 298 : cpositive[j]++;
797 298 : ncross++;
798 298 : fgPairs[i*nn+j]=100;
799 298 : }
800 : else
801 0 : AliError(Form("Index out of range: ip=%d, in=%d",ip,in));
802 298 : }
803 354 : }
804 298 : }
805 :
806 274 : if (!ncross) {
807 16 : delete [] cnegative;
808 16 : delete [] cused1;
809 16 : delete [] negativepair;
810 16 : delete [] cpositive;
811 16 : delete [] cused2;
812 16 : delete [] positivepair;
813 8 : return;
814 : }
815 : //why not to allocate memorey here? if(!clusters) clusters = new TClonesArray("AliITSRecPoint", ncross);
816 :
817 : /* //
818 : // try to recover points out of but close to the module boundaries
819 : //
820 : for (Int_t i=0; i<np; i++) {
821 : Float_t yp=pos[i].GetY();
822 : if ( (pos[i].GetQ()>0) && (pos[i].GetQ()<3) ) continue;
823 : for (Int_t j=0; j<nn; j++) {
824 : if ( (neg[j].GetQ()>0) && (neg[j].GetQ()<3) ) continue;
825 : // if both 1Dclusters have an other cross continue
826 : if (cpositive[j]&&cnegative[i]) continue;
827 : Float_t yn=neg[j].GetY();
828 :
829 : Float_t xt, zt;
830 : seg->GetPadCxz(yn, yp, xt, zt);
831 :
832 : if (TMath::Abs(xt)<hwSSD+0.1)
833 : if (TMath::Abs(zt)<hlSSD+0.15) {
834 : // tag 1Dcluster (eventually will produce low quality recpoint)
835 : if (cnegative[i]==0) pos[i].SetNd(100); // not available pair
836 : if (cpositive[j]==0) neg[j].SetNd(100); // not available pair
837 : Int_t in = i*10+cnegative[i];
838 : Int_t ip = j*10+cpositive[j];
839 : if ((in < 10*np) && (ip < 10*nn)) {
840 : negativepair[in] =j; //index
841 : positivepair[ip] =i;
842 : cnegative[i]++; //counters
843 : cpositive[j]++;
844 : fgPairs[i*nn+j]=100;
845 : }
846 : else
847 : AliError(Form("Index out of range: ip=%d, in=%d",ip,in));
848 : }
849 : }
850 : }
851 : */
852 :
853 : //
854 266 : Float_t lp[6];
855 266 : Int_t milab[10];
856 : Double_t ratio;
857 :
858 :
859 266 : if(repa->GetUseChargeMatchingInClusterFinderSSD()==kTRUE) {
860 :
861 :
862 : //
863 : // sign gold tracks
864 : //
865 1112 : for (Int_t ip=0;ip<np;ip++){
866 : Float_t xbest=1000,zbest=1000,qbest=0;
867 : //
868 : // select gold clusters
869 572 : if ( (cnegative[ip]==1) && cpositive[negativepair[10*ip]]==1){
870 274 : Float_t yp=pos[ip].GetY();
871 274 : Int_t j = negativepair[10*ip];
872 :
873 274 : if( (pos[ip].GetQ()==0) && (neg[j].GetQ() ==0) ) {
874 : // both bad, hence continue;
875 : // mark both as used (to avoid recover at the end)
876 0 : cused1[ip]++;
877 0 : cused2[j]++;
878 0 : continue;
879 : }
880 :
881 274 : ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
882 : //cout<<"ratio="<<ratio<<endl;
883 :
884 : // charge matching (note that if posQ or negQ is 0 -> ratio=1 and the following condition is met
885 274 : if (TMath::Abs(ratio)>0.2) continue; // note: 0.2=3xsigma_ratio calculated in cosmics tests
886 :
887 : //
888 274 : Float_t yn=neg[j].GetY();
889 :
890 274 : Float_t xt, zt;
891 274 : seg->GetPadCxz(yn, yp, xt, zt);
892 :
893 274 : xbest=xt; zbest=zt;
894 :
895 :
896 274 : qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
897 548 : if( (pos[ip].GetQ()==0)||(neg[j].GetQ()==0)) qbest*=2; // in case of bad strips on one side keep all charge from the other one
898 :
899 : {
900 274 : Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
901 274 : mT2L->MasterToLocal(loc,trk);
902 274 : lp[0]=trk[1];
903 274 : lp[1]=trk[2];
904 274 : }
905 274 : lp[4]=qbest; //Q
906 6028 : for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
907 2192 : for (Int_t ilab=0;ilab<3;ilab++){
908 822 : milab[ilab] = pos[ip].GetLabel(ilab);
909 822 : milab[ilab+3] = neg[j].GetLabel(ilab);
910 : }
911 : //
912 274 : CheckLabels2(milab);
913 274 : milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
914 274 : Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
915 :
916 274 : lp[2]=4.968e-06; // 0.00223*0.00223; //SigmaY2
917 274 : lp[3]=0.012; // 0.110*0.110; //SigmaZ2
918 : // out-of-diagonal element of covariance matrix
919 536 : if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
920 286 : else if ( (info[0]>1) && (info[1]>1) ) {
921 58 : lp[2]=2.63e-06; // 0.0016*0.0016; //SigmaY2
922 58 : lp[3]=0.0065; // 0.08*0.08; //SigmaZ2
923 58 : lp[5]=-6.48e-05;
924 58 : }
925 : else {
926 118 : lp[2]=4.80e-06; // 0.00219*0.00219
927 118 : lp[3]=0.0093; // 0.0964*0.0964;
928 118 : if (info[0]==1) {
929 66 : lp[5]=-0.00014;
930 66 : }
931 : else {
932 52 : lp[2]=2.79e-06; // 0.0017*0.0017;
933 52 : lp[3]=0.00935; // 0.967*0.967;
934 52 : lp[5]=-4.32e-05;
935 : }
936 : }
937 :
938 : AliITSRecPoint * cl2;
939 274 : cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
940 :
941 274 : cl2->SetChargeRatio(ratio);
942 274 : cl2->SetType(1);
943 274 : fgPairs[ip*nn+j]=1;
944 :
945 274 : if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
946 14 : cl2->SetType(2);
947 14 : fgPairs[ip*nn+j]=2;
948 14 : }
949 :
950 274 : if(pos[ip].GetQ()==0) cl2->SetType(3);
951 274 : if(neg[j].GetQ()==0) cl2->SetType(4);
952 :
953 274 : cused1[ip]++;
954 274 : cused2[j]++;
955 :
956 274 : ncl++;
957 274 : }
958 290 : }
959 :
960 1112 : for (Int_t ip=0;ip<np;ip++){
961 : Float_t xbest=1000,zbest=1000,qbest=0;
962 : //
963 : //
964 : // select "silber" cluster
965 572 : if ( cnegative[ip]==1 && cpositive[negativepair[10*ip]]==2){
966 : Int_t in = negativepair[10*ip];
967 8 : Int_t ip2 = positivepair[10*in];
968 16 : if (ip2==ip) ip2 = positivepair[10*in+1];
969 8 : Float_t pcharge = pos[ip].GetQ()+pos[ip2].GetQ();
970 :
971 :
972 :
973 8 : ratio = (pcharge-neg[in].GetQ())/(pcharge+neg[in].GetQ());
974 8 : if ( (TMath::Abs(ratio)<0.2) && (pcharge!=0) ) {
975 : //if ( (TMath::Abs(pcharge-neg[in].GetQ())<30) && (pcharge!=0) ) { //
976 :
977 : //
978 : // add first pair
979 4 : if ( (fgPairs[ip*nn+in]==100)&&(pos[ip].GetQ() ) ) { //
980 :
981 2 : Float_t yp=pos[ip].GetY();
982 2 : Float_t yn=neg[in].GetY();
983 :
984 2 : Float_t xt, zt;
985 2 : seg->GetPadCxz(yn, yp, xt, zt);
986 :
987 2 : xbest=xt; zbest=zt;
988 :
989 2 : qbest =pos[ip].GetQ();
990 2 : Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
991 2 : mT2L->MasterToLocal(loc,trk);
992 2 : lp[0]=trk[1];
993 2 : lp[1]=trk[2];
994 :
995 2 : lp[4]=qbest; //Q
996 44 : for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
997 16 : for (Int_t ilab=0;ilab<3;ilab++){
998 6 : milab[ilab] = pos[ip].GetLabel(ilab);
999 6 : milab[ilab+3] = neg[in].GetLabel(ilab);
1000 : }
1001 : //
1002 2 : CheckLabels2(milab);
1003 2 : ratio = (pos[ip].GetQ()-neg[in].GetQ())/(pos[ip].GetQ()+neg[in].GetQ());
1004 2 : milab[3]=(((ip<<10) + in)<<10) + idet; // pos|neg|det
1005 2 : Int_t info[3] = {pos[ip].GetNd(),neg[in].GetNd(),fNlayer[fModule]};
1006 :
1007 2 : lp[2]=4.968e-06; // 0.00223*0.00223; //SigmaY2
1008 2 : lp[3]=0.012; // 0.110*0.110; //SigmaZ2
1009 : // out-of-diagonal element of covariance matrix
1010 6 : if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1011 0 : else if ( (info[0]>1) && (info[1]>1) ) {
1012 0 : lp[2]=2.63e-06; // 0.0016*0.0016; //SigmaY2
1013 0 : lp[3]=0.0065; // 0.08*0.08; //SigmaZ2
1014 0 : lp[5]=-6.48e-05;
1015 0 : }
1016 : else {
1017 0 : lp[2]=4.80e-06; // 0.00219*0.00219
1018 0 : lp[3]=0.0093; // 0.0964*0.0964;
1019 0 : if (info[0]==1) {
1020 0 : lp[5]=-0.00014;
1021 0 : }
1022 : else {
1023 0 : lp[2]=2.79e-06; // 0.0017*0.0017;
1024 0 : lp[3]=0.00935; // 0.967*0.967;
1025 0 : lp[5]=-4.32e-05;
1026 : }
1027 : }
1028 :
1029 : AliITSRecPoint * cl2;
1030 2 : cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1031 2 : cl2->SetChargeRatio(ratio);
1032 2 : cl2->SetType(5);
1033 2 : fgPairs[ip*nn+in] = 5;
1034 2 : if ((pos[ip].GetNd()+neg[in].GetNd())>6){ //multi cluster
1035 0 : cl2->SetType(6);
1036 0 : fgPairs[ip*nn+in] = 6;
1037 0 : }
1038 2 : ncl++;
1039 2 : }
1040 :
1041 :
1042 : //
1043 : // add second pair
1044 :
1045 : // if (!(cused1[ip2] || cused2[in])){ //
1046 4 : if ( (fgPairs[ip2*nn+in]==100) && (pos[ip2].GetQ()) ) {
1047 :
1048 2 : Float_t yp=pos[ip2].GetY();
1049 2 : Float_t yn=neg[in].GetY();
1050 :
1051 2 : Float_t xt, zt;
1052 2 : seg->GetPadCxz(yn, yp, xt, zt);
1053 :
1054 2 : xbest=xt; zbest=zt;
1055 :
1056 2 : qbest =pos[ip2].GetQ();
1057 :
1058 2 : Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
1059 2 : mT2L->MasterToLocal(loc,trk);
1060 2 : lp[0]=trk[1];
1061 2 : lp[1]=trk[2];
1062 :
1063 2 : lp[4]=qbest; //Q
1064 44 : for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1065 16 : for (Int_t ilab=0;ilab<3;ilab++){
1066 6 : milab[ilab] = pos[ip2].GetLabel(ilab);
1067 6 : milab[ilab+3] = neg[in].GetLabel(ilab);
1068 : }
1069 : //
1070 2 : CheckLabels2(milab);
1071 2 : ratio = (pos[ip2].GetQ()-neg[in].GetQ())/(pos[ip2].GetQ()+neg[in].GetQ());
1072 2 : milab[3]=(((ip2<<10) + in)<<10) + idet; // pos|neg|det
1073 2 : Int_t info[3] = {pos[ip2].GetNd(),neg[in].GetNd(),fNlayer[fModule]};
1074 :
1075 2 : lp[2]=4.968e-06; // 0.00223*0.00223; //SigmaY2
1076 2 : lp[3]=0.012; // 0.110*0.110; //SigmaZ2
1077 : // out-of-diagonal element of covariance matrix
1078 6 : if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1079 0 : else if ( (info[0]>1) && (info[1]>1) ) {
1080 0 : lp[2]=2.63e-06; // 0.0016*0.0016; //SigmaY2
1081 0 : lp[3]=0.0065; // 0.08*0.08; //SigmaZ2
1082 0 : lp[5]=-6.48e-05;
1083 0 : }
1084 : else {
1085 0 : lp[2]=4.80e-06; // 0.00219*0.00219
1086 0 : lp[3]=0.0093; // 0.0964*0.0964;
1087 0 : if (info[0]==1) {
1088 0 : lp[5]=-0.00014;
1089 0 : }
1090 : else {
1091 0 : lp[2]=2.79e-06; // 0.0017*0.0017;
1092 0 : lp[3]=0.00935; // 0.967*0.967;
1093 0 : lp[5]=-4.32e-05;
1094 : }
1095 : }
1096 :
1097 : AliITSRecPoint * cl2;
1098 2 : cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1099 :
1100 2 : cl2->SetChargeRatio(ratio);
1101 2 : cl2->SetType(5);
1102 2 : fgPairs[ip2*nn+in] =5;
1103 2 : if ((pos[ip2].GetNd()+neg[in].GetNd())>6){ //multi cluster
1104 0 : cl2->SetType(6);
1105 0 : fgPairs[ip2*nn+in] =6;
1106 0 : }
1107 2 : ncl++;
1108 2 : }
1109 :
1110 2 : cused1[ip]++;
1111 2 : cused1[ip2]++;
1112 2 : cused2[in]++;
1113 :
1114 2 : } // charge matching condition
1115 :
1116 8 : } // 2 Pside cross 1 Nside
1117 : } // loop over Pside clusters
1118 :
1119 :
1120 :
1121 : //
1122 1112 : for (Int_t jn=0;jn<nn;jn++){
1123 290 : if (cused2[jn]) continue;
1124 : Float_t xbest=1000,zbest=1000,qbest=0;
1125 : // select "silber" cluster
1126 22 : if ( cpositive[jn]==1 && cnegative[positivepair[10*jn]]==2){
1127 : Int_t ip = positivepair[10*jn];
1128 8 : Int_t jn2 = negativepair[10*ip];
1129 8 : if (jn2==jn) jn2 = negativepair[10*ip+1];
1130 8 : Float_t pcharge = neg[jn].GetQ()+neg[jn2].GetQ();
1131 : //
1132 :
1133 :
1134 8 : ratio = (pcharge-pos[ip].GetQ())/(pcharge+pos[ip].GetQ());
1135 8 : if ( (TMath::Abs(ratio)<0.2) && (pcharge!=0) ) {
1136 :
1137 : /*
1138 : if ( (TMath::Abs(pcharge-pos[ip].GetQ())<30) && // charge matching
1139 : (pcharge!=0) ) { // reject combinations of bad strips
1140 : */
1141 :
1142 :
1143 : //
1144 : // add first pair
1145 : // if (!(cused1[ip]||cused2[jn])){
1146 0 : if ( (fgPairs[ip*nn+jn]==100) && (neg[jn].GetQ()) ) { //
1147 :
1148 0 : Float_t yn=neg[jn].GetY();
1149 0 : Float_t yp=pos[ip].GetY();
1150 :
1151 0 : Float_t xt, zt;
1152 0 : seg->GetPadCxz(yn, yp, xt, zt);
1153 :
1154 0 : xbest=xt; zbest=zt;
1155 :
1156 0 : qbest =neg[jn].GetQ();
1157 :
1158 : {
1159 0 : Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
1160 0 : mT2L->MasterToLocal(loc,trk);
1161 0 : lp[0]=trk[1];
1162 0 : lp[1]=trk[2];
1163 0 : }
1164 :
1165 0 : lp[4]=qbest; //Q
1166 0 : for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1167 0 : for (Int_t ilab=0;ilab<3;ilab++){
1168 0 : milab[ilab] = pos[ip].GetLabel(ilab);
1169 0 : milab[ilab+3] = neg[jn].GetLabel(ilab);
1170 : }
1171 : //
1172 0 : CheckLabels2(milab);
1173 0 : ratio = (pos[ip].GetQ()-neg[jn].GetQ())/(pos[ip].GetQ()+neg[jn].GetQ());
1174 0 : milab[3]=(((ip<<10) + jn)<<10) + idet; // pos|neg|det
1175 0 : Int_t info[3] = {pos[ip].GetNd(),neg[jn].GetNd(),fNlayer[fModule]};
1176 :
1177 0 : lp[2]=4.968e-06; // 0.00223*0.00223; //SigmaY2
1178 0 : lp[3]=0.012; // 0.110*0.110; //SigmaZ2
1179 : // out-of-diagonal element of covariance matrix
1180 0 : if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1181 0 : else if ( (info[0]>1) && (info[1]>1) ) {
1182 0 : lp[2]=2.63e-06; // 0.0016*0.0016; //SigmaY2
1183 0 : lp[3]=0.0065; // 0.08*0.08; //SigmaZ2
1184 0 : lp[5]=-6.48e-05;
1185 0 : }
1186 : else {
1187 0 : lp[2]=4.80e-06; // 0.00219*0.00219
1188 0 : lp[3]=0.0093; // 0.0964*0.0964;
1189 0 : if (info[0]==1) {
1190 0 : lp[5]=-0.00014;
1191 0 : }
1192 : else {
1193 0 : lp[2]=2.79e-06; // 0.0017*0.0017;
1194 0 : lp[3]=0.00935; // 0.967*0.967;
1195 0 : lp[5]=-4.32e-05;
1196 : }
1197 : }
1198 :
1199 : AliITSRecPoint * cl2;
1200 0 : cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1201 :
1202 0 : cl2->SetChargeRatio(ratio);
1203 0 : cl2->SetType(7);
1204 0 : fgPairs[ip*nn+jn] =7;
1205 0 : if ((pos[ip].GetNd()+neg[jn].GetNd())>6){ //multi cluster
1206 0 : cl2->SetType(8);
1207 0 : fgPairs[ip*nn+jn]=8;
1208 0 : }
1209 0 : ncl++;
1210 0 : }
1211 : //
1212 : // add second pair
1213 : // if (!(cused1[ip]||cused2[jn2])){
1214 0 : if ( (fgPairs[ip*nn+jn2]==100)&&(neg[jn2].GetQ() ) ) { //
1215 :
1216 0 : Float_t yn=neg[jn2].GetY();
1217 0 : Double_t yp=pos[ip].GetY();
1218 :
1219 0 : Float_t xt, zt;
1220 0 : seg->GetPadCxz(yn, yp, xt, zt);
1221 :
1222 0 : xbest=xt; zbest=zt;
1223 :
1224 0 : qbest =neg[jn2].GetQ();
1225 :
1226 : {
1227 0 : Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
1228 0 : mT2L->MasterToLocal(loc,trk);
1229 0 : lp[0]=trk[1];
1230 0 : lp[1]=trk[2];
1231 0 : }
1232 :
1233 0 : lp[4]=qbest; //Q
1234 0 : for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1235 0 : for (Int_t ilab=0;ilab<3;ilab++){
1236 0 : milab[ilab] = pos[ip].GetLabel(ilab);
1237 0 : milab[ilab+3] = neg[jn2].GetLabel(ilab);
1238 : }
1239 : //
1240 0 : CheckLabels2(milab);
1241 0 : ratio = (pos[ip].GetQ()-neg[jn2].GetQ())/(pos[ip].GetQ()+neg[jn2].GetQ());
1242 0 : milab[3]=(((ip<<10) + jn2)<<10) + idet; // pos|neg|det
1243 0 : Int_t info[3] = {pos[ip].GetNd(),neg[jn2].GetNd(),fNlayer[fModule]};
1244 :
1245 0 : lp[2]=4.968e-06; // 0.00223*0.00223; //SigmaY2
1246 0 : lp[3]=0.012; // 0.110*0.110; //SigmaZ2
1247 : // out-of-diagonal element of covariance matrix
1248 0 : if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1249 0 : else if ( (info[0]>1) && (info[1]>1) ) {
1250 0 : lp[2]=2.63e-06; // 0.0016*0.0016; //SigmaY2
1251 0 : lp[3]=0.0065; // 0.08*0.08; //SigmaZ2
1252 0 : lp[5]=-6.48e-05;
1253 0 : }
1254 : else {
1255 0 : lp[2]=4.80e-06; // 0.00219*0.00219
1256 0 : lp[3]=0.0093; // 0.0964*0.0964;
1257 0 : if (info[0]==1) {
1258 0 : lp[5]=-0.00014;
1259 0 : }
1260 : else {
1261 0 : lp[2]=2.79e-06; // 0.0017*0.0017;
1262 0 : lp[3]=0.00935; // 0.967*0.967;
1263 0 : lp[5]=-4.32e-05;
1264 : }
1265 : }
1266 :
1267 : AliITSRecPoint * cl2;
1268 0 : cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1269 :
1270 :
1271 0 : cl2->SetChargeRatio(ratio);
1272 0 : fgPairs[ip*nn+jn2]=7;
1273 0 : cl2->SetType(7);
1274 0 : if ((pos[ip].GetNd()+neg[jn2].GetNd())>6){ //multi cluster
1275 0 : cl2->SetType(8);
1276 0 : fgPairs[ip*nn+jn2]=8;
1277 0 : }
1278 0 : ncl++;
1279 0 : }
1280 0 : cused1[ip]++;
1281 0 : cused2[jn]++;
1282 0 : cused2[jn2]++;
1283 :
1284 0 : } // charge matching condition
1285 :
1286 8 : } // 2 Nside cross 1 Pside
1287 14 : } // loop over Pside clusters
1288 :
1289 :
1290 :
1291 1112 : for (Int_t ip=0;ip<np;ip++){
1292 :
1293 290 : if(cused1[ip]) continue;
1294 :
1295 :
1296 : Float_t xbest=1000,zbest=1000,qbest=0;
1297 : //
1298 : // 2x2 clusters
1299 : //
1300 18 : if ( (cnegative[ip]==2) && cpositive[negativepair[10*ip]]==2){
1301 : Float_t minchargediff =4.;
1302 : // Float_t minchargeratio =0.2;
1303 :
1304 : Int_t j=-1;
1305 36 : for (Int_t di=0;di<cnegative[ip];di++){
1306 12 : Int_t jc = negativepair[ip*10+di];
1307 12 : Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
1308 12 : ratio = (pos[ip].GetQ()-neg[jc].GetQ())/(pos[ip].GetQ()+neg[jc].GetQ());
1309 : //if (TMath::Abs(chargedif)<minchargediff){
1310 12 : if (TMath::Abs(ratio)<0.2){
1311 : j =jc;
1312 10 : minchargediff = TMath::Abs(chargedif);
1313 : // minchargeratio = TMath::Abs(ratio);
1314 10 : }
1315 : }
1316 6 : if (j<0) continue; // not proper cluster
1317 :
1318 :
1319 : Int_t count =0;
1320 36 : for (Int_t di=0;di<cnegative[ip];di++){
1321 12 : Int_t jc = negativepair[ip*10+di];
1322 12 : Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
1323 20 : if (TMath::Abs(chargedif)<minchargediff+3.) count++;
1324 : }
1325 8 : if (count>1) continue; // more than one "proper" cluster for positive
1326 : //
1327 :
1328 : count =0;
1329 16 : for (Int_t dj=0;dj<cpositive[j];dj++){
1330 4 : Int_t ic = positivepair[j*10+dj];
1331 4 : Float_t chargedif = pos[ic].GetQ()-neg[j].GetQ();
1332 8 : if (TMath::Abs(chargedif)<minchargediff+3.) count++;
1333 : }
1334 4 : if (count>1) continue; // more than one "proper" cluster for negative
1335 :
1336 : Int_t jp = 0;
1337 :
1338 : count =0;
1339 16 : for (Int_t dj=0;dj<cnegative[jp];dj++){
1340 4 : Int_t ic = positivepair[jp*10+dj];
1341 4 : Float_t chargedif = pos[ic].GetQ()-neg[jp].GetQ();
1342 8 : if (TMath::Abs(chargedif)<minchargediff+4.) count++;
1343 : }
1344 4 : if (count>1) continue;
1345 4 : if (fgPairs[ip*nn+j]<100) continue;
1346 : //
1347 :
1348 :
1349 :
1350 : //almost gold clusters
1351 4 : Float_t yp=pos[ip].GetY();
1352 4 : Float_t yn=neg[j].GetY();
1353 4 : Float_t xt, zt;
1354 4 : seg->GetPadCxz(yn, yp, xt, zt);
1355 4 : xbest=xt; zbest=zt;
1356 4 : qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
1357 : {
1358 4 : Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
1359 4 : mT2L->MasterToLocal(loc,trk);
1360 4 : lp[0]=trk[1];
1361 4 : lp[1]=trk[2];
1362 4 : }
1363 4 : lp[4]=qbest; //Q
1364 88 : for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1365 32 : for (Int_t ilab=0;ilab<3;ilab++){
1366 12 : milab[ilab] = pos[ip].GetLabel(ilab);
1367 12 : milab[ilab+3] = neg[j].GetLabel(ilab);
1368 : }
1369 : //
1370 4 : CheckLabels2(milab);
1371 4 : if ((neg[j].GetQ()==0)&&(pos[ip].GetQ()==0)) continue; // reject crosses of bad strips!!
1372 4 : ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
1373 4 : milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
1374 4 : Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
1375 :
1376 4 : lp[2]=4.968e-06; // 0.00223*0.00223; //SigmaY2
1377 4 : lp[3]=0.012; // 0.110*0.110; //SigmaZ2
1378 : // out-of-diagonal element of covariance matrix
1379 10 : if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1380 2 : else if ( (info[0]>1) && (info[1]>1) ) {
1381 0 : lp[2]=2.63e-06; // 0.0016*0.0016; //SigmaY2
1382 0 : lp[3]=0.0065; // 0.08*0.08; //SigmaZ2
1383 0 : lp[5]=-6.48e-05;
1384 0 : }
1385 : else {
1386 2 : lp[2]=4.80e-06; // 0.00219*0.00219
1387 2 : lp[3]=0.0093; // 0.0964*0.0964;
1388 2 : if (info[0]==1) {
1389 2 : lp[5]=-0.00014;
1390 2 : }
1391 : else {
1392 0 : lp[2]=2.79e-06; // 0.0017*0.0017;
1393 0 : lp[3]=0.00935; // 0.967*0.967;
1394 0 : lp[5]=-4.32e-05;
1395 : }
1396 : }
1397 :
1398 : AliITSRecPoint * cl2;
1399 4 : cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1400 :
1401 4 : cl2->SetChargeRatio(ratio);
1402 4 : cl2->SetType(10);
1403 4 : fgPairs[ip*nn+j]=10;
1404 4 : if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
1405 0 : cl2->SetType(11);
1406 0 : fgPairs[ip*nn+j]=11;
1407 0 : }
1408 4 : cused1[ip]++;
1409 4 : cused2[j]++;
1410 4 : ncl++;
1411 :
1412 8 : } // 2X2
1413 10 : } // loop over Pside 1Dclusters
1414 :
1415 :
1416 :
1417 1112 : for (Int_t ip=0;ip<np;ip++){
1418 :
1419 290 : if(cused1[ip]) continue;
1420 :
1421 :
1422 : Float_t xbest=1000,zbest=1000,qbest=0;
1423 : //
1424 : // manyxmany clusters
1425 : //
1426 16 : if ( (cnegative[ip]<5) && cpositive[negativepair[10*ip]]<5){
1427 : Float_t minchargediff =4.;
1428 : Int_t j=-1;
1429 36 : for (Int_t di=0;di<cnegative[ip];di++){
1430 10 : Int_t jc = negativepair[ip*10+di];
1431 10 : Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
1432 10 : if (TMath::Abs(chargedif)<minchargediff){
1433 : j =jc;
1434 4 : minchargediff = TMath::Abs(chargedif);
1435 4 : }
1436 : }
1437 13 : if (j<0) continue; // not proper cluster
1438 :
1439 : Int_t count =0;
1440 16 : for (Int_t di=0;di<cnegative[ip];di++){
1441 5 : Int_t jc = negativepair[ip*10+di];
1442 5 : Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
1443 10 : if (TMath::Abs(chargedif)<minchargediff+3.) count++;
1444 : }
1445 5 : if (count>1) continue; // more than one "proper" cluster for positive
1446 : //
1447 :
1448 : count =0;
1449 6 : for (Int_t dj=0;dj<cpositive[j];dj++){
1450 2 : Int_t ic = positivepair[j*10+dj];
1451 2 : Float_t chargedif = pos[ic].GetQ()-neg[j].GetQ();
1452 3 : if (TMath::Abs(chargedif)<minchargediff+3.) count++;
1453 : }
1454 1 : if (count>1) continue; // more than one "proper" cluster for negative
1455 :
1456 : Int_t jp = 0;
1457 :
1458 : count =0;
1459 4 : for (Int_t dj=0;dj<cnegative[jp];dj++){
1460 1 : Int_t ic = positivepair[jp*10+dj];
1461 1 : Float_t chargedif = pos[ic].GetQ()-neg[jp].GetQ();
1462 2 : if (TMath::Abs(chargedif)<minchargediff+4.) count++;
1463 : }
1464 1 : if (count>1) continue;
1465 1 : if (fgPairs[ip*nn+j]<100) continue;
1466 : //
1467 :
1468 : //almost gold clusters
1469 1 : Float_t yp=pos[ip].GetY();
1470 1 : Float_t yn=neg[j].GetY();
1471 :
1472 :
1473 1 : Float_t xt, zt;
1474 1 : seg->GetPadCxz(yn, yp, xt, zt);
1475 :
1476 1 : xbest=xt; zbest=zt;
1477 :
1478 1 : qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
1479 :
1480 : {
1481 1 : Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
1482 1 : mT2L->MasterToLocal(loc,trk);
1483 1 : lp[0]=trk[1];
1484 1 : lp[1]=trk[2];
1485 1 : }
1486 1 : lp[4]=qbest; //Q
1487 22 : for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1488 8 : for (Int_t ilab=0;ilab<3;ilab++){
1489 3 : milab[ilab] = pos[ip].GetLabel(ilab);
1490 3 : milab[ilab+3] = neg[j].GetLabel(ilab);
1491 : }
1492 : //
1493 1 : CheckLabels2(milab);
1494 1 : if ((neg[j].GetQ()==0)&&(pos[ip].GetQ()==0)) continue; // reject crosses of bad strips!!
1495 1 : ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
1496 1 : milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
1497 1 : Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
1498 :
1499 1 : lp[2]=4.968e-06; // 0.00223*0.00223; //SigmaY2
1500 1 : lp[3]=0.012; // 0.110*0.110; //SigmaZ2
1501 : // out-of-diagonal element of covariance matrix
1502 2 : if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1503 1 : else if ( (info[0]>1) && (info[1]>1) ) {
1504 0 : lp[2]=2.63e-06; // 0.0016*0.0016; //SigmaY2
1505 0 : lp[3]=0.0065; // 0.08*0.08; //SigmaZ2
1506 0 : lp[5]=-6.48e-05;
1507 0 : }
1508 : else {
1509 1 : lp[2]=4.80e-06; // 0.00219*0.00219
1510 1 : lp[3]=0.0093; // 0.0964*0.0964;
1511 1 : if (info[0]==1) {
1512 1 : lp[5]=-0.00014;
1513 1 : }
1514 : else {
1515 0 : lp[2]=2.79e-06; // 0.0017*0.0017;
1516 0 : lp[3]=0.00935; // 0.967*0.967;
1517 0 : lp[5]=-4.32e-05;
1518 : }
1519 : }
1520 : //
1521 1 : if (fRawID2ClusID) { // set rawID <-> clusterID correspondence for embedding
1522 : const int kMaxRefRW = 200;
1523 0 : UInt_t nrefsRW,refsRW[kMaxRefRW];
1524 0 : nrefsRW = fRawIDRef[0].GetReferences(j,refsRW,kMaxRefRW); // n-side
1525 0 : for (int ir=nrefsRW;ir--;) {
1526 0 : int rwid = (int)refsRW[ir];
1527 0 : if (fRawID2ClusID->GetSize()<=rwid) fRawID2ClusID->Set( (rwid+10)<<1 );
1528 0 : (*fRawID2ClusID)[rwid] = fNClusters+1; // RS: store clID+1 as a reference to the cluster
1529 : }
1530 : //
1531 0 : nrefsRW = fRawIDRef[1].GetReferences(ip,refsRW,kMaxRefRW); // p-side
1532 0 : for (int ir=nrefsRW;ir--;) {
1533 0 : int rwid = (int)refsRW[ir];
1534 0 : if (fRawID2ClusID->GetSize()<=rwid) fRawID2ClusID->Set( (rwid+10)<<1 );
1535 0 : (*fRawID2ClusID)[rwid] = fNClusters+1; // RS: store clID+1 as a reference to the cluster
1536 : }
1537 : //
1538 0 : milab[0] = fNClusters+1; // RS: assign id as cluster label
1539 0 : }
1540 :
1541 : AliITSRecPoint * cl2;
1542 1 : cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1543 :
1544 1 : cl2->SetChargeRatio(ratio);
1545 1 : cl2->SetType(12);
1546 1 : fgPairs[ip*nn+j]=12;
1547 1 : if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
1548 0 : cl2->SetType(13);
1549 0 : fgPairs[ip*nn+j]=13;
1550 0 : }
1551 1 : cused1[ip]++;
1552 1 : cused2[j]++;
1553 1 : ncl++;
1554 1 : fNClusters++;
1555 :
1556 2 : } // manyXmany
1557 1 : } // loop over Pside 1Dclusters
1558 :
1559 266 : } // use charge matching
1560 :
1561 :
1562 : // recover all the other crosses
1563 : //
1564 1112 : for (Int_t i=0; i<np; i++) {
1565 : Float_t xbest=1000,zbest=1000,qbest=0;
1566 290 : Float_t yp=pos[i].GetY();
1567 580 : if ((pos[i].GetQ()>0)&&(pos[i].GetQ()<3)) continue;
1568 1272 : for (Int_t j=0; j<nn; j++) {
1569 : // for (Int_t di = 0;di<cpositive[i];di++){
1570 : // Int_t j = negativepair[10*i+di];
1571 692 : if ((neg[j].GetQ()>0)&&(neg[j].GetQ()<3)) continue;
1572 :
1573 346 : if ((neg[j].GetQ()==0)&&(pos[i].GetQ()==0)) continue; // reject crosses of bad strips!!
1574 :
1575 368 : if (cused2[j]||cused1[i]) continue;
1576 20 : if (fgPairs[i*nn+j]>0 &&fgPairs[i*nn+j]<100) continue;
1577 11 : ratio = (pos[i].GetQ()-neg[j].GetQ())/(pos[i].GetQ()+neg[j].GetQ());
1578 11 : Float_t yn=neg[j].GetY();
1579 :
1580 11 : Float_t xt, zt;
1581 11 : seg->GetPadCxz(yn, yp, xt, zt);
1582 :
1583 11 : if (TMath::Abs(xt)<hwSSD)
1584 11 : if (TMath::Abs(zt)<hlSSD) {
1585 9 : xbest=xt; zbest=zt;
1586 :
1587 9 : qbest=0.5*(pos[i].GetQ()+neg[j].GetQ());
1588 :
1589 : {
1590 9 : Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
1591 9 : mT2L->MasterToLocal(loc,trk);
1592 9 : lp[0]=trk[1];
1593 9 : lp[1]=trk[2];
1594 9 : }
1595 9 : lp[4]=qbest; //Q
1596 198 : for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1597 72 : for (Int_t ilab=0;ilab<3;ilab++){
1598 27 : milab[ilab] = pos[i].GetLabel(ilab);
1599 27 : milab[ilab+3] = neg[j].GetLabel(ilab);
1600 : }
1601 : //
1602 9 : CheckLabels2(milab);
1603 9 : milab[3]=(((i<<10) + j)<<10) + idet; // pos|neg|det
1604 9 : Int_t info[3] = {pos[i].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
1605 :
1606 9 : lp[2]=4.968e-06; // 0.00223*0.00223; //SigmaY2
1607 9 : lp[3]=0.012; // 0.110*0.110; //SigmaZ2
1608 : // out-of-diagonal element of covariance matrix
1609 20 : if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1610 7 : else if ( (info[0]>1) && (info[1]>1) ) {
1611 0 : lp[2]=2.63e-06; // 0.0016*0.0016; //SigmaY2
1612 0 : lp[3]=0.0065; // 0.08*0.08; //SigmaZ2
1613 0 : lp[5]=-6.48e-05;
1614 0 : }
1615 : else {
1616 5 : lp[2]=4.80e-06; // 0.00219*0.00219
1617 5 : lp[3]=0.0093; // 0.0964*0.0964;
1618 5 : if (info[0]==1) {
1619 3 : lp[5]=-0.00014;
1620 3 : }
1621 : else {
1622 2 : lp[2]=2.79e-06; // 0.0017*0.0017;
1623 2 : lp[3]=0.00935; // 0.967*0.967;
1624 2 : lp[5]=-4.32e-05;
1625 : }
1626 : }
1627 :
1628 : AliITSRecPoint * cl2;
1629 9 : cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1630 :
1631 9 : cl2->SetChargeRatio(ratio);
1632 9 : cl2->SetType(100+cpositive[j]+cnegative[i]);
1633 :
1634 9 : if(pos[i].GetQ()==0) cl2->SetType(200+cpositive[j]+cnegative[i]);
1635 9 : if(neg[j].GetQ()==0) cl2->SetType(300+cpositive[j]+cnegative[i]);
1636 9 : ncl++;
1637 9 : }
1638 11 : }
1639 290 : }
1640 :
1641 :
1642 :
1643 266 : if(repa->GetUseBadChannelsInClusterFinderSSD()==kTRUE) {
1644 :
1645 : //---------------------------------------------------------
1646 : // recover crosses of good 1D clusters with bad strips on the other side
1647 : // Note1: at first iteration skip modules with a bad side (or almost), (would produce too many fake!)
1648 : // Note2: for modules with a bad side see below
1649 :
1650 0 : AliITSCalibrationSSD* cal = (AliITSCalibrationSSD*)GetResp(fModule);
1651 : Int_t countPbad=0, countNbad=0;
1652 0 : for(Int_t ib=0; ib<768; ib++) {
1653 0 : if(cal->IsPChannelBad(ib)) countPbad++;
1654 0 : if(cal->IsNChannelBad(ib)) countNbad++;
1655 : }
1656 : // AliInfo(Form("module %d has %d P- and %d N-bad strips",fModule,countPbad,countNbad));
1657 :
1658 0 : if( (countPbad<100) && (countNbad<100) ) { // no bad side!!
1659 :
1660 0 : for (Int_t i=0; i<np; i++) { // loop over Nside 1Dclusters with no crosses
1661 0 : if(cnegative[i]) continue; // if intersecting Pside clusters continue;
1662 :
1663 : // for(Int_t ib=0; ib<768; ib++) { // loop over all Pstrips
1664 0 : for(Int_t ib=15; ib<753; ib++) { // loop over all Pstrips
1665 :
1666 0 : if(cal->IsPChannelBad(ib)) { // check if strips is bad
1667 0 : Float_t yN=pos[i].GetY();
1668 0 : Float_t xt, zt;
1669 0 : seg->GetPadCxz(1.*ib, yN, xt, zt);
1670 :
1671 : //----------
1672 : // bad Pstrip is crossing the Nside 1Dcluster -> create recpoint
1673 : //
1674 0 : if ( (TMath::Abs(xt)<hwSSD) && (TMath::Abs(zt)<hlSSD) ) {
1675 0 : Double_t loc[3]={xt,0.,zt},trk[3]={0.,0.,0.};
1676 0 : mT2L->MasterToLocal(loc,trk);
1677 0 : lp[0]=trk[1];
1678 0 : lp[1]=trk[2];
1679 0 : lp[4]=pos[i].GetQ(); //Q
1680 0 : for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1681 0 : for (Int_t ilab=0;ilab<3;ilab++) milab[ilab] = pos[i].GetLabel(ilab);
1682 0 : CheckLabels2(milab);
1683 0 : milab[3]=( (i<<10) << 10 ) + idet; // pos|neg|det
1684 0 : Int_t info[3] = {pos[i].GetNd(),0,fNlayer[fModule]};
1685 :
1686 0 : lp[2]=4.968e-06; // 0.00223*0.00223; //SigmaY2
1687 0 : lp[3]=0.012; // 0.110*0.110; //SigmaZ2
1688 0 : lp[5]=-0.00012; // out-of-diagonal element of covariance matrix
1689 0 : if (info[0]>1) {
1690 0 : lp[2]=4.80e-06;
1691 0 : lp[3]=0.0093;
1692 0 : lp[5]=0.00014;
1693 0 : }
1694 :
1695 : AliITSRecPoint * cl2;
1696 0 : cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1697 0 : cl2->SetChargeRatio(1.);
1698 0 : cl2->SetType(50);
1699 0 : ncl++;
1700 0 : } // cross is within the detector
1701 : //
1702 : //--------------
1703 :
1704 0 : } // bad Pstrip
1705 :
1706 : } // end loop over Pstrips
1707 :
1708 0 : } // end loop over Nside 1D clusters
1709 :
1710 0 : for (Int_t j=0; j<nn; j++) { // loop over Pside 1D clusters with no crosses
1711 0 : if(cpositive[j]) continue;
1712 :
1713 : // for(Int_t ib=0; ib<768; ib++) { // loop over all Nside strips
1714 0 : for(Int_t ib=15; ib<753; ib++) { // loop over all Nside strips
1715 :
1716 0 : if(cal->IsNChannelBad(ib)) { // check if strip is bad
1717 0 : Float_t yP=neg[j].GetY();
1718 0 : Float_t xt, zt;
1719 0 : seg->GetPadCxz(yP, 1.*ib, xt, zt);
1720 :
1721 : //----------
1722 : // bad Nstrip is crossing the Pside 1Dcluster -> create recpoint
1723 : //
1724 0 : if ( (TMath::Abs(xt)<hwSSD) && (TMath::Abs(zt)<hlSSD) ) {
1725 0 : Double_t loc[3]={xt,0.,zt},trk[3]={0.,0.,0.};
1726 0 : mT2L->MasterToLocal(loc,trk);
1727 0 : lp[0]=trk[1];
1728 0 : lp[1]=trk[2];
1729 0 : lp[4]=neg[j].GetQ(); //Q
1730 0 : for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1731 0 : for (Int_t ilab=0;ilab<3;ilab++) milab[ilab] = neg[j].GetLabel(ilab);
1732 0 : CheckLabels2(milab);
1733 0 : milab[3]=( j << 10 ) + idet; // pos|neg|det
1734 0 : Int_t info[3]={0,(Int_t)neg[j].GetNd(),fNlayer[fModule]};
1735 :
1736 0 : lp[2]=4.968e-06; // 0.00223*0.00223; //SigmaY2
1737 0 : lp[3]=0.012; // 0.110*0.110; //SigmaZ2
1738 0 : lp[5]=-0.00012; // out-of-diagonal element of covariance matrix
1739 0 : if (info[0]>1) {
1740 0 : lp[2]=2.79e-06;
1741 0 : lp[3]=0.00935;
1742 0 : lp[5]=-4.32e-05;
1743 0 : }
1744 :
1745 : AliITSRecPoint * cl2;
1746 0 : cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1747 0 : cl2->SetChargeRatio(1.);
1748 0 : cl2->SetType(60);
1749 0 : ncl++;
1750 0 : } // cross is within the detector
1751 : //
1752 : //--------------
1753 :
1754 0 : } // bad Nstrip
1755 : } // end loop over Nstrips
1756 0 : } // end loop over Pside 1D clusters
1757 :
1758 0 : } // no bad sides
1759 :
1760 : //---------------------------------------------------------
1761 :
1762 0 : else if( (countPbad>700) && (countNbad<100) ) { // bad Pside!!
1763 :
1764 0 : for (Int_t i=0; i<np; i++) { // loop over Nside 1Dclusters with no crosses
1765 0 : if(cnegative[i]) continue; // if intersecting Pside clusters continue;
1766 :
1767 0 : Float_t xt, zt;
1768 0 : Float_t yN=pos[i].GetY();
1769 : Float_t yP=0.;
1770 0 : if (seg->GetLayer()==5) yP = yN + (7.6/1.9);
1771 0 : else yP = yN - (7.6/1.9);
1772 0 : seg->GetPadCxz(yP, yN, xt, zt);
1773 :
1774 0 : if ( (TMath::Abs(xt)<hwSSD) && (TMath::Abs(zt)<hlSSD) ) {
1775 0 : Double_t loc[3]={xt,0.,zt},trk[3]={0.,0.,0.};
1776 0 : mT2L->MasterToLocal(loc,trk);
1777 0 : lp[0]=trk[1];
1778 0 : lp[1]=trk[2];
1779 0 : lp[4]=pos[i].GetQ(); //Q
1780 0 : for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1781 0 : for (Int_t ilab=0;ilab<3;ilab++) milab[ilab] = pos[i].GetLabel(ilab);
1782 0 : CheckLabels2(milab);
1783 0 : milab[3]=( (i<<10) << 10 ) + idet; // pos|neg|det
1784 0 : Int_t info[3] = {(Int_t)pos[i].GetNd(),0,fNlayer[fModule]};
1785 :
1786 0 : lp[2]=0.00098; // 0.031*0.031; //SigmaY2
1787 0 : lp[3]=1.329; // 1.15*1.15; //SigmaZ2
1788 0 : lp[5]=-0.0359;
1789 0 : if(info[0]>1) lp[2]=0.00097;
1790 :
1791 : AliITSRecPoint * cl2;
1792 0 : cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1793 0 : cl2->SetChargeRatio(1.);
1794 0 : cl2->SetType(70);
1795 0 : ncl++;
1796 0 : } // cross is within the detector
1797 : //
1798 : //--------------
1799 :
1800 0 : } // end loop over Nside 1D clusters
1801 :
1802 0 : } // bad Pside module
1803 :
1804 0 : else if( (countNbad>700) && (countPbad<100) ) { // bad Nside!!
1805 :
1806 0 : for (Int_t j=0; j<nn; j++) { // loop over Pside 1D clusters with no crosses
1807 0 : if(cpositive[j]) continue;
1808 :
1809 0 : Float_t xt, zt;
1810 0 : Float_t yP=neg[j].GetY();
1811 : Float_t yN=0.;
1812 0 : if (seg->GetLayer()==5) yN = yP - (7.6/1.9);
1813 0 : else yN = yP + (7.6/1.9);
1814 0 : seg->GetPadCxz(yP, yN, xt, zt);
1815 :
1816 0 : if ( (TMath::Abs(xt)<hwSSD) && (TMath::Abs(zt)<hlSSD) ) {
1817 0 : Double_t loc[3]={xt,0.,zt},trk[3]={0.,0.,0.};
1818 0 : mT2L->MasterToLocal(loc,trk);
1819 0 : lp[0]=trk[1];
1820 0 : lp[1]=trk[2];
1821 0 : lp[4]=neg[j].GetQ(); //Q
1822 0 : for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1823 0 : for (Int_t ilab=0;ilab<3;ilab++) milab[ilab] = neg[j].GetLabel(ilab);
1824 0 : CheckLabels2(milab);
1825 0 : milab[3]=( j << 10 ) + idet; // pos|neg|det
1826 0 : Int_t info[3] = {0,(Int_t)neg[j].GetNd(),fNlayer[fModule]};
1827 :
1828 0 : lp[2]=7.27e-05; // 0.0085*0.0085; //SigmaY2
1829 0 : lp[3]=1.33; // 1.15*1.15; //SigmaZ2
1830 0 : lp[5]=0.00931;
1831 0 : if(info[1]>1) lp[2]=6.91e-05;
1832 :
1833 : AliITSRecPoint * cl2;
1834 0 : cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1835 0 : cl2->SetChargeRatio(1.);
1836 0 : cl2->SetType(80);
1837 0 : ncl++;
1838 0 : } // cross is within the detector
1839 : //
1840 : //--------------
1841 :
1842 0 : } // end loop over Pside 1D clusters
1843 :
1844 0 : } // bad Nside module
1845 :
1846 : //---------------------------------------------------------
1847 :
1848 0 : } // use bad channels
1849 :
1850 : //cout<<ncl<<" clusters for this module"<<endl;
1851 :
1852 532 : delete [] cnegative;
1853 532 : delete [] cused1;
1854 532 : delete [] negativepair;
1855 532 : delete [] cpositive;
1856 532 : delete [] cused2;
1857 532 : delete [] positivepair;
1858 :
1859 540 : }
|
