1 /**************************************************************************
2 * This file is property of and copyright by the ALICE HLT Project *
3 * All rights reserved. *
6 * Artur Szostak <artursz@iafrica.com> *
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 **************************************************************************/
20 /// @file AliAnalysisTaskLinkToMC.cxx
21 /// @author Artur Szostak <artursz@iafrica.com>
23 /// @brief Implementation of the AliAnalysisTaskLinkToMC task.
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
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
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
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
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,
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
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.
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"
93 #include "TObjArray.h"
101 ClassImp(AliAnalysisTaskLinkToMC)
104 void AliAnalysisTaskLinkToMC::HardCutLimitX(Double_t value)
106 /// Sets the absolute maximum distance in centimetres between the cluster
107 /// and track reference that is allowed in the X direction.
111 fHardCutLimitX = value;
115 AliWarning("Setting value for fHardCutLimitX which is negative."
116 " Will set it to a positive value."
118 fHardCutLimitX = TMath::Abs(value);
123 void AliAnalysisTaskLinkToMC::HardCutLimitY(Double_t value)
125 /// Sets the absolute maximum distance in centimetres between the cluster
126 /// and track reference that is allowed in the Y direction.
130 fHardCutLimitY = value;
134 AliWarning("Setting value for fHardCutLimitY which is negative."
135 " Will set it to a positive value."
137 fHardCutLimitY = TMath::Abs(value);
142 void AliAnalysisTaskLinkToMC::HardCutLimitZ(Double_t value)
144 /// Sets the absolute maximum distance in centimetres between the cluster
145 /// and track reference that is allowed in the Z direction.
149 fHardCutLimitZ = value;
153 AliWarning("Setting value for fHardCutLimitZ which is negative."
154 " Will set it to a positive value."
156 fHardCutLimitZ = TMath::Abs(value);
161 Double_t AliAnalysisTaskLinkToMC::SigmaCut() const
163 /// Returns the maximum sigma value to allow for each X and Y dimension
164 /// between clusters and track references.
166 return TMath::Sqrt(fSigmaCut2);
170 Bool_t AliAnalysisTaskLinkToMC::ChamberMustMatch(Int_t chamber) const
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.
176 if (chamber >= 1 and chamber <= 14)
178 return fChamberMustMatch[chamber-1];
182 AliError(Form("The chamber number %d is invalid. Expected a value"
183 " in the range [1..14]. Will return a value of 'false'.",
191 void AliAnalysisTaskLinkToMC::ChamberMustMatch(Int_t chamber, Bool_t value)
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.
198 if (chamber >= 1 and chamber <= 14)
200 fChamberMustMatch[chamber-1] = value;
204 AliError(Form("The chamber number %d is invalid. Expected a value"
205 " in the range [1..14]. Will not set any fChamberMustMatch flag.",
212 Bool_t AliAnalysisTaskLinkToMC::StationMustMatch(Int_t station) const
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.
219 if (station >= 1 and station <= 7)
221 return fStationMustMatch[station-1];
225 AliError(Form("The station number %d is invalid. Expected a value"
226 " in the range [1..7]. Will return a value of 'false'.",
234 void AliAnalysisTaskLinkToMC::StationMustMatch(Int_t station, Bool_t value)
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.
242 if (station >= 1 and station <= 7)
244 fStationMustMatch[station-1] = value;
248 AliError(Form("The station number %d is invalid. Expected a value"
249 " in the range [1..7]. Will not set any fStationMustMatch flag.",
256 AliAnalysisTaskLinkToMC::AliAnalysisTaskLinkToMC() :
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
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.
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).
282 // Do not force any matching for particular chambers by default.
283 for (Int_t i = 0; i < 14; i++)
285 fChamberMustMatch[i] = false;
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++)
290 fStationMustMatch[i] = true;
292 for (Int_t i = 3; i < 7; i++)
294 fStationMustMatch[i] = false;
297 DefineOutput(1, TList::Class());
301 AliAnalysisTaskLinkToMC::AliAnalysisTaskLinkToMC(const char* name) :
302 AliAnalysisTaskSE(name),
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
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.
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
325 /// \param name The name of the task.
327 // Do not force any matching for particular chambers by default.
328 for (Int_t i = 0; i < 14; i++)
330 fChamberMustMatch[i] = false;
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++)
335 fStationMustMatch[i] = true;
337 for (Int_t i = 3; i < 7; i++)
339 fStationMustMatch[i] = false;
342 DefineOutput(1, TList::Class());
346 void AliAnalysisTaskLinkToMC::UserCreateOutputObjects()
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.
363 Double_t maxY = -2.4;
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);
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);
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);
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);
395 void AliAnalysisTaskLinkToMC::UserExec(Option_t* /*option*/)
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.
412 (fFindableHist == NULL or fFoundHistMC == NULL or fFoundHist == NULL or fFakeHist == NULL)
415 AliError("Histograms have not been created.");
419 if (fMakeHists and MCEvent() != NULL) FillHistsFromMC();
421 if (InputEvent() != NULL)
423 if (InputEvent()->IsA() == AliESDEvent::Class())
425 if (MCEvent() != NULL)
428 FindLinksToMC(links);
429 if (AODEvent() != NULL) CreateAODTracks(links);
430 if (fMakeHists) FillHistsFromLinks(links);
434 AliError("To process the ESD event we must have the Monte Carlo data also.");
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()
445 if (fMakeHists) PostData(1, fHistos);
449 void AliAnalysisTaskLinkToMC::Terminate(Option_t* /*option*/)
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.
457 gStyle->SetPalette(1);
459 fFindableHist->DrawCopy("colz");
461 fFoundHistMC->DrawCopy("colz");
463 fFoundHist->DrawCopy("colz");
465 fFakeHist->DrawCopy("colz");
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()) ));
474 bool AliAnalysisTaskLinkToMC::IsFindable(AliMCParticle* track) const
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.
487 // Select only muons.
488 if (TMath::Abs(track->Particle()->GetPdgCode()) != 13) return false;
490 // Build hit mask for chambers.
492 for (Int_t i = 0; i < 14; i++)
496 for (Int_t i = 0; i < track->GetNumberOfTrackReferences(); i++)
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)
503 hitOnCh[chamberId] = true;
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];
514 return (hitOnSt1 and hitOnSt2 and hitOnSt3 and hitOnSt4 and hitOnSt5);
518 void AliAnalysisTaskLinkToMC::FillHistsFromMC()
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.
524 for (Int_t i = 0; i < MCEvent()->GetNumberOfTracks(); i++)
526 AliMCParticle* mcTrack = (AliMCParticle*) MCEvent()->GetTrack(i);
528 // Select only reconstructible tracks to fill into findable hist.
529 if (IsFindable(mcTrack))
531 fFindableHist->Fill(mcTrack->Pt(), mcTrack->Y());
537 void AliAnalysisTaskLinkToMC::FillHistsFromLinks(TMap& links)
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.
548 TIter itmap(links.GetTable());
550 while ( (pair = static_cast<TPair*>(itmap())) != NULL )
552 if (pair->Value() != NULL)
554 AliESDMuonTrack* esdTrack = static_cast<AliESDMuonTrack*>(pair->Key());
555 fFoundHist->Fill(esdTrack->Pt(), esdTrack->Y());
557 AliMCParticle* mcTrack = static_cast<AliMCParticle*>(pair->Value());
558 fFoundHistMC->Fill(mcTrack->Pt(), mcTrack->Y());
562 AliESDMuonTrack* esdTrack = static_cast<AliESDMuonTrack*>(pair->Key());
563 fFakeHist->Fill(esdTrack->Pt(), esdTrack->Y());
569 void AliAnalysisTaskLinkToMC::CreateAODTracks(TMap& links)
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.
575 // ESD Muon Filter analysis task executed for each event
576 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
579 AliError("Cannot get input event");
583 // Define arrays for muons
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.;
592 AliAODHeader* header = AODEvent()->GetHeader();
593 AliAODTrack *aodTrack = 0x0;
594 AliESDMuonTrack *esdMuTrack = 0x0;
596 // Access to the AOD container of tracks
597 TClonesArray &tracks = *(AODEvent()->GetTracks());
598 Int_t jTracks = tracks.GetEntriesFast();
600 // Read primary vertex from AOD event
601 AliAODVertex *primary = AODEvent()->GetPrimaryVertex();
602 if (primary != NULL) primary->Print();
604 // Loop on muon tracks to fill the AOD track branch
605 Int_t nMuTracks = esd->GetNumberOfMuonTracks();
606 printf("Number of Muon Tracks=%d\n",nMuTracks);
608 // Update number of positive and negative tracks from AOD event (M.G.)
609 Int_t nPosTracks = header->GetRefMultiplicityPos();
610 Int_t nNegTracks = header->GetRefMultiplicityNeg();
612 for (Int_t nMuTrack = 0; nMuTrack < nMuTracks; ++nMuTrack)
614 esdMuTrack = esd->GetMuonTrack(nMuTrack);
616 //if (!esdMuTrack->ContainTrackerData()) continue;
618 UInt_t selectInfo = 0;
620 p[0] = esdMuTrack->Px();
621 p[1] = esdMuTrack->Py();
622 p[2] = esdMuTrack->Pz();
624 pos[0] = esdMuTrack->GetNonBendingCoor();
625 pos[1] = esdMuTrack->GetBendingCoor();
626 pos[2] = esdMuTrack->GetZ();
629 TPair* pair = static_cast<TPair*>( links.FindObject(esdMuTrack) );
630 if (pair != NULL and pair->Value() != NULL)
632 AliMCParticle* mcTrack = static_cast<AliMCParticle*>(pair->Value());
633 label = mcTrack->GetLabel();
636 AliWarning("The MC track was set to -1.");
640 aodTrack = new(tracks[jTracks++]) AliAODTrack(
641 esdMuTrack->GetUniqueID(), // ID
644 kTRUE, // cartesian coordinate system
647 0x0, // covariance matrix
648 esdMuTrack->Charge(), // charge
651 primary, // primary vertex
652 kFALSE, // used for vertex fit?
653 kFALSE, // used for primary vertex fit?
654 AliAODTrack::kPrimary,// track type
658 aodTrack->SetXYAtDCA(esdMuTrack->GetNonBendingCoorAtDCA(), esdMuTrack->GetBendingCoorAtDCA());
659 aodTrack->SetPxPyPzAtDCA(esdMuTrack->PxAtDCA(), esdMuTrack->PyAtDCA(), esdMuTrack->PzAtDCA());
660 aodTrack->ConvertAliPIDtoAODPID();
661 aodTrack->SetChi2perNDF(esdMuTrack->GetChi2() / (2.*esdMuTrack->GetNHit() - 5.));
662 aodTrack->SetChi2MatchTrigger(esdMuTrack->GetChi2MatchTrigger());
663 aodTrack->SetHitsPatternInTrigCh(esdMuTrack->GetHitsPatternInTrigCh());
664 aodTrack->SetMuonClusterMap(esdMuTrack->GetMuonClusterMap());
665 aodTrack->SetMatchTrigger(esdMuTrack->GetMatchTrigger());
667 if (primary != NULL) primary->AddDaughter(aodTrack);
669 if (esdMuTrack->Charge() > 0) nPosTracks++;
673 header->SetRefMultiplicity(jTracks);
674 header->SetRefMultiplicityPos(nPosTracks);
675 header->SetRefMultiplicityNeg(nNegTracks);
677 // Note: we do not have to call PostData for the AOD because the AliAnalysisTaskSE
678 // base class already does so just after the UserExec.
682 void AliAnalysisTaskLinkToMC::FindLinksToMC(TMap& links)
684 /// Finds all MC tracks that correspond to their ESD tracks and creates a
685 /// mapping between the ESD tracks and MC tracks which is filled into the
686 /// 'links' TMap object.
687 /// \param links This will be filled with the ESD and MC track pairs which
688 /// have been matched. Thus, for each ESD track which is set as the TMap key,
689 /// the corresponding MC track is identified by the TMap value. If no
690 /// MC track could be matched then the MC track value is set to NULL.
692 // Algorithm description for connecting reconstructed ESD muon tracks to MC:
694 // Build ESD track list
695 // Build MC track list
696 // while ESD track list not empty and MC track list not empty:
697 // Find ESD/MC track pair with best compatibility.
698 // if no pair is compatible then break loop.
699 // Store found pair in return list.
700 // Remove pair from ESD and MC track list.
701 // Mark remaining tracks in ESD track list as fakes.
703 // To find ESD/MC track pair with best compatibility:
704 // bestQuality := 1e100
705 // bestESDTrack := nil
706 // bestMCTrack := nil
707 // for each ESD track:
708 // for each MC track:
709 // chi^2 / numOfClustersMatched := calculate compatibility
710 // if tracks are compatible and chi^2 / numOfClustersMatched < bestQuality then:
711 // bestQuality := chi^2 / noOfClustersMatched
712 // bestESDTrack := current ESD track
713 // bestMCTrack := current MC track
715 // To calculate compatibility (ESD track, MC track):
717 // numberOfClustersAssociated := 0
718 // for each cluster in ESD track:
719 // find nearest 2 track refs from MC track
720 // if 2 track refs found then find average of them into MChit
721 // else MChit := nearest ref track.
722 // k := calculate chi^2 distance between ESD cluster and MChit
723 // if k <= maxValue then:
724 // chi^2 := chi^2 + k
725 // numberOfClustersAssociated := numberOfClustersAssociated + 1
727 // chi^2 := chi^2 + maxAllowedChi2
728 // if numberOfClustersAssociated >= minGoodClusters then:
729 // return (compatible, chi^2)
731 // return (not compatible, 1e100)
733 AliESDEvent* esd = static_cast<AliESDEvent*>(InputEvent());
737 for (Int_t i = 0; i < esd->GetNumberOfMuonTracks(); i++)
739 esdTracks.Add(esd->GetMuonTrack(i));
741 for (Int_t i = 0; i < MCEvent()->GetNumberOfTracks(); i++)
743 mcTracks.Add(MCEvent()->GetTrack(i));
746 while (not esdTracks.IsEmpty() and not mcTracks.IsEmpty())
748 Double_t bestQuality = 1e100;
749 Int_t bestNoOfClusters = 0;
750 AliESDMuonTrack* bestEsdTrack = NULL;
751 AliMCParticle* bestMcTrack = NULL;
752 TIter itesd(&esdTracks);
753 TIter itmc(&mcTracks);
754 AliESDMuonTrack* esdTrack = NULL;
755 AliMCParticle* mcTrack = NULL;
757 // Find the ESD and MC track pair that matches the best.
758 while ( (esdTrack = static_cast<AliESDMuonTrack*>(itesd())) != NULL )
760 while ( (mcTrack = static_cast<AliMCParticle*>(itmc())) != NULL )
762 Double_t chi2 = 1e100;
763 bool tracksCompatible = false;
764 Int_t noOfClusters = 0;
765 CalculateTrackCompatibility(esdTrack, mcTrack, tracksCompatible, chi2, noOfClusters);
766 Double_t quality = noOfClusters > 0 ? chi2 / Double_t(noOfClusters) : 1e100;
767 if (tracksCompatible and quality < bestQuality and noOfClusters >= bestNoOfClusters)
769 bestQuality = quality;
770 bestNoOfClusters = noOfClusters;
771 bestEsdTrack = esdTrack;
772 bestMcTrack = mcTrack;
778 // If no track pair could be matched then we are done because we
779 // will not be able to match anything more next time.
780 if (bestMcTrack == NULL) break;
782 links.Add(bestEsdTrack, bestMcTrack);
784 esdTracks.Remove(bestEsdTrack);
785 mcTracks.Remove(bestMcTrack);
788 // Add all remaining ESD tracks which could not be matched to any MC track
789 // to the mapping but with NULL for the MC track.
790 TIter itesd(&esdTracks);
791 AliESDMuonTrack* esdTrack = NULL;
792 while ( (esdTrack = static_cast<AliESDMuonTrack*>(itesd())) != NULL )
794 links.Add(esdTrack, NULL);
799 void AliAnalysisTaskLinkToMC::CalculateTrackCompatibility(
800 AliESDMuonTrack* esdTrack, AliMCParticle* mcTrack,
801 bool& tracksCompatible, Double_t& chi2, Int_t& noOfClusters
804 /// Calculates the compatibility between a ESD and MC track pair.
805 /// For each cluster of the ESD track, two track references are found
806 /// which are spatially closest to it. If the distance between the cluster
807 /// and track reference is larger than the HardCutLimit*() parameters then
808 /// the track reference is not used. The average coordinate is calculated
809 /// from the two track references found and the chi-square between the cluster
810 /// and average coordinate is calculated. If only one track reference is
811 /// found then it is used for the calculation without taking an average.
812 /// If no track references are found then the cluster is not matched and
813 /// we continue to the next cluster.
814 /// If the chi-squared value for the cluster / track reference pair is
815 /// larger than that allowed by the SigmaCut() parameter then the cluster
816 /// is also considered not to be matched.
817 /// The tracks are compatible if the ESD track has a minimum number of
818 /// clusters matched with MC track references from the MC track.
819 /// [in] \param esdTrack The ESD track we are trying to match.
820 /// [in] \param mcTrack The MC track we are trying to match the ESD track to.
821 /// [out] \param tracksCompatible Filled with true if the given esdTrack
822 /// and mcTrack are compatible.
823 /// [out] \param chi2 Filled with the total chi-squared calculated for all
824 /// matched clusters.
825 /// [out] \param noOfClusters Filled with the number of clusters that were
826 /// matched to at least one track reference.
830 tracksCompatible = false;
831 Int_t noOfClustersSt12 = 0;
832 Int_t noOfClustersSt45 = 0;
833 bool chamberMatched[14] = {
834 false, false, false, false, false, false, false,
835 false, false, false, false, false, false, false
837 bool stationMatched[7] = {
838 false, false, false, false, false, false, false
841 for (Int_t i = 0; i < esdTrack->GetNClusters(); i++)
843 AliESDMuonCluster* cluster = static_cast<AliESDMuonCluster*>( esdTrack->GetClusters()[i] );
844 Double_t varX = cluster->GetErrX2();
845 Double_t varY = cluster->GetErrY2();
846 // If the variance is zero then use the default one or just 1
847 // if the default is also not set.
848 if (varX == 0) varX = fVarX != 0 ? fVarX : 1.;
849 if (varY == 0) varY = fVarY != 0 ? fVarY : 1.;
851 Double_t chi2A = 1e100;
852 AliTrackReference* refA = NULL;
853 Double_t chi2B = 1e100;
854 AliTrackReference* refB = NULL;
856 for (Int_t j = 0; j < mcTrack->GetNumberOfTrackReferences(); j++)
858 AliTrackReference* ref = mcTrack->GetTrackReference(j);
860 if (ref->DetectorId() != AliTrackReference::kMUON) continue;
861 if (ref->UserId() != cluster->GetDetElemId()) continue;
863 Double_t dx = ref->X() - cluster->GetX();
864 if (TMath::Abs(dx) > fHardCutLimitX) continue;
865 Double_t dy = ref->Y() - cluster->GetY();
866 if (TMath::Abs(dy) > fHardCutLimitY) continue;
869 Double_t dz = ref->Z() - cluster->GetZ();
870 if (TMath::Abs(dz) > fHardCutLimitZ) continue;
873 Double_t chi2sum = dx*dx / varX + dy*dy / varY;
876 // Copy track reference A to B and set new value for A
882 else if (chi2sum < chi2B)
884 // Track reference A still has smaller chi2 so just replace B.
890 Double_t x = 0, y = 0;
891 if (refA != NULL and refB != NULL)
893 // Average the track references if 2 were found.
894 x = (refA->X() + refB->X()) * 0.5;
895 y = (refA->Y() + refB->Y()) * 0.5;
897 else if (refA != NULL)
904 continue; // no track reference hit found.
907 Double_t dx = x - cluster->GetX();
908 Double_t dy = y - cluster->GetY();
910 Double_t chi2x = dx*dx / varX;
911 if (chi2x > fSigmaCut2) continue;
912 Double_t chi2y = dy*dy / varY;
913 if (chi2y > fSigmaCut2) continue;
915 chi2 += chi2x + chi2y;
917 if (cluster->GetChamberId() >= 0 and cluster->GetChamberId() < 4)
921 if (cluster->GetChamberId() >= 6 and cluster->GetChamberId() < 10)
926 if (0 <= cluster->GetChamberId() and cluster->GetChamberId() < 14)
928 chamberMatched[cluster->GetChamberId()] = true;
929 stationMatched[cluster->GetChamberId()/2] = true;
933 // Need to check that all the chambers and stations that had to match we
934 // actually found matching track references.
936 bool forcedChambersMatched = true;
937 for (Int_t i = 0; i < 14; i++)
939 if (fChamberMustMatch[i] and not chamberMatched[i])
941 forcedChambersMatched = false;
944 bool forcedStationsMatched = true;
945 for (Int_t i = 0; i < 7; i++)
947 if (fStationMustMatch[i] and not stationMatched[i])
949 forcedStationsMatched = false;
953 if (noOfClusters >= fMinClusters and
954 noOfClustersSt12 >= fMinClustersInSt12 and
955 noOfClustersSt45 >= fMinClustersInSt45 and
956 forcedChambersMatched and forcedStationsMatched
959 tracksCompatible = true;