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