]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG1/TPC/AliMCTrackingTestTask.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
add calculation of epsilon using (1-x)+x*ncoll as a weight
[u/mrichter/AliRoot.git] / PWG1 / TPC / AliMCTrackingTestTask.cxx
CommitLineData
7cc34f08 1//
2// This class is the task to check the
3// Propagation and Update method used in the
4// 1. AliExternalTrackParam
5// 2. AliTracker
6//
7// Pure Monte-Carlo data used, not influence of detectors
8
9// Input - TParticle + Array of track references - (Points alogn track trajectories)
10// Output - Trees with track references - no histograms
11// MC tree - test for material budget correction
12// see function ProcessRefTracker
13// MCupdate tree - test for correctness of propagation and update
14// see function AliMCTrackingTestTask::FitTrackRefs
15//
16// Principle - Creates AliExternalTrackParam form 1 Track Refernece -
17// Propagate it to other
18// Magnetic field and the geometry has to be created before using it
19
20//
21//
22//
23
24//
25// ROOT includes
26#include <TChain.h>
27#include <TMath.h>
28#include <TVectorD.h>
29#include <TSystem.h>
30#include <TFile.h>
31
32// ALIROOT includes
33#include <TTreeStream.h>
34#include <AliAnalysisManager.h>
35#include <AliESDInputHandler.h>
36#include "AliStack.h"
37#include "AliMCEvent.h"
38#include "AliMCEventHandler.h"
39
40#include <AliESD.h>
41#include "AliMCTrackingTestTask.h"
42#include "AliGenInfoMaker.h"
43#include "AliHelix.h"
44
45//
46#include "AliMCInfo.h"
7cc34f08 47#include "AliESDRecInfo.h"
48#include "AliTPCParamSR.h"
49#include "AliTracker.h"
50#include "AliTPCseed.h"
51
52// STL includes
53#include <iostream>
54
55using namespace std;
56
57ClassImp(AliMCTrackingTestTask)
58
59//________________________________________________________________________
60AliMCTrackingTestTask::AliMCTrackingTestTask() :
61 AliAnalysisTask(),
62 fMCinfo(0), //! MC event handler
63 fESD(0),
64 fDebugStreamer(0),
65 fStreamLevel(0),
66 fDebugLevel(0),
67 fDebugOutputPath()
68{
69 //
70 // Default constructor (should not be used)
71 //
72}
73
74AliMCTrackingTestTask::AliMCTrackingTestTask(const AliMCTrackingTestTask& /*info*/) :
75 AliAnalysisTask(),
76 fMCinfo(0), //! MC event handler
77 fESD(0),
78 //
79 fDebugStreamer(0),
80 fStreamLevel(0),
81 fDebugLevel(),
82 fDebugOutputPath()
83{
84 //
85 // Default constructor
86 //
87}
88
89
90
91//________________________________________________________________________
92AliMCTrackingTestTask::AliMCTrackingTestTask(const char *name) :
93 AliAnalysisTask(name, "AliMCTrackingTestTask"),
94 fMCinfo(0), //! MC event handler
95 fESD(0),
96 fDebugStreamer(0),
97 fStreamLevel(0),
98 fDebugLevel(0),
99 fDebugOutputPath()
100{
101 //
102 // Normal constructor
103 //
104 // Input slot #0 works with a TChain
105 DefineInput(0, TChain::Class());
106 // Output slot #0 writes into a TList
107 DefineOutput(0, TObjArray::Class());
108 //
109 //
110}
111
112AliMCTrackingTestTask::~AliMCTrackingTestTask(){
113 //
114 //
115 //
116 if (fDebugLevel>0) printf("AliMCTrackingTestTask::~AliMCTrackingTestTask\n");
117 if (fDebugStreamer) delete fDebugStreamer;
118 fDebugStreamer=0;
119}
120
121
122//________________________________________________________________________
123void AliMCTrackingTestTask::ConnectInputData(Option_t *)
124{
125 //
126 // Connect the input data
127 //
128 if(fDebugLevel>3)
129 cout << "AnalysisTaskTPCCluster::ConnectInputData()" << endl;
130
131 TTree* tree=dynamic_cast<TTree*>(GetInputData(0));
132 if (!tree) {
133 //Printf("ERROR: Could not read chain from input slot 0");
134 }
135 else {
136 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
137 if (!esdH) {
138 //Printf("ERROR: Could not get ESDInputHandler");
139 }
140 else {
141 fESD = esdH->GetEvent();
142 //Printf("*** CONNECTED NEW EVENT ****");
143 }
144 }
145 AliMCEventHandler* mcinfo = (AliMCEventHandler*) (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
146 mcinfo->SetReadTR(kTRUE);
147
148 fMCinfo = mcinfo->MCEvent();
149
150
151}
152
153
154
155
156
157
158//________________________________________________________________________
159void AliMCTrackingTestTask::CreateOutputObjects()
160{
161 //
162 // Connect the output objects
163 //
164 if(fDebugLevel>3)
165 cout << "AnalysisTaskTPCCluster::CreateOutputObjects()" << endl;
166
167}
168
169
170//________________________________________________________________________
171void AliMCTrackingTestTask::Exec(Option_t *) {
172 //
173 // Execute analysis for current event
174 //
175
176 if(fDebugLevel>3)
177 cout << "AliMCTrackingTestTask::Exec()" << endl;
178
179
180 // If MC has been connected
181
182 if (!fMCinfo){
183 cout << "Not MC info\n" << endl;
184 }else{
185 ProcessMCInfo();
186 //mcinfo->Print();
187 //DumpInfo();
188 }
189 //
190}
191
192
193
194
195//________________________________________________________________________
196void AliMCTrackingTestTask::Terminate(Option_t *) {
197 //
198 // Terminate loop
199 //
200 if(fDebugLevel>3)
201 printf("AliMCTrackingTestTask: Terminate() \n");
202 //
203 if (fDebugLevel>0) printf("AliMCtrackingTestTask::Terminate\n");
204 if (fDebugStreamer) delete fDebugStreamer;
205 fDebugStreamer = 0;
206 return;
207}
208
209
210
211TTreeSRedirector *AliMCTrackingTestTask::GetDebugStreamer(){
212 //
213 // Get Debug streamer
214 // In case debug streamer not yet initialized and StreamLevel>0 create new one
215 //
216 if (fStreamLevel==0) return 0;
217 if (fDebugStreamer) return fDebugStreamer;
218 TString dsName;
219 dsName=GetName();
220 dsName+="Debug.root";
221 dsName.ReplaceAll(" ","");
222 fDebugStreamer = new TTreeSRedirector(dsName.Data());
223 return fDebugStreamer;
224}
225
226
227
228
229
230
231AliExternalTrackParam * AliMCTrackingTestTask::MakeTrack(const AliTrackReference* ref, TParticle*part)
232{
233 //
234 // Make track out of the track ref
235 // part - TParticle used to determine chargr
236 // the covariance matrix - equal 0 - starting from ideal MC position
237 Double_t xyz[3]={ref->X(),ref->Y(),ref->Z()};
238 Double_t pxyz[3]={ref->Px(),ref->Py(),ref->Pz()};
239 Double_t cv[21];
240 for (Int_t i=0; i<21;i++) cv[i]=0;
241 if (!part->GetPDG()) return 0;
242 AliExternalTrackParam * param = new AliExternalTrackParam(xyz,pxyz,cv,TMath::Nint(part->GetPDG()->Charge()/3.));
243 return param;
244}
245
246Bool_t AliMCTrackingTestTask::PropagateToPoint(AliExternalTrackParam *param, Double_t *xyz, Double_t mass, Float_t step){
247 //
248 // Propagate track to point xyz using
249 // AliTracker::PropagateTo functionality
250 //
251 // param - track parameters
252 // xyz - position to propagate
253 // mass - particle mass
254 // step - step to be used
255 Double_t radius=TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
256 AliTracker::PropagateTrackTo(param, radius+step, mass, step, kTRUE,0.99);
257 AliTracker::PropagateTrackTo(param, radius+0.5, mass, step*0.1, kTRUE,0.99);
258 Double_t sxyz[3]={0,0,0};
259 AliESDVertex vertex(xyz,sxyz);
260 Bool_t isOK = param->PropagateToDCA(&vertex,AliTracker::GetBz(),10);
261 return isOK;
262}
263
264
265void AliMCTrackingTestTask::ProcessMCInfo(){
266 //
267 //
268 //
269 //
270 TParticle * particle= new TParticle;
271 TClonesArray * trefs = new TClonesArray("AliTrackReference");
272 const Double_t kPcut=0.1;
273 //
274 //
275 // Process tracks
276 //
277 Int_t npart = fMCinfo->GetNumberOfTracks();
278 if (npart==0) return;
279 Double_t vertex[4]={0,0,0,0};
280 fMCinfo->GetParticleAndTR(0, particle, trefs);
281 if (particle){
282 vertex[0]=particle->Vx();
283 vertex[1]=particle->Vy();
284 vertex[2]=particle->Vz();
285 vertex[3]=particle->R();
286 }
287 //
288 //
289
290 for (Int_t ipart=0;ipart<npart;ipart++){
291 Int_t status = fMCinfo->GetParticleAndTR(ipart, particle, trefs);
292 if (status<0) continue;
293 if (!particle) continue;
294 if (!trefs) continue;
295 Int_t nref = trefs->GetEntries();
296 if (nref<5) continue;
297 FitTrackRefs(particle,trefs);
298 AliTrackReference * tpcIn=0;
299 AliTrackReference * tpcOut=0;
300 AliTrackReference * trdIn=0;
301 AliTrackReference * trdOut=0;
302 AliTrackReference * itsIn=0;
303 AliTrackReference * itsOut=0;
304 Double_t rmax=0;
305 Double_t rmin=1000;
306 for (Int_t iref=0;iref<nref;iref++){
307 AliTrackReference * ref = (AliTrackReference*)trefs->At(iref);
308 if (!ref) continue;
309
310 Float_t dir = ref->X()*ref->Px()+ref->Y()*ref->Py();
311
312 if (dir<0) break; // oposite direction - looping track - return back
313 if (ref->P()<kPcut) continue;
314 if (ref->R()<rmax) break;
315 //if (ref->R()<rmin) break;
316 //TPC
317 if (ref->DetectorId()==AliTrackReference::kTPC){
318 if (!tpcIn) {
319 tpcIn = ref;
320 }else{
321 if (ref->R()>tpcIn->R()) tpcOut = ref;
322 }
323 }
324 //ITS
325 if (ref->DetectorId()==AliTrackReference::kITS){
326 if (!itsIn) {
327 itsIn = ref;
328 }else{
329 if (ref->R()>itsIn->R()) itsOut = ref;
330 }
331 }
332 //TRD
333 if (ref->DetectorId()==AliTrackReference::kTRD){
334 if (!trdIn) {
335 trdIn = ref;
336 }else{
337 if (ref->R()>trdIn->R()) trdOut = ref;
338 }
339 }
340 if (ref->R()<rmin) rmin=ref->R();
341 if (ref->R()>rmax) rmax=ref->R();
342 }
343 if (tpcIn && tpcOut) {
344 ProcessRefTracker(tpcIn,tpcOut,particle,1);
345 ProcessRefTracker(tpcIn,tpcOut,particle,3);
346 }
347 if (itsIn && itsOut) ProcessRefTracker(itsIn,itsOut,particle,0);
348 if (trdIn && trdOut) ProcessRefTracker(trdIn,trdOut,particle,2);
349 }
350}
351
352
353
354void AliMCTrackingTestTask::ProcessRefTracker(AliTrackReference* refIn, AliTrackReference* refOut, TParticle*part,Int_t type){
355 //
356 // Test propagation from In to out
357 //
358 AliExternalTrackParam *param = 0;
359 AliExternalTrackParam *paramMC = 0;
360 Double_t xyzIn[3]={refIn->X(),refIn->Y(), refIn->Z()};
361 Double_t mass = part->GetMass();
362 Double_t step=1;
363 //
364 param=MakeTrack(refOut,part);
365 paramMC=MakeTrack(refOut,part);
366 if (!param) return;
367 if (type<3) PropagateToPoint(param,xyzIn, mass, step);
368 if (type==3) {
369 AliTPCseed seed;
370 seed.Set(param->GetX(),param->GetAlpha(),param->GetParameter(),param->GetCovariance());
371 Float_t alpha= TMath::ATan2(refIn->Y(),refIn->X());
26a35193 372 if(seed.Rotate(alpha-seed.GetAlpha())==kFALSE) return;
7cc34f08 373 seed.SetMass(mass);
374 for (Float_t xlayer= seed.GetX(); xlayer>refIn->R(); xlayer-=step){
375 seed.PropagateTo(xlayer);
376 }
377 seed.PropagateTo(refIn->R());
378 param->Set(seed.GetX(),seed.GetAlpha(),seed.GetParameter(),seed.GetCovariance());
379 }
380 TTreeSRedirector *pcstream = GetDebugStreamer();
381 TVectorD gpos(3);
382 TVectorD gmom(3);
383 param->GetXYZ(gpos.GetMatrixArray());
384 param->GetPxPyPz(gmom.GetMatrixArray());
385 if (pcstream){
386 (*pcstream)<<"MC"<<
387 "type="<<type<<
388 "step="<<step<<
389 "refIn.="<<refIn<<
390 "refOut.="<<refOut<<
391 "p.="<<part<<
392 "par.="<<param<<
393 "parMC.="<<paramMC<<
394 "gpos.="<<&gpos<<
395 "gmom.="<<&gmom<<
396 "\n";
397 }
398}
399
400
401void AliMCTrackingTestTask::FitTrackRefs(TParticle * part, TClonesArray * trefs){
402 //
403 //
404 //
405 //
406 const Int_t kMinRefs=6;
407 Int_t nrefs = trefs->GetEntries();
408 if (nrefs<kMinRefs) return; // we should have enough references
409 Int_t iref0 =-1;
410 Int_t iref1 =-1;
411
412 for (Int_t iref=0; iref<nrefs; iref++){
413 AliTrackReference * ref = (AliTrackReference*)trefs->At(iref);
414 if (!ref) continue;
415 Float_t dir = ref->X()*ref->Px()+ref->Y()*ref->Py();
416 if (dir<0) break;
417 if (ref->DetectorId()!=AliTrackReference::kTPC) continue;
418 if (iref0<0) iref0 = iref;
419 iref1 = iref;
420 }
421 if (iref1-iref0<kMinRefs) return;
c6d45264 422 Double_t covar[15];
423 for (Int_t icov=0; icov<15; icov++) covar[icov]=0;
7cc34f08 424 covar[0]=1;
425 covar[2]=1;
426 covar[5]=1;
427 covar[9]=1;
428 covar[14]=1;
429
430 AliTrackReference * refIn = (AliTrackReference*)trefs->At(iref0);
431 AliTrackReference * refOut = (AliTrackReference*)trefs->At(iref1);
432 AliExternalTrackParam *paramPropagate= MakeTrack(refIn,part);
433 AliExternalTrackParam *paramUpdate = MakeTrack(refIn,part);
434 paramUpdate->AddCovariance(covar);
435 Double_t mass = part->GetMass();
436 Double_t charge = part->GetPDG()->Charge()/3.;
437/*
438 Float_t alphaIn= TMath::ATan2(refIn->Y(),refIn->X());
439 Float_t radiusIn= refIn->R();
440 Float_t alphaOut= TMath::ATan2(refOut->Y(),refOut->X());
441 Float_t radiusOut= refOut->R();
442*/
443 Bool_t isOKP=kTRUE;
444 Bool_t isOKU=kTRUE;
445 AliMagF * field = (AliMagF*) TGeoGlobalMagField::Instance()->GetField();
446 for (Int_t iref = iref0; iref<=iref1; iref++){
447 AliTrackReference * ref = (AliTrackReference*)trefs->At(iref);
448 Float_t alphaC= TMath::ATan2(ref->Y(),ref->X());
449 Double_t pos[3] = {ref->X(), ref->Y(), ref->Z()};
450 Double_t mag[3];
451 field->Field(pos,mag);
452 isOKP&=paramPropagate->Rotate(alphaC);
453 isOKU&=paramUpdate->Rotate(alphaC);
454 for (Float_t xref= paramPropagate->GetX(); xref<ref->R(); xref++){
455 isOKP&=paramPropagate->PropagateTo(xref, mag[2]);
456 isOKU&=paramUpdate->PropagateTo(xref, mag[2]);
457 }
458 isOKP&=paramPropagate->PropagateTo(ref->R(), mag[2]);
459 isOKU&=paramUpdate->PropagateTo(ref->R(), mag[2]);
460 Double_t clpos[2] = {0, ref->Z()};
461 Double_t clcov[3] = { 0.005,0,0.005};
462 isOKU&= paramUpdate->Update(clpos, clcov);
463 }
464 TTreeSRedirector *pcstream = GetDebugStreamer();
465 if (pcstream){
466 TVectorD gposU(3);
467 TVectorD gmomU(3);
468 TVectorD gposP(3);
469 TVectorD gmomP(3);
470 paramUpdate->GetXYZ(gposU.GetMatrixArray());
471 paramUpdate->GetPxPyPz(gmomU.GetMatrixArray());
472 paramPropagate->GetXYZ(gposP.GetMatrixArray());
473 paramPropagate->GetPxPyPz(gmomP.GetMatrixArray());
474
475 (*pcstream)<<"MCupdate"<<
476 "isOKU="<<isOKU<<
477 "isOKP="<<isOKP<<
478 "m="<<mass<<
479 "q="<<charge<<
480 "part.="<<part<<
481 "refIn.="<<refIn<<
482 "refOut.="<<refOut<<
483 "pP.="<<paramPropagate<<
484 "pU.="<<paramUpdate<<
485 "gposU.="<<&gposU<<
486 "gmomU.="<<&gmomU<<
487 "gposP.="<<&gposP<<
488 "gmomP.="<<&gmomP<<
489 "\n";
490 }
491}
492
493
494
495
496void AliMCTrackingTestTask::FinishTaskOutput()
497{
498 //
499 // According description in AliAnalisysTask this method is call
500 // on the slaves before sending data
501 //
502 Terminate("slave");
503 gSystem->Exec("pwd");
504 RegisterDebugOutput();
505
506}
507
508
509void AliMCTrackingTestTask::RegisterDebugOutput(){
510 //
511 //
512 //
513 //
514 // store - copy debug output to the destination position
515 // currently ONLY for local copy
516 TString dsName;
517 dsName=GetName();
518 dsName+="Debug.root";
519 dsName.ReplaceAll(" ","");
520 TString dsName2=fDebugOutputPath.Data();
521 gSystem->MakeDirectory(dsName2.Data());
522 dsName2+=gSystem->HostName();
523 gSystem->MakeDirectory(dsName2.Data());
524 dsName2+="/";
525 dsName2+=gSystem->BaseName(gSystem->pwd());
526 dsName2+="/";
527 gSystem->MakeDirectory(dsName2.Data());
528 dsName2+=dsName;
529 AliInfo(Form("copy %s\t%s\n",dsName.Data(),dsName2.Data()));
530 printf("copy %s\t%s\n",dsName.Data(),dsName2.Data());
531 TFile::Cp(dsName.Data(),dsName2.Data());
532}
533
534
535/*
536 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
537 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
538 AliXRDPROOFtoolkit tool;
539 TChain * chain = tool.MakeChain("mctracking.txt","MC",0,100);
540 chain->Lookup();
541 //
542 //
543 chain->SetAlias("pdg","(p.fPdgCode)");
544 chain->SetAlias("dPRec","(refOut.P()-par.P())/refIn.P()");
545 chain->SetAlias("dPMC","(refOut.P()-refIn.P())/refIn->P()");
546 chain->SetAlias("dPtRec","(refOut.Pt()-par.Pt())/refIn.Pt()");
547 chain->SetAlias("dPtMC","(refOut.Pt()-refIn.Pt())/refIn->Pt()");
548
549
550 // ITS
551 chain->Draw("-sqrt(dPRec):-sqrt(dPMC)","abs(pdg)!=11&&type==0&&p.Pt()<0.5&&abs(p.R())<1&&abs(refOut.fZ)<220");
552 htemp->SetYTitle("#sqrt{#DeltaP_{rec}/P}");
553 htemp->SetXTitle("#sqrt{#DeltaP_{mc}/P}");
554 gPad->SaveAs("picLoss/dPcorr_ITS_step1.gif");
555 gPad->SaveAs("picLoss/dPcorr_ITS_step1.eps");
556 // TPC
557 chain->Draw("-sqrt(dPRec):-sqrt(dPMC)","abs(pdg)!=11&&type==1&&p.Pt()<0.5&&abs(p.R())<1&&abs(refOut.fZ)<220");
558 htemp->SetYTitle("#sqrt{#DeltaP_{rec}/P}");
559 htemp->SetXTitle("#sqrt{#DeltaP_{mc}/P}");
560 gPad->SaveAs("picLoss/dPcorr_TPC_step1.gif");
561 gPad->SaveAs("picLoss/dPcorr_TPC_step1.eps");
562 //
563 // TPC
564 chain->Draw("-sqrt(dPRec):-sqrt(dPMC)","abs(pdg)!=11&&type==3&&p.Pt()<0.5&&abs(p.R())<1&&abs(refOut.fZ)<220");
565 htemp->SetYTitle("#sqrt{#DeltaP_{rec}/P}");
566 htemp->SetXTitle("#sqrt{#DeltaP_{mc}/P}");
567 gPad->SaveAs("picLoss/dPcorr_TPCseed_step1.gif");
568 gPad->SaveAs("picLoss/dPcorr_TPCseed_step1.eps");
569
570
571 // TRD
572 chain->Draw("-sqrt(dPRec):-sqrt(dPMC)","abs(pdg)!=11&&type==2&&p.Pt()<0.5&&abs(p.R())<1&&abs(refOut.fZ)<220");
573 htemp->SetYTitle("#sqrt{#DeltaP_{rec}/P}");
574 htemp->SetXTitle("#sqrt{#DeltaP_{mc}/P}");
575 gPad->SaveAs("picLoss/dPcorr_TRD_step1.gif");
576 gPad->SaveAs("picLoss/dPcorr_TRD_step1.eps");
577
578 //
579 //
580 //
581 chain->Draw("(par.Pt()-refIn.Pt())/refIn.Pt()>>his(100,-0.02,0.02)","abs(pdg)!=11&&type==3&&p.Pt()<0.5&&abs(p.R())<1&&abs(refOut.fZ)<220");
582 his->SetXTitle("(P_{trec}-P_{tmc})/P_{tmc}");
583 gPad->SaveAs("picLoss/dPtcorr_TPCseed_step1_1D.eps");
584 gPad->SaveAs("picLoss/dPtcorr_TPCseed_step1_1D.gif");
585
586 chain->Draw("(par.P()-refIn.P())/refIn.P()>>his(100,-0.02,0.02)","abs(pdg)!=11&&type==3&&p.Pt()<0.5&&abs(p.R())<1&&abs(refOut.fZ)<220");
587 his->SetXTitle("(P_{rec}-P_{mc})/P_{mc}");
588 gPad->SaveAs("picLoss/dPcorr_TPCseed_step1_1D.eps");
589 gPad->SaveAs("picLoss/dPcorr_TPCseed_step1_1D.gif");
590*/