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