EMCAL/AliEMCALTracker.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 //                       Class AliEMCALTracker
  17 //                      -----------------------
  18 // Implementation of the track matching method between barrel tracks and
  19 // EMCAL clusters.
  20 // Besides algorithm implementation, some cuts are required to be set
  21 // in order to define, for each track, an acceptance window where clusters
  22 // are searched to find best match (if any).
  23 // The class accepts as input an ESD container, and works directly on it,
  24 // simply setting, for each of its tracks, the fEMCALindex flag, for each
  25 // track which is matched to a cluster.
  26 // In order to use method, one must launch PropagateBack().
  27 //
  28 // ------------------------------------------------------------------------
  29 // author: A. Pulvirenti (alberto.pulvirenti@ct.infn.it)
  30 //=========================================================================
  31
  32 #include <Riostream.h>
  33 #include <iomanip>
  34
  35 #include <TFile.h>
  36 #include <TTree.h>
  37 #include <TList.h>
  38 #include <TString.h>
  39 #include <TVector3.h>
  40 #include <TClonesArray.h>
  41 #include <TGeoMatrix.h>
  42
  43 #include "AliLog.h"
  44 #include "AliESDEvent.h"
  45 #include "AliESDtrack.h"
  46 #include "AliESDCaloCluster.h"
  47 #include "AliEMCALRecPoint.h"
  48 #include "AliRunLoader.h"
  49 #include "AliEMCALTrack.h"
  50 #include "AliEMCALLoader.h"
  51 #include "AliEMCALGeometry.h"
  52 #include "AliEMCALReconstructor.h"
  53 #include "AliEMCALRecParam.h"
  54 #include "AliCDBEntry.h"
  55 #include "AliCDBManager.h"
  56
  57 #include "AliEMCALTracker.h"
  58
  59 ClassImp(AliEMCALTracker)
  60
  61 //
  62 //------------------------------------------------------------------------------
  63 //
  64 AliEMCALTracker::AliEMCALTracker()
  65   : AliTracker(),
  66     fNPropSteps(0),
  67     fTrackCorrMode(kTrackCorrNone),
  68     fCutX(50.0),
  69     fCutY(50.0),
  70     fCutZ(50.0),
  71     fCutAlphaMin(-200.0),
  72     fCutAlphaMax(200.0),
  73     fCutAngle(100.0),
  74     fMaxDist(10.0),
  75     fRho(1.0),
  76     fX0(1.0),
  77     fTracks(0),
  78     fClusters(0),
  79     fMatches(0),
  80     fGeom(0)
  81 {
  82         //
  83         // Default constructor.
  84         // Initializes al simple data members to default values,
  85         // and all collections to NULL.
  86         // Output file name is set to a default value.
  87         //
  88
  89         InitParameters();
  90 }
  91 //
  92 //------------------------------------------------------------------------------
  93 //
  94 AliEMCALTracker::AliEMCALTracker(const AliEMCALTracker& copy)
  95   : AliTracker(),
  96     fNPropSteps(copy.fNPropSteps),
  97     fTrackCorrMode(copy.fTrackCorrMode),
  98     fCutX(copy.fCutX),
  99     fCutY(copy.fCutY),
 100     fCutZ(copy.fCutZ),
 101     fCutAlphaMin(copy.fCutAlphaMin),
 102     fCutAlphaMax(copy.fCutAlphaMax),
 103     fCutAngle(copy.fCutAngle),
 104     fMaxDist(copy.fMaxDist),
 105     fRho(copy.fRho),
 106     fX0(copy.fX0),
 107     fTracks((TObjArray*)copy.fTracks->Clone()),
 108     fClusters((TObjArray*)copy.fClusters->Clone()),
 109     fMatches((TList*)copy.fMatches->Clone()),
 110     fGeom(copy.fGeom)
 111 {
 112         //
 113         // Copy constructor
 114         // Besides copying all parameters, duplicates all collections.
 115         //
 116 }
 117 //
 118 //------------------------------------------------------------------------------
 119 //
 120 AliEMCALTracker& AliEMCALTracker::operator=(const AliEMCALTracker& copy)
 121 {
 122         //
 123         // Assignment operator.
 124         // Besides copying all parameters, duplicates all collections.
 125         //
 126
 127         fCutX = copy.fCutX;
 128         fCutY = copy.fCutY;
 129         fCutZ = copy.fCutZ;
 130         fCutAlphaMin = copy.fCutAlphaMin;
 131         fCutAlphaMax = copy.fCutAlphaMax;
 132         fCutAngle = copy.fCutAngle;
 133         fMaxDist = copy.fMaxDist;
 134
 135         fTracks = (TObjArray*)copy.fTracks->Clone();
 136         fClusters = (TObjArray*)copy.fClusters->Clone();
 137         fMatches = (TList*)copy.fMatches->Clone();
 138
 139         fGeom = copy.fGeom;
 140
 141         return (*this);
 142 }
 143 //
 144 //------------------------------------------------------------------------------
 145 //
 146 void AliEMCALTracker::InitParameters()
 147 {
 148         //
 149         // Retrieve initialization parameters
 150         //
 151
 152   // Check if the instance of AliEMCALRecParam exists,
 153   const AliEMCALRecParam* recParam = AliEMCALReconstructor::GetRecParam();
 154   if(!recParam){
 155     AliFatal("Reconstruction parameters for EMCAL not set!");
 156   }
 157
 158   fCutX =  recParam->GetTrkCutX();
 159   fCutY =  recParam->GetTrkCutY();
 160   fCutZ =  recParam->GetTrkCutZ();
 161   fMaxDist =  recParam->GetTrkCutR();
 162   fCutAngle =  recParam->GetTrkCutAngle();
 163   fCutAlphaMin =  recParam->GetTrkCutAlphaMin();
 164   fCutAlphaMax =  recParam->GetTrkCutAlphaMax();
 165
 166 }
 167 //
 168 //------------------------------------------------------------------------------
 169 //
 170 TTree* AliEMCALTracker::SearchTrueMatches()
 171 {
 172   //Search through the list of
 173   //track match candidates and clusters
 174   //and look for true matches
 175   //
 176   //
 177         if (!fClusters) return 0;
 178         if (fClusters->IsEmpty()) return 0;
 179         if (!fTracks) return 0;
 180         if (fTracks->IsEmpty()) return 0;
 181
 182         TTree *outTree = new TTree("tree", "True matches from event");
 183         Int_t indexT, indexC, label;
 184         outTree->Branch("indexC", &indexC, "indexC/I");
 185         outTree->Branch("indexT", &indexT, "indexT/I");
 186         outTree->Branch("label",  &label , "label/I");
 187
 188         Double_t dist;
 189         Int_t ic, nClusters = (Int_t)fClusters->GetEntries();
 190         Int_t it, nTracks = fTracks->GetEntries();
 191
 192         for (ic = 0; ic < nClusters; ic++) {
 193                 AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(ic);
 194                 label = cluster->Label();
 195                 indexC = cluster->Index();
 196                 for (it = 0; it < nTracks; it++) {
 197                         AliEMCALTrack *track = (AliEMCALTrack*)fTracks->At(it);
 198                         if (TMath::Abs(track->GetSeedLabel()) != label) continue;
 199                         dist = CheckPair(track, cluster);
 200                         if (dist <= fMaxDist) {
 201                                 indexT = track->GetSeedIndex();
 202                                 outTree->Fill();
 203                         }
 204                 }
 205         }
 206
 207         return outTree;
 208 }
 209 //
 210 //------------------------------------------------------------------------------
 211 //
 212 void AliEMCALTracker::Clear(Option_t* option)
 213 {
 214         //
 215         // Clearing method
 216         // Deletes all objects in arrays and the arrays themselves
 217         //
 218
 219         TString opt(option);
 220         Bool_t clearTracks = opt.Contains("TRACKS");
 221         Bool_t clearClusters = opt.Contains("CLUSTERS");
 222         Bool_t clearMatches = opt.Contains("MATCHES");
 223         if (opt.Contains("ALL")) {
 224                 clearTracks = kTRUE;
 225                 clearClusters = kTRUE;
 226                 clearMatches = kTRUE;
 227         }
 228
 229         if (fTracks != 0x0 && clearTracks) {
 230            fTracks->Delete();
 231            delete fTracks;
 232            fTracks = 0;
 233         }
 234         if (fClusters != 0x0 && clearClusters) {
 235            fClusters->Delete();
 236            delete fClusters;
 237            fClusters = 0;
 238         }
 239         if (fMatches != 0x0 && clearMatches) {
 240            fMatches->Delete();
 241            delete fMatches;
 242            fMatches = 0;
 243         }
 244 }
 245 //
 246 //------------------------------------------------------------------------------
 247 //
 248 Int_t AliEMCALTracker::LoadClusters(TTree *cTree)
 249 {
 250         //
 251         // Load EMCAL clusters in the form of AliEMCALRecPoint,
 252         // from simulation temporary files.
 253         // (When included in reconstruction chain, this method is used automatically)
 254         //
 255
 256         Clear("CLUSTERS");
 257
 258         TBranch *branch = cTree->GetBranch("EMCALECARP");
 259         if (!branch) {
 260                 AliError("Can't get the branch with the EMCAL clusters");
 261                 return 1;
 262         }
 263
 264         TClonesArray *clusters = new TClonesArray("AliEMCALRecPoint", 1000);
 265         branch->SetAddress(&clusters);
 266
 267         cTree->GetEvent(0);
 268         Int_t nClusters = (Int_t)clusters->GetEntries();
 269         fClusters = new TObjArray(0);
 270         for (Int_t i = 0; i < nClusters; i++) {
 271                 AliEMCALRecPoint *cluster = (AliEMCALRecPoint*)clusters->At(i);
 272                 if (!cluster) continue;
 273                 if (cluster->GetClusterType() != AliESDCaloCluster::kEMCALClusterv1) continue;
 274                 AliEMCALMatchCluster *matchCluster = new AliEMCALMatchCluster(i, cluster);
 275                 fClusters->AddLast(matchCluster);
 276         }
 277
 278         branch->SetAddress(0);
 279         clusters->Delete();
 280         delete clusters;
 281         if (fClusters->IsEmpty())
 282            AliWarning("No clusters collected");
 283
 284         AliInfo(Form("Collected %d clusters (RecPoints)", fClusters->GetEntries()));
 285
 286         return 0;
 287 }
 288 //
 289 //------------------------------------------------------------------------------
 290 //
 291 Int_t AliEMCALTracker::LoadClusters(AliESDEvent *esd)
 292 {
 293         //
 294         // Load EMCAL clusters in the form of AliESDCaloClusters,
 295         // from an AliESD object.
 296         //
 297
 298         // make sure that tracks/clusters collections are empty
 299         Clear("CLUSTERS");
 300
 301         Int_t start = esd->GetFirstEMCALCluster();
 302         Int_t nClustersEMC = esd->GetNumberOfEMCALClusters();
 303         Int_t end = start + nClustersEMC;
 304
 305         fClusters = new TObjArray(0);
 306
 307         Int_t i;
 308         for (i = start; i < end; i++) {
 309                 AliESDCaloCluster *cluster = esd->GetCaloCluster(i);
 310                 if (!cluster) continue;
 311                 if (cluster->GetClusterType() != AliESDCaloCluster::kEMCALClusterv1) continue;
 312                 AliEMCALMatchCluster *matchCluster = new AliEMCALMatchCluster(i, cluster);
 313                 fClusters->AddLast(matchCluster);
 314         }
 315         if (fClusters->IsEmpty())
 316            AliWarning("No clusters collected");
 317
 318         AliInfo(Form("Collected %d clusters from ESD", fClusters->GetEntries()));
 319
 320         return 0;
 321 }
 322 //
 323 //------------------------------------------------------------------------------
 324 //
 325 Int_t AliEMCALTracker::LoadTracks(AliESDEvent *esd)
 326 {
 327         //
 328         // Load ESD tracks.
 329         //
 330
 331         Clear("TRACKS");
 332
 333         Int_t nTracks = esd->GetNumberOfTracks();
 334         fTracks = new TObjArray(0);
 335
 336         Int_t i, j;
 337         Bool_t isKink;
 338         Double_t alpha;
 339         for (i = 0; i < nTracks; i++) {
 340                 AliESDtrack *esdTrack = esd->GetTrack(i);
 341                 // set by default the value corresponding to "no match"
 342                 esdTrack->SetEMCALcluster(kUnmatched);
 343 //              if (esdTrack->GetLabel() < 0) continue;
 344 //              if (!(esdTrack->GetStatus() & AliESDtrack::kTOFout)) continue;
 345                 isKink = kFALSE;
 346                 for (j = 0; j < 3; j++) {
 347                         if (esdTrack->GetKinkIndex(j) != 0) isKink = kTRUE;
 348                 }
 349                 if (isKink) continue;
 350                 AliEMCALTrack *track = new AliEMCALTrack(*esdTrack);
 351                 track->SetMass(0.13957018);
 352                 // check alpha and reject the tracks which fall outside EMCAL acceptance
 353                 alpha = track->GetAlpha() * TMath::RadToDeg();
 354                 if (alpha >  -155.0 && alpha < 67.0) {
 355                         delete track;
 356                         continue;
 357                 }
 358 //              if (!PropagateToEMCAL(track)) {
 359 //                      delete track;
 360 //                      continue;
 361 //              }
 362                 track->SetSeedIndex(i);
 363                 track->SetSeedLabel(esdTrack->GetLabel());
 364                 fTracks->AddLast(track);
 365         }
 366         if (fTracks->IsEmpty()) {
 367                 AliWarning("No tracks collected");
 368         }
 369
 370         AliInfo(Form("Collected %d tracks", fTracks->GetEntries()));
 371
 372         return 0;
 373 }
 374 //
 375 //------------------------------------------------------------------------------
 376 //
 377 Int_t AliEMCALTracker::PropagateBack(AliESDEvent* esd)
 378 {
 379         //
 380         // Main operation method.
 381         // Gets external AliESD containing tracks to be matched.
 382         // After executing match finding, stores in the same ESD object all infos
 383         // and releases the object for further reconstruction steps.
 384         //
 385         //
 386         // Note: should always return 0=OK, because otherwise all tracking
 387         // is aborted for this event
 388
 389         if (!esd) {
 390                 AliError("NULL ESD passed");
 391                 return 1;
 392         }
 393
 394         // step 1:
 395         // if cluster array is empty, cluster are collected
 396         // from the passed ESD, and work is done with ESDCaloClusters
 397         Int_t okLoadClusters, nClusters;
 398         if (!fClusters || (fClusters && fClusters->IsEmpty())) {
 399                 okLoadClusters = LoadClusters(esd);
 400         }
 401         nClusters = fClusters->GetEntries();
 402
 403         // step 2:
 404         // collect ESD tracks
 405         Int_t okLoadTracks = LoadTracks(esd), nTracks;
 406         if (okLoadTracks) return 3;
 407         nTracks = fTracks->GetEntries();
 408
 409         // step 3:
 410         // each track is propagated to the "R" position of each cluster.
 411         // The closest cluster is assigned as match.
 412         // IF no clusters lie within the maximum allowed distance, no matches are assigned.
 413         Int_t nMatches = CreateMatches();
 414         if (!nMatches) {
 415                 AliInfo(Form("#clusters = %d -- #tracks = %d --> No good matches found.", nClusters, nTracks));
 416                 return 0;
 417         }
 418         else {
 419                 AliInfo(Form("#clusters = %d -- #tracks = %d --> Found %d matches.", nClusters, nTracks, nMatches));
 420         }
 421
 422         // step 4:
 423         // when more than 1 track share the same matched cluster, only the closest one is kept.
 424         Int_t nRemoved = SolveCompetitions();
 425         AliInfo(Form("Removed %d duplicate matches", nRemoved));
 426         if (nRemoved >= nMatches) {
 427                 AliError("Removed ALL matches! Check the algorithm or data. Nothing to save");
 428                 return 5;
 429         }
 430
 431         // step 5:
 432         // save obtained information setting the 'fEMCALindex' field of AliESDtrack object
 433         Int_t nSaved = 0, trackID, nGood = 0, nFake = 0;
 434         TListIter iter(fMatches);
 435         AliEMCALMatch *match = 0;
 436         while ( (match = (AliEMCALMatch*)iter.Next()) ) {
 437                 if (!match->CanBeSaved()) continue;
 438                 AliEMCALTrack *track = (AliEMCALTrack*)fTracks->At(match->GetIndexT());
 439                 AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(match->GetIndexC());
 440                 trackID = track->GetSeedIndex();
 441                 AliESDtrack *esdTrack = esd->GetTrack(trackID);
 442                 if (!esdTrack) continue;
 443                 if (TMath::Abs(esdTrack->GetLabel()) == cluster->Label()) {
 444                         esdTrack->SetEMCALcluster(cluster->Index());
 445                         nGood++;
 446                 }
 447                 else {
 448                         esdTrack->SetEMCALcluster(-cluster->Index());
 449                         nFake++;
 450                 }
 451                 nSaved++;
 452         }
 453         /*
 454         AliEMCALTrack *track = 0;
 455         TObjArrayIter tracks(fTracks);
 456         while ( (track = (AliEMCALTrack*)tracks.Next()) ) {
 457                 trackID = track->GetSeedIndex();
 458                 clusterID = track->GetMatchedClusterIndex();
 459                 AliESDtrack *esdTrack = esd->GetTrack(trackID);
 460                 if (!esdTrack) continue;
 461                 if (clusterID < 0) {
 462                         esdTrack->SetEMCALcluster(kUnmatched);
 463                 }
 464                 else {
 465                         AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(clusterID);
 466                         if (!cluster) continue;
 467                         if (esdTrack->GetLabel() == cluster->Label()) {
 468                                 nGood++;
 469                                 esdTrack->SetEMCALcluster(cluster->Index());
 470                         }
 471                         else {
 472                                 esdTrack->SetEMCALcluster(-cluster->Index());
 473                         }
 474                         nSaved++;
 475                 }
 476         }
 477         */
 478         AliInfo(Form("Saved %d matches [%d good + %d fake]", nSaved, nGood, nFake));
 479
 480         return 0;
 481 }
 482 //
 483 //------------------------------------------------------------------------------
 484 //
 485 void AliEMCALTracker::SetTrackCorrectionMode(Option_t *option)
 486 {
 487         //
 488         // Set track correction mode
 489         // gest the choice in string format and converts into
 490         // internal enum
 491         //
 492
 493         TString opt(option);
 494         opt.ToUpper();
 495
 496         if (!opt.CompareTo("NONE")) {
 497                 fTrackCorrMode = kTrackCorrNone;
 498         }
 499         else if (!opt.CompareTo("MMB")) {
 500                 fTrackCorrMode = kTrackCorrMMB;
 501         }
 502         else if (!opt.CompareTo("FIXED")) {
 503                 fTrackCorrMode = kTrackCorrFixed;
 504         }
 505         else {
 506                 cerr << "E-AliEMCALTracker::SetTrackCorrectionMode '" << option << "': Unrecognized option" << endl;
 507         }
 508 }
 509
 510 //
 511 //------------------------------------------------------------------------------
 512 //
 513 Double_t AliEMCALTracker::AngleDiff(Double_t angle1, Double_t angle2)
 514 {
 515         //
 516         // [PRIVATE]
 517         // Given two angles in radiants, it converts them in the range 0-2pi
 518         // then computes their true difference, i.e. if the difference a1-a2
 519         // results to be larger than 180 degrees, it returns 360 - diff.
 520         //
 521
 522         if (angle1 < 0.0) angle1 += TMath::TwoPi();
 523         if (angle1 > TMath::TwoPi()) angle1 -= TMath::TwoPi();
 524         if (angle2 < 0.0) angle2 += TMath::TwoPi();
 525         if (angle2 > TMath::TwoPi()) angle2 -= TMath::TwoPi();
 526
 527         Double_t diff = TMath::Abs(angle1 - angle2);
 528         if (diff > TMath::Pi()) diff = TMath::TwoPi() - diff;
 529
 530         if (angle2 > angle1) diff = -diff;
 531
 532         return diff;
 533 }
 534 //
 535 //------------------------------------------------------------------------------
 536 //
 537 Double_t AliEMCALTracker::CheckPair
 538 (AliEMCALTrack *track, AliEMCALMatchCluster *cl)
 539 {
 540         //
 541         // Given a track and a cluster,
 542         // propagates the first to the radius of the second.
 543         // Then, checks the propagation point against all cuts.
 544         // If at least a cut is not passed, a valuer equal to
 545         // twice the maximum allowed distance is passed (so the value returned
 546         // will not be taken into account when creating matches)
 547         //
 548
 549         // TEMP
 550         Bool_t isTrue = kFALSE;
 551 //      if (tr->GetSeedLabel() == cl->Label()) {
 552 //              isTrue = kTRUE;
 553 //      }
 554
 555         // copy track into temporary variable
 556         AliEMCALTrack *tr = new AliEMCALTrack(*track);
 557
 558         Double_t distance = 2.0 * fMaxDist;
 559
 560         // check against cut on difference 'alpha - phi'
 561         Double_t phi = TMath::ATan2(cl->Y(), cl->X());
 562         phi = AngleDiff(phi, tr->GetAlpha());
 563         if (phi < fCutAlphaMin || phi > fCutAlphaMax){
 564           delete tr;
 565           return distance;
 566         }
 567
 568         // try to propagate to cluster radius
 569         // (return the 'distance' value if it fails)
 570         Double_t pos[3], &x = pos[0], &y = pos[1], &z = pos[2];
 571         Double_t x0, rho;
 572         tr->GetXYZ(pos);
 573         Double_t rt = TMath::Sqrt(x*x + y*y);
 574         Double_t rc = TMath::Sqrt(cl->X()*cl->X() + cl->Y()*cl->Y());
 575
 576         if (fTrackCorrMode == kTrackCorrMMB) {
 577                 Double_t pos1[3], pos2[3], param[6];
 578                 pos1[0] = x;
 579                 pos1[1] = y;
 580                 pos1[2] = z;
 581                 pos2[0] = cl->X();
 582                 pos2[1] = cl->Y();
 583                 pos2[2] = cl->Z();
 584                 MeanMaterialBudget(pos1, pos2, param);
 585                 rho = param[0]*param[4];
 586                 x0 = param[1];
 587         }
 588         else if (fTrackCorrMode == kTrackCorrFixed) {
 589                 rho = fRho;
 590                 x0 = fX0;
 591         }
 592         else {
 593                 rho = 0.0;
 594                 x0 = 0.0;
 595         }
 596         if (fNPropSteps) {
 597                 Int_t i;
 598                 Double_t r;
 599                 cout.setf(ios::fixed);
 600                 cout.precision(5);
 601                 if (isTrue) cout << "Init : " << rt << ' ' << x << ' ' << y << ' ' << z << endl;
 602                 for (i = 0; i < fNPropSteps; i++) {
 603                         r = rt + (rc - rt) * ((Double_t)(i+1)/(Double_t)fNPropSteps);
 604                         if (!tr->PropagateTo(r, x0, rho)){
 605                           delete tr;
 606                           return distance;
 607                         }
 608                         tr->GetXYZ(pos);
 609                         if (isTrue) cout << "Step : " << r << ' ' << x << ' ' << y << ' ' << z << endl;
 610                 }
 611                 if (isTrue) cout << "Clstr: " << rc << ' ' << cl->X() << ' ' << cl->Y() << ' ' << cl->Z() << endl;
 612         }
 613         else {
 614                 // when no steps are used, no correction makes sense
 615                 //if (!tr->PropagateTo(rc, 0.0, 0.0)) return distance;
 616                 if (!tr->PropagateToGlobal(cl->X(), cl->Y(), cl->Z(), 0.0, 0.0)){
 617                   delete tr;
 618                   return distance;
 619                 }
 620                 /*
 621                 Bool_t propOK = kFALSE;
 622                 cout << "START" << endl;
 623                 Double_t dist, rCHK, bestDist = 10000000.0;
 624                 for (Double_t rTMP = rc; rTMP> rc*0.95; rTMP -= 0.1) {
 625                         if (!tr->PropagateTo(rTMP)) continue;
 626                         propOK = kTRUE;
 627                         tr->GetXYZ(pos);
 628                         rCHK = TMath::Sqrt(x*x + y*y);
 629                         dist = TMath::Abs(rCHK - rc);
 630                         cout << rCHK << " vs. " << rc << endl;
 631
 632                         if (TMath::Abs(rCHK - rc) < 0.01) break;
 633                 }
 634                 cout << "STOP" << endl;
 635                 if (!propOK) return distance;
 636                 */
 637         }
 638
 639         // get global propagation of track at end of propagation
 640         tr->GetXYZ(pos);
 641
 642         // check angle cut
 643         TVector3 vc(cl->X(), cl->Y(), cl->Z());
 644         TVector3 vt(x, y, z);
 645         Double_t angle = TMath::Abs(vc.Angle(vt)) * TMath::RadToDeg();
 646         // check: where is the track?
 647         Double_t r, phiT, phiC;
 648         r = TMath::Sqrt(pos[0]*pos[0] + pos[1]*pos[1]);
 649         phiT = TMath::ATan2(pos[1], pos[0]) * TMath::RadToDeg();
 650         phiC = vc.Phi() * TMath::RadToDeg();
 651         //cout << "Propagated R, phiT, phiC = " << r << ' ' << phiT << ' ' << phiC << endl;
 652
 653         if (angle > fCutAngle) {
 654                 //cout << "angle" << endl;
 655           delete tr;
 656                 return distance;
 657         }
 658
 659         // compute differences wr to each coordinate
 660         x -= cl->X();
 661         if (TMath::Abs(x) > fCutX) {
 662                 //cout << "cut X" << endl;
 663           delete tr;
 664                 return distance;
 665         }
 666         y -= cl->Y();
 667         if (TMath::Abs(y) > fCutY) {
 668                 //cout << "cut Y" << endl;
 669           delete tr;
 670                 return distance;
 671         }
 672         z -= cl->Z();
 673         if (TMath::Abs(z) > fCutZ) {
 674                 //cout << "cut Z" << endl;
 675           delete tr;
 676                 return distance;
 677         }
 678
 679         // compute true distance
 680         distance = TMath::Sqrt(x*x + y*y + z*z);
 681         //Double_t temp = CheckPairV2(tr, cl);
 682         //if (temp < distance) return temp; else
 683
 684         // delete temporary object
 685         delete tr;
 686
 687         return distance;
 688 }
 689 //
 690 //------------------------------------------------------------------------------
 691 //
 692 Double_t AliEMCALTracker::CheckPairV2
 693 (AliEMCALTrack *tr, AliEMCALMatchCluster *cl)
 694 {
 695         //
 696         // Given a track and a cluster,
 697         // propagates the first to the radius of the second.
 698         // Then, checks the propagation point against all cuts.
 699         // If at least a cut is not passed, a valuer equal to
 700         // twice the maximum allowed distance is passed (so the value returned
 701         // will not be taken into account when creating matches)
 702         //
 703
 704         // TEMP
 705 //      Bool_t isTrue = kFALSE;
 706 //      if (tr->GetSeedLabel() == cl->Label()) {
 707 //              isTrue = kTRUE;
 708 //              cout << "TRUE MATCH!!!" << endl;
 709 //      }
 710
 711         Double_t distance = 2.0 * fMaxDist;
 712
 713         Double_t x0, rho;
 714         if (fTrackCorrMode == kTrackCorrMMB) {
 715                 Double_t pos1[3], pos2[3], param[6];
 716                 tr->GetXYZ(pos1);
 717 //              pos1[0] = x;
 718 //              pos1[1] = y;
 719 //              pos1[2] = z;
 720                 pos2[0] = cl->X();
 721                 pos2[1] = cl->Y();
 722                 pos2[2] = cl->Z();
 723                 MeanMaterialBudget(pos1, pos2, param);
 724                 rho = param[0]*param[4];
 725                 x0 = param[1];
 726         }
 727         else if (fTrackCorrMode == kTrackCorrFixed) {
 728                 rho = fRho;
 729                 x0 = fX0;
 730         }
 731         else {
 732                 rho = 0.0;
 733                 x0 = 0.0;
 734         }
 735
 736         // check against cut on difference 'alpha - phi'
 737         Double_t phi = TMath::ATan2(cl->Y(), cl->X());
 738         phi = AngleDiff(phi, tr->GetAlpha());
 739         if (phi < fCutAlphaMin || phi > fCutAlphaMax) return distance;
 740
 741         // get cluster position and put them into a vector
 742         TVector3 vc(cl->X(), cl->Y(), cl->Z());
 743         // rotate the vector in order to put all clusters on a plane intersecting
 744         // vertically the X axis; the angle depends on the sector
 745         Double_t clusterRot, clusterPhi = vc.Phi() * TMath::RadToDeg();
 746         if (clusterPhi < 0.0) clusterPhi += 360.0;
 747         if (clusterPhi < 100.0) {
 748                 clusterRot = -90.0;
 749         }
 750         else if (clusterPhi < 120.0) {
 751                 clusterRot = -110.0;
 752         }
 753         else if (clusterPhi < 140.0) {
 754                 clusterRot = -130.0;
 755         }
 756         else if (clusterPhi < 160.0) {
 757                 clusterRot = -150.0;
 758         }
 759         else if (clusterPhi < 180.0) {
 760                 clusterRot = -170.0;
 761         }
 762         else {
 763                 clusterRot = -190.0;
 764         }
 765         vc.RotateZ(clusterRot * TMath::DegToRad());
 766         // generate a track from the ESD track selected
 767         AliEMCALTrack *track = new AliEMCALTrack(*tr);
 768         // compute the 'phi' coordinate of the intersection point to
 769         // the EMCAL surface
 770         Double_t x = vc.X();
 771         Double_t y;
 772         track->GetYAt(vc.X(), track->GetBz(), y);
 773         Double_t tmp = x*TMath::Cos(track->GetAlpha()) - y*TMath::Sin(track->GetAlpha());
 774         y = x*TMath::Sin(track->GetAlpha()) + y*TMath::Cos(track->GetAlpha());
 775         x = tmp;
 776         Double_t trackPhi = TMath::ATan2(y, x) * TMath::RadToDeg();
 777         // compute phi difference
 778         Double_t dphi = trackPhi - clusterPhi;
 779         if (TMath::Abs(dphi) > 180.0) {
 780                 dphi = 360.0 - TMath::Abs(dphi);
 781                 if (clusterPhi > trackPhi) dphi = -dphi;
 782         }
 783         // propagate track to the X position of rotated cluster
 784         // and get the vector of X, Y, Z in the local ref. frame of the track
 785         track->PropagateTo(vc.X(), x0, rho);
 786         TVector3 vt(track->GetX(), track->GetY(), track->GetZ());
 787         vt.RotateZ((clusterPhi - trackPhi) * TMath::DegToRad());
 788         TVector3 vdiff = vt-vc;
 789
 790         // compute differences wr to each coordinate
 791         delete track;
 792         if (vdiff.X() > fCutX) return distance;
 793         if (vdiff.Y() > fCutY) return distance;
 794         if (vdiff.Z() > fCutZ) return distance;
 795
 796         // compute true distance
 797         distance = vdiff.Mag();
 798         return distance;
 799 }
 800 //
 801 //------------------------------------------------------------------------------
 802 //
 803 Double_t AliEMCALTracker::CheckPairV3
 804 (AliEMCALTrack *track, AliEMCALMatchCluster *cl)
 805 {
 806         //
 807         // Given a track and a cluster,
 808         // propagates the first to the radius of the second.
 809         // Then, checks the propagation point against all cuts.
 810         // If at least a cut is not passed, a valuer equal to
 811         // twice the maximum allowed distance is passed (so the value returned
 812         // will not be taken into account when creating matches)
 813         //
 814
 815         AliEMCALTrack tr(*track);
 816
 817         Int_t    sector;
 818         Double_t distance = 2.0 * fMaxDist;
 819         Double_t dx, dy, dz;
 820         Double_t phi, alpha, slope, tgtXnum, tgtXden, sectorWidth = 20.0 * TMath::DegToRad();
 821         Double_t xcurr, xprop, param[6] = {0., 0., 0., 0., 0., 0.}, x0, rho, bz;
 822         Double_t x[3], x1[3], x2[3];
 823
 824         // get initial track position
 825         xcurr = tr.GetX();
 826
 827         // evaluate the EMCAL sector number
 828         phi = cl->Phi();
 829         if (phi < 0.0) phi += TMath::TwoPi();
 830         sector = (Int_t)(phi / sectorWidth);
 831         alpha = ((Double_t)sector + 0.5) * sectorWidth;
 832         // evaluate the corresponding X for track propagation
 833         slope = TMath::Tan(alpha - 0.5*TMath::Pi());
 834         tgtXnum = cl->Y() - slope * cl->X();
 835         tgtXden = TMath::Sqrt(1.0 + slope*slope);
 836         xprop = TMath::Abs(tgtXnum / tgtXden);
 837
 838         // propagate by small steps
 839         tr.GetXYZ(x1);
 840         bz = tr.GetBz();
 841         if (!tr.GetXYZAt(xprop, bz, x2)) return distance;
 842         //AliKalmanTrack::MeanMaterialBudget(x1, x2, param);
 843         rho = param[0]*param[4];
 844         x0 = param[1];
 845         if (!tr.PropagateTo(xprop, x0, rho)) return distance;
 846         //if (!tr.PropagateTo(xprop, 0.0, 0.0)) return distance;
 847
 848         // get propagated position at the end
 849         tr.GetXYZ(x);
 850         dx = TMath::Abs(x[0] - cl->X());
 851         dy = TMath::Abs(x[1] - cl->Y());
 852         dz = TMath::Abs(x[2] - cl->Z());
 853         if (dx > fCutX || dy > fCutY || dz > fCutZ) return distance;
 854
 855         distance = TMath::Sqrt(dx*dx + dy*dy + dz*dz);
 856
 857         return distance;
 858 }
 859 //
 860 //------------------------------------------------------------------------------
 861 //
 862 Bool_t AliEMCALTracker::PropagateToEMCAL(AliEMCALTrack *tr)
 863 {
 864         //
 865         // Propagates the track to the proximity of the EMCAL surface
 866         //
 867
 868         Double_t xcurr, xtemp, xprop = 438.0, step = 10.0, param[6], x0, rho, bz;
 869         Double_t x1[3], x2[3];
 870
 871         // get initial track position
 872         xcurr = tr->GetX();
 873
 874         // propagate by small steps
 875         for (xtemp = xcurr + step; xtemp < xprop; xtemp += step) {
 876                 // to compute material budget, take current position and
 877                 // propagated hypothesis without energy loss
 878                 tr->GetXYZ(x1);
 879                 bz = tr->GetBz();
 880                 if (!tr->GetXYZAt(xtemp, bz, x2)) return kFALSE;
 881                 MeanMaterialBudget(x1, x2, param);
 882                 rho = param[0]*param[4];
 883                 x0 = param[1];
 884                 if (!tr->PropagateTo(xtemp, x0, rho)) return kFALSE;
 885         }
 886
 887         return kTRUE;
 888 }
 889 //
 890 //------------------------------------------------------------------------------
 891 //
 892 Int_t AliEMCALTracker::CreateMatches()
 893 {
 894         //
 895         // Creation of matches between tracks and clusters.
 896         // For each ESD track collected by ReadESD(), an AliEMCALTrack is made.
 897         // If it finds a cluster close enough to its propagation to EMCAL,
 898         // which passes all cuts, its index is stored.
 899         // If many clusters are found which satisfy the criteria described above,
 900         // only the closest one is stored.
 901         // At this level, it is possible that two tracks share the same cluster.
 902         //
 903
 904         // if matches collection is already present, it is deleted
 905         if (fMatches) {
 906                 fMatches->Delete();
 907                 delete fMatches;
 908         }
 909         fMatches = new TList;
 910
 911         // initialize counters and indexes
 912         Int_t count = 0;
 913         Int_t ic, nClusters = (Int_t)fClusters->GetEntries();
 914         Int_t it, nTracks = fTracks->GetEntries();
 915
 916         // external loop on clusters, internal loop on tracks
 917         Double_t dist;
 918         for (ic = 0; ic < nClusters; ic++) {
 919                 AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(ic);
 920                 for (it = 0; it < nTracks; it++) {
 921                         AliEMCALTrack *track = (AliEMCALTrack*)fTracks->At(it);
 922                         dist = CheckPair(track, cluster);
 923                         //cout << dist << endl;
 924                         if (dist <= fMaxDist) {
 925                                 AliEMCALMatch *candidate = new AliEMCALMatch;
 926                                 candidate->SetIndexT(it);
 927                                 candidate->SetIndexC(ic);
 928                                 candidate->SetDistance(dist);
 929                                 candidate->SetPt(track->GetSignedPt());
 930                                 fMatches->Add(candidate);
 931                                 count++;
 932                         }
 933                 }
 934         }
 935
 936         /*
 937         // loop on clusters and tracks
 938         Int_t icBest;
 939         Double_t dist, distBest;
 940         for (it = 0; it < nTracks; it++) {
 941                 AliEMCALTrack *track = (AliEMCALTrack*)fTracks->At(it);
 942                 if (!track) continue;
 943                 icBest = -1;
 944                 distBest = fMaxDist;
 945                 for (ic = 0; ic < nClusters; ic++) {
 946                         AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(ic);
 947                         if (!cluster) continue;
 948                         dist = CheckPair(track, cluster);
 949                         if (dist < distBest) {
 950                                 distBest = dist;
 951                                 icBest = ic;
 952                         }
 953                 }
 954                 if (icBest >= 0) {
 955                         track->SetMatchedClusterIndex(icBest);
 956                         track->SetMatchedClusterDist(distBest);
 957                         count++;
 958                 }
 959                 else {
 960                         track->SetMatchedClusterIndex(-1);
 961                 }
 962         }
 963         */
 964
 965         return count;
 966 }
 967 //
 968 //------------------------------------------------------------------------------
 969 //
 970 Int_t AliEMCALTracker::SolveCompetitions()
 971 {
 972         //
 973         // Match selector.
 974         // The match list is sorted from the best to the worst match, w.r. to the
 975         // distance between track prolongation and cluster position.
 976         // Based on this criterion, starting from the first (best) match, a flag
 977         // is set to both the involved track and cluster, and all matches containing
 978         // an already used track or cluster are removed, leaving only the best match
 979         // for each cluster.
 980         //
 981
 982         // sort matches with respect to track-cluster distance
 983         fMatches->Sort(kSortAscending);
 984
 985         // keep track of eliminated matches
 986         Int_t count = 0;
 987
 988         // initialize flags to check repetitions
 989         Int_t ic, nClusters = (Int_t)fClusters->GetEntries();
 990         Int_t it, nTracks = fTracks->GetEntries();
 991         Bool_t *usedC = new Bool_t[nClusters];
 992         Bool_t *usedT = new Bool_t[nTracks];
 993         for (ic = 0; ic < nClusters; ic++) usedC[ic] = kFALSE;
 994         for (it = 0; it < nTracks; it++) usedT[it] = kFALSE;
 995
 996         // loop on matches
 997         TListIter iter(fMatches);
 998         AliEMCALMatch *match = 0;
 999         while ( (match = (AliEMCALMatch*)iter.Next()) ) {
1000                 ic = match->GetIndexC();
1001                 it = match->GetIndexT();
1002                 if (!usedT[it] && !usedC[ic]) {
1003                         usedT[it] = kTRUE;
1004                         usedC[ic] = kTRUE;
1005                         match->CanBeSaved() = kTRUE;
1006                 }
1007                 else {
1008                         count++;
1009                 }
1010         }
1011
1012         /*
1013         Int_t it1, it2, nTracks = (Int_t)fTracks->GetEntries();
1014         AliEMCALTrack *track1 = 0, *track2 = 0;
1015         for (it1 = 0; it1 < nTracks; it1++) {
1016                 track1 = (AliEMCALTrack*)fTracks->At(it1);
1017                 if (!track1) continue;
1018                 if (track1->GetMatchedClusterIndex() < 0) continue;
1019                 for (it2 = it1+1; it2 < nTracks; it2++) {
1020                         track2 = (AliEMCALTrack*)fTracks->At(it2);
1021                         if (!track2) continue;
1022                         if (track2->GetMatchedClusterIndex() < 0) continue;
1023                         if (track1->GetMatchedClusterIndex() != track2->GetMatchedClusterIndex()) continue;
1024                         count++;
1025                         if (track1->GetMatchedClusterDist() < track2->GetMatchedClusterDist()) {
1026                                 track2->SetMatchedClusterIndex(-1);
1027                         }
1028                         else if (track2->GetMatchedClusterDist() < track1->GetMatchedClusterDist()) {
1029                                 track1->SetMatchedClusterIndex(-1);
1030                         }
1031                 }
1032         }
1033         */
1034
1035         delete [] usedC;
1036         delete [] usedT;
1037
1038         return count;
1039 }
1040 //
1041 //------------------------------------------------------------------------------
1042 //
1043 void AliEMCALTracker::UnloadClusters()
1044 {
1045         //
1046         // Free memory from all arrays
1047         // This method is called after the local tracking step
1048         // so we can safely delete everything
1049         //
1050
1051         Clear();
1052 }
1053 //
1054 //------------------------------------------------------------------------------
1055 //
1056 AliEMCALTracker::AliEMCALMatchCluster::AliEMCALMatchCluster(Int_t index, AliEMCALRecPoint *recPoint)
1057   : fIndex(index),
1058     fLabel(recPoint->GetPrimaryIndex()),
1059     fX(0.),
1060     fY(0.),
1061     fZ(0.)
1062 {
1063         //
1064         // Translates an AliEMCALRecPoint object into the internal format.
1065         // Index of passed cluster in its native array must be specified.
1066         //
1067         TVector3 clpos;
1068         recPoint->GetGlobalPosition(clpos);
1069
1070         fX = clpos.X();
1071         fY = clpos.Y();
1072         fZ = clpos.Z();
1073 }
1074 //
1075 //------------------------------------------------------------------------------
1076 //
1077 AliEMCALTracker::AliEMCALMatchCluster::AliEMCALMatchCluster(Int_t index, AliESDCaloCluster *caloCluster)
1078   : fIndex(index),
1079     fLabel(caloCluster->GetLabel()),
1080     fX(0.),
1081     fY(0.),
1082     fZ(0.)
1083 {
1084         //
1085         // Translates an AliESDCaloCluster object into the internal format.
1086         // Index of passed cluster in its native array must be specified.
1087         //
1088         Float_t clpos[3];
1089         caloCluster->GetPosition(clpos);
1090
1091         fX = (Double_t)clpos[0];
1092         fY = (Double_t)clpos[1];
1093         fZ = (Double_t)clpos[2];
1094 }
1095 //
1096 //------------------------------------------------------------------------------
1097 //
1098 Int_t AliEMCALTracker::AliEMCALMatch::Compare(const TObject *obj) const
1099 {
1100         //
1101         // Tracks compared wrt their distance from matched point
1102         //
1103
1104         AliEMCALTracker::AliEMCALMatch *that = (AliEMCALTracker::AliEMCALMatch*)obj;
1105
1106         Double_t thisDist = fPt;//fDistance;
1107         Double_t thatDist = that->fPt;//that->GetDistance();
1108
1109         if (thisDist > thatDist) return 1;
1110         else if (thisDist < thatDist) return -1;
1111         return 0;
1112 }
1113
1114 AliEMCALTracker::AliEMCALMatch::AliEMCALMatch()
1115   : TObject(),
1116     fCanBeSaved(kFALSE),
1117     fIndexC(0),
1118     fIndexT(0),
1119     fDistance(0.),
1120     fPt(0.)
1121 {
1122   //default constructor
1123
1124 }
1125
1126 AliEMCALTracker::AliEMCALMatch::AliEMCALMatch(const AliEMCALMatch& copy)
1127   : TObject(),
1128     fCanBeSaved(copy.fCanBeSaved),
1129     fIndexC(copy.fIndexC),
1130     fIndexT(copy.fIndexT),
1131     fDistance(copy.fDistance),
1132     fPt(copy.fPt)
1133 {
1134   //copy ctor
1135 }