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 <TMath.h>
  38 #include <TList.h>
  39 #include <TString.h>
  40 #include <TVector3.h>
  41 #include <TClonesArray.h>
  42
  43 #include "AliLog.h"
  44 #include "AliESD.h"
  45 #include "AliESDtrack.h"
  46 #include "AliKalmanTrack.h"
  47 #include "AliEMCALRecPoint.h"
  48 #include "AliRunLoader.h"
  49 #include "AliEMCALTrack.h"
  50 #include "AliEMCALLoader.h"
  51 #include "AliEMCALTracker.h"
  52
  53 ClassImp(AliEMCALTracker)
  54 //
  55 //------------------------------------------------------------------------------
  56 //
  57 AliEMCALTracker::AliEMCALTracker()
  58   : AliTracker(),
  59     fNPropSteps(0),
  60     fTrackCorrMode(kTrackCorrNone),
  61     fCutX(50.0),
  62     fCutY(50.0),
  63     fCutZ(50.0),
  64     fCutAlphaMin(-200.0),
  65     fCutAlphaMax(200.0),
  66     fCutAngle(100.0),
  67     fMaxDist(100.0),
  68     fRho(1.0),
  69     fX0(1.0),
  70     fTracks(0),
  71     fClusters(0),
  72     fMatches(0)
  73 {
  74         //
  75         // Default constructor.
  76         // Initializes al simple data members to default values,
  77         // and all collections to NULL.
  78         // Output file name is set to a default value.
  79         //
  80 }
  81 //
  82 //------------------------------------------------------------------------------
  83 //
  84 AliEMCALTracker::AliEMCALTracker(const AliEMCALTracker& copy)
  85   : AliTracker(),
  86     fNPropSteps(copy.fNPropSteps),
  87     fTrackCorrMode(copy.fTrackCorrMode),
  88     fCutX(copy.fCutX),
  89     fCutY(copy.fCutY),
  90     fCutZ(copy.fCutZ),
  91     fCutAlphaMin(copy.fCutAlphaMin),
  92     fCutAlphaMax(copy.fCutAlphaMax),
  93     fCutAngle(copy.fCutAngle),
  94     fMaxDist(copy.fMaxDist),
  95     fRho(copy.fRho),
  96     fX0(copy.fX0),
  97     fTracks((TObjArray*)copy.fTracks->Clone()),
  98     fClusters((TObjArray*)copy.fClusters->Clone()),
  99     fMatches((TList*)copy.fMatches->Clone())
 100 {
 101         //
 102         // Copy constructor
 103         // Besides copying all parameters, duplicates all collections.
 104         //
 105 }
 106 //
 107 //------------------------------------------------------------------------------
 108 //
 109 AliEMCALTracker& AliEMCALTracker::operator=(const AliEMCALTracker& copy)
 110 {
 111         //
 112         // Assignment operator.
 113         // Besides copying all parameters, duplicates all collections.
 114         //
 115
 116         fCutX = copy.fCutX;
 117         fCutY = copy.fCutY;
 118         fCutZ = copy.fCutZ;
 119         fCutAlphaMin = copy.fCutAlphaMin;
 120         fCutAlphaMax = copy.fCutAlphaMax;
 121         fCutAngle = copy.fCutAngle;
 122         fMaxDist = copy.fMaxDist;
 123
 124         fTracks = (TObjArray*)copy.fTracks->Clone();
 125         fClusters = (TObjArray*)copy.fClusters->Clone();
 126         fMatches = (TList*)copy.fMatches->Clone();
 127
 128         return (*this);
 129 }
 130 //
 131 //------------------------------------------------------------------------------
 132 //
 133 void AliEMCALTracker::Clear(Option_t* option)
 134 {
 135         //
 136         // Clearing method
 137         // Clears all specified arrays and the containers themselves.
 138         //
 139
 140         TString opt(option);
 141         Bool_t clearTracks = opt.Contains("TRACKS");
 142         Bool_t clearClusters = opt.Contains("CLUSTERS");
 143         Bool_t clearMatches = opt.Contains("MATCHES");
 144         if (opt.Contains("ALL")) {
 145                 clearTracks = kTRUE;
 146                 clearClusters = kTRUE;
 147                 clearMatches = kTRUE;
 148         }
 149
 150         if (fTracks != 0x0 && clearTracks) {
 151                 if (!fTracks->IsEmpty()) fTracks->Delete();
 152                 delete fTracks;
 153                 fTracks = 0;
 154         }
 155         if (fClusters != 0x0 && clearClusters) {
 156                 if (!fClusters->IsEmpty()) fClusters->Delete();
 157                 delete fClusters;
 158                 fClusters = 0;
 159         }
 160         if (fMatches != 0x0 && clearMatches) {
 161                 if (!fMatches->IsEmpty()) fMatches->Delete();
 162                 delete fMatches;
 163                 fMatches = 0;
 164         }
 165 }
 166 //
 167 //------------------------------------------------------------------------------
 168 //
 169 Int_t AliEMCALTracker::LoadClusters(TTree *cTree)
 170 {
 171         //
 172         // Load EMCAL clusters in the form of AliEMCALRecPoint,
 173         // from simulation temporary files.
 174         // (When included in reconstruction chain, this method is used automatically)
 175         //
 176
 177         Clear("CLUSTERS");
 178
 179         TBranch *branch = cTree->GetBranch("EMCALECARP");
 180         if (!branch) {
 181                 AliError("can't get the branch with the EMCAL clusters");
 182                 return 1;
 183         }
 184
 185         TClonesArray dummy("AliEMCALRecPoint", 10000);
 186         TClonesArray *clusters = &dummy;
 187         branch->SetAddress(&clusters);
 188         Int_t nClusters = (Int_t)clusters->GetEntries();
 189
 190         cTree->GetEvent(0);
 191         fClusters = new TObjArray(0);
 192         for (Int_t i = 0; i < nClusters; i++) {
 193                 AliEMCALRecPoint *cluster = (AliEMCALRecPoint*)clusters->At(i);
 194                 if (!cluster) continue;
 195                 if (cluster->GetClusterType() != AliESDCaloCluster::kClusterv1) continue;
 196                 AliEMCALMatchCluster *matchCluster = new AliEMCALMatchCluster(i, cluster);
 197                 fClusters->AddLast(matchCluster);
 198         }
 199         if (fClusters->IsEmpty()) {
 200                 AliError("No clusters collected");
 201                 return 1;
 202         }
 203
 204         AliInfo(Form("Collected %d clusters", fClusters->GetEntries()));
 205
 206         return 0;
 207 }
 208 //
 209 //------------------------------------------------------------------------------
 210 //
 211 Int_t AliEMCALTracker::LoadClusters(AliESD *esd)
 212 {
 213         //
 214         // Load EMCAL clusters in the form of AliESDCaloClusters,
 215         // from an AliESD object.
 216         //
 217
 218         // make sure that tracks/clusters collections are empty
 219         Clear("CLUSTERS");
 220
 221         Int_t start = esd->GetFirstEMCALCluster();
 222         Int_t nClustersEMC = esd->GetNumberOfEMCALClusters();
 223         Int_t end = start + nClustersEMC;
 224
 225         fClusters = new TObjArray(0);
 226
 227         Int_t i;
 228         for (i = start; i < end; i++) {
 229                 AliESDCaloCluster *cluster = esd->GetCaloCluster(i);
 230                 if (!cluster) continue;
 231                 if (cluster->GetClusterType() != AliESDCaloCluster::kClusterv1) continue;
 232                 AliEMCALMatchCluster *matchCluster = new AliEMCALMatchCluster(i, cluster);
 233                 fClusters->AddLast(matchCluster);
 234         }
 235         if (fClusters->IsEmpty()) {
 236                 AliError("No clusters collected");
 237                 return 1;
 238         }
 239
 240         AliInfo(Form("Collected %d clusters", fClusters->GetEntries()));
 241
 242         return 0;
 243 }
 244 //
 245 //------------------------------------------------------------------------------
 246 //
 247 Int_t AliEMCALTracker::LoadTracks(AliESD *esd)
 248 {
 249         //
 250         // Load ESD tracks.
 251         //
 252
 253         Clear("TRACKS");
 254
 255         Int_t nTracks = esd->GetNumberOfTracks();
 256         fTracks = new TObjArray(0);
 257
 258         Int_t i, j;
 259         Bool_t isKink;
 260         Double_t alpha;
 261         for (i = 0; i < nTracks; i++) {
 262                 AliESDtrack *esdTrack = esd->GetTrack(i);
 263                 // set by default the value corresponding to "no match"
 264                 esdTrack->SetEMCALcluster(-99999);
 265 //              if (esdTrack->GetLabel() < 0) continue;
 266 //              if (!(esdTrack->GetStatus() & AliESDtrack::kTOFout)) continue;
 267                 isKink = kFALSE;
 268                 for (j = 0; j < 3; j++) {
 269                         if (esdTrack->GetKinkIndex(j) > 0) isKink = kTRUE;
 270                 }
 271                 if (isKink) continue;
 272                 AliEMCALTrack *track = new AliEMCALTrack(*esdTrack);
 273                 // check alpha and reject the tracks which fall outside EMCAL acceptance
 274                 alpha = track->GetAlpha() * TMath::RadToDeg();
 275                 if (alpha >  -155.0 && alpha < 67.0) {
 276                         delete track;
 277                         continue;
 278                 }
 279                 track->SetSeedIndex(i);
 280                 track->SetSeedLabel(esdTrack->GetLabel());
 281                 fTracks->AddLast(track);
 282         }
 283         if (fTracks->IsEmpty()) {
 284                 AliError("No tracks collected");
 285                 return 1;
 286         }
 287
 288         AliInfo(Form("Collected %d tracks", fTracks->GetEntries()));
 289
 290         return 0;
 291 }
 292 //
 293 //------------------------------------------------------------------------------
 294 //
 295 Int_t AliEMCALTracker::PropagateBack(AliESD* esd)
 296 {
 297         //
 298         // Main operation method.
 299         // Gets external AliESD containing tracks to be matched.
 300         // After executing match finding, stores in the same ESD object all infos
 301         // and releases the object for further reconstruction steps.
 302         //
 303
 304         if (!esd) {
 305                 AliError("NULL ESD passed");
 306                 return 1;
 307         }
 308
 309         // step 1:
 310         // if cluster array is empty, cluster are collected
 311         // from the passed ESD, and work is done with ESDCaloClusters
 312         Int_t okLoadClusters;
 313         if (!fClusters || (fClusters && fClusters->IsEmpty())) {
 314                 AliInfo("Cluster array is empty. Loading clusters...");
 315                 okLoadClusters = LoadClusters(esd);
 316                 if (okLoadClusters) return 2;
 317         }
 318
 319         // step 2:
 320         // collect ESD tracks
 321         Int_t okLoadTracks = LoadTracks(esd);
 322         if (okLoadTracks) return 3;
 323
 324         // step 3:
 325         // each track is propagated to the "R" position of each cluster.
 326         // The closest cluster is assigned as match.
 327         // IF no clusters lie within the maximum allowed distance, no matches are assigned.
 328         Int_t nMatches = CreateMatches();
 329         if (!nMatches) {
 330                 AliInfo("No good matches found.");
 331                 return 4;
 332         }
 333         else {
 334                 AliInfo(Form("Found %d matches", nMatches));
 335         }
 336
 337         // step 4:
 338         // when more than 1 track share the same matched cluster, only the closest one is kept.
 339         Int_t nRemoved = SolveCompetitions();
 340         AliInfo(Form("Removed %d duplicate matches", nRemoved));
 341         if (nRemoved >= nMatches) {
 342                 AliError("Removed ALL matches! Check the algorithm or data. Nothing to save");
 343                 return 5;
 344         }
 345
 346         // step 5:
 347         // save obtained information setting the 'fEMCALindex' field of AliESDtrack object
 348         Int_t nSaved = 0, trackID;
 349         TListIter iter(fMatches);
 350         AliEMCALMatch *match = 0;
 351         while ( (match = (AliEMCALMatch*)iter.Next()) ) {
 352                 if (!match->CanBeSaved()) continue;
 353                 AliEMCALTrack *track = (AliEMCALTrack*)fTracks->At(match->GetIndexT());
 354                 AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(match->GetIndexC());
 355                 trackID = track->GetSeedIndex();
 356                 AliESDtrack *esdTrack = esd->GetTrack(trackID);
 357                 if (!esdTrack) continue;
 358                 if (esdTrack->GetLabel() == cluster->Label()) {
 359                         esdTrack->SetEMCALcluster(cluster->Index());
 360                 }
 361                 else {
 362                         esdTrack->SetEMCALcluster(-cluster->Index());
 363                 }
 364                 nSaved++;
 365         }
 366         /*
 367         AliEMCALTrack *track = 0;
 368         TObjArrayIter tracks(fTracks);
 369         while ( (track = (AliEMCALTrack*)tracks.Next()) ) {
 370                 trackID = track->GetSeedIndex();
 371                 clusterID = track->GetMatchedClusterIndex();
 372                 AliESDtrack *esdTrack = esd->GetTrack(trackID);
 373                 if (!esdTrack) continue;
 374                 if (clusterID < 0) {
 375                         esdTrack->SetEMCALcluster(-99999);
 376                 }
 377                 else {
 378                         AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(clusterID);
 379                         if (!cluster) continue;
 380                         if (esdTrack->GetLabel() == cluster->Label()) {
 381                                 nGood++;
 382                                 esdTrack->SetEMCALcluster(cluster->Index());
 383                         }
 384                         else {
 385                                 esdTrack->SetEMCALcluster(-cluster->Index());
 386                         }
 387                         nSaved++;
 388                 }
 389         }
 390         */
 391         AliInfo(Form("Saved %d matches", nSaved));
 392
 393         return 0;
 394 }
 395 //
 396 //------------------------------------------------------------------------------
 397 //
 398 void AliEMCALTracker::SetTrackCorrectionMode(Option_t *option)
 399 {
 400         //
 401         // Set track correction mode
 402         // gest the choice in string format and converts into
 403         // internal enum
 404         //
 405
 406         TString opt(option);
 407         opt.ToUpper();
 408
 409         if (!opt.CompareTo("NONE")) {
 410                 fTrackCorrMode = kTrackCorrNone;
 411         }
 412         else if (!opt.CompareTo("MMB")) {
 413                 fTrackCorrMode = kTrackCorrMMB;
 414         }
 415         else if (!opt.CompareTo("FIXED")) {
 416                 fTrackCorrMode = kTrackCorrFixed;
 417         }
 418         else {
 419                 cerr << "E-AliEMCALTracker::SetTrackCorrectionMode '" << option << "': Unrecognized option" << endl;
 420         }
 421 }
 422
 423 //
 424 //------------------------------------------------------------------------------
 425 //
 426 Double_t AliEMCALTracker::AngleDiff(Double_t angle1, Double_t angle2)
 427 {
 428         //
 429         // [PRIVATE]
 430         // Given two angles in radiants, it converts them in the range 0-2pi
 431         // then computes their true difference, i.e. if the difference a1-a2
 432         // results to be larger than 180 degrees, it returns 360 - diff.
 433         //
 434
 435         if (angle1 < 0.0) angle1 += TMath::TwoPi();
 436         if (angle1 > TMath::TwoPi()) angle1 -= TMath::TwoPi();
 437         if (angle2 < 0.0) angle2 += TMath::TwoPi();
 438         if (angle2 > TMath::TwoPi()) angle2 -= TMath::TwoPi();
 439
 440         Double_t diff = TMath::Abs(angle1 - angle2);
 441         if (diff > TMath::Pi()) diff = TMath::TwoPi() - diff;
 442
 443         if (angle2 > angle1) diff = -diff;
 444
 445         return diff;
 446 }
 447 //
 448 //------------------------------------------------------------------------------
 449 //
 450 Double_t AliEMCALTracker::CheckPair
 451 (AliEMCALTrack *track, AliEMCALMatchCluster *cl)
 452 {
 453         //
 454         // Given a track and a cluster,
 455         // propagates the first to the radius of the second.
 456         // Then, checks the propagation point against all cuts.
 457         // If at least a cut is not passed, a valuer equal to
 458         // twice the maximum allowed distance is passed (so the value returned
 459         // will not be taken into account when creating matches)
 460         //
 461
 462         // TEMP
 463         Bool_t isTrue = kFALSE;
 464 //      if (tr->GetSeedLabel() == cl->Label()) {
 465 //              isTrue = kTRUE;
 466 //              cout << "TRUE MATCH!!!" << endl;
 467 //      }
 468
 469         // copy track into temporary variable
 470         AliEMCALTrack *tr = new AliEMCALTrack(*track);
 471
 472         Double_t distance = 2.0 * fMaxDist;
 473
 474         // check against cut on difference 'alpha - phi'
 475         Double_t phi = TMath::ATan2(cl->Y(), cl->X());
 476         phi = AngleDiff(phi, tr->GetAlpha());
 477         if (phi < fCutAlphaMin || phi > fCutAlphaMax) return distance;
 478
 479         // try to propagate to cluster radius
 480         // (return the 'distance' value if it fails)
 481         Double_t pos[3], &x = pos[0], &y = pos[1], &z = pos[2];
 482         Double_t x0, rho;
 483         tr->GetXYZ(pos);
 484         Double_t rt = TMath::Sqrt(x*x + y*y);
 485         Double_t rc = TMath::Sqrt(cl->X()*cl->X() + cl->Y()*cl->Y());
 486
 487         if (fTrackCorrMode == kTrackCorrMMB) {
 488                 Double_t pos1[3], pos2[3], param[6];
 489                 pos1[0] = x;
 490                 pos1[1] = y;
 491                 pos1[2] = z;
 492                 pos2[0] = cl->X();
 493                 pos2[1] = cl->Y();
 494                 pos2[2] = cl->Z();
 495                 AliKalmanTrack::MeanMaterialBudget(pos1, pos2, param);
 496                 rho = param[0];
 497                 x0 = param[1];
 498         }
 499         else if (fTrackCorrMode == kTrackCorrFixed) {
 500                 rho = fRho;
 501                 x0 = fX0;
 502         }
 503         else {
 504                 rho = 0.0;
 505                 x0 = 0.0;
 506         }
 507         if (fNPropSteps) {
 508                 Int_t i;
 509                 Double_t r;
 510                 cout.setf(ios::fixed);
 511                 cout.precision(5);
 512                 if (isTrue) cout << "Init : " << rt << ' ' << x << ' ' << y << ' ' << z << endl;
 513                 for (i = 0; i < fNPropSteps; i++) {
 514                         r = rt + (rc - rt) * ((Double_t)(i+1)/(Double_t)fNPropSteps);
 515                         if (!tr->PropagateTo(r, x0, rho)) return distance;
 516                         tr->GetXYZ(pos);
 517                         if (isTrue) cout << "Step : " << r << ' ' << x << ' ' << y << ' ' << z << endl;
 518                 }
 519                 if (isTrue) cout << "Clstr: " << rc << ' ' << cl->X() << ' ' << cl->Y() << ' ' << cl->Z() << endl;
 520         }
 521         else {
 522                 // when no steps are used, no correction makes sense
 523                 if (!tr->PropagateTo(rc, 0.0, 0.0)) return distance;
 524         }
 525
 526         // get global propagation of track at end of propagation
 527         tr->GetXYZ(pos);
 528
 529         // check angle cut
 530         TVector3 vc(cl->X(), cl->Y(), cl->Z());
 531         TVector3 vt(x, y, z);
 532         Double_t angle = TMath::Abs(vc.Angle(vt)) * TMath::RadToDeg();
 533         if (angle > fCutAngle) {
 534                 //cout << "angle" << endl;
 535                 return distance;
 536         }
 537
 538         // compute differences wr to each coordinate
 539         x -= cl->X();
 540         if (x > fCutX) {
 541                 //cout << "cut X" << endl;
 542                 return distance;
 543         }
 544         y -= cl->Y();
 545         if (y > fCutY) {
 546                 //cout << "cut Y" << endl;
 547                 return distance;
 548         }
 549         z -= cl->Z();
 550         if (z > fCutZ) {
 551                 //cout << "cut Z" << endl;
 552                 return distance;
 553         }
 554
 555         // compute true distance
 556         distance = TMath::Sqrt(x*x + y*y + z*z);
 557         //Double_t temp = CheckPairV2(tr, cl);
 558         //if (temp < distance) return temp; else
 559
 560         // delete temporary object
 561         delete tr;
 562
 563         return distance;
 564 }
 565 //
 566 //------------------------------------------------------------------------------
 567 //
 568 Double_t AliEMCALTracker::CheckPairV2
 569 (AliEMCALTrack *tr, AliEMCALMatchCluster *cl)
 570 {
 571         //
 572         // Given a track and a cluster,
 573         // propagates the first to the radius of the second.
 574         // Then, checks the propagation point against all cuts.
 575         // If at least a cut is not passed, a valuer equal to
 576         // twice the maximum allowed distance is passed (so the value returned
 577         // will not be taken into account when creating matches)
 578         //
 579
 580         // TEMP
 581 //      Bool_t isTrue = kFALSE;
 582 //      if (tr->GetSeedLabel() == cl->Label()) {
 583 //              isTrue = kTRUE;
 584 //              cout << "TRUE MATCH!!!" << endl;
 585 //      }
 586
 587         Double_t distance = 2.0 * fMaxDist;
 588
 589         Double_t x0, rho;
 590         if (fTrackCorrMode == kTrackCorrMMB) {
 591                 Double_t pos1[3], pos2[3], param[6];
 592                 tr->GetXYZ(pos1);
 593 //              pos1[0] = x;
 594 //              pos1[1] = y;
 595 //              pos1[2] = z;
 596                 pos2[0] = cl->X();
 597                 pos2[1] = cl->Y();
 598                 pos2[2] = cl->Z();
 599                 AliKalmanTrack::MeanMaterialBudget(pos1, pos2, param);
 600                 rho = param[0];
 601                 x0 = param[1];
 602         }
 603         else if (fTrackCorrMode == kTrackCorrFixed) {
 604                 rho = fRho;
 605                 x0 = fX0;
 606         }
 607         else {
 608                 rho = 0.0;
 609                 x0 = 0.0;
 610         }
 611
 612         // check against cut on difference 'alpha - phi'
 613         Double_t phi = TMath::ATan2(cl->Y(), cl->X());
 614         phi = AngleDiff(phi, tr->GetAlpha());
 615         if (phi < fCutAlphaMin || phi > fCutAlphaMax) return distance;
 616
 617         // get cluster position and put them into a vector
 618         TVector3 vc(cl->X(), cl->Y(), cl->Z());
 619         // rotate the vector in order to put all clusters on a plane intersecting
 620         // vertically the X axis; the angle depends on the sector
 621         Double_t clusterRot, clusterPhi = vc.Phi() * TMath::RadToDeg();
 622         if (clusterPhi < 0.0) clusterPhi += 360.0;
 623         if (clusterPhi < 100.0) {
 624                 clusterRot = -90.0;
 625         }
 626         else if (clusterPhi < 120.0) {
 627                 clusterRot = -110.0;
 628         }
 629         else if (clusterPhi < 140.0) {
 630                 clusterRot = -130.0;
 631         }
 632         else if (clusterPhi < 160.0) {
 633                 clusterRot = -150.0;
 634         }
 635         else if (clusterPhi < 180.0) {
 636                 clusterRot = -170.0;
 637         }
 638         else {
 639                 clusterRot = -190.0;
 640         }
 641         vc.RotateZ(clusterRot * TMath::DegToRad());
 642         // generate a track from the ESD track selected
 643         AliEMCALTrack *track = new AliEMCALTrack(*tr);
 644         // compute the 'phi' coordinate of the intersection point to
 645         // the EMCAL surface
 646         Double_t x = vc.X();
 647         Double_t y;
 648         track->GetYAt(vc.X(), track->GetBz(), y);
 649         Double_t tmp = x*TMath::Cos(track->GetAlpha()) - y*TMath::Sin(track->GetAlpha());
 650         y = x*TMath::Sin(track->GetAlpha()) + y*TMath::Cos(track->GetAlpha());
 651         x = tmp;
 652         Double_t trackPhi = TMath::ATan2(y, x) * TMath::RadToDeg();
 653         // compute phi difference
 654         Double_t dphi = trackPhi - clusterPhi;
 655         if (TMath::Abs(dphi) > 180.0) {
 656                 dphi = 360.0 - TMath::Abs(dphi);
 657                 if (clusterPhi > trackPhi) dphi = -dphi;
 658         }
 659         // propagate track to the X position of rotated cluster
 660         // and get the vector of X, Y, Z in the local ref. frame of the track
 661         track->PropagateTo(vc.X(), x0, rho);
 662         TVector3 vt(track->GetX(), track->GetY(), track->GetZ());
 663         vt.RotateZ((clusterPhi - trackPhi) * TMath::DegToRad());
 664         TVector3 vdiff = vt-vc;
 665
 666         // compute differences wr to each coordinate
 667         if (vdiff.X() > fCutX) return distance;
 668         if (vdiff.Y() > fCutY) return distance;
 669         if (vdiff.Z() > fCutZ) return distance;
 670
 671         // compute true distance
 672         distance = vdiff.Mag();
 673         return distance;
 674 }
 675 //
 676 //------------------------------------------------------------------------------
 677 //
 678 Int_t AliEMCALTracker::CreateMatches()
 679 {
 680         //
 681         // Creation of matches between tracks and clusters.
 682         // For each ESD track collected by ReadESD(), an AliEMCALTrack is made.
 683         // If it finds a cluster close enough to its propagation to EMCAL,
 684         // which passes all cuts, its index is stored.
 685         // If many clusters are found which satisfy the criteria described above,
 686         // only the closest one is stored.
 687         // At this level, it is possible that two tracks share the same cluster.
 688         //
 689
 690         // if matches collection is already present, it is deleted
 691         if (fMatches) {
 692                 fMatches->Delete();
 693                 delete fMatches;
 694         }
 695         fMatches = new TList;
 696
 697         // initialize counters and indexes
 698         Int_t count = 0;
 699         Int_t ic, nClusters = (Int_t)fClusters->GetEntries();
 700         Int_t it, nTracks = fTracks->GetEntries();
 701
 702         // external loop on clusters, internal loop on tracks
 703         Double_t dist;
 704         for (ic = 0; ic < nClusters; ic++) {
 705                 cout << "\rMatching cluster " << ic+1 << " of " << nClusters << flush;
 706                 AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(ic);
 707                 for (it = 0; it < nTracks; it++) {
 708                         AliEMCALTrack *track = (AliEMCALTrack*)fTracks->At(it);
 709                         dist = CheckPair(track, cluster);
 710                         if (dist <= fMaxDist) {
 711                                 AliEMCALMatch *candidate = new AliEMCALMatch;
 712                                 candidate->SetIndexT(it);
 713                                 candidate->SetIndexC(ic);
 714                                 candidate->SetDistance(dist);
 715                                 fMatches->Add(candidate);
 716                                 count++;
 717                         }
 718                 }
 719         }
 720         cout << endl;
 721
 722         /*
 723         // loop on clusters and tracks
 724         Int_t icBest;
 725         Double_t dist, distBest;
 726         for (it = 0; it < nTracks; it++) {
 727                 AliEMCALTrack *track = (AliEMCALTrack*)fTracks->At(it);
 728                 if (!track) continue;
 729                 icBest = -1;
 730                 distBest = fMaxDist;
 731                 for (ic = 0; ic < nClusters; ic++) {
 732                         AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(ic);
 733                         if (!cluster) continue;
 734                         dist = CheckPair(track, cluster);
 735                         if (dist < distBest) {
 736                                 distBest = dist;
 737                                 icBest = ic;
 738                         }
 739                 }
 740                 if (icBest >= 0) {
 741                         track->SetMatchedClusterIndex(icBest);
 742                         track->SetMatchedClusterDist(distBest);
 743                         count++;
 744                 }
 745                 else {
 746                         track->SetMatchedClusterIndex(-1);
 747                 }
 748         }
 749         */
 750
 751         return count;
 752 }
 753 //
 754 //------------------------------------------------------------------------------
 755 //
 756 Int_t AliEMCALTracker::SolveCompetitions()
 757 {
 758         //
 759         // Match selector.
 760         // The match list is sorted from the best to the worst match, w.r. to the
 761         // distance between track prolongation and cluster position.
 762         // Based on this criterion, starting from the first (best) match, a flag
 763         // is set to both the involved track and cluster, and all matches containing
 764         // an already used track or cluster are removed, leaving only the best match
 765         // for each cluster.
 766         //
 767
 768         // sort matches with respect to track-cluster distance
 769         fMatches->Sort(kSortAscending);
 770
 771         // keep track of eliminated matches
 772         Int_t count = 0;
 773
 774         // initialize flags to check repetitions
 775         Int_t ic, nClusters = (Int_t)fClusters->GetEntries();
 776         Int_t it, nTracks = fTracks->GetEntries();
 777         Bool_t *usedC = new Bool_t[nClusters];
 778         Bool_t *usedT = new Bool_t[nTracks];
 779         for (ic = 0; ic < nClusters; ic++) usedC[ic] = kFALSE;
 780         for (it = 0; it < nTracks; it++) usedT[it] = kFALSE;
 781
 782         // loop on matches
 783         TListIter iter(fMatches);
 784         AliEMCALMatch *match = 0;
 785         while ( (match = (AliEMCALMatch*)iter.Next()) ) {
 786                 ic = match->GetIndexC();
 787                 it = match->GetIndexT();
 788                 if (!usedT[it] && !usedC[ic]) {
 789                         usedT[it] = kTRUE;
 790                         usedC[ic] = kTRUE;
 791                         match->CanBeSaved() = kTRUE;
 792                 }
 793                 else {
 794                         count++;
 795                 }
 796         }
 797
 798         /*
 799         Int_t it1, it2, nTracks = (Int_t)fTracks->GetEntries();
 800         AliEMCALTrack *track1 = 0, *track2 = 0;
 801         for (it1 = 0; it1 < nTracks; it1++) {
 802                 track1 = (AliEMCALTrack*)fTracks->At(it1);
 803                 if (!track1) continue;
 804                 if (track1->GetMatchedClusterIndex() < 0) continue;
 805                 for (it2 = it1+1; it2 < nTracks; it2++) {
 806                         track2 = (AliEMCALTrack*)fTracks->At(it2);
 807                         if (!track2) continue;
 808                         if (track2->GetMatchedClusterIndex() < 0) continue;
 809                         if (track1->GetMatchedClusterIndex() != track2->GetMatchedClusterIndex()) continue;
 810                         count++;
 811                         if (track1->GetMatchedClusterDist() < track2->GetMatchedClusterDist()) {
 812                                 track2->SetMatchedClusterIndex(-1);
 813                         }
 814                         else if (track2->GetMatchedClusterDist() < track1->GetMatchedClusterDist()) {
 815                                 track1->SetMatchedClusterIndex(-1);
 816                         }
 817                 }
 818         }
 819         */
 820
 821         return count;
 822 }
 823 //
 824 //------------------------------------------------------------------------------
 825 //
 826 void AliEMCALTracker::UnloadClusters()
 827 {
 828         //
 829         // Free memory from clusters
 830         //
 831
 832         Clear("CLUSTERS");
 833 }
 834 //
 835 //------------------------------------------------------------------------------
 836 //
 837 AliEMCALTracker::AliEMCALMatchCluster::AliEMCALMatchCluster(Int_t index, AliEMCALRecPoint *recPoint)
 838   : fIndex(index),
 839     fLabel(recPoint->GetPrimaryIndex()),
 840     fX(0.),
 841     fY(0.),
 842     fZ(0.)
 843 {
 844         //
 845         // Translates an AliEMCALRecPoint object into the internal format.
 846         // Index of passed cluster in its native array must be specified.
 847         //
 848         TVector3 clpos;
 849         recPoint->GetGlobalPosition(clpos);
 850
 851         fX = clpos.X();
 852         fY = clpos.Y();
 853         fZ = clpos.Z();
 854 }
 855 //
 856 //------------------------------------------------------------------------------
 857 //
 858 AliEMCALTracker::AliEMCALMatchCluster::AliEMCALMatchCluster(Int_t index, AliESDCaloCluster *caloCluster)
 859   : fIndex(index),
 860     fLabel(caloCluster->GetPrimaryIndex()),
 861     fX(0.),
 862     fY(0.),
 863     fZ(0.)
 864 {
 865         //
 866         // Translates an AliESDCaloCluster object into the internal format.
 867         // Index of passed cluster in its native array must be specified.
 868         //
 869         Float_t clpos[3];
 870         caloCluster->GetGlobalPosition(clpos);
 871
 872         fX = (Double_t)clpos[0];
 873         fY = (Double_t)clpos[1];
 874         fZ = (Double_t)clpos[2];
 875 }
 876 //
 877 //------------------------------------------------------------------------------
 878 //
 879 Int_t AliEMCALTracker::AliEMCALMatch::Compare(const TObject *obj) const
 880 {
 881         //
 882         // Tracks compared wrt their distance from matched point
 883         //
 884
 885         AliEMCALTracker::AliEMCALMatch *that = (AliEMCALTracker::AliEMCALMatch*)obj;
 886
 887         Double_t thisDist = fDistance;
 888         Double_t thatDist = that->GetDistance();
 889
 890         if (thisDist > thatDist) return 1;
 891         else if (thisDist < thatDist) return -1;
 892         return 0;
 893 }
 894
 895 AliEMCALTracker::AliEMCALMatch::AliEMCALMatch()
 896   : TObject(),
 897     fCanBeSaved(kFALSE),
 898     fIndexC(0),
 899     fIndexT(0),
 900     fDistance(0.)
 901 {
 902   //default constructor
 903
 904 }
 905
 906 AliEMCALTracker::AliEMCALMatch::AliEMCALMatch(const AliEMCALMatch& copy)
 907   : TObject(),
 908     fCanBeSaved(copy.fCanBeSaved),
 909     fIndexC(copy.fIndexC),
 910     fIndexT(copy.fIndexT),
 911     fDistance(copy.fDistance)
 912 {
 913   //copy ctor
 914 }