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"
41 Short_t *AliITSClusterFinderV2SSD::fgPairs = 0x0;
42 Int_t AliITSClusterFinderV2SSD::fgPairsSize = 0;
44 ClassImp(AliITSClusterFinderV2SSD)
47 AliITSClusterFinderV2SSD::AliITSClusterFinderV2SSD(AliITSDetTypeRec* dettyp):AliITSClusterFinderV2(dettyp),
48 fLastSSD1(AliITSgeomTGeo::GetModuleIndex(6,1,1)-1),
60 //______________________________________________________________________
61 AliITSClusterFinderV2SSD::AliITSClusterFinderV2SSD(const AliITSClusterFinderV2SSD &cf) : AliITSClusterFinderV2(cf), fLastSSD1(cf.fLastSSD1),
62 fYpitchSSD(cf.fYpitchSSD),
71 //______________________________________________________________________
72 AliITSClusterFinderV2SSD& AliITSClusterFinderV2SSD::operator=(const AliITSClusterFinderV2SSD& cf ){
73 // Assignment operator
75 this->~AliITSClusterFinderV2SSD();
76 new(this) AliITSClusterFinderV2SSD(cf);
81 void AliITSClusterFinderV2SSD::FindRawClusters(Int_t mod){
85 FindClustersSSD(fDigits);
89 void AliITSClusterFinderV2SSD::FindClustersSSD(TClonesArray *alldigits) {
90 //------------------------------------------------------------
91 // Actual SSD cluster finder
92 //------------------------------------------------------------
94 static AliITSRecoParam *repa = NULL;
98 repa = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
100 repa = AliITSRecoParam::GetHighFluxParam();
101 AliWarning("Using default AliITSRecoParam class");
105 AliITSCalibrationSSD* cal = (AliITSCalibrationSSD*)GetResp(fModule);
108 Int_t smaxall=alldigits->GetEntriesFast();
109 if (smaxall==0) return;
110 // TObjArray *digits = new TObjArray;
112 for (Int_t i=0;i<smaxall; i++){
113 AliITSdigitSSD *d=(AliITSdigitSSD*)alldigits->UncheckedAt(i);
115 if(d->IsSideP()) gain = cal->GetGainP(d->GetStripNumber());
116 else gain = cal->GetGainN(d->GetStripNumber());
118 Float_t q=gain*d->GetSignal(); // calibration brings mip peaks around 120 (in ADC units)
119 q=cal->ADCToKeV(q); // converts the charge in KeV from ADC units
120 //Float_t q=d->GetSignal()/4.29;// temp. fix (for PID purposed - normalis. to be checked)
121 d->SetSignal(Int_t(q));
123 if (d->GetSignal()<3) continue;
126 Int_t smax = digits.GetEntriesFast();
129 const Int_t kMax=1000;
131 Ali1Dcluster pos[kMax], neg[kMax];
132 Float_t y=0., q=0., qmax=0.;
133 Int_t lab[4]={-2,-2,-2,-2};
135 AliITSdigitSSD *d=(AliITSdigitSSD*)digits.UncheckedAt(0);
137 y += d->GetCoord2()*d->GetSignal();
139 lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2);
140 Int_t curr=d->GetCoord2();
141 Int_t flag=d->GetCoord1();
146 for (Int_t ilab=0;ilab<10;ilab++){
149 milab[0]=d->GetTrack(0); milab[1]=d->GetTrack(1); milab[2]=d->GetTrack(2);
151 for (Int_t s=1; s<smax; s++) {
152 d=(AliITSdigitSSD*)digits.UncheckedAt(s);
153 Int_t strip=d->GetCoord2();
154 if ((strip-curr) > 1 || flag!=d->GetCoord1()) {
159 c[*n].SetLabels(milab);
161 if(repa->GetUseUnfoldingInClusterFinderSSD()==kTRUE) {
163 //Split suspiciously big cluster
165 c[*n].SetY(y/q-0.25*nd);
169 Error("FindClustersSSD","Too many 1D clusters !");
172 c[*n].SetY(y/q+0.25*nd);
175 c[*n].SetLabels(milab);
182 Error("FindClustersSSD","Too many 1D clusters !");
187 lab[0]=lab[1]=lab[2]=-2;
189 for (Int_t ilab=0;ilab<10;ilab++){
193 if (flag!=d->GetCoord1()) { n=&np; c=pos; }
197 y += d->GetCoord2()*d->GetSignal();
199 if (d->GetSignal()>qmax) {
201 lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2);
203 for (Int_t ilab=0;ilab<10;ilab++) {
204 if (d->GetTrack(ilab)>=0) AddLabel(milab, (d->GetTrack(ilab)));
211 c[*n].SetLabels(lab);
213 if(repa->GetUseUnfoldingInClusterFinderSSD()==kTRUE) {
215 //Split suspiciously big cluster
217 c[*n].SetY(y/q-0.25*nd);
221 Error("FindClustersSSD","Too many 1D clusters !");
224 c[*n].SetY(y/q+0.25*nd);
227 c[*n].SetLabels(lab);
233 Error("FindClustersSSD","Too many 1D clusters !");
237 FindClustersSSD(neg, nn, pos, np);
241 void AliITSClusterFinderV2SSD::RawdataToClusters(AliRawReader* rawReader,TClonesArray** clusters){
243 //------------------------------------------------------------
244 // This function creates ITS clusters from raw data
245 //------------------------------------------------------------
247 AliITSRawStreamSSD inputSSD(rawReader);
248 FindClustersSSD(&inputSSD,clusters);
252 void AliITSClusterFinderV2SSD::FindClustersSSD(AliITSRawStreamSSD* input,
253 TClonesArray** clusters)
255 //------------------------------------------------------------
256 // Actual SSD cluster finder for raw data
257 //------------------------------------------------------------
259 static AliITSRecoParam *repa = NULL;
261 repa = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
263 repa = AliITSRecoParam::GetHighFluxParam();
264 AliWarning("Using default AliITSRecoParam class");
268 Int_t nClustersSSD = 0;
269 const Int_t kMax = 1000;
270 Ali1Dcluster clusters1D[2][kMax];
271 Int_t nClusters[2] = {0, 0};
272 Int_t lab[3]={-2,-2,-2};
278 // Float_t pedestal=0.;
280 AliITSCalibrationSSD* cal=NULL;
282 Int_t matrix[12][1536];
289 Int_t osignal = 65535;
293 // read raw data input stream
296 // reset signal matrix
297 for(Int_t i=0; i<12; i++) { for(Int_t j=0; j<1536; j++) { matrix[i][j] = 65535;} }
301 matrix[oadc][ostrip] = osignal; // recover data from previous occurence of input->Next()
304 // buffer data for ddl=iddl and ad=iad
307 next = input->Next();
308 if((!next)&&(input->flag)) continue;
309 Int_t ddl=input->GetDDL();
310 Int_t ad=input->GetAD();
311 Int_t adc = input->GetADC(); adc = (adc<6)? adc : adc - 2;
312 Int_t strip = input->GetStrip();
313 if(input->GetSideFlag()) strip=1535-strip;
314 Int_t signal = input->GetSignal();
316 if((ddl==iddl)&&(ad==iad)) {n++; matrix[adc][strip] = signal;}
317 else {oddl=iddl; oad=iad; oadc = adc; ostrip = strip; osignal=signal; iddl=ddl; iad=ad; break;}
319 if(!next) {oddl=iddl; oad=iad; oadc = adc; ostrip = strip; osignal=signal; iddl=ddl; iad=ad; break;}
324 if(!next && oddl<0) break;
326 if(n==0) continue; // first occurence
330 for(Int_t iadc=0; iadc<12; iadc++) { // loop over ADC index for ddl=oddl and ad=oad
332 Int_t iimod = (oad - 1) * 12 + iadc;
333 Int_t iModule = AliITSRawStreamSSD::GetModuleNumber(oddl,iimod);
334 if(iModule==-1) continue;
335 cal = (AliITSCalibrationSSD*)GetResp(iModule);
340 for(Int_t istrip=0; istrip<768; istrip++) { // P-side
341 Int_t signal = matrix[iadc][istrip];
342 pedestal = cal->GetPedestalP(istrip);
343 matrix[iadc][istrip]=signal-(Int_t)pedestal;
349 for(Int_t l=0; l<6; l++) {
351 for(Int_t n=20; n<108; n++) cmode+=matrix[iadc][l*128+n];
353 for(Int_t n=0; n<128; n++) matrix[iadc][l*128+n]-=(Int_t)cmode;
359 for(istrip=0; istrip<768; istrip++) { // P-side
361 Int_t signal = TMath::Abs(matrix[iadc][istrip]);
364 noise = cal->GetNoiseP(istrip); if(noise<1.) signal = 65535;
365 if(signal<3*noise) signal = 65535; // in case ZS was not done in hw do it now
367 // if(cal->IsPChannelB ad(istrip)) signal=0;
370 gain = cal->GetGainP(istrip);
371 signal = (Int_t) ( signal * gain ); // signal is corrected for gain
372 signal = (Int_t) cal->ADCToKeV( signal ); // signal is converted in KeV
374 q += signal; // add digit to current cluster
375 y += istrip * signal;
382 if ( ( (nDigits==1) && ( (q==0) || (q>5*oldnoise)) ) || (nDigits>1) ) {
384 Ali1Dcluster& cluster = clusters1D[0][nClusters[0]++];
386 if(q!=0) cluster.SetY(y/q);
387 else cluster.SetY(istrip-1);
390 cluster.SetNd(nDigits);
391 cluster.SetLabels(lab);
393 if(repa->GetUseUnfoldingInClusterFinderSSD()==kTRUE) {
395 //Split suspiciously big cluster
396 if (nDigits > 4&&nDigits < 25) {
397 if(q!=0) cluster.SetY(y/q - 0.25*nDigits);
398 else cluster.SetY(istrip-1 - 0.25*nDigits);
400 if (nClusters[0] == kMax) {
401 Error("FindClustersSSD", "Too many 1D clusters !");
404 Ali1Dcluster& cluster2 = clusters1D[0][nClusters[0]++];
405 if(q!=0) cluster2.SetY(y/q + 0.25*nDigits);
406 else cluster2.SetY(istrip-1 + 0.25*nDigits);
407 cluster2.SetQ(0.5*q);
408 cluster2.SetNd(nDigits);
409 cluster2.SetLabels(lab);
419 } // loop over strip on P-side
421 // if last strip does have signal
424 if ( ( (nDigits==1) && ( (q==0) || (q>5*oldnoise)) ) || (nDigits>1) ) {
426 Ali1Dcluster& cluster = clusters1D[0][nClusters[0]++];
428 if(q!=0) cluster.SetY(y/q);
429 else cluster.SetY(istrip-1);
432 cluster.SetNd(nDigits);
433 cluster.SetLabels(lab);
435 if(repa->GetUseUnfoldingInClusterFinderSSD()==kTRUE) {
437 //Split suspiciously big cluster
438 if (nDigits > 4&&nDigits < 25) {
439 if(q!=0) cluster.SetY(y/q - 0.25*nDigits);
440 else cluster.SetY(istrip-1 - 0.25*nDigits);
442 if (nClusters[0] == kMax) {
443 Error("FindClustersSSD", "Too many 1D clusters !");
446 Ali1Dcluster& cluster2 = clusters1D[0][nClusters[0]++];
447 if(q!=0) cluster2.SetY(y/q + 0.25*nDigits);
448 else cluster2.SetY(istrip-1 + 0.25*nDigits);
449 cluster2.SetQ(0.5*q);
450 cluster2.SetNd(nDigits);
451 cluster2.SetLabels(lab);
462 for(Int_t istrip=768; istrip<1536; istrip++) { // P-side
463 Int_t signal = matrix[iadc][istrip];
464 pedestal = cal->GetPedestalN(1535-istrip);
465 matrix[iadc][istrip]=signal-(Int_t)pedestal;
470 for(Int_t l=6; l<12; l++) {
472 for(Int_t n=20; n<108; n++) cmode+=matrix[iadc][l*128+n];
474 for(Int_t n=0; n<128; n++) matrix[iadc][l*128+n]-=(Int_t)cmode;
481 for(Int_t iistrip=768; iistrip<1536; iistrip++) { // N-side
483 Int_t signal = TMath::Abs(matrix[iadc][iistrip]);
484 strip = 1535-iistrip;
487 noise = cal->GetNoiseN(strip); if(noise<1.) signal=65535;
489 // if(cal->IsNChannelBad(strip)) signal=0;
491 if(signal<3*noise) signal = 65535; // in case ZS was not done in hw do it now
494 gain = cal->GetGainN(strip);
495 signal = (Int_t) ( signal * gain); // signal is corrected for gain
496 signal = (Int_t) cal->ADCToKeV( signal ); // signal is converted in KeV
498 // add digit to current cluster
507 if ( ( (nDigits==1) && ( (q==0) || (q>5*oldnoise)) ) || (nDigits>1) ) {
509 Ali1Dcluster& cluster = clusters1D[1][nClusters[1]++];
511 if(q!=0) cluster.SetY(y/q);
512 else cluster.SetY(strip+1);
515 cluster.SetNd(nDigits);
516 cluster.SetLabels(lab);
518 if(repa->GetUseUnfoldingInClusterFinderSSD()==kTRUE) {
520 //Split suspiciously big cluster
521 if (nDigits > 4&&nDigits < 25) {
522 cluster.SetY(y/q - 0.25*nDigits);
524 if (nClusters[1] == kMax) {
525 Error("FindClustersSSD", "Too many 1D clusters !");
528 Ali1Dcluster& cluster2 = clusters1D[1][nClusters[1]++];
529 cluster2.SetY(y/q + 0.25*nDigits);
530 cluster2.SetQ(0.5*q);
531 cluster2.SetNd(nDigits);
532 cluster2.SetLabels(lab);
542 } // loop over strips on N-side
546 if ( ( (nDigits==1) && ( (q==0) || (q>5*oldnoise)) ) || (nDigits>1) ) {
548 Ali1Dcluster& cluster = clusters1D[1][nClusters[1]++];
550 if(q!=0) cluster.SetY(y/q);
551 else cluster.SetY(strip+1);
554 cluster.SetNd(nDigits);
555 cluster.SetLabels(lab);
557 if(repa->GetUseUnfoldingInClusterFinderSSD()==kTRUE) {
559 //Split suspiciously big cluster
560 if (nDigits > 4&&nDigits < 25) {
561 if(q!=0) cluster.SetY(y/q - 0.25*nDigits);
562 else cluster.SetY(strip+1 - 0.25*nDigits);
564 if (nClusters[1] == kMax) {
565 Error("FindClustersSSD", "Too many 1D clusters !");
568 Ali1Dcluster& cluster2 = clusters1D[1][nClusters[1]++];
569 if(q!=0) cluster2.SetY(y/q + 0.25*nDigits);
570 else cluster2.SetY(strip+1 + 0.25*nDigits);
571 cluster2.SetQ(0.5*q);
572 cluster2.SetNd(nDigits);
573 cluster2.SetLabels(lab);
584 if((nClusters[0])&&(nClusters[1])) {
586 clusters[iModule] = new TClonesArray("AliITSRecPoint");
588 FindClustersSSD(&clusters1D[0][0], nClusters[0],
589 &clusters1D[1][0], nClusters[1], clusters[iModule]);
590 Int_t nClustersn = clusters[iModule]->GetEntriesFast();
591 nClustersSSD += nClustersn;
594 nClusters[0] = nClusters[1] = 0;
603 Info("FindClustersSSD", "found clusters in ITS SSD: %d", nClustersSSD);
606 void AliITSClusterFinderV2SSD::
607 FindClustersSSD(Ali1Dcluster* neg, Int_t nn,
608 Ali1Dcluster* pos, Int_t np,
609 TClonesArray *clusters) {
610 //------------------------------------------------------------
611 // Actual SSD cluster finder
612 //------------------------------------------------------------
614 const TGeoHMatrix *mT2L=AliITSgeomTGeo::GetTracking2LocalMatrix(fModule);
616 TClonesArray &cl=*clusters;
618 // Float_t tanp=fTanP, tann=fTanN;
619 const Float_t kStartXzero=3.64325;
620 const Float_t kDeltaXzero5or6=0.02239;
621 const Float_t kDeltaZ5to6=7.6/7.0;
625 if (fModule>fLastSSD1) is6=1;
626 Int_t idet=fNdet[fModule];
630 Int_t negativepair[30000];
631 Int_t cnegative[3000];
633 Int_t positivepair[30000];
634 Int_t cpositive[3000];
636 for (Int_t i=0;i<3000;i++) {cnegative[i]=0; cused1[i]=0;}
637 for (Int_t i=0;i<3000;i++) {cpositive[i]=0; cused2[i]=0;}
638 for (Int_t i=0;i<30000;i++) {negativepair[i]=0; positivepair[i]=0;}
640 if ((np*nn) > fgPairsSize) {
642 if (fgPairs) delete [] fgPairs;
643 fgPairsSize = 4*np*nn;
644 fgPairs = new Short_t[fgPairsSize];
646 memset(fgPairs,0,sizeof(Short_t)*np*nn);
649 // find available pairs
651 for (Int_t i=0; i<np; i++) {
652 Float_t yp=pos[i].GetY();
653 if ( (pos[i].GetQ()>0) && (pos[i].GetQ()<3) ) continue;
654 for (Int_t j=0; j<nn; j++) {
655 if ( (neg[j].GetQ()>0) && (neg[j].GetQ()<3) ) continue;
656 Float_t yn=neg[j].GetY();
658 Float_t zt=1.9*(yp-yn)/7.0-is6*kDeltaZ5to6;
659 Float_t yt=is6*(kStartXzero-(285*yp + 1045*yn)/140000.0) + kDeltaXzero5or6;
661 if (TMath::Abs(yt)<fHwSSD+0.01)
662 if (TMath::Abs(zt)<fHlSSD+0.01*(neg[j].GetNd()+pos[i].GetNd())) {
663 negativepair[i*10+cnegative[i]] =j; //index
664 positivepair[j*10+cpositive[j]] =i;
665 cnegative[i]++; //counters
673 // try to recover points out of but close to the module boundaries
675 for (Int_t i=0; i<np; i++) {
676 Float_t yp=pos[i].GetY();
677 if ( (pos[i].GetQ()>0) && (pos[i].GetQ()<3) ) continue;
678 for (Int_t j=0; j<nn; j++) {
679 if ( (neg[j].GetQ()>0) && (neg[j].GetQ()<3) ) continue;
680 // if both 1Dclusters have an other cross continue
681 if (cpositive[j]&&cnegative[i]) continue;
682 Float_t yn=neg[j].GetY();
684 Float_t zt=1.9*(yp-yn)/7.0-is6*kDeltaZ5to6;
685 Float_t yt=is6*(kStartXzero-(285*yp + 1045*yn)/140000.0) + kDeltaXzero5or6;
687 if (TMath::Abs(yt)<fHwSSD+0.1)
688 if (TMath::Abs(zt)<fHlSSD+0.15) {
689 // tag 1Dcluster (eventually will produce low quality recpoint)
690 if (cnegative[i]==0) pos[i].SetNd(100); // not available pair
691 if (cpositive[j]==0) neg[j].SetNd(100); // not available pair
692 negativepair[i*10+cnegative[i]] =j; //index
693 positivepair[j*10+cpositive[j]] =i;
694 cnegative[i]++; //counters
707 static AliITSRecoParam *repa = NULL;
709 repa = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
711 repa = AliITSRecoParam::GetHighFluxParam();
712 AliWarning("Using default AliITSRecoParam class");
716 if(repa->GetUseChargeMatchingInClusterFinderSSD()==kTRUE) {
722 for (Int_t ip=0;ip<np;ip++){
723 Float_t ybest=1000,zbest=1000,qbest=0;
725 // select gold clusters
726 if ( (cnegative[ip]==1) && cpositive[negativepair[10*ip]]==1){
727 Float_t yp=pos[ip].GetY();
728 Int_t j = negativepair[10*ip];
730 if( (pos[ip].GetQ()==0) && (neg[j].GetQ() ==0) ) {
731 // both bad, hence continue;
732 // mark both as used (to avoid recover at the end)
738 ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
740 // charge matching (note that if posQ or negQ is 0 -> ratio=1 and the following condition is met
741 if (TMath::Abs(ratio)>0.2) continue; // note: 0.2=3xsigma_ratio calculated in cosmics tests
744 Float_t yn=neg[j].GetY();
746 Float_t zt=1.9*(yp-yn)/7.0-is6*kDeltaZ5to6;
747 Float_t yt=is6*(kStartXzero-(285*yp + 1045*yn)/140000.0) + kDeltaXzero5or6;
752 qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
753 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
756 Double_t loc[3]={ybest,0.,zbest},trk[3]={0.,0.,0.};
757 mT2L->MasterToLocal(loc,trk);
761 lp[2]=0.0025*0.0025; //SigmaY2
762 lp[3]=0.110*0.110; //SigmaZ2
765 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
766 for (Int_t ilab=0;ilab<3;ilab++){
767 milab[ilab] = pos[ip].GetLabel(ilab);
768 milab[ilab+3] = neg[j].GetLabel(ilab);
772 milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
773 Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
774 AliITSRecPoint * cl2;
776 if(clusters){ // Note clusters != 0 when method is called for rawdata
779 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
781 cl2->SetChargeRatio(ratio);
785 if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
790 if(pos[ip].GetQ()==0) cl2->SetType(3);
791 if(neg[ip].GetQ()==0) cl2->SetType(4);
797 else{ // Note clusters == 0 when method is called for digits
799 cl2 = new AliITSRecPoint(milab,lp,info);
801 cl2->SetChargeRatio(ratio);
805 if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
810 if(pos[ip].GetQ()==0) cl2->SetType(3);
811 if(neg[ip].GetQ()==0) cl2->SetType(4);
816 fDetTypeRec->AddRecPoint(*cl2);
822 for (Int_t ip=0;ip<np;ip++){
823 Float_t ybest=1000,zbest=1000,qbest=0;
826 // select "silber" cluster
827 if ( cnegative[ip]==1 && cpositive[negativepair[10*ip]]==2){
828 Int_t in = negativepair[10*ip];
829 Int_t ip2 = positivepair[10*in];
830 if (ip2==ip) ip2 = positivepair[10*in+1];
831 Float_t pcharge = pos[ip].GetQ()+pos[ip2].GetQ();
833 if ( (TMath::Abs(pcharge-neg[in].GetQ())<30) && (pcharge!=0) ) { //
837 if ( (fgPairs[ip*nn+in]==100)&&(pos[ip].GetQ() ) ) { //
839 Float_t yp=pos[ip].GetY();
840 Float_t yn=neg[in].GetY();
842 Float_t zt=1.9*(yp-yn)/7.0-is6*kDeltaZ5to6;
843 Float_t yt=is6*(kStartXzero-(285*yp + 1045*yn)/140000.0) + kDeltaXzero5or6;
847 qbest =pos[ip].GetQ();
848 Double_t loc[3]={ybest,0.,zbest},trk[3]={0.,0.,0.};
849 mT2L->MasterToLocal(loc,trk);
853 lp[2]=0.0025*0.0025; //SigmaY2
854 lp[3]=0.110*0.110; //SigmaZ2
857 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
858 for (Int_t ilab=0;ilab<3;ilab++){
859 milab[ilab] = pos[ip].GetLabel(ilab);
860 milab[ilab+3] = neg[in].GetLabel(ilab);
864 ratio = (pos[ip].GetQ()-neg[in].GetQ())/(pos[ip].GetQ()+neg[in].GetQ());
865 milab[3]=(((ip<<10) + in)<<10) + idet; // pos|neg|det
866 Int_t info[3] = {pos[ip].GetNd(),neg[in].GetNd(),fNlayer[fModule]};
868 AliITSRecPoint * cl2;
871 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
872 cl2->SetChargeRatio(ratio);
874 fgPairs[ip*nn+in] = 5;
875 if ((pos[ip].GetNd()+neg[in].GetNd())>6){ //multi cluster
877 fgPairs[ip*nn+in] = 6;
881 cl2 = new AliITSRecPoint(milab,lp,info);
882 cl2->SetChargeRatio(ratio);
884 fgPairs[ip*nn+in] = 5;
885 if ((pos[ip].GetNd()+neg[in].GetNd())>6){ //multi cluster
887 fgPairs[ip*nn+in] = 6;
890 fDetTypeRec->AddRecPoint(*cl2);
899 // if (!(cused1[ip2] || cused2[in])){ //
900 if ( (fgPairs[ip2*nn+in]==100) && (pos[ip2].GetQ()) ) {
902 Float_t yp=pos[ip2].GetY();
903 Float_t yn=neg[in].GetY();
905 Float_t zt=1.9*(yp-yn)/7.0-is6*kDeltaZ5to6;
906 Float_t yt=is6*(kStartXzero-(285*yp + 1045*yn)/140000.0) + kDeltaXzero5or6;
909 qbest =pos[ip2].GetQ();
911 Double_t loc[3]={ybest,0.,zbest},trk[3]={0.,0.,0.};
912 mT2L->MasterToLocal(loc,trk);
916 lp[2]=0.0025*0.0025; //SigmaY2
917 lp[3]=0.110*0.110; //SigmaZ2
920 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
921 for (Int_t ilab=0;ilab<3;ilab++){
922 milab[ilab] = pos[ip2].GetLabel(ilab);
923 milab[ilab+3] = neg[in].GetLabel(ilab);
927 ratio = (pos[ip2].GetQ()-neg[in].GetQ())/(pos[ip2].GetQ()+neg[in].GetQ());
928 milab[3]=(((ip2<<10) + in)<<10) + idet; // pos|neg|det
929 Int_t info[3] = {pos[ip2].GetNd(),neg[in].GetNd(),fNlayer[fModule]};
931 AliITSRecPoint * cl2;
933 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
935 cl2->SetChargeRatio(ratio);
937 fgPairs[ip2*nn+in] =5;
938 if ((pos[ip2].GetNd()+neg[in].GetNd())>6){ //multi cluster
940 fgPairs[ip2*nn+in] =6;
944 cl2 = new AliITSRecPoint(milab,lp,info);
945 cl2->SetChargeRatio(ratio);
947 fgPairs[ip2*nn+in] =5;
948 if ((pos[ip2].GetNd()+neg[in].GetNd())>6){ //multi cluster
950 fgPairs[ip2*nn+in] =6;
953 fDetTypeRec->AddRecPoint(*cl2);
962 } // charge matching condition
964 } // 2 Pside cross 1 Nside
965 } // loop over Pside clusters
970 for (Int_t jn=0;jn<nn;jn++){
971 if (cused2[jn]) continue;
972 Float_t ybest=1000,zbest=1000,qbest=0;
973 // select "silber" cluster
974 if ( cpositive[jn]==1 && cnegative[positivepair[10*jn]]==2){
975 Int_t ip = positivepair[10*jn];
976 Int_t jn2 = negativepair[10*ip];
977 if (jn2==jn) jn2 = negativepair[10*ip+1];
978 Float_t pcharge = neg[jn].GetQ()+neg[jn2].GetQ();
981 if ( (TMath::Abs(pcharge-pos[ip].GetQ())<30) && // charge matching
982 (pcharge!=0) ) { // reject combinations of bad strips
986 // if (!(cused1[ip]||cused2[jn])){
987 if ( (fgPairs[ip*nn+jn]==100) && (neg[jn].GetQ()) ) { //
989 Float_t yn=neg[jn].GetY();
990 Float_t yp=pos[ip].GetY();
992 Float_t zt=1.9*(yp-yn)/7.0-is6*kDeltaZ5to6;
993 Float_t yt=is6*(kStartXzero-(285*yp + 1045*yn)/140000.0) + kDeltaXzero5or6;
998 qbest =neg[jn].GetQ();
1001 Double_t loc[3]={ybest,0.,zbest},trk[3]={0.,0.,0.};
1002 mT2L->MasterToLocal(loc,trk);
1006 lp[2]=0.0025*0.0025; //SigmaY2
1007 lp[3]=0.110*0.110; //SigmaZ2
1010 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1011 for (Int_t ilab=0;ilab<3;ilab++){
1012 milab[ilab] = pos[ip].GetLabel(ilab);
1013 milab[ilab+3] = neg[jn].GetLabel(ilab);
1016 CheckLabels2(milab);
1017 ratio = (pos[ip].GetQ()-neg[jn].GetQ())/(pos[ip].GetQ()+neg[jn].GetQ());
1018 milab[3]=(((ip<<10) + jn)<<10) + idet; // pos|neg|det
1019 Int_t info[3] = {pos[ip].GetNd(),neg[jn].GetNd(),fNlayer[fModule]};
1021 AliITSRecPoint * cl2;
1023 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
1025 cl2->SetChargeRatio(ratio);
1027 fgPairs[ip*nn+jn] =7;
1028 if ((pos[ip].GetNd()+neg[jn].GetNd())>6){ //multi cluster
1030 fgPairs[ip*nn+jn]=8;
1035 cl2 = new AliITSRecPoint(milab,lp,info);
1036 cl2->SetChargeRatio(ratio);
1038 fgPairs[ip*nn+jn] =7;
1039 if ((pos[ip].GetNd()+neg[jn].GetNd())>6){ //multi cluster
1041 fgPairs[ip*nn+jn]=8;
1044 fDetTypeRec->AddRecPoint(*cl2);
1050 // if (!(cused1[ip]||cused2[jn2])){
1051 if ( (fgPairs[ip*nn+jn2]==100)&&(neg[jn2].GetQ() ) ) { //
1053 Float_t yn=neg[jn2].GetY();
1054 Double_t yp=pos[ip].GetY();
1055 Double_t zt=1.9*(yp-yn)/7.0-is6*kDeltaZ5to6;
1056 Double_t yt=is6*(kStartXzero-(285*yp + 1045*yn)/140000.0) + kDeltaXzero5or6;
1058 ybest =yt; zbest=zt;
1061 qbest =neg[jn2].GetQ();
1064 Double_t loc[3]={ybest,0.,zbest},trk[3]={0.,0.,0.};
1065 mT2L->MasterToLocal(loc,trk);
1069 lp[2]=0.0025*0.0025; //SigmaY2
1070 lp[3]=0.110*0.110; //SigmaZ2
1073 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1074 for (Int_t ilab=0;ilab<3;ilab++){
1075 milab[ilab] = pos[ip].GetLabel(ilab);
1076 milab[ilab+3] = neg[jn2].GetLabel(ilab);
1079 CheckLabels2(milab);
1080 ratio = (pos[ip].GetQ()-neg[jn2].GetQ())/(pos[ip].GetQ()+neg[jn2].GetQ());
1081 milab[3]=(((ip<<10) + jn2)<<10) + idet; // pos|neg|det
1082 Int_t info[3] = {pos[ip].GetNd(),neg[jn2].GetNd(),fNlayer[fModule]};
1083 AliITSRecPoint * cl2;
1085 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
1088 cl2->SetChargeRatio(ratio);
1089 fgPairs[ip*nn+jn2]=7;
1091 if ((pos[ip].GetNd()+neg[jn2].GetNd())>6){ //multi cluster
1093 fgPairs[ip*nn+jn2]=8;
1098 cl2 = new AliITSRecPoint(milab,lp,info);
1099 cl2->SetChargeRatio(ratio);
1100 fgPairs[ip*nn+jn2]=7;
1102 if ((pos[ip].GetNd()+neg[jn2].GetNd())>6){ //multi cluster
1104 fgPairs[ip*nn+jn2]=8;
1107 fDetTypeRec->AddRecPoint(*cl2);
1116 } // charge matching condition
1118 } // 2 Nside cross 1 Pside
1119 } // loop over Pside clusters
1123 for (Int_t ip=0;ip<np;ip++){
1124 Float_t ybest=1000,zbest=1000,qbest=0;
1128 if ( (cnegative[ip]<5) && cpositive[negativepair[10*ip]]<5){
1129 Float_t minchargediff =4.;
1131 for (Int_t di=0;di<cnegative[ip];di++){
1132 Int_t jc = negativepair[ip*10+di];
1133 Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
1134 if (TMath::Abs(chargedif)<minchargediff){
1136 minchargediff = TMath::Abs(chargedif);
1139 if (j<0) continue; // not proper cluster
1142 for (Int_t di=0;di<cnegative[ip];di++){
1143 Int_t jc = negativepair[ip*10+di];
1144 Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
1145 if (TMath::Abs(chargedif)<minchargediff+3.) count++;
1147 if (count>1) continue; // more than one "proper" cluster for positive
1151 for (Int_t dj=0;dj<cpositive[j];dj++){
1152 Int_t ic = positivepair[j*10+dj];
1153 Float_t chargedif = pos[ic].GetQ()-neg[j].GetQ();
1154 if (TMath::Abs(chargedif)<minchargediff+3.) count++;
1156 if (count>1) continue; // more than one "proper" cluster for negative
1161 for (Int_t dj=0;dj<cnegative[jp];dj++){
1162 Int_t ic = positivepair[jp*10+dj];
1163 Float_t chargedif = pos[ic].GetQ()-neg[jp].GetQ();
1164 if (TMath::Abs(chargedif)<minchargediff+4.) count++;
1166 if (count>1) continue;
1167 if (fgPairs[ip*nn+j]<100) continue;
1170 //almost gold clusters
1171 Float_t yp=pos[ip].GetY();
1172 Float_t yn=neg[j].GetY();
1174 Float_t zt=1.9*(yp-yn)/7.0-is6*kDeltaZ5to6;
1175 Float_t yt=is6*(kStartXzero-(285*yp + 1045*yn)/140000.0) + kDeltaXzero5or6;
1178 qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
1181 Double_t loc[3]={ybest,0.,zbest},trk[3]={0.,0.,0.};
1182 mT2L->MasterToLocal(loc,trk);
1186 lp[2]=0.0025*0.0025; //SigmaY2
1187 lp[3]=0.110*0.110; //SigmaZ2
1189 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1190 for (Int_t ilab=0;ilab<3;ilab++){
1191 milab[ilab] = pos[ip].GetLabel(ilab);
1192 milab[ilab+3] = neg[j].GetLabel(ilab);
1195 CheckLabels2(milab);
1196 if ((neg[j].GetQ()==0)&&(pos[ip].GetQ()==0)) continue; // reject crosses of bad strips!!
1197 ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
1198 milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
1199 Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
1200 AliITSRecPoint * cl2;
1202 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
1204 cl2->SetChargeRatio(ratio);
1206 fgPairs[ip*nn+j]=10;
1207 if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
1209 fgPairs[ip*nn+j]=11;
1215 cl2 = new AliITSRecPoint(milab,lp,info);
1216 cl2->SetChargeRatio(ratio);
1218 fgPairs[ip*nn+j]=10;
1219 if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
1221 fgPairs[ip*nn+j]=11;
1226 fDetTypeRec->AddRecPoint(*cl2);
1231 } // loop over Pside 1Dclusters
1234 } // use charge matching
1237 // recover all the other crosses
1239 for (Int_t i=0; i<np; i++) {
1240 Float_t ybest=1000,zbest=1000,qbest=0;
1241 Float_t yp=pos[i].GetY();
1242 if ((pos[i].GetQ()>0)&&(pos[i].GetQ()<3)) continue;
1243 for (Int_t j=0; j<nn; j++) {
1244 // for (Int_t di = 0;di<cpositive[i];di++){
1245 // Int_t j = negativepair[10*i+di];
1246 if ((neg[j].GetQ()>0)&&(neg[j].GetQ()<3)) continue;
1248 if ((neg[j].GetQ()==0)&&(pos[i].GetQ()==0)) continue; // reject crosses of bad strips!!
1250 if (cused2[j]||cused1[i]) continue;
1251 if (fgPairs[i*nn+j]>0 &&fgPairs[i*nn+j]<100) continue;
1252 ratio = (pos[i].GetQ()-neg[j].GetQ())/(pos[i].GetQ()+neg[j].GetQ());
1253 Float_t yn=neg[j].GetY();
1255 Float_t zt=1.9*(yp-yn)/7.0-is6*kDeltaZ5to6;
1256 Float_t yt=is6*(kStartXzero-(285*yp + 1045*yn)/140000.0) + kDeltaXzero5or6;
1258 if (TMath::Abs(yt)<fHwSSD+0.01)
1259 if (TMath::Abs(zt)<fHlSSD+0.01*(neg[j].GetNd()+pos[i].GetNd())) {
1261 qbest=0.5*(pos[i].GetQ()+neg[j].GetQ());
1264 Double_t loc[3]={ybest,0.,zbest},trk[3]={0.,0.,0.};
1265 mT2L->MasterToLocal(loc,trk);
1269 lp[2]=0.0025*0.0025; //SigmaY2
1270 lp[3]=0.110*0.110; //SigmaZ2
1273 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1274 for (Int_t ilab=0;ilab<3;ilab++){
1275 milab[ilab] = pos[i].GetLabel(ilab);
1276 milab[ilab+3] = neg[j].GetLabel(ilab);
1279 CheckLabels2(milab);
1280 milab[3]=(((i<<10) + j)<<10) + idet; // pos|neg|det
1281 Int_t info[3] = {pos[i].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
1282 AliITSRecPoint * cl2;
1284 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
1286 cl2->SetChargeRatio(ratio);
1287 cl2->SetType(100+cpositive[j]+cnegative[i]);
1289 if(pos[i].GetQ()==0) cl2->SetType(200+cpositive[j]+cnegative[i]);
1290 if(neg[j].GetQ()==0) cl2->SetType(300+cpositive[j]+cnegative[i]);
1294 cl2 = new AliITSRecPoint(milab,lp,info);
1295 cl2->SetChargeRatio(ratio);
1296 cl2->SetType(100+cpositive[j]+cnegative[i]);
1298 if(pos[i].GetQ()==0) cl2->SetType(200+cpositive[j]+cnegative[i]);
1299 if(neg[j].GetQ()==0) cl2->SetType(300+cpositive[j]+cnegative[i]);
1301 fDetTypeRec->AddRecPoint(*cl2);