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