55dd806b7787822a2efa664599f645e152f1a68d
[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 //
17 // Add the muon tracks to the generic AOD track branch during the 
18 // filtering of the ESD. 
19 //
20 // Authors: R. Arnaldi 5/5/08 and L. Aphecetche January 2011
21 //
22 // Note that we :
23 //   - completely disable all the branches that are not required by (most) the muon analyses,
24 //     e.g. cascades, v0s, kinks, jets, etc...
25 //   - filter severely the tracks (keep only muon tracks) and vertices (keep only primary -including
26 //     pile-up - vertices) branches 
27 // 
28 // (see AddFilteredAOD method)
29 //
30
31 #include "AliAnalysisTaskESDMuonFilter.h"
32
33 #include "AliAODDimuon.h"
34 #include "AliAODEvent.h"
35 #include "AliAODHandler.h"
36 #include "AliAODMCParticle.h"
37 #include "AliAODMuonReplicator.h"
38 #include "AliAODVertex.h"
39 #include "AliAnalysisFilter.h"
40 #include "AliAnalysisManager.h"
41 #include "AliCodeTimer.h"
42 #include "AliESDEvent.h"
43 #include "AliESDInputHandler.h"
44 #include "AliESDMuonTrack.h"
45 #include "AliESDVertex.h"
46 #include "AliESDtrack.h"
47 #include "AliLog.h"
48 #include "AliMCEvent.h"
49 #include "AliMCEventHandler.h"
50 #include "AliMultiplicity.h"
51 #include "AliStack.h"
52 #include <TChain.h>
53 #include <TFile.h>
54 #include <TParticle.h>
55
56 ClassImp(AliAnalysisTaskESDMuonFilter)
57 ClassImp(AliAnalysisNonMuonTrackCuts)
58
59 ////////////////////////////////////////////////////////////////////////
60
61 AliAnalysisNonMuonTrackCuts::AliAnalysisNonMuonTrackCuts()
62 {
63   // default ctor 
64 }
65
66 Bool_t AliAnalysisNonMuonTrackCuts::IsSelected(TObject* obj)
67 {
68   // Returns true if the object is a muon track
69   AliAODTrack* track = dynamic_cast<AliAODTrack*>(obj);
70   if (track && track->IsMuonTrack()) return kTRUE;
71   return kFALSE;
72 }
73
74 AliAnalysisNonPrimaryVertices::AliAnalysisNonPrimaryVertices()
75 {
76   // default ctor   
77 }
78
79 Bool_t AliAnalysisNonPrimaryVertices::IsSelected(TObject* obj)
80 {
81   // Returns true if the object is a primary vertex
82
83   AliAODVertex* vertex = dynamic_cast<AliAODVertex*>(obj);
84   if (vertex)
85   {
86     if ( vertex->GetType() == AliAODVertex::kPrimary ||
87         vertex->GetType() == AliAODVertex::kMainSPD ||
88         vertex->GetType() == AliAODVertex::kPileupSPD ||
89         vertex->GetType() == AliAODVertex::kPileupTracks ||
90         vertex->GetType() == AliAODVertex::kMainTPC )
91     {
92       return kTRUE;
93     }
94   }
95   
96 //  enum AODVtx_t {kUndef=-1, kPrimary, kKink, kV0, kCascade, kMulti, kMainSPD, kPileupSPD, kPileupTracks,kMainTPC};
97
98   return kFALSE;
99   
100 }
101
102 AliAnalysisTaskESDMuonFilter::AliAnalysisTaskESDMuonFilter(Bool_t onlyMuon, Bool_t keepAllEvents):
103   AliAnalysisTaskSE(),
104   fTrackFilter(0x0),
105   fEnableMuonAOD(kFALSE),
106   fEnableDimuonAOD(kFALSE),
107   fOnlyMuon(onlyMuon),
108   fKeepAllEvents(keepAllEvents)
109 {
110   // Default constructor
111 }
112
113 AliAnalysisTaskESDMuonFilter::AliAnalysisTaskESDMuonFilter(const char* name, Bool_t onlyMuon, Bool_t keepAllEvents):
114   AliAnalysisTaskSE(name),
115   fTrackFilter(0x0),
116   fEnableMuonAOD(kFALSE),
117   fEnableDimuonAOD(kFALSE),
118   fOnlyMuon(onlyMuon),
119   fKeepAllEvents(keepAllEvents)
120 {
121   // Constructor
122 }
123
124 //______________________________________________________________________________
125 void AliAnalysisTaskESDMuonFilter::UserCreateOutputObjects()
126 {
127   // Create the output container
128   if (fTrackFilter) OutputTree()->GetUserInfo()->Add(fTrackFilter);
129 }
130
131 //______________________________________________________________________________
132 void AliAnalysisTaskESDMuonFilter::PrintTask(Option_t *option, Int_t indent) const
133 {
134   // Specify how we are configured
135   
136   AliAnalysisTaskSE::PrintTask(option,indent);
137   
138   TString spaces(' ',indent+3);
139   
140   if ( fOnlyMuon ) 
141   {
142     cout << spaces.Data() << "Keep only muon information " << endl;        
143   }
144   else 
145   {
146     cout << spaces.Data() << "Keep all information from standard AOD" << endl;
147   }
148
149   if ( fKeepAllEvents ) 
150   {
151     cout << spaces.Data() << "Keep all events, regardless of number of muons" << endl;    
152   }
153   else 
154   {
155     cout << spaces.Data() << "Keep only events with at least one muon" << endl;
156   }
157 }
158
159 //______________________________________________________________________________
160 void AliAnalysisTaskESDMuonFilter::AddFilteredAOD(const char* aodfilename, const char* title)
161 {
162   // Add an output filtered and replicated aod
163   
164   AliAODHandler *aodH = (AliAODHandler*)((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
165   if (!aodH) Fatal("UserCreateOutputObjects", "No AOD handler");
166
167   AliAODExtension* ext = aodH->AddFilteredAOD(aodfilename,title);
168
169   if (!ext) return;
170   
171   if ( fOnlyMuon ) 
172   {    
173     AliAODBranchReplicator* murep = new AliAODMuonReplicator("MuonReplicator",
174                                                              "remove non muon tracks and non primary or pileup vertices",
175                                                              new AliAnalysisNonMuonTrackCuts,
176                                                              new AliAnalysisNonPrimaryVertices);
177     
178     ext->DropUnspecifiedBranches(); // all branches not part of a FilterBranch call (below) will be dropped
179     
180     ext->FilterBranch("tracks",murep);    
181     ext->FilterBranch("vertices",murep);  
182     ext->FilterBranch("dimuons",murep);
183   }  
184 }
185
186 //______________________________________________________________________________
187 void AliAnalysisTaskESDMuonFilter::Init()
188 {
189   // Initialization
190   if(fEnableMuonAOD) AddFilteredAOD("AliAOD.Muons.root", "MuonEvents");
191   if(fEnableDimuonAOD) AddFilteredAOD("AliAOD.Dimuons.root", "DimuonEvents");    
192 }
193
194
195 //______________________________________________________________________________
196 void AliAnalysisTaskESDMuonFilter::UserExec(Option_t */*option*/)
197 {
198   // Execute analysis for current event                                     
199   
200   Long64_t ientry = Entry();
201   if(fDebug)printf("Muon Filter: Analysing event # %5d\n", (Int_t) ientry);
202   
203   ConvertESDtoAOD();
204 }
205
206 //______________________________________________________________________________
207 void AliAnalysisTaskESDMuonFilter::ConvertESDtoAOD() 
208 {
209   // ESD Muon Filter analysis task executed for each event
210   
211   AliCodeTimerAuto("",0);
212   
213   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
214   
215   // Fetch Stack for debuggging if available 
216   AliStack *pStack = 0;
217   AliMCEventHandler *mcH = 0;
218   if(MCEvent()){
219     pStack = MCEvent()->Stack();
220     mcH = (AliMCEventHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler()); 
221   }
222     
223   // Define arrays for muons
224   Double_t pos[3];
225   Double_t p[3];
226   Double_t pid[10];
227   
228   // has to be changed once the muon pid is provided by the ESD
229   for (Int_t i = 0; i < 10; pid[i++] = 0.) {}
230   pid[AliAODTrack::kMuon]=1.;
231   
232   AliAODHeader* header = AODEvent()->GetHeader();
233   AliAODTrack *aodTrack = 0x0;
234   AliESDMuonTrack *esdMuTrack = 0x0;
235   
236   // Access to the AOD container of tracks
237   TClonesArray &tracks = *(AODEvent()->GetTracks());
238   Int_t jTracks = header->GetRefMultiplicity();
239   
240   // Read primary vertex from AOD event 
241   AliAODVertex *primary = AODEvent()->GetPrimaryVertex();
242   if(fDebug)primary->Print();
243   
244   // Loop on muon tracks to fill the AOD track branch
245   Int_t nMuTracks = esd->GetNumberOfMuonTracks();
246
247   for (Int_t iTrack=0; iTrack<nMuTracks; ++iTrack) esd->GetMuonTrack(iTrack)->SetESDEvent(esd);
248   
249   // Update number of positive and negative tracks from AOD event (M.G.)
250   Int_t nPosTracks = header->GetRefMultiplicityPos();
251   Int_t nNegTracks = header->GetRefMultiplicityNeg();
252   
253   // Access to the AOD container of dimuons
254   TClonesArray &dimuons = *(AODEvent()->GetDimuons());
255   AliAODDimuon *aodDimuon = 0x0;
256   
257   Int_t nMuons=0;
258   Int_t nDimuons=0;
259   Int_t jDimuons=0;
260   Int_t nMuonTrack[100];
261   
262   for(int imuon=0;imuon<100;imuon++) nMuonTrack[imuon]=0;
263   
264   for (Int_t nMuTrack = 0; nMuTrack < nMuTracks; ++nMuTrack)
265   {
266     esdMuTrack = esd->GetMuonTrack(nMuTrack);
267     
268     if (!esdMuTrack->ContainTrackerData()) continue;
269     
270     UInt_t selectInfo = 0;
271     // Track selection
272     if (fTrackFilter) {
273         selectInfo = fTrackFilter->IsSelected(esdMuTrack);
274         if (!selectInfo) {
275           continue;
276         }  
277     }
278     
279     p[0] = esdMuTrack->Px(); 
280     p[1] = esdMuTrack->Py(); 
281     p[2] = esdMuTrack->Pz();
282     
283     pos[0] = esdMuTrack->GetNonBendingCoor(); 
284     pos[1] = esdMuTrack->GetBendingCoor(); 
285     pos[2] = esdMuTrack->GetZ();
286     
287     if(mcH)mcH->SelectParticle(esdMuTrack->GetLabel());
288     
289     aodTrack = new(tracks[jTracks++]) AliAODTrack(esdMuTrack->GetUniqueID(), // ID
290                                                   esdMuTrack->GetLabel(), // label
291                                                   p, // momentum
292                                                   kTRUE, // cartesian coordinate system
293                                                   pos, // position
294                                                   kFALSE, // isDCA
295                                                   0x0, // covariance matrix
296                                                   esdMuTrack->Charge(), // charge
297                                                   0, // ITSClusterMap
298                                                   pid, // pid
299                                                   primary, // primary vertex
300                                                   kFALSE, // used for vertex fit?
301                                                   kFALSE, // used for primary vertex fit?
302                                                   AliAODTrack::kPrimary,// track type
303                                                   selectInfo); 
304     
305     aodTrack->SetXYAtDCA(esdMuTrack->GetNonBendingCoorAtDCA(), esdMuTrack->GetBendingCoorAtDCA());
306     aodTrack->SetPxPyPzAtDCA(esdMuTrack->PxAtDCA(), esdMuTrack->PyAtDCA(), esdMuTrack->PzAtDCA());
307     aodTrack->SetRAtAbsorberEnd(esdMuTrack->GetRAtAbsorberEnd());
308     aodTrack->ConvertAliPIDtoAODPID();
309     aodTrack->SetChi2perNDF(esdMuTrack->GetChi2() / (2.*esdMuTrack->GetNHit() - 5.));
310     aodTrack->SetChi2MatchTrigger(esdMuTrack->GetChi2MatchTrigger());
311     aodTrack->SetHitsPatternInTrigCh(esdMuTrack->GetHitsPatternInTrigCh());
312     aodTrack->SetMuonClusterMap(esdMuTrack->GetMuonClusterMap());
313     aodTrack->SetMatchTrigger(esdMuTrack->GetMatchTrigger());
314     aodTrack->Connected(esdMuTrack->IsConnected());
315     primary->AddDaughter(aodTrack);
316     
317     if (esdMuTrack->Charge() > 0) nPosTracks++;
318     else nNegTracks++;
319     
320     nMuonTrack[nMuons]= jTracks-1.;
321     ++nMuons;
322   }
323   
324   if(nMuons>=2) 
325   { 
326     for(int i=0;i<nMuons;i++){
327       Int_t index0 = nMuonTrack[i];
328       for(int j=i+1;j<nMuons;j++){
329         Int_t index1 = nMuonTrack[j];
330         aodDimuon = new(dimuons[jDimuons++]) AliAODDimuon(tracks.At(index0),tracks.At(index1));
331         ++nDimuons;
332         if (fDebug > 1){
333           AliAODDimuon *dimuon0 = (AliAODDimuon*)dimuons.At(jDimuons-1);
334           printf("Dimuon: mass = %f, px=%f, py=%f, pz=%f\n",dimuon0->M(),dimuon0->Px(),dimuon0->Py(),dimuon0->Pz());  
335           AliAODTrack  *mu0 = (AliAODTrack*) dimuon0->GetMu(0);
336           AliAODTrack  *mu1 = (AliAODTrack*) dimuon0->GetMu(1);
337           printf("Muon0 px=%f py=%f pz=%f\n",mu0->Px(),mu0->Py(),mu0->Pz());
338           printf("Muon1 px=%f py=%f pz=%f\n",mu1->Px(),mu1->Py(),mu1->Pz());
339         }  
340       }
341     }
342   }
343   
344   
345   header->SetRefMultiplicity(jTracks); 
346   header->SetRefMultiplicityPos(nPosTracks);
347   header->SetRefMultiplicityNeg(nNegTracks);
348   header->SetNumberOfMuons(nMuons);
349   header->SetNumberOfDimuons(nDimuons);
350   
351   if ( fEnableMuonAOD && ( (nMuons>0) || fKeepAllEvents ) )
352   {
353     AliAODExtension *extMuons = dynamic_cast<AliAODHandler*>
354     ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler())->GetFilteredAOD("AliAOD.Muons.root");
355     extMuons->SelectEvent();
356   }
357   
358   if ( fEnableDimuonAOD && ( (nMuons>1) || fKeepAllEvents )  )
359   {
360     AliAODExtension *extDimuons = dynamic_cast<AliAODHandler*>
361     ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler())->GetFilteredAOD("AliAOD.Dimuons.root");
362     extDimuons->SelectEvent();
363   }
364   
365 }
366
367 void AliAnalysisTaskESDMuonFilter::Terminate(Option_t */*option*/)
368 {
369   // Terminate analysis
370   //
371   if (fDebug > 1) printf("AnalysisESDfilter: Terminate() \n");
372 }
373
374 void  AliAnalysisTaskESDMuonFilter::PrintMCInfo(AliStack *pStack,Int_t label)
375 {
376   // print mc info
377   if(!pStack)return;
378   label = TMath::Abs(label);
379   TParticle *part = pStack->Particle(label);
380   Printf("########################");
381   Printf("%s:%d %d UniqueID %d PDG %d P %3.3f",(char*)__FILE__,__LINE__,label,part->GetUniqueID(),part->GetPdgCode(),part->P());
382   part->Print();
383   TParticle* mother = part;
384   Int_t imo = part->GetFirstMother();
385   Int_t nprim = pStack->GetNprimary();
386   //  while((imo >= nprim) && (mother->GetUniqueID() == 4)) {
387   while((imo >= nprim)) {
388     mother =  pStack->Particle(imo);
389     Printf("Mother %s:%d Label %d UniqueID %d PDG %d P %3.3f",(char*)__FILE__,__LINE__,imo,mother->GetUniqueID(),mother->GetPdgCode(),mother->P());
390     mother->Print();
391     imo =  mother->GetFirstMother();
392   }
393   Printf("########################");
394 }