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 */