don't lie in the log!
[u/mrichter/AliRoot.git] / PWGPP / TPC / AliMCTrackingTestTask.cxx
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 #include <TList.h>
32 #include <TTree.h>
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>
42 #include "AliTrackComparison.h"
43 #include "AliMCTrackingTestTask.h"
44 #include "AliGenInfoMaker.h"
45 #include "AliHelix.h"
46 #include "AliTrackPointArray.h"
47 #include "AliESDCaloCluster.h"
48
49 //
50 #include "AliMCInfo.h"
51 #include "AliComparisonObject.h"
52 #include "AliESDRecInfo.h"
53 #include "AliTPCParamSR.h"
54 #include "AliTracker.h"
55 #include "AliTPCseed.h"
56
57 #include "AliCDBManager.h"
58 #include "AliGRPManager.h"
59 #include "AliGeomManager.h"
60
61 // STL includes
62 #include <iostream>
63
64 using namespace std;
65
66 ClassImp(AliMCTrackingTestTask)
67
68 #define USE_STREAMER 0
69 #define DEBUG 0
70
71 //________________________________________________________________________
72 AliMCTrackingTestTask::AliMCTrackingTestTask() : 
73   AliAnalysisTask(), 
74   fMCinfo(0),     //! MC event handler
75   fESD(0),
76   fCurrentRun(-1),
77   fDebugStreamer(0),
78   fStreamLevel(0),
79   fDebugLevel(0),
80   fDebugOutputPath(),
81   fOutList(NULL),
82   fPitList(NULL),
83   fCompList(NULL)
84 {
85   //
86   // Default constructor (should not be used)
87   //
88 }
89
90 AliMCTrackingTestTask::AliMCTrackingTestTask(const AliMCTrackingTestTask& /*info*/) : 
91   AliAnalysisTask(), 
92   fMCinfo(0),     //! MC event handler
93   fESD(0),
94   fCurrentRun(-1),
95   fDebugStreamer(0),
96   fStreamLevel(0),
97   fDebugLevel(),
98   fDebugOutputPath(),
99   fOutList(NULL),
100   fPitList(NULL),
101   fCompList(NULL)
102 {
103   //
104   // Default constructor 
105   //
106 }
107
108
109
110 //________________________________________________________________________
111 AliMCTrackingTestTask::AliMCTrackingTestTask(const char *name) : 
112   AliAnalysisTask(name, "AliMCTrackingTestTask"), 
113   fMCinfo(0),     //! MC event handler
114   fESD(0),
115   fCurrentRun(-1),
116   fDebugStreamer(0),
117   fStreamLevel(0),
118   fDebugLevel(0),
119   fDebugOutputPath(),
120   fOutList(NULL),
121   fPitList(NULL),
122   fCompList(NULL)
123 {
124   //
125   // Normal constructor
126   //
127   // Input slot #0 works with a TChain
128   DefineInput(0, TChain::Class());
129   DefineOutput(0, TList::Class());
130
131   // create the list for comparison objects
132   fCompList = new TList;
133   //
134   //
135 }
136
137 AliMCTrackingTestTask::~AliMCTrackingTestTask(){
138   //
139   //
140   //
141   if (fDebugLevel>0)  printf("AliMCTrackingTestTask::~AliMCTrackingTestTask\n");
142   if (fDebugStreamer) delete fDebugStreamer;
143   fDebugStreamer=0;
144   if (fOutList)    delete fOutList;   fOutList   = 0;
145   if (fCompList)   delete fCompList;  fCompList = 0; 
146 }
147
148
149 //_____________________________________________________________________________
150 Bool_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 }
163 //________________________________________________________________________
164 void 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();
190 }
191
192 //________________________________________________________________________
193 void AliMCTrackingTestTask::CreateOutputObjects() 
194 {
195   //
196   // Connect the output objects
197   //
198   if(fDebugLevel>3)
199     cout << "AnalysisTaskTPCCluster::CreateOutputObjects()" << endl;
200   
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);  
220 }
221
222
223 //________________________________________________________________________
224 void AliMCTrackingTestTask::Exec(Option_t *) {
225   //
226   // Execute analysis for current event 
227   //
228
229 #if DEBUG
230   printf("New event!\n");
231 #endif
232
233   if(fDebugLevel>3)
234     cout << "AliMCTrackingTestTask::Exec()" << endl;
235     
236
237
238   // If MC has been connected   
239   if (!fMCinfo){
240     cout << "Not MC info\n" << endl;
241   }else{
242     ProcessMCInfo();
243     // fMCinfo->Print();
244     //DumpInfo();
245   }
246
247   PostData(0, fOutList);
248 }      
249
250
251 //________________________________________________________________________
252 void 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;
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
271   return;
272 }
273
274
275 //________________________________________________________________________
276 TTreeSRedirector *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";
286
287   printf(" get debug streamer \n");
288
289   dsName.ReplaceAll(" ","");
290   fDebugStreamer = new TTreeSRedirector(dsName.Data());
291   return fDebugStreamer;
292 }
293
294 //________________________________________________________________________
295 AliExternalTrackParam * 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
301   if (!ref) return 0x0;
302   if (!part) return 0x0;
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
312 //________________________________________________________________________
313 Bool_t  AliMCTrackingTestTask::PropagateToPoint(AliExternalTrackParam *param, Double_t *xyz, Double_t mass, Float_t step){
314   // 
315   // Propagate track to point xyz using 
316   // AliTracker::PropagateToBxByBz functionality
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]);
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
329   return isOK;
330 }
331
332 //________________________________________________________________________
333 void  AliMCTrackingTestTask::ProcessMCInfo(){
334   //
335   //
336   //
337    //
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();
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
367
368   for (Int_t ipart=0;ipart<npart;ipart++){
369     Int_t status = fMCinfo->GetParticleAndTR(ipart, particle, trefs);
370     if (status<0 || !particle || !trefs) continue;
371     Int_t nref = trefs->GetEntries();
372     if (nref<5) continue;
373     FitTrackRefs(particle,trefs);
374
375     AliTrackReference * tpcIn=0;
376     AliTrackReference * tpcOut=0;
377     AliTrackReference * trdIn=0;
378     AliTrackReference * trdOut=0;
379     AliTrackReference * itsIn=0;
380     AliTrackReference * itsOut=0;
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
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       }      
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
453       if (ref->R()<rmin) rmin=ref->R();
454       if (ref->R()>rmax) rmax=ref->R();
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       }
499     }
500
501     // -----------------------------------------------
502     if (tpcIn && tpcOut) {
503       ProcessRefTracker(tpcIn,tpcOut,particle,1);
504       ProcessRefTracker(tpcIn,tpcOut,particle,3);
505     }
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
599   }
600
601   if (trefs) trefs->Clear("C");
602   //delete particle;
603   //delete tpcIn;
604
605 }
606
607 void AliMCTrackingTestTask::ProcessRefTracker(AliTrackReference* refIn,  AliTrackReference* refOut, TParticle*part,Int_t type){
608   //
609   // Test propagation from In to out
610   //
611
612 #if DEBUG
613   printf("ProcessRefTracker\n");
614 #endif
615
616   AliExternalTrackParam *param = 0;
617   AliExternalTrackParam *paramMC = 0;
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()};
622   Double_t mass = part->GetMass();
623   Double_t step=1;
624   //
625   param=MakeTrack(refIn,part);
626   paramMC=MakeTrack(refOut,part);
627   paramDebug=MakeTrack(refIn,part);
628
629   if (!param) return;
630   if (type!=3) PropagateToPoint(param,xyzOut, mass, step);
631   //
632 #if 0
633   /*
634     if (type==3) {
635     AliTPCseed seed;
636     seed.Set(param->GetX(),param->GetAlpha(),param->GetParameter(),param->GetCovariance());
637     Float_t alpha= TMath::ATan2(refIn->Y(),refIn->X());
638     seed.Rotate(alpha-seed.GetAlpha());
639     seed.SetMass(mass);
640     for (Float_t xlayer= seed.GetX(); xlayer<refOut->R(); xlayer+=step){
641       seed.PropagateTo(xlayer);
642     }
643     seed.PropagateTo(refOut->R());
644     param->Set(seed.GetX(),seed.GetAlpha(),seed.GetParameter(),seed.GetCovariance());
645   }
646 */
647 #endif
648 #if USE_STREAMER
649   TTreeSRedirector *pcstream = GetDebugStreamer();
650   TVectorD gpos(3);
651   TVectorD gmom(3);
652   Bool_t isOK=kTRUE;
653   isOK&=param->Rotate(paramMC->GetAlpha());
654   isOK&=param->PropagateTo(paramMC->GetX(),AliTracker::GetBz());
655   param->GetXYZ(gpos.GetMatrixArray());
656   param->GetPxPyPz(gmom.GetMatrixArray());
657   if (pcstream){
658     (*pcstream)<<"MC"<<
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
669       "\n";
670   }
671 #endif
672   delete param;
673   delete paramMC;
674   delete paramDebug;
675 }
676
677  
678 Bool_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);    
691 }
692
693
694 void  AliMCTrackingTestTask::FitTrackRefs(TParticle * part, TClonesArray * trefs){
695   //
696   //
697   //
698   //
699
700 #if DEBUG
701   printf("FitTrackRefs\n");
702 #endif
703
704   if (!trefs) return;
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;
721   Double_t covar[15];
722   for (Int_t icov=0; icov<15; icov++) covar[icov]=0;
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);
730   //AliTrackReference * refOut = (AliTrackReference*)trefs->At(iref1);
731   AliExternalTrackParam *paramPropagate= MakeTrack(refIn,part);
732   AliExternalTrackParam *paramUpdate   = MakeTrack(refIn,part);
733   paramUpdate->AddCovariance(covar);
734   //Double_t mass = part->GetMass();
735   //Double_t charge = part->GetPDG()->Charge()/3.;
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   }
763 #if USE_STREAMER
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    }
791 #endif
792   delete paramPropagate;
793   delete paramUpdate;
794 }
795
796
797
798
799 void 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
812 void 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 */