]> 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
prtotect against missing QA projections in mass production
[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>
40d962ff 31#include <TList.h>
32#include <TTree.h>
7cc34f08 33// ALIROOT includes
34#include <TTreeStream.h>
35#include <AliAnalysisManager.h>
36#include <AliESDInputHandler.h>
37#include "AliStack.h"
38#include "AliMCEvent.h"
39#include "AliMCEventHandler.h"
40
41#include <AliESD.h>
40d962ff 42#include "AliTrackComparison.h"
7cc34f08 43#include "AliMCTrackingTestTask.h"
44#include "AliGenInfoMaker.h"
45#include "AliHelix.h"
40d962ff 46#include "AliTrackPointArray.h"
47#include "AliESDCaloCluster.h"
7cc34f08 48
49//
50#include "AliMCInfo.h"
40d962ff 51#include "AliComparisonObject.h"
7cc34f08 52#include "AliESDRecInfo.h"
53#include "AliTPCParamSR.h"
54#include "AliTracker.h"
55#include "AliTPCseed.h"
56
40d962ff 57#include "AliCDBManager.h"
58#include "AliGRPManager.h"
59#include "AliGeomManager.h"
60
7cc34f08 61// STL includes
62#include <iostream>
63
64using namespace std;
65
66ClassImp(AliMCTrackingTestTask)
67
40d962ff 68#define USE_STREAMER 0
69#define DEBUG 0
70
7cc34f08 71//________________________________________________________________________
72AliMCTrackingTestTask::AliMCTrackingTestTask() :
73 AliAnalysisTask(),
74 fMCinfo(0), //! MC event handler
75 fESD(0),
40d962ff 76 fCurrentRun(-1),
7cc34f08 77 fDebugStreamer(0),
78 fStreamLevel(0),
79 fDebugLevel(0),
40d962ff 80 fDebugOutputPath(),
81 fOutList(NULL),
82 fPitList(NULL),
83 fCompList(NULL)
7cc34f08 84{
85 //
86 // Default constructor (should not be used)
87 //
88}
89
90AliMCTrackingTestTask::AliMCTrackingTestTask(const AliMCTrackingTestTask& /*info*/) :
91 AliAnalysisTask(),
92 fMCinfo(0), //! MC event handler
93 fESD(0),
40d962ff 94 fCurrentRun(-1),
7cc34f08 95 fDebugStreamer(0),
96 fStreamLevel(0),
97 fDebugLevel(),
40d962ff 98 fDebugOutputPath(),
99 fOutList(NULL),
100 fPitList(NULL),
101 fCompList(NULL)
7cc34f08 102{
103 //
104 // Default constructor
105 //
106}
107
108
109
110//________________________________________________________________________
111AliMCTrackingTestTask::AliMCTrackingTestTask(const char *name) :
112 AliAnalysisTask(name, "AliMCTrackingTestTask"),
113 fMCinfo(0), //! MC event handler
114 fESD(0),
40d962ff 115 fCurrentRun(-1),
7cc34f08 116 fDebugStreamer(0),
117 fStreamLevel(0),
118 fDebugLevel(0),
40d962ff 119 fDebugOutputPath(),
120 fOutList(NULL),
121 fPitList(NULL),
122 fCompList(NULL)
7cc34f08 123{
124 //
125 // Normal constructor
126 //
127 // Input slot #0 works with a TChain
128 DefineInput(0, TChain::Class());
40d962ff 129 DefineOutput(0, TList::Class());
130
131 // create the list for comparison objects
132 fCompList = new TList;
7cc34f08 133 //
134 //
135}
136
137AliMCTrackingTestTask::~AliMCTrackingTestTask(){
138 //
139 //
140 //
141 if (fDebugLevel>0) printf("AliMCTrackingTestTask::~AliMCTrackingTestTask\n");
142 if (fDebugStreamer) delete fDebugStreamer;
143 fDebugStreamer=0;
40d962ff 144 if (fOutList) delete fOutList; fOutList = 0;
145 if (fCompList) delete fCompList; fCompList = 0;
7cc34f08 146}
147
148
40d962ff 149//_____________________________________________________________________________
150Bool_t AliMCTrackingTestTask::AddComparisonObject(AliTrackComparison *cObj)
151{
152 // add comparison object to the list
153 if(cObj == 0) {
154 Printf("ERROR: Could not add comparison object");
155 return kFALSE;
156 }
157
158 // add object to the list
159 fCompList->AddLast(cObj);
160
161 return kTRUE;
162}
7cc34f08 163//________________________________________________________________________
164void AliMCTrackingTestTask::ConnectInputData(Option_t *)
165{
166 //
167 // Connect the input data
168 //
169 if(fDebugLevel>3)
170 cout << "AnalysisTaskTPCCluster::ConnectInputData()" << endl;
171
172 TTree* tree=dynamic_cast<TTree*>(GetInputData(0));
173 if (!tree) {
174 //Printf("ERROR: Could not read chain from input slot 0");
175 }
176 else {
177 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
178 if (!esdH) {
179 //Printf("ERROR: Could not get ESDInputHandler");
180 }
181 else {
182 fESD = esdH->GetEvent();
183 //Printf("*** CONNECTED NEW EVENT ****");
184 }
185 }
186 AliMCEventHandler* mcinfo = (AliMCEventHandler*) (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
187 mcinfo->SetReadTR(kTRUE);
188
189 fMCinfo = mcinfo->MCEvent();
7cc34f08 190}
191
7cc34f08 192//________________________________________________________________________
193void AliMCTrackingTestTask::CreateOutputObjects()
194{
195 //
196 // Connect the output objects
197 //
198 if(fDebugLevel>3)
199 cout << "AnalysisTaskTPCCluster::CreateOutputObjects()" << endl;
200
40d962ff 201 if (!fOutList)
202 fOutList = new TList;
203 fOutList->SetOwner(kTRUE);
204 fPitList = fOutList->MakeIterator();
205
206 // create output list
207
208 // add comparison objects to the output
209 AliTrackComparison *cObj=0;
210 Int_t count=0;
211 TIterator *pitCompList = fCompList->MakeIterator();
212 pitCompList->Reset();
213 while(( cObj = (AliTrackComparison*)pitCompList->Next()) != NULL) {
214 fOutList->Add(cObj);
215 count++;
216 }
217 Printf("UserCreateOutputObjects(): Number of output comparison objects: %d \n", count);
218
219 PostData(0, fOutList);
7cc34f08 220}
221
222
223//________________________________________________________________________
224void AliMCTrackingTestTask::Exec(Option_t *) {
225 //
226 // Execute analysis for current event
227 //
228
40d962ff 229#if DEBUG
230 printf("New event!\n");
231#endif
232
7cc34f08 233 if(fDebugLevel>3)
234 cout << "AliMCTrackingTestTask::Exec()" << endl;
235
236
7cc34f08 237
40d962ff 238 // If MC has been connected
7cc34f08 239 if (!fMCinfo){
240 cout << "Not MC info\n" << endl;
241 }else{
242 ProcessMCInfo();
40d962ff 243 // fMCinfo->Print();
7cc34f08 244 //DumpInfo();
245 }
7cc34f08 246
40d962ff 247 PostData(0, fOutList);
248}
7cc34f08 249
250
251//________________________________________________________________________
252void AliMCTrackingTestTask::Terminate(Option_t *) {
253 //
254 // Terminate loop
255 //
256 if(fDebugLevel>3)
257 printf("AliMCTrackingTestTask: Terminate() \n");
258 //
259 if (fDebugLevel>0) printf("AliMCtrackingTestTask::Terminate\n");
260 if (fDebugStreamer) delete fDebugStreamer;
261 fDebugStreamer = 0;
40d962ff 262
263// AliTrackComparison *cObj=0;
264// TIterator *pitCompList = fCompList->MakeIterator();
265// pitCompList->Reset();
266// while(( cObj = (AliTrackComparison*)pitCompList->Next()) != NULL) {
267// for(Int_t i=0; i<6; i++)
268// cObj->MakeDistortionMap(438,i);
269// }
270
7cc34f08 271 return;
272}
273
274
40d962ff 275//________________________________________________________________________
7cc34f08 276TTreeSRedirector *AliMCTrackingTestTask::GetDebugStreamer(){
277 //
278 // Get Debug streamer
279 // In case debug streamer not yet initialized and StreamLevel>0 create new one
280 //
281 if (fStreamLevel==0) return 0;
282 if (fDebugStreamer) return fDebugStreamer;
283 TString dsName;
284 dsName=GetName();
285 dsName+="Debug.root";
40d962ff 286
287 printf(" get debug streamer \n");
288
7cc34f08 289 dsName.ReplaceAll(" ","");
290 fDebugStreamer = new TTreeSRedirector(dsName.Data());
291 return fDebugStreamer;
292}
293
40d962ff 294//________________________________________________________________________
7cc34f08 295AliExternalTrackParam * AliMCTrackingTestTask::MakeTrack(const AliTrackReference* ref, TParticle*part)
296{
297 //
298 // Make track out of the track ref
299 // part - TParticle used to determine chargr
300 // the covariance matrix - equal 0 - starting from ideal MC position
2a3fb9f6 301 if (!ref) return 0x0;
302 if (!part) return 0x0;
7cc34f08 303 Double_t xyz[3]={ref->X(),ref->Y(),ref->Z()};
304 Double_t pxyz[3]={ref->Px(),ref->Py(),ref->Pz()};
305 Double_t cv[21];
306 for (Int_t i=0; i<21;i++) cv[i]=0;
307 if (!part->GetPDG()) return 0;
308 AliExternalTrackParam * param = new AliExternalTrackParam(xyz,pxyz,cv,TMath::Nint(part->GetPDG()->Charge()/3.));
309 return param;
310}
311
40d962ff 312//________________________________________________________________________
7cc34f08 313Bool_t AliMCTrackingTestTask::PropagateToPoint(AliExternalTrackParam *param, Double_t *xyz, Double_t mass, Float_t step){
314 //
315 // Propagate track to point xyz using
40d962ff 316 // AliTracker::PropagateToBxByBz functionality
7cc34f08 317 //
318 // param - track parameters
319 // xyz - position to propagate
320 // mass - particle mass
321 // step - step to be used
322 Double_t radius=TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
40d962ff 323 Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
324
325 AliTracker::PropagateTrackToBxByBz(param, radius+step, mass, step, kTRUE,0.99,-1);
326 param->Rotate(alpha);
327 Bool_t isOK = param->PropagateTo(radius,AliTracker::GetBz());
328
7cc34f08 329 return isOK;
330}
331
40d962ff 332//________________________________________________________________________
7cc34f08 333void AliMCTrackingTestTask::ProcessMCInfo(){
334 //
335 //
336 //
337 //
40d962ff 338#if DEBUG
339 printf("ProcessMCInfo\n");
340#endif
341
342 const AliESDVertex *pVertex = fESD->GetPrimaryVertex();
343 if(!pVertex) AliError("No primary vertex found!\n");
344 Double_t vPos[3];
345 pVertex->GetXYZ(vPos);
346
347 TParticle * particle= new TParticle();
7cc34f08 348 TClonesArray * trefs = new TClonesArray("AliTrackReference");
349 const Double_t kPcut=0.1;
350 //
351 //
352 // Process tracks
353 //
354 Int_t npart = fMCinfo->GetNumberOfTracks();
355 if (npart==0) return;
356 Double_t vertex[4]={0,0,0,0};
357 fMCinfo->GetParticleAndTR(0, particle, trefs);
358 if (particle){
359 vertex[0]=particle->Vx();
360 vertex[1]=particle->Vy();
361 vertex[2]=particle->Vz();
362 vertex[3]=particle->R();
363 }
364 //
365 //
366
40d962ff 367
7cc34f08 368 for (Int_t ipart=0;ipart<npart;ipart++){
369 Int_t status = fMCinfo->GetParticleAndTR(ipart, particle, trefs);
40d962ff 370 if (status<0 || !particle || !trefs) continue;
7cc34f08 371 Int_t nref = trefs->GetEntries();
372 if (nref<5) continue;
373 FitTrackRefs(particle,trefs);
40d962ff 374
7cc34f08 375 AliTrackReference * tpcIn=0;
376 AliTrackReference * tpcOut=0;
377 AliTrackReference * trdIn=0;
378 AliTrackReference * trdOut=0;
379 AliTrackReference * itsIn=0;
380 AliTrackReference * itsOut=0;
40d962ff 381
382 AliTrackReference * tofIn=0;
383 AliTrackReference * tofOut=0;
384 AliTrackReference * hmpidIn=0;
385 AliTrackReference * hmpidOut=0;
386 AliTrackReference * emcalIn=0;
387 AliTrackReference * emcalOut=0;
388
7cc34f08 389 Double_t rmax=0;
390 Double_t rmin=1000;
391 for (Int_t iref=0;iref<nref;iref++){
392 AliTrackReference * ref = (AliTrackReference*)trefs->At(iref);
393 if (!ref) continue;
394
395 Float_t dir = ref->X()*ref->Px()+ref->Y()*ref->Py();
396
397 if (dir<0) break; // oposite direction - looping track - return back
398 if (ref->P()<kPcut) continue;
399 if (ref->R()<rmax) break;
400 //if (ref->R()<rmin) break;
401 //TPC
402 if (ref->DetectorId()==AliTrackReference::kTPC){
403 if (!tpcIn) {
404 tpcIn = ref;
405 }else{
406 if (ref->R()>tpcIn->R()) tpcOut = ref;
407 }
408 }
409 //ITS
410 if (ref->DetectorId()==AliTrackReference::kITS){
411 if (!itsIn) {
412 itsIn = ref;
413 }else{
414 if (ref->R()>itsIn->R()) itsOut = ref;
415 }
416 }
417 //TRD
418 if (ref->DetectorId()==AliTrackReference::kTRD){
419 if (!trdIn) {
420 trdIn = ref;
421 }else{
422 if (ref->R()>trdIn->R()) trdOut = ref;
423 }
424 }
40d962ff 425 //TOF
426 if (ref->DetectorId()==AliTrackReference::kTOF){
427 if (!tofIn) {
428 tofIn = ref;
429 }else{
430 if (ref->R()>tofIn->R()) tofOut = ref;
431 }
432 }
433
434 //HMPID
435 if (ref->DetectorId()==AliTrackReference::kHMPID){
436 if (!hmpidIn) {
437 hmpidIn = ref;
438 }else{
439 if (ref->R()>hmpidIn->R()) hmpidOut = ref;
440 }
441 }
442
443
444 //EMCAL
445 if (ref->DetectorId()==AliTrackReference::kEMCAL){
446 if (!emcalIn) {
447 emcalIn = ref;
448 }else{
449 if (ref->R()>emcalIn->R()) emcalOut = ref;
450 }
451 }
452
7cc34f08 453 if (ref->R()<rmin) rmin=ref->R();
454 if (ref->R()>rmax) rmax=ref->R();
40d962ff 455 } // end ref loop
456
457 // -----------------------------------------------
458 // check for gammas
459
460 // electron
461 if ( TMath::Abs(particle->GetPdgCode()) == 11 ) {
462
463 // findable
464 if (IsFindable(ipart,70)) {
465
466 // is from gamma conversion
467 Int_t motherId = particle->GetFirstMother();
468 if (motherId > 0) {
469 if (motherId < npart) {
470 TParticle* mother = (fMCinfo->Stack())->Particle(motherId);
471 if (mother && TMath::Abs(mother->GetPdgCode()) == 22) {
472 Int_t nDaughters = mother->GetNDaughters();
473
474 for (Int_t idx=0; idx<nDaughters; ++idx) {
475 Int_t daughterId = mother->GetDaughter(idx);
476 if ( daughterId == ipart || daughterId >= npart )
477 continue;
478
479 TParticle* daughter = (fMCinfo->Stack())->Particle(daughterId);
480 if (daughter && TMath::Abs(daughter->GetPdgCode()) == 11) {
481 //Bool_t findable = IsFindable(daughterId,70);
482#if USE_STREAMER
483 TTreeSRedirector *pcstream = GetDebugStreamer();
484 if (pcstream){
485 (*pcstream)<<"MCgamma"<<
486 "triggerP="<<particle<< // trigger electron
487 "motherP="<<mother<< // mother gamma
488 "daughterP="<<daughter<< // daughter electron
489 "isFindable="<<findable<< // 2nd is findable
490 "\n";
491 }
492#endif
493 }
494 }
495 }
496 }
497 }
498 }
7cc34f08 499 }
40d962ff 500
501 // -----------------------------------------------
7cc34f08 502 if (tpcIn && tpcOut) {
503 ProcessRefTracker(tpcIn,tpcOut,particle,1);
504 ProcessRefTracker(tpcIn,tpcOut,particle,3);
505 }
40d962ff 506 if (itsIn && itsOut) ProcessRefTracker(itsIn, itsOut, particle, 0);
507 if (trdIn && trdOut) ProcessRefTracker(trdIn, trdOut, particle, 2);
508
509 if (tpcOut && trdIn) ProcessRefTracker(tpcOut,trdIn, particle, 4);
510 if (tpcOut && tofIn) ProcessRefTracker(tpcOut,tofIn, particle, 5);
511 if (tpcOut && hmpidIn) ProcessRefTracker(tpcOut,hmpidIn,particle, 6);
512 if (tpcOut && emcalIn) ProcessRefTracker(tpcOut,emcalIn,particle, 7);
513
514 if (tpcOut && trdOut) ProcessRefTracker(tpcOut,trdOut, particle, 8);
515 if (trdIn && tofIn) ProcessRefTracker(trdIn, tofIn, particle, 9);
516 if (tofIn && tofOut) ProcessRefTracker(tofIn, tofOut, particle,10);
517
518
519 // -----------------------------------------------
520 //Test for new base tracking class
521
522 AliTrackComparison *cObj=0;
523 AliTrackPoint *point=new AliTrackPoint();
524 AliESDCaloCluster *cluster=0;
525 AliExternalTrackParam *tpcOuter=0;
526
527 Double_t eMax=0;
528 Int_t clsIndex=-1;
529
530 for(Int_t iCluster=0; iCluster<fESD->GetNumberOfCaloClusters(); iCluster++)
531 {
532 cluster = fESD->GetCaloCluster(iCluster);
533 if(!cluster->IsEMCAL()) continue;
534 if(cluster->GetLabel()==ipart && cluster->E()>eMax)
535 {
536 clsIndex=iCluster;
537 eMax=cluster->E();
538 }
539 }
540
541 if(clsIndex>-1)
542 {
543 if(cluster) cluster=0;
544 cluster = fESD->GetCaloCluster(clsIndex);
545 Float_t clusterPos[3];
546 cluster->GetPosition(clusterPos);
547 point->SetXYZ(clusterPos[0],clusterPos[1],clusterPos[2],0);
548 }
549
550 if(tpcOut)
551 tpcOuter = MakeTrack(tpcOut,particle);
552
553
554 Double_t mass = particle->GetMass();
555 Int_t charge = TMath::Nint(particle->GetPDG()->Charge()/3.);
556 fPitList->Reset();
557 while(( cObj = (AliTrackComparison *)fPitList->Next()) != NULL) {
558 TString objName(cObj->GetName());
559 if(!objName.CompareTo("TPCOutToEMCalInElecCls"))
560 {
561 if(TMath::Abs(particle->GetPdgCode())==11 && tpcOuter && point && cluster && emcalIn)
562 {
563 printf("\n\nTPCOutToEMCalInElecCls: ");
564 cout<<cObj->AddTracks(tpcOuter,point,mass,cluster->E(),vPos)<<endl;
565 }
566 }
567
568 if(!objName.CompareTo("TPCOutToEMCalInElec"))
569 {
570 if(TMath::Abs(particle->GetPdgCode())==11 && tpcOut && emcalIn)
571 {
572 printf("TPCOutToEMCalInElec: ");
573 cout<<cObj->AddTracks(tpcOut,emcalIn,mass,charge)<<endl;
574 }
575 }
576
577 if(!objName.CompareTo("TPCOutToEMCalInPion"))
578 {
579 if(TMath::Abs(particle->GetPdgCode())==211 && tpcOut && emcalIn)
580 cObj->AddTracks(tpcOut,emcalIn,mass,charge);
581 }
582
583 if(!objName.CompareTo("TPCOutToTOFIn"))
584 {
585 if(tpcOut && tofIn)
586 cObj->AddTracks(tpcOut,tofIn,mass,charge);
587 }
588
589 if(!objName.CompareTo("TPCOutToHMPIDIn"))
590 {
591 if(tpcOut && hmpidIn)
592 cObj->AddTracks(tpcOut,hmpidIn,mass,charge);
593 }
594 }
595 //End of the test for new base tracking class
596 // -----------------------------------------------
597 delete point;
598
7cc34f08 599 }
7cc34f08 600
2a3fb9f6 601 if (trefs) trefs->Clear("C");
40d962ff 602 //delete particle;
603 //delete tpcIn;
7cc34f08 604
40d962ff 605}
7cc34f08 606
607void AliMCTrackingTestTask::ProcessRefTracker(AliTrackReference* refIn, AliTrackReference* refOut, TParticle*part,Int_t type){
608 //
609 // Test propagation from In to out
610 //
40d962ff 611
612#if DEBUG
613 printf("ProcessRefTracker\n");
614#endif
615
7cc34f08 616 AliExternalTrackParam *param = 0;
617 AliExternalTrackParam *paramMC = 0;
40d962ff 618 AliExternalTrackParam *paramDebug = 0;
619
620 // Double_t xyzIn[3]={refIn->X(),refIn->Y(), refIn->Z()};
621 Double_t xyzOut[3]={refOut->X(),refOut->Y(), refOut->Z()};
7cc34f08 622 Double_t mass = part->GetMass();
623 Double_t step=1;
624 //
40d962ff 625 param=MakeTrack(refIn,part);
7cc34f08 626 paramMC=MakeTrack(refOut,part);
40d962ff 627 paramDebug=MakeTrack(refIn,part);
628
7cc34f08 629 if (!param) return;
40d962ff 630 if (type!=3) PropagateToPoint(param,xyzOut, mass, step);
631 //
632#if 0
633 /*
634 if (type==3) {
7cc34f08 635 AliTPCseed seed;
636 seed.Set(param->GetX(),param->GetAlpha(),param->GetParameter(),param->GetCovariance());
637 Float_t alpha= TMath::ATan2(refIn->Y(),refIn->X());
40d962ff 638 seed.Rotate(alpha-seed.GetAlpha());
7cc34f08 639 seed.SetMass(mass);
40d962ff 640 for (Float_t xlayer= seed.GetX(); xlayer<refOut->R(); xlayer+=step){
7cc34f08 641 seed.PropagateTo(xlayer);
642 }
40d962ff 643 seed.PropagateTo(refOut->R());
7cc34f08 644 param->Set(seed.GetX(),seed.GetAlpha(),seed.GetParameter(),seed.GetCovariance());
645 }
40d962ff 646*/
647#endif
648#if USE_STREAMER
7cc34f08 649 TTreeSRedirector *pcstream = GetDebugStreamer();
650 TVectorD gpos(3);
651 TVectorD gmom(3);
40d962ff 652 Bool_t isOK=kTRUE;
653 isOK&=param->Rotate(paramMC->GetAlpha());
654 isOK&=param->PropagateTo(paramMC->GetX(),AliTracker::GetBz());
7cc34f08 655 param->GetXYZ(gpos.GetMatrixArray());
656 param->GetPxPyPz(gmom.GetMatrixArray());
657 if (pcstream){
658 (*pcstream)<<"MC"<<
40d962ff 659 "isOK="<<isOK<<
660 "type="<<type<< // detector matching type
661 "step="<<step<< // propagation step length
662 "refIn.="<<refIn<< // starting track refernce
663 "refOut.="<<refOut<< // outer track reference
664 "p.="<<part<< // particle desription (TParticle)
665 "par.="<<param<< // AliExternalTrackParam create at starting point propagated to outer track ref radius
666 "parMC.="<<paramMC<< // AliExternalTrackParam created at the outer point
667 "gpos.="<<&gpos<< // global position
668 "gmom.="<<&gmom<< // global momenta
7cc34f08 669 "\n";
670 }
40d962ff 671#endif
672 delete param;
673 delete paramMC;
674 delete paramDebug;
675}
676
677
678Bool_t AliMCTrackingTestTask::IsFindable(Int_t label, Float_t minTrackLength ) {
679 //
680 // Find findable tracks
681 //
682
683 AliMCParticle *mcParticle = (AliMCParticle*) fMCinfo->GetTrack(label);
684 if(!mcParticle) return kFALSE;
685
686 Int_t counter;
687 Float_t tpcTrackLength = mcParticle->GetTPCTrackLength(AliTracker::GetBz(),0.05,counter,3.0);
688 //printf("tpcTrackLength %f \n", tpcTrackLength);
689
690 return (tpcTrackLength>minTrackLength);
7cc34f08 691}
692
693
694void AliMCTrackingTestTask::FitTrackRefs(TParticle * part, TClonesArray * trefs){
695 //
696 //
697 //
698 //
40d962ff 699
700#if DEBUG
701 printf("FitTrackRefs\n");
702#endif
703
2a3fb9f6 704 if (!trefs) return;
7cc34f08 705 const Int_t kMinRefs=6;
706 Int_t nrefs = trefs->GetEntries();
707 if (nrefs<kMinRefs) return; // we should have enough references
708 Int_t iref0 =-1;
709 Int_t iref1 =-1;
710
711 for (Int_t iref=0; iref<nrefs; iref++){
712 AliTrackReference * ref = (AliTrackReference*)trefs->At(iref);
713 if (!ref) continue;
714 Float_t dir = ref->X()*ref->Px()+ref->Y()*ref->Py();
715 if (dir<0) break;
716 if (ref->DetectorId()!=AliTrackReference::kTPC) continue;
717 if (iref0<0) iref0 = iref;
718 iref1 = iref;
719 }
720 if (iref1-iref0<kMinRefs) return;
c6d45264 721 Double_t covar[15];
722 for (Int_t icov=0; icov<15; icov++) covar[icov]=0;
7cc34f08 723 covar[0]=1;
724 covar[2]=1;
725 covar[5]=1;
726 covar[9]=1;
727 covar[14]=1;
728
729 AliTrackReference * refIn = (AliTrackReference*)trefs->At(iref0);
40d962ff 730 //AliTrackReference * refOut = (AliTrackReference*)trefs->At(iref1);
7cc34f08 731 AliExternalTrackParam *paramPropagate= MakeTrack(refIn,part);
732 AliExternalTrackParam *paramUpdate = MakeTrack(refIn,part);
733 paramUpdate->AddCovariance(covar);
40d962ff 734 //Double_t mass = part->GetMass();
735 //Double_t charge = part->GetPDG()->Charge()/3.;
7cc34f08 736/*
737 Float_t alphaIn= TMath::ATan2(refIn->Y(),refIn->X());
738 Float_t radiusIn= refIn->R();
739 Float_t alphaOut= TMath::ATan2(refOut->Y(),refOut->X());
740 Float_t radiusOut= refOut->R();
741*/
742 Bool_t isOKP=kTRUE;
743 Bool_t isOKU=kTRUE;
744 AliMagF * field = (AliMagF*) TGeoGlobalMagField::Instance()->GetField();
745 for (Int_t iref = iref0; iref<=iref1; iref++){
746 AliTrackReference * ref = (AliTrackReference*)trefs->At(iref);
747 Float_t alphaC= TMath::ATan2(ref->Y(),ref->X());
748 Double_t pos[3] = {ref->X(), ref->Y(), ref->Z()};
749 Double_t mag[3];
750 field->Field(pos,mag);
751 isOKP&=paramPropagate->Rotate(alphaC);
752 isOKU&=paramUpdate->Rotate(alphaC);
753 for (Float_t xref= paramPropagate->GetX(); xref<ref->R(); xref++){
754 isOKP&=paramPropagate->PropagateTo(xref, mag[2]);
755 isOKU&=paramUpdate->PropagateTo(xref, mag[2]);
756 }
757 isOKP&=paramPropagate->PropagateTo(ref->R(), mag[2]);
758 isOKU&=paramUpdate->PropagateTo(ref->R(), mag[2]);
759 Double_t clpos[2] = {0, ref->Z()};
760 Double_t clcov[3] = { 0.005,0,0.005};
761 isOKU&= paramUpdate->Update(clpos, clcov);
762 }
40d962ff 763#if USE_STREAMER
7cc34f08 764 TTreeSRedirector *pcstream = GetDebugStreamer();
765 if (pcstream){
766 TVectorD gposU(3);
767 TVectorD gmomU(3);
768 TVectorD gposP(3);
769 TVectorD gmomP(3);
770 paramUpdate->GetXYZ(gposU.GetMatrixArray());
771 paramUpdate->GetPxPyPz(gmomU.GetMatrixArray());
772 paramPropagate->GetXYZ(gposP.GetMatrixArray());
773 paramPropagate->GetPxPyPz(gmomP.GetMatrixArray());
774
775 (*pcstream)<<"MCupdate"<<
776 "isOKU="<<isOKU<<
777 "isOKP="<<isOKP<<
778 "m="<<mass<<
779 "q="<<charge<<
780 "part.="<<part<<
781 "refIn.="<<refIn<<
782 "refOut.="<<refOut<<
783 "pP.="<<paramPropagate<<
784 "pU.="<<paramUpdate<<
785 "gposU.="<<&gposU<<
786 "gmomU.="<<&gmomU<<
787 "gposP.="<<&gposP<<
788 "gmomP.="<<&gmomP<<
789 "\n";
790 }
40d962ff 791#endif
792 delete paramPropagate;
793 delete paramUpdate;
7cc34f08 794}
795
796
797
798
799void AliMCTrackingTestTask::FinishTaskOutput()
800{
801 //
802 // According description in AliAnalisysTask this method is call
803 // on the slaves before sending data
804 //
805 Terminate("slave");
806 gSystem->Exec("pwd");
807 RegisterDebugOutput();
808
809}
810
811
812void AliMCTrackingTestTask::RegisterDebugOutput(){
813 //
814 //
815 //
816 //
817 // store - copy debug output to the destination position
818 // currently ONLY for local copy
819 TString dsName;
820 dsName=GetName();
821 dsName+="Debug.root";
822 dsName.ReplaceAll(" ","");
823 TString dsName2=fDebugOutputPath.Data();
824 gSystem->MakeDirectory(dsName2.Data());
825 dsName2+=gSystem->HostName();
826 gSystem->MakeDirectory(dsName2.Data());
827 dsName2+="/";
828 dsName2+=gSystem->BaseName(gSystem->pwd());
829 dsName2+="/";
830 gSystem->MakeDirectory(dsName2.Data());
831 dsName2+=dsName;
832 AliInfo(Form("copy %s\t%s\n",dsName.Data(),dsName2.Data()));
833 printf("copy %s\t%s\n",dsName.Data(),dsName2.Data());
834 TFile::Cp(dsName.Data(),dsName2.Data());
835}
836
837
838/*
839 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
840 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
841 AliXRDPROOFtoolkit tool;
842 TChain * chain = tool.MakeChain("mctracking.txt","MC",0,100);
843 chain->Lookup();
844 //
845 //
846 chain->SetAlias("pdg","(p.fPdgCode)");
847 chain->SetAlias("dPRec","(refOut.P()-par.P())/refIn.P()");
848 chain->SetAlias("dPMC","(refOut.P()-refIn.P())/refIn->P()");
849 chain->SetAlias("dPtRec","(refOut.Pt()-par.Pt())/refIn.Pt()");
850 chain->SetAlias("dPtMC","(refOut.Pt()-refIn.Pt())/refIn->Pt()");
851
852
853 // ITS
854 chain->Draw("-sqrt(dPRec):-sqrt(dPMC)","abs(pdg)!=11&&type==0&&p.Pt()<0.5&&abs(p.R())<1&&abs(refOut.fZ)<220");
855 htemp->SetYTitle("#sqrt{#DeltaP_{rec}/P}");
856 htemp->SetXTitle("#sqrt{#DeltaP_{mc}/P}");
857 gPad->SaveAs("picLoss/dPcorr_ITS_step1.gif");
858 gPad->SaveAs("picLoss/dPcorr_ITS_step1.eps");
859 // TPC
860 chain->Draw("-sqrt(dPRec):-sqrt(dPMC)","abs(pdg)!=11&&type==1&&p.Pt()<0.5&&abs(p.R())<1&&abs(refOut.fZ)<220");
861 htemp->SetYTitle("#sqrt{#DeltaP_{rec}/P}");
862 htemp->SetXTitle("#sqrt{#DeltaP_{mc}/P}");
863 gPad->SaveAs("picLoss/dPcorr_TPC_step1.gif");
864 gPad->SaveAs("picLoss/dPcorr_TPC_step1.eps");
865 //
866 // TPC
867 chain->Draw("-sqrt(dPRec):-sqrt(dPMC)","abs(pdg)!=11&&type==3&&p.Pt()<0.5&&abs(p.R())<1&&abs(refOut.fZ)<220");
868 htemp->SetYTitle("#sqrt{#DeltaP_{rec}/P}");
869 htemp->SetXTitle("#sqrt{#DeltaP_{mc}/P}");
870 gPad->SaveAs("picLoss/dPcorr_TPCseed_step1.gif");
871 gPad->SaveAs("picLoss/dPcorr_TPCseed_step1.eps");
872
873
874 // TRD
875 chain->Draw("-sqrt(dPRec):-sqrt(dPMC)","abs(pdg)!=11&&type==2&&p.Pt()<0.5&&abs(p.R())<1&&abs(refOut.fZ)<220");
876 htemp->SetYTitle("#sqrt{#DeltaP_{rec}/P}");
877 htemp->SetXTitle("#sqrt{#DeltaP_{mc}/P}");
878 gPad->SaveAs("picLoss/dPcorr_TRD_step1.gif");
879 gPad->SaveAs("picLoss/dPcorr_TRD_step1.eps");
880
881 //
882 //
883 //
884 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");
885 his->SetXTitle("(P_{trec}-P_{tmc})/P_{tmc}");
886 gPad->SaveAs("picLoss/dPtcorr_TPCseed_step1_1D.eps");
887 gPad->SaveAs("picLoss/dPtcorr_TPCseed_step1_1D.gif");
888
889 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");
890 his->SetXTitle("(P_{rec}-P_{mc})/P_{mc}");
891 gPad->SaveAs("picLoss/dPcorr_TPCseed_step1_1D.eps");
892 gPad->SaveAs("picLoss/dPcorr_TPCseed_step1_1D.gif");
893*/