Update for AliAODHeader. (R. Arnaldi)
[u/mrichter/AliRoot.git] / PWG3 / muon / AliAnalysisTaskESDMuonFilter.cxx
CommitLineData
1d1ea3ce 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
aba11173 18
1d1ea3ce 19#include <TChain.h>
20#include <TFile.h>
fc3a4c45 21#include <TParticle.h>
1d1ea3ce 22
23#include "AliAnalysisTaskESDMuonFilter.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 "AliESDVertex.h"
33#include "AliMultiplicity.h"
34#include "AliLog.h"
fc3a4c45 35#include "AliStack.h"
36#include "AliMCEvent.h"
37#include "AliMCEventHandler.h"
38#include "AliAODMCParticle.h"
866d8d78 39#include "AliAODDimuon.h"
1d1ea3ce 40
41ClassImp(AliAnalysisTaskESDMuonFilter)
42
43////////////////////////////////////////////////////////////////////////
44
45AliAnalysisTaskESDMuonFilter::AliAnalysisTaskESDMuonFilter():
aba11173 46 AliAnalysisTaskSE(),
f785fc59 47 fTrackFilter(0x0),
4a497116 48 fEnableMuonAOD(kFALSE),
49 fEnableDimuonAOD(kFALSE)
1d1ea3ce 50{
51 // Default constructor
52}
53
54AliAnalysisTaskESDMuonFilter::AliAnalysisTaskESDMuonFilter(const char* name):
aba11173 55 AliAnalysisTaskSE(name),
f785fc59 56 fTrackFilter(0x0),
4a497116 57 fEnableMuonAOD(kFALSE),
58 fEnableDimuonAOD(kFALSE)
1d1ea3ce 59{
60 // Constructor
61}
62
63void AliAnalysisTaskESDMuonFilter::UserCreateOutputObjects()
64{
65 // Create the output container
aba11173 66 if (fTrackFilter) OutputTree()->GetUserInfo()->Add(fTrackFilter);
1d1ea3ce 67}
68
69void AliAnalysisTaskESDMuonFilter::Init()
70{
71 // Initialization
aba11173 72 if (fDebug > 1) AliInfo("Init() \n");
f785fc59 73 // From Andrei
74 AliAODHandler *aodH = (AliAODHandler*)((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
75 if (!aodH) Fatal("UserCreateOutputObjects", "No AOD handler. Aborting.");
76 if(fEnableMuonAOD)aodH->AddFilteredAOD("AliAOD.Muons.root", "MuonEvents");
866d8d78 77 if(fEnableDimuonAOD)aodH->AddFilteredAOD("AliAOD.Dimuons.root", "DimuonEvents");
1d1ea3ce 78}
79
80
81void AliAnalysisTaskESDMuonFilter::UserExec(Option_t */*option*/)
82{
83 // Execute analysis for current event
84 Long64_t ientry = Entry();
edee4062 85 if(fDebug)printf("Muon Filter: Analysing event # %5d\n", (Int_t) ientry);
aba11173 86
1d1ea3ce 87 ConvertESDtoAOD();
88}
89
90void AliAnalysisTaskESDMuonFilter::ConvertESDtoAOD()
91{
92 // ESD Muon Filter analysis task executed for each event
93 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
aba11173 94
fc3a4c45 95 // Fetch Stack for debuggging if available
96 AliStack *pStack = 0;
97 AliMCEventHandler *mcH = 0;
98 if(MCEvent()){
99 pStack = MCEvent()->Stack();
100 mcH = (AliMCEventHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
101 }
102
1d1ea3ce 103 // Define arrays for muons
104 Double_t pos[3];
105 Double_t p[3];
106 Double_t pid[10];
aba11173 107
108 // has to be changed once the muon pid is provided by the ESD
a6e0ebfe 109 for (Int_t i = 0; i < 10; pid[i++] = 0.) {}
aba11173 110 pid[AliAODTrack::kMuon]=1.;
111
1d1ea3ce 112 AliAODHeader* header = AODEvent()->GetHeader();
113 AliAODTrack *aodTrack = 0x0;
aba11173 114 AliESDMuonTrack *esdMuTrack = 0x0;
115
1d1ea3ce 116 // Access to the AOD container of tracks
117 TClonesArray &tracks = *(AODEvent()->GetTracks());
118 Int_t jTracks = header->GetRefMultiplicity();
119
120 // Read primary vertex from AOD event
aba11173 121 AliAODVertex *primary = AODEvent()->GetPrimaryVertex();
edee4062 122 if(fDebug)primary->Print();
aba11173 123
1d1ea3ce 124 // Loop on muon tracks to fill the AOD track branch
125 Int_t nMuTracks = esd->GetNumberOfMuonTracks();
866d8d78 126
daf100e1 127 for (Int_t iTrack=0; iTrack<nMuTracks; ++iTrack) esd->GetMuonTrack(iTrack)->SetESDEvent(esd);
aba11173 128
129 // Update number of positive and negative tracks from AOD event (M.G.)
1d1ea3ce 130 Int_t nPosTracks = header->GetRefMultiplicityPos();
131 Int_t nNegTracks = header->GetRefMultiplicityNeg();
aba11173 132
866d8d78 133 // Access to the AOD container of dimuons
134 TClonesArray &dimuons = *(AODEvent()->GetDimuons());
135 AliAODDimuon *aodDimuon = 0x0;
136
f785fc59 137 Bool_t MuonsExist = kFALSE;
866d8d78 138 Bool_t DimuonsExist = kFALSE;
866d8d78 139 Int_t nMuons=0;
29e7e51b 140 Int_t nDimuons=0;
866d8d78 141 Int_t jDimuons=0;
87ac315c 142 Int_t nMuonTrack[10];
143
144 for(int imuon=0;imuon<10;imuon++) nMuonTrack[imuon]=0;
f785fc59 145
1d1ea3ce 146 for (Int_t nMuTrack = 0; nMuTrack < nMuTracks; ++nMuTrack) {
aba11173 147 esdMuTrack = esd->GetMuonTrack(nMuTrack);
148
149 if (!esdMuTrack->ContainTrackerData()) continue;
daf100e1 150
151 MuonsExist = kTRUE;
152
aba11173 153 UInt_t selectInfo = 0;
154 // Track selection
155 if (fTrackFilter) {
156 selectInfo = fTrackFilter->IsSelected(esdMuTrack);
157 if (!selectInfo) {
158 continue;
159 }
87ac315c 160 }
aba11173 161
162 p[0] = esdMuTrack->Px();
163 p[1] = esdMuTrack->Py();
164 p[2] = esdMuTrack->Pz();
165
166 pos[0] = esdMuTrack->GetNonBendingCoor();
167 pos[1] = esdMuTrack->GetBendingCoor();
168 pos[2] = esdMuTrack->GetZ();
169
fc3a4c45 170 if(mcH)mcH->SelectParticle(esdMuTrack->GetLabel());
171
aba11173 172 aodTrack = new(tracks[jTracks++]) AliAODTrack(esdMuTrack->GetUniqueID(), // ID
2e2d0c44 173 esdMuTrack->GetLabel(), // label
aba11173 174 p, // momentum
175 kTRUE, // cartesian coordinate system
176 pos, // position
177 kFALSE, // isDCA
178 0x0, // covariance matrix
179 esdMuTrack->Charge(), // charge
180 0, // ITSClusterMap
181 pid, // pid
182 primary, // primary vertex
183 kFALSE, // used for vertex fit?
184 kFALSE, // used for primary vertex fit?
185 AliAODTrack::kPrimary,// track type
186 selectInfo);
187
188 aodTrack->SetXYAtDCA(esdMuTrack->GetNonBendingCoorAtDCA(), esdMuTrack->GetBendingCoorAtDCA());
189 aodTrack->SetPxPyPzAtDCA(esdMuTrack->PxAtDCA(), esdMuTrack->PyAtDCA(), esdMuTrack->PzAtDCA());
190 aodTrack->ConvertAliPIDtoAODPID();
191 aodTrack->SetChi2perNDF(esdMuTrack->GetChi2() / (2.*esdMuTrack->GetNHit() - 5.));
192 aodTrack->SetChi2MatchTrigger(esdMuTrack->GetChi2MatchTrigger());
193 aodTrack->SetHitsPatternInTrigCh(esdMuTrack->GetHitsPatternInTrigCh());
194 aodTrack->SetMuonClusterMap(esdMuTrack->GetMuonClusterMap());
195 aodTrack->SetMatchTrigger(esdMuTrack->GetMatchTrigger());
5c15a68b 196 aodTrack->Connected(esdMuTrack->IsConnected());
aba11173 197 primary->AddDaughter(aodTrack);
198
199 if (esdMuTrack->Charge() > 0) nPosTracks++;
200 else nNegTracks++;
866d8d78 201
87ac315c 202 nMuonTrack[nMuons]= jTracks-1.;
866d8d78 203 nMuons++;
87ac315c 204 }
205
206 if(nMuons>=2) DimuonsExist = kTRUE;
866d8d78 207 if(DimuonsExist) {
87ac315c 208 for(int i=0;i<nMuons;i++){
209 Int_t index0 = nMuonTrack[i];
210 for(int j=i+1;j<nMuons;j++){
211 Int_t index1 = nMuonTrack[j];
212 aodDimuon = new(dimuons[jDimuons++]) AliAODDimuon(tracks.At(index0),tracks.At(index1));
29e7e51b 213 nDimuons++;
87ac315c 214 if (fDebug > 1){
87ac315c 215 AliAODDimuon *dimuon0 = (AliAODDimuon*)dimuons.At(jDimuons-1);
216 printf("Dimuon: mass = %f, px=%f, py=%f, pz=%f\n",dimuon0->M(),dimuon0->Px(),dimuon0->Py(),dimuon0->Pz());
217 AliAODTrack *mu0 = (AliAODTrack*) dimuon0->GetMu(0);
218 AliAODTrack *mu1 = (AliAODTrack*) dimuon0->GetMu(1);
219 printf("Muon0 px=%f py=%f pz=%f\n",mu0->Px(),mu0->Py(),mu0->Pz());
220 printf("Muon1 px=%f py=%f pz=%f\n",mu1->Px(),mu1->Py(),mu1->Pz());
221 }
222 }
223 }
866d8d78 224 }
225
aba11173 226
1d1ea3ce 227 header->SetRefMultiplicity(jTracks);
228 header->SetRefMultiplicityPos(nPosTracks);
229 header->SetRefMultiplicityNeg(nNegTracks);
29e7e51b 230 header->SetNumberOfMuons(nMuons);
231 header->SetNumberOfDimuons(nDimuons);
866d8d78 232
f785fc59 233 // From Andrei
234 if(fEnableMuonAOD && MuonsExist){
235 AliAODExtension *extMuons = dynamic_cast<AliAODHandler*>
236 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler())->GetFilteredAOD("AliAOD.Muons.root");
237 extMuons->SelectEvent();
866d8d78 238 }
239
240 if(fEnableDimuonAOD && DimuonsExist){
241 AliAODExtension *extDimuons = dynamic_cast<AliAODHandler*>
242 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler())->GetFilteredAOD("AliAOD.Dimuons.root");
243 extDimuons->SelectEvent();
f785fc59 244 }
245
1d1ea3ce 246}
247
248void AliAnalysisTaskESDMuonFilter::Terminate(Option_t */*option*/)
249{
aba11173 250 // Terminate analysis
251 //
252 if (fDebug > 1) printf("AnalysisESDfilter: Terminate() \n");
1d1ea3ce 253}
aba11173 254
fc3a4c45 255void AliAnalysisTaskESDMuonFilter::PrintMCInfo(AliStack *pStack,Int_t label){
256 if(!pStack)return;
257 label = TMath::Abs(label);
258 TParticle *part = pStack->Particle(label);
259 Printf("########################");
260 Printf("%s:%d %d UniqueID %d PDG %d P %3.3f",(char*)__FILE__,__LINE__,label,part->GetUniqueID(),part->GetPdgCode(),part->P());
261 part->Print();
262 TParticle* mother = part;
263 Int_t imo = part->GetFirstMother();
264 Int_t nprim = pStack->GetNprimary();
265 // while((imo >= nprim) && (mother->GetUniqueID() == 4)) {
266 while((imo >= nprim)) {
267 mother = pStack->Particle(imo);
268 Printf("Mother %s:%d Label %d UniqueID %d PDG %d P %3.3f",(char*)__FILE__,__LINE__,imo,mother->GetUniqueID(),mother->GetPdgCode(),mother->P());
269 mother->Print();
270 imo = mother->GetFirstMother();
271 }
272 Printf("########################");
273}