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