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