1 /**************************************************************************
2 * Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 ////////////////////////////////////////////////////////////////////////////
19 // Implementation of the ITS clusterer V2 class //
21 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch //
22 // Last revision: 13-05-09 Enrico Fragiacomo //
23 // enrico.fragiacomo@ts.infn.it //
25 ///////////////////////////////////////////////////////////////////////////
27 #include <Riostream.h>
30 #include "AliITSClusterFinderV2SSD.h"
31 #include "AliITSRecPoint.h"
32 #include "AliITSgeomTGeo.h"
33 #include "AliITSDetTypeRec.h"
34 #include "AliRawReader.h"
35 #include "AliITSRawStreamSSD.h"
36 #include <TClonesArray.h>
37 #include "AliITSdigitSSD.h"
38 #include "AliITSReconstructor.h"
39 #include "AliITSCalibrationSSD.h"
40 #include "AliITSsegmentationSSD.h"
42 Short_t *AliITSClusterFinderV2SSD::fgPairs = 0x0;
43 Int_t AliITSClusterFinderV2SSD::fgPairsSize = 0;
44 const Float_t AliITSClusterFinderV2SSD::fgkThreshold = 5.;
46 const Float_t AliITSClusterFinderV2SSD::fgkCosmic2008StripShifts[16][9] =
47 {{-0.35,-0.35,-0.35,-0.35,-0.35,-0.35,-0.35,-0.35,-0.35}, // DDL 512
48 {-0.35,-0.35,-0.35,-0.35,-0.35,-0.35,-0.35,-0.35,-0.35}, // DDL 513
49 {-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15}, // DDL 514
50 {-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15}, // DDL 515
51 { 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00}, // DDL 516
52 { 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00}, // DDL 517
53 {-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15}, // DDL 518
54 {-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15}, // DDL 519
55 {-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.25,-0.15}, // DDL 520
56 {-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15}, // DDL 521
57 {-0.10,-0.10,-0.10,-0.40,-0.40,-0.40,-0.10,-0.10,-0.45}, // DDL 522
58 {-0.10,-0.10,-0.10,-0.35,-0.35,-0.35,-0.10,-0.35,-0.50}, // DDL 523
59 { 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00}, // DDL 524
60 { 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00}, // DDL 525
61 { 0.35, 0.35, 0.35, 0.35, 0.35, 0.35, 0.35, 0.35, 0.35}, // DDL 526
62 { 0.45, 0.45, 0.45, 0.45, 0.45, 0.45, 0.45, 0.45, 0.45}}; // DDL 527
64 ClassImp(AliITSClusterFinderV2SSD)
67 AliITSClusterFinderV2SSD::AliITSClusterFinderV2SSD(AliITSDetTypeRec* dettyp):AliITSClusterFinder(dettyp),
68 fLastSSD1(AliITSgeomTGeo::GetModuleIndex(6,1,1)-1)
74 //______________________________________________________________________
75 AliITSClusterFinderV2SSD::AliITSClusterFinderV2SSD(const AliITSClusterFinderV2SSD &cf) : AliITSClusterFinder(cf), fLastSSD1(cf.fLastSSD1)
80 //______________________________________________________________________
81 AliITSClusterFinderV2SSD& AliITSClusterFinderV2SSD::operator=(const AliITSClusterFinderV2SSD& cf ){
82 // Assignment operator
84 this->~AliITSClusterFinderV2SSD();
85 new(this) AliITSClusterFinderV2SSD(cf);
90 void AliITSClusterFinderV2SSD::FindRawClusters(Int_t mod){
94 FindClustersSSD(fDigits);
98 void AliITSClusterFinderV2SSD::FindClustersSSD(TClonesArray *alldigits) {
99 //------------------------------------------------------------
100 // Actual SSD cluster finder
101 //------------------------------------------------------------
102 Int_t smaxall=alldigits->GetEntriesFast();
103 if (smaxall==0) return;
106 //---------------------------------------
107 // load recoparam and calibration
109 static AliITSRecoParam *repa = NULL;
111 repa = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
113 repa = AliITSRecoParam::GetHighFluxParam();
114 AliWarning("Using default AliITSRecoParam class");
118 AliITSCalibrationSSD* cal = (AliITSCalibrationSSD*)GetResp(fModule);
121 //---------------------------------------
124 //------------------------------------
125 // fill the digits array with zero-suppression condition
126 // Signal is converted in KeV
129 for (Int_t i=0;i<smaxall; i++){
130 AliITSdigitSSD *d=(AliITSdigitSSD*)alldigits->UncheckedAt(i);
132 if(d->IsSideP()) noise = cal->GetNoiseP(d->GetStripNumber());
133 else noise = cal->GetNoiseN(d->GetStripNumber());
134 if (d->GetSignal()<3.*noise) continue;
136 if(d->IsSideP()) gain = cal->GetGainP(d->GetStripNumber());
137 else gain = cal->GetGainN(d->GetStripNumber());
139 Float_t q=gain*d->GetSignal(); //
140 q=cal->ADCToKeV(q); // converts the charge in KeV from ADC units
141 d->SetSignal(Int_t(q));
145 Int_t smax = digits.GetEntriesFast();
147 //------------------------------------
150 const Int_t kMax=1000;
152 Ali1Dcluster pos[kMax], neg[kMax];
153 Float_t y=0., q=0., qmax=0.;
154 Int_t lab[4]={-2,-2,-2,-2};
158 cout<<"-----------------------------"<<endl;
159 cout<<"this is module "<<fModule;
164 //--------------------------------------------------------
165 // start 1D-clustering from the first digit in the digits array
167 AliITSdigitSSD *d=(AliITSdigitSSD*)digits.UncheckedAt(0);
169 y += d->GetCoord2()*d->GetSignal();
171 lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2);
174 noise = cal->GetNoiseP(d->GetStripNumber());
175 gain = cal->GetGainP(d->GetStripNumber());
178 noise = cal->GetNoiseN(d->GetStripNumber());
179 gain = cal->GetGainN(d->GetStripNumber());
182 noise=cal->ADCToKeV(noise); // converts noise in KeV from ADC units
184 if(qmax>fgkThreshold*noise) flag5=1; // seed for the cluster
187 cout<<d->GetSignal()<<" "<<noise<<" "<<flag5<<" "<<
188 d->GetCoord1()<<" "<<d->GetCoord2()<<endl;
191 Int_t curr=d->GetCoord2();
192 Int_t flag=d->GetCoord1();
194 // Note: the first side which will be processed is supposed to be the
195 // P-side which is neg
198 if(flag) {n=&np; c=pos;} // in case we have only Nstrips (P was bad!)
202 for (Int_t ilab=0;ilab<10;ilab++){
205 milab[0]=d->GetTrack(0); milab[1]=d->GetTrack(1); milab[2]=d->GetTrack(2);
208 //----------------------------------------------------------
209 // search for neighboring digits
211 for (Int_t s=1; s<smax; s++) {
212 d=(AliITSdigitSSD*)digits.UncheckedAt(s);
213 Int_t strip=d->GetCoord2();
215 // if digits is not a neighbour or side did not change
216 // and at least one of the previous digits met the seed condition
217 // then creates a new 1D cluster
218 if ( ( ((strip-curr) > 1) || (flag!=d->GetCoord1()) ) ) {
221 //cout<<"here1"<<endl;
226 c[*n].SetLabels(milab);
228 if(repa->GetUseUnfoldingInClusterFinderSSD()==kTRUE) {
229 // Note: fUseUnfoldingInClusterFinderSSD=kFALSE by default in RecoParam
231 //Split suspiciously big cluster
233 c[*n].SetY(y/q-0.25*nd);
237 Error("FindClustersSSD","Too many 1D clusters !");
240 c[*n].SetY(y/q+0.25*nd);
243 c[*n].SetLabels(milab);
250 Error("FindClustersSSD","Too many 1D clusters !");
260 lab[0]=lab[1]=lab[2]=-2;
261 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
263 // if side changed from P to N, switch to pos 1D clusters
264 // (if for some reason the side changed from N to P then do the opposite)
265 if (flag!=d->GetCoord1())
266 { if(!flag) {n=&np; c=pos;} else {n=&nn; c=neg;} }
268 } // end create new 1D cluster from previous neighboring digits
270 // continues adding digits to the previous cluster
271 // or start a new one
274 y += d->GetCoord2()*d->GetSignal();
278 noise = cal->GetNoiseP(d->GetStripNumber());
279 gain = cal->GetGainP(d->GetStripNumber());
282 noise = cal->GetNoiseN(d->GetStripNumber());
283 gain = cal->GetGainN(d->GetStripNumber());
286 noise=cal->ADCToKeV(noise); // converts the charge in KeV from ADC units
288 if(d->GetSignal()>fgkThreshold*noise) flag5=1;
291 cout<<d->GetSignal()<<" "<<noise<<" "<<flag5<<" "<<
292 d->GetCoord1()<<" "<<d->GetCoord2()<<endl;
295 if (d->GetSignal()>qmax) {
297 lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2);
299 for (Int_t ilab=0;ilab<10;ilab++) {
300 if (d->GetTrack(ilab)>=0) AddLabel(milab, (d->GetTrack(ilab)));
305 } // loop over digits, no more digits in the digits array
308 // add the last 1D cluster
311 // cout<<"here2"<<endl;
316 c[*n].SetLabels(lab);
318 if(repa->GetUseUnfoldingInClusterFinderSSD()==kTRUE) {
320 //Split suspiciously big cluster
322 c[*n].SetY(y/q-0.25*nd);
326 Error("FindClustersSSD","Too many 1D clusters !");
329 c[*n].SetY(y/q+0.25*nd);
332 c[*n].SetLabels(lab);
338 Error("FindClustersSSD","Too many 1D clusters !");
342 } // if flag5 last 1D cluster added
345 //------------------------------------------------------
346 // call FindClustersSSD to pair neg and pos 1D clusters
347 // and create recpoints from the crosses
348 // Note1: neg are Pside and pos are Nside!!
349 // Note2: if there are no Pside digits nn=0 (bad strips!!) (same for Nside)
351 // cout<<nn<<" Pside and "<<np<<" Nside clusters"<<endl;
352 FindClustersSSD(neg, nn, pos, np);
354 //-----------------------------------------------------
359 void AliITSClusterFinderV2SSD::RawdataToClusters(AliRawReader* rawReader,TClonesArray** clusters){
361 //------------------------------------------------------------
362 // This function creates ITS clusters from raw data
363 //------------------------------------------------------------
365 AliITSRawStreamSSD inputSSD(rawReader);
366 FindClustersSSD(&inputSSD,clusters);
370 void AliITSClusterFinderV2SSD::FindClustersSSD(AliITSRawStreamSSD* input,
371 TClonesArray** clusters)
373 //------------------------------------------------------------
374 // Actual SSD cluster finder for raw data
375 //------------------------------------------------------------
377 static AliITSRecoParam *repa = NULL;
379 repa = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
381 repa = AliITSRecoParam::GetHighFluxParam();
382 AliWarning("Using default AliITSRecoParam class");
386 Int_t nClustersSSD = 0;
387 const Int_t kMax = 1000;
388 Ali1Dcluster clusters1D[2][kMax];
389 Int_t nClusters[2] = {0, 0};
390 Int_t lab[3]={-2,-2,-2};
396 // Float_t pedestal=0.;
398 AliITSCalibrationSSD* cal=NULL;
400 Int_t matrix[12][1536];
407 Int_t osignal = 65535;
411 // read raw data input stream
414 // reset signal matrix
415 for(Int_t i=0; i<12; i++) { for(Int_t j=0; j<1536; j++) { matrix[i][j] = 65535;} }
417 if((osignal!=65535)&&(ostrip>0)&&(ostrip<1536)) {
419 matrix[oadc][ostrip] = osignal; // recover data from previous occurence of input->Next()
422 // buffer data for ddl=iddl and ad=iad
425 next = input->Next();
426 if((!next)&&(input->flag)) continue;
427 Int_t ddl=input->GetDDL();
428 Int_t ad=input->GetAD();
429 Int_t adc = input->GetADC(); adc = (adc<6)? adc : adc - 2;
430 Int_t strip = input->GetStrip();
431 if(input->GetSideFlag()) strip=1535-strip;
432 Int_t signal = input->GetSignal();
434 if((ddl==iddl)&&(ad==iad)&&(strip>0)&&(strip<1536)) {n++; matrix[adc][strip] = signal;}
435 else {if ((strip<1536) && (strip>0)) {oddl=iddl; oad=iad; oadc = adc; ostrip = strip; osignal=signal; iddl=ddl; iad=ad; break;}}
437 if(!next) {oddl=iddl; oad=iad; oadc = adc; ostrip = strip; osignal=signal; iddl=ddl; iad=ad; break;}
442 if(!next && oddl<0) break;
444 if(n==0) continue; // first occurence
448 if (repa->GetUseCosmicRunShiftsSSD()) { // Special condition for 2007/2008 cosmic data
449 dStrip = fgkCosmic2008StripShifts[oddl][oad-1];
451 if (TMath::Abs(dStrip) > 1.5)
452 AliError(Form("Indexing error ? oddl = %d, dStrip %f\n",oddl,dStrip));
454 for(Int_t iadc=0; iadc<12; iadc++) { // loop over ADC index for ddl=oddl and ad=oad
456 Int_t iimod = (oad - 1) * 12 + iadc;
457 Int_t iModule = AliITSRawStreamSSD::GetModuleNumber(oddl,iimod);
458 if(iModule==-1) continue;
459 cal = (AliITSCalibrationSSD*)GetResp(iModule);
465 for(Int_t istrip=0; istrip<768; istrip++) { // P-side
466 Int_t signal = matrix[iadc][istrip];
467 pedestal = cal->GetPedestalP(istrip);
468 matrix[iadc][istrip]=signal-(Int_t)pedestal;
474 for(Int_t l=0; l<6; l++) {
476 for(Int_t n=20; n<108; n++) cmode+=matrix[iadc][l*128+n];
478 for(Int_t n=0; n<128; n++) matrix[iadc][l*128+n]-=(Int_t)cmode;
484 for(istrip=0; istrip<768; istrip++) { // P-side
486 Int_t signal = TMath::Abs(matrix[iadc][istrip]);
489 noise = cal->GetNoiseP(istrip); if(noise<1.) signal = 65535;
490 if(signal<3*noise) signal = 65535; // in case ZS was not done in hw do it now
492 // if(cal->IsPChannelBad(istrip)) signal=0;
495 gain = cal->GetGainP(istrip);
496 signal = (Int_t) ( signal * gain ); // signal is corrected for gain
497 if(signal>fgkThreshold*noise) flag5=1;
498 signal = (Int_t) cal->ADCToKeV( signal ); // signal is converted in KeV
500 q += signal; // add digit to current cluster
501 y += istrip * signal;
508 if ( (nDigits>0) && flag5 ) {
510 Ali1Dcluster& cluster = clusters1D[0][nClusters[0]++];
512 if(q!=0) cluster.SetY(y/q + dStrip);
513 else cluster.SetY(istrip + dStrip -1);
516 cluster.SetNd(nDigits);
517 cluster.SetLabels(lab);
519 if(repa->GetUseUnfoldingInClusterFinderSSD()==kTRUE) {
521 //Split suspiciously big cluster
522 if (nDigits > 4&&nDigits < 25) {
523 if(q!=0) cluster.SetY(y/q + dStrip - 0.25*nDigits);
524 else cluster.SetY(istrip-1 + dStrip - 0.25*nDigits);
526 if (nClusters[0] == kMax) {
527 Error("FindClustersSSD", "Too many 1D clusters !");
530 Ali1Dcluster& cluster2 = clusters1D[0][nClusters[0]++];
531 if(q!=0) cluster2.SetY(y/q + dStrip + 0.25*nDigits);
532 else cluster2.SetY(istrip-1 + dStrip + 0.25*nDigits);
533 cluster2.SetQ(0.5*q);
534 cluster2.SetNd(nDigits);
535 cluster2.SetLabels(lab);
546 } // loop over strip on P-side
548 // if last strip does have signal
551 if ( (nDigits>0) && flag5 ) {
553 Ali1Dcluster& cluster = clusters1D[0][nClusters[0]++];
555 if(q!=0) cluster.SetY(y/q + dStrip);
556 else cluster.SetY(istrip - 1 + dStrip);
559 cluster.SetNd(nDigits);
560 cluster.SetLabels(lab);
562 if(repa->GetUseUnfoldingInClusterFinderSSD()==kTRUE) {
564 //Split suspiciously big cluster
565 if (nDigits > 4&&nDigits < 25) {
566 if(q!=0) cluster.SetY(y/q + dStrip - 0.25*nDigits);
567 else cluster.SetY(istrip-1 + dStrip - 0.25*nDigits);
569 if (nClusters[0] == kMax) {
570 Error("FindClustersSSD", "Too many 1D clusters !");
573 Ali1Dcluster& cluster2 = clusters1D[0][nClusters[0]++];
574 if(q!=0) cluster2.SetY(y/q + dStrip + 0.25*nDigits);
575 else cluster2.SetY(istrip-1 + dStrip + 0.25*nDigits);
576 cluster2.SetQ(0.5*q);
577 cluster2.SetNd(nDigits);
578 cluster2.SetLabels(lab);
590 for(Int_t istrip=768; istrip<1536; istrip++) { // P-side
591 Int_t signal = matrix[iadc][istrip];
592 pedestal = cal->GetPedestalN(1535-istrip);
593 matrix[iadc][istrip]=signal-(Int_t)pedestal;
598 for(Int_t l=6; l<12; l++) {
600 for(Int_t n=20; n<108; n++) cmode+=matrix[iadc][l*128+n];
602 for(Int_t n=0; n<128; n++) matrix[iadc][l*128+n]-=(Int_t)cmode;
609 for(Int_t iistrip=768; iistrip<1536; iistrip++) { // N-side
611 Int_t signal = TMath::Abs(matrix[iadc][iistrip]);
612 strip = 1535-iistrip;
615 noise = cal->GetNoiseN(strip); if(noise<1.) signal=65535;
617 // if(cal->IsNChannelBad(strip)) signal=0;
619 if(signal<3*noise) signal = 65535; // in case ZS was not done in hw do it now
622 gain = cal->GetGainN(strip);
623 signal = (Int_t) ( signal * gain); // signal is corrected for gain
624 if(signal>fgkThreshold*noise) flag5=1;
625 signal = (Int_t) cal->ADCToKeV( signal ); // signal is converted in KeV
627 // add digit to current cluster
636 if ( (nDigits>0) && flag5 ) {
638 Ali1Dcluster& cluster = clusters1D[1][nClusters[1]++];
640 if(q!=0) cluster.SetY(y/q - dStrip);
641 else cluster.SetY(strip+1 - dStrip);
644 cluster.SetNd(nDigits);
645 cluster.SetLabels(lab);
647 if(repa->GetUseUnfoldingInClusterFinderSSD()==kTRUE) {
649 //Split suspiciously big cluster
650 if (nDigits > 4&&nDigits < 25) {
651 cluster.SetY(y/q - dStrip - 0.25*nDigits);
653 if (nClusters[1] == kMax) {
654 Error("FindClustersSSD", "Too many 1D clusters !");
657 Ali1Dcluster& cluster2 = clusters1D[1][nClusters[1]++];
658 cluster2.SetY(y/q - dStrip + 0.25*nDigits);
659 cluster2.SetQ(0.5*q);
660 cluster2.SetNd(nDigits);
661 cluster2.SetLabels(lab);
672 } // loop over strips on N-side
676 if ( (nDigits>0) && flag5 ) {
678 Ali1Dcluster& cluster = clusters1D[1][nClusters[1]++];
680 if(q!=0) cluster.SetY(y/q - dStrip);
681 else cluster.SetY(strip - dStrip + 1);
684 cluster.SetNd(nDigits);
685 cluster.SetLabels(lab);
687 if(repa->GetUseUnfoldingInClusterFinderSSD()==kTRUE) {
689 //Split suspiciously big cluster
690 if (nDigits > 4&&nDigits < 25) {
691 if(q!=0) cluster.SetY(y/q - dStrip - 0.25*nDigits);
692 else cluster.SetY(strip+1 - dStrip - 0.25*nDigits);
694 if (nClusters[1] == kMax) {
695 Error("FindClustersSSD", "Too many 1D clusters !");
698 Ali1Dcluster& cluster2 = clusters1D[1][nClusters[1]++];
699 if(q!=0) cluster2.SetY(y/q - dStrip + 0.25*nDigits);
700 else cluster2.SetY(strip+1 - dStrip + 0.25*nDigits);
701 cluster2.SetQ(0.5*q);
702 cluster2.SetNd(nDigits);
703 cluster2.SetLabels(lab);
715 if((nClusters[0])&&(nClusters[1])) {
717 clusters[iModule] = new TClonesArray("AliITSRecPoint");
719 FindClustersSSD(&clusters1D[0][0], nClusters[0],
720 &clusters1D[1][0], nClusters[1], clusters[iModule]);
721 Int_t nClustersn = clusters[iModule]->GetEntriesFast();
722 nClustersSSD += nClustersn;
725 nClusters[0] = nClusters[1] = 0;
734 Info("FindClustersSSD", "found clusters in ITS SSD: %d", nClustersSSD);
737 void AliITSClusterFinderV2SSD::
738 FindClustersSSD(Ali1Dcluster* neg, Int_t nn,
739 Ali1Dcluster* pos, Int_t np,
740 TClonesArray *clusters) {
741 //------------------------------------------------------------
742 // Actual SSD cluster finder
743 //------------------------------------------------------------
745 const TGeoHMatrix *mT2L=AliITSgeomTGeo::GetTracking2LocalMatrix(fModule);
747 TClonesArray &cl=*clusters;
749 AliITSsegmentationSSD *seg = dynamic_cast<AliITSsegmentationSSD*>(fDetTypeRec->GetSegmentationModel(2));
750 if (fModule>fLastSSD1)
755 Float_t hwSSD = seg->Dx()*1e-4/2;
756 Float_t hlSSD = seg->Dz()*1e-4/2;
758 Int_t idet=fNdet[fModule];
762 Int_t *cnegative = new Int_t[np];
763 Int_t *cused1 = new Int_t[np];
764 Int_t *negativepair = new Int_t[10*np];
765 Int_t *cpositive = new Int_t[nn];
766 Int_t *cused2 = new Int_t[nn];
767 Int_t *positivepair = new Int_t[10*nn];
768 for (Int_t i=0;i<np;i++) {cnegative[i]=0; cused1[i]=0;}
769 for (Int_t i=0;i<nn;i++) {cpositive[i]=0; cused2[i]=0;}
770 for (Int_t i=0;i<10*np;i++) {negativepair[i]=0;}
771 for (Int_t i=0;i<10*nn;i++) {positivepair[i]=0;}
773 if ((np*nn) > fgPairsSize) {
775 if (fgPairs) delete [] fgPairs;
776 fgPairsSize = 4*np*nn;
777 fgPairs = new Short_t[fgPairsSize];
779 memset(fgPairs,0,sizeof(Short_t)*np*nn);
782 // find available pairs
784 for (Int_t i=0; i<np; i++) {
785 Float_t yp=pos[i].GetY();
786 if ( (pos[i].GetQ()>0) && (pos[i].GetQ()<3) ) continue;
787 for (Int_t j=0; j<nn; j++) {
788 if ( (neg[j].GetQ()>0) && (neg[j].GetQ()<3) ) continue;
789 Float_t yn=neg[j].GetY();
792 seg->GetPadCxz(yn, yp, xt, zt);
793 //cout<<yn<<" "<<yp<<" "<<xt<<" "<<zt<<endl;
795 if (TMath::Abs(xt)<hwSSD+0.01)
796 if (TMath::Abs(zt)<hlSSD+0.01*(neg[j].GetNd()+pos[i].GetNd())) {
797 negativepair[i*10+cnegative[i]] =j; //index
798 positivepair[j*10+cpositive[j]] =i;
799 cnegative[i]++; //counters
807 // try to recover points out of but close to the module boundaries
809 for (Int_t i=0; i<np; i++) {
810 Float_t yp=pos[i].GetY();
811 if ( (pos[i].GetQ()>0) && (pos[i].GetQ()<3) ) continue;
812 for (Int_t j=0; j<nn; j++) {
813 if ( (neg[j].GetQ()>0) && (neg[j].GetQ()<3) ) continue;
814 // if both 1Dclusters have an other cross continue
815 if (cpositive[j]&&cnegative[i]) continue;
816 Float_t yn=neg[j].GetY();
819 seg->GetPadCxz(yn, yp, xt, zt);
821 if (TMath::Abs(xt)<hwSSD+0.1)
822 if (TMath::Abs(zt)<hlSSD+0.15) {
823 // tag 1Dcluster (eventually will produce low quality recpoint)
824 if (cnegative[i]==0) pos[i].SetNd(100); // not available pair
825 if (cpositive[j]==0) neg[j].SetNd(100); // not available pair
826 negativepair[i*10+cnegative[i]] =j; //index
827 positivepair[j*10+cpositive[j]] =i;
828 cnegative[i]++; //counters
841 static AliITSRecoParam *repa = NULL;
843 repa = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
845 repa = AliITSRecoParam::GetHighFluxParam();
846 AliWarning("Using default AliITSRecoParam class");
850 if(repa->GetUseChargeMatchingInClusterFinderSSD()==kTRUE) {
856 for (Int_t ip=0;ip<np;ip++){
857 Float_t xbest=1000,zbest=1000,qbest=0;
859 // select gold clusters
860 if ( (cnegative[ip]==1) && cpositive[negativepair[10*ip]]==1){
861 Float_t yp=pos[ip].GetY();
862 Int_t j = negativepair[10*ip];
864 if( (pos[ip].GetQ()==0) && (neg[j].GetQ() ==0) ) {
865 // both bad, hence continue;
866 // mark both as used (to avoid recover at the end)
872 ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
873 //cout<<"ratio="<<ratio<<endl;
875 // charge matching (note that if posQ or negQ is 0 -> ratio=1 and the following condition is met
876 if (TMath::Abs(ratio)>0.2) continue; // note: 0.2=3xsigma_ratio calculated in cosmics tests
879 Float_t yn=neg[j].GetY();
882 seg->GetPadCxz(yn, yp, xt, zt);
887 qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
888 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
891 Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
892 mT2L->MasterToLocal(loc,trk);
897 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
898 for (Int_t ilab=0;ilab<3;ilab++){
899 milab[ilab] = pos[ip].GetLabel(ilab);
900 milab[ilab+3] = neg[j].GetLabel(ilab);
904 milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
905 Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
907 lp[2]=0.0022*0.0022; //SigmaY2
908 lp[3]=0.110*0.110; //SigmaZ2
909 // out-of-diagonal element of covariance matrix
910 if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
911 else if ( (info[0]>1) && (info[1]>1) ) {
912 lp[2]=0.0016*0.0016; //SigmaY2
913 lp[3]=0.08*0.08; //SigmaZ2
918 if (info[0]==1) { lp[5]=-0.00014;}
919 else { lp[2]=0.0017*0.0017; lp[5]=-0.00004;}
922 AliITSRecPoint * cl2;
924 if(clusters){ // Note clusters != 0 when method is called for rawdata
927 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
929 cl2->SetChargeRatio(ratio);
933 if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
938 if(pos[ip].GetQ()==0) cl2->SetType(3);
939 if(neg[j].GetQ()==0) cl2->SetType(4);
945 else{ // Note clusters == 0 when method is called for digits
947 cl2 = new AliITSRecPoint(milab,lp,info);
949 cl2->SetChargeRatio(ratio);
953 if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
958 if(pos[ip].GetQ()==0) cl2->SetType(3);
959 if(neg[j].GetQ()==0) cl2->SetType(4);
964 fDetTypeRec->AddRecPoint(*cl2);
970 for (Int_t ip=0;ip<np;ip++){
971 Float_t xbest=1000,zbest=1000,qbest=0;
974 // select "silber" cluster
975 if ( cnegative[ip]==1 && cpositive[negativepair[10*ip]]==2){
976 Int_t in = negativepair[10*ip];
977 Int_t ip2 = positivepair[10*in];
978 if (ip2==ip) ip2 = positivepair[10*in+1];
979 Float_t pcharge = pos[ip].GetQ()+pos[ip2].GetQ();
983 ratio = (pcharge-neg[in].GetQ())/(pcharge+neg[in].GetQ());
984 if ( (TMath::Abs(ratio)<0.2) && (pcharge!=0) ) {
985 //if ( (TMath::Abs(pcharge-neg[in].GetQ())<30) && (pcharge!=0) ) { //
989 if ( (fgPairs[ip*nn+in]==100)&&(pos[ip].GetQ() ) ) { //
991 Float_t yp=pos[ip].GetY();
992 Float_t yn=neg[in].GetY();
995 seg->GetPadCxz(yn, yp, xt, zt);
999 qbest =pos[ip].GetQ();
1000 Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
1001 mT2L->MasterToLocal(loc,trk);
1006 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1007 for (Int_t ilab=0;ilab<3;ilab++){
1008 milab[ilab] = pos[ip].GetLabel(ilab);
1009 milab[ilab+3] = neg[in].GetLabel(ilab);
1012 CheckLabels2(milab);
1013 ratio = (pos[ip].GetQ()-neg[in].GetQ())/(pos[ip].GetQ()+neg[in].GetQ());
1014 milab[3]=(((ip<<10) + in)<<10) + idet; // pos|neg|det
1015 Int_t info[3] = {pos[ip].GetNd(),neg[in].GetNd(),fNlayer[fModule]};
1017 lp[2]=0.0022*0.0022; //SigmaY2
1018 lp[3]=0.110*0.110; //SigmaZ2
1019 // out-of-diagonal element of covariance matrix
1020 if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1021 else if ( (info[0]>1) && (info[1]>1) ) {
1022 lp[2]=0.0016*0.0016; //SigmaY2
1023 lp[3]=0.08*0.08; //SigmaZ2
1027 lp[3]=0.0093*0.0093;
1028 if (info[0]==1) { lp[5]=-0.00014;}
1029 else { lp[2]=0.0017*0.0017; lp[5]=-0.00004;}
1032 AliITSRecPoint * cl2;
1035 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
1036 cl2->SetChargeRatio(ratio);
1038 fgPairs[ip*nn+in] = 5;
1039 if ((pos[ip].GetNd()+neg[in].GetNd())>6){ //multi cluster
1041 fgPairs[ip*nn+in] = 6;
1045 cl2 = new AliITSRecPoint(milab,lp,info);
1046 cl2->SetChargeRatio(ratio);
1048 fgPairs[ip*nn+in] = 5;
1049 if ((pos[ip].GetNd()+neg[in].GetNd())>6){ //multi cluster
1051 fgPairs[ip*nn+in] = 6;
1054 fDetTypeRec->AddRecPoint(*cl2);
1063 // if (!(cused1[ip2] || cused2[in])){ //
1064 if ( (fgPairs[ip2*nn+in]==100) && (pos[ip2].GetQ()) ) {
1066 Float_t yp=pos[ip2].GetY();
1067 Float_t yn=neg[in].GetY();
1070 seg->GetPadCxz(yn, yp, xt, zt);
1074 qbest =pos[ip2].GetQ();
1076 Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
1077 mT2L->MasterToLocal(loc,trk);
1082 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1083 for (Int_t ilab=0;ilab<3;ilab++){
1084 milab[ilab] = pos[ip2].GetLabel(ilab);
1085 milab[ilab+3] = neg[in].GetLabel(ilab);
1088 CheckLabels2(milab);
1089 ratio = (pos[ip2].GetQ()-neg[in].GetQ())/(pos[ip2].GetQ()+neg[in].GetQ());
1090 milab[3]=(((ip2<<10) + in)<<10) + idet; // pos|neg|det
1091 Int_t info[3] = {pos[ip2].GetNd(),neg[in].GetNd(),fNlayer[fModule]};
1093 lp[2]=0.0022*0.0022; //SigmaY2
1094 lp[3]=0.110*0.110; //SigmaZ2
1095 // out-of-diagonal element of covariance matrix
1096 if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1097 else if ( (info[0]>1) && (info[1]>1) ) {
1098 lp[2]=0.0016*0.0016; //SigmaY2
1099 lp[3]=0.08*0.08; //SigmaZ2
1103 lp[3]=0.0093*0.0093;
1104 if (info[0]==1) { lp[5]=-0.00014;}
1105 else { lp[2]=0.0017*0.0017; lp[5]=-0.00004;}
1108 AliITSRecPoint * cl2;
1110 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
1112 cl2->SetChargeRatio(ratio);
1114 fgPairs[ip2*nn+in] =5;
1115 if ((pos[ip2].GetNd()+neg[in].GetNd())>6){ //multi cluster
1117 fgPairs[ip2*nn+in] =6;
1121 cl2 = new AliITSRecPoint(milab,lp,info);
1122 cl2->SetChargeRatio(ratio);
1124 fgPairs[ip2*nn+in] =5;
1125 if ((pos[ip2].GetNd()+neg[in].GetNd())>6){ //multi cluster
1127 fgPairs[ip2*nn+in] =6;
1130 fDetTypeRec->AddRecPoint(*cl2);
1139 } // charge matching condition
1141 } // 2 Pside cross 1 Nside
1142 } // loop over Pside clusters
1147 for (Int_t jn=0;jn<nn;jn++){
1148 if (cused2[jn]) continue;
1149 Float_t xbest=1000,zbest=1000,qbest=0;
1150 // select "silber" cluster
1151 if ( cpositive[jn]==1 && cnegative[positivepair[10*jn]]==2){
1152 Int_t ip = positivepair[10*jn];
1153 Int_t jn2 = negativepair[10*ip];
1154 if (jn2==jn) jn2 = negativepair[10*ip+1];
1155 Float_t pcharge = neg[jn].GetQ()+neg[jn2].GetQ();
1159 ratio = (pcharge-pos[ip].GetQ())/(pcharge+pos[ip].GetQ());
1160 if ( (TMath::Abs(ratio)<0.2) && (pcharge!=0) ) {
1163 if ( (TMath::Abs(pcharge-pos[ip].GetQ())<30) && // charge matching
1164 (pcharge!=0) ) { // reject combinations of bad strips
1170 // if (!(cused1[ip]||cused2[jn])){
1171 if ( (fgPairs[ip*nn+jn]==100) && (neg[jn].GetQ()) ) { //
1173 Float_t yn=neg[jn].GetY();
1174 Float_t yp=pos[ip].GetY();
1177 seg->GetPadCxz(yn, yp, xt, zt);
1181 qbest =neg[jn].GetQ();
1184 Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
1185 mT2L->MasterToLocal(loc,trk);
1191 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1192 for (Int_t ilab=0;ilab<3;ilab++){
1193 milab[ilab] = pos[ip].GetLabel(ilab);
1194 milab[ilab+3] = neg[jn].GetLabel(ilab);
1197 CheckLabels2(milab);
1198 ratio = (pos[ip].GetQ()-neg[jn].GetQ())/(pos[ip].GetQ()+neg[jn].GetQ());
1199 milab[3]=(((ip<<10) + jn)<<10) + idet; // pos|neg|det
1200 Int_t info[3] = {pos[ip].GetNd(),neg[jn].GetNd(),fNlayer[fModule]};
1202 lp[2]=0.0022*0.0022; //SigmaY2
1203 lp[3]=0.110*0.110; //SigmaZ2
1204 // out-of-diagonal element of covariance matrix
1205 if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1206 else if ( (info[0]>1) && (info[1]>1) ) {
1207 lp[2]=0.0016*0.0016; //SigmaY2
1208 lp[3]=0.08*0.08; //SigmaZ2
1212 lp[3]=0.0093*0.0093;
1213 if (info[0]==1) { lp[5]=-0.00014;}
1214 else { lp[2]=0.0017*0.0017; lp[5]=-0.00004;}
1217 AliITSRecPoint * cl2;
1219 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
1221 cl2->SetChargeRatio(ratio);
1223 fgPairs[ip*nn+jn] =7;
1224 if ((pos[ip].GetNd()+neg[jn].GetNd())>6){ //multi cluster
1226 fgPairs[ip*nn+jn]=8;
1231 cl2 = new AliITSRecPoint(milab,lp,info);
1232 cl2->SetChargeRatio(ratio);
1234 fgPairs[ip*nn+jn] =7;
1235 if ((pos[ip].GetNd()+neg[jn].GetNd())>6){ //multi cluster
1237 fgPairs[ip*nn+jn]=8;
1240 fDetTypeRec->AddRecPoint(*cl2);
1246 // if (!(cused1[ip]||cused2[jn2])){
1247 if ( (fgPairs[ip*nn+jn2]==100)&&(neg[jn2].GetQ() ) ) { //
1249 Float_t yn=neg[jn2].GetY();
1250 Double_t yp=pos[ip].GetY();
1253 seg->GetPadCxz(yn, yp, xt, zt);
1257 qbest =neg[jn2].GetQ();
1260 Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
1261 mT2L->MasterToLocal(loc,trk);
1267 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1268 for (Int_t ilab=0;ilab<3;ilab++){
1269 milab[ilab] = pos[ip].GetLabel(ilab);
1270 milab[ilab+3] = neg[jn2].GetLabel(ilab);
1273 CheckLabels2(milab);
1274 ratio = (pos[ip].GetQ()-neg[jn2].GetQ())/(pos[ip].GetQ()+neg[jn2].GetQ());
1275 milab[3]=(((ip<<10) + jn2)<<10) + idet; // pos|neg|det
1276 Int_t info[3] = {pos[ip].GetNd(),neg[jn2].GetNd(),fNlayer[fModule]};
1278 lp[2]=0.0022*0.0022; //SigmaY2
1279 lp[3]=0.110*0.110; //SigmaZ2
1280 // out-of-diagonal element of covariance matrix
1281 if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1282 else if ( (info[0]>1) && (info[1]>1) ) {
1283 lp[2]=0.0016*0.0016; //SigmaY2
1284 lp[3]=0.08*0.08; //SigmaZ2
1288 lp[3]=0.0093*0.0093;
1289 if (info[0]==1) { lp[5]=-0.00014;}
1290 else { lp[2]=0.0017*0.0017; lp[5]=-0.00004;}
1293 AliITSRecPoint * cl2;
1295 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
1298 cl2->SetChargeRatio(ratio);
1299 fgPairs[ip*nn+jn2]=7;
1301 if ((pos[ip].GetNd()+neg[jn2].GetNd())>6){ //multi cluster
1303 fgPairs[ip*nn+jn2]=8;
1308 cl2 = new AliITSRecPoint(milab,lp,info);
1309 cl2->SetChargeRatio(ratio);
1310 fgPairs[ip*nn+jn2]=7;
1312 if ((pos[ip].GetNd()+neg[jn2].GetNd())>6){ //multi cluster
1314 fgPairs[ip*nn+jn2]=8;
1317 fDetTypeRec->AddRecPoint(*cl2);
1326 } // charge matching condition
1328 } // 2 Nside cross 1 Pside
1329 } // loop over Pside clusters
1333 for (Int_t ip=0;ip<np;ip++){
1335 if(cused1[ip]) continue;
1338 Float_t xbest=1000,zbest=1000,qbest=0;
1342 if ( (cnegative[ip]==2) && cpositive[negativepair[10*ip]]==2){
1343 Float_t minchargediff =4.;
1344 Float_t minchargeratio =0.2;
1347 for (Int_t di=0;di<cnegative[ip];di++){
1348 Int_t jc = negativepair[ip*10+di];
1349 Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
1350 ratio = (pos[ip].GetQ()-neg[jc].GetQ())/(pos[ip].GetQ()+neg[jc].GetQ());
1351 //if (TMath::Abs(chargedif)<minchargediff){
1352 if (TMath::Abs(ratio)<0.2){
1354 minchargediff = TMath::Abs(chargedif);
1355 minchargeratio = TMath::Abs(ratio);
1358 if (j<0) continue; // not proper cluster
1362 for (Int_t di=0;di<cnegative[ip];di++){
1363 Int_t jc = negativepair[ip*10+di];
1364 Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
1365 if (TMath::Abs(chargedif)<minchargediff+3.) count++;
1367 if (count>1) continue; // more than one "proper" cluster for positive
1371 for (Int_t dj=0;dj<cpositive[j];dj++){
1372 Int_t ic = positivepair[j*10+dj];
1373 Float_t chargedif = pos[ic].GetQ()-neg[j].GetQ();
1374 if (TMath::Abs(chargedif)<minchargediff+3.) count++;
1376 if (count>1) continue; // more than one "proper" cluster for negative
1381 for (Int_t dj=0;dj<cnegative[jp];dj++){
1382 Int_t ic = positivepair[jp*10+dj];
1383 Float_t chargedif = pos[ic].GetQ()-neg[jp].GetQ();
1384 if (TMath::Abs(chargedif)<minchargediff+4.) count++;
1386 if (count>1) continue;
1387 if (fgPairs[ip*nn+j]<100) continue;
1392 //almost gold clusters
1393 Float_t yp=pos[ip].GetY();
1394 Float_t yn=neg[j].GetY();
1396 seg->GetPadCxz(yn, yp, xt, zt);
1398 qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
1400 Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
1401 mT2L->MasterToLocal(loc,trk);
1406 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1407 for (Int_t ilab=0;ilab<3;ilab++){
1408 milab[ilab] = pos[ip].GetLabel(ilab);
1409 milab[ilab+3] = neg[j].GetLabel(ilab);
1412 CheckLabels2(milab);
1413 if ((neg[j].GetQ()==0)&&(pos[ip].GetQ()==0)) continue; // reject crosses of bad strips!!
1414 ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
1415 milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
1416 Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
1418 lp[2]=0.0022*0.0022; //SigmaY2
1419 lp[3]=0.110*0.110; //SigmaZ2
1420 // out-of-diagonal element of covariance matrix
1421 if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1422 else if ( (info[0]>1) && (info[1]>1) ) {
1423 lp[2]=0.0016*0.0016; //SigmaY2
1424 lp[3]=0.08*0.08; //SigmaZ2
1428 lp[3]=0.0093*0.0093;
1429 if (info[0]==1) { lp[5]=-0.00014;}
1430 else { lp[2]=0.0017*0.0017; lp[5]=-0.00004;}
1433 AliITSRecPoint * cl2;
1435 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
1437 cl2->SetChargeRatio(ratio);
1439 fgPairs[ip*nn+j]=10;
1440 if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
1442 fgPairs[ip*nn+j]=11;
1448 cl2 = new AliITSRecPoint(milab,lp,info);
1449 cl2->SetChargeRatio(ratio);
1451 fgPairs[ip*nn+j]=10;
1452 if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
1454 fgPairs[ip*nn+j]=11;
1459 fDetTypeRec->AddRecPoint(*cl2);
1464 } // loop over Pside 1Dclusters
1468 for (Int_t ip=0;ip<np;ip++){
1470 if(cused1[ip]) continue;
1473 Float_t xbest=1000,zbest=1000,qbest=0;
1475 // manyxmany clusters
1477 if ( (cnegative[ip]<5) && cpositive[negativepair[10*ip]]<5){
1478 Float_t minchargediff =4.;
1480 for (Int_t di=0;di<cnegative[ip];di++){
1481 Int_t jc = negativepair[ip*10+di];
1482 Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
1483 if (TMath::Abs(chargedif)<minchargediff){
1485 minchargediff = TMath::Abs(chargedif);
1488 if (j<0) continue; // not proper cluster
1491 for (Int_t di=0;di<cnegative[ip];di++){
1492 Int_t jc = negativepair[ip*10+di];
1493 Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
1494 if (TMath::Abs(chargedif)<minchargediff+3.) count++;
1496 if (count>1) continue; // more than one "proper" cluster for positive
1500 for (Int_t dj=0;dj<cpositive[j];dj++){
1501 Int_t ic = positivepair[j*10+dj];
1502 Float_t chargedif = pos[ic].GetQ()-neg[j].GetQ();
1503 if (TMath::Abs(chargedif)<minchargediff+3.) count++;
1505 if (count>1) continue; // more than one "proper" cluster for negative
1510 for (Int_t dj=0;dj<cnegative[jp];dj++){
1511 Int_t ic = positivepair[jp*10+dj];
1512 Float_t chargedif = pos[ic].GetQ()-neg[jp].GetQ();
1513 if (TMath::Abs(chargedif)<minchargediff+4.) count++;
1515 if (count>1) continue;
1516 if (fgPairs[ip*nn+j]<100) continue;
1519 //almost gold clusters
1520 Float_t yp=pos[ip].GetY();
1521 Float_t yn=neg[j].GetY();
1525 seg->GetPadCxz(yn, yp, xt, zt);
1529 qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
1532 Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
1533 mT2L->MasterToLocal(loc,trk);
1538 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1539 for (Int_t ilab=0;ilab<3;ilab++){
1540 milab[ilab] = pos[ip].GetLabel(ilab);
1541 milab[ilab+3] = neg[j].GetLabel(ilab);
1544 CheckLabels2(milab);
1545 if ((neg[j].GetQ()==0)&&(pos[ip].GetQ()==0)) continue; // reject crosses of bad strips!!
1546 ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
1547 milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
1548 Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
1550 lp[2]=0.0022*0.0022; //SigmaY2
1551 lp[3]=0.110*0.110; //SigmaZ2
1552 // out-of-diagonal element of covariance matrix
1553 if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1554 else if ( (info[0]>1) && (info[1]>1) ) {
1555 lp[2]=0.0016*0.0016; //SigmaY2
1556 lp[3]=0.08*0.08; //SigmaZ2
1560 lp[3]=0.0093*0.0093;
1561 if (info[0]==1) { lp[5]=-0.00014;}
1562 else { lp[2]=0.0017*0.0017; lp[5]=-0.00004;}
1565 AliITSRecPoint * cl2;
1567 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
1569 cl2->SetChargeRatio(ratio);
1571 fgPairs[ip*nn+j]=12;
1572 if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
1574 fgPairs[ip*nn+j]=13;
1580 cl2 = new AliITSRecPoint(milab,lp,info);
1581 cl2->SetChargeRatio(ratio);
1583 fgPairs[ip*nn+j]=12;
1584 if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
1586 fgPairs[ip*nn+j]=13;
1591 fDetTypeRec->AddRecPoint(*cl2);
1596 } // loop over Pside 1Dclusters
1598 } // use charge matching
1601 // recover all the other crosses
1603 for (Int_t i=0; i<np; i++) {
1604 Float_t xbest=1000,zbest=1000,qbest=0;
1605 Float_t yp=pos[i].GetY();
1606 if ((pos[i].GetQ()>0)&&(pos[i].GetQ()<3)) continue;
1607 for (Int_t j=0; j<nn; j++) {
1608 // for (Int_t di = 0;di<cpositive[i];di++){
1609 // Int_t j = negativepair[10*i+di];
1610 if ((neg[j].GetQ()>0)&&(neg[j].GetQ()<3)) continue;
1612 if ((neg[j].GetQ()==0)&&(pos[i].GetQ()==0)) continue; // reject crosses of bad strips!!
1614 if (cused2[j]||cused1[i]) continue;
1615 if (fgPairs[i*nn+j]>0 &&fgPairs[i*nn+j]<100) continue;
1616 ratio = (pos[i].GetQ()-neg[j].GetQ())/(pos[i].GetQ()+neg[j].GetQ());
1617 Float_t yn=neg[j].GetY();
1620 seg->GetPadCxz(yn, yp, xt, zt);
1622 if (TMath::Abs(xt)<hwSSD+0.01)
1623 if (TMath::Abs(zt)<hlSSD+0.01*(neg[j].GetNd()+pos[i].GetNd())) {
1626 qbest=0.5*(pos[i].GetQ()+neg[j].GetQ());
1629 Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
1630 mT2L->MasterToLocal(loc,trk);
1635 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1636 for (Int_t ilab=0;ilab<3;ilab++){
1637 milab[ilab] = pos[i].GetLabel(ilab);
1638 milab[ilab+3] = neg[j].GetLabel(ilab);
1641 CheckLabels2(milab);
1642 milab[3]=(((i<<10) + j)<<10) + idet; // pos|neg|det
1643 Int_t info[3] = {pos[i].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
1645 lp[2]=0.0022*0.0022; //SigmaY2
1646 lp[3]=0.110*0.110; //SigmaZ2
1647 // out-of-diagonal element of covariance matrix
1648 if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1649 else if ( (info[0]>1) && (info[1]>1) ) {
1650 lp[2]=0.0016*0.0016; //SigmaY2
1651 lp[3]=0.08*0.08; //SigmaZ2
1655 lp[3]=0.0093*0.0093;
1656 if (info[0]==1) { lp[5]=-0.00014;}
1657 else { lp[2]=0.0017*0.0017; lp[5]=-0.00004;}
1660 AliITSRecPoint * cl2;
1662 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
1664 cl2->SetChargeRatio(ratio);
1665 cl2->SetType(100+cpositive[j]+cnegative[i]);
1667 if(pos[i].GetQ()==0) cl2->SetType(200+cpositive[j]+cnegative[i]);
1668 if(neg[j].GetQ()==0) cl2->SetType(300+cpositive[j]+cnegative[i]);
1672 cl2 = new AliITSRecPoint(milab,lp,info);
1673 cl2->SetChargeRatio(ratio);
1674 cl2->SetType(100+cpositive[j]+cnegative[i]);
1676 if(pos[i].GetQ()==0) cl2->SetType(200+cpositive[j]+cnegative[i]);
1677 if(neg[j].GetQ()==0) cl2->SetType(300+cpositive[j]+cnegative[i]);
1679 fDetTypeRec->AddRecPoint(*cl2);
1687 //---------------------------------------------------------
1688 // recover crosses of good 1D clusters with bad strips on the other side
1689 // Note1: at first iteration skip modules with a bad side (or almost), (would produce too many fake!)
1690 // Note2: for modules with a bad side see below
1692 AliITSCalibrationSSD* cal = (AliITSCalibrationSSD*)GetResp(fModule);
1693 Int_t countPbad=0, countNbad=0;
1694 for(Int_t ib=0; ib<768; ib++) {
1695 if(cal->IsPChannelBad(ib)) countPbad++;
1696 if(cal->IsNChannelBad(ib)) countNbad++;
1698 // AliInfo(Form("module %d has %d P- and %d N-bad strips",fModule,countPbad,countNbad));
1700 if( (countPbad<100) && (countNbad<100) ) { // no bad side!!
1702 for (Int_t i=0; i<np; i++) { // loop over Nside 1Dclusters with no crosses
1703 if(cnegative[i]) continue; // if intersecting Pside clusters continue;
1705 // for(Int_t ib=0; ib<768; ib++) { // loop over all Pstrips
1706 for(Int_t ib=15; ib<753; ib++) { // loop over all Pstrips
1708 if(cal->IsPChannelBad(ib)) { // check if strips is bad
1709 Float_t yN=pos[i].GetY();
1711 seg->GetPadCxz(1.*ib, yN, xt, zt);
1714 // bad Pstrip is crossing the Nside 1Dcluster -> create recpoint
1716 if ( (TMath::Abs(xt)<hwSSD+0.01) && (TMath::Abs(zt)<hlSSD+0.01) ) {
1717 Double_t loc[3]={xt,0.,zt},trk[3]={0.,0.,0.};
1718 mT2L->MasterToLocal(loc,trk);
1721 lp[4]=pos[i].GetQ(); //Q
1722 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1723 for (Int_t ilab=0;ilab<3;ilab++) milab[ilab] = pos[i].GetLabel(ilab);
1724 CheckLabels2(milab);
1725 milab[3]=( (i<<10) << 10 ) + idet; // pos|neg|det
1726 Int_t info[3] = {pos[i].GetNd(),0,fNlayer[fModule]};
1728 // out-of-diagonal element of covariance matrix
1729 if (info[0]==1) lp[5]=0.0065;
1732 lp[2]=0.0022*0.0022; //SigmaY2
1733 lp[3]=0.110*0.110; //SigmaZ2
1734 lp[5]=-0.00012; // out-of-diagonal element of covariance matrix
1736 AliITSRecPoint * cl2;
1738 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
1739 cl2->SetChargeRatio(1.);
1743 cl2 = new AliITSRecPoint(milab,lp,info);
1744 cl2->SetChargeRatio(1.);
1746 fDetTypeRec->AddRecPoint(*cl2);
1749 } // cross is within the detector
1755 } // end loop over Pstrips
1757 } // end loop over Nside 1D clusters
1759 for (Int_t j=0; j<nn; j++) { // loop over Pside 1D clusters with no crosses
1760 if(cpositive[j]) continue;
1762 // for(Int_t ib=0; ib<768; ib++) { // loop over all Nside strips
1763 for(Int_t ib=15; ib<753; ib++) { // loop over all Nside strips
1765 if(cal->IsNChannelBad(ib)) { // check if strip is bad
1766 Float_t yP=neg[j].GetY();
1768 seg->GetPadCxz(yP, 1.*ib, xt, zt);
1771 // bad Nstrip is crossing the Pside 1Dcluster -> create recpoint
1773 if ( (TMath::Abs(xt)<hwSSD+0.01) && (TMath::Abs(zt)<hlSSD+0.01) ) {
1774 Double_t loc[3]={xt,0.,zt},trk[3]={0.,0.,0.};
1775 mT2L->MasterToLocal(loc,trk);
1778 lp[4]=neg[j].GetQ(); //Q
1779 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1780 for (Int_t ilab=0;ilab<3;ilab++) milab[ilab] = neg[j].GetLabel(ilab);
1781 CheckLabels2(milab);
1782 milab[3]=( j << 10 ) + idet; // pos|neg|det
1783 Int_t info[3] = {0.,neg[j].GetNd(),fNlayer[fModule]};
1785 lp[2]=0.0022*0.0022; //SigmaY2
1786 lp[3]=0.110*0.110; //SigmaZ2
1787 lp[5]=-0.00012; // out-of-diagonal element of covariance matrix
1789 AliITSRecPoint * cl2;
1791 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
1792 cl2->SetChargeRatio(1.);
1796 cl2 = new AliITSRecPoint(milab,lp,info);
1797 cl2->SetChargeRatio(1.);
1799 fDetTypeRec->AddRecPoint(*cl2);
1802 } // cross is within the detector
1807 } // end loop over Nstrips
1808 } // end loop over Pside 1D clusters
1812 //---------------------------------------------------------
1814 else if( (countPbad>700) && (countNbad<100) ) { // bad Pside!!
1816 for (Int_t i=0; i<np; i++) { // loop over Nside 1Dclusters with no crosses
1817 if(cnegative[i]) continue; // if intersecting Pside clusters continue;
1820 Float_t yN=pos[i].GetY();
1822 if (seg->GetLayer()==5) yP = yN + (7.6/1.9);
1823 else yP = yN - (7.6/1.9);
1824 seg->GetPadCxz(yP, yN, xt, zt);
1826 if ( (TMath::Abs(xt)<hwSSD+0.01) && (TMath::Abs(zt)<hlSSD+0.01) ) {
1827 Double_t loc[3]={xt,0.,zt},trk[3]={0.,0.,0.};
1828 mT2L->MasterToLocal(loc,trk);
1831 lp[4]=pos[i].GetQ(); //Q
1832 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1833 for (Int_t ilab=0;ilab<3;ilab++) milab[ilab] = pos[i].GetLabel(ilab);
1834 CheckLabels2(milab);
1835 milab[3]=( (i<<10) << 10 ) + idet; // pos|neg|det
1836 Int_t info[3] = {pos[i].GetNd(),0.,fNlayer[fModule]};
1838 lp[2]=0.031*0.031; //SigmaY2
1839 lp[3]=1.15*1.15; //SigmaZ2
1842 AliITSRecPoint * cl2;
1844 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
1845 cl2->SetChargeRatio(1.);
1849 cl2 = new AliITSRecPoint(milab,lp,info);
1850 cl2->SetChargeRatio(1.);
1852 fDetTypeRec->AddRecPoint(*cl2);
1855 } // cross is within the detector
1859 } // end loop over Nside 1D clusters
1861 } // bad Pside module
1863 else if( (countNbad>700) && (countPbad<100) ) { // bad Nside!!
1865 for (Int_t j=0; j<nn; j++) { // loop over Pside 1D clusters with no crosses
1866 if(cpositive[j]) continue;
1869 Float_t yP=neg[j].GetY();
1871 if (seg->GetLayer()==5) yN = yP - (7.6/1.9);
1872 else yN = yP + (7.6/1.9);
1873 seg->GetPadCxz(yP, yN, xt, zt);
1875 if ( (TMath::Abs(xt)<hwSSD+0.01) && (TMath::Abs(zt)<hlSSD+0.01) ) {
1876 Double_t loc[3]={xt,0.,zt},trk[3]={0.,0.,0.};
1877 mT2L->MasterToLocal(loc,trk);
1880 lp[4]=neg[j].GetQ(); //Q
1881 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1882 for (Int_t ilab=0;ilab<3;ilab++) milab[ilab] = neg[j].GetLabel(ilab);
1883 CheckLabels2(milab);
1884 milab[3]=( j << 10 ) + idet; // pos|neg|det
1885 Int_t info[3] = {0.,neg[j].GetNd(),fNlayer[fModule]};
1887 lp[2]=0.0085*0.0085; //SigmaY2
1888 lp[3]=1.15*1.15; //SigmaZ2
1891 AliITSRecPoint * cl2;
1893 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
1894 cl2->SetChargeRatio(1.);
1898 cl2 = new AliITSRecPoint(milab,lp,info);
1899 cl2->SetChargeRatio(1.);
1901 fDetTypeRec->AddRecPoint(*cl2);
1904 } // cross is within the detector
1908 } // end loop over Pside 1D clusters
1910 } // bad Nside module
1912 //---------------------------------------------------------
1914 //cout<<ncl<<" clusters for this module"<<endl;
1916 delete [] cnegative;
1918 delete [] negativepair;
1919 delete [] cpositive;
1921 delete [] positivepair;