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