PWG1/AliMaterialBudget.cxx

   1
   2
   3
   4
   5 //
   6 //
   7
   8 // This class estiamtes the material budget of the inner detectors in ALICE based
   9 // on the "upper/lower track matching"-method of the ALICE TPC.
  10
  11 //
  12 //
  13 //
  14
  15 //
  16
  17
  18
  19
  20 // ROOT includes
  21 #include <TChain.h>
  22 #include <TMath.h>
  23 #include <TVectorD.h>
  24 #include <TSystem.h>
  25 #include <TFile.h>
  26 #include "TGeoGlobalMagField.h"
  27
  28 // ALIROOT includes
  29 #include <TTreeStream.h>
  30 #include <AliAnalysisManager.h>
  31 #include <AliESDInputHandler.h>
  32 #include "AliStack.h"
  33 #include "AliMCEvent.h"
  34 #include "AliMCEventHandler.h"
  35
  36 #include <AliESD.h>
  37 #include "AliMaterialBudget.h"
  38 #include "AliGenInfoMaker.h"
  39 #include "AliHelix.h"
  40
  41 //
  42 #include "AliMCInfo.h"
  43 #include "AliComparisonObject.h"
  44 #include "AliESDRecInfo.h"
  45 #include "AliTPCParamSR.h"
  46 #include "AliTracker.h"
  47 #include "AliTPCseed.h"
  48
  49 // STL includes
  50 #include <iostream>
  51
  52 using namespace std;
  53
  54 ClassImp(AliMaterialBudget)
  55
  56 //________________________________________________________________________
  57 AliMaterialBudget::AliMaterialBudget() :
  58   AliAnalysisTask(),
  59   fMCinfo(0),     //! MC event handler
  60   fESD(0),
  61   fDebugStreamer(0),
  62   fStreamLevel(0),
  63   fDebugLevel(0),
  64   fDebugOutputPath(),
  65   fListHist(0),
  66   fHistMult(0),
  67   fCutMaxD(5),        // maximal distance in rfi ditection
  68   fCutMaxDz(40),      // maximal distance in z ditection
  69   fCutTheta(0.03),    // maximal distan theta
  70   fCutMinDir(-0.99)   // direction vector products
  71 {
  72   //
  73   // Default constructor (should not be used)
  74   //
  75 }
  76
  77 AliMaterialBudget::AliMaterialBudget(const AliMaterialBudget& /*info*/) :
  78   AliAnalysisTask(),
  79   fMCinfo(0),     //! MC event handler
  80   fESD(0),
  81   //
  82   fDebugStreamer(0),
  83   fStreamLevel(0),
  84   fDebugLevel(),
  85   fDebugOutputPath(),
  86   fListHist(0),
  87   fHistMult(0),
  88   fCutMaxD(5),        // maximal distance in rfi ditection
  89   fCutMaxDz(40),      // maximal distance in z ditection
  90   fCutTheta(0.03),    // maximal distan theta
  91   fCutMinDir(-0.99)   // direction vector products
  92 {
  93   //
  94   // Default constructor
  95   //
  96 }
  97
  98
  99
 100 //________________________________________________________________________
 101 AliMaterialBudget::AliMaterialBudget(const char *name) :
 102   AliAnalysisTask(name, "AliMaterialBudget"),
 103   fMCinfo(0),     //! MC event handler
 104   fESD(0),
 105   fDebugStreamer(0),
 106   fStreamLevel(0),
 107   fDebugLevel(0),
 108   fDebugOutputPath(),
 109   fListHist(0),
 110   fHistMult(0),
 111   fCutMaxD(5),        // maximal distance in rfi ditection
 112   fCutMaxDz(40),      // maximal distance in z ditection
 113   fCutTheta(0.03),    // maximal distan theta
 114   fCutMinDir(-0.99)   // direction vector products
 115 {
 116   //
 117   // Normal constructor
 118   //
 119   // Input slot #0 works with a TChain
 120   DefineInput(0, TChain::Class());
 121   // Output slot #0 writes into a TList
 122   DefineOutput(0, TList::Class());
 123   //
 124   //
 125 }
 126
 127 AliMaterialBudget::~AliMaterialBudget(){
 128   //
 129   //
 130   //
 131   if (fDebugLevel>0)  printf("AliMaterialBudget::~AliMaterialBudget\n");
 132   if (fDebugStreamer) delete fDebugStreamer;
 133   fDebugStreamer=0;
 134 }
 135
 136
 137 //________________________________________________________________________
 138 void AliMaterialBudget::ConnectInputData(Option_t *)
 139 {
 140   //
 141   // Connect the input data
 142   //
 143   if(fDebugLevel>3)
 144     cout << "AnalysisTaskTPCCluster::ConnectInputData()" << endl;
 145
 146   TTree* tree=dynamic_cast<TTree*>(GetInputData(0));
 147   if (!tree) {
 148     //Printf("ERROR: Could not read chain from input slot 0");
 149   }
 150   else {
 151     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 152     if (!esdH) {
 153       //Printf("ERROR: Could not get ESDInputHandler");
 154     }
 155     else {
 156       esdH->SetActiveBranches("ESDfriend");
 157       fESD = esdH->GetEvent();
 158       //Printf("*** CONNECTED NEW EVENT ****");
 159     }
 160   }
 161   AliMCEventHandler* mcinfo = (AliMCEventHandler*) (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 162   mcinfo->SetReadTR(kTRUE);
 163
 164   fMCinfo = mcinfo->MCEvent();
 165
 166
 167 }
 168
 169
 170
 171
 172
 173
 174 //________________________________________________________________________
 175 void AliMaterialBudget::CreateOutputObjects()
 176 {
 177   //
 178   // Connect the output objects
 179   //
 180   if(fDebugLevel>3)
 181     cout << "AnalysisTaskTPCCluster::CreateOutputObjects()" << endl;
 182   //
 183   fListHist = new TList();
 184   fListHist->SetOwner();
 185   //
 186   fHistMult = new TH1F("HistMult", "Number of Tracks per Event; number of tracks per event; number of tracks",501,-0.5,500.5);
 187   fListHist->Add(fHistMult);
 188
 189
 190 }
 191
 192
 193 //________________________________________________________________________
 194 void AliMaterialBudget::Exec(Option_t *) {
 195   //
 196   // Execute analysis for current event
 197   //
 198
 199   if(fDebugLevel>3)
 200     cout << "AliMaterialBudget::Exec()" << endl;
 201
 202   fHistMult->Fill(fESD->GetNumberOfTracks());
 203   //FindPairs(fESD); // nearly everything takes place in find pairs...
 204
 205   // If MC has been connected
 206
 207   if (!fMCinfo){
 208     cout << "Not MC info\n" << endl;
 209   }else{
 210     ProcessMCInfo();
 211     //mcinfo->Print();
 212     //DumpInfo();
 213   }
 214   //
 215   PostData(0, fListHist);
 216 }
 217
 218
 219
 220
 221 //________________________________________________________________________
 222 void AliMaterialBudget::Terminate(Option_t *) {
 223     //
 224     // Terminate loop
 225     //
 226   if(fDebugLevel>3)
 227     printf("AliMaterialBudget: Terminate() \n");
 228   //
 229   if (fDebugLevel>0) printf("AliMCtrackingTestTask::Terminate\n");
 230   if (fDebugStreamer) delete fDebugStreamer;
 231   fDebugStreamer = 0;
 232   return;
 233 }
 234
 235
 236
 237 TTreeSRedirector *AliMaterialBudget::GetDebugStreamer(){
 238   //
 239   // Get Debug streamer
 240   // In case debug streamer not yet initialized and StreamLevel>0 create new one
 241   //
 242   if (fStreamLevel==0) return 0;
 243   if (fDebugStreamer) return fDebugStreamer;
 244   TString dsName;
 245   dsName=GetName();
 246   dsName+="Debug.root";
 247   dsName.ReplaceAll(" ","");
 248   fDebugStreamer = new TTreeSRedirector(dsName.Data());
 249   return fDebugStreamer;
 250 }
 251
 252
 253
 254
 255
 256
 257 AliExternalTrackParam * AliMaterialBudget::MakeTrack(const AliTrackReference* ref, TParticle*part)
 258 {
 259   //
 260   // Make track out of the track ref
 261   // part - TParticle used to determine chargr
 262   // the covariance matrix - equal 0 - starting from ideal MC position
 263   if (!part->GetPDG()) return 0;
 264   Double_t xyz[3]={ref->X(),ref->Y(),ref->Z()};
 265   Double_t pxyz[3]={ref->Px(),ref->Py(),ref->Pz()};
 266   Double_t charge = TMath::Nint(part->GetPDG()->Charge()/3.);
 267   if (ref->X()*ref->Px()+ref->Y()*ref->Py() <0){
 268     pxyz[0]*=-1;
 269     pxyz[1]*=-1;
 270     pxyz[2]*=-1;
 271     charge*=-1.;
 272   }
 273   Double_t cv[21];
 274   for (Int_t i=0; i<21;i++) cv[i]=0;
 275   AliExternalTrackParam * param = new AliExternalTrackParam(xyz,pxyz,cv,charge);
 276   return param;
 277 }
 278
 279 Bool_t  AliMaterialBudget::PropagateToPoint(AliExternalTrackParam *param, Double_t *xyz, Double_t mass, Float_t step){
 280   //
 281   // Propagate track to point xyz using
 282   // AliTracker::PropagateTo functionality
 283   //
 284   //  param - track parameters
 285   //  xyz   - position to propagate
 286   //  mass  - particle mass
 287   //  step  - step to be used
 288   Double_t radius=TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
 289   AliTracker::PropagateTrackToBxByBz(param, radius+step, mass, step, kTRUE,0.99);
 290   AliTracker::PropagateTrackToBxByBz(param, radius+0.5, mass, step*0.1, kTRUE,0.99);
 291   Double_t sxyz[3]={0,0,0};
 292   AliESDVertex vertex(xyz,sxyz);
 293   Bool_t isOK = param->PropagateToDCA(&vertex,AliTracker::GetBz(),10);
 294   return isOK;
 295 }
 296
 297
 298 void  AliMaterialBudget::ProcessMCInfo(){
 299   //
 300   //
 301   //
 302   //
 303   TParticle * particle= new TParticle;
 304   TClonesArray * trefs = new TClonesArray("AliTrackReference");
 305   const Double_t kPcut=0.05;
 306   const Double_t kMinDrITS = 2.;   // minimal distance between references
 307   const Double_t kMinDrTRD = 8.;   // minimal distance between references
 308   const Double_t kMinDrTOF = 10.;   // minimal distance between references
 309   //
 310   //
 311   // Process tracks
 312   //
 313   Int_t npart = fMCinfo->GetNumberOfTracks();
 314   if (npart==0) return;
 315   Double_t vertex[4]={0,0,0,0};
 316   fMCinfo->GetParticleAndTR(0, particle, trefs);
 317   if (particle){
 318     vertex[0]=particle->Vx();
 319     vertex[1]=particle->Vy();
 320     vertex[2]=particle->Vz();
 321     vertex[3]=particle->R();
 322   }
 323   //
 324   //
 325   AliTrackReference dummy,*pdummy= &dummy;
 326   AliExternalTrackParam epdummy,*pepdummy= &epdummy;
 327   Int_t nRefITS =0;
 328   Int_t nRefTPC =0;
 329   Int_t nRefTRD =0;
 330   Int_t nRefTOF =0;
 331   AliTrackReference * refITS0, *refITS1;
 332   AliTrackReference * refTPC0, *refTPC1;
 333   AliTrackReference * refTPCIn0, *refTPCIn1;
 334   AliTrackReference * refTRD0, *refTRD1;
 335   AliTrackReference * refTOF0, *refTOF1;
 336   AliTrackReference *refMinR;
 337   //
 338   for (Int_t ipart=0;ipart<npart;ipart++){
 339     Int_t status = fMCinfo->GetParticleAndTR(ipart, particle, trefs);
 340     AliMCParticle * pp = fMCinfo->GetTrack(ipart);
 341     if (!pp) continue;
 342     if (particle->P()<kPcut) continue;
 343     Double_t mass = particle->GetMass();
 344
 345     // RESET
 346     nRefITS =0;
 347     nRefTPC =0;
 348     nRefTRD =0;
 349     nRefTOF =0;
 350     refITS0=pdummy;    refITS1=pdummy;
 351     refTPC0=pdummy;    refTPC1=pdummy;
 352     refTPCIn0=pdummy;    refTPCIn1=pdummy;
 353     refTRD0=pdummy;    refTRD1=pdummy;
 354     refTOF0=pdummy;    refTOF1=pdummy;
 355     refMinR = pdummy;
 356     //
 357     Int_t nref = pp->GetNumberOfTrackReferences();
 358     if (nref==0) continue;
 359     for (Int_t iref = 0; iref < nref; iref++) {
 360       AliTrackReference *ref = pp->GetTrackReference(iref);
 361       if (ref->DetectorId()==AliTrackReference::kDisappeared) continue;
 362       //      if (ref.Px()*particle.Px()+ref.Py()*particle.Py()<0) break; // returning track
 363       //
 364       if (ref->DetectorId()==AliTrackReference::kITS){
 365         if (TMath::Abs(ref->R()-refITS1->R())>kMinDrITS) {
 366           refITS1 = ref;
 367           nRefITS++;
 368         }
 369         if (refITS0==pdummy) refITS0=ref;
 370       }
 371       if (ref->DetectorId()==AliTrackReference::kTPC){
 372         nRefTPC++;
 373         refTPC1 = ref;
 374         if (refTPC0==pdummy) refTPC0=ref;
 375       }
 376       if (ref->DetectorId()==AliTrackReference::kTRD){
 377         if (TMath::Abs(ref->R()-refTRD1->R())>kMinDrTRD) {
 378           refTRD1 = ref;
 379           nRefTRD++;
 380         }
 381         if (refTRD0==pdummy) refTRD0=ref;
 382       }
 383       if (ref->DetectorId()==AliTrackReference::kTOF){
 384         if (TMath::Abs(ref->X()-refTOF1->X()) + TMath::Abs(ref->Y()-refTOF1->Y()) + TMath::Abs(ref->Z()-refTOF1->Z())>kMinDrTOF) {
 385           refTOF1 = ref;
 386           nRefTOF++;
 387         }
 388         if (refTOF0==pdummy) refTOF0=ref;
 389       }
 390       //
 391       // "find inner track ref"
 392       if (ref->DetectorId()==AliTrackReference::kTPC){
 393         if (ref->Px()*ref->X()+ref->Py()*ref->Y()<0){
 394           //  track in
 395           if (refTPCIn0 == pdummy) refTPCIn0=ref;
 396           if (refTPCIn0 != pdummy &&  refTPCIn0->R()>ref->R())
 397             refTPCIn0=ref;
 398         }
 399         if (ref->Px()*ref->X()+ref->Py()*ref->Y()>0){
 400           //  track in
 401           if (refTPCIn1 == pdummy) refTPCIn1=ref;
 402           if (refTPCIn1 != pdummy &&  refTPCIn1->R()>ref->R())
 403             refTPCIn1=ref;
 404         }
 405       }
 406
 407
 408       if (refMinR==pdummy && ref->P()>0  ){
 409         refMinR=ref;
 410       }
 411       if (refMinR->R()>ref->R() && ref->P()>0 ) refMinR=ref;
 412
 413     }
 414     //
 415     AliExternalTrackParam * trackMC = pepdummy;
 416     //track0->GetDZ(0,0,0,bz,dvertex0)
 417     Float_t dist[2]={0,0};
 418     AliMagF* field = (AliMagF*) TGeoGlobalMagField::Instance()->GetField();
 419     Double_t esdfield= fESD->GetMagneticField();
 420     Double_t xyz[3]={0,0,0};
 421     Double_t bxyz[3]={0,0,0};
 422     field->Field(xyz,bxyz);
 423     if (refMinR->P()>0) {
 424       trackMC = MakeTrack(refMinR,particle);
 425       trackMC->GetDZ(0,0,0,bxyz[2],dist);
 426     }
 427     Double_t alphaTOF0 = TMath::ATan2(refTOF0->Y(),refTOF0->X());
 428     Double_t alphaTOF1 = TMath::ATan2(refTOF1->Y(),refTOF1->X());
 429     Int_t dsecTOF   = TMath::Nint(180*(alphaTOF0-alphaTOF1)/(TMath::Pi()*20.)-0.5);
 430     //
 431     // make the two different TPC tracks and propagate them to their DCA
 432     //
 433     Double_t dP = 0;
 434     Bool_t statusProp = kFALSE;
 435     Double_t dY = 0;
 436     Double_t dZ = 0;
 437     AliExternalTrackParam * track0  = pepdummy;
 438     AliExternalTrackParam * track1  = pepdummy;
 439     AliExternalTrackParam * otrack0 = pepdummy;
 440     AliExternalTrackParam * otrack1 = pepdummy;
 441     if (refTPCIn0!=pdummy && refTPCIn1!=pdummy) {
 442       track0 = MakeTrack(refTPCIn0,particle);
 443       track1 = MakeTrack(refTPCIn1,particle);
 444       otrack0 = MakeTrack(refTPCIn0,particle);
 445       otrack1 = MakeTrack(refTPCIn1,particle);
 446       dP = track0->P() - track1->P(); // momentum loss
 447       statusProp = AliMaterialBudget::PropagateCosmicToDCA(track0,track1,mass);
 448       if (statusProp) {
 449         dY = track0->GetY() - track1->GetY();
 450         dZ = track0->GetZ() - track1->GetZ();
 451       }
 452     }
 453     //
 454     TTreeSRedirector *pcstream = GetDebugStreamer();
 455     if (pcstream){
 456       char name[100];
 457       for (Int_t id=0; id<3;id++){
 458
 459         //      if (id==0) sprintf(name,"mcAll"); // all tracks: inconvenient to cut if on is only interest in tracks which reach the TPC
 460         if (id==0) continue; // require TPC
 461         if (id==1) sprintf(name,"mcITS");
 462         if (id==2) sprintf(name,"mcTPC");
 463         if (id==1&& nRefITS==0) continue;
 464         if (id==2&& nRefTPC==0) continue;
 465
 466         (*pcstream)<<name<<
 467           "ipart="<<ipart<<
 468           "p.="<<particle<<
 469           "mass="<<mass<<
 470           "tbfield="<<bxyz[2]<<
 471           "esdbfield="<<esdfield<<
 472           // counter
 473           "nref="<<nref<<
 474           "nRefITS="<<nRefITS<<
 475           "nRefTPC="<<nRefTPC<<
 476           "nRefTRD="<<nRefTRD<<
 477           "nRefTOF="<<nRefTOF<<
 478           //references
 479           "refMinR.="<<refMinR<<
 480           "refITS0.="<<refITS0<<
 481           "refITS1.="<<refITS1<<
 482           "refTPC0.="<<refTPC0<<
 483           "refTPC1.="<<refTPC1<<
 484           "refTPCIn0.="<<refTPCIn0<<
 485           "refTPCIn1.="<<refTPCIn1<<
 486           "refTRD0.="<<refTRD0<<
 487           "refTRD1.="<<refTRD1<<
 488           "refTOF0.="<<refTOF0<<
 489           "refTOF1.="<<refTOF1<<
 490           //trigger variables
 491           "dsecTOF="<<dsecTOF<<    // delta TOF sectors
 492           "aTOF0="<<alphaTOF0<<
 493           "aTOF1="<<alphaTOF1<<
 494           //
 495           // track
 496           "dr="<<dist[0]<<
 497           "dz="<<dist[1]<<
 498           //
 499           // "two" TPC tracks
 500           "status="<<statusProp<<
 501           "dP="<<dP<<
 502           "otrack0.="<<otrack0<<
 503           "otrack1.="<<otrack1<<
 504           "track0.="<<track0<<
 505           "track1.="<<track1<<
 506           "dY="<<dY<<
 507           "dZ="<<dZ<<
 508           "\n";
 509       }
 510     }
 511   }
 512 }
 513
 514
 515
 516 void AliMaterialBudget::ProcessRefTracker(AliTrackReference* refIn,  AliTrackReference* refOut, TParticle*part,Int_t type){
 517   //
 518   // Test propagation from In to out
 519   //
 520   AliExternalTrackParam *param = 0;
 521   AliExternalTrackParam *paramMC = 0;
 522   Double_t xyzIn[3]={refIn->X(),refIn->Y(), refIn->Z()};
 523   Double_t mass = part->GetMass();
 524   Double_t step=1;
 525   //
 526   param=MakeTrack(refOut,part);
 527   paramMC=MakeTrack(refOut,part);
 528   if (!param) return;
 529   if (type<3) PropagateToPoint(param,xyzIn, mass, step);
 530   if (type==3) {
 531     AliTPCseed seed;
 532     seed.Set(param->GetX(),param->GetAlpha(),param->GetParameter(),param->GetCovariance());
 533     Float_t alpha= TMath::ATan2(refIn->Y(),refIn->X());
 534     seed.Rotate(alpha-seed.GetAlpha());
 535     seed.SetMass(mass);
 536     for (Float_t xlayer= seed.GetX(); xlayer>refIn->R(); xlayer-=step){
 537       seed.PropagateTo(xlayer);
 538     }
 539     seed.PropagateTo(refIn->R());
 540     param->Set(seed.GetX(),seed.GetAlpha(),seed.GetParameter(),seed.GetCovariance());
 541   }
 542   TTreeSRedirector *pcstream = GetDebugStreamer();
 543   TVectorD gpos(3);
 544   TVectorD gmom(3);
 545   param->GetXYZ(gpos.GetMatrixArray());
 546   param->GetPxPyPz(gmom.GetMatrixArray());
 547   if (pcstream){
 548     (*pcstream)<<"MC"<<
 549       "type="<<type<<
 550       "step="<<step<<
 551       "refIn.="<<refIn<<
 552       "refOut.="<<refOut<<
 553       "p.="<<part<<
 554       "par.="<<param<<
 555       "parMC.="<<paramMC<<
 556       "gpos.="<<&gpos<<
 557       "gmom.="<<&gmom<<
 558       "\n";
 559   }
 560 }
 561
 562
 563 void  AliMaterialBudget::FitTrackRefs(TParticle * part, TClonesArray * trefs){
 564   //
 565   //
 566   //
 567   //
 568   const Int_t kMinRefs=6;
 569   Int_t nrefs = trefs->GetEntries();
 570   if (nrefs<kMinRefs) return; // we should have enough references
 571   Int_t iref0 =-1;
 572   Int_t iref1 =-1;
 573
 574   for (Int_t iref=0; iref<nrefs; iref++){
 575     AliTrackReference * ref = (AliTrackReference*)trefs->At(iref);
 576     if (!ref) continue;
 577     Float_t dir = ref->X()*ref->Px()+ref->Y()*ref->Py();
 578     if (dir<0) break;
 579     if (ref->DetectorId()!=AliTrackReference::kTPC) continue;
 580     if (iref0<0) iref0 = iref;
 581     iref1 = iref;
 582   }
 583   if (iref1-iref0<kMinRefs) return;
 584   Double_t covar[14];
 585   for (Int_t icov=0; icov<14; icov++) covar[icov]=0;
 586   covar[0]=1;
 587   covar[2]=1;
 588   covar[5]=1;
 589   covar[9]=1;
 590   covar[14]=1;
 591
 592   AliTrackReference * refIn = (AliTrackReference*)trefs->At(iref0);
 593   AliTrackReference * refOut = (AliTrackReference*)trefs->At(iref1);
 594   AliExternalTrackParam *paramPropagate= MakeTrack(refIn,part);
 595   AliExternalTrackParam *paramUpdate   = MakeTrack(refIn,part);
 596   paramUpdate->AddCovariance(covar);
 597   Double_t mass = part->GetMass();
 598   Double_t charge = part->GetPDG()->Charge()/3.;
 599 /*
 600   Float_t alphaIn= TMath::ATan2(refIn->Y(),refIn->X());
 601   Float_t radiusIn= refIn->R();
 602   Float_t alphaOut= TMath::ATan2(refOut->Y(),refOut->X());
 603   Float_t radiusOut= refOut->R();
 604 */
 605   Bool_t isOKP=kTRUE;
 606   Bool_t isOKU=kTRUE;
 607   AliMagF * field = (AliMagF*) TGeoGlobalMagField::Instance()->GetField();
 608   for (Int_t iref = iref0; iref<=iref1; iref++){
 609     AliTrackReference * ref = (AliTrackReference*)trefs->At(iref);
 610     Float_t alphaC= TMath::ATan2(ref->Y(),ref->X());
 611     Double_t pos[3] = {ref->X(), ref->Y(), ref->Z()};
 612     Double_t mag[3];
 613     field->Field(pos,mag);
 614     isOKP&=paramPropagate->Rotate(alphaC);
 615     isOKU&=paramUpdate->Rotate(alphaC);
 616     for (Float_t xref= paramPropagate->GetX(); xref<ref->R(); xref++){
 617       isOKP&=paramPropagate->PropagateTo(xref, mag[2]);
 618       isOKU&=paramUpdate->PropagateTo(xref, mag[2]);
 619     }
 620     isOKP&=paramPropagate->PropagateTo(ref->R(), mag[2]);
 621     isOKU&=paramUpdate->PropagateTo(ref->R(), mag[2]);
 622     Double_t clpos[2] = {0, ref->Z()};
 623     Double_t clcov[3] = { 0.005,0,0.005};
 624     isOKU&= paramUpdate->Update(clpos, clcov);
 625   }
 626   TTreeSRedirector *pcstream = GetDebugStreamer();
 627   if (pcstream){
 628     TVectorD gposU(3);
 629     TVectorD gmomU(3);
 630     TVectorD gposP(3);
 631     TVectorD gmomP(3);
 632     paramUpdate->GetXYZ(gposU.GetMatrixArray());
 633     paramUpdate->GetPxPyPz(gmomU.GetMatrixArray());
 634     paramPropagate->GetXYZ(gposP.GetMatrixArray());
 635     paramPropagate->GetPxPyPz(gmomP.GetMatrixArray());
 636
 637      (*pcstream)<<"MCupdate"<<
 638        "isOKU="<<isOKU<<
 639        "isOKP="<<isOKP<<
 640        "m="<<mass<<
 641        "q="<<charge<<
 642        "part.="<<part<<
 643        "refIn.="<<refIn<<
 644        "refOut.="<<refOut<<
 645        "pP.="<<paramPropagate<<
 646        "pU.="<<paramUpdate<<
 647        "gposU.="<<&gposU<<
 648        "gmomU.="<<&gmomU<<
 649        "gposP.="<<&gposP<<
 650        "gmomP.="<<&gmomP<<
 651        "\n";
 652    }
 653 }
 654
 655
 656
 657 void AliMaterialBudget::FindPairs(AliESDEvent * event) {
 658   //
 659   // This function matches the cosmic tracks and calculates the energy loss in the material.
 660   // If accessible the "true" energy loss is determined with the MC track references.
 661   //
 662   //
 663   // Find cosmic pairs
 664   //
 665   // Track0 is choosen in upper TPC part
 666   // Track1 is choosen in lower TPC part
 667   //
 668   if(fDebugLevel>3)
 669     cout << "AliMaterialBudget::FindPairs()" << endl;
 670
 671
 672   AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
 673   Int_t ntracks=event->GetNumberOfTracks();
 674   TObjArray  tpcSeeds(ntracks);
 675   if (ntracks==0) return;
 676   Double_t vtxx[3]={0,0,0};
 677   Double_t svtxx[3]={0.000001,0.000001,100.};
 678   AliESDVertex vtx(vtxx,svtxx);
 679   //
 680   //track loop
 681   //
 682   for (Int_t i=0;i<ntracks;++i) {
 683    AliESDtrack *track = event->GetTrack(i);
 684      const AliExternalTrackParam * trackIn = track->GetInnerParam();
 685    const AliExternalTrackParam * trackOut = track->GetOuterParam();
 686    if (!trackIn) continue;
 687    if (!trackOut) continue;
 688    AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
 689    TObject *calibObject;
 690    AliTPCseed *seed = 0;
 691    for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
 692      if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
 693    }
 694    if (seed) tpcSeeds.AddAt(seed,i);
 695   }
 696
 697   if (ntracks<2) return;
 698   //
 699   // Find pairs
 700   //
 701   for (Int_t i=0;i<ntracks;++i) {
 702     AliESDtrack *track0 = event->GetTrack(i);
 703     // track0 - choosen upper part
 704     if (!track0) continue;
 705     if (!track0->GetOuterParam()) continue;
 706     if (track0->GetOuterParam()->GetAlpha()<0) continue;
 707     Double_t dir0[3];
 708     track0->GetDirection(dir0);
 709     for (Int_t j=0;j<ntracks;++j) {
 710       if (i==j) continue;
 711       AliESDtrack *track1 = event->GetTrack(j);
 712       //track 1 lower part
 713       if (!track1) continue;
 714       if (!track1->GetOuterParam()) continue;
 715       if (track1->GetOuterParam()->GetAlpha()>0) continue;
 716       //
 717       Double_t dir1[3];
 718       track1->GetDirection(dir1);
 719
 720       AliTPCseed * seed0 = (AliTPCseed*) tpcSeeds.At(i);
 721       AliTPCseed * seed1 = (AliTPCseed*) tpcSeeds.At(j);
 722       if (! seed0) continue;
 723       if (! seed1) continue;
 724       //
 725       Float_t dir = (dir0[0]*dir1[0] + dir0[1]*dir1[1] + dir0[2]*dir1[2]);
 726       Float_t d0  = track0->GetLinearD(0,0);
 727       Float_t d1  = track1->GetLinearD(0,0);
 728       //
 729       // conservative cuts - convergence to be guarantied
 730       // applying before track propagation
 731       if (TMath::Abs(d0+d1)>fCutMaxD) continue;   // distance to the 0,0
 732       if (dir>fCutMinDir) continue;               // direction vector product
 733       Float_t bz = AliTracker::GetBz();
 734       Float_t dvertex0[2];   //distance to 0,0
 735       Float_t dvertex1[2];   //distance to 0,0
 736       track0->GetDZ(0,0,0,bz,dvertex0);
 737       track1->GetDZ(0,0,0,bz,dvertex1);
 738       if (TMath::Abs(dvertex0[1])>250) continue;
 739       if (TMath::Abs(dvertex1[1])>250) continue;
 740       //
 741       //
 742       //
 743       Float_t dmax = TMath::Max(TMath::Abs(d0),TMath::Abs(d1));
 744       AliExternalTrackParam param0(*track0);
 745       AliExternalTrackParam param1(*track1);
 746       //
 747       // Propagate using Magnetic field and correct fo material budget
 748       //
 749       AliTracker::PropagateTrackToBxByBz(&param0,dmax+1,0.0005,3,kTRUE);
 750       AliTracker::PropagateTrackToBxByBz(&param1,dmax+1,0.0005,3,kTRUE);
 751       //
 752       // Propagate rest to the 0,0 DCA - z should be ignored
 753       //
 754       Bool_t b0 = param0.PropagateToDCA(&vtx,bz,1000);
 755       Bool_t b1 = param1.PropagateToDCA(&vtx,bz,1000);
 756       //
 757       param0.GetDZ(0,0,0,bz,dvertex0);
 758       param1.GetDZ(0,0,0,bz,dvertex1);
 759       if (TMath::Abs(param0.GetZ()-param1.GetZ())>fCutMaxDz) continue;
 760       //
 761       Double_t xyz0[3];//,pxyz0[3];
 762       Double_t xyz1[3];//,pxyz1[3];
 763       param0.GetXYZ(xyz0);
 764       param1.GetXYZ(xyz1);
 765       Bool_t isPair = IsPair(&param0,&param1);
 766       //
 767       // HERE WE WILL PUT THE ACCESS TO THE MC TRACKS AND MATCH THESE !!!!
 768       //
 769       Int_t label0 = TMath::Abs(track0->GetLabel());
 770       AliMCParticle *mcParticle0 = fMCinfo->GetTrack(label0);
 771       TParticle *particle0 = mcParticle0->Particle();
 772       AliTrackReference *ref0 = GetFirstTPCTrackRef(mcParticle0); // get the first TPC track reference
 773       if (!ref0) continue;
 774       AliExternalTrackParam *paramMC0 = 0;
 775       paramMC0 = MakeTrack(ref0, particle0);
 776       //
 777       Int_t label1 = TMath::Abs(track1->GetLabel());
 778       AliMCParticle *mcParticle1 = fMCinfo->GetTrack(label1);
 779       TParticle *particle1 = mcParticle1->Particle();
 780       AliTrackReference *ref1 = GetFirstTPCTrackRef(mcParticle1); // get the first TPC track reference
 781       if (!ref1) continue;
 782       AliExternalTrackParam *paramMC1 = 0;
 783       paramMC1 = MakeTrack(ref1, particle1);
 784       //
 785       // ACCESS TOF INFORMATION
 786       Int_t nTrackRefTOF0 = 0;
 787       Int_t nTrackRefITS0 = 0;
 788       AliTrackReference * refLastTOF0 = 0;
 789       AliTrackReference * refFirstTOF0 = GetAllTOFinfo(mcParticle0, nTrackRefTOF0, nTrackRefITS0);
 790       Float_t alphaTOF0 = 0;
 791       if (refFirstTOF0) alphaTOF0 = refFirstTOF0->Alpha();
 792       //
 793       Int_t nTrackRefTOF1 = 0;
 794       Int_t nTrackRefITS1 = 0;
 795       AliTrackReference  * refLastTOF1 = 0;
 796       AliTrackReference  * refFirstTOF1 =GetAllTOFinfo(mcParticle1, nTrackRefTOF1, nTrackRefITS1);
 797       Float_t alphaTOF1 = 0;
 798       if (refFirstTOF1) alphaTOF1 = refFirstTOF1->Alpha();
 799       //cout << " STATUS "<<nTrackRefTOF0<<" "<<refFirstTOF0<<" "<<refLastTOF0<<" " <<refFirstTOF1<<" "<<refLastTOF1<<endl;
 800       //
 801       //
 802       //
 803       if (fStreamLevel>0){
 804         TTreeSRedirector * cstream =  GetDebugStreamer();
 805         AliExternalTrackParam *ip0 = (AliExternalTrackParam *)track0->GetInnerParam();
 806         AliExternalTrackParam *ip1 = (AliExternalTrackParam *)track1->GetInnerParam();
 807         AliExternalTrackParam *op0 = (AliExternalTrackParam *)track0->GetOuterParam();
 808         AliExternalTrackParam *op1 = (AliExternalTrackParam *)track1->GetOuterParam();
 809         //
 810         //
 811         //
 812         if (cstream) {
 813           (*cstream) << "Track0" <<
 814             "dir="<<dir<<               //  direction
 815             "OK="<<isPair<<             //  will be accepted
 816             "b0="<<b0<<                 //  propagate status
 817             "b1="<<b1<<                 //  propagate status
 818             //
 819             "Particle.="<<particle0<<             // TParticle with generated momentum
 820             "NTrackRefTOF0="<<nTrackRefTOF0<<     // Number of TOF track references upper
 821             "NTrackRefTOF1="<<nTrackRefTOF1<<     // Number of TOF track references lower
 822             "NTrackRefITS0="<<nTrackRefITS0<<     // Number of ITS track references upper
 823             "NTrackRefITS1="<<nTrackRefITS1<<     // Number of ITS track references lower
 824             "Alpha0="<<alphaTOF0<<                   // alpha upper
 825             "Alpha1="<<alphaTOF1<<                   // alpha lower
 826             "RefFirstTOF0.="<<refFirstTOF0<<        // first tof reference upper
 827             "RefLastTOF0.="<<refLastTOF0<<          // last tof reference upper
 828             "RefFirstTOF1.="<<refFirstTOF1<<        // first tof reference lower
 829             "RefLastTOF1.="<<refLastTOF1<<          // last tof reference lower
 830             //
 831             "Orig0.=" << track0 <<      //  original track  0
 832             "Orig1.=" << track1 <<      //  original track  1
 833             "Tr0.="<<&param0<<          //  track propagated to the DCA 0,0
 834             "Tr1.="<<&param1<<          //  track propagated to the DCA 0,0
 835             "Ip0.="<<ip0<<              //  inner param - upper
 836             "Ip1.="<<ip1<<              //  inner param - lower
 837             "Op0.="<<op0<<              //  outer param - upper
 838             "Op1.="<<op1<<              //  outer param - lower
 839             //
 840             "paramTrackRef0.="<<paramMC0<< //  "ideal" MC true track parameters from track references - upper
 841             "paramTrackRef1.="<<paramMC1<< //  "ideal" MC true track parameters from track references - lower
 842             //
 843             "v00="<<dvertex0[0]<<       //  distance using kalman
 844             "v01="<<dvertex0[1]<<       //
 845             "v10="<<dvertex1[0]<<       //
 846             "v11="<<dvertex1[1]<<       //
 847             "d0="<<d0<<                 //  linear distance to 0,0
 848             "d1="<<d1<<                 //  linear distance to 0,0
 849             //
 850             //
 851             //
 852             "x00="<<xyz0[0]<<           // global position close to vertex
 853             "x01="<<xyz0[1]<<
 854             "x02="<<xyz0[2]<<
 855             //
 856             "x10="<<xyz1[0]<<           // global position close to vertex
 857             "x11="<<xyz1[1]<<
 858             "x12="<<xyz1[2]<<
 859             //
 860             "dir00="<<dir0[0]<<           // direction upper
 861             "dir01="<<dir0[1]<<
 862             "dir02="<<dir0[2]<<
 863             //
 864             "dir10="<<dir1[0]<<           // direction lower
 865             "dir11="<<dir1[1]<<
 866             "dir12="<<dir1[2]<<
 867             //
 868             //
 869             "Seed0.=" << seed0 <<       //  original seed 0
 870             "Seed1.=" << seed1 <<       //  original seed 1
 871             //
 872             "\n";
 873         }
 874       }
 875       delete paramMC0;
 876       delete paramMC1;
 877     }
 878   }
 879
 880   return;
 881
 882
 883 }
 884
 885 Bool_t  AliMaterialBudget::IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1){
 886   //
 887   //
 888   /*
 889   // 0. Same direction - OPOSITE  - cutDir +cutT
 890   TCut cutDir("cutDir","dir<-0.99")
 891   // 1.
 892   TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03")
 893   //
 894   // 2. The same rphi
 895   TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5")
 896   //
 897   //
 898   //
 899   TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
 900   // 1/Pt diff cut
 901   */
 902   const Double_t *p0 = tr0->GetParameter();
 903   const Double_t *p1 = tr1->GetParameter();
 904   if (TMath::Abs(p0[3]+p1[3])>fCutTheta) return kFALSE;
 905   if (TMath::Abs(p0[1]-p1[1])>fCutMaxDz) return kFALSE;
 906   if (TMath::Abs(p0[0]+p1[0])>fCutMaxD)  return kFALSE;
 907
 908   Double_t d0[3], d1[3];
 909   tr0->GetDirection(d0);
 910   tr1->GetDirection(d1);
 911   if (d0[0]*d1[0] + d0[1]*d1[1] + d0[2]*d1[2] >fCutMinDir) return kFALSE;
 912   //
 913   return kTRUE;
 914 }
 915
 916
 917 AliTrackReference * AliMaterialBudget::GetFirstTPCTrackRef(AliMCParticle *mcParticle)
 918 {
 919   // return first TPC track reference
 920   if(!mcParticle) return 0;
 921
 922   // find first track reference
 923   // check direction to select proper reference point for looping tracks
 924   Int_t nTrackRef = mcParticle->GetNumberOfTrackReferences();
 925   AliTrackReference *ref = 0;
 926   AliTrackReference *refIn = 0;
 927   for (Int_t iref = 0; iref < nTrackRef; iref++) {
 928     ref = mcParticle->GetTrackReference(iref);
 929     if(ref && (ref->DetectorId()==AliTrackReference::kTPC))
 930     {
 931       //Float_t dir = ref->X()*ref->Px()+ref->Y()*ref->Py();
 932       //if(dir < 0.) break;
 933
 934       refIn = ref;
 935       break;
 936     }
 937   }
 938
 939   return refIn;
 940 }
 941
 942
 943 AliTrackReference *  AliMaterialBudget::GetAllTOFinfo(AliMCParticle *mcParticle, Int_t &nTrackRef, Int_t &nTrackRefITS, Int_t retValue) {
 944   //
 945   //
 946   //
 947
 948   if(!mcParticle) return 0;
 949   Int_t counter = 0;
 950   nTrackRef = 0;
 951   nTrackRefITS = 0;
 952   AliTrackReference *ref = 0;
 953   for (Int_t iref = 0; iref < mcParticle->GetNumberOfTrackReferences(); iref++) {
 954     ref = mcParticle->GetTrackReference(iref);
 955     if(ref && (ref->DetectorId()==AliTrackReference::kTOF)) {
 956       counter = iref;
 957       nTrackRef++;
 958     }
 959     if(ref && (ref->DetectorId()==AliTrackReference::kITS)) nTrackRefITS++;
 960   }
 961   if (nTrackRef ==0) return 0;
 962   if (retValue == 0) return mcParticle->GetTrackReference(counter - nTrackRef +1);
 963   return mcParticle->GetTrackReference(counter);
 964
 965 }
 966
 967
 968 void AliMaterialBudget::FinishTaskOutput()
 969 {
 970   //
 971   // According description in AliAnalisysTask this method is call
 972   // on the slaves before sending data
 973   //
 974   Terminate("slave");
 975   gSystem->Exec("pwd");
 976   RegisterDebugOutput();
 977
 978 }
 979
 980
 981 void AliMaterialBudget::RegisterDebugOutput(){
 982   //
 983   //
 984   //
 985   //
 986   // store  - copy debug output to the destination position
 987   // currently ONLY for local copy
 988   TString dsName;
 989   dsName=GetName();
 990   dsName+="Debug.root";
 991   dsName.ReplaceAll(" ","");
 992   TString dsName2=fDebugOutputPath.Data();
 993   gSystem->MakeDirectory(dsName2.Data());
 994   dsName2+=gSystem->HostName();
 995   gSystem->MakeDirectory(dsName2.Data());
 996   dsName2+="/";
 997   dsName2+=gSystem->BaseName(gSystem->pwd());
 998   dsName2+="/";
 999   gSystem->MakeDirectory(dsName2.Data());
1000   dsName2+=dsName;
1001   AliInfo(Form("copy %s\t%s\n",dsName.Data(),dsName2.Data()));
1002   printf("copy %s\t%s\n",dsName.Data(),dsName2.Data());
1003   TFile::Cp(dsName.Data(),dsName2.Data());
1004 }
1005
1006
1007 Bool_t AliMaterialBudget::PropagateCosmicToDCA(AliExternalTrackParam *param0, AliExternalTrackParam *param1, Double_t mass){
1008   //
1009   // param0 - upper part of cosmic track
1010   // param1 - lower part of cosmic track
1011   //
1012   // 0. Propagate both tracks to DCA to (0,0,0)
1013   // 1. After propagation to DCA rotate track param1 to corrdinate system of track1 <-> rotate param0 to coordinate system of param 1 ????
1014   // 2. Propagate track 1 to refernce x from track0
1015   //
1016
1017   // step 0.
1018
1019   Float_t d0  = param0->GetLinearD(0,0);
1020   Float_t d1  = param1->GetLinearD(0,0);
1021   Float_t dmax = TMath::Max(TMath::Abs(d0),TMath::Abs(d1));
1022   //
1023   // propagate in the beginning taking all material into account
1024   //
1025   AliTracker::PropagateTrackToBxByBz(param0,dmax+1.,mass,0.5,kTRUE,0.99,-1.);
1026   AliTracker::PropagateTrackToBxByBz(param1,dmax+1.,mass,0.5,kTRUE,0.99,1.);
1027   //
1028   Double_t vtxx[3]={0,0,0};
1029   Double_t svtxx[3]={0.000001,0.000001,100.};
1030   AliESDVertex vtx(vtxx,svtxx);
1031   //
1032   Bool_t b0 = param0->PropagateToDCA(&vtx,AliTracker::GetBz(),1000);
1033   Bool_t b1 = param1->PropagateToDCA(&vtx,AliTracker::GetBz(),1000);
1034
1035   if (!(b0 && b1)) return 0;
1036
1037   // step 1.
1038
1039   Float_t dAlpha = param0->GetAlpha();
1040   param1->Rotate(dAlpha);
1041
1042   // step 2.
1043
1044   Float_t refX = param0->GetX();
1045   param1->PropagateTo(refX,AliTracker::GetBz());
1046
1047   return kTRUE;
1048
1049
1050 }
1051
1052