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