TRD/AliTRDtrackerV1.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, 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 //                                                                           //
  20 //  Track finder                                                             //
  21 //                                                                           //
  22 //  Authors:                                                                 //
  23 //    Alex Bercuci <A.Bercuci@gsi.de>                                        //
  24 //    Markus Fasel <M.Fasel@gsi.de>                                          //
  25 //                                                                           //
  26 ///////////////////////////////////////////////////////////////////////////////
  27
  28 #include <Riostream.h>
  29 #include <stdio.h>
  30 #include <string.h>
  31
  32 #include <TBranch.h>
  33 #include <TFile.h>
  34 #include <TGraph.h>
  35 #include <TH1D.h>
  36 #include <TH2D.h>
  37 #include <TLinearFitter.h>
  38 #include <TObjArray.h>
  39 #include <TROOT.h>
  40 #include <TTree.h>
  41 #include <TClonesArray.h>
  42 #include <TRandom.h>
  43 #include <TTreeStream.h>
  44
  45 #include "AliLog.h"
  46 #include "AliESDEvent.h"
  47 #include "AliAlignObj.h"
  48 #include "AliRieman.h"
  49 #include "AliTrackPointArray.h"
  50
  51 #include "AliTRDtracker.h"
  52 #include "AliTRDtrackerV1.h"
  53 #include "AliTRDgeometry.h"
  54 #include "AliTRDpadPlane.h"
  55 #include "AliTRDgeometry.h"
  56 #include "AliTRDcluster.h"
  57 #include "AliTRDtrack.h"
  58 #include "AliTRDseed.h"
  59 #include "AliTRDcalibDB.h"
  60 #include "AliTRDCommonParam.h"
  61 #include "AliTRDReconstructor.h"
  62 #include "AliTRDCalibraFillHisto.h"
  63 #include "AliTRDtrackerFitter.h"
  64 #include "AliTRDstackLayer.h"
  65 #include "AliTRDrecoParam.h"
  66 #include "AliTRDseedV1.h"
  67
  68 #define DEBUG
  69
  70 ClassImp(AliTRDtrackerV1)
  71 Double_t AliTRDtrackerV1::fgTopologicQA[kNConfigs] = {
  72                 0.1112, 0.1112, 0.1112, 0.0786, 0.0786,
  73                 0.0786, 0.0786, 0.0579, 0.0579, 0.0474,
  74                 0.0474, 0.0408, 0.0335, 0.0335, 0.0335
  75 };
  76
  77 //____________________________________________________________________
  78 AliTRDtrackerV1::AliTRDtrackerV1(AliTRDrecoParam *p)
  79   :AliTRDtracker()
  80   ,fSieveSeeding(0)
  81   ,fRecoParam(p)
  82   ,fFitter(0x0)
  83   ,fDebugStreamerV1(0x0)
  84 {
  85   //
  86   // Default constructor. Nothing is initialized.
  87   //
  88
  89 }
  90
  91 //____________________________________________________________________
  92 AliTRDtrackerV1::AliTRDtrackerV1(const TFile *in, AliTRDrecoParam *p)
  93   :AliTRDtracker(in)
  94   ,fSieveSeeding(0)
  95   ,fRecoParam(p)
  96   ,fFitter(0x0)
  97   ,fDebugStreamerV1(0x0)
  98 {
  99   //
 100   // Standard constructor.
 101   // Setting of the geometry file, debug output (if enabled)
 102   // and initilize fitter helper.
 103   //
 104
 105         fFitter = new AliTRDtrackerFitter();
 106
 107 #ifdef DEBUG
 108         fDebugStreamerV1 = new TTreeSRedirector("TRDdebug.root");
 109         fFitter->SetDebugStream(fDebugStreamerV1);
 110 #endif
 111
 112 }
 113
 114 //____________________________________________________________________
 115 AliTRDtrackerV1::~AliTRDtrackerV1()
 116 {
 117   //
 118   // Destructor
 119   //
 120
 121         if(fDebugStreamerV1) delete fDebugStreamerV1;
 122         if(fFitter) delete fFitter;
 123         if(fRecoParam) delete fRecoParam;
 124
 125 }
 126
 127 //____________________________________________________________________
 128 Int_t AliTRDtrackerV1::Clusters2Tracks(AliESDEvent *esd)
 129 {
 130   //
 131   // Steering stand alone tracking for full TRD detector
 132   //
 133   // Parameters :
 134   //   esd     : The ESD event. On output it contains
 135   //             the ESD tracks found in TRD.
 136   //
 137   // Output :
 138   //   Number of tracks found in the TRD detector.
 139   //
 140   // Detailed description
 141   // 1. Launch individual SM trackers.
 142   //    See AliTRDtrackerV1::Clusters2TracksSM() for details.
 143   //
 144
 145         if(!fRecoParam){
 146                 AliError("Reconstruction configuration not initialized. Call first AliTRDtrackerV1::SetRecoParam().");
 147                 return 0;
 148         }
 149
 150         //AliInfo("Start Track Finder ...");
 151         Int_t ntracks = 0;
 152         for(int ism=0; ism<AliTRDtracker::kTrackingSectors; ism++){
 153                         //AliInfo(Form("Processing supermodule %i ...", ism));
 154                         ntracks += Clusters2TracksSM(fTrSec[ism], esd);
 155         }
 156         AliInfo(Form("Found %d TRD tracks.", ntracks));
 157         return ntracks;
 158 }
 159
 160 //____________________________________________________________________
 161 Int_t AliTRDtrackerV1::Clusters2TracksSM(AliTRDtracker::AliTRDtrackingSector *sector
 162                                        , AliESDEvent *esd)
 163 {
 164   //
 165   // Steer tracking for one SM.
 166   //
 167   // Parameters :
 168   //   sector  : Array of (SM) propagation layers containing clusters
 169   //   esd     : The current ESD event. On output it contains the also
 170   //             the ESD (TRD) tracks found in this SM.
 171   //
 172   // Output :
 173   //   Number of tracks found in this TRD supermodule.
 174   //
 175   // Detailed description
 176   //
 177   // 1. Unpack AliTRDpropagationLayers objects for each stack.
 178   // 2. Launch stack tracking.
 179   //    See AliTRDtrackerV1::Clusters2TracksStack() for details.
 180   // 3. Pack results in the ESD event.
 181   //
 182
 183         AliTRDpadPlane *pp = 0x0;
 184
 185         // allocate space for esd tracks in this SM
 186         TClonesArray esdTrackList("AliESDtrack", 2*kMaxTracksStack);
 187         esdTrackList.SetOwner();
 188         AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
 189         Int_t nTimeBins = cal->GetNumberOfTimeBins();
 190         const Int_t kFindable = Int_t(fRecoParam->GetFindableClusters()*6.*nTimeBins);
 191
 192         Int_t ntracks = 0;
 193         Int_t nClStack = 0;
 194         for(int istack = 0; istack<AliTRDpropagationLayer::kZones; istack++){
 195                 AliTRDstackLayer stackLayer[kNPlanes*kNTimeBins];
 196
 197                 nClStack = 0;
 198                 //AliInfo(Form("Processing stack %i ...",istack));
 199                 //AliInfo("Building stack propagation layers ...");
 200                 for(int ilayer=0; ilayer<kNPlanes*nTimeBins; ilayer++){
 201                         pp = fGeom->GetPadPlane((Int_t)(ilayer/nTimeBins), istack);
 202                         Double_t stacklength = (pp->GetNrows() - 2) * pp->GetLengthIPad()
 203                                              + 2 * pp->GetLengthOPad() + 2 * pp->GetLengthRim();
 204                         //Debug
 205                         Double_t z0  = fGeom->GetRow0((Int_t)(ilayer/nTimeBins),istack,0);
 206                         const AliTRDpropagationLayer kSMlayer(*(sector->GetLayer(ilayer)));
 207                         stackLayer[ilayer] = kSMlayer;
 208 #ifdef DEBUG
 209                         stackLayer[ilayer].SetDebugStream(fDebugStreamerV1);
 210 #endif
 211                         stackLayer[ilayer].SetRange(z0 - stacklength, stacklength);
 212                         stackLayer[ilayer].SetSector(sector->GetSector());
 213                         stackLayer[ilayer].SetStackNr(istack);
 214                         stackLayer[ilayer].SetNRows(pp->GetNrows());
 215                         stackLayer[ilayer].SetRecoParam(fRecoParam);
 216                         stackLayer[ilayer].BuildIndices();
 217                         nClStack += stackLayer[ilayer].GetNClusters();
 218                 }
 219                 //AliInfo(Form("Finish building stack propagation layers. nClusters %d.", nClStack));
 220                 if(nClStack < kFindable) continue;
 221                 ntracks += Clusters2TracksStack(&stackLayer[0], &esdTrackList);
 222         }
 223         //AliInfo(Form("Found %d tracks in SM", ntracks));
 224
 225         for(int itrack=0; itrack<ntracks; itrack++)
 226           esd->AddTrack((AliESDtrack*)esdTrackList[itrack]);
 227
 228         return ntracks;
 229 }
 230
 231 //____________________________________________________________________
 232 Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDstackLayer *layer
 233                                           , TClonesArray *esdTrackList)
 234 {
 235   //
 236   // Make tracks in one TRD stack.
 237   //
 238   // Parameters :
 239   //   layer  : Array of stack propagation layers containing clusters
 240   //   esdTrackList  : Array of ESD tracks found by the stand alone tracker.
 241   //                   On exit the tracks found in this stack are appended.
 242   //
 243   // Output :
 244   //   Number of tracks found in this stack.
 245   //
 246   // Detailed description
 247   //
 248   // 1. Find the 3 most useful seeding chambers. See BuildSeedingConfigs() for details.
 249   // 2. Steer AliTRDtrackerV1::MakeSeeds() for 3 seeding layer configurations.
 250   //    See AliTRDtrackerV1::MakeSeeds() for more details.
 251   // 3. Arrange track candidates in decreasing order of their quality
 252   // 4. Classify tracks in 5 categories according to:
 253   //    a) number of layers crossed
 254   //    b) track quality
 255   // 5. Sign clusters by tracks in decreasing order of track quality
 256   // 6. Build AliTRDtrack out of seeding tracklets
 257   // 7. Cook MC label
 258   // 8. Build ESD track and register it to the output list
 259   //
 260
 261         AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
 262         Int_t nTimeBins = cal->GetNumberOfTimeBins();
 263         AliTRDseedV1 sseed[kMaxTracksStack*6]; // to be initialized
 264         Int_t pars[3]; // MakeSeeds parameters
 265
 266         //Double_t alpha = AliTRDgeometry::GetAlpha();
 267         //Double_t shift = .5 * alpha;
 268         Int_t configs[kNConfigs];
 269
 270         // Build initial seeding configurations
 271         Double_t quality = BuildSeedingConfigs(layer, configs);
 272 #ifdef DEBUG
 273                 if(AliTRDReconstructor::StreamLevel() > 1)
 274                   AliInfo(Form("Plane config %d %d %d Quality %f"
 275                               , configs[0], configs[1], configs[2], quality));
 276 #endif
 277
 278         // Initialize contors
 279         Int_t ntracks,      // number of TRD track candidates
 280               ntracks1,     // number of registered TRD tracks/iter
 281               ntracks2 = 0; // number of all registered TRD tracks in stack
 282         fSieveSeeding = 0;
 283         do{
 284                 // Loop over seeding configurations
 285                 ntracks = 0; ntracks1 = 0;
 286                 for (Int_t iconf = 0; iconf<3; iconf++) {
 287                         pars[0] = configs[iconf];
 288                         pars[1] = layer->GetStackNr();
 289                         pars[2] = ntracks;
 290                         ntracks = MakeSeeds(layer, &sseed[6*ntracks], pars);
 291                         if(ntracks == kMaxTracksStack) break;
 292                 }
 293 #ifdef DEBUG
 294                 if(AliTRDReconstructor::StreamLevel() > 1)
 295                   AliInfo(Form("Candidate TRD tracks %d in stack %d.", ntracks, pars[1]));
 296 #endif
 297                 if(!ntracks) break;
 298
 299                 // Sort the seeds according to their quality
 300                 Int_t sort[kMaxTracksStack];
 301                 TMath::Sort(ntracks, fTrackQuality, sort, kTRUE);
 302
 303                 // Initialize number of tracks so far and logic switches
 304                 Int_t ntracks0 = esdTrackList->GetEntriesFast();
 305                 Bool_t signedTrack[kMaxTracksStack];
 306                 Bool_t fakeTrack[kMaxTracksStack];
 307                 for (Int_t i=0; i<ntracks; i++){
 308                         signedTrack[i] = kFALSE;
 309                         fakeTrack[i] = kFALSE;
 310                 }
 311                 //AliInfo("Selecting track candidates ...");
 312
 313                 // Sieve clusters in decreasing order of track quality
 314                 Double_t trackParams[7];
 315 //              AliTRDseedV1 *lseed = 0x0;
 316                 Int_t jSieve = 0, candidates;
 317                 do{
 318                         //AliInfo(Form("\t\tITER = %i ", jSieve));
 319
 320                         // Check track candidates
 321                         candidates = 0;
 322                         for (Int_t itrack = 0; itrack < ntracks; itrack++) {
 323                                 Int_t trackIndex = sort[itrack];
 324                                 if (signedTrack[trackIndex] || fakeTrack[trackIndex]) continue;
 325
 326
 327                                 // Calculate track parameters from tracklets seeds
 328                                 Int_t labelsall[1000];
 329                                 Int_t nlabelsall = 0;
 330                                 Int_t naccepted  = 0;
 331                                 Int_t ncl        = 0;
 332                                 Int_t nused      = 0;
 333                                 Int_t nlayers    = 0;
 334                                 Int_t findable   = 0;
 335                                 for (Int_t jLayer = 0; jLayer < kNPlanes; jLayer++) {
 336                                         Int_t jseed = kNPlanes*trackIndex+jLayer;
 337                                         if (TMath::Abs(sseed[jseed].GetYref(0) / sseed[jseed].GetX0()) < 0.15)
 338                                           findable++;
 339
 340                                         if(!sseed[jseed].IsOK()) continue;
 341                                         sseed[jseed].UpdateUsed();
 342                                         ncl   += sseed[jseed].GetN2();
 343                                         nused += sseed[jseed].GetNUsed();
 344                                         nlayers++;
 345
 346                                         // Cooking label
 347                                         for (Int_t itime = 0; itime < nTimeBins; itime++) {
 348                                                 if(!sseed[jseed].IsUsable(itime)) continue;
 349                                                 naccepted++;
 350                                                 Int_t tindex = 0, ilab = 0;
 351                                                 while(ilab<3 && (tindex = sseed[jseed].GetClusters(itime)->GetLabel(ilab)) >= 0){
 352                                                         labelsall[nlabelsall++] = tindex;
 353                                                         ilab++;
 354                                                 }
 355                                         }
 356                                 }
 357                                 // Filter duplicated tracks
 358                                 if (nused > 30){
 359                                         //printf("Skip nused %d\n", nused);
 360                                         fakeTrack[trackIndex] = kTRUE;
 361                                         continue;
 362                                 }
 363                                 if (Float_t(nused)/ncl >= .25){
 364                                         //printf("Skip nused/ncl >= .25\n");
 365                                         fakeTrack[trackIndex] = kTRUE;
 366                                         continue;
 367                                 }
 368
 369                                 // Classify tracks
 370                                 Bool_t skip = kFALSE;
 371                                 switch(jSieve){
 372                                 case 0:
 373                                         if(nlayers < 6) {skip = kTRUE; break;}
 374                                         if(TMath::Log(1.E-9+fTrackQuality[trackIndex]) < -5.){skip = kTRUE; break;}
 375                                         break;
 376
 377                                 case 1:
 378                                         if(nlayers < findable){skip = kTRUE; break;}
 379                                         if(TMath::Log(1.E-9+fTrackQuality[trackIndex]) < -4.){skip = kTRUE; break;}
 380                                         break;
 381
 382                                 case 2:
 383                                         if ((nlayers == findable) || (nlayers == 6)) { skip = kTRUE; break;}
 384                                         if (TMath::Log(1.E-9+fTrackQuality[trackIndex]) < -6.0){skip = kTRUE; break;}
 385                                         break;
 386
 387                                 case 3:
 388                                         if (TMath::Log(1.E-9+fTrackQuality[trackIndex]) < -5.){skip = kTRUE; break;}
 389                                         break;
 390
 391                                 case 4:
 392                                         if (nlayers == 3){skip = kTRUE; break;}
 393                                         if (TMath::Log(1.E-9+fTrackQuality[trackIndex]) - nused/(nlayers-3.0) < -15.0){skip = kTRUE; break;}
 394                                         break;
 395                                 }
 396                                 if(skip){
 397                                         candidates++;
 398                                         //printf("REJECTED : %d [%d] nlayers %d trackQuality = %e nused %d\n", itrack, trackIndex, nlayers, fTrackQuality[trackIndex], nused);
 399                                         continue;
 400                                 }
 401                                 signedTrack[trackIndex] = kTRUE;
 402
 403
 404                                 // Build track label - what happens if measured data ???
 405                                 Int_t labels[1000];
 406                                 Int_t outlab[1000];
 407                                 Int_t nlab = 0;
 408                                 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
 409                                         Int_t jseed = kNPlanes*trackIndex+iLayer;
 410                                         if(!sseed[jseed].IsOK()) continue;
 411                                         for(int ilab=0; ilab<2; ilab++){
 412                                                 if(sseed[jseed].GetLabels(ilab) < 0) continue;
 413                                                 labels[nlab] = sseed[jseed].GetLabels(ilab);
 414                                                 nlab++;
 415                                         }
 416                                 }
 417                                 Freq(nlab,labels,outlab,kFALSE);
 418                                 Int_t   label     = outlab[0];
 419                                 Int_t   frequency = outlab[1];
 420                                 Freq(nlabelsall,labelsall,outlab,kFALSE);
 421                                 Int_t   label1    = outlab[0];
 422                                 Int_t   label2    = outlab[2];
 423                                 Float_t fakeratio = (naccepted - outlab[1]) / Float_t(naccepted);
 424
 425
 426                                 // Sign clusters
 427                                 AliTRDcluster *cl = 0x0; Int_t clusterIndex = -1;
 428                                 for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
 429                                         Int_t jseed = kNPlanes*trackIndex+jLayer;
 430                                         if(!sseed[jseed].IsOK()) continue;
 431                                         if(TMath::Abs(sseed[jseed].GetYfit(1) - sseed[jseed].GetYfit(1)) >= .2) continue; // check this condition with Marian
 432                                         sseed[jseed].UseClusters();
 433                                         if(!cl){
 434                                                 Int_t ic = 0;
 435                                                 while(!(cl = sseed[jseed].GetClusters(ic))) ic++;
 436                                                 clusterIndex =  sseed[jseed].GetIndexes(ic);
 437                                         }
 438                                 }
 439                                 if(!cl) continue;
 440
 441
 442                                 // Build track parameters
 443                                 AliTRDseedV1 *lseed =&sseed[trackIndex*6];
 444                                 Int_t idx = 0;
 445                                 while(idx<3 && !lseed->IsOK()) {
 446                                         idx++;
 447                                         lseed++;
 448                                 }
 449                                 Double_t cR = lseed->GetC();
 450                                 trackParams[1] = lseed->GetYref(0);
 451                                 trackParams[2] = lseed->GetZref(0);
 452                                 trackParams[3] = lseed->GetX0() * cR - TMath::Sin(TMath::ATan(lseed->GetYref(1)));
 453                                 trackParams[4] = lseed->GetZref(1) / TMath::Sqrt(1. + lseed->GetYref(1) * lseed->GetYref(1));
 454                                 trackParams[5] = cR;
 455                                 trackParams[0] = lseed->GetX0();
 456                                 trackParams[6] = layer[0].GetSector();/* *alpha+shift;  // Supermodule*/
 457
 458 #ifdef DEBUG
 459                                 if(AliTRDReconstructor::StreamLevel() > 1) printf("Track %d [%d] nlayers %d trackQuality = %e nused %d, yref = %3.3f\n", itrack, trackIndex, nlayers, fTrackQuality[trackIndex], nused, trackParams[1]);
 460
 461                                 if(AliTRDReconstructor::StreamLevel() >= 1){
 462                                         Int_t sector = layer[0].GetSector();
 463                                         Int_t nclusters = 0;
 464                                         AliTRDseedV1 *dseed[6];
 465                                         for(int is=0; is<6; is++){
 466                                                 dseed[is] = new AliTRDseedV1(sseed[trackIndex*6+is], kTRUE);
 467                                                 nclusters += sseed[is].GetN2();
 468                                                 //for(int ic=0; ic<30; ic++) if(sseed[trackIndex*6+is].GetClusters(ic)) printf("l[%d] tb[%d] cptr[%p]\n", is, ic, sseed[trackIndex*6+is].GetClusters(ic));
 469                                         }
 470                                         //Int_t eventNrInFile = esd->GetEventNumberInFile();
 471                                         //AliInfo(Form("Number of clusters %d.", nclusters));
 472                                         TTreeSRedirector &cstreamer = *fDebugStreamerV1;
 473                                         cstreamer << "Clusters2TracksStack"
 474                                                 << "Iter="      << fSieveSeeding
 475                                                 << "Like="      << fTrackQuality[trackIndex]
 476                                                 << "S0.="       << dseed[0]
 477                                                 << "S1.="       << dseed[1]
 478                                                 << "S2.="       << dseed[2]
 479                                                 << "S3.="       << dseed[3]
 480                                                 << "S4.="       << dseed[4]
 481                                                 << "S5.="       << dseed[5]
 482                                                 << "p0=" << trackParams[0]
 483                                                 << "p1=" << trackParams[1]
 484                                                 << "p2=" << trackParams[2]
 485                                                 << "p3=" << trackParams[3]
 486                                                 << "p4=" << trackParams[4]
 487                                                 << "p5=" << trackParams[5]
 488                                                 << "p6=" << trackParams[6]
 489                                                 << "Sector="    << sector
 490                                                 << "Stack="     << pars[1]
 491                                                 << "Label="     << label
 492                                                 << "Label1="    << label1
 493                                                 << "Label2="    << label2
 494                                                 << "FakeRatio=" << fakeratio
 495                                                 << "Freq="      << frequency
 496                                                 << "Ncl="       << ncl
 497                                                 << "NLayers="   << nlayers
 498                                                 << "Findable="  << findable
 499                                                 << "NUsed="     << nused
 500                                                 << "\n";
 501                                         //???for(int is=0; is<6; is++) delete dseed[is];
 502                                 }
 503 #endif
 504
 505                                 AliTRDtrack *track = AliTRDtrackerV1::RegisterSeed(&sseed[trackIndex*kNPlanes], trackParams);
 506                                 if(!track){
 507                                         AliWarning("Fail to build a TRD Track.");
 508                                         continue;
 509                                 }
 510                                 AliESDtrack esdTrack;
 511                                 esdTrack.UpdateTrackParams(track, AliESDtrack::kTRDout);
 512                                 esdTrack.SetLabel(track->GetLabel());
 513                                 new ((*esdTrackList)[ntracks0++]) AliESDtrack(esdTrack);
 514                                 ntracks1++;
 515                         }
 516
 517                         jSieve++;
 518                 } while(jSieve<5 && candidates); // end track candidates sieve
 519                 if(!ntracks1) break;
 520
 521                 // increment counters
 522                 ntracks2 += ntracks1;
 523                 fSieveSeeding++;
 524
 525                 // Rebuild plane configurations and indices taking only unused clusters into account
 526                 quality = BuildSeedingConfigs(layer, configs);
 527                 //if(quality < fRecoParam->GetPlaneQualityThreshold()) break;
 528
 529                 for(Int_t il = 0; il < kNPlanes * nTimeBins; il++) layer[il].BuildIndices(fSieveSeeding);
 530
 531 #ifdef DEBUG
 532                                 if(AliTRDReconstructor::StreamLevel() > 1) AliInfo(Form("Sieve level %d Plane config %d %d %d Quality %f", fSieveSeeding, configs[0], configs[1], configs[2], quality));
 533 #endif
 534         } while(fSieveSeeding<10); // end stack clusters sieve
 535
 536
 537
 538         //AliInfo(Form("Registered TRD tracks %d in stack %d.", ntracks2, pars[1]));
 539
 540         return ntracks2;
 541 }
 542
 543 //___________________________________________________________________
 544 Double_t AliTRDtrackerV1::BuildSeedingConfigs(AliTRDstackLayer *layers
 545                                             , Int_t *configs)
 546 {
 547   //
 548   // Assign probabilities to chambers according to their
 549   // capability of producing seeds.
 550   //
 551   // Parameters :
 552   //
 553   //   layers : Array of stack propagation layers for all 6 chambers in one stack
 554   //   configs : On exit array of configuration indexes (see GetSeedingConfig()
 555   // for details) in the decreasing order of their seeding probabilities.
 556   //
 557   // Output :
 558   //
 559   //  Return top configuration quality
 560   //
 561   // Detailed description:
 562   //
 563   // To each chamber seeding configuration (see GetSeedingConfig() for
 564   // the list of all configurations) one defines 2 quality factors:
 565   //  - an apriori topological quality (see GetSeedingConfig() for details) and
 566   //  - a data quality based on the uniformity of the distribution of
 567   //    clusters over the x range (time bins population). See CookChamberQA() for details.
 568   // The overall chamber quality is given by the product of this 2 contributions.
 569   //
 570
 571         AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
 572         Int_t nTimeBins = cal->GetNumberOfTimeBins();
 573
 574         Double_t chamberQA[kNPlanes];
 575         for(int iplane=0; iplane<kNPlanes; iplane++){
 576                 chamberQA[iplane] = CookPlaneQA(&layers[iplane*nTimeBins]);
 577                 //printf("chamberQA[%d] = %f\n", iplane, chamberQA[iplane]);
 578         }
 579
 580         Double_t tconfig[kNConfigs];
 581         Int_t planes[4];
 582         for(int iconf=0; iconf<kNConfigs; iconf++){
 583                 GetSeedingConfig(iconf, planes);
 584                 tconfig[iconf] = fgTopologicQA[iconf];
 585                 for(int iplane=0; iplane<4; iplane++) tconfig[iconf] *= chamberQA[planes[iplane]];
 586         }
 587
 588         TMath::Sort(kNConfigs, tconfig, configs, kTRUE);
 589         return tconfig[configs[0]];
 590 }
 591
 592 //____________________________________________________________________
 593 Int_t AliTRDtrackerV1::MakeSeeds(AliTRDstackLayer *layers
 594                                , AliTRDseedV1 *sseed
 595                                , Int_t *ipar)
 596 {
 597   //
 598   // Make tracklet seeds in the TRD stack.
 599   //
 600   // Parameters :
 601   //   layers : Array of stack propagation layers containing clusters
 602   //   sseed  : Array of empty tracklet seeds. On exit they are filled.
 603   //   ipar   : Control parameters:
 604   //       ipar[0] -> seeding chambers configuration
 605   //       ipar[1] -> stack index
 606   //       ipar[2] -> number of track candidates found so far
 607   //
 608   // Output :
 609   //   Number of tracks candidates found.
 610   //
 611   // Detailed description
 612   //
 613   // The following steps are performed:
 614   // 1. Select seeding layers from seeding chambers
 615   // 2. Select seeding clusters from the seeding AliTRDpropagationLayerStack.
 616   //   The clusters are taken from layer 3, layer 0, layer 1 and layer 2, in
 617   //   this order. The parameters controling the range of accepted clusters in
 618   //   layer 0, 1, and 2 are defined in AliTRDstackLayer::BuildCond().
 619   // 3. Helix fit of the cluster set. (see AliTRDtrackerFitter::FitRieman(AliTRDcluster**))
 620   // 4. Initialize seeding tracklets in the seeding chambers.
 621   // 5. Filter 0.
 622   //   Chi2 in the Y direction less than threshold ... (1./(3. - sLayer))
 623   //   Chi2 in the Z direction less than threshold ... (1./(3. - sLayer))
 624   // 6. Attach clusters to seeding tracklets and find linear approximation of
 625   //   the tracklet (see AliTRDseedV1::AttachClustersIter()). The number of used
 626   //   clusters used by current seeds should not exceed ... (25).
 627   // 7. Filter 1.
 628   //   All 4 seeding tracklets should be correctly constructed (see
 629   //   AliTRDseedV1::AttachClustersIter())
 630   // 8. Helix fit of the seeding tracklets
 631   // 9. Filter 2.
 632   //   Likelihood calculation of the fit. (See AliTRDtrackerV1::CookLikelihood() for details)
 633   // 10. Extrapolation of the helix fit to the other 2 chambers:
 634   //    a) Initialization of extrapolation tracklet with fit parameters
 635   //    b) Helix fit of tracklets
 636   //    c) Attach clusters and linear interpolation to extrapolated tracklets
 637   //    d) Helix fit of tracklets
 638   // 11. Improve seeding tracklets quality by reassigning clusters.
 639   //      See AliTRDtrackerV1::ImproveSeedQuality() for details.
 640   // 12. Helix fit of all 6 seeding tracklets and chi2 calculation
 641   // 13. Hyperplane fit and track quality calculation. See AliTRDtrackerFitter::FitHyperplane() for details.
 642   // 14. Cooking labels for tracklets. Should be done only for MC
 643   // 15. Register seeds.
 644   //
 645
 646         AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
 647         Int_t nTimeBins = cal->GetNumberOfTimeBins();
 648         AliTRDcluster *c[4] = {0x0, 0x0, 0x0, 0x0}; // initilize seeding clusters
 649         AliTRDseedV1 *cseed = &sseed[0]; // initialize tracklets for first track
 650         Int_t ncl, mcl; // working variable for looping over clusters
 651         Int_t index[AliTRDstackLayer::kMaxClustersLayer], jndex[AliTRDstackLayer::kMaxClustersLayer];
 652         // chi2 storage
 653         // chi2[0] = tracklet chi2 on the Z direction
 654         // chi2[1] = tracklet chi2 on the R direction
 655         Double_t chi2[4];
 656
 657
 658         // this should be data member of AliTRDtrack
 659         Double_t seedQuality[kMaxTracksStack];
 660
 661         // unpack control parameters
 662         Int_t config  = ipar[0];
 663         Int_t istack  = ipar[1];
 664         Int_t ntracks = ipar[2];
 665         Int_t planes[kNSeedPlanes]; GetSeedingConfig(config, planes);
 666 #ifdef DEBUG
 667                 if(AliTRDReconstructor::StreamLevel() > 1) AliInfo(Form("Making seeds Stack[%d] Config[%d] Tracks[%d]...", istack, config, ntracks));
 668 #endif
 669
 670         // Init chambers geometry
 671         Double_t hL[kNPlanes];       // Tilting angle
 672         Float_t padlength[kNPlanes]; // pad lenghts
 673         AliTRDpadPlane *pp;
 674         for(int il=0; il<kNPlanes; il++){
 675                 pp = fGeom->GetPadPlane(il, istack); // istack has to be imported
 676                 hL[il]        = TMath::Tan(-TMath::DegToRad()*pp->GetTiltingAngle());
 677                 padlength[il] = 10.; //pp->GetLengthIPad();
 678         }
 679
 680         Double_t cond0[4], cond1[4], cond2[4];
 681         // make seeding layers (to be moved in Clusters2TracksStack)
 682         AliTRDstackLayer *layer[] = {0x0, 0x0, 0x0, 0x0};
 683         for(int isl=0; isl<kNSeedPlanes; isl++) layer[isl] = MakeSeedingLayer(&layers[planes[isl] * nTimeBins], planes[isl]);
 684
 685
 686         // Start finding seeds
 687         Int_t icl = 0;
 688         while((c[3] = (*layer[3])[icl++])){
 689                 if(!c[3]) continue;
 690                 layer[0]->BuildCond(c[3], cond0, 0);
 691                 layer[0]->GetClusters(cond0, index, ncl);
 692                 Int_t jcl = 0;
 693                 while(jcl<ncl) {
 694                         c[0] = (*layer[0])[index[jcl++]];
 695                         if(!c[0]) continue;
 696                         Double_t dx    = c[3]->GetX() - c[0]->GetX();
 697                         Double_t theta = (c[3]->GetZ() - c[0]->GetZ())/dx;
 698                         Double_t phi   = (c[3]->GetY() - c[0]->GetY())/dx;
 699                         layer[1]->BuildCond(c[0], cond1, 1, theta, phi);
 700                         layer[1]->GetClusters(cond1, jndex, mcl);
 701
 702                         Int_t kcl = 0;
 703                         while(kcl<mcl) {
 704                                 c[1] = (*layer[1])[jndex[kcl++]];
 705                                 if(!c[1]) continue;
 706                                 layer[2]->BuildCond(c[1], cond2, 2, theta, phi);
 707                                 c[2] = layer[2]->GetNearestCluster(cond2);
 708                                 if(!c[2]) continue;
 709
 710                                 //AliInfo("Seeding clusters found. Building seeds ...");
 711                                 //for(Int_t i = 0; i < kNSeedPlanes; i++) printf("%i. coordinates: x = %3.3f, y = %3.3f, z = %3.3f\n", i, c[i]->GetX(), c[i]->GetY(), c[i]->GetZ());
 712                                 for (Int_t il = 0; il < 6; il++) cseed[il].Reset();
 713
 714                                 fFitter->Reset();
 715
 716                                 fFitter->FitRieman(c, kNSeedPlanes);
 717
 718                                 chi2[0] = 0.; chi2[1] = 0.;
 719                                 AliTRDseedV1 *tseed = 0x0;
 720                                 for(int iLayer=0; iLayer<kNSeedPlanes; iLayer++){
 721                                         tseed = &cseed[planes[iLayer]];
 722                                         tseed->SetRecoParam(fRecoParam);
 723                                         tseed->SetLayer(planes[iLayer]);
 724                                         tseed->SetTilt(hL[planes[iLayer]]);
 725                                         tseed->SetPadLength(TMath::Sqrt(c[iLayer]->GetSigmaZ2()*12));
 726                                         tseed->SetX0(layer[iLayer]->GetX());
 727
 728                                         tseed->Update(fFitter->GetRiemanFitter());
 729                                         chi2[0] += tseed->GetChi2Z(c[iLayer]->GetZ());
 730                                         chi2[1] += tseed->GetChi2Y(c[iLayer]->GetY());
 731                                 }
 732
 733                                 Bool_t isFake = kFALSE;
 734                                 if (c[0]->GetLabel(0) != c[3]->GetLabel(0)) isFake = kTRUE;
 735                                 if (c[1]->GetLabel(0) != c[3]->GetLabel(0)) isFake = kTRUE;
 736                                 if (c[2]->GetLabel(0) != c[3]->GetLabel(0)) isFake = kTRUE;
 737 #ifdef DEBUG
 738                                 if(AliTRDReconstructor::StreamLevel() >= 2){
 739                                         Float_t yref[4], ycluster[4];
 740                                         for(int il=0; il<4; il++){
 741                                                 tseed = &cseed[planes[il]];
 742                                                 yref[il] = tseed->GetYref(0);
 743                                                 ycluster[il] = c[il]->GetY();
 744                                         }
 745                                         Float_t threshold = .5;//1./(3. - sLayer);
 746                                         Int_t ll = c[3]->GetLabel(0);
 747                                         TTreeSRedirector &cs0 = *fDebugStreamerV1;
 748                                                         cs0 << "MakeSeeds0"
 749                                                         <<"isFake=" << isFake
 750                                                         <<"label=" << ll
 751                                                         <<"threshold=" << threshold
 752                                                         <<"chi2=" << chi2[1]
 753                                                         <<"yref0="<<yref[0]
 754                                                         <<"yref1="<<yref[1]
 755                                                         <<"yref2="<<yref[2]
 756                                                         <<"yref3="<<yref[3]
 757                                                         <<"ycluster0="<<ycluster[0]
 758                                                         <<"ycluster1="<<ycluster[1]
 759                                                         <<"ycluster2="<<ycluster[2]
 760                                                         <<"ycluster3="<<ycluster[3]
 761                                                         <<"\n";
 762                                 }
 763 #endif
 764
 765                                 if(chi2[0] > fRecoParam->GetChi2Z()/*7./(3. - sLayer)*//*iter*/){
 766                                         //AliInfo(Form("Failed chi2 filter on chi2Z [%f].", chi2[0]));
 767                                         continue;
 768                                 }
 769                                 if(chi2[1] > fRecoParam->GetChi2Y()/*1./(3. - sLayer)*//*iter*/){
 770                                         //AliInfo(Form("Failed chi2 filter on chi2Y [%f].", chi2[1]));
 771                                         continue;
 772                                 }
 773                                 //AliInfo("Passed chi2 filter.");
 774
 775 #ifdef DEBUG
 776                                 if(AliTRDReconstructor::StreamLevel() >= 2){
 777                                         Float_t minmax[2] = { -100.0,  100.0 };
 778                                         for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
 779                                                 Float_t max = c[iLayer]->GetZ() + cseed[planes[iLayer]].GetPadLength() * 0.5 + 1.0 - cseed[planes[iLayer]].GetZref(0);
 780                                                 if (max < minmax[1]) minmax[1] = max;
 781                                                 Float_t min = c[iLayer]->GetZ()-cseed[planes[iLayer]].GetPadLength() * 0.5 - 1.0 - cseed[planes[iLayer]].GetZref(0);
 782                                                 if (min > minmax[0]) minmax[0] = min;
 783                                         }
 784                                         Double_t xpos[4];
 785                                         for(Int_t l = 0; l < kNSeedPlanes; l++) xpos[l] = layer[l]->GetX();
 786                                         TTreeSRedirector &cstreamer = *fDebugStreamerV1;
 787                                                         cstreamer << "MakeSeeds1"
 788                                                 << "isFake=" << isFake
 789                                                 << "config="   << config
 790                                                 << "Cl0.="   << c[0]
 791                                                 << "Cl1.="   << c[1]
 792                                                 << "Cl2.="   << c[2]
 793                                                 << "Cl3.="   << c[3]
 794                                                 << "X0="     << xpos[0] //layer[sLayer]->GetX()
 795                                                 << "X1="     << xpos[1] //layer[sLayer + 1]->GetX()
 796                                                 << "X2="     << xpos[2] //layer[sLayer + 2]->GetX()
 797                                                 << "X3="     << xpos[3] //layer[sLayer + 3]->GetX()
 798                                                 << "Y2exp="  << cond2[0]
 799                                                 << "Z2exp="  << cond2[1]
 800                                                 << "Chi2R="  << chi2[0]
 801                                                 << "Chi2Z="  << chi2[1]
 802                                                 << "Seed0.=" << &cseed[planes[0]]
 803                                                 << "Seed1.=" << &cseed[planes[1]]
 804                                                 << "Seed2.=" << &cseed[planes[2]]
 805                                                 << "Seed3.=" << &cseed[planes[3]]
 806                                                 << "Zmin="   << minmax[0]
 807                                                 << "Zmax="   << minmax[1]
 808                                                 << "\n" ;
 809                                 }
 810 #endif
 811                                 // try attaching clusters to tracklets
 812                                 Int_t nUsedCl = 0;
 813                                 Int_t nlayers = 0;
 814                                 for(int iLayer=0; iLayer<kNSeedPlanes; iLayer++){
 815                                         AliTRDseedV1 tseed = cseed[planes[iLayer]];
 816                                         if(!tseed.AttachClustersIter(&layers[planes[iLayer]*nTimeBins], 5., kFALSE, c[iLayer])) continue;
 817                                         cseed[planes[iLayer]] = tseed;
 818                                         nUsedCl += cseed[planes[iLayer]].GetNUsed();
 819                                         if(nUsedCl > 25) break;
 820                                         nlayers++;
 821                                 }
 822                                 if(nlayers < kNSeedPlanes){
 823                                         //AliInfo("Failed updating all seeds.");
 824                                         continue;
 825                                 }
 826                                 // fit tracklets and cook likelihood
 827                                 chi2[0] = 0.; chi2[1] = 0.;
 828                                 fFitter->FitRieman(&cseed[0], &planes[0]);
 829                                 AliRieman *rim = fFitter->GetRiemanFitter();
 830                                 for(int iLayer=0; iLayer<4; iLayer++){
 831                                         cseed[planes[iLayer]].Update(rim);
 832                                         chi2[0] += (Float_t)cseed[planes[iLayer]].GetChi2Z();
 833                                         chi2[1] += cseed[planes[iLayer]].GetChi2Y();
 834                                 }
 835                                 Double_t chi2r = chi2[1], chi2z = chi2[0];
 836                                 Double_t like = CookLikelihood(&cseed[0], planes, chi2); // to be checked
 837                                 if (TMath::Log(1.E-9 + like) < fRecoParam->GetTrackLikelihood()){
 838                                         //AliInfo(Form("Failed likelihood %f[%e].", TMath::Log(1.E-9 + like), like));
 839                                         continue;
 840                                 }
 841                                 //AliInfo(Form("Passed likelihood %f[%e].", TMath::Log(1.E-9 + like), like));
 842
 843
 844                                 // book preliminary results
 845                                 seedQuality[ntracks] = like;
 846                                 fSeedLayer[ntracks]  = config;/*sLayer;*/
 847
 848                                 // attach clusters to the extrapolation seeds
 849                                 Int_t lextrap[2];
 850                                 GetExtrapolationConfig(config, lextrap);
 851                                 Int_t nusedf   = 0; // debug value
 852                                 for(int iLayer=0; iLayer<2; iLayer++){
 853                                         Int_t jLayer = lextrap[iLayer];
 854
 855                                         // prepare extrapolated seed
 856                                         cseed[jLayer].Reset();
 857                                         cseed[jLayer].SetRecoParam(fRecoParam);
 858                                         cseed[jLayer].SetLayer(jLayer);
 859                                         cseed[jLayer].SetTilt(hL[jLayer]);
 860                                         cseed[jLayer].SetX0(layers[jLayer * nTimeBins + (nTimeBins/2)].GetX());  // ????????
 861                                         //cseed[jLayer].SetPadLength(??????????);
 862                                         cseed[jLayer].Update(rim);
 863
 864                                         AliTRDcluster *cd = FindSeedingCluster(&layers[jLayer*nTimeBins], &cseed[jLayer]);
 865                                         if(cd == 0x0) continue;
 866 //                                      if(!cd) continue;
 867                                         cseed[jLayer].SetPadLength(TMath::Sqrt(cd->GetSigmaZ2() * 12.));
 868                                         cseed[jLayer].SetX0(cd->GetX());        // reference defined by a seedingLayer which is defined by the x-coordinate of the layers inside
 869
 870                                         // fit extrapolated seed
 871                                         AliTRDseedV1::FitRiemanTilt(cseed, kTRUE);
 872                                         if ((jLayer == 0) && !(cseed[1].IsOK())) continue;
 873                                         if ((jLayer == 5) && !(cseed[4].IsOK())) continue;
 874                                         AliTRDseedV1 tseed = cseed[jLayer];
 875                                         if(!tseed.AttachClustersIter(&layers[jLayer*nTimeBins], 1000.)) continue;
 876                                         cseed[jLayer] = tseed;
 877                                         nusedf += cseed[jLayer].GetNUsed(); // debug value
 878                                         AliTRDseedV1::FitRiemanTilt(cseed, kTRUE);
 879                                 }
 880                                 //AliInfo("Extrapolation done.");
 881
 882                                 ImproveSeedQuality(layers, cseed);
 883                                 //AliInfo("Improve seed quality done.");
 884
 885                                 nlayers   = 0;
 886                                 Int_t nclusters = 0;
 887                                 Int_t findable  = 0;
 888                                 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
 889                                         if (TMath::Abs(cseed[iLayer].GetYref(0) / cseed[iLayer].GetX0()) < 0.15) findable++;
 890                                         if (!cseed[iLayer].IsOK()) continue;
 891                                         nclusters += cseed[iLayer].GetN2();
 892                                         nlayers++;
 893                                 }
 894                                 if (nlayers < 3){
 895                                         //AliInfo("Failed quality check on seeds.");
 896                                         continue;
 897                                 }
 898
 899                                 // fit full track and cook likelihoods
 900                                 fFitter->FitRieman(&cseed[0]);
 901                                 Double_t chi2ZF = 0., chi2RF = 0.;
 902                                 for(int ilayer=0; ilayer<6; ilayer++){
 903                                         cseed[ilayer].Update(fFitter->GetRiemanFitter());
 904                                         if (!cseed[ilayer].IsOK()) continue;
 905                                         //tchi2 = cseed[ilayer].GetChi2Z();
 906                                         //printf("layer %d chi2 %e\n", ilayer, tchi2);
 907                                         chi2ZF += cseed[ilayer].GetChi2Z();
 908                                         chi2RF += cseed[ilayer].GetChi2Y();
 909                                 }
 910                                 chi2ZF /= TMath::Max((nlayers - 3.), 1.);
 911                                 chi2RF /= TMath::Max((nlayers - 3.), 1.);
 912
 913                                 // do the final track fitting
 914                                 fFitter->SetLayers(nlayers);
 915 #ifdef DEBUG
 916                                 fFitter->SetDebugStream(fDebugStreamerV1);
 917 #endif
 918                                 fTrackQuality[ntracks] = fFitter->FitHyperplane(&cseed[0], chi2ZF, GetZ());
 919                                 Double_t param[3];
 920                                 Double_t chi2[2];
 921                                 fFitter->GetHyperplaneFitResults(param);
 922                                 fFitter->GetHyperplaneFitChi2(chi2);
 923                                 //AliInfo("Hyperplane fit done\n");
 924
 925                                 // finalize tracklets
 926                                 Int_t labels[12];
 927                                 Int_t outlab[24];
 928                                 Int_t nlab = 0;
 929                                 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
 930                                         if (!cseed[iLayer].IsOK()) continue;
 931
 932                                         if (cseed[iLayer].GetLabels(0) >= 0) {
 933                                                 labels[nlab] = cseed[iLayer].GetLabels(0);
 934                                                 nlab++;
 935                                         }
 936
 937                                         if (cseed[iLayer].GetLabels(1) >= 0) {
 938                                                 labels[nlab] = cseed[iLayer].GetLabels(1);
 939                                                 nlab++;
 940                                         }
 941                                 }
 942                                 Freq(nlab,labels,outlab,kFALSE);
 943                                 Int_t label     = outlab[0];
 944                                 Int_t frequency = outlab[1];
 945                                 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
 946                                         cseed[iLayer].SetFreq(frequency);
 947                                         cseed[iLayer].SetC(param[1]/*cR*/);
 948                                         cseed[iLayer].SetCC(param[0]/*cC*/);
 949                                         cseed[iLayer].SetChi2(chi2[0]);
 950                                         cseed[iLayer].SetChi2Z(chi2ZF);
 951                                 }
 952
 953 #ifdef DEBUG
 954                                 if(AliTRDReconstructor::StreamLevel() >= 2){
 955                                         Double_t curv = (fFitter->GetRiemanFitter())->GetC();
 956                                         TTreeSRedirector &cstreamer = *fDebugStreamerV1;
 957                                         cstreamer << "MakeSeeds2"
 958                                                 << "C="       << curv
 959                                                 << "Chi2R="   << chi2r
 960                                                 << "Chi2Z="   << chi2z
 961                                                 << "Chi2TR="  << chi2[0]
 962                                                 << "Chi2TC="  << chi2[1]
 963                                                 << "Chi2RF="  << chi2RF
 964                                                 << "Chi2ZF="  << chi2ZF
 965                                                 << "Ncl="     << nclusters
 966                                                 << "Nlayers=" << nlayers
 967                                                 << "NUsedS="  << nUsedCl
 968                                                 << "NUsed="   << nusedf
 969                                                 << "Findable" << findable
 970                                                 << "Like="    << like
 971                                                 << "S0.="     << &cseed[0]
 972                                                 << "S1.="     << &cseed[1]
 973                                                 << "S2.="     << &cseed[2]
 974                                                 << "S3.="     << &cseed[3]
 975                                                 << "S4.="     << &cseed[4]
 976                                                 << "S5.="     << &cseed[5]
 977                                                 << "Label="   << label
 978                                                 << "Freq="    << frequency
 979                                                 << "\n";
 980                                 }
 981 #endif
 982
 983                                 ntracks++;
 984                                 if(ntracks == kMaxTracksStack){
 985                                         AliWarning(Form("Number of seeds reached maximum allowed (%d) in stack.", kMaxTracksStack));
 986                                         return ntracks;
 987                                 }
 988                                 cseed += 6;
 989                         }
 990                 }
 991         }
 992         for(int isl=0; isl<4; isl++) delete layer[isl];
 993
 994         return ntracks;
 995 }
 996
 997 //_____________________________________________________________________________
 998 AliTRDtrack *AliTRDtrackerV1::RegisterSeed(AliTRDseedV1 *seeds, Double_t *params)
 999 {
1000   //
1001   // Build a TRD track out of tracklet candidates
1002   //
1003   // Parameters :
1004   //   seeds  : array of tracklets
1005   //   params : track parameters (see MakeSeeds() function body for a detailed description)
1006   //
1007   // Output :
1008   //   The TRD track.
1009   //
1010   // Detailed description
1011   //
1012   // To be discussed with Marian !!
1013   //
1014
1015         AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1016         Int_t nTimeBins = cal->GetNumberOfTimeBins();
1017
1018   Double_t alpha = AliTRDgeometry::GetAlpha();
1019   Double_t shift = AliTRDgeometry::GetAlpha()/2.0;
1020   Double_t c[15];
1021
1022   c[ 0] = 0.2;
1023   c[ 1] = 0.0; c[ 2] = 2.0;
1024   c[ 3] = 0.0; c[ 4] = 0.0; c[ 5] = 0.02;
1025   c[ 6] = 0.0; c[ 7] = 0.0; c[ 8] = 0.0;  c[ 9] = 0.1;
1026   c[10] = 0.0; c[11] = 0.0; c[12] = 0.0;  c[13] = 0.0; c[14] = params[5]*params[5]*0.01;
1027
1028   Int_t index = 0;
1029   AliTRDcluster *cl = 0;
1030
1031   for (Int_t ilayer = 0; ilayer < 6; ilayer++) {
1032     if (seeds[ilayer].IsOK()) {
1033       for (Int_t itime = nTimeBins - 1; itime > 0; itime--) {
1034         if (seeds[ilayer].GetIndexes(itime) > 0) {
1035           index = seeds[ilayer].GetIndexes(itime);
1036           cl    = seeds[ilayer].GetClusters(itime);
1037           break;
1038         }
1039       }
1040     }
1041     if (index > 0) {
1042       break;
1043     }
1044   }
1045   if (cl == 0) return 0;
1046   AliTRDtrack *track = new AliTRDtrack(cl
1047                                       ,index
1048                                       ,&params[1]
1049                                       ,c
1050                                       ,params[0]
1051                                       ,params[6]*alpha+shift);
1052         // SetCluster(cl, 0); // A. Bercuci
1053         track->PropagateTo(params[0]-5.0);
1054   track->ResetCovariance(1);
1055   Int_t rc = FollowBackProlongation(*track);
1056   if (rc < 30) {
1057     delete track;
1058     track = 0;
1059   }
1060   else {
1061     track->CookdEdx();
1062     track->CookdEdxTimBin(-1);
1063     CookLabel(track,0.9);
1064   }
1065
1066   return track;
1067 }
1068
1069 //____________________________________________________________________
1070 void AliTRDtrackerV1::ImproveSeedQuality(AliTRDstackLayer *layers
1071                                        , AliTRDseedV1 *cseed)
1072 {
1073   //
1074   // Sort tracklets according to "quality" and try to "improve" the first 4 worst
1075   //
1076   // Parameters :
1077   //  layers : Array of propagation layers for a stack/supermodule
1078   //  cseed  : Array of 6 seeding tracklets which has to be improved
1079   //
1080   // Output :
1081   //   cssed : Improved seeds
1082   //
1083   // Detailed description
1084   //
1085   // Iterative procedure in which new clusters are searched for each
1086   // tracklet seed such that the seed quality (see AliTRDseed::GetQuality())
1087   // can be maximized. If some optimization is found the old seeds are replaced.
1088   //
1089
1090         AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1091         Int_t nTimeBins = cal->GetNumberOfTimeBins();
1092
1093         // make a local working copy
1094         AliTRDseedV1 bseed[6];
1095         for (Int_t jLayer = 0; jLayer < 6; jLayer++) bseed[jLayer] = cseed[jLayer];
1096
1097
1098         Float_t lastquality = 10000.0;
1099         Float_t lastchi2    = 10000.0;
1100         Float_t chi2        =  1000.0;
1101
1102         for (Int_t iter = 0; iter < 4; iter++) {
1103                 Float_t sumquality = 0.0;
1104                 Float_t squality[6];
1105                 Int_t   sortindexes[6];
1106
1107                 for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
1108                         squality[jLayer]  = bseed[jLayer].IsOK() ? bseed[jLayer].GetQuality(kTRUE) : -1.;
1109                         sumquality +=squality[jLayer];
1110                 }
1111                 if ((sumquality >= lastquality) || (chi2       >     lastchi2)) break;
1112
1113
1114                 lastquality = sumquality;
1115                 lastchi2    = chi2;
1116                 if (iter > 0) for (Int_t jLayer = 0; jLayer < 6; jLayer++) cseed[jLayer] = bseed[jLayer];
1117
1118
1119                 TMath::Sort(6, squality, sortindexes, kFALSE);
1120                 for (Int_t jLayer = 5; jLayer > 1; jLayer--) {
1121                         Int_t bLayer = sortindexes[jLayer];
1122                         bseed[bLayer].AttachClustersIter(&layers[bLayer*nTimeBins], squality[bLayer], kTRUE);
1123                 }
1124
1125                 chi2 = AliTRDseedV1::FitRiemanTilt(bseed,kTRUE);
1126         } // Loop: iter
1127 }
1128
1129 //____________________________________________________________________
1130 Double_t  AliTRDtrackerV1::CookPlaneQA(AliTRDstackLayer *layers)
1131 {
1132   //
1133   // Calculate plane quality for seeding.
1134   //
1135   //
1136   // Parameters :
1137   //   layers : Array of propagation layers for this plane.
1138   //
1139   // Output :
1140   //   plane quality factor for seeding
1141   //
1142   // Detailed description
1143   //
1144   // The quality of the plane for seeding is higher if:
1145   //  1. the average timebin population is closer to an integer number
1146   //  2. the distribution of clusters/timebin is closer to a uniform distribution.
1147   //    - the slope of the first derivative of a parabolic fit is small or
1148   //    - the slope of a linear fit is small
1149   //
1150
1151         AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1152         Int_t nTimeBins = cal->GetNumberOfTimeBins();
1153
1154 //      Double_t x;
1155 //      TLinearFitter fitter(1, "pol1");
1156 //      fitter.ClearPoints();
1157         Int_t ncl = 0;
1158         Int_t nused = 0;
1159         Int_t nClLayer;
1160         for(int itb=0; itb<nTimeBins; itb++){
1161                 //x = layer[itb].GetX();
1162                 //printf("x[%d] = %f nCls %d\n", itb, x, layer[itb].GetNClusters());
1163                 //if(!layer[itb].GetNClusters()) continue;
1164                 //fitter.AddPoint(&x, layer[itb].GetNClusters(), 1.);
1165                 nClLayer = layers[itb].GetNClusters();
1166                 ncl += nClLayer;
1167                 for(Int_t incl = 0; incl < nClLayer; incl++)
1168                         if((layers[itb].GetCluster(incl))->IsUsed()) nused++;
1169         }
1170
1171         // calculate the deviation of the mean number of clusters from the
1172         // closest integer values
1173         Float_t nclMed = float(ncl-nused)/nTimeBins;
1174         Int_t ncli = Int_t(nclMed);
1175         Float_t nclDev = TMath::Abs(nclMed - TMath::Max(ncli, 1));
1176         nclDev -= (nclDev>.5) && ncli ? .5 : 0.;
1177         /*Double_t quality = */ return TMath::Exp(-2.*nclDev);
1178
1179 //      // get slope of the derivative
1180 //      if(!fitter.Eval()) return quality;
1181 //      fitter.PrintResults(3);
1182 //      Double_t a = fitter.GetParameter(1);
1183 //
1184 //      printf("nclDev(%f)  a(%f)\n", nclDev, a);
1185 //      return quality*TMath::Exp(-a);
1186 }
1187
1188 //____________________________________________________________________
1189 Double_t AliTRDtrackerV1::CookLikelihood(AliTRDseedV1 *cseed
1190                                        , Int_t planes[4]
1191                                        , Double_t *chi2)
1192 {
1193   //
1194   // Calculate the probability of this track candidate.
1195   //
1196   // Parameters :
1197   //   cseeds : array of candidate tracklets
1198   //   planes : array of seeding planes (see seeding configuration)
1199   //   chi2   : chi2 values (on the Z and Y direction) from the rieman fit of the track.
1200   //
1201   // Output :
1202   //   likelihood value
1203   //
1204   // Detailed description
1205   //
1206   // The track quality is estimated based on the following 4 criteria:
1207   //  1. precision of the rieman fit on the Y direction (likea)
1208   //  2. chi2 on the Y direction (likechi2y)
1209   //  3. chi2 on the Z direction (likechi2z)
1210   //  4. number of attached clusters compared to a reference value
1211   //     (see AliTRDrecoParam::fkFindable) (likeN)
1212   //
1213   // The distributions for each type of probabilities are given below as of
1214   // (date). They have to be checked to assure consistency of estimation.
1215   //
1216
1217         AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1218         Int_t nTimeBins = cal->GetNumberOfTimeBins();
1219         // ratio of the total number of clusters/track which are expected to be found by the tracker.
1220         Float_t fgFindable = fRecoParam->GetFindableClusters();
1221
1222
1223         Int_t nclusters = 0;
1224         Double_t sumda = 0.;
1225         for(UChar_t ilayer = 0; ilayer < 4; ilayer++){
1226                 Int_t jlayer = planes[ilayer];
1227                 nclusters += cseed[jlayer].GetN2();
1228                 sumda += TMath::Abs(cseed[jlayer].GetYfitR(1) - cseed[jlayer].GetYref(1));
1229         }
1230         Double_t likea     = TMath::Exp(-sumda*10.6);
1231         Double_t likechi2y  = 0.0000000001;
1232         if (chi2[1] < 0.5) likechi2y += TMath::Exp(-TMath::Sqrt(chi2[1]) * 7.73);
1233         Double_t likechi2z = TMath::Exp(-chi2[0] * 0.088) / TMath::Exp(-chi2[0] * 0.019);
1234         Int_t enc = Int_t(fgFindable*4.*nTimeBins);     // Expected Number Of Clusters, normally 72
1235         Double_t likeN     = TMath::Exp(-(enc - nclusters) * 0.19);
1236
1237         Double_t like      = likea * likechi2y * likechi2z * likeN;
1238
1239 #ifdef DEBUG
1240         //AliInfo(Form("sumda(%f) chi2[0](%f) chi2[1](%f) likea(%f) likechi2y(%f) likechi2z(%f) nclusters(%d) likeN(%f)", sumda, chi2[0], chi2[1], likea, likechi2y, likechi2z, nclusters, likeN));
1241         if(AliTRDReconstructor::StreamLevel() >= 2){
1242                 TTreeSRedirector &cstreamer = *fDebugStreamerV1;
1243                 cstreamer << "CookLikelihood"
1244                         << "sumda="     << sumda
1245                         << "chi0="      << chi2[0]
1246                         << "chi1="      << chi2[1]
1247                         << "likea="     << likea
1248                         << "likechi2y=" << likechi2y
1249                         << "likechi2z=" << likechi2z
1250                         << "nclusters=" << nclusters
1251                         << "likeN="     << likeN
1252                         << "like="      << like
1253                         << "\n";
1254         }
1255 #endif
1256
1257         return like;
1258 }
1259
1260 //___________________________________________________________________
1261 void AliTRDtrackerV1::GetMeanCLStack(AliTRDstackLayer *layers
1262                                    , Int_t *planes
1263                                    , Double_t *params)
1264 {
1265   //
1266   // Determines the Mean number of clusters per layer.
1267   // Needed to determine good Seeding Layers
1268   //
1269   // Parameters:
1270   //    - Array of AliTRDstackLayers
1271   //    - Container for the params
1272   //
1273   // Detailed description
1274   //
1275   // Two Iterations:
1276   // In the first Iteration the mean is calculted using all layers.
1277   // After this, all layers outside the 1-sigma-region are rejected.
1278   // Then the mean value and the standard-deviation are calculted a second
1279   // time in order to select all layers in the 1-sigma-region as good-candidates.
1280   //
1281
1282         AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1283         Int_t nTimeBins = cal->GetNumberOfTimeBins();
1284
1285         Float_t mean = 0, stdev = 0;
1286         Double_t ncl[kNTimeBins*kNSeedPlanes], mcl[kNTimeBins*kNSeedPlanes];
1287         Int_t position = 0;
1288         memset(ncl, 0, sizeof(Int_t)*kNTimeBins*kNSeedPlanes);
1289         memset(mcl, 0, sizeof(Int_t)*kNTimeBins*kNSeedPlanes);
1290         Int_t nused = 0;
1291         for(Int_t ipl = 0; ipl < kNSeedPlanes; ipl++){
1292                 for(Int_t ils = 0; ils < nTimeBins; ils++){
1293                         position = planes[ipl]*nTimeBins + ils;
1294                         ncl[ipl * nTimeBins + ils] = layers[position].GetNClusters();
1295                         nused = 0;
1296                         for(Int_t icl = 0; icl < ncl[ipl * nTimeBins + ils]; icl++)
1297                                 if((layers[position].GetCluster(icl))->IsUsed()) nused++;
1298                         ncl[ipl * nTimeBins + ils] -= nused;
1299                 }
1300         }
1301         // Declaration of quartils:
1302         //Double_t qvals[3] = {0.0, 0.0, 0.0};
1303         //Double_t qprop[3] = {0.16667, 0.5, 0.83333};
1304         // Iterations
1305         Int_t counter;
1306         Double_t *array;
1307         Int_t *limit;
1308         Int_t nLayers = nTimeBins * kNSeedPlanes;
1309         for(Int_t iter = 0; iter < 2; iter++){
1310                 array = (iter == 0) ? &ncl[0] : &mcl[0];
1311                 limit = (iter == 0) ? &nLayers : &counter;
1312                 counter = 0;
1313                 if(iter == 1){
1314                         for(Int_t i = 0; i < nTimeBins *kNSeedPlanes; i++){
1315                                 if((ncl[i] >  mean + stdev) || (ncl[i] <  mean - stdev)) continue; // Outside 1-sigma region
1316 //                              if((ncl[i] >  qvals[2]) || (ncl[i] <  qvals[0])) continue; // Outside 1-sigma region
1317                                 if(ncl[i] == 0) continue;                                                // 0-Layers also rejected
1318                                 mcl[counter] = ncl[i];
1319                                 counter++;
1320                         }
1321                 }
1322                 if(*limit == 0) break;
1323                 printf("Limit = %d\n", *limit);
1324                 //using quartils instead of mean and RMS
1325 //              TMath::Quantiles(*limit,3,array,qvals,qprop,kFALSE);
1326                 mean = TMath::Median(*limit, array, 0x0);
1327                 stdev  = TMath::RMS(*limit, array);
1328         }
1329 //      printf("Quantiles: 0.16667 = %3.3f, 0.5 = %3.3f, 0.83333 = %3.3f\n", qvals[0],qvals[1],qvals[2]);
1330 //      memcpy(params,qvals,sizeof(Double_t)*3);
1331         params[1] = (Double_t)TMath::Nint(mean);
1332         params[0] = (Double_t)TMath::Nint(mean - stdev);
1333         params[2] = (Double_t)TMath::Nint(mean + stdev);
1334
1335 }
1336
1337 //___________________________________________________________________
1338 Int_t AliTRDtrackerV1::GetSeedingLayers(AliTRDstackLayer *layers
1339                                       , Double_t *params)
1340 {
1341   //
1342   // Algorithm to find optimal seeding layer
1343   // Layers inside one sigma region (given by Quantiles) are sorted
1344   // according to their difference.
1345   // All layers outside are sorted according t
1346   //
1347   // Parameters:
1348   //     - Array of AliTRDstackLayers (in the current plane !!!)
1349   //     - Container for the Indices of the seeding Layer candidates
1350   //
1351   // Output:
1352   //     - Number of Layers inside the 1-sigma-region
1353   //
1354   // The optimal seeding layer should contain the mean number of
1355   // custers in the layers in one chamber.
1356   //
1357
1358         //printf("Params: %3.3f, %3.3f, %3.3f\n", params[0], params[1], params[2]);
1359         AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1360         const Int_t kMaxClustersLayer = AliTRDstackLayer::kMaxClustersLayer;
1361         Int_t nTimeBins = cal->GetNumberOfTimeBins();
1362         Int_t ncl[kNTimeBins], indices[kNTimeBins], bins[kMaxClustersLayer];
1363         memset(ncl, 0, sizeof(Int_t)*kNTimeBins);
1364         memset(indices, 0, sizeof(Int_t)*kNTimeBins);
1365         memset(bins, 0, sizeof(Int_t)*kMaxClustersLayer);
1366         Int_t nused = 0;
1367         for(Int_t ils = 0; ils < nTimeBins; ils++){
1368                 ncl[ils] = layers[ils].GetNClusters();
1369                 nused = 0;
1370                 for(Int_t icl = 0; icl < ncl[ils]; icl++)
1371                         if((layers[ils].GetCluster(icl))->IsUsed()) nused++;
1372                 ncl[ils] -= nused;
1373         }
1374
1375         Float_t mean = params[1];
1376         for(Int_t ils = 0; ils < nTimeBins; ils++){
1377                 memmove(indices + bins[ncl[ils]+1] + 1, indices + bins[ncl[ils]+1], sizeof(Int_t)*(nTimeBins - ils));
1378                 indices[bins[ncl[ils]+1]] = ils;
1379                 for(Int_t i = ncl[ils]+1; i < kMaxClustersLayer; i++)
1380                         bins[i]++;
1381         }
1382
1383         //for(Int_t i = 0; i < nTimeBins; i++) printf("Bin %d = %d\n", i, bins[i]);
1384         Int_t sbin = -1;
1385         Int_t nElements;
1386         Int_t position = 0;
1387         TRandom *r = new TRandom();
1388         Int_t iter = 0;
1389         while(1){
1390                 while(sbin < (Int_t)params[0] || sbin > (Int_t)params[2]){
1391                         // Randomly selecting one bin
1392                         sbin = (Int_t)r->Poisson(mean);
1393                 }
1394                 printf("Bin = %d\n",sbin);
1395                 //Randomly selecting one Layer in the bin
1396                 nElements = bins[sbin + 1] - bins[sbin];
1397                 printf("nElements = %d\n", nElements);
1398                 if(iter == 5){
1399                         position = (Int_t)(gRandom->Rndm()*(nTimeBins-1));
1400                         break;
1401                 }
1402                 else if(nElements==0){
1403                         iter++;
1404                         continue;
1405                 }
1406                 position = (Int_t)(gRandom->Rndm()*(nElements-1)) + bins[sbin];
1407                 break;
1408         }
1409         delete r;
1410         return indices[position];
1411 }
1412
1413 //____________________________________________________________________
1414 AliTRDcluster *AliTRDtrackerV1::FindSeedingCluster(AliTRDstackLayer *layers
1415                                                  , AliTRDseedV1/*AliRieman*/ *reference)
1416 {
1417   //
1418   // Finds a seeding Cluster for the extrapolation chamber.
1419   //
1420   // The seeding cluster should be as close as possible to the assumed
1421   // track which is represented by a Rieman fit.
1422   // Therefore the selecting criterion is the minimum distance between
1423   // the best fitting cluster and the Reference which is derived from
1424   // the AliTRDseed. Because all layers are assumed to be equally good
1425   // a linear search is performed.
1426   //
1427   // Imput parameters: - layers: array of AliTRDstackLayers (in one chamber!!!)
1428   //                   - sfit: the reference
1429   //
1430   // Output:           - the best seeding cluster
1431   //
1432
1433         AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1434         Int_t nTimeBins = cal->GetNumberOfTimeBins();
1435
1436         // distances as squared distances
1437         Int_t index = 0;
1438         Float_t ypos = 0.0, zpos = 0.0, distance = 0.0, nearestDistance =100000.0;
1439         ypos = reference->GetYref(0);
1440         zpos = reference->GetZref(0);
1441         AliTRDcluster *currentBest = 0x0, *temp = 0x0;
1442         for(Int_t ils = 0; ils < nTimeBins; ils++){
1443                 // Reference positions
1444 //              ypos = reference->GetYat(layers[ils].GetX());
1445 //              zpos = reference->GetZat(layers[ils].GetX());
1446                 index = layers[ils].SearchNearestCluster(ypos, zpos, fRecoParam->GetRoad2y(), fRecoParam->GetRoad2z());
1447                 if(index == -1) continue;
1448                 temp = layers[ils].GetCluster(index);
1449                 if(!temp) continue;
1450                 distance = (temp->GetY() - ypos) * (temp->GetY() - ypos) + (temp->GetZ() - zpos) * (temp->GetZ() - zpos);
1451                 if(distance < nearestDistance){
1452                         nearestDistance = distance;
1453                         currentBest = temp;
1454                 }
1455         }
1456         return currentBest;
1457 }
1458
1459 //____________________________________________________________________
1460 AliTRDstackLayer *AliTRDtrackerV1::MakeSeedingLayer(AliTRDstackLayer *layers
1461                                                   , Int_t plane)
1462 {
1463   //
1464   // Creates a seeding layer
1465   //
1466
1467         // constants
1468         const Int_t kMaxRows = 16;
1469         const Int_t kMaxCols = 144;
1470         const Int_t kMaxPads = 2304;
1471
1472         // Get the calculation
1473         AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1474         Int_t nTimeBins = cal->GetNumberOfTimeBins();
1475
1476         // Get the geometrical data of the chamber
1477         AliTRDpadPlane *pp = fGeom->GetPadPlane(plane, layers[0].GetStackNr());
1478         Int_t nCols = pp->GetNcols();
1479         Float_t ymin = TMath::Min(pp->GetCol0(), pp->GetColEnd());
1480         Float_t ymax = TMath::Max(pp->GetCol0(), pp->GetColEnd());
1481         Float_t zmin = TMath::Min(pp->GetRow0(), pp->GetRowEnd());
1482         Float_t zmax = TMath::Max(pp->GetRow0(), pp->GetRowEnd());
1483         Int_t nRows = pp->GetNrows();
1484         Float_t binlength = (ymax - ymin)/nCols;
1485         //AliInfo(Form("ymin(%f) ymax(%f) zmin(%f) zmax(%f) nRows(%d) binlength(%f)", ymin, ymax, zmin, zmax, nRows, binlength));
1486
1487         // Fill the histogram
1488         Int_t arrpos;
1489         Float_t ypos;
1490         Int_t irow, nClusters;
1491         Int_t *histogram[kMaxRows];                                                                                     // 2D-Histogram
1492         Int_t hvals[kMaxPads];  memset(hvals, 0, sizeof(Int_t)*kMaxPads);
1493         Float_t *sigmas[kMaxRows];
1494         Float_t svals[kMaxPads];        memset(svals, 0, sizeof(Float_t)*kMaxPads);
1495         AliTRDcluster *c = 0x0;
1496         for(Int_t irs = 0; irs < kMaxRows; irs++){
1497                 histogram[irs] = &hvals[irs*kMaxCols];
1498                 sigmas[irs] = &svals[irs*kMaxCols];
1499         }
1500         for(Int_t iTime = 0; iTime < nTimeBins; iTime++){
1501                 nClusters = layers[iTime].GetNClusters();
1502                 for(Int_t incl = 0; incl < nClusters; incl++){
1503                         c = layers[iTime].GetCluster(incl);
1504                         ypos = c->GetY();
1505                         if(ypos > ymax && ypos < ymin) continue;
1506                         irow = pp->GetPadRowNumber(c->GetZ());                          // Zbin
1507                         if(irow < 0)continue;
1508                         arrpos = static_cast<Int_t>((ypos - ymin)/binlength);
1509                         if(ypos == ymax) arrpos = nCols - 1;
1510                         histogram[irow][arrpos]++;
1511                         sigmas[irow][arrpos] += c->GetSigmaZ2();
1512                 }
1513         }
1514
1515 // Now I have everything in the histogram, do the selection
1516 //      printf("Starting the analysis\n");
1517         //Int_t nPads = nCols * nRows;
1518         // This is what we are interested in: The center of gravity of the best candidates
1519         Float_t cogyvals[kMaxPads]; memset(cogyvals, 0, sizeof(Float_t)*kMaxPads);
1520         Float_t cogzvals[kMaxPads]; memset(cogzvals, 0, sizeof(Float_t)*kMaxPads);
1521         Float_t *cogy[kMaxRows];
1522         Float_t *cogz[kMaxRows];
1523         // Lookup-Table storing coordinates according ti the bins
1524         Float_t yLengths[kMaxCols];
1525         Float_t zLengths[kMaxRows];
1526         for(Int_t icnt = 0; icnt < nCols; icnt++){
1527                 yLengths[icnt] = pp->GetColPos(nCols - 1 - icnt) + binlength/2;
1528         }
1529         for(Int_t icnt = 0; icnt < nRows; icnt++){
1530                 zLengths[icnt] = pp->GetRowPos(icnt) - pp->GetRowSize(icnt)/2;
1531         }
1532
1533         // A bitfield is used to mask the pads as usable
1534         Short_t mask[kMaxCols]; memset(mask, 0 ,sizeof(Short_t) * kMaxCols);//bool mvals[kMaxPads];
1535         for(UChar_t icount = 0; icount < nRows; icount++){
1536                 cogy[icount] = &cogyvals[icount*kMaxCols];
1537                 cogz[icount] = &cogzvals[icount*kMaxCols];
1538         }
1539         // In this array the array position of the best candidates will be stored
1540         Int_t cand[kMaxTracksStack];
1541         Float_t sigcands[kMaxTracksStack];
1542
1543         // helper variables
1544         Int_t indices[kMaxPads]; memset(indices, 0, sizeof(Int_t)*kMaxPads);
1545         Int_t nCandidates = 0;
1546         Float_t norm, cogv;
1547         // histogram filled -> Select best bins
1548         TMath::Sort(kMaxPads, hvals, indices);                  // bins storing a 0 should not matter
1549         // Set Threshold
1550         Int_t maximum = hvals[indices[0]];      // best
1551         Int_t threshold = static_cast<UChar_t>(maximum * fRecoParam->GetFindableClusters());
1552         Int_t col, row, lower, lower1, upper, upper1;
1553         for(Int_t ib = 0; ib < kMaxPads; ib++){
1554                 if(nCandidates >= kMaxTracksStack){
1555                         AliWarning(Form("Number of seed candidates %d exceeded maximum allowed per stack %d", nCandidates, kMaxTracksStack));
1556                         break;
1557                 }
1558                 // Positions
1559                 row = indices[ib]/nCols;
1560                 col = indices[ib]%nCols;
1561                 // here will be the threshold condition:
1562                 if((mask[col] & (1 << row)) != 0) continue;             // Pad is masked: continue
1563                 if(histogram[row][col] < TMath::Max(threshold, 1)){     // of course at least one cluster is needed
1564                         break;                  // number of clusters below threshold: break;
1565                 }
1566                 // passing: Mark the neighbors
1567                 lower  = TMath::Max(col - 1, 0); upper  = TMath::Min(col + 2, nCols);
1568                 lower1 = TMath::Max(row - 1, 0); upper1 = TMath::Min(row + 2, nCols);
1569                 for(Int_t ic = lower; ic < upper; ++ic)
1570                         for(Int_t ir = lower1; ir < upper1; ++ir){
1571                                 if(ic == col && ir == row) continue;
1572                                 mask[ic] |= (1 << ir);
1573                         }
1574                 // Storing the position in an array
1575                 // testing for neigboring
1576                 cogv = 0;
1577                 norm = 0;
1578                 lower = TMath::Max(col - 1,0);
1579                 upper = TMath::Min(col + 2, nCols);
1580                 for(Int_t inb = lower; inb < upper; ++inb){
1581                         cogv += yLengths[inb] * histogram[row][inb];
1582                         norm += histogram[row][inb];
1583                 }
1584                 cogy[row][col] = cogv / norm;
1585                 cogv = 0; norm = 0;
1586                 lower = TMath::Max(row - 1, 0);
1587                 upper = TMath::Min(row + 2, nRows);
1588                 for(Int_t inb = lower; inb < upper; ++inb){
1589                         cogv += zLengths[inb] * histogram[inb][col];
1590                         norm += histogram[inb][col];
1591                 }
1592                 cogz[row][col] = cogv /  norm;
1593                 // passed the filter
1594                 cand[nCandidates] = row*kMaxCols + col; // store the position of a passig candidate into an Array
1595                 sigcands[nCandidates] = sigmas[row][col] / histogram[row][col]; // never be a floating point exeption
1596                 // Analysis output
1597                 nCandidates++;
1598         }
1599         AliTRDstackLayer *fakeLayer = new AliTRDstackLayer(layers[0].GetZ0(), layers[0].GetDZ0(), layers[0].GetStackNr());
1600         fakeLayer->SetX((TMath::Abs(layers[nTimeBins-1].GetX() + layers[0].GetX()))/2);
1601         fakeLayer->SetSector(layers[0].GetSector());
1602         AliTRDcluster **fakeClusters = 0x0;
1603         UInt_t *fakeIndices = 0x0;
1604         if(nCandidates){
1605                 fakeClusters = new AliTRDcluster*[nCandidates];
1606                 fakeIndices = new UInt_t[nCandidates];
1607                 UInt_t fakeIndex = 0;
1608                 for(Int_t ican = 0; ican < nCandidates; ican++){
1609                         fakeClusters[ican] = new AliTRDcluster();
1610                         fakeClusters[ican]->SetX(fakeLayer->GetX());
1611                         fakeClusters[ican]->SetY(cogyvals[cand[ican]]);
1612                         fakeClusters[ican]->SetZ(cogzvals[cand[ican]]);
1613                         fakeClusters[ican]->SetSigmaZ2(sigcands[ican]);
1614                         fakeIndices[ican] = fakeIndex++;// fantasy number
1615                 }
1616         }
1617         fakeLayer->SetRecoParam(fRecoParam);
1618         fakeLayer->SetClustersArray(fakeClusters, nCandidates);
1619         fakeLayer->SetIndexArray(fakeIndices);
1620         fakeLayer->SetNRows(nRows);
1621         fakeLayer->BuildIndices();
1622         //fakeLayer->PrintClusters();
1623
1624 #ifdef DEBUG
1625         if(AliTRDReconstructor::StreamLevel() >= 3){
1626                 TMatrixT<double> hist(nRows, nCols);
1627                 for(Int_t i = 0; i < nRows; i++)
1628                         for(Int_t j = 0; j < nCols; j++)
1629                                 hist(i,j) = histogram[i][j];
1630                 TTreeSRedirector &cstreamer = *fDebugStreamerV1;
1631                 cstreamer << "MakeSeedingLayer"
1632                         << "Iteration="  << fSieveSeeding
1633                         << "plane="      << plane
1634                         << "ymin="       << ymin
1635                         << "ymax="       << ymax
1636                         << "zmin="       << zmin
1637                         << "zmax="       << zmax
1638                         << "L.="         << fakeLayer
1639                         << "Histogram.=" << &hist
1640                         << "\n";
1641         }
1642 #endif
1643         return fakeLayer;
1644 }
1645
1646 //____________________________________________________________________
1647 void AliTRDtrackerV1::GetSeedingConfig(Int_t iconfig
1648                                      , Int_t planes[4]) const
1649 {
1650   //
1651   // Map seeding configurations to detector planes.
1652   //
1653   // Parameters :
1654   //   iconfig : configuration index
1655   //   planes  : member planes of this configuration. On input empty.
1656   //
1657   // Output :
1658   //   planes : contains the planes which are defining the configuration
1659   //
1660   // Detailed description
1661   //
1662   // Here is the list of seeding planes configurations together with
1663   // their topological classification:
1664   //
1665   //  0 - 5432 TQ 0
1666   //  1 - 4321 TQ 0
1667   //  2 - 3210 TQ 0
1668   //  3 - 5321 TQ 1
1669   //  4 - 4210 TQ 1
1670   //  5 - 5431 TQ 1
1671   //  6 - 4320 TQ 1
1672   //  7 - 5430 TQ 2
1673   //  8 - 5210 TQ 2
1674   //  9 - 5421 TQ 3
1675   // 10 - 4310 TQ 3
1676   // 11 - 5410 TQ 4
1677   // 12 - 5420 TQ 5
1678   // 13 - 5320 TQ 5
1679   // 14 - 5310 TQ 5
1680   //
1681   // The topologic quality is modeled as follows:
1682   // 1. The general model is define by the equation:
1683   //  p(conf) = exp(-conf/2)
1684   // 2. According to the topologic classification, configurations from the same
1685   //    class are assigned the agerage value over the model values.
1686   // 3. Quality values are normalized.
1687   //
1688   // The topologic quality distribution as function of configuration is given below:
1689   //Begin_Html
1690   // <img src="gif/topologicQA.gif">
1691   //End_Html
1692   //
1693
1694         switch(iconfig){
1695         case 0: // 5432 TQ 0
1696                 planes[0] = 2;
1697                 planes[1] = 3;
1698                 planes[2] = 4;
1699                 planes[3] = 5;
1700                 break;
1701         case 1: // 4321 TQ 0
1702                 planes[0] = 1;
1703                 planes[1] = 2;
1704                 planes[2] = 3;
1705                 planes[3] = 4;
1706                 break;
1707         case 2: // 3210 TQ 0
1708                 planes[0] = 0;
1709                 planes[1] = 1;
1710                 planes[2] = 2;
1711                 planes[3] = 3;
1712                 break;
1713         case 3: // 5321 TQ 1
1714                 planes[0] = 1;
1715                 planes[1] = 2;
1716                 planes[2] = 3;
1717                 planes[3] = 5;
1718                 break;
1719         case 4: // 4210 TQ 1
1720                 planes[0] = 0;
1721                 planes[1] = 1;
1722                 planes[2] = 2;
1723                 planes[3] = 4;
1724                 break;
1725         case 5: // 5431 TQ 1
1726                 planes[0] = 1;
1727                 planes[1] = 3;
1728                 planes[2] = 4;
1729                 planes[3] = 5;
1730                 break;
1731         case 6: // 4320 TQ 1
1732                 planes[0] = 0;
1733                 planes[1] = 2;
1734                 planes[2] = 3;
1735                 planes[3] = 4;
1736                 break;
1737         case 7: // 5430 TQ 2
1738                 planes[0] = 0;
1739                 planes[1] = 3;
1740                 planes[2] = 4;
1741                 planes[3] = 5;
1742                 break;
1743         case 8: // 5210 TQ 2
1744                 planes[0] = 0;
1745                 planes[1] = 1;
1746                 planes[2] = 2;
1747                 planes[3] = 5;
1748                 break;
1749         case 9: // 5421 TQ 3
1750                 planes[0] = 1;
1751                 planes[1] = 2;
1752                 planes[2] = 4;
1753                 planes[3] = 5;
1754                 break;
1755         case 10: // 4310 TQ 3
1756                 planes[0] = 0;
1757                 planes[1] = 1;
1758                 planes[2] = 3;
1759                 planes[3] = 4;
1760                 break;
1761         case 11: // 5410 TQ 4
1762                 planes[0] = 0;
1763                 planes[1] = 1;
1764                 planes[2] = 4;
1765                 planes[3] = 5;
1766                 break;
1767         case 12: // 5420 TQ 5
1768                 planes[0] = 0;
1769                 planes[1] = 2;
1770                 planes[2] = 4;
1771                 planes[3] = 5;
1772                 break;
1773         case 13: // 5320 TQ 5
1774                 planes[0] = 0;
1775                 planes[1] = 2;
1776                 planes[2] = 3;
1777                 planes[3] = 5;
1778                 break;
1779         case 14: // 5310 TQ 5
1780                 planes[0] = 0;
1781                 planes[1] = 1;
1782                 planes[2] = 3;
1783                 planes[3] = 5;
1784                 break;
1785         }
1786 }
1787
1788 //____________________________________________________________________
1789 void AliTRDtrackerV1::GetExtrapolationConfig(Int_t iconfig
1790                                            , Int_t planes[2]) const
1791 {
1792   //
1793   // Returns the extrapolation planes for a seeding configuration.
1794   //
1795   // Parameters :
1796   //   iconfig : configuration index
1797   //   planes  : planes which are not in this configuration. On input empty.
1798   //
1799   // Output :
1800   //   planes : contains the planes which are not in the configuration
1801   //
1802   // Detailed description
1803   //
1804
1805         switch(iconfig){
1806         case 0: // 5432 TQ 0
1807                 planes[0] = 1;
1808                 planes[1] = 0;
1809                 break;
1810         case 1: // 4321 TQ 0
1811                 planes[0] = 5;
1812                 planes[1] = 0;
1813                 break;
1814         case 2: // 3210 TQ 0
1815                 planes[0] = 4;
1816                 planes[1] = 5;
1817                 break;
1818         case 3: // 5321 TQ 1
1819                 planes[0] = 4;
1820                 planes[1] = 0;
1821                 break;
1822         case 4: // 4210 TQ 1
1823                 planes[0] = 5;
1824                 planes[1] = 3;
1825                 break;
1826         case 5: // 5431 TQ 1
1827                 planes[0] = 2;
1828                 planes[1] = 0;
1829                 break;
1830         case 6: // 4320 TQ 1
1831                 planes[0] = 5;
1832                 planes[1] = 1;
1833                 break;
1834         case 7: // 5430 TQ 2
1835                 planes[0] = 2;
1836                 planes[1] = 1;
1837                 break;
1838         case 8: // 5210 TQ 2
1839                 planes[0] = 4;
1840                 planes[1] = 3;
1841                 break;
1842         case 9: // 5421 TQ 3
1843                 planes[0] = 3;
1844                 planes[1] = 0;
1845                 break;
1846         case 10: // 4310 TQ 3
1847                 planes[0] = 5;
1848                 planes[1] = 2;
1849                 break;
1850         case 11: // 5410 TQ 4
1851                 planes[0] = 3;
1852                 planes[1] = 2;
1853                 break;
1854         case 12: // 5420 TQ 5
1855                 planes[0] = 3;
1856                 planes[1] = 1;
1857                 break;
1858         case 13: // 5320 TQ 5
1859                 planes[0] = 4;
1860                 planes[1] = 1;
1861                 break;
1862         case 14: // 5310 TQ 5
1863                 planes[0] = 4;
1864                 planes[1] = 2;
1865                 break;
1866         }
1867 }