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