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 // Revised: Enrico Fragiacomo, enrico.fragiacomo@ts.infn.it //
23 // Revised 23/06/08: Marco Bregant
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;
45 ClassImp(AliITSClusterFinderV2SSD)
48 AliITSClusterFinderV2SSD::AliITSClusterFinderV2SSD(AliITSDetTypeRec* dettyp):AliITSClusterFinderV2(dettyp),
49 fLastSSD1(AliITSgeomTGeo::GetModuleIndex(6,1,1)-1)
55 //______________________________________________________________________
56 AliITSClusterFinderV2SSD::AliITSClusterFinderV2SSD(const AliITSClusterFinderV2SSD &cf) : AliITSClusterFinderV2(cf), fLastSSD1(cf.fLastSSD1)
61 //______________________________________________________________________
62 AliITSClusterFinderV2SSD& AliITSClusterFinderV2SSD::operator=(const AliITSClusterFinderV2SSD& cf ){
63 // Assignment operator
65 this->~AliITSClusterFinderV2SSD();
66 new(this) AliITSClusterFinderV2SSD(cf);
71 void AliITSClusterFinderV2SSD::FindRawClusters(Int_t mod){
75 FindClustersSSD(fDigits);
79 void AliITSClusterFinderV2SSD::FindClustersSSD(TClonesArray *alldigits) {
80 //------------------------------------------------------------
81 // Actual SSD cluster finder
82 //------------------------------------------------------------
84 static AliITSRecoParam *repa = NULL;
88 repa = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
90 repa = AliITSRecoParam::GetHighFluxParam();
91 AliWarning("Using default AliITSRecoParam class");
95 AliITSCalibrationSSD* cal = (AliITSCalibrationSSD*)GetResp(fModule);
98 Int_t smaxall=alldigits->GetEntriesFast();
99 if (smaxall==0) return;
100 // TObjArray *digits = new TObjArray;
102 for (Int_t i=0;i<smaxall; i++){
103 AliITSdigitSSD *d=(AliITSdigitSSD*)alldigits->UncheckedAt(i);
105 if(d->IsSideP()) gain = cal->GetGainP(d->GetStripNumber());
106 else gain = cal->GetGainN(d->GetStripNumber());
108 Float_t q=gain*d->GetSignal(); // calibration brings mip peaks around 120 (in ADC units)
109 q=cal->ADCToKeV(q); // converts the charge in KeV from ADC units
110 d->SetSignal(Int_t(q));
112 if (d->GetSignal()<3) continue;
115 Int_t smax = digits.GetEntriesFast();
118 const Int_t kMax=1000;
120 Ali1Dcluster pos[kMax], neg[kMax];
121 Float_t y=0., q=0., qmax=0.;
122 Int_t lab[4]={-2,-2,-2,-2};
124 AliITSdigitSSD *d=(AliITSdigitSSD*)digits.UncheckedAt(0);
126 y += d->GetCoord2()*d->GetSignal();
128 lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2);
129 Int_t curr=d->GetCoord2();
130 Int_t flag=d->GetCoord1();
135 for (Int_t ilab=0;ilab<10;ilab++){
138 milab[0]=d->GetTrack(0); milab[1]=d->GetTrack(1); milab[2]=d->GetTrack(2);
140 for (Int_t s=1; s<smax; s++) {
141 d=(AliITSdigitSSD*)digits.UncheckedAt(s);
142 Int_t strip=d->GetCoord2();
143 if ((strip-curr) > 1 || flag!=d->GetCoord1()) {
148 c[*n].SetLabels(milab);
150 if(repa->GetUseUnfoldingInClusterFinderSSD()==kTRUE) {
152 //Split suspiciously big cluster
154 c[*n].SetY(y/q-0.25*nd);
158 Error("FindClustersSSD","Too many 1D clusters !");
161 c[*n].SetY(y/q+0.25*nd);
164 c[*n].SetLabels(milab);
171 Error("FindClustersSSD","Too many 1D clusters !");
176 lab[0]=lab[1]=lab[2]=-2;
178 for (Int_t ilab=0;ilab<10;ilab++){
182 if (flag!=d->GetCoord1()) { n=&np; c=pos; }
186 y += d->GetCoord2()*d->GetSignal();
188 if (d->GetSignal()>qmax) {
190 lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2);
192 for (Int_t ilab=0;ilab<10;ilab++) {
193 if (d->GetTrack(ilab)>=0) AddLabel(milab, (d->GetTrack(ilab)));
200 c[*n].SetLabels(lab);
202 if(repa->GetUseUnfoldingInClusterFinderSSD()==kTRUE) {
204 //Split suspiciously big cluster
206 c[*n].SetY(y/q-0.25*nd);
210 Error("FindClustersSSD","Too many 1D clusters !");
213 c[*n].SetY(y/q+0.25*nd);
216 c[*n].SetLabels(lab);
222 Error("FindClustersSSD","Too many 1D clusters !");
226 FindClustersSSD(neg, nn, pos, np);
230 void AliITSClusterFinderV2SSD::RawdataToClusters(AliRawReader* rawReader,TClonesArray** clusters){
232 //------------------------------------------------------------
233 // This function creates ITS clusters from raw data
234 //------------------------------------------------------------
236 AliITSRawStreamSSD inputSSD(rawReader);
237 FindClustersSSD(&inputSSD,clusters);
241 void AliITSClusterFinderV2SSD::FindClustersSSD(AliITSRawStreamSSD* input,
242 TClonesArray** clusters)
244 //------------------------------------------------------------
245 // Actual SSD cluster finder for raw data
246 //------------------------------------------------------------
248 static AliITSRecoParam *repa = NULL;
250 repa = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
252 repa = AliITSRecoParam::GetHighFluxParam();
253 AliWarning("Using default AliITSRecoParam class");
257 Int_t nClustersSSD = 0;
258 const Int_t kMax = 1000;
259 Ali1Dcluster clusters1D[2][kMax];
260 Int_t nClusters[2] = {0, 0};
261 Int_t lab[3]={-2,-2,-2};
267 // Float_t pedestal=0.;
269 AliITSCalibrationSSD* cal=NULL;
271 Int_t matrix[12][1536];
278 Int_t osignal = 65535;
282 // read raw data input stream
285 // reset signal matrix
286 for(Int_t i=0; i<12; i++) { for(Int_t j=0; j<1536; j++) { matrix[i][j] = 65535;} }
288 if((osignal!=65535)&&(ostrip>0)&&(ostrip<1536)) {
290 matrix[oadc][ostrip] = osignal; // recover data from previous occurence of input->Next()
293 // buffer data for ddl=iddl and ad=iad
296 next = input->Next();
297 if((!next)&&(input->flag)) continue;
298 Int_t ddl=input->GetDDL();
299 Int_t ad=input->GetAD();
300 Int_t adc = input->GetADC(); adc = (adc<6)? adc : adc - 2;
301 Int_t strip = input->GetStrip();
302 if(input->GetSideFlag()) strip=1535-strip;
303 Int_t signal = input->GetSignal();
305 if((ddl==iddl)&&(ad==iad)&&(strip>0)&&(strip<1536)) {n++; matrix[adc][strip] = signal;}
306 else {oddl=iddl; oad=iad; oadc = adc; ostrip = strip; osignal=signal; iddl=ddl; iad=ad; break;}
308 if(!next) {oddl=iddl; oad=iad; oadc = adc; ostrip = strip; osignal=signal; iddl=ddl; iad=ad; break;}
313 if(!next && oddl<0) break;
315 if(n==0) continue; // first occurence
319 for(Int_t iadc=0; iadc<12; iadc++) { // loop over ADC index for ddl=oddl and ad=oad
321 Int_t iimod = (oad - 1) * 12 + iadc;
322 Int_t iModule = AliITSRawStreamSSD::GetModuleNumber(oddl,iimod);
323 if(iModule==-1) continue;
324 cal = (AliITSCalibrationSSD*)GetResp(iModule);
329 for(Int_t istrip=0; istrip<768; istrip++) { // P-side
330 Int_t signal = matrix[iadc][istrip];
331 pedestal = cal->GetPedestalP(istrip);
332 matrix[iadc][istrip]=signal-(Int_t)pedestal;
338 for(Int_t l=0; l<6; l++) {
340 for(Int_t n=20; n<108; n++) cmode+=matrix[iadc][l*128+n];
342 for(Int_t n=0; n<128; n++) matrix[iadc][l*128+n]-=(Int_t)cmode;
348 for(istrip=0; istrip<768; istrip++) { // P-side
350 Int_t signal = TMath::Abs(matrix[iadc][istrip]);
353 noise = cal->GetNoiseP(istrip); if(noise<1.) signal = 65535;
354 if(signal<3*noise) signal = 65535; // in case ZS was not done in hw do it now
356 // if(cal->IsPChannelBad(istrip)) signal=0;
359 gain = cal->GetGainP(istrip);
360 signal = (Int_t) ( signal * gain ); // signal is corrected for gain
361 signal = (Int_t) cal->ADCToKeV( signal ); // signal is converted in KeV
363 q += signal; // add digit to current cluster
364 y += istrip * signal;
371 if ( ( (nDigits==1) && ( (q==0) || (q>5*oldnoise)) ) || (nDigits>1) ) {
373 Ali1Dcluster& cluster = clusters1D[0][nClusters[0]++];
375 if(q!=0) cluster.SetY(y/q);
376 else cluster.SetY(istrip-1);
379 cluster.SetNd(nDigits);
380 cluster.SetLabels(lab);
382 if(repa->GetUseUnfoldingInClusterFinderSSD()==kTRUE) {
384 //Split suspiciously big cluster
385 if (nDigits > 4&&nDigits < 25) {
386 if(q!=0) cluster.SetY(y/q - 0.25*nDigits);
387 else cluster.SetY(istrip-1 - 0.25*nDigits);
389 if (nClusters[0] == kMax) {
390 Error("FindClustersSSD", "Too many 1D clusters !");
393 Ali1Dcluster& cluster2 = clusters1D[0][nClusters[0]++];
394 if(q!=0) cluster2.SetY(y/q + 0.25*nDigits);
395 else cluster2.SetY(istrip-1 + 0.25*nDigits);
396 cluster2.SetQ(0.5*q);
397 cluster2.SetNd(nDigits);
398 cluster2.SetLabels(lab);
408 } // loop over strip on P-side
410 // if last strip does have signal
413 if ( ( (nDigits==1) && ( (q==0) || (q>5*oldnoise)) ) || (nDigits>1) ) {
415 Ali1Dcluster& cluster = clusters1D[0][nClusters[0]++];
417 if(q!=0) cluster.SetY(y/q);
418 else cluster.SetY(istrip-1);
421 cluster.SetNd(nDigits);
422 cluster.SetLabels(lab);
424 if(repa->GetUseUnfoldingInClusterFinderSSD()==kTRUE) {
426 //Split suspiciously big cluster
427 if (nDigits > 4&&nDigits < 25) {
428 if(q!=0) cluster.SetY(y/q - 0.25*nDigits);
429 else cluster.SetY(istrip-1 - 0.25*nDigits);
431 if (nClusters[0] == kMax) {
432 Error("FindClustersSSD", "Too many 1D clusters !");
435 Ali1Dcluster& cluster2 = clusters1D[0][nClusters[0]++];
436 if(q!=0) cluster2.SetY(y/q + 0.25*nDigits);
437 else cluster2.SetY(istrip-1 + 0.25*nDigits);
438 cluster2.SetQ(0.5*q);
439 cluster2.SetNd(nDigits);
440 cluster2.SetLabels(lab);
451 for(Int_t istrip=768; istrip<1536; istrip++) { // P-side
452 Int_t signal = matrix[iadc][istrip];
453 pedestal = cal->GetPedestalN(1535-istrip);
454 matrix[iadc][istrip]=signal-(Int_t)pedestal;
459 for(Int_t l=6; l<12; l++) {
461 for(Int_t n=20; n<108; n++) cmode+=matrix[iadc][l*128+n];
463 for(Int_t n=0; n<128; n++) matrix[iadc][l*128+n]-=(Int_t)cmode;
470 for(Int_t iistrip=768; iistrip<1536; iistrip++) { // N-side
472 Int_t signal = TMath::Abs(matrix[iadc][iistrip]);
473 strip = 1535-iistrip;
476 noise = cal->GetNoiseN(strip); if(noise<1.) signal=65535;
478 // if(cal->IsNChannelBad(strip)) signal=0;
480 if(signal<3*noise) signal = 65535; // in case ZS was not done in hw do it now
483 gain = cal->GetGainN(strip);
484 signal = (Int_t) ( signal * gain); // signal is corrected for gain
485 signal = (Int_t) cal->ADCToKeV( signal ); // signal is converted in KeV
487 // add digit to current cluster
496 if ( ( (nDigits==1) && ( (q==0) || (q>5*oldnoise)) ) || (nDigits>1) ) {
498 Ali1Dcluster& cluster = clusters1D[1][nClusters[1]++];
500 if(q!=0) cluster.SetY(y/q);
501 else cluster.SetY(strip+1);
504 cluster.SetNd(nDigits);
505 cluster.SetLabels(lab);
507 if(repa->GetUseUnfoldingInClusterFinderSSD()==kTRUE) {
509 //Split suspiciously big cluster
510 if (nDigits > 4&&nDigits < 25) {
511 cluster.SetY(y/q - 0.25*nDigits);
513 if (nClusters[1] == kMax) {
514 Error("FindClustersSSD", "Too many 1D clusters !");
517 Ali1Dcluster& cluster2 = clusters1D[1][nClusters[1]++];
518 cluster2.SetY(y/q + 0.25*nDigits);
519 cluster2.SetQ(0.5*q);
520 cluster2.SetNd(nDigits);
521 cluster2.SetLabels(lab);
531 } // loop over strips on N-side
535 if ( ( (nDigits==1) && ( (q==0) || (q>5*oldnoise)) ) || (nDigits>1) ) {
537 Ali1Dcluster& cluster = clusters1D[1][nClusters[1]++];
539 if(q!=0) cluster.SetY(y/q);
540 else cluster.SetY(strip+1);
543 cluster.SetNd(nDigits);
544 cluster.SetLabels(lab);
546 if(repa->GetUseUnfoldingInClusterFinderSSD()==kTRUE) {
548 //Split suspiciously big cluster
549 if (nDigits > 4&&nDigits < 25) {
550 if(q!=0) cluster.SetY(y/q - 0.25*nDigits);
551 else cluster.SetY(strip+1 - 0.25*nDigits);
553 if (nClusters[1] == kMax) {
554 Error("FindClustersSSD", "Too many 1D clusters !");
557 Ali1Dcluster& cluster2 = clusters1D[1][nClusters[1]++];
558 if(q!=0) cluster2.SetY(y/q + 0.25*nDigits);
559 else cluster2.SetY(strip+1 + 0.25*nDigits);
560 cluster2.SetQ(0.5*q);
561 cluster2.SetNd(nDigits);
562 cluster2.SetLabels(lab);
573 if((nClusters[0])&&(nClusters[1])) {
575 clusters[iModule] = new TClonesArray("AliITSRecPoint");
577 FindClustersSSD(&clusters1D[0][0], nClusters[0],
578 &clusters1D[1][0], nClusters[1], clusters[iModule]);
579 Int_t nClustersn = clusters[iModule]->GetEntriesFast();
580 nClustersSSD += nClustersn;
583 nClusters[0] = nClusters[1] = 0;
592 Info("FindClustersSSD", "found clusters in ITS SSD: %d", nClustersSSD);
595 void AliITSClusterFinderV2SSD::
596 FindClustersSSD(Ali1Dcluster* neg, Int_t nn,
597 Ali1Dcluster* pos, Int_t np,
598 TClonesArray *clusters) {
599 //------------------------------------------------------------
600 // Actual SSD cluster finder
601 //------------------------------------------------------------
603 const TGeoHMatrix *mT2L=AliITSgeomTGeo::GetTracking2LocalMatrix(fModule);
605 TClonesArray &cl=*clusters;
607 AliITSsegmentationSSD *seg = dynamic_cast<AliITSsegmentationSSD*>(fDetTypeRec->GetSegmentationModel(2));
608 if (fModule>fLastSSD1)
613 Float_t hwSSD = seg->Dx()*1e-4/2;
614 Float_t hlSSD = seg->Dz()*1e-4/2;
616 Int_t idet=fNdet[fModule];
620 Int_t *cnegative = new Int_t[np];
621 Int_t *cused1 = new Int_t[np];
622 Int_t *negativepair = new Int_t[10*np];
623 Int_t *cpositive = new Int_t[nn];
624 Int_t *cused2 = new Int_t[nn];
625 Int_t *positivepair = new Int_t[10*nn];
626 for (Int_t i=0;i<np;i++) {cnegative[i]=0; cused1[i]=0;}
627 for (Int_t i=0;i<nn;i++) {cpositive[i]=0; cused2[i]=0;}
628 for (Int_t i=0;i<10*np;i++) {negativepair[i]=0;}
629 for (Int_t i=0;i<10*nn;i++) {positivepair[i]=0;}
631 if ((np*nn) > fgPairsSize) {
633 if (fgPairs) delete [] fgPairs;
634 fgPairsSize = 4*np*nn;
635 fgPairs = new Short_t[fgPairsSize];
637 memset(fgPairs,0,sizeof(Short_t)*np*nn);
640 // find available pairs
642 for (Int_t i=0; i<np; i++) {
643 Float_t yp=pos[i].GetY();
644 if ( (pos[i].GetQ()>0) && (pos[i].GetQ()<3) ) continue;
645 for (Int_t j=0; j<nn; j++) {
646 if ( (neg[j].GetQ()>0) && (neg[j].GetQ()<3) ) continue;
647 Float_t yn=neg[j].GetY();
650 seg->GetPadCxz(yn, yp, xt, zt);
652 if (TMath::Abs(xt)<hwSSD+0.01)
653 if (TMath::Abs(zt)<hlSSD+0.01*(neg[j].GetNd()+pos[i].GetNd())) {
654 negativepair[i*10+cnegative[i]] =j; //index
655 positivepair[j*10+cpositive[j]] =i;
656 cnegative[i]++; //counters
664 // try to recover points out of but close to the module boundaries
666 for (Int_t i=0; i<np; i++) {
667 Float_t yp=pos[i].GetY();
668 if ( (pos[i].GetQ()>0) && (pos[i].GetQ()<3) ) continue;
669 for (Int_t j=0; j<nn; j++) {
670 if ( (neg[j].GetQ()>0) && (neg[j].GetQ()<3) ) continue;
671 // if both 1Dclusters have an other cross continue
672 if (cpositive[j]&&cnegative[i]) continue;
673 Float_t yn=neg[j].GetY();
676 seg->GetPadCxz(yn, yp, xt, zt);
678 if (TMath::Abs(xt)<hwSSD+0.1)
679 if (TMath::Abs(zt)<hlSSD+0.15) {
680 // tag 1Dcluster (eventually will produce low quality recpoint)
681 if (cnegative[i]==0) pos[i].SetNd(100); // not available pair
682 if (cpositive[j]==0) neg[j].SetNd(100); // not available pair
683 negativepair[i*10+cnegative[i]] =j; //index
684 positivepair[j*10+cpositive[j]] =i;
685 cnegative[i]++; //counters
698 static AliITSRecoParam *repa = NULL;
700 repa = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
702 repa = AliITSRecoParam::GetHighFluxParam();
703 AliWarning("Using default AliITSRecoParam class");
707 if(repa->GetUseChargeMatchingInClusterFinderSSD()==kTRUE) {
713 for (Int_t ip=0;ip<np;ip++){
714 Float_t xbest=1000,zbest=1000,qbest=0;
716 // select gold clusters
717 if ( (cnegative[ip]==1) && cpositive[negativepair[10*ip]]==1){
718 Float_t yp=pos[ip].GetY();
719 Int_t j = negativepair[10*ip];
721 if( (pos[ip].GetQ()==0) && (neg[j].GetQ() ==0) ) {
722 // both bad, hence continue;
723 // mark both as used (to avoid recover at the end)
729 ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
731 // charge matching (note that if posQ or negQ is 0 -> ratio=1 and the following condition is met
732 if (TMath::Abs(ratio)>0.2) continue; // note: 0.2=3xsigma_ratio calculated in cosmics tests
735 Float_t yn=neg[j].GetY();
738 seg->GetPadCxz(yn, yp, xt, zt);
743 qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
744 if( (pos[ip].GetQ()==0)||(neg[ip].GetQ()==0)) qbest*=2; // in case of bad strips on one side keep all charge from the other one
747 Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
748 mT2L->MasterToLocal(loc,trk);
752 lp[2]=0.0025*0.0025; //SigmaY2
753 lp[3]=0.110*0.110; //SigmaZ2
756 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
757 for (Int_t ilab=0;ilab<3;ilab++){
758 milab[ilab] = pos[ip].GetLabel(ilab);
759 milab[ilab+3] = neg[j].GetLabel(ilab);
763 milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
764 Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
765 AliITSRecPoint * cl2;
767 if(clusters){ // Note clusters != 0 when method is called for rawdata
770 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
772 cl2->SetChargeRatio(ratio);
776 if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
781 if(pos[ip].GetQ()==0) cl2->SetType(3);
782 if(neg[ip].GetQ()==0) cl2->SetType(4);
788 else{ // Note clusters == 0 when method is called for digits
790 cl2 = new AliITSRecPoint(milab,lp,info);
792 cl2->SetChargeRatio(ratio);
796 if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
801 if(pos[ip].GetQ()==0) cl2->SetType(3);
802 if(neg[ip].GetQ()==0) cl2->SetType(4);
807 fDetTypeRec->AddRecPoint(*cl2);
813 for (Int_t ip=0;ip<np;ip++){
814 Float_t xbest=1000,zbest=1000,qbest=0;
817 // select "silber" cluster
818 if ( cnegative[ip]==1 && cpositive[negativepair[10*ip]]==2){
819 Int_t in = negativepair[10*ip];
820 Int_t ip2 = positivepair[10*in];
821 if (ip2==ip) ip2 = positivepair[10*in+1];
822 Float_t pcharge = pos[ip].GetQ()+pos[ip2].GetQ();
824 if ( (TMath::Abs(pcharge-neg[in].GetQ())<30) && (pcharge!=0) ) { //
828 if ( (fgPairs[ip*nn+in]==100)&&(pos[ip].GetQ() ) ) { //
830 Float_t yp=pos[ip].GetY();
831 Float_t yn=neg[in].GetY();
834 seg->GetPadCxz(yn, yp, xt, zt);
838 qbest =pos[ip].GetQ();
839 Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
840 mT2L->MasterToLocal(loc,trk);
844 lp[2]=0.0025*0.0025; //SigmaY2
845 lp[3]=0.110*0.110; //SigmaZ2
848 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
849 for (Int_t ilab=0;ilab<3;ilab++){
850 milab[ilab] = pos[ip].GetLabel(ilab);
851 milab[ilab+3] = neg[in].GetLabel(ilab);
855 ratio = (pos[ip].GetQ()-neg[in].GetQ())/(pos[ip].GetQ()+neg[in].GetQ());
856 milab[3]=(((ip<<10) + in)<<10) + idet; // pos|neg|det
857 Int_t info[3] = {pos[ip].GetNd(),neg[in].GetNd(),fNlayer[fModule]};
859 AliITSRecPoint * cl2;
862 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
863 cl2->SetChargeRatio(ratio);
865 fgPairs[ip*nn+in] = 5;
866 if ((pos[ip].GetNd()+neg[in].GetNd())>6){ //multi cluster
868 fgPairs[ip*nn+in] = 6;
872 cl2 = new AliITSRecPoint(milab,lp,info);
873 cl2->SetChargeRatio(ratio);
875 fgPairs[ip*nn+in] = 5;
876 if ((pos[ip].GetNd()+neg[in].GetNd())>6){ //multi cluster
878 fgPairs[ip*nn+in] = 6;
881 fDetTypeRec->AddRecPoint(*cl2);
890 // if (!(cused1[ip2] || cused2[in])){ //
891 if ( (fgPairs[ip2*nn+in]==100) && (pos[ip2].GetQ()) ) {
893 Float_t yp=pos[ip2].GetY();
894 Float_t yn=neg[in].GetY();
897 seg->GetPadCxz(yn, yp, xt, zt);
901 qbest =pos[ip2].GetQ();
903 Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
904 mT2L->MasterToLocal(loc,trk);
908 lp[2]=0.0025*0.0025; //SigmaY2
909 lp[3]=0.110*0.110; //SigmaZ2
912 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
913 for (Int_t ilab=0;ilab<3;ilab++){
914 milab[ilab] = pos[ip2].GetLabel(ilab);
915 milab[ilab+3] = neg[in].GetLabel(ilab);
919 ratio = (pos[ip2].GetQ()-neg[in].GetQ())/(pos[ip2].GetQ()+neg[in].GetQ());
920 milab[3]=(((ip2<<10) + in)<<10) + idet; // pos|neg|det
921 Int_t info[3] = {pos[ip2].GetNd(),neg[in].GetNd(),fNlayer[fModule]};
923 AliITSRecPoint * cl2;
925 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
927 cl2->SetChargeRatio(ratio);
929 fgPairs[ip2*nn+in] =5;
930 if ((pos[ip2].GetNd()+neg[in].GetNd())>6){ //multi cluster
932 fgPairs[ip2*nn+in] =6;
936 cl2 = new AliITSRecPoint(milab,lp,info);
937 cl2->SetChargeRatio(ratio);
939 fgPairs[ip2*nn+in] =5;
940 if ((pos[ip2].GetNd()+neg[in].GetNd())>6){ //multi cluster
942 fgPairs[ip2*nn+in] =6;
945 fDetTypeRec->AddRecPoint(*cl2);
954 } // charge matching condition
956 } // 2 Pside cross 1 Nside
957 } // loop over Pside clusters
962 for (Int_t jn=0;jn<nn;jn++){
963 if (cused2[jn]) continue;
964 Float_t xbest=1000,zbest=1000,qbest=0;
965 // select "silber" cluster
966 if ( cpositive[jn]==1 && cnegative[positivepair[10*jn]]==2){
967 Int_t ip = positivepair[10*jn];
968 Int_t jn2 = negativepair[10*ip];
969 if (jn2==jn) jn2 = negativepair[10*ip+1];
970 Float_t pcharge = neg[jn].GetQ()+neg[jn2].GetQ();
973 if ( (TMath::Abs(pcharge-pos[ip].GetQ())<30) && // charge matching
974 (pcharge!=0) ) { // reject combinations of bad strips
978 // if (!(cused1[ip]||cused2[jn])){
979 if ( (fgPairs[ip*nn+jn]==100) && (neg[jn].GetQ()) ) { //
981 Float_t yn=neg[jn].GetY();
982 Float_t yp=pos[ip].GetY();
985 seg->GetPadCxz(yn, yp, xt, zt);
989 qbest =neg[jn].GetQ();
992 Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
993 mT2L->MasterToLocal(loc,trk);
997 lp[2]=0.0025*0.0025; //SigmaY2
998 lp[3]=0.110*0.110; //SigmaZ2
1001 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1002 for (Int_t ilab=0;ilab<3;ilab++){
1003 milab[ilab] = pos[ip].GetLabel(ilab);
1004 milab[ilab+3] = neg[jn].GetLabel(ilab);
1007 CheckLabels2(milab);
1008 ratio = (pos[ip].GetQ()-neg[jn].GetQ())/(pos[ip].GetQ()+neg[jn].GetQ());
1009 milab[3]=(((ip<<10) + jn)<<10) + idet; // pos|neg|det
1010 Int_t info[3] = {pos[ip].GetNd(),neg[jn].GetNd(),fNlayer[fModule]};
1012 AliITSRecPoint * cl2;
1014 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
1016 cl2->SetChargeRatio(ratio);
1018 fgPairs[ip*nn+jn] =7;
1019 if ((pos[ip].GetNd()+neg[jn].GetNd())>6){ //multi cluster
1021 fgPairs[ip*nn+jn]=8;
1026 cl2 = new AliITSRecPoint(milab,lp,info);
1027 cl2->SetChargeRatio(ratio);
1029 fgPairs[ip*nn+jn] =7;
1030 if ((pos[ip].GetNd()+neg[jn].GetNd())>6){ //multi cluster
1032 fgPairs[ip*nn+jn]=8;
1035 fDetTypeRec->AddRecPoint(*cl2);
1041 // if (!(cused1[ip]||cused2[jn2])){
1042 if ( (fgPairs[ip*nn+jn2]==100)&&(neg[jn2].GetQ() ) ) { //
1044 Float_t yn=neg[jn2].GetY();
1045 Double_t yp=pos[ip].GetY();
1048 seg->GetPadCxz(yn, yp, xt, zt);
1052 qbest =neg[jn2].GetQ();
1055 Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
1056 mT2L->MasterToLocal(loc,trk);
1060 lp[2]=0.0025*0.0025; //SigmaY2
1061 lp[3]=0.110*0.110; //SigmaZ2
1064 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1065 for (Int_t ilab=0;ilab<3;ilab++){
1066 milab[ilab] = pos[ip].GetLabel(ilab);
1067 milab[ilab+3] = neg[jn2].GetLabel(ilab);
1070 CheckLabels2(milab);
1071 ratio = (pos[ip].GetQ()-neg[jn2].GetQ())/(pos[ip].GetQ()+neg[jn2].GetQ());
1072 milab[3]=(((ip<<10) + jn2)<<10) + idet; // pos|neg|det
1073 Int_t info[3] = {pos[ip].GetNd(),neg[jn2].GetNd(),fNlayer[fModule]};
1074 AliITSRecPoint * cl2;
1076 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
1079 cl2->SetChargeRatio(ratio);
1080 fgPairs[ip*nn+jn2]=7;
1082 if ((pos[ip].GetNd()+neg[jn2].GetNd())>6){ //multi cluster
1084 fgPairs[ip*nn+jn2]=8;
1089 cl2 = new AliITSRecPoint(milab,lp,info);
1090 cl2->SetChargeRatio(ratio);
1091 fgPairs[ip*nn+jn2]=7;
1093 if ((pos[ip].GetNd()+neg[jn2].GetNd())>6){ //multi cluster
1095 fgPairs[ip*nn+jn2]=8;
1098 fDetTypeRec->AddRecPoint(*cl2);
1107 } // charge matching condition
1109 } // 2 Nside cross 1 Pside
1110 } // loop over Pside clusters
1114 for (Int_t ip=0;ip<np;ip++){
1115 Float_t xbest=1000,zbest=1000,qbest=0;
1119 if ( (cnegative[ip]<5) && cpositive[negativepair[10*ip]]<5){
1120 Float_t minchargediff =4.;
1122 for (Int_t di=0;di<cnegative[ip];di++){
1123 Int_t jc = negativepair[ip*10+di];
1124 Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
1125 if (TMath::Abs(chargedif)<minchargediff){
1127 minchargediff = TMath::Abs(chargedif);
1130 if (j<0) continue; // not proper cluster
1133 for (Int_t di=0;di<cnegative[ip];di++){
1134 Int_t jc = negativepair[ip*10+di];
1135 Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
1136 if (TMath::Abs(chargedif)<minchargediff+3.) count++;
1138 if (count>1) continue; // more than one "proper" cluster for positive
1142 for (Int_t dj=0;dj<cpositive[j];dj++){
1143 Int_t ic = positivepair[j*10+dj];
1144 Float_t chargedif = pos[ic].GetQ()-neg[j].GetQ();
1145 if (TMath::Abs(chargedif)<minchargediff+3.) count++;
1147 if (count>1) continue; // more than one "proper" cluster for negative
1152 for (Int_t dj=0;dj<cnegative[jp];dj++){
1153 Int_t ic = positivepair[jp*10+dj];
1154 Float_t chargedif = pos[ic].GetQ()-neg[jp].GetQ();
1155 if (TMath::Abs(chargedif)<minchargediff+4.) count++;
1157 if (count>1) continue;
1158 if (fgPairs[ip*nn+j]<100) continue;
1161 //almost gold clusters
1162 Float_t yp=pos[ip].GetY();
1163 Float_t yn=neg[j].GetY();
1167 seg->GetPadCxz(yn, yp, xt, zt);
1171 qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
1174 Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
1175 mT2L->MasterToLocal(loc,trk);
1179 lp[2]=0.0025*0.0025; //SigmaY2
1180 lp[3]=0.110*0.110; //SigmaZ2
1182 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1183 for (Int_t ilab=0;ilab<3;ilab++){
1184 milab[ilab] = pos[ip].GetLabel(ilab);
1185 milab[ilab+3] = neg[j].GetLabel(ilab);
1188 CheckLabels2(milab);
1189 if ((neg[j].GetQ()==0)&&(pos[ip].GetQ()==0)) continue; // reject crosses of bad strips!!
1190 ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
1191 milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
1192 Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
1193 AliITSRecPoint * cl2;
1195 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
1197 cl2->SetChargeRatio(ratio);
1199 fgPairs[ip*nn+j]=10;
1200 if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
1202 fgPairs[ip*nn+j]=11;
1208 cl2 = new AliITSRecPoint(milab,lp,info);
1209 cl2->SetChargeRatio(ratio);
1211 fgPairs[ip*nn+j]=10;
1212 if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
1214 fgPairs[ip*nn+j]=11;
1219 fDetTypeRec->AddRecPoint(*cl2);
1224 } // loop over Pside 1Dclusters
1227 } // use charge matching
1230 // recover all the other crosses
1232 for (Int_t i=0; i<np; i++) {
1233 Float_t xbest=1000,zbest=1000,qbest=0;
1234 Float_t yp=pos[i].GetY();
1235 if ((pos[i].GetQ()>0)&&(pos[i].GetQ()<3)) continue;
1236 for (Int_t j=0; j<nn; j++) {
1237 // for (Int_t di = 0;di<cpositive[i];di++){
1238 // Int_t j = negativepair[10*i+di];
1239 if ((neg[j].GetQ()>0)&&(neg[j].GetQ()<3)) continue;
1241 if ((neg[j].GetQ()==0)&&(pos[i].GetQ()==0)) continue; // reject crosses of bad strips!!
1243 if (cused2[j]||cused1[i]) continue;
1244 if (fgPairs[i*nn+j]>0 &&fgPairs[i*nn+j]<100) continue;
1245 ratio = (pos[i].GetQ()-neg[j].GetQ())/(pos[i].GetQ()+neg[j].GetQ());
1246 Float_t yn=neg[j].GetY();
1249 seg->GetPadCxz(yn, yp, xt, zt);
1251 if (TMath::Abs(xt)<hwSSD+0.01)
1252 if (TMath::Abs(zt)<hlSSD+0.01*(neg[j].GetNd()+pos[i].GetNd())) {
1255 qbest=0.5*(pos[i].GetQ()+neg[j].GetQ());
1258 Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
1259 mT2L->MasterToLocal(loc,trk);
1263 lp[2]=0.0025*0.0025; //SigmaY2
1264 lp[3]=0.110*0.110; //SigmaZ2
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[i].GetLabel(ilab);
1270 milab[ilab+3] = neg[j].GetLabel(ilab);
1273 CheckLabels2(milab);
1274 milab[3]=(((i<<10) + j)<<10) + idet; // pos|neg|det
1275 Int_t info[3] = {pos[i].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
1276 AliITSRecPoint * cl2;
1278 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
1280 cl2->SetChargeRatio(ratio);
1281 cl2->SetType(100+cpositive[j]+cnegative[i]);
1283 if(pos[i].GetQ()==0) cl2->SetType(200+cpositive[j]+cnegative[i]);
1284 if(neg[j].GetQ()==0) cl2->SetType(300+cpositive[j]+cnegative[i]);
1288 cl2 = new AliITSRecPoint(milab,lp,info);
1289 cl2->SetChargeRatio(ratio);
1290 cl2->SetType(100+cpositive[j]+cnegative[i]);
1292 if(pos[i].GetQ()==0) cl2->SetType(200+cpositive[j]+cnegative[i]);
1293 if(neg[j].GetQ()==0) cl2->SetType(300+cpositive[j]+cnegative[i]);
1295 fDetTypeRec->AddRecPoint(*cl2);
1301 delete [] cnegative;
1303 delete [] negativepair;
1304 delete [] cpositive;
1306 delete [] positivepair;