QA ref defaut storage setter in sim and rec
[u/mrichter/AliRoot.git] / PWG1 / AliGenInfoTask.cxx
CommitLineData
13db274d 1//
2// This class is the task for connecting together
3// MC information and the RC information
4//
5// The task is a wrapper over two components
6// AliGenInfoMaker
7// AliESDRecInfoMaker.h
8
9// ROOT includes
10#include <TChain.h>
11#include <TMath.h>
12
13// ALIROOT includes
9c3fd353 14#include <TTreeStream.h>
13db274d 15#include <AliAnalysisManager.h>
16#include <AliESDInputHandler.h>
17#include "AliStack.h"
18#include "AliMCEvent.h"
19#include "AliMCEventHandler.h"
20
21#include <AliESD.h>
22#include "AliGenInfoTask.h"
23#include "AliGenInfoMaker.h"
24#include "AliHelix.h"
25
9c3fd353 26//
27#include "AliMCInfo.h"
28#include "AliComparisonObject.h"
29#include "AliESDRecInfo.h"
30#include "AliTPCParamSR.h"
31
13db274d 32// STL includes
33#include <iostream>
34
35using namespace std;
36
37ClassImp(AliGenInfoTask)
38
39//________________________________________________________________________
40AliGenInfoTask::AliGenInfoTask() :
cd875161 41 AliAnalysisTask(),
9c3fd353 42 fMCinfo(0), //! MC event handler
43 fESD(0),
44 fCompList(0), //array of comparison objects
45 fGenTracksArray(0), //clones array with filtered particles
46 fGenKinkArray(0), //clones array with filtered Kinks
47 fGenV0Array(0), //clones array with filtered V0s
48 fRecTracksArray(0), //clones array with filtered particles
49 fDebugStreamer(0),
50 fStreamLevel(0),
51 fDebugLevel(0)
13db274d 52{
53 //
54 // Default constructor (should not be used)
55 //
13db274d 56}
57
cd875161 58AliGenInfoTask::AliGenInfoTask(const AliGenInfoTask& /*info*/) :
59 AliAnalysisTask(),
9c3fd353 60 fMCinfo(0), //! MC event handler
61 fESD(0),
62 fCompList(0),
63 fGenTracksArray(0), //clones array with filtered particles
64 fGenKinkArray(0), //clones array with filtered Kinks
65 fGenV0Array(0), //clones array with filtered V0s
66 fRecTracksArray(0), //clones array with filtered particles
67 //
68 fDebugStreamer(0),
69 fStreamLevel(0),
70 fDebugLevel()
cd875161 71{
72 //
73 // Default constructor
74 //
cd875161 75}
76
77
78
13db274d 79//________________________________________________________________________
80AliGenInfoTask::AliGenInfoTask(const char *name) :
81 AliAnalysisTask(name, "AliGenInfoTask"),
9c3fd353 82 fMCinfo(0), //! MC event handler
83 fESD(0),
84 fCompList(0),
85 fGenTracksArray(0), //clones array with filtered particles
86 fGenKinkArray(0), //clones array with filtered Kinks
87 fGenV0Array(0), //clones array with filtered V0s
88 fRecTracksArray(0), //clones array with filtered particles
89 fDebugStreamer(0),
90 fStreamLevel(0),
91 fDebugLevel(0)
13db274d 92{
93 //
94 // Normal constructor
95 //
13db274d 96 // Input slot #0 works with a TChain
97 DefineInput(0, TChain::Class());
98 // Output slot #0 writes into a TList
9c3fd353 99 DefineOutput(0, TObjArray::Class());
100 //
101 //
102 fCompList = new TObjArray;
103}
104
105AliGenInfoTask::~AliGenInfoTask(){
106 //
107 //
108 //
109 if (fDebugLevel>0) printf("AliGenInfoTask::~AliGenInfoTask\n");
110 if (fDebugStreamer) delete fDebugStreamer;
111 fDebugStreamer=0;
112 if(fCompList) delete fCompList;
113 fCompList =0;
13db274d 114}
115
9c3fd353 116
13db274d 117//________________________________________________________________________
118void AliGenInfoTask::ConnectInputData(Option_t *)
119{
120 //
121 // Connect the input data
122 //
9c3fd353 123 if(fDebugLevel>3)
13db274d 124 cout << "AnalysisTaskTPCCluster::ConnectInputData()" << endl;
125
9c3fd353 126 TTree* tree=dynamic_cast<TTree*>(GetInputData(0));
127 if (!tree) {
128 //Printf("ERROR: Could not read chain from input slot 0");
129 }
130 else {
131 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
132 if (!esdH) {
133 //Printf("ERROR: Could not get ESDInputHandler");
134 }
135 else {
136 fESD = esdH->GetEvent();
137 //Printf("*** CONNECTED NEW EVENT ****");
138 }
139 }
140 AliMCEventHandler* mcinfo = (AliMCEventHandler*) (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
141 mcinfo->SetReadTR(kTRUE);
142
143 fMCinfo = mcinfo->MCEvent();
13db274d 144
13db274d 145
9c3fd353 146}
147
148
149//_____________________________________________________________________________
150Bool_t AliGenInfoTask::AddComparisonObject(AliComparisonObject *pObj)
151{
152 // add comparison object to the list
153 if(pObj == 0) {
154 Printf("ERROR: Could not add comparison object");
155 return kFALSE;
13db274d 156 }
9c3fd353 157 // add object to the list
158 fCompList->AddLast(pObj);
159 return kTRUE;
13db274d 160}
161
9c3fd353 162
163
164
13db274d 165//________________________________________________________________________
166void AliGenInfoTask::CreateOutputObjects()
167{
168 //
169 // Connect the output objects
170 //
9c3fd353 171 if(fDebugLevel>3)
13db274d 172 cout << "AnalysisTaskTPCCluster::CreateOutputObjects()" << endl;
173
13db274d 174}
175
176
177//________________________________________________________________________
178void AliGenInfoTask::Exec(Option_t *) {
179 //
180 // Execute analysis for current event
13db274d 181 //
182
9c3fd353 183 if(fDebugLevel>3)
13db274d 184 cout << "AliGenInfoTask::Exec()" << endl;
9c3fd353 185
13db274d 186
13db274d 187 // If MC has been connected
9c3fd353 188 if (fGenTracksArray) fGenTracksArray->Delete();
189 if (fRecTracksArray) fRecTracksArray->Delete();
190
191 if (!fMCinfo){
13db274d 192 cout << "Not MC info\n" << endl;
9c3fd353 193 }else{
194 //mcinfo->Print();
195 ProcessMCInfo();
196 ProcessESDInfo();
197 DumpInfo();
198 ProcessComparison();
13db274d 199 }
9c3fd353 200 //
201}
13db274d 202
203
204
9c3fd353 205
206//________________________________________________________________________
207void AliGenInfoTask::Terminate(Option_t *) {
208 //
209 // Terminate loop
210 //
211 if(fDebugLevel>3)
212 printf("AliGenInfoTask: Terminate() \n");
213 //
214 if (fDebugLevel>0) printf("AliTPCcalibBase::Terminate\n");
215 if (fDebugStreamer) delete fDebugStreamer;
216 fDebugStreamer = 0;
217 return;
218}
219
220
221
222TTreeSRedirector *AliGenInfoTask::GetDebugStreamer(){
223 //
224 // Get Debug streamer
225 // In case debug streamer not yet initialized and StreamLevel>0 create new one
226 //
227 if (fStreamLevel==0) return 0;
228 if (fDebugStreamer) return fDebugStreamer;
229 TString dsName;
230 dsName=GetName();
231 dsName+="Debug.root";
232 dsName.ReplaceAll(" ","");
233 fDebugStreamer = new TTreeSRedirector(dsName.Data());
234 return fDebugStreamer;
235}
236
237
238
239AliMCInfo* AliGenInfoTask::GetTrack(Int_t index, Bool_t force){
240 //
241 // Get the MC info for given track
242 //
243 if (!fGenTracksArray) fGenTracksArray = new TClonesArray("AliMCInfo",1000);
244 if (index>fGenTracksArray->GetEntriesFast()) fGenTracksArray->Expand(index*2+10);
245 AliMCInfo * info = (AliMCInfo*)fGenTracksArray->At(index);
246 if (!force) return info;
247 if (!info){
248 info = new ((*fGenTracksArray)[index]) AliMCInfo;
13db274d 249 }
9c3fd353 250 return info;
251}
252
253AliESDRecInfo* AliGenInfoTask::GetRecTrack(Int_t index, Bool_t force){
254 //
255 // Get the MC info for given track
256 //
257 if (!fRecTracksArray) fRecTracksArray = new TClonesArray("AliESDRecInfo",1000);
258 if (index>fRecTracksArray->GetEntriesFast()) fRecTracksArray->Expand(index*2+10);
259 AliESDRecInfo * info = (AliESDRecInfo*)fRecTracksArray->At(index);
260 if (!force) return info;
261 if (!info){
262 info = new ((*fRecTracksArray)[index]) AliESDRecInfo;
13db274d 263 }
9c3fd353 264 return info;
265}
13db274d 266
9c3fd353 267
268
269
270void AliGenInfoTask::ProcessMCInfo(){
271 //
272 // Dump information from MC to the array
273 //
274 //
275 TParticle * particle= new TParticle;
276 TClonesArray * trefs = new TClonesArray("AliTrackReference");
13db274d 277 //
13db274d 278 //
9c3fd353 279 // Process tracks
280 //
281 Int_t npart = fMCinfo->GetNumberOfTracks();
282 if (npart==0) return;
283 Double_t vertex[4]={0,0,0,0};
284 fMCinfo->GetParticleAndTR(0, particle, trefs);
285 if (particle){
286 vertex[0]=particle->Vx();
287 vertex[1]=particle->Vy();
288 vertex[2]=particle->Vz();
289 vertex[3]=particle->R();
290 }
13db274d 291
9c3fd353 292 for (Int_t ipart=0;ipart<npart;ipart++){
293 Int_t status = fMCinfo->GetParticleAndTR(ipart, particle, trefs);
294 if (status<0) continue;
295 if (!particle) continue;
296 if (!trefs) continue;
297 if (!AcceptParticle(particle)) continue;
298 //if (trefs->GetEntries()<1) continue;
299 AliMCInfo * mcinfo = GetTrack(ipart,kTRUE);
300 mcinfo->Update(particle,trefs,vertex,ipart);
301 //
302 TTreeSRedirector *pcstream = GetDebugStreamer();
303 if (pcstream){
304 (*pcstream)<<"MC"<<
305 "p.="<<particle<<
306 "MC.="<<mcinfo<<
307 "\n";
308 }
309 }
13db274d 310}
311
9c3fd353 312void AliGenInfoTask::ProcessESDInfo(){
313 //
314 //
315 //
316 static AliTPCParamSR param;
317 //
318 //
319 if (!fESD) return;
320 Int_t ntracks = fESD->GetNumberOfTracks();
321 for (Int_t itrack=0; itrack<ntracks; itrack++){
322 AliESDtrack *track = fESD->GetTrack(itrack);
323 Int_t label = TMath::Abs(track->GetLabel());
324 AliMCInfo * mcinfo = GetTrack(label,kFALSE);
325 if (!mcinfo) continue;
326 AliESDRecInfo *recInfo= GetRecTrack(label,kTRUE);
327 recInfo->AddESDtrack(track,mcinfo);
328 recInfo->Update(mcinfo,&param,kTRUE);
329 }
330 //
331 //
332 //
333 Int_t ntracksMC = fMCinfo->GetNumberOfTracks();
334 for (Int_t imc=0; imc<ntracksMC; imc++){
335 AliMCInfo * mcinfo = GetTrack(imc,kFALSE);
336 if (!mcinfo) continue;
337 AliESDRecInfo *recInfo= GetRecTrack(imc,kFALSE);
338 if (recInfo) continue;
339 if (mcinfo->GetNTPCRef()<2) continue;
13db274d 340 //
13db274d 341 //
9c3fd353 342 for (Int_t itrack=0; itrack<ntracks; itrack++){
343 AliESDtrack *track = fESD->GetTrack(itrack);
344 Int_t label = TMath::Abs(track->GetLabel());
345 if (label!=mcinfo->GetLabel()) continue;
346
347 AliMCInfo * mcinfo2 = GetTrack(label,kFALSE);
348 if (!mcinfo2) continue;
349 AliESDRecInfo *recInfo= GetRecTrack(label,kTRUE);
350 recInfo->AddESDtrack(track,mcinfo2);
351 recInfo->Update(mcinfo2,&param,kTRUE);
352 }
353 }
354
355
356
357
358}
359
360
361void AliGenInfoTask::ProcessComparison(){
362 //
363 //
364 //
365 static AliESDRecInfo dummy;
366 Int_t npart = fMCinfo->GetNumberOfTracks();
367 for (Int_t ipart=0;ipart<npart;ipart++){
368 AliMCInfo * mcinfo = GetTrack(ipart,kFALSE);
369 if (!mcinfo) continue;
370 AliESDRecInfo *recInfo= GetRecTrack(ipart,kFALSE);
371 if (!recInfo) recInfo=&dummy;
372 //
373 for (Int_t icomp = 0; icomp<fCompList->GetEntries(); icomp++){
374 AliComparisonObject *pObj= (AliComparisonObject *)fCompList->At(icomp);
375 if (pObj){
376 pObj->Exec(mcinfo,recInfo);
377 }
378 }
379 }
380 PostData(0, fCompList);
381}
382
383void AliGenInfoTask::DumpInfo(){
384 //
385 //
386 //
387 static AliESDRecInfo dummy;
388 Int_t npart = fMCinfo->GetNumberOfTracks();
389 for (Int_t ipart=0;ipart<npart;ipart++){
390 AliMCInfo * mcinfo = GetTrack(ipart,kFALSE);
391 if (!mcinfo) continue;
392 AliESDRecInfo *recInfo= GetRecTrack(ipart,kFALSE);
393 if (!recInfo) recInfo=&dummy;
394 TTreeSRedirector *pcstream = GetDebugStreamer();
395 if (pcstream){
396 (*pcstream)<<"RC"<<
397 "MC.="<<mcinfo<<
398 "RC.="<<recInfo<<
399 "\n";
400 }
401 }
13db274d 402}
9c3fd353 403
404
405
406Bool_t AliGenInfoTask::AcceptParticle(TParticle *part){
407 //
408 /*
409 MC cuts
410 TCut cutPt("p.Pt()>0.1");
411 TCut cutZ("abs(p.Vz())<250");
412 TCut cutR("abs(p.R())<250");
413 //
414 */
415 //
416 //
417 if (part->Pt()<0.1) return kFALSE;
418 if (TMath::Abs(part->Vz())>250) return kFALSE;
419 if (part->R()>360) return kFALSE;
420 if (part->GetPDG()){
421 if (part->GetPDG()->Charge()==0) return kFALSE;
422 }
423 return kTRUE;
424}
425