PWG3/muon/AliAnalysisTaskLinkToMC.cxx

   1 /**************************************************************************
   2  * This file is property of and copyright by the ALICE HLT Project        *
   3  * All rights reserved.                                                   *
   4  *                                                                        *
   5  * Primary Authors:                                                       *
   6  *   Artur Szostak <artursz@iafrica.com>                                  *
   7  *                                                                        *
   8  * Permission to use, copy, modify and distribute this software and its   *
   9  * documentation strictly for non-commercial purposes is hereby granted   *
  10  * without fee, provided that the above copyright notice appears in all   *
  11  * copies and that both the copyright notice and this permission notice   *
  12  * appear in the supporting documentation. The authors make no claims     *
  13  * about the suitability of this software for any purpose. It is          *
  14  * provided "as is" without express or implied warranty.                  *
  15  **************************************************************************/
  16
  17 ///
  18 /// @file   AliAnalysisTaskLinkToMC.cxx
  19 /// @author Artur Szostak <artursz@iafrica.com>
  20 /// @date   27 Oct 2008
  21 /// @brief  Implementation of the AliAnalysisTaskLinkToMC task.
  22 ///
  23 /// The AliAnalysisTaskLinkToMC connects reconstructed tracks found in the ESD
  24 /// to their corresponding Monte Carlo (MC) tracks via the proximity of their
  25 /// hit coordinates. Each ESD track is copied over to an AOD track and put into
  26 /// the AOD event container, but the AOD track label is filled such that the
  27 /// AOD track label is the same as the MC track label.
  28 /// If a MC track could not be found for a given ESD track then -1 is filled
  29 /// into the AOD track label.
  30 /// The corresponding MC track is found by correlating the reconstructed track
  31 /// clusters with the MC track references. For each cluster the closest track
  32 /// reference coordinate is found and a chi-squared value is computed for each
  33 /// X and Y dimension independently. The chi-square's are summed up and the MC
  34 /// track with the smallest total chi-squared is chosen as the best match to
  35 /// the ESD track.
  36 ///
  37 /// The following cuts control the conditions under which we consider a MC track
  38 /// to match a given ESD track:
  39 /// 1) SigmaCut(Double_t value) - This is maximum distance expressed in standard
  40 ///      deviations that is allowed between a reconstructed cluster coordinate
  41 ///      and the MC track reference coordinate. The condition is expressed as:
  42 ///        (x_{cluster} - x_{trackref})^2 / (sigma_{x})^2 < #sigma_{cut}
  43 ///      This conditions is applied to each X and Y dimension respectively and
  44 ///      cluster / track reference pairs that do not pass this cut are considered
  45 ///      not to match.
  46 /// 2) MinClusters(Int_t value) - Indicates the minimum number of clusters that
  47 ///      must be matched to a track reference from a MC track for the ESD and
  48 ///      MC track to be considered to match.
  49 /// The following are absolute cuts which are used mainly as a speed optimisation
  50 /// so as to skip over cluster / track reference pairs that are far away from
  51 /// each other early.
  52 /// 3) HardCutLimitX(Double_t value) - The absolute maximum distance in centimetres
  53 ///      allowed between a cluster's X coordinate and a track reference's X coordinate.
  54 ///      If a cluster / track reference pair does not pass this condition then the
  55 ///      pair is considered not to match.
  56 /// 4) HardCutLimitY(Double_t value) - Similar to HardCutLimitX but applied to the
  57 ///      Y dimension of the cluster and track reference coordinate.
  58 /// 5) HardCutLimitZ(Double_t value) - Similar to HardCutLimitX but applied to the
  59 ///      Z dimension.
  60 ///
  61 /// The analysis task also optionally generates histograms which show can be used
  62 /// to check the quality of the correlation between ESD and MC tracks.
  63 /// The histograms are generated in 2D for the pT and rapidity variable,
  64 /// and are called:
  65 /// 1) findableTracksHist - Filled with the MC tracks which should in principle
  66 ///      be findable by the reconstruction software. A findable track is defined
  67 ///      by the IsFindable method.
  68 /// 2) foundTracksHistMC - Filled with the MC information for each muon track that
  69 ///      was found in the ESD and for which we could match its corresponding MC track.
  70 /// 3) foundTracksHist - Filled just like the foundTracksHistMC histogram, but with
  71 ///      the momentum information taken from the ESD track and not the corresponding
  72 ///      MC track.
  73 /// 4) fakeTracksHist - Filled with all the ESD tracks for which we could not find
  74 ///      its matching MC track and assume it is a fake track.
  75 /// These histograms are filled by default, but can be disabled by using the
  76 /// GenerateHistograms(false) method call.
  77 ///
  78
  79 #include "AliAnalysisTaskLinkToMC.h"
  80 #include "AliESDEvent.h"
  81 #include "AliESDMuonTrack.h"
  82 #include "AliESDMuonCluster.h"
  83 #include "AliMCEvent.h"
  84 #include "AliMCParticle.h"
  85 #include "AliAODVertex.h"
  86 #include "AliAODEvent.h"
  87 #include "AliAODTrack.h"
  88 #include "AliTrackReference.h"
  89 #include "AliStack.h"
  90 #include "AliLog.h"
  91 #include "TObjArray.h"
  92 #include "TCanvas.h"
  93 #include "TMath.h"
  94 #include "TMap.h"
  95 #include "TList.h"
  96 #include "TH2D.h"
  97 #include "TStyle.h"
  98
  99 ClassImp(AliAnalysisTaskLinkToMC)
 100
 101
 102 void AliAnalysisTaskLinkToMC::HardCutLimitX(Double_t value)
 103 {
 104         /// Sets the absolute maximum distance in centimetres between the cluster
 105         /// and track reference that is allowed in the X direction.
 106
 107         if (value >= 0)
 108         {
 109                 fHardCutLimitX = value;
 110         }
 111         else
 112         {
 113                 AliWarning("Setting value for fHardCutLimitX which is negative."
 114                         " Will set it to a positive value."
 115                         );
 116                 fHardCutLimitX = TMath::Abs(value);
 117         }
 118 }
 119
 120
 121 void AliAnalysisTaskLinkToMC::HardCutLimitY(Double_t value)
 122 {
 123         /// Sets the absolute maximum distance in centimetres between the cluster
 124         /// and track reference that is allowed in the Y direction.
 125
 126         if (value >= 0)
 127         {
 128                 fHardCutLimitY = value;
 129         }
 130         else
 131         {
 132                 AliWarning("Setting value for fHardCutLimitY which is negative."
 133                         " Will set it to a positive value."
 134                         );
 135                 fHardCutLimitY = TMath::Abs(value);
 136         }
 137 }
 138
 139
 140 void AliAnalysisTaskLinkToMC::HardCutLimitZ(Double_t value)
 141 {
 142         /// Sets the absolute maximum distance in centimetres between the cluster
 143         /// and track reference that is allowed in the Z direction.
 144
 145         if (value >= 0)
 146         {
 147                 fHardCutLimitZ = value;
 148         }
 149         else
 150         {
 151                 AliWarning("Setting value for fHardCutLimitZ which is negative."
 152                         " Will set it to a positive value."
 153                         );
 154                 fHardCutLimitZ = TMath::Abs(value);
 155         }
 156 }
 157
 158
 159 Double_t AliAnalysisTaskLinkToMC::SigmaCut() const
 160 {
 161         /// Returns the maximum sigma value to allow for each X and Y dimension
 162         /// between clusters and track references.
 163
 164         return TMath::Sqrt(fSigmaCut2);
 165 }
 166
 167
 168 Bool_t AliAnalysisTaskLinkToMC::ChamberMustMatch(Int_t chamber) const
 169 {
 170         /// Returns the flag indicating if the ESD track cluster on the given
 171         /// chamber must match a corresponding track reference to be valid.
 172         /// \param chamber The chamber number must be in the range 1 to 14.
 173
 174         if (chamber >= 1 and chamber <= 14)
 175         {
 176                 return fChamberMustMatch[chamber-1];
 177         }
 178         else
 179         {
 180                 AliError(Form("The chamber number %d is invalid. Expected a value"
 181                         " in the range [1..14]. Will return a value of 'false'.",
 182                         chamber
 183                         ));
 184                 return false;
 185         }
 186 }
 187
 188
 189 void AliAnalysisTaskLinkToMC::ChamberMustMatch(Int_t chamber, Bool_t value)
 190 {
 191         /// Set the flag indicating if the ESD track cluster on the given chamber
 192         /// must match a corresponding track reference to be valid.
 193         /// \param chamber The chamber number must be in the range 1 to 14.
 194         /// \param value  The new value to set to.
 195
 196         if (chamber >= 1 and chamber <= 14)
 197         {
 198                 fChamberMustMatch[chamber-1] = value;
 199         }
 200         else
 201         {
 202                 AliError(Form("The chamber number %d is invalid. Expected a value"
 203                         " in the range [1..14]. Will not set any fChamberMustMatch flag.",
 204                         chamber
 205                         ));
 206         }
 207 }
 208
 209
 210 Bool_t AliAnalysisTaskLinkToMC::StationMustMatch(Int_t station) const
 211 {
 212         /// Returns the flag indicating if at least one ESD track cluster on the
 213         /// given station must match a corresponding track reference to be valid.
 214         /// \param station The station number must be in the range 1..7. Station
 215         ///     1 to 5 are the tracking stations while 6 is MT1 and 7 is MT2.
 216
 217         if (station >= 1 and station <= 7)
 218         {
 219                 return fStationMustMatch[station-1];
 220         }
 221         else
 222         {
 223                 AliError(Form("The station number %d is invalid. Expected a value"
 224                         " in the range [1..7]. Will return a value of 'false'.",
 225                         station
 226                         ));
 227                 return false;
 228         }
 229 }
 230
 231
 232 void AliAnalysisTaskLinkToMC::StationMustMatch(Int_t station, Bool_t value)
 233 {
 234         /// Sets the flag indicating if at least one ESD track cluster on the
 235         /// given station must match a corresponding track reference to be valid.
 236         /// \param station The station number must be in the range 1..7. Station
 237         ///     1 to 5 are the tracking stations while 6 is MT1 and 7 is MT2.
 238         /// \param value  The new value to set to.
 239
 240         if (station >= 1 and station <= 7)
 241         {
 242                 fStationMustMatch[station-1] = value;
 243         }
 244         else
 245         {
 246                 AliError(Form("The station number %d is invalid. Expected a value"
 247                         " in the range [1..7]. Will not set any fStationMustMatch flag.",
 248                         station
 249                         ));
 250         }
 251 }
 252
 253
 254 AliAnalysisTaskLinkToMC::AliAnalysisTaskLinkToMC() :
 255         AliAnalysisTaskSE(),
 256         fHistos(NULL),
 257         fFindableHist(NULL),
 258         fFoundHistMC(NULL),
 259         fFoundHist(NULL),
 260         fFakeHist(NULL),
 261         fHardCutLimitX(4.), // cm
 262         fHardCutLimitY(4.), // cm
 263         fHardCutLimitZ(4.), // cm
 264         fVarX(0.144*0.144), // cm^2
 265         fVarY(0.01*0.01), // cm^2
 266         fSigmaCut2(5*5),  // sigma^2
 267         fMinClusters(6),
 268         fMinClustersInSt12(0),  // no minimum requirement for stations 1 and 2.
 269         fMinClustersInSt45(3),  // 3 hits required on station 4 and 5 by default.
 270         fMakeHists(true),  // Generate histograms for monitoring by default
 271         fUseZCoordinate(false)  // Just use detElemId to match in the Z direction.
 272 {
 273         /// Default constructor.
 274         /// Sets input slot 0 to TChain input for ESDs, output slot 0 to TTree
 275         /// for AODs and output slot 1 to TList for the list of output histograms.
 276         /// The output histograms generated are to cross check the quality of the
 277         /// correlation and their generation can be turned off with
 278         /// AliAnalysisTaskLinkToMC::GenerateHistograms(false).
 279
 280         // Do not force any matching for particular chambers by default.
 281         for (Int_t i = 0; i < 14; i++)
 282         {
 283                 fChamberMustMatch[i] = false;
 284         }
 285         // Only force a minimum of 1 hit to be matched on stations 1, 2 and 3 by default.
 286         for (Int_t i = 0; i < 3; i++)
 287         {
 288                 fStationMustMatch[i] = true;
 289         }
 290         for (Int_t i = 3; i < 7; i++)
 291         {
 292                 fStationMustMatch[i] = false;
 293         }
 294
 295         DefineOutput(1, TList::Class());
 296 }
 297
 298
 299 AliAnalysisTaskLinkToMC::AliAnalysisTaskLinkToMC(const char* name) :
 300         AliAnalysisTaskSE(name),
 301         fHistos(NULL),
 302         fFindableHist(NULL),
 303         fFoundHistMC(NULL),
 304         fFoundHist(NULL),
 305         fFakeHist(NULL),
 306         fHardCutLimitX(4.), // cm
 307         fHardCutLimitY(4.), // cm
 308         fHardCutLimitZ(4.), // cm
 309         fVarX(0.144*0.144), // cm^2
 310         fVarY(0.01*0.01), // cm^2
 311         fSigmaCut2(5*5),  // sigma^2
 312         fMinClusters(6),
 313         fMinClustersInSt12(0),  // no minimum requirement for stations 1 and 2.
 314         fMinClustersInSt45(3),  // 3 hits required on station 4 and 5 by default.
 315         fMakeHists(true),  // Generate histograms for monitoring by default
 316         fUseZCoordinate(false)  // Just use detElemId to match in the Z direction.
 317 {
 318         /// Constructor with a task name.
 319         /// Sets input slot 0 to TChain input for ESDs, output slot 0 to TTree
 320         /// for AODs and output slot 1 to TList for the list of output histograms
 321         /// storing histograms generated to cross check the quility of the
 322         /// correlation.
 323         /// \param name  The name of the task.
 324
 325         // Do not force any matching for particular chambers by default.
 326         for (Int_t i = 0; i < 14; i++)
 327         {
 328                 fChamberMustMatch[i] = false;
 329         }
 330         // Only force a minimum of 1 hit to be matched on stations 1, 2 and 3 by default.
 331         for (Int_t i = 0; i < 3; i++)
 332         {
 333                 fStationMustMatch[i] = true;
 334         }
 335         for (Int_t i = 3; i < 7; i++)
 336         {
 337                 fStationMustMatch[i] = false;
 338         }
 339
 340         DefineOutput(1, TList::Class());
 341 }
 342
 343
 344 void AliAnalysisTaskLinkToMC::UserCreateOutputObjects()
 345 {
 346         /// Creates the output histograms containing findable tracks, found tracks
 347         /// and fake tracks. The binning range for all histograms is 30 bins along
 348         /// the pT dimension (X direction) from 0 to 20 GeV/c and 30 bins in the
 349         /// rapidity dimension (Y direction) form -4 to -2.4 units.
 350
 351         fHistos = new TList;
 352
 353         char name[1024];
 354         char title[1024];
 355
 356         Int_t nBinsX = 30;
 357         Double_t minX = 0.;
 358         Double_t maxX = 20.;
 359         Int_t nBinsY = 30;
 360         Double_t minY = -4.;
 361         Double_t maxY = -2.4;
 362
 363         sprintf(name, "findableTracksHist");
 364         sprintf(title, "Findable tracks");
 365         fFindableHist = new TH2D(name, title, nBinsX, minX, maxX, nBinsY, minY, maxY);
 366         fFindableHist->SetXTitle("p_{T} [GeV/c]");
 367         fFindableHist->SetYTitle("Rapidity (Y)");
 368         fHistos->Add(fFindableHist);
 369
 370         sprintf(name, "foundTracksHistMC");
 371         sprintf(title, "Found tracks (filled with Monte Carlo kinematic information)");
 372         fFoundHistMC = new TH2D(name, title, nBinsX, minX, maxX, nBinsY, minY, maxY);
 373         fFoundHistMC->SetXTitle("p_{T} [GeV/c]");
 374         fFoundHistMC->SetYTitle("Rapidity (Y)");
 375         fHistos->Add(fFoundHistMC);
 376
 377         sprintf(name, "foundTracksHist");
 378         sprintf(title, "Found tracks (filled with reconstructed kinematic information)");
 379         fFoundHist = new TH2D(name, title, nBinsX, minX, maxX, nBinsY, minY, maxY);
 380         fFoundHist->SetXTitle("p_{T} [GeV/c]");
 381         fFoundHist->SetYTitle("Rapidity (Y)");
 382         fHistos->Add(fFoundHist);
 383
 384         sprintf(name, "fakeTracksHist");
 385         sprintf(title, "Fake tracks");
 386         fFakeHist = new TH2D(name, title, nBinsX, minX, maxX, nBinsY, minY, maxY);
 387         fFakeHist->SetXTitle("p_{T} [GeV/c]");
 388         fFakeHist->SetYTitle("Rapidity (Y)");
 389         fHistos->Add(fFakeHist);
 390 }
 391
 392
 393 void AliAnalysisTaskLinkToMC::UserExec(Option_t* /*option*/)
 394 {
 395         /// This is the top level method that performs the work of the task.
 396         /// It will fill the histograms for quality checking with the Monte
 397         /// Carlo (MC) information if they are to be generated.
 398         /// The input event is then processed which must be of type AliESDEvent.
 399         /// For all the ESD muon tracks we try to find the corresponding MC
 400         /// track and build a TMap for the mapping.
 401         /// If the AOD output event is available then it is filled with AOD
 402         /// tracks which copy the information form the ESD track and the AOD
 403         /// track label is filled with the corresponding MC track label.
 404         /// If a corresponding MC track could not be found then the AOD track
 405         /// label is filled with -1.
 406         /// Finally the checking histograms are further filled with the ESD
 407         /// information based on the ESD to MC track map.
 408
 409         if (fMakeHists and
 410             (fFindableHist == NULL or fFoundHistMC == NULL or fFoundHist == NULL or fFakeHist == NULL)
 411            )
 412         {
 413                 AliError("Histograms have not been created.");
 414                 return;
 415         }
 416
 417         if (fMakeHists and MCEvent() != NULL) FillHistsFromMC();
 418
 419         if (InputEvent() != NULL)
 420         {
 421                 if (InputEvent()->IsA() == AliESDEvent::Class())
 422                 {
 423                         if (MCEvent() != NULL)
 424                         {
 425                                 TMap links;
 426                                 FindLinksToMC(links);
 427                                 if (AODEvent() != NULL) CreateAODTracks(links);
 428                                 if (fMakeHists) FillHistsFromLinks(links);
 429                         }
 430                         else
 431                         {
 432                                 AliError("To process the ESD event we must have the Monte Carlo data also.");
 433                         }
 434                 }
 435                 else
 436                 {
 437                         AliError(Form("The input event is of type \"%s\". Do not know how to handle it so it will be skipped.",
 438                                         InputEvent()->ClassName()
 439                                 ));
 440                 }
 441         }
 442
 443         if (fMakeHists) PostData(1, fHistos);
 444 }
 445
 446
 447 void AliAnalysisTaskLinkToMC::Terminate(Option_t* /*option*/)
 448 {
 449         /// The terminate method will draw the histograms filled by this task if
 450         /// they were filled by setting appropriate flag with the
 451         /// GenerateHistograms(true) method, which is the default.
 452
 453         if (fMakeHists)
 454         {
 455                 gStyle->SetPalette(1);
 456                 new TCanvas;
 457                 fFindableHist->DrawCopy("colz");
 458                 new TCanvas;
 459                 fFoundHistMC->DrawCopy("colz");
 460                 new TCanvas;
 461                 fFoundHist->DrawCopy("colz");
 462                 new TCanvas;
 463                 fFakeHist->DrawCopy("colz");
 464
 465                 AliInfo(Form("Findable tracks = %d", Int_t(fFindableHist->GetEntries()) ));
 466                 AliInfo(Form("Found tracks    = %d", Int_t(fFoundHistMC->GetEntries()) ));
 467                 AliInfo(Form("Fake tracks     = %d", Int_t(fFakeHist->GetEntries()) ));
 468         }
 469 }
 470
 471
 472 bool AliAnalysisTaskLinkToMC::IsFindable(AliMCParticle* track) const
 473 {
 474         /// This method is used to identify if a Monte Carlo (MC) track is in principle
 475         /// findable in the muon spectrometer.
 476         /// This means the MC track must be a muon particle and leave at least one track
 477         /// reference point on either chamber in each tracking station.
 478         /// \param track  The MC track to check.
 479         /// \returns  true if the track is findable by the muon spectrometer
 480         ///    and false otherwise.
 481         /// \note  The definition of track find-ability does not affect the correlation
 482         ///    between ESD and MC tracks nor the generation of the AOD output tracks,
 483         ///    only the histogram of findable tracks is effected if it is generated at all.
 484
 485         // Select only muons.
 486         if (TMath::Abs(track->Particle()->GetPdgCode()) != 13) return false;
 487
 488         // Build hit mask for chambers.
 489         bool hitOnCh[14];
 490         for (Int_t i = 0; i < 14; i++)
 491         {
 492                 hitOnCh[i] = false;
 493         }
 494         for (Int_t i = 0; i < track->GetNumberOfTrackReferences(); i++)
 495         {
 496                 AliTrackReference* ref = track->GetTrackReference(i);
 497                 if (ref->DetectorId() != AliTrackReference::kMUON) continue;
 498                 Int_t chamberId = ref->UserId() / 100 - 1;
 499                 if (0 <= chamberId and chamberId < 14)
 500                 {
 501                         hitOnCh[chamberId] = true;
 502                 }
 503         }
 504
 505         // Check that enough hits are available on all tracking chambers.
 506         bool hitOnSt1 = hitOnCh[0] or hitOnCh[1];
 507         bool hitOnSt2 = hitOnCh[2] or hitOnCh[3];
 508         bool hitOnSt3 = hitOnCh[4] or hitOnCh[5];
 509         bool hitOnSt4 = hitOnCh[6] or hitOnCh[7];
 510         bool hitOnSt5 = hitOnCh[8] or hitOnCh[9];
 511
 512         return (hitOnSt1 and hitOnSt2 and hitOnSt3 and hitOnSt4 and hitOnSt5);
 513 }
 514
 515
 516 void AliAnalysisTaskLinkToMC::FillHistsFromMC()
 517 {
 518         /// Fills the histogram containing findable tracks from the MC event information.
 519         /// Only MC tracks for which the IsFindable method returns true are added
 520         /// to the histogram.
 521
 522         for (Int_t i = 0; i < MCEvent()->GetNumberOfTracks(); i++)
 523         {
 524                 AliMCParticle* mcTrack = (AliMCParticle*) MCEvent()->GetTrack(i);
 525
 526                 // Select only reconstructible tracks to fill into findable hist.
 527                 if (IsFindable(mcTrack))
 528                 {
 529                         fFindableHist->Fill(mcTrack->Pt(), mcTrack->Y());
 530                 }
 531         }
 532 }
 533
 534
 535 void AliAnalysisTaskLinkToMC::FillHistsFromLinks(TMap& links)
 536 {
 537         /// Fills the histograms containing found tracks and fake tracks.
 538         /// Found tracks are all those that have a ESD muon track in the ESD event
 539         /// and for which a corresponding MC track could be found.
 540         /// Fake tracks are ESD muon tracks for which no corresponding MC track
 541         /// could be matched so it is considered a fake track.
 542         /// \param links  This contains the mapping between ESD muon tracks and
 543         ///    MC tracks. ESD tracks for which no MC track could be matched should
 544         ///    have the "value" of the key/value pair in the mapping set to NULL.
 545
 546         TIter itmap(links.GetTable());
 547         TPair* pair = NULL;
 548         while ( (pair = static_cast<TPair*>(itmap())) != NULL )
 549         {
 550                 if (pair->Value() != NULL)
 551                 {
 552                         AliESDMuonTrack* esdTrack = static_cast<AliESDMuonTrack*>(pair->Key());
 553                         fFoundHist->Fill(esdTrack->Pt(), esdTrack->Y());
 554
 555                         AliMCParticle* mcTrack = static_cast<AliMCParticle*>(pair->Value());
 556                         fFoundHistMC->Fill(mcTrack->Pt(), mcTrack->Y());
 557                 }
 558                 else
 559                 {
 560                         AliESDMuonTrack* esdTrack = static_cast<AliESDMuonTrack*>(pair->Key());
 561                         fFakeHist->Fill(esdTrack->Pt(), esdTrack->Y());
 562                 }
 563         }
 564 }
 565
 566
 567 void AliAnalysisTaskLinkToMC::CreateAODTracks(TMap& links)
 568 {
 569         /// This code is copied from AliAnalysisTaskESDMuonFilter::ConvertESDtoAOD with
 570         /// modifications to fill the MC track label using the ESD to MC track mapping.
 571         /// \param links  This is the mapping between ESD tracks and MC tracks.
 572
 573         // ESD Muon Filter analysis task executed for each event
 574         AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
 575
 576         // Define arrays for muons
 577         Double_t pos[3];
 578         Double_t p[3];
 579         Double_t pid[10];
 580
 581         // has to be changed once the muon pid is provided by the ESD
 582         for (Int_t i = 0; i < 10; pid[i++] = 0.) {}
 583         pid[AliAODTrack::kMuon]=1.;
 584
 585         AliAODHeader* header = AODEvent()->GetHeader();
 586         AliAODTrack *aodTrack = 0x0;
 587         AliESDMuonTrack *esdMuTrack = 0x0;
 588
 589         // Access to the AOD container of tracks
 590         TClonesArray &tracks = *(AODEvent()->GetTracks());
 591         Int_t jTracks = tracks.GetEntriesFast();
 592
 593         // Read primary vertex from AOD event
 594         AliAODVertex *primary = AODEvent()->GetPrimaryVertex();
 595         if (primary != NULL) primary->Print();
 596
 597         // Loop on muon tracks to fill the AOD track branch
 598         Int_t nMuTracks = esd->GetNumberOfMuonTracks();
 599         printf("Number of Muon Tracks=%d\n",nMuTracks);
 600
 601         // Update number of positive and negative tracks from AOD event (M.G.)
 602         Int_t nPosTracks = header->GetRefMultiplicityPos();
 603         Int_t nNegTracks = header->GetRefMultiplicityNeg();
 604
 605         for (Int_t nMuTrack = 0; nMuTrack < nMuTracks; ++nMuTrack)
 606         {
 607                 esdMuTrack = esd->GetMuonTrack(nMuTrack);
 608
 609                 //if (!esdMuTrack->ContainTrackerData()) continue;
 610
 611                 UInt_t selectInfo = 0;
 612
 613                 p[0] = esdMuTrack->Px();
 614                 p[1] = esdMuTrack->Py();
 615                 p[2] = esdMuTrack->Pz();
 616
 617                 pos[0] = esdMuTrack->GetNonBendingCoor();
 618                 pos[1] = esdMuTrack->GetBendingCoor();
 619                 pos[2] = esdMuTrack->GetZ();
 620
 621                 Int_t label = -1;
 622                 TPair* pair = static_cast<TPair*>( links.FindObject(esdMuTrack) );
 623                 if (pair != NULL and pair->Value() != NULL)
 624                 {
 625                         AliMCParticle* mcTrack = static_cast<AliMCParticle*>(pair->Value());
 626                         label = mcTrack->GetLabel();
 627                         if (label == -1)
 628                         {
 629                                 AliWarning("The MC track was set to -1.");
 630                         }
 631                 }
 632
 633                 aodTrack = new(tracks[jTracks++]) AliAODTrack(
 634                                 esdMuTrack->GetUniqueID(), // ID
 635                                 label, // label
 636                                 p, // momentum
 637                                 kTRUE, // cartesian coordinate system
 638                                 pos, // position
 639                                 kFALSE, // isDCA
 640                                 0x0, // covariance matrix
 641                                 esdMuTrack->Charge(), // charge
 642                                 0, // ITSClusterMap
 643                                 pid, // pid
 644                                 primary, // primary vertex
 645                                 kFALSE, // used for vertex fit?
 646                                 kFALSE, // used for primary vertex fit?
 647                                 AliAODTrack::kPrimary,// track type
 648                                 selectInfo
 649                         );
 650
 651                 aodTrack->SetXYAtDCA(esdMuTrack->GetNonBendingCoorAtDCA(), esdMuTrack->GetBendingCoorAtDCA());
 652                 aodTrack->SetPxPyPzAtDCA(esdMuTrack->PxAtDCA(), esdMuTrack->PyAtDCA(), esdMuTrack->PzAtDCA());
 653                 aodTrack->ConvertAliPIDtoAODPID();
 654                 aodTrack->SetChi2perNDF(esdMuTrack->GetChi2() / (2.*esdMuTrack->GetNHit() - 5.));
 655                 aodTrack->SetChi2MatchTrigger(esdMuTrack->GetChi2MatchTrigger());
 656                 aodTrack->SetHitsPatternInTrigCh(esdMuTrack->GetHitsPatternInTrigCh());
 657                 aodTrack->SetMuonClusterMap(esdMuTrack->GetMuonClusterMap());
 658                 aodTrack->SetMatchTrigger(esdMuTrack->GetMatchTrigger());
 659
 660                 if (primary != NULL) primary->AddDaughter(aodTrack);
 661
 662                 if (esdMuTrack->Charge() > 0) nPosTracks++;
 663                 else nNegTracks++;
 664         }
 665
 666         header->SetRefMultiplicity(jTracks);
 667         header->SetRefMultiplicityPos(nPosTracks);
 668         header->SetRefMultiplicityNeg(nNegTracks);
 669
 670         // Note: we do not have to call PostData for the AOD because the AliAnalysisTaskSE
 671         // base class already does so just after the UserExec.
 672 }
 673
 674
 675 void AliAnalysisTaskLinkToMC::FindLinksToMC(TMap& links)
 676 {
 677         /// Finds all MC tracks that correspond to their ESD tracks and creates a
 678         /// mapping between the ESD tracks and MC tracks which is filled into the
 679         /// 'links' TMap object.
 680         /// \param links  This will be filled with the ESD and MC track pairs which
 681         ///    have been matched. Thus, for each ESD track which is set as the TMap key,
 682         ///    the corresponding MC track is identified by the TMap value. If no
 683         ///    MC track could be matched then the MC track value is set to NULL.
 684
 685         // Algorithm description for connecting reconstructed ESD muon tracks to MC:
 686         //
 687         // Build ESD track list
 688         // Build MC track list
 689         // while ESD track list not empty and MC track list not empty:
 690         //   Find ESD/MC track pair with best compatibility.
 691         //   if no pair is compatible then break loop.
 692         //   Store found pair in return list.
 693         //   Remove pair from ESD and MC track list.
 694         // Mark remaining tracks in ESD track list as fakes.
 695         //
 696         // To find ESD/MC track pair with best compatibility:
 697         //   bestQuality := 1e100
 698         //   bestESDTrack := nil
 699         //   bestMCTrack := nil
 700         //   for each ESD track:
 701         //     for each MC track:
 702         //       chi^2 / numOfClustersMatched := calculate compatibility
 703         //       if tracks are compatible and chi^2 / numOfClustersMatched < bestQuality then:
 704         //         bestQuality := chi^2 / noOfClustersMatched
 705         //         bestESDTrack := current ESD track
 706         //         bestMCTrack := current MC track
 707         //
 708         // To calculate compatibility (ESD track, MC track):
 709         //   chi^2 := 1e100
 710         //   numberOfClustersAssociated := 0
 711         //   for each cluster in ESD track:
 712         //     find nearest 2 track refs from MC track
 713         //     if 2 track refs found then find average of them into MChit
 714         //     else MChit := nearest ref track.
 715         //     k := calculate chi^2 distance between ESD cluster and MChit
 716         //     if k <= maxValue then:
 717         //       chi^2 := chi^2 + k
 718         //       numberOfClustersAssociated := numberOfClustersAssociated + 1
 719         //     else:
 720         //       chi^2 := chi^2 + maxAllowedChi2
 721         //   if numberOfClustersAssociated >= minGoodClusters then:
 722         //     return (compatible, chi^2)
 723         //   else:
 724         //     return (not compatible, 1e100)
 725
 726         AliESDEvent* esd = static_cast<AliESDEvent*>(InputEvent());
 727
 728         TObjArray esdTracks;
 729         TObjArray mcTracks;
 730         for (Int_t i = 0; i < esd->GetNumberOfMuonTracks(); i++)
 731         {
 732                 esdTracks.Add(esd->GetMuonTrack(i));
 733         }
 734         for (Int_t i = 0; i < MCEvent()->GetNumberOfTracks(); i++)
 735         {
 736                 mcTracks.Add(MCEvent()->GetTrack(i));
 737         }
 738
 739         while (not esdTracks.IsEmpty() and not mcTracks.IsEmpty())
 740         {
 741                 Double_t bestQuality = 1e100;
 742                 Int_t bestNoOfClusters = 0;
 743                 AliESDMuonTrack* bestEsdTrack = NULL;
 744                 AliMCParticle* bestMcTrack = NULL;
 745                 TIter itesd(&esdTracks);
 746                 TIter itmc(&mcTracks);
 747                 AliESDMuonTrack* esdTrack = NULL;
 748                 AliMCParticle* mcTrack = NULL;
 749
 750                 // Find the ESD and MC track pair that matches the best.
 751                 while ( (esdTrack = static_cast<AliESDMuonTrack*>(itesd())) != NULL )
 752                 {
 753                         while ( (mcTrack = static_cast<AliMCParticle*>(itmc())) != NULL )
 754                         {
 755                                 Double_t chi2 = 1e100;
 756                                 bool tracksCompatible = false;
 757                                 Int_t noOfClusters = 0;
 758                                 CalculateTrackCompatibility(esdTrack, mcTrack, tracksCompatible, chi2, noOfClusters);
 759                                 Double_t quality = noOfClusters > 0 ? chi2 / Double_t(noOfClusters) : 1e100;
 760                                 if (tracksCompatible and quality < bestQuality and noOfClusters >= bestNoOfClusters)
 761                                 {
 762                                         bestQuality = quality;
 763                                         bestNoOfClusters = noOfClusters;
 764                                         bestEsdTrack = esdTrack;
 765                                         bestMcTrack = mcTrack;
 766                                 }
 767                         }
 768                         itmc.Reset();
 769                 }
 770
 771                 // If no track pair could be matched then we are done because we
 772                 // will not be able to match anything more next time.
 773                 if (bestMcTrack == NULL) break;
 774
 775                 links.Add(bestEsdTrack, bestMcTrack);
 776
 777                 esdTracks.Remove(bestEsdTrack);
 778                 mcTracks.Remove(bestMcTrack);
 779         }
 780
 781         // Add all remaining ESD tracks which could not be matched to any MC track
 782         // to the mapping but with NULL for the MC track.
 783         TIter itesd(&esdTracks);
 784         AliESDMuonTrack* esdTrack = NULL;
 785         while ( (esdTrack = static_cast<AliESDMuonTrack*>(itesd())) != NULL )
 786         {
 787                 links.Add(esdTrack, NULL);
 788         }
 789 }
 790
 791
 792 void AliAnalysisTaskLinkToMC::CalculateTrackCompatibility(
 793                 AliESDMuonTrack* esdTrack, AliMCParticle* mcTrack,
 794                 bool& tracksCompatible, Double_t& chi2, Int_t& noOfClusters
 795         )
 796 {
 797         /// Calculates the compatibility between a ESD and MC track pair.
 798         /// For each cluster of the ESD track, two track references are found
 799         /// which are spatially closest to it. If the distance between the cluster
 800         /// and track reference is larger than the HardCutLimit*() parameters then
 801         /// the track reference is not used. The average coordinate is calculated
 802         /// from the two track references found and the chi-square between the cluster
 803         /// and average coordinate is calculated. If only one track reference is
 804         /// found then it is used for the calculation without taking an average.
 805         /// If no track references are found then the cluster is not matched and
 806         /// we continue to the next cluster.
 807         /// If the chi-squared value for the cluster / track reference pair is
 808         /// larger than that allowed by the SigmaCut() parameter then the cluster
 809         /// is also considered not to be matched.
 810         /// The tracks are compatible if the ESD track has a minimum number of
 811         /// clusters matched with MC track references from the MC track.
 812         /// [in] \param esdTrack  The ESD track we are trying to match.
 813         /// [in] \param mcTrack  The MC track we are trying to match the ESD track to.
 814         /// [out] \param tracksCompatible  Filled with true if the given esdTrack
 815         ///     and mcTrack are compatible.
 816         /// [out] \param chi2  Filled with the total chi-squared calculated for all
 817         ///     matched clusters.
 818         /// [out] \param noOfClusters  Filled with the number of clusters that were
 819         ///     matched to at least one track reference.
 820
 821         chi2 = 0;
 822         noOfClusters = 0;
 823         tracksCompatible = false;
 824         Int_t noOfClustersSt12 = 0;
 825         Int_t noOfClustersSt45 = 0;
 826         bool chamberMatched[14] = {
 827                         false, false, false, false, false, false, false,
 828                         false, false, false, false, false, false, false
 829                 };
 830         bool stationMatched[7] = {
 831                         false, false, false, false, false, false, false
 832                 };
 833
 834         for (Int_t i = 0; i < esdTrack->GetNClusters(); i++)
 835         {
 836                 AliESDMuonCluster* cluster = static_cast<AliESDMuonCluster*>( esdTrack->GetClusters()[i] );
 837                 Double_t varX = cluster->GetErrX2();
 838                 Double_t varY = cluster->GetErrY2();
 839                 // If the variance is zero then use the default one or just 1
 840                 // if the default is also not set.
 841                 if (varX == 0) varX = fVarX != 0 ? fVarX : 1.;
 842                 if (varY == 0) varY = fVarY != 0 ? fVarY : 1.;
 843
 844                 Double_t chi2A = 1e100;
 845                 AliTrackReference* refA = NULL;
 846                 Double_t chi2B = 1e100;
 847                 AliTrackReference* refB = NULL;
 848
 849                 for (Int_t j = 0; j < mcTrack->GetNumberOfTrackReferences(); j++)
 850                 {
 851                         AliTrackReference* ref = mcTrack->GetTrackReference(j);
 852
 853                         if (ref->DetectorId() != AliTrackReference::kMUON) continue;
 854                         if (ref->UserId() != cluster->GetDetElemId()) continue;
 855
 856                         Double_t dx = ref->X() - cluster->GetX();
 857                         if (TMath::Abs(dx) > fHardCutLimitX) continue;
 858                         Double_t dy = ref->Y() - cluster->GetY();
 859                         if (TMath::Abs(dy) > fHardCutLimitY) continue;
 860                         if (fUseZCoordinate)
 861                         {
 862                                 Double_t dz = ref->Z() - cluster->GetZ();
 863                                 if (TMath::Abs(dz) > fHardCutLimitZ) continue;
 864                         }
 865
 866                         Double_t chi2sum = dx*dx / varX + dy*dy / varY;
 867                         if (chi2sum < chi2A)
 868                         {
 869                                 // Copy track reference A to B and set new value for A
 870                                 chi2B = chi2A;
 871                                 refB = refA;
 872                                 chi2A = chi2sum;
 873                                 refA = ref;
 874                         }
 875                         else if (chi2sum < chi2B)
 876                         {
 877                                 // Track reference A still has smaller chi2 so just replace B.
 878                                 chi2B = chi2sum;
 879                                 refB = ref;
 880                         }
 881                 }
 882
 883                 Double_t x = 0, y = 0;
 884                 if (refA != NULL and refB != NULL)
 885                 {
 886                         // Average the track references if 2 were found.
 887                         x = (refA->X() + refB->X()) * 0.5;
 888                         y = (refA->Y() + refB->Y()) * 0.5;
 889                 }
 890                 else if (refA != NULL)
 891                 {
 892                         x = refA->X();
 893                         y = refA->Y();
 894                 }
 895                 else
 896                 {
 897                         continue;  // no track reference hit found.
 898                 }
 899
 900                 Double_t dx = x - cluster->GetX();
 901                 Double_t dy = y - cluster->GetY();
 902
 903                 Double_t chi2x = dx*dx / varX;
 904                 if (chi2x > fSigmaCut2) continue;
 905                 Double_t chi2y = dy*dy / varY;
 906                 if (chi2y > fSigmaCut2) continue;
 907
 908                 chi2 += chi2x + chi2y;
 909                 noOfClusters++;
 910                 if (cluster->GetChamberId() >= 0 and cluster->GetChamberId() < 4)
 911                 {
 912                         noOfClustersSt12++;
 913                 }
 914                 if (cluster->GetChamberId() >= 6 and cluster->GetChamberId() < 10)
 915                 {
 916                         noOfClustersSt45++;
 917                 }
 918
 919                 if (0 <= cluster->GetChamberId() and cluster->GetChamberId() < 14)
 920                 {
 921                         chamberMatched[cluster->GetChamberId()] = true;
 922                         stationMatched[cluster->GetChamberId()/2] = true;
 923                 }
 924         }
 925
 926         // Need to check that all the chambers and stations that had to match we
 927         // actually found matching track references.
 928
 929         bool forcedChambersMatched = true;
 930         for (Int_t i = 0; i < 14; i++)
 931         {
 932                 if (fChamberMustMatch[i] and not chamberMatched[i])
 933                 {
 934                         forcedChambersMatched = false;
 935                 }
 936         }
 937         bool forcedStationsMatched = true;
 938         for (Int_t i = 0; i < 7; i++)
 939         {
 940                 if (fStationMustMatch[i] and not stationMatched[i])
 941                 {
 942                         forcedStationsMatched = false;
 943                 }
 944         }
 945
 946         if (noOfClusters >= fMinClusters and
 947             noOfClustersSt12 >= fMinClustersInSt12 and
 948             noOfClustersSt45 >= fMinClustersInSt45 and
 949             forcedChambersMatched and forcedStationsMatched
 950            )
 951         {
 952                 tracksCompatible = true;
 953         }
 954 }
 955