Reseting the MC label to a default value of -1 before evaluating them (Philippe P...
[u/mrichter/AliRoot.git] / PWG3 / muondep / AliAnalysisTaskESDMCLabelAddition.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 #include <TChain.h>
17 #include <TFile.h>
18 #include <TParticle.h>
19
20 #include "AliAnalysisTaskESDMCLabelAddition.h"
21 #include "AliAnalysisManager.h"
22 #include "AliESDEvent.h"
23 #include "AliAODEvent.h"
24 #include "AliESDInputHandler.h"
25 #include "AliAODHandler.h"
26 #include "AliAnalysisFilter.h"
27 #include "AliESDtrack.h"
28 #include "AliESDMuonTrack.h"
29 #include "AliESDMuonCluster.h"
30 #include "AliESDVertex.h"
31 #include "AliMultiplicity.h"
32 #include "AliLog.h"
33 #include "AliStack.h"
34 #include "AliMCEvent.h"
35 #include "AliMCEventHandler.h"
36 #include "AliAODMCParticle.h"
37
38 #include "AliMUONRecoCheck.h"
39 #include "AliMUONESDInterface.h"
40 #include "AliMUONTrack.h"
41 #include "AliMUONTrackParam.h"
42 #include "AliMUONVTrackStore.h"
43 #include "AliMUONVCluster.h"
44 #include "AliMUONVClusterStore.h"
45
46 ClassImp(AliAnalysisTaskESDMCLabelAddition)
47
48 // sigma cut applied to match a reconstructed cluster with a trackref
49 const Double_t AliAnalysisTaskESDMCLabelAddition::fgkSigmaCut = 10.;
50
51 //----------------------------------------------------------------------
52 AliAnalysisTaskESDMCLabelAddition::AliAnalysisTaskESDMCLabelAddition():
53   AliAnalysisTaskSE()
54 {
55   // Default constructor
56 }
57
58
59 //----------------------------------------------------------------------
60 AliAnalysisTaskESDMCLabelAddition::AliAnalysisTaskESDMCLabelAddition(const char* name):
61   AliAnalysisTaskSE(name)
62 {
63   // Constructor
64 }
65
66
67 //----------------------------------------------------------------------
68 void AliAnalysisTaskESDMCLabelAddition::UserCreateOutputObjects()
69 {
70 }
71
72
73 //----------------------------------------------------------------------
74 void AliAnalysisTaskESDMCLabelAddition::Init()
75 {
76   AliDebug(2, "Init()");
77 }
78
79
80 //----------------------------------------------------------------------
81 void AliAnalysisTaskESDMCLabelAddition::UserExec(Option_t */*option*/)
82 {
83   // Execute analysis for current event                                     
84   Long64_t ientry = Entry();
85   AliDebug(1, Form("MCLabel Addition: Analysing event # %5d\n",(Int_t) ientry)); 
86   AddMCLabel();
87 }
88
89
90 //----------------------------------------------------------------------
91 void AliAnalysisTaskESDMCLabelAddition::AddMCLabel() 
92 {
93   // Load ESD event
94   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
95   
96   // Load MC event 
97   AliMCEventHandler *mcH = 0;
98   if(MCEvent()) mcH = (AliMCEventHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler()); 
99   
100   // Get reference tracks
101   AliMUONRecoCheck rc(esd,mcH);
102   AliMUONVTrackStore* trackRefStore = rc.TrackRefs(-1);
103   
104   // Loop over reconstructed tracks
105   AliESDMuonTrack *esdTrack = 0x0;
106   Int_t nMuTracks = esd->GetNumberOfMuonTracks();
107   for (Int_t nMuTrack = 0; nMuTrack < nMuTracks; ++nMuTrack) {
108     
109     esdTrack = esd->GetMuonTrack(nMuTrack);
110     esdTrack->SetLabel(-1);
111     
112     // skip ghosts
113     if (!esdTrack->ContainTrackerData()) continue;
114     
115     // try to match the reconstructed track with a simulated one
116     AliMUONTrack* matchedTrackRef = MatchWithTrackRef(*esdTrack, *trackRefStore);
117     
118     if (matchedTrackRef) {
119       
120       // set the MC label
121       esdTrack->SetLabel(matchedTrackRef->GetUniqueID());
122        
123       // remove already matched trackRefs
124       trackRefStore->Remove(*matchedTrackRef);
125       
126     }
127     
128   }
129   
130 }
131
132
133 //----------------------------------------------------------------------
134 void AliAnalysisTaskESDMCLabelAddition::Terminate(Option_t */*option*/)
135 {
136   // Terminate analysis
137   //
138   AliDebug(2, "Terminate()");
139 }
140
141
142 //----------------------------------------------------------------------
143 AliMUONTrack* AliAnalysisTaskESDMCLabelAddition::ESDToMUON(AliESDMuonTrack &esdTrack)
144 {
145   /// Convert an ESD track into a MUON track with dummy parameters (just to fill the clusters).
146   /// It is the responsability of the user to delete the track afterward
147   
148   AliMUONTrack *track = new AliMUONTrack();
149   
150   // ckeck whether the ESD track contains clusters
151   if(!esdTrack.ClustersStored()) return track;
152   
153   // track parameters at first cluster
154   AliMUONTrackParam param;
155   AliMUONESDInterface::GetParamAtFirstCluster(esdTrack, param);
156   AliMUONESDInterface::GetParamCov(esdTrack, param);
157   
158   // create empty cluster
159   AliMUONVClusterStore* cStore = AliMUONESDInterface::NewClusterStore();
160   AliMUONVCluster* cluster = cStore->CreateCluster(0,0,0);
161   
162   // loop over ESD clusters
163   AliESDMuonCluster *esdCluster = (AliESDMuonCluster*) esdTrack.GetClusters().First();
164   while (esdCluster) {
165     
166     // copy cluster information
167     AliMUONESDInterface::ESDToMUON(*esdCluster, *cluster);
168     
169     // only set the Z parameter to avoid error in the AddTrackParamAtCluster(...) method
170     param.SetZ(cluster->GetZ());
171     
172     // add common track parameters at current cluster
173     track->AddTrackParamAtCluster(param, *cluster, kTRUE);
174     
175     esdCluster = (AliESDMuonCluster*) esdTrack.GetClusters().After(esdCluster);
176   }
177   
178   // clean memory
179   delete cluster;
180   delete cStore;
181   
182   return track;
183   
184 }
185
186
187 //----------------------------------------------------------------------
188 AliMUONTrack* AliAnalysisTaskESDMCLabelAddition::MatchWithTrackRef(AliESDMuonTrack &esdTrack,
189                                                                    AliMUONVTrackStore &trackRefStore)
190 {
191   /// Return the trackRef matched with the reconstructed track and the fraction of matched clusters
192   
193   AliMUONTrack *matchedTrackRef = 0x0;
194   
195   // convert ESD track to MUON track
196   AliMUONTrack *track = ESDToMUON(esdTrack);
197   
198   // look for the corresponding simulated track if any
199   TIter next(trackRefStore.CreateIterator());
200   AliMUONTrack* trackRef;
201   while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) {
202     
203     // check compatibility
204     if (TrackMatched(*track, *trackRef)) {
205       matchedTrackRef = trackRef;
206       break;
207     }
208     
209   }
210   
211   // clean memory
212   delete track;
213   
214   return matchedTrackRef;
215   
216 }
217
218
219 //----------------------------------------------------------------------
220 Bool_t AliAnalysisTaskESDMCLabelAddition::TrackMatched(AliMUONTrack &track, AliMUONTrack &trackRef)
221 {
222   /// Try to match 2 tracks
223   
224   Bool_t compTrack[10];
225   Int_t nMatchClusters = track.CompatibleTrack(&trackRef, fgkSigmaCut, compTrack);
226   Double_t fractionOfMatchCluster = ((Double_t)nMatchClusters) / ((Double_t)track.GetNClusters());
227   
228   if ((compTrack[0] || compTrack[1] || compTrack[2] || compTrack[3]) && // at least 1 cluster matched in st 1 & 2
229       (compTrack[6] || compTrack[7] || compTrack[8] || compTrack[9]) && // at least 1 cluster matched in st 4 & 5
230       fractionOfMatchCluster > 0.5) return kTRUE;                       // more than 50% of clusters matched
231   else return kFALSE;
232   
233 }
234