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 "AliITSClusterFinderV2SSD.h"
29 #include <Riostream.h>
30 #include <TGeoGlobalMagField.h>
34 #include "AliITSRecPoint.h"
35 #include "AliITSRecPointContainer.h"
36 #include "AliITSgeomTGeo.h"
37 #include "AliITSDetTypeRec.h"
38 #include "AliRawReader.h"
39 #include "AliITSRawStreamSSD.h"
40 #include <TClonesArray.h>
41 #include <TCollection.h>
42 #include "AliITSdigitSSD.h"
43 #include "AliITSReconstructor.h"
44 #include "AliITSCalibrationSSD.h"
45 #include "AliITSsegmentationSSD.h"
47 Short_t *AliITSClusterFinderV2SSD::fgPairs = 0x0;
48 Int_t AliITSClusterFinderV2SSD::fgPairsSize = 0;
49 const Float_t AliITSClusterFinderV2SSD::fgkThreshold = 5.;
51 const Float_t AliITSClusterFinderV2SSD::fgkCosmic2008StripShifts[16][9] =
52 {{-0.35,-0.35,-0.35,-0.35,-0.35,-0.35,-0.35,-0.35,-0.35}, // DDL 512
53 {-0.35,-0.35,-0.35,-0.35,-0.35,-0.35,-0.35,-0.35,-0.35}, // DDL 513
54 {-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15}, // DDL 514
55 {-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15}, // DDL 515
56 { 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00}, // DDL 516
57 { 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00}, // DDL 517
58 {-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15}, // DDL 518
59 {-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15}, // DDL 519
60 {-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.25,-0.15}, // DDL 520
61 {-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15}, // DDL 521
62 {-0.10,-0.10,-0.10,-0.40,-0.40,-0.40,-0.10,-0.10,-0.45}, // DDL 522
63 {-0.10,-0.10,-0.10,-0.35,-0.35,-0.35,-0.10,-0.35,-0.50}, // DDL 523
64 { 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00}, // DDL 524
65 { 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00}, // DDL 525
66 { 0.35, 0.35, 0.35, 0.35, 0.35, 0.35, 0.35, 0.35, 0.35}, // DDL 526
67 { 0.45, 0.45, 0.45, 0.45, 0.45, 0.45, 0.45, 0.45, 0.45}}; // DDL 527
69 ClassImp(AliITSClusterFinderV2SSD)
72 AliITSClusterFinderV2SSD::AliITSClusterFinderV2SSD(AliITSDetTypeRec* dettyp):AliITSClusterFinder(dettyp),fLastSSD1(AliITSgeomTGeo::GetModuleIndex(6,1,1)-1), fLorentzShiftP(0), fLorentzShiftN(0)
75 static AliITSRecoParam *repa = NULL;
77 repa = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
79 repa = AliITSRecoParam::GetHighFluxParam();
80 AliWarning("Using default AliITSRecoParam class");
84 if (repa->GetCorrectLorentzAngleSSD()) {
85 AliMagF* field = dynamic_cast<AliMagF*>(TGeoGlobalMagField::Instance()->GetField());
87 AliError("Cannot get magnetic field from TGeoGlobalMagField");
90 Float_t Bfield = field->SolenoidField();
91 // NB: spatial shift has opposite sign for lay 5 and 6, but strip numbering also changes direction, so no sign-change
92 // Shift due to ExB on drift N-side, units: strip width
93 fLorentzShiftP = -repa->GetTanLorentzAngleElectronsSSD() * 150.e-4/95.e-4 * Bfield / 5.0;
94 // Shift due to ExB on drift P-side, units: strip width
95 fLorentzShiftN = -repa->GetTanLorentzAngleHolesSSD() * 150.e-4/95.e-4 * Bfield / 5.0;
96 AliDebug(1,Form("Bfield %f Lorentz Shift P-side %f N-side %f",Bfield,fLorentzShiftN,fLorentzShiftP));
101 //______________________________________________________________________
102 AliITSClusterFinderV2SSD::AliITSClusterFinderV2SSD(const AliITSClusterFinderV2SSD &cf) : AliITSClusterFinder(cf), fLastSSD1(cf.fLastSSD1), fLorentzShiftP(cf.fLorentzShiftP), fLorentzShiftN(cf.fLorentzShiftN)
107 //______________________________________________________________________
108 AliITSClusterFinderV2SSD& AliITSClusterFinderV2SSD::operator=(const AliITSClusterFinderV2SSD& cf ){
109 // Assignment operator
111 this->~AliITSClusterFinderV2SSD();
112 new(this) AliITSClusterFinderV2SSD(cf);
117 void AliITSClusterFinderV2SSD::FindRawClusters(Int_t mod){
121 FindClustersSSD(fDigits);
125 void AliITSClusterFinderV2SSD::FindClustersSSD(TClonesArray *alldigits) {
126 //------------------------------------------------------------
127 // Actual SSD cluster finder
128 //------------------------------------------------------------
129 Int_t smaxall=alldigits->GetEntriesFast();
130 if (smaxall==0) return;
133 //---------------------------------------
134 // load recoparam and calibration
136 static AliITSRecoParam *repa = NULL;
138 repa = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
140 repa = AliITSRecoParam::GetHighFluxParam();
141 AliWarning("Using default AliITSRecoParam class");
145 AliITSCalibrationSSD* cal = (AliITSCalibrationSSD*)GetResp(fModule);
148 //---------------------------------------
151 //------------------------------------
152 // fill the digits array with zero-suppression condition
153 // Signal is converted in KeV
156 for (Int_t i=0;i<smaxall; i++){
157 AliITSdigitSSD *d=(AliITSdigitSSD*)alldigits->UncheckedAt(i);
159 if(d->IsSideP()) noise = cal->GetNoiseP(d->GetStripNumber());
160 else noise = cal->GetNoiseN(d->GetStripNumber());
161 if (d->GetSignal()<3.*noise) continue;
163 if(d->IsSideP()) gain = cal->GetGainP(d->GetStripNumber());
164 else gain = cal->GetGainN(d->GetStripNumber());
166 Float_t q=gain*d->GetSignal(); //
167 q=cal->ADCToKeV(q); // converts the charge in KeV from ADC units
168 d->SetSignal(Int_t(q));
172 Int_t smax = digits.GetEntriesFast();
174 //------------------------------------
177 const Int_t kMax=1000;
179 Ali1Dcluster pos[kMax], neg[kMax];
180 Float_t y=0., q=0., qmax=0.;
181 Int_t lab[4]={-2,-2,-2,-2};
185 cout<<"-----------------------------"<<endl;
186 cout<<"this is module "<<fModule;
191 if (fModule>fLastSSD1)
194 //--------------------------------------------------------
195 // start 1D-clustering from the first digit in the digits array
197 AliITSdigitSSD *d=(AliITSdigitSSD*)digits.UncheckedAt(0);
199 y += d->GetCoord2()*d->GetSignal();
201 lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2);
204 noise = cal->GetNoiseP(d->GetStripNumber());
205 gain = cal->GetGainP(d->GetStripNumber());
208 noise = cal->GetNoiseN(d->GetStripNumber());
209 gain = cal->GetGainN(d->GetStripNumber());
212 noise=cal->ADCToKeV(noise); // converts noise in KeV from ADC units
214 if(qmax>fgkThreshold*noise) flag5=1; // seed for the cluster
217 cout<<d->GetSignal()<<" "<<noise<<" "<<flag5<<" "<<
218 d->GetCoord1()<<" "<<d->GetCoord2()<<endl;
221 Int_t curr=d->GetCoord2();
222 Int_t flag=d->GetCoord1();
224 // Note: the first side which will be processed is supposed to be the
225 // P-side which is neg
228 if(flag) {n=&np; c=pos;} // in case we have only Nstrips (P was bad!)
232 for (Int_t ilab=0;ilab<10;ilab++){
235 milab[0]=d->GetTrack(0); milab[1]=d->GetTrack(1); milab[2]=d->GetTrack(2);
238 //----------------------------------------------------------
239 // search for neighboring digits
241 for (Int_t s=1; s<smax; s++) {
242 d=(AliITSdigitSSD*)digits.UncheckedAt(s);
243 Int_t strip=d->GetCoord2();
245 // if digits is not a neighbour or side did not change
246 // and at least one of the previous digits met the seed condition
247 // then creates a new 1D cluster
248 if ( ( ((strip-curr) > 1) || (flag!=d->GetCoord1()) ) ) {
251 //cout<<"here1"<<endl;
252 Float_t dLorentz = 0;
253 if (!flag) { // P-side is neg clust
254 dLorentz = fLorentzShiftN;
256 else { // N-side is p clust
257 dLorentz = fLorentzShiftP;
259 c[*n].SetY(y/q+dLorentz);
263 c[*n].SetLabels(milab);
265 if(repa->GetUseUnfoldingInClusterFinderSSD()==kTRUE) {
266 // Note: fUseUnfoldingInClusterFinderSSD=kFALSE by default in RecoParam
268 //Split suspiciously big cluster
270 c[*n].SetY(y/q-0.25*nd+dLorentz);
274 Error("FindClustersSSD","Too many 1D clusters !");
277 c[*n].SetY(y/q+0.25*nd+dLorentz);
280 c[*n].SetLabels(milab);
287 Error("FindClustersSSD","Too many 1D clusters !");
297 lab[0]=lab[1]=lab[2]=-2;
298 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
300 // if side changed from P to N, switch to pos 1D clusters
301 // (if for some reason the side changed from N to P then do the opposite)
302 if (flag!=d->GetCoord1())
303 { if(!flag) {n=&np; c=pos;} else {n=&nn; c=neg;} }
305 } // end create new 1D cluster from previous neighboring digits
307 // continues adding digits to the previous cluster
308 // or start a new one
311 y += d->GetCoord2()*d->GetSignal();
315 noise = cal->GetNoiseP(d->GetStripNumber());
316 gain = cal->GetGainP(d->GetStripNumber());
319 noise = cal->GetNoiseN(d->GetStripNumber());
320 gain = cal->GetGainN(d->GetStripNumber());
323 noise=cal->ADCToKeV(noise); // converts the charge in KeV from ADC units
325 if(d->GetSignal()>fgkThreshold*noise) flag5=1;
328 cout<<d->GetSignal()<<" "<<noise<<" "<<flag5<<" "<<
329 d->GetCoord1()<<" "<<d->GetCoord2()<<endl;
332 if (d->GetSignal()>qmax) {
334 lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2);
336 for (Int_t ilab=0;ilab<10;ilab++) {
337 if (d->GetTrack(ilab)>=0) AddLabel(milab, (d->GetTrack(ilab)));
342 } // loop over digits, no more digits in the digits array
345 // add the last 1D cluster
348 // cout<<"here2"<<endl;
349 Float_t dLorentz = 0;
350 if (!flag) { // P-side is neg clust
351 dLorentz = fLorentzShiftN;
353 else { // N-side is p clust
354 dLorentz = fLorentzShiftP;
357 c[*n].SetY(y/q + dLorentz);
360 c[*n].SetLabels(lab);
362 if(repa->GetUseUnfoldingInClusterFinderSSD()==kTRUE) {
364 //Split suspiciously big cluster
366 c[*n].SetY(y/q-0.25*nd + dLorentz);
370 Error("FindClustersSSD","Too many 1D clusters !");
373 c[*n].SetY(y/q+0.25*nd + dLorentz);
376 c[*n].SetLabels(lab);
382 Error("FindClustersSSD","Too many 1D clusters !");
386 } // if flag5 last 1D cluster added
389 //------------------------------------------------------
390 // call FindClustersSSD to pair neg and pos 1D clusters
391 // and create recpoints from the crosses
392 // Note1: neg are Pside and pos are Nside!!
393 // Note2: if there are no Pside digits nn=0 (bad strips!!) (same for Nside)
395 // cout<<nn<<" Pside and "<<np<<" Nside clusters"<<endl;
397 AliITSRecPointContainer* rpc = AliITSRecPointContainer::Instance();
399 TClonesArray* clusters = rpc->UncheckedGetClusters(fModule);
401 FindClustersSSD(neg, nn, pos, np, clusters);
404 while ((irp = (AliITSRecPoint*)itr.Next())) fDetTypeRec->AddRecPoint(*irp);
406 //-----------------------------------------------------
410 void AliITSClusterFinderV2SSD::RawdataToClusters(AliRawReader* rawReader){
412 //------------------------------------------------------------
413 // This function creates ITS clusters from raw data
414 //------------------------------------------------------------
416 AliITSRawStreamSSD inputSSD(rawReader);
417 FindClustersSSD(&inputSSD);
422 void AliITSClusterFinderV2SSD::FindClustersSSD(AliITSRawStreamSSD* input)
424 //------------------------------------------------------------
425 // Actual SSD cluster finder for raw data
426 //------------------------------------------------------------
428 AliITSRecPointContainer* rpc = AliITSRecPointContainer::Instance();
429 static AliITSRecoParam *repa = NULL;
431 repa = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
433 repa = AliITSRecoParam::GetHighFluxParam();
434 AliWarning("Using default AliITSRecoParam class");
437 Int_t nClustersSSD = 0;
438 const Int_t kNADC = 12;
439 const Int_t kMaxADCClusters = 1000;
441 Int_t strips[kNADC][2][kMaxADCClusters][2]; // [ADC],[side],[istrip], [0]=istrip [1]=signal
442 Int_t nStrips[kNADC][2];
444 for( int i=0; i<kNADC; i++ ){
453 //* Loop over modules DDL+AD
458 bool next = input->Next();
461 //* Continue if corrupted input
464 if( (!next)&&(input->flag) ){
465 AliWarning(Form("ITSClustersFinderSSD: Corrupted data: warning from RawReader"));
469 Int_t newDDL = input->GetDDL();
470 Int_t newAD = input->GetAD();
473 if( newDDL<0 || newDDL>15 ){
474 AliWarning(Form("ITSClustersFinderSSD: Corrupted data: wrong DDL number (%d)",newDDL));
478 if( newAD<1 || newAD>9 ){
479 AliWarning(Form("ITSClustersFinderSSD: Corrupted data: wrong AD number (%d)",newAD));
484 bool newModule = ( !next || ddl!= newDDL || ad!=newAD );
486 if( newModule && ddl>=0 && ad>=0 ){
489 //* Reconstruct the previous block of 12 modules --- actual clusterfinder
492 for( int adc = 0; adc<kNADC; adc++ ){
496 Ali1Dcluster clusters1D[2][kMaxADCClusters]; // per ADC, per side
497 Int_t nClusters1D[2] = {0,0};
498 //int nstat[2] = {0,0};
499 fModule = AliITSRawStreamSSD::GetModuleNumber(ddl, (ad - 1) * 12 + adc );
502 // AliWarning(Form("ITSClustersFinderSSD: Corrupted data: module (ddl %d ad %d adc %d) not found in the map",ddl,ad,adc));
503 //CM channels are always present even everything is suppressed
508 if (fModule>fLastSSD1)
511 AliITSCalibrationSSD* cal = (AliITSCalibrationSSD*)fDetTypeRec->GetCalibrationModel(fModule);
513 AliWarning(Form("ITSClustersFinderSSD: No calibration found for module (ddl %d ad %d adc %d)",ddl,ad,adc));
519 if( repa->GetUseCosmicRunShiftsSSD()) { // Special condition for 2007/2008 cosmic data
520 dStrip = fgkCosmic2008StripShifts[ddl][ad-1];
521 if (TMath::Abs(dStrip) > 1.5){
522 AliWarning(Form("Indexing error in Cosmic calibration: ddl = %d, dStrip %f\n",ddl,dStrip));
527 for( int side=0; side<=1; side++ ){
529 Int_t lab[3]={-2,-2,-2};
536 Float_t dLorentz = 0;
537 if (side==0) { // P-side is neg clust
538 dLorentz = fLorentzShiftN;
540 else { // N-side is pos clust
541 dLorentz = fLorentzShiftP;
544 Int_t n = nStrips[adc][side];
545 for( int istr = 0; istr<n+1; istr++ ){
549 Float_t signal=0.0, noise=0.0, gain=0.0;
552 strip = strips[adc][side][istr][0];
553 signal = strips[adc][side][istr][1];
555 //cout<<"strip "<<adc<<" / "<<side<<": "<<strip<<endl;
558 noise = side ?cal->GetNoiseN(strip) :cal->GetNoiseP(strip);
559 gain = side ?cal->GetGainN(strip) :cal->GetGainP(strip);
560 stripOK = ( noise>=1. && signal>=3.0*noise
561 //&& !cal->IsPChannelBad(strip)
564 } else stripOK = 0; // end of data
566 bool newCluster = ( (abs(strip-ostrip)!=1) || !stripOK );
570 //* Store the previous cluster
572 if( nDigits>0 && q>0 && snFlag ){
574 if (nClusters1D[side] >= kMaxADCClusters-1 ) {
575 AliWarning("HLT ClustersFinderSSD: Too many 1D clusters !");
578 Ali1Dcluster &cluster = clusters1D[side][nClusters1D[side]++];
579 cluster.SetY( y / q + dStrip + dLorentz);
581 cluster.SetNd(nDigits);
582 cluster.SetLabels(lab);
583 //cout<<"cluster 1D side "<<side<<": y= "<<y<<" q= "<<q<<" d="<<dStrip<<" Y="<<cluster.GetY()<<endl;
584 //Split suspiciously big cluster
586 if( repa->GetUseUnfoldingInClusterFinderSSD()
587 && nDigits > 4 && nDigits < 25
589 cluster.SetY(y/q + dStrip - 0.25*nDigits + dLorentz);
591 Ali1Dcluster& cluster2 = clusters1D[side][nClusters1D[side]++];
592 cluster2.SetY(y/q + dStrip + 0.25*nDigits + dLorentz);
593 cluster2.SetQ(0.5*q);
594 cluster2.SetNd(nDigits);
595 cluster2.SetLabels(lab);
603 } //* End store the previous cluster
605 if( stripOK ){ // add new signal to the cluster
606 // signal = (Int_t) (signal * gain); // signal is corrected for gain
607 if( signal>fgkThreshold*noise) snFlag = 1;
608 signal = signal * gain; // signal is corrected for gain
609 // if( cal ) signal = (Int_t) cal->ADCToKeV( signal ); // signal is converted in KeV
610 if( cal ) signal = cal->ADCToKeV( signal ); // signal is converted in KeV
611 q += signal; // add digit to current cluster
618 } //* end loop over strips
620 } //* end loop over ADC sides
624 if( nClusters1D[0] && nClusters1D[1] && fModule>=0 ){
625 TClonesArray* clusters = rpc->UncheckedGetClusters(fModule);
626 FindClustersSSD( clusters1D[0], nClusters1D[0], clusters1D[1], nClusters1D[1], clusters);
627 Int_t nClustersn = clusters->GetEntriesFast();
628 nClustersSSD += nClustersn;
631 //cout<<"SG: "<<ddl<<" "<<ad<<" "<<adc<<": strips "<<nstat[0]<<"+"<<nstat[1]<<", clusters 1D= "<<nClusters1D[0]<<" + "<<nClusters1D[1]<<", 2D= "<<clusters.size()<<endl;
633 }//* end loop over adc
635 }//* end of reconstruction of previous block of 12 modules
640 //* Clean up arrays and set new module
643 for( int i=0; i<kNADC; i++ ){
653 //* Exit main loop when there is no more input
659 //* Fill the current strip information
662 Int_t adc = input->GetADC();
663 if( adc<0 || adc>=kNADC+2 || (adc>5&&adc<8) ){
664 AliWarning(Form("HLT ClustersFinderSSD: Corrupted data: wrong adc number (%d)", adc));
668 if( adc>7 ) adc-= 2; // shift ADC numbers 8-13 to 6-11
670 Bool_t side = input->GetSideFlag();
671 Int_t strip = input->GetStrip();
672 Int_t signal = input->GetSignal();
675 //cout<<"SSD: "<<ddl<<" "<<ad<<" "<<adc<<" "<<side<<" "<<strip<<" : "<<signal<<endl;
678 AliWarning(Form("HLT ClustersFinderSSD: Corrupted data: wrong strip number (ddl %d ad %d adc %d side %d, strip %d",
679 ddl, ad, adc, side,strip) );
682 if (strip < 0) continue;
684 int &n = nStrips[adc][side];
686 Int_t oldStrip = strips[adc][side][n-1][0];
688 if( strip==oldStrip ){
689 AliWarning(Form("HLT ClustersFinderSSD: Corrupted data: duplicated signal: ddl %d ad %d adc %d, side %d, strip %d",
690 ddl, ad, adc, side, strip ));
694 strips[adc][side][n][0] = strip;
695 strips[adc][side][n][1] = signal;
698 //cout<<"SSD: "<<input->GetDDL()<<" "<<input->GetAD()<<" "
699 //<<input->GetADC()<<" "<<input->GetSideFlag()<<" "<<((int)input->GetStrip())<<" "<<strip<<" : "<<input->GetSignal()<<endl;
701 } //* End main loop over the input
703 AliDebug(1,Form("found clusters in ITS SSD: %d", nClustersSSD));
707 void AliITSClusterFinderV2SSD::
708 FindClustersSSD(Ali1Dcluster* neg, Int_t nn,
709 Ali1Dcluster* pos, Int_t np,
710 TClonesArray *clusters) {
711 //------------------------------------------------------------
712 // Actual SSD cluster finder
713 //------------------------------------------------------------
715 const TGeoHMatrix *mT2L=AliITSgeomTGeo::GetTracking2LocalMatrix(fModule);
717 //---------------------------------------
720 static AliITSRecoParam *repa = NULL;
722 repa = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
724 repa = AliITSRecoParam::GetHighFluxParam();
725 AliWarning("Using default AliITSRecoParam class");
729 // TClonesArray &cl=*clusters;
731 AliITSsegmentationSSD *seg = dynamic_cast<AliITSsegmentationSSD*>(fDetTypeRec->GetSegmentationModel(2));
732 if (fModule>fLastSSD1)
737 Float_t hwSSD = seg->Dx()*1e-4/2;
738 Float_t hlSSD = seg->Dz()*1e-4/2;
740 Int_t idet=fNdet[fModule];
744 Int_t *cnegative = new Int_t[np];
745 Int_t *cused1 = new Int_t[np];
746 Int_t *negativepair = new Int_t[10*np];
747 Int_t *cpositive = new Int_t[nn];
748 Int_t *cused2 = new Int_t[nn];
749 Int_t *positivepair = new Int_t[10*nn];
750 for (Int_t i=0;i<np;i++) {cnegative[i]=0; cused1[i]=0;}
751 for (Int_t i=0;i<nn;i++) {cpositive[i]=0; cused2[i]=0;}
752 for (Int_t i=0;i<10*np;i++) {negativepair[i]=0;}
753 for (Int_t i=0;i<10*nn;i++) {positivepair[i]=0;}
755 if ((np*nn) > fgPairsSize) {
757 if (fgPairs) delete [] fgPairs;
758 fgPairsSize = 4*np*nn;
759 fgPairs = new Short_t[fgPairsSize];
761 memset(fgPairs,0,sizeof(Short_t)*np*nn);
764 // find available pairs
767 for (Int_t i=0; i<np; i++) {
768 Float_t yp=pos[i].GetY();
769 if ( (pos[i].GetQ()>0) && (pos[i].GetQ()<3) ) continue;
770 for (Int_t j=0; j<nn; j++) {
771 if ( (neg[j].GetQ()>0) && (neg[j].GetQ()<3) ) continue;
772 Float_t yn=neg[j].GetY();
775 seg->GetPadCxz(yn, yp, xt, zt);
776 //cout<<yn<<" "<<yp<<" "<<xt<<" "<<zt<<endl;
778 if (TMath::Abs(xt)<hwSSD)
779 if (TMath::Abs(zt)<hlSSD) {
780 Int_t in = i*10+cnegative[i];
781 Int_t ip = j*10+cpositive[j];
782 if ((in < 10*np) && (ip < 10*nn)) {
783 negativepair[in] =j; //index
785 cnegative[i]++; //counters
791 AliError(Form("Index out of range: ip=%d, in=%d",ip,in));
799 delete [] negativepair;
802 delete [] positivepair;
805 //why not to allocate memorey here? if(!clusters) clusters = new TClonesArray("AliITSRecPoint", ncross);
808 // try to recover points out of but close to the module boundaries
810 for (Int_t i=0; i<np; i++) {
811 Float_t yp=pos[i].GetY();
812 if ( (pos[i].GetQ()>0) && (pos[i].GetQ()<3) ) continue;
813 for (Int_t j=0; j<nn; j++) {
814 if ( (neg[j].GetQ()>0) && (neg[j].GetQ()<3) ) continue;
815 // if both 1Dclusters have an other cross continue
816 if (cpositive[j]&&cnegative[i]) continue;
817 Float_t yn=neg[j].GetY();
820 seg->GetPadCxz(yn, yp, xt, zt);
822 if (TMath::Abs(xt)<hwSSD+0.1)
823 if (TMath::Abs(zt)<hlSSD+0.15) {
824 // tag 1Dcluster (eventually will produce low quality recpoint)
825 if (cnegative[i]==0) pos[i].SetNd(100); // not available pair
826 if (cpositive[j]==0) neg[j].SetNd(100); // not available pair
827 Int_t in = i*10+cnegative[i];
828 Int_t ip = j*10+cpositive[j];
829 if ((in < 10*np) && (ip < 10*nn)) {
830 negativepair[in] =j; //index
832 cnegative[i]++; //counters
837 AliError(Form("Index out of range: ip=%d, in=%d",ip,in));
849 if(repa->GetUseChargeMatchingInClusterFinderSSD()==kTRUE) {
855 for (Int_t ip=0;ip<np;ip++){
856 Float_t xbest=1000,zbest=1000,qbest=0;
858 // select gold clusters
859 if ( (cnegative[ip]==1) && cpositive[negativepair[10*ip]]==1){
860 Float_t yp=pos[ip].GetY();
861 Int_t j = negativepair[10*ip];
863 if( (pos[ip].GetQ()==0) && (neg[j].GetQ() ==0) ) {
864 // both bad, hence continue;
865 // mark both as used (to avoid recover at the end)
871 ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
872 //cout<<"ratio="<<ratio<<endl;
874 // charge matching (note that if posQ or negQ is 0 -> ratio=1 and the following condition is met
875 if (TMath::Abs(ratio)>0.2) continue; // note: 0.2=3xsigma_ratio calculated in cosmics tests
878 Float_t yn=neg[j].GetY();
881 seg->GetPadCxz(yn, yp, xt, zt);
886 qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
887 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
890 Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
891 mT2L->MasterToLocal(loc,trk);
896 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
897 for (Int_t ilab=0;ilab<3;ilab++){
898 milab[ilab] = pos[ip].GetLabel(ilab);
899 milab[ilab+3] = neg[j].GetLabel(ilab);
903 milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
904 Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
906 lp[2]=4.968e-06; // 0.00223*0.00223; //SigmaY2
907 lp[3]=0.012; // 0.110*0.110; //SigmaZ2
908 // out-of-diagonal element of covariance matrix
909 if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
910 else if ( (info[0]>1) && (info[1]>1) ) {
911 lp[2]=2.63e-06; // 0.0016*0.0016; //SigmaY2
912 lp[3]=0.0065; // 0.08*0.08; //SigmaZ2
916 lp[2]=4.80e-06; // 0.00219*0.00219
917 lp[3]=0.0093; // 0.0964*0.0964;
922 lp[2]=2.79e-06; // 0.0017*0.0017;
923 lp[3]=0.00935; // 0.967*0.967;
928 AliITSRecPoint * cl2;
929 cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
931 cl2->SetChargeRatio(ratio);
935 if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
940 if(pos[ip].GetQ()==0) cl2->SetType(3);
941 if(neg[j].GetQ()==0) cl2->SetType(4);
950 for (Int_t ip=0;ip<np;ip++){
951 Float_t xbest=1000,zbest=1000,qbest=0;
954 // select "silber" cluster
955 if ( cnegative[ip]==1 && cpositive[negativepair[10*ip]]==2){
956 Int_t in = negativepair[10*ip];
957 Int_t ip2 = positivepair[10*in];
958 if (ip2==ip) ip2 = positivepair[10*in+1];
959 Float_t pcharge = pos[ip].GetQ()+pos[ip2].GetQ();
963 ratio = (pcharge-neg[in].GetQ())/(pcharge+neg[in].GetQ());
964 if ( (TMath::Abs(ratio)<0.2) && (pcharge!=0) ) {
965 //if ( (TMath::Abs(pcharge-neg[in].GetQ())<30) && (pcharge!=0) ) { //
969 if ( (fgPairs[ip*nn+in]==100)&&(pos[ip].GetQ() ) ) { //
971 Float_t yp=pos[ip].GetY();
972 Float_t yn=neg[in].GetY();
975 seg->GetPadCxz(yn, yp, xt, zt);
979 qbest =pos[ip].GetQ();
980 Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
981 mT2L->MasterToLocal(loc,trk);
986 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
987 for (Int_t ilab=0;ilab<3;ilab++){
988 milab[ilab] = pos[ip].GetLabel(ilab);
989 milab[ilab+3] = neg[in].GetLabel(ilab);
993 ratio = (pos[ip].GetQ()-neg[in].GetQ())/(pos[ip].GetQ()+neg[in].GetQ());
994 milab[3]=(((ip<<10) + in)<<10) + idet; // pos|neg|det
995 Int_t info[3] = {pos[ip].GetNd(),neg[in].GetNd(),fNlayer[fModule]};
997 lp[2]=4.968e-06; // 0.00223*0.00223; //SigmaY2
998 lp[3]=0.012; // 0.110*0.110; //SigmaZ2
999 // out-of-diagonal element of covariance matrix
1000 if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1001 else if ( (info[0]>1) && (info[1]>1) ) {
1002 lp[2]=2.63e-06; // 0.0016*0.0016; //SigmaY2
1003 lp[3]=0.0065; // 0.08*0.08; //SigmaZ2
1007 lp[2]=4.80e-06; // 0.00219*0.00219
1008 lp[3]=0.0093; // 0.0964*0.0964;
1013 lp[2]=2.79e-06; // 0.0017*0.0017;
1014 lp[3]=0.00935; // 0.967*0.967;
1019 AliITSRecPoint * cl2;
1020 cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1021 cl2->SetChargeRatio(ratio);
1023 fgPairs[ip*nn+in] = 5;
1024 if ((pos[ip].GetNd()+neg[in].GetNd())>6){ //multi cluster
1026 fgPairs[ip*nn+in] = 6;
1035 // if (!(cused1[ip2] || cused2[in])){ //
1036 if ( (fgPairs[ip2*nn+in]==100) && (pos[ip2].GetQ()) ) {
1038 Float_t yp=pos[ip2].GetY();
1039 Float_t yn=neg[in].GetY();
1042 seg->GetPadCxz(yn, yp, xt, zt);
1046 qbest =pos[ip2].GetQ();
1048 Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
1049 mT2L->MasterToLocal(loc,trk);
1054 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1055 for (Int_t ilab=0;ilab<3;ilab++){
1056 milab[ilab] = pos[ip2].GetLabel(ilab);
1057 milab[ilab+3] = neg[in].GetLabel(ilab);
1060 CheckLabels2(milab);
1061 ratio = (pos[ip2].GetQ()-neg[in].GetQ())/(pos[ip2].GetQ()+neg[in].GetQ());
1062 milab[3]=(((ip2<<10) + in)<<10) + idet; // pos|neg|det
1063 Int_t info[3] = {pos[ip2].GetNd(),neg[in].GetNd(),fNlayer[fModule]};
1065 lp[2]=4.968e-06; // 0.00223*0.00223; //SigmaY2
1066 lp[3]=0.012; // 0.110*0.110; //SigmaZ2
1067 // out-of-diagonal element of covariance matrix
1068 if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1069 else if ( (info[0]>1) && (info[1]>1) ) {
1070 lp[2]=2.63e-06; // 0.0016*0.0016; //SigmaY2
1071 lp[3]=0.0065; // 0.08*0.08; //SigmaZ2
1075 lp[2]=4.80e-06; // 0.00219*0.00219
1076 lp[3]=0.0093; // 0.0964*0.0964;
1081 lp[2]=2.79e-06; // 0.0017*0.0017;
1082 lp[3]=0.00935; // 0.967*0.967;
1087 AliITSRecPoint * cl2;
1088 cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1090 cl2->SetChargeRatio(ratio);
1092 fgPairs[ip2*nn+in] =5;
1093 if ((pos[ip2].GetNd()+neg[in].GetNd())>6){ //multi cluster
1095 fgPairs[ip2*nn+in] =6;
1104 } // charge matching condition
1106 } // 2 Pside cross 1 Nside
1107 } // loop over Pside clusters
1112 for (Int_t jn=0;jn<nn;jn++){
1113 if (cused2[jn]) continue;
1114 Float_t xbest=1000,zbest=1000,qbest=0;
1115 // select "silber" cluster
1116 if ( cpositive[jn]==1 && cnegative[positivepair[10*jn]]==2){
1117 Int_t ip = positivepair[10*jn];
1118 Int_t jn2 = negativepair[10*ip];
1119 if (jn2==jn) jn2 = negativepair[10*ip+1];
1120 Float_t pcharge = neg[jn].GetQ()+neg[jn2].GetQ();
1124 ratio = (pcharge-pos[ip].GetQ())/(pcharge+pos[ip].GetQ());
1125 if ( (TMath::Abs(ratio)<0.2) && (pcharge!=0) ) {
1128 if ( (TMath::Abs(pcharge-pos[ip].GetQ())<30) && // charge matching
1129 (pcharge!=0) ) { // reject combinations of bad strips
1135 // if (!(cused1[ip]||cused2[jn])){
1136 if ( (fgPairs[ip*nn+jn]==100) && (neg[jn].GetQ()) ) { //
1138 Float_t yn=neg[jn].GetY();
1139 Float_t yp=pos[ip].GetY();
1142 seg->GetPadCxz(yn, yp, xt, zt);
1146 qbest =neg[jn].GetQ();
1149 Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
1150 mT2L->MasterToLocal(loc,trk);
1156 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1157 for (Int_t ilab=0;ilab<3;ilab++){
1158 milab[ilab] = pos[ip].GetLabel(ilab);
1159 milab[ilab+3] = neg[jn].GetLabel(ilab);
1162 CheckLabels2(milab);
1163 ratio = (pos[ip].GetQ()-neg[jn].GetQ())/(pos[ip].GetQ()+neg[jn].GetQ());
1164 milab[3]=(((ip<<10) + jn)<<10) + idet; // pos|neg|det
1165 Int_t info[3] = {pos[ip].GetNd(),neg[jn].GetNd(),fNlayer[fModule]};
1167 lp[2]=4.968e-06; // 0.00223*0.00223; //SigmaY2
1168 lp[3]=0.012; // 0.110*0.110; //SigmaZ2
1169 // out-of-diagonal element of covariance matrix
1170 if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1171 else if ( (info[0]>1) && (info[1]>1) ) {
1172 lp[2]=2.63e-06; // 0.0016*0.0016; //SigmaY2
1173 lp[3]=0.0065; // 0.08*0.08; //SigmaZ2
1177 lp[2]=4.80e-06; // 0.00219*0.00219
1178 lp[3]=0.0093; // 0.0964*0.0964;
1183 lp[2]=2.79e-06; // 0.0017*0.0017;
1184 lp[3]=0.00935; // 0.967*0.967;
1189 AliITSRecPoint * cl2;
1190 cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1192 cl2->SetChargeRatio(ratio);
1194 fgPairs[ip*nn+jn] =7;
1195 if ((pos[ip].GetNd()+neg[jn].GetNd())>6){ //multi cluster
1197 fgPairs[ip*nn+jn]=8;
1203 // if (!(cused1[ip]||cused2[jn2])){
1204 if ( (fgPairs[ip*nn+jn2]==100)&&(neg[jn2].GetQ() ) ) { //
1206 Float_t yn=neg[jn2].GetY();
1207 Double_t yp=pos[ip].GetY();
1210 seg->GetPadCxz(yn, yp, xt, zt);
1214 qbest =neg[jn2].GetQ();
1217 Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
1218 mT2L->MasterToLocal(loc,trk);
1224 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1225 for (Int_t ilab=0;ilab<3;ilab++){
1226 milab[ilab] = pos[ip].GetLabel(ilab);
1227 milab[ilab+3] = neg[jn2].GetLabel(ilab);
1230 CheckLabels2(milab);
1231 ratio = (pos[ip].GetQ()-neg[jn2].GetQ())/(pos[ip].GetQ()+neg[jn2].GetQ());
1232 milab[3]=(((ip<<10) + jn2)<<10) + idet; // pos|neg|det
1233 Int_t info[3] = {pos[ip].GetNd(),neg[jn2].GetNd(),fNlayer[fModule]};
1235 lp[2]=4.968e-06; // 0.00223*0.00223; //SigmaY2
1236 lp[3]=0.012; // 0.110*0.110; //SigmaZ2
1237 // out-of-diagonal element of covariance matrix
1238 if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1239 else if ( (info[0]>1) && (info[1]>1) ) {
1240 lp[2]=2.63e-06; // 0.0016*0.0016; //SigmaY2
1241 lp[3]=0.0065; // 0.08*0.08; //SigmaZ2
1245 lp[2]=4.80e-06; // 0.00219*0.00219
1246 lp[3]=0.0093; // 0.0964*0.0964;
1251 lp[2]=2.79e-06; // 0.0017*0.0017;
1252 lp[3]=0.00935; // 0.967*0.967;
1257 AliITSRecPoint * cl2;
1258 cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1261 cl2->SetChargeRatio(ratio);
1262 fgPairs[ip*nn+jn2]=7;
1264 if ((pos[ip].GetNd()+neg[jn2].GetNd())>6){ //multi cluster
1266 fgPairs[ip*nn+jn2]=8;
1274 } // charge matching condition
1276 } // 2 Nside cross 1 Pside
1277 } // loop over Pside clusters
1281 for (Int_t ip=0;ip<np;ip++){
1283 if(cused1[ip]) continue;
1286 Float_t xbest=1000,zbest=1000,qbest=0;
1290 if ( (cnegative[ip]==2) && cpositive[negativepair[10*ip]]==2){
1291 Float_t minchargediff =4.;
1292 Float_t minchargeratio =0.2;
1295 for (Int_t di=0;di<cnegative[ip];di++){
1296 Int_t jc = negativepair[ip*10+di];
1297 Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
1298 ratio = (pos[ip].GetQ()-neg[jc].GetQ())/(pos[ip].GetQ()+neg[jc].GetQ());
1299 //if (TMath::Abs(chargedif)<minchargediff){
1300 if (TMath::Abs(ratio)<0.2){
1302 minchargediff = TMath::Abs(chargedif);
1303 minchargeratio = TMath::Abs(ratio);
1306 if (j<0) continue; // not proper cluster
1310 for (Int_t di=0;di<cnegative[ip];di++){
1311 Int_t jc = negativepair[ip*10+di];
1312 Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
1313 if (TMath::Abs(chargedif)<minchargediff+3.) count++;
1315 if (count>1) continue; // more than one "proper" cluster for positive
1319 for (Int_t dj=0;dj<cpositive[j];dj++){
1320 Int_t ic = positivepair[j*10+dj];
1321 Float_t chargedif = pos[ic].GetQ()-neg[j].GetQ();
1322 if (TMath::Abs(chargedif)<minchargediff+3.) count++;
1324 if (count>1) continue; // more than one "proper" cluster for negative
1329 for (Int_t dj=0;dj<cnegative[jp];dj++){
1330 Int_t ic = positivepair[jp*10+dj];
1331 Float_t chargedif = pos[ic].GetQ()-neg[jp].GetQ();
1332 if (TMath::Abs(chargedif)<minchargediff+4.) count++;
1334 if (count>1) continue;
1335 if (fgPairs[ip*nn+j]<100) continue;
1340 //almost gold clusters
1341 Float_t yp=pos[ip].GetY();
1342 Float_t yn=neg[j].GetY();
1344 seg->GetPadCxz(yn, yp, xt, zt);
1346 qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
1348 Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
1349 mT2L->MasterToLocal(loc,trk);
1354 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1355 for (Int_t ilab=0;ilab<3;ilab++){
1356 milab[ilab] = pos[ip].GetLabel(ilab);
1357 milab[ilab+3] = neg[j].GetLabel(ilab);
1360 CheckLabels2(milab);
1361 if ((neg[j].GetQ()==0)&&(pos[ip].GetQ()==0)) continue; // reject crosses of bad strips!!
1362 ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
1363 milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
1364 Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
1366 lp[2]=4.968e-06; // 0.00223*0.00223; //SigmaY2
1367 lp[3]=0.012; // 0.110*0.110; //SigmaZ2
1368 // out-of-diagonal element of covariance matrix
1369 if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1370 else if ( (info[0]>1) && (info[1]>1) ) {
1371 lp[2]=2.63e-06; // 0.0016*0.0016; //SigmaY2
1372 lp[3]=0.0065; // 0.08*0.08; //SigmaZ2
1376 lp[2]=4.80e-06; // 0.00219*0.00219
1377 lp[3]=0.0093; // 0.0964*0.0964;
1382 lp[2]=2.79e-06; // 0.0017*0.0017;
1383 lp[3]=0.00935; // 0.967*0.967;
1388 AliITSRecPoint * cl2;
1389 cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1391 cl2->SetChargeRatio(ratio);
1393 fgPairs[ip*nn+j]=10;
1394 if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
1396 fgPairs[ip*nn+j]=11;
1403 } // loop over Pside 1Dclusters
1407 for (Int_t ip=0;ip<np;ip++){
1409 if(cused1[ip]) continue;
1412 Float_t xbest=1000,zbest=1000,qbest=0;
1414 // manyxmany clusters
1416 if ( (cnegative[ip]<5) && cpositive[negativepair[10*ip]]<5){
1417 Float_t minchargediff =4.;
1419 for (Int_t di=0;di<cnegative[ip];di++){
1420 Int_t jc = negativepair[ip*10+di];
1421 Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
1422 if (TMath::Abs(chargedif)<minchargediff){
1424 minchargediff = TMath::Abs(chargedif);
1427 if (j<0) continue; // not proper cluster
1430 for (Int_t di=0;di<cnegative[ip];di++){
1431 Int_t jc = negativepair[ip*10+di];
1432 Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
1433 if (TMath::Abs(chargedif)<minchargediff+3.) count++;
1435 if (count>1) continue; // more than one "proper" cluster for positive
1439 for (Int_t dj=0;dj<cpositive[j];dj++){
1440 Int_t ic = positivepair[j*10+dj];
1441 Float_t chargedif = pos[ic].GetQ()-neg[j].GetQ();
1442 if (TMath::Abs(chargedif)<minchargediff+3.) count++;
1444 if (count>1) continue; // more than one "proper" cluster for negative
1449 for (Int_t dj=0;dj<cnegative[jp];dj++){
1450 Int_t ic = positivepair[jp*10+dj];
1451 Float_t chargedif = pos[ic].GetQ()-neg[jp].GetQ();
1452 if (TMath::Abs(chargedif)<minchargediff+4.) count++;
1454 if (count>1) continue;
1455 if (fgPairs[ip*nn+j]<100) continue;
1458 //almost gold clusters
1459 Float_t yp=pos[ip].GetY();
1460 Float_t yn=neg[j].GetY();
1464 seg->GetPadCxz(yn, yp, xt, zt);
1468 qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
1471 Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
1472 mT2L->MasterToLocal(loc,trk);
1477 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1478 for (Int_t ilab=0;ilab<3;ilab++){
1479 milab[ilab] = pos[ip].GetLabel(ilab);
1480 milab[ilab+3] = neg[j].GetLabel(ilab);
1483 CheckLabels2(milab);
1484 if ((neg[j].GetQ()==0)&&(pos[ip].GetQ()==0)) continue; // reject crosses of bad strips!!
1485 ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
1486 milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
1487 Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
1489 lp[2]=4.968e-06; // 0.00223*0.00223; //SigmaY2
1490 lp[3]=0.012; // 0.110*0.110; //SigmaZ2
1491 // out-of-diagonal element of covariance matrix
1492 if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1493 else if ( (info[0]>1) && (info[1]>1) ) {
1494 lp[2]=2.63e-06; // 0.0016*0.0016; //SigmaY2
1495 lp[3]=0.0065; // 0.08*0.08; //SigmaZ2
1499 lp[2]=4.80e-06; // 0.00219*0.00219
1500 lp[3]=0.0093; // 0.0964*0.0964;
1505 lp[2]=2.79e-06; // 0.0017*0.0017;
1506 lp[3]=0.00935; // 0.967*0.967;
1511 AliITSRecPoint * cl2;
1512 cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1514 cl2->SetChargeRatio(ratio);
1516 fgPairs[ip*nn+j]=12;
1517 if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
1519 fgPairs[ip*nn+j]=13;
1526 } // loop over Pside 1Dclusters
1528 } // use charge matching
1531 // recover all the other crosses
1533 for (Int_t i=0; i<np; i++) {
1534 Float_t xbest=1000,zbest=1000,qbest=0;
1535 Float_t yp=pos[i].GetY();
1536 if ((pos[i].GetQ()>0)&&(pos[i].GetQ()<3)) continue;
1537 for (Int_t j=0; j<nn; j++) {
1538 // for (Int_t di = 0;di<cpositive[i];di++){
1539 // Int_t j = negativepair[10*i+di];
1540 if ((neg[j].GetQ()>0)&&(neg[j].GetQ()<3)) continue;
1542 if ((neg[j].GetQ()==0)&&(pos[i].GetQ()==0)) continue; // reject crosses of bad strips!!
1544 if (cused2[j]||cused1[i]) continue;
1545 if (fgPairs[i*nn+j]>0 &&fgPairs[i*nn+j]<100) continue;
1546 ratio = (pos[i].GetQ()-neg[j].GetQ())/(pos[i].GetQ()+neg[j].GetQ());
1547 Float_t yn=neg[j].GetY();
1550 seg->GetPadCxz(yn, yp, xt, zt);
1552 if (TMath::Abs(xt)<hwSSD)
1553 if (TMath::Abs(zt)<hlSSD) {
1556 qbest=0.5*(pos[i].GetQ()+neg[j].GetQ());
1559 Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
1560 mT2L->MasterToLocal(loc,trk);
1565 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1566 for (Int_t ilab=0;ilab<3;ilab++){
1567 milab[ilab] = pos[i].GetLabel(ilab);
1568 milab[ilab+3] = neg[j].GetLabel(ilab);
1571 CheckLabels2(milab);
1572 milab[3]=(((i<<10) + j)<<10) + idet; // pos|neg|det
1573 Int_t info[3] = {pos[i].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
1575 lp[2]=4.968e-06; // 0.00223*0.00223; //SigmaY2
1576 lp[3]=0.012; // 0.110*0.110; //SigmaZ2
1577 // out-of-diagonal element of covariance matrix
1578 if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1579 else if ( (info[0]>1) && (info[1]>1) ) {
1580 lp[2]=2.63e-06; // 0.0016*0.0016; //SigmaY2
1581 lp[3]=0.0065; // 0.08*0.08; //SigmaZ2
1585 lp[2]=4.80e-06; // 0.00219*0.00219
1586 lp[3]=0.0093; // 0.0964*0.0964;
1591 lp[2]=2.79e-06; // 0.0017*0.0017;
1592 lp[3]=0.00935; // 0.967*0.967;
1597 AliITSRecPoint * cl2;
1598 cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1600 cl2->SetChargeRatio(ratio);
1601 cl2->SetType(100+cpositive[j]+cnegative[i]);
1603 if(pos[i].GetQ()==0) cl2->SetType(200+cpositive[j]+cnegative[i]);
1604 if(neg[j].GetQ()==0) cl2->SetType(300+cpositive[j]+cnegative[i]);
1612 if(repa->GetUseBadChannelsInClusterFinderSSD()==kTRUE) {
1614 //---------------------------------------------------------
1615 // recover crosses of good 1D clusters with bad strips on the other side
1616 // Note1: at first iteration skip modules with a bad side (or almost), (would produce too many fake!)
1617 // Note2: for modules with a bad side see below
1619 AliITSCalibrationSSD* cal = (AliITSCalibrationSSD*)GetResp(fModule);
1620 Int_t countPbad=0, countNbad=0;
1621 for(Int_t ib=0; ib<768; ib++) {
1622 if(cal->IsPChannelBad(ib)) countPbad++;
1623 if(cal->IsNChannelBad(ib)) countNbad++;
1625 // AliInfo(Form("module %d has %d P- and %d N-bad strips",fModule,countPbad,countNbad));
1627 if( (countPbad<100) && (countNbad<100) ) { // no bad side!!
1629 for (Int_t i=0; i<np; i++) { // loop over Nside 1Dclusters with no crosses
1630 if(cnegative[i]) continue; // if intersecting Pside clusters continue;
1632 // for(Int_t ib=0; ib<768; ib++) { // loop over all Pstrips
1633 for(Int_t ib=15; ib<753; ib++) { // loop over all Pstrips
1635 if(cal->IsPChannelBad(ib)) { // check if strips is bad
1636 Float_t yN=pos[i].GetY();
1638 seg->GetPadCxz(1.*ib, yN, xt, zt);
1641 // bad Pstrip is crossing the Nside 1Dcluster -> create recpoint
1643 if ( (TMath::Abs(xt)<hwSSD) && (TMath::Abs(zt)<hlSSD) ) {
1644 Double_t loc[3]={xt,0.,zt},trk[3]={0.,0.,0.};
1645 mT2L->MasterToLocal(loc,trk);
1648 lp[4]=pos[i].GetQ(); //Q
1649 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1650 for (Int_t ilab=0;ilab<3;ilab++) milab[ilab] = pos[i].GetLabel(ilab);
1651 CheckLabels2(milab);
1652 milab[3]=( (i<<10) << 10 ) + idet; // pos|neg|det
1653 Int_t info[3] = {pos[i].GetNd(),0,fNlayer[fModule]};
1655 lp[2]=4.968e-06; // 0.00223*0.00223; //SigmaY2
1656 lp[3]=0.012; // 0.110*0.110; //SigmaZ2
1657 lp[5]=-0.00012; // out-of-diagonal element of covariance matrix
1664 AliITSRecPoint * cl2;
1665 cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1666 cl2->SetChargeRatio(1.);
1669 } // cross is within the detector
1675 } // end loop over Pstrips
1677 } // end loop over Nside 1D clusters
1679 for (Int_t j=0; j<nn; j++) { // loop over Pside 1D clusters with no crosses
1680 if(cpositive[j]) continue;
1682 // for(Int_t ib=0; ib<768; ib++) { // loop over all Nside strips
1683 for(Int_t ib=15; ib<753; ib++) { // loop over all Nside strips
1685 if(cal->IsNChannelBad(ib)) { // check if strip is bad
1686 Float_t yP=neg[j].GetY();
1688 seg->GetPadCxz(yP, 1.*ib, xt, zt);
1691 // bad Nstrip is crossing the Pside 1Dcluster -> create recpoint
1693 if ( (TMath::Abs(xt)<hwSSD) && (TMath::Abs(zt)<hlSSD) ) {
1694 Double_t loc[3]={xt,0.,zt},trk[3]={0.,0.,0.};
1695 mT2L->MasterToLocal(loc,trk);
1698 lp[4]=neg[j].GetQ(); //Q
1699 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1700 for (Int_t ilab=0;ilab<3;ilab++) milab[ilab] = neg[j].GetLabel(ilab);
1701 CheckLabels2(milab);
1702 milab[3]=( j << 10 ) + idet; // pos|neg|det
1703 Int_t info[3]={0,(Int_t)neg[j].GetNd(),fNlayer[fModule]};
1705 lp[2]=4.968e-06; // 0.00223*0.00223; //SigmaY2
1706 lp[3]=0.012; // 0.110*0.110; //SigmaZ2
1707 lp[5]=-0.00012; // out-of-diagonal element of covariance matrix
1714 AliITSRecPoint * cl2;
1715 cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1716 cl2->SetChargeRatio(1.);
1719 } // cross is within the detector
1724 } // end loop over Nstrips
1725 } // end loop over Pside 1D clusters
1729 //---------------------------------------------------------
1731 else if( (countPbad>700) && (countNbad<100) ) { // bad Pside!!
1733 for (Int_t i=0; i<np; i++) { // loop over Nside 1Dclusters with no crosses
1734 if(cnegative[i]) continue; // if intersecting Pside clusters continue;
1737 Float_t yN=pos[i].GetY();
1739 if (seg->GetLayer()==5) yP = yN + (7.6/1.9);
1740 else yP = yN - (7.6/1.9);
1741 seg->GetPadCxz(yP, yN, xt, zt);
1743 if ( (TMath::Abs(xt)<hwSSD) && (TMath::Abs(zt)<hlSSD) ) {
1744 Double_t loc[3]={xt,0.,zt},trk[3]={0.,0.,0.};
1745 mT2L->MasterToLocal(loc,trk);
1748 lp[4]=pos[i].GetQ(); //Q
1749 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1750 for (Int_t ilab=0;ilab<3;ilab++) milab[ilab] = pos[i].GetLabel(ilab);
1751 CheckLabels2(milab);
1752 milab[3]=( (i<<10) << 10 ) + idet; // pos|neg|det
1753 Int_t info[3] = {(Int_t)pos[i].GetNd(),0,fNlayer[fModule]};
1755 lp[2]=0.00098; // 0.031*0.031; //SigmaY2
1756 lp[3]=1.329; // 1.15*1.15; //SigmaZ2
1758 if(info[0]>1) lp[2]=0.00097;
1760 AliITSRecPoint * cl2;
1761 cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1762 cl2->SetChargeRatio(1.);
1765 } // cross is within the detector
1769 } // end loop over Nside 1D clusters
1771 } // bad Pside module
1773 else if( (countNbad>700) && (countPbad<100) ) { // bad Nside!!
1775 for (Int_t j=0; j<nn; j++) { // loop over Pside 1D clusters with no crosses
1776 if(cpositive[j]) continue;
1779 Float_t yP=neg[j].GetY();
1781 if (seg->GetLayer()==5) yN = yP - (7.6/1.9);
1782 else yN = yP + (7.6/1.9);
1783 seg->GetPadCxz(yP, yN, xt, zt);
1785 if ( (TMath::Abs(xt)<hwSSD) && (TMath::Abs(zt)<hlSSD) ) {
1786 Double_t loc[3]={xt,0.,zt},trk[3]={0.,0.,0.};
1787 mT2L->MasterToLocal(loc,trk);
1790 lp[4]=neg[j].GetQ(); //Q
1791 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1792 for (Int_t ilab=0;ilab<3;ilab++) milab[ilab] = neg[j].GetLabel(ilab);
1793 CheckLabels2(milab);
1794 milab[3]=( j << 10 ) + idet; // pos|neg|det
1795 Int_t info[3] = {0,(Int_t)neg[j].GetNd(),fNlayer[fModule]};
1797 lp[2]=7.27e-05; // 0.0085*0.0085; //SigmaY2
1798 lp[3]=1.33; // 1.15*1.15; //SigmaZ2
1800 if(info[1]>1) lp[2]=6.91e-05;
1802 AliITSRecPoint * cl2;
1803 cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1804 cl2->SetChargeRatio(1.);
1807 } // cross is within the detector
1811 } // end loop over Pside 1D clusters
1813 } // bad Nside module
1815 //---------------------------------------------------------
1817 } // use bad channels
1819 //cout<<ncl<<" clusters for this module"<<endl;
1821 delete [] cnegative;
1823 delete [] negativepair;
1824 delete [] cpositive;
1826 delete [] positivepair;