c4a2bdaf03d1a8975e9db79e61f1167f144a7638
[u/mrichter/AliRoot.git] / PWG3 / muondep / AliAnalysisTaskRecoCheck.cxx
1 #include "TChain.h"
2 #include "TClonesArray.h"
3 #include "TLorentzVector.h"
4 #include "TIterator.h"
5
6 #include "AliStack.h"
7 #include "AliMagFMaps.h"
8 #include "AliTracker.h"
9
10 #include "AliAnalysisTask.h"
11 #include "AliAnalysisManager.h"
12
13 #include "AliESDEvent.h"
14 #include "AliESDInputHandler.h"
15 #include "AliESDMuonTrack.h"
16
17 #include "AliMCEventHandler.h"
18 #include "AliMCEvent.h"
19
20 #include "AliAnalysisTaskRecoCheck.h"
21
22 #include "AliMUONTrackLight.h"
23 #include "AliMUONPairLight.h"
24 #include "AliMUONVTrackStore.h"
25 #include "AliMUONTrack.h"
26 #include "AliMUONRecoCheck.h"
27 #include "AliMUONTrackExtrap.h"
28
29 // analysis task for decoding reconstructed tracks and kinematics (AliMUONRecoCheck)
30 // Authors: Bogdan Vulpescu
31
32 ClassImp(AliAnalysisTaskRecoCheck)
33
34 //________________________________________________________________________
35 AliAnalysisTaskRecoCheck::AliAnalysisTaskRecoCheck(const char *name) 
36   : AliAnalysisTask(name, ""), fESDEvent(0), fTree(0), fArray1Mu(0), fArray2Mu(0), fL3Current(30000.0)
37 {
38   // Constructor
39
40   // Define input and output slots here
41   // Input slot #0 works with a TChain
42   DefineInput(0, TChain::Class());
43   // Output slot #0 writes into a TTree container
44   DefineOutput(0, TTree::Class());
45
46 }
47
48 //________________________________________________________________________
49 void AliAnalysisTaskRecoCheck::ConnectInputData(Option_t *) 
50 {
51   // Connect ESD or AOD here
52   // Called once
53
54   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
55   if (!tree) {
56     Printf("ERROR: Could not read chain from input slot 0");
57   } else {
58     tree->SetBranchStatus("*", kFALSE);
59     tree->SetBranchStatus("MuonTracks.*", kTRUE);
60     tree->SetBranchStatus("AliESDHeader.*", kTRUE);
61
62     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
63
64     if (!esdH) {
65       Printf("ERROR: Could not get ESDInputHandler");
66     } else
67       fESDEvent = esdH->GetEvent();
68   }
69   
70   // calculate the filed map in the L3 magnet using the current value
71   AliMagFMaps* field = 0x0;
72   if (fL3Current == 30000.0) {
73     field = new AliMagFMaps("Maps","Maps", 1, 1., 10., AliMagFMaps::k5kG);
74   }
75
76   // set the tracker field map
77   AliTracker::SetFieldMap(field, kFALSE);
78   AliMUONTrackExtrap::SetField(AliTracker::GetFieldMap());
79
80 }
81
82 //________________________________________________________________________
83 void AliAnalysisTaskRecoCheck::CreateOutputObjects() 
84 {
85   // Create output tree with branches of arrays
86   // Called once
87
88   fArray1Mu = new TClonesArray("AliMUONTrackLight",0); 
89   fArray2Mu = new TClonesArray("AliMUONPairLight",0); 
90   fTree = new TTree("tree","tree",0); 
91   fTree->Branch("muons",  &fArray1Mu); 
92   fTree->Branch("dimuons",&fArray2Mu); 
93   fTree->SetDirectory(0);
94   
95 }
96
97 //________________________________________________________________________
98 void AliAnalysisTaskRecoCheck::Exec(Option_t *) 
99 {
100   // Main loop
101   // Called for each event
102   
103   // Process MC truth, therefore we receive the AliAnalysisManager and ask it for the AliMCEventHandler
104   // This handler can return the current MC event
105   
106   Int_t nTracks = fESDEvent->GetNumberOfMuonTracks();
107   if (nTracks == 0) return;
108   
109   AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
110   if (!eventHandler) {
111     Printf("ERROR: Could not retrieve MC event handler");
112     return;
113   }
114
115   AliMCEvent* mcEvent = eventHandler->MCEvent();
116   if (!mcEvent) {
117      Printf("ERROR: Could not retrieve MC event");
118      return;
119   }
120
121   Int_t eventNumber = fESDEvent->GetEventNumberInFile();
122
123   fArray1Mu->Clear();
124   fArray2Mu->Clear();
125
126   TLorentzVector v; 
127
128   AliMUONRecoCheck *rc = new AliMUONRecoCheck(fESDEvent, eventHandler);
129
130   AliMUONVTrackStore* trackRefs = rc->TrackRefs(eventNumber);
131   AliMUONVTrackStore* recoTracks = rc->ReconstructedTracks(eventNumber);
132   AliStack* pstack = mcEvent->Stack();
133   
134   TIter next(recoTracks->CreateIterator());
135   AliMUONTrack* trackReco = NULL;
136   
137   Int_t nTrackReco = recoTracks->GetSize();
138   Int_t nTracksESD = fESDEvent->GetNumberOfMuonTracks();
139   //printf("Tracks in recoTracks (%d) and in ESD (%d).\n", nTrackReco, nTracksESD); 
140
141   if (nTrackReco != nTracksESD) printf ("Tracks in recoTracks (%d) and in ESD (%d) do not match!\n", nTrackReco, nTracksESD);
142   
143   Int_t nreftracks = 0;
144   Int_t itrRec = 0;
145   while ( (trackReco = static_cast<AliMUONTrack*>(next())) != NULL ) {
146     // assign parameters concerning the reconstructed tracks
147     AliMUONTrackLight muLight;
148     
149     muLight.FillFromESD(fESDEvent->GetMuonTrack(itrRec));
150     
151     // find the reference track and store further information   
152     TParticle *part = muLight.FindRefTrack(trackReco, trackRefs, pstack);
153     if (part) { 
154       v.SetPxPyPzE(part->Px(), part->Py(), part->Pz(), part->Energy());
155       muLight.SetPGen(v); 
156       muLight.FillMuonHistory(pstack, part);
157       //muLight.PrintInfo("A");
158       //store the referenced track in the muonArray:
159       TClonesArray &muons = *fArray1Mu;
160       new (muons[nreftracks++]) AliMUONTrackLight(muLight);
161     } 
162     
163     itrRec++;
164   }  // end reco track
165
166   // now loop over muon pairs to build dimuons
167   Int_t nmuons = fArray1Mu->GetEntriesFast(); 
168   Int_t ndimuons = 0; 
169   for(Int_t itrRec1 = 0; itrRec1 < nmuons-1; itrRec1++) { 
170     AliMUONTrackLight* mu1 = (AliMUONTrackLight*) fArray1Mu->At(itrRec1); 
171     for(Int_t itrRec2 = itrRec1+1; itrRec2 < nmuons; itrRec2++){
172       AliMUONTrackLight* mu2 = (AliMUONTrackLight*) fArray1Mu->At(itrRec2); 
173       AliMUONPairLight dimuLight;
174       dimuLight.SetMuons(*mu1, *mu2);
175       //dimuLight.PrintInfo("A");
176       TClonesArray &dimuons = *fArray2Mu;
177       new (dimuons[ndimuons++]) AliMUONPairLight(dimuLight);
178     }
179   }
180     
181   delete rc;
182
183   // Post output data.
184   if (nreftracks != 0) {
185     fTree->Fill();
186   }
187   
188   PostData(0, fTree);
189   
190   
191 }      
192
193 //________________________________________________________________________
194 void AliAnalysisTaskRecoCheck::Terminate(Option_t *) 
195 {
196   // the end
197
198 }