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