PWGPP/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
 303   const Double_t kPcut=0.05;
 304   const Double_t kMinDrITS = 2.;   // minimal distance between references
 305   const Double_t kMinDrTRD = 8.;   // minimal distance between references
 306   const Double_t kMinDrTOF = 10.;   // minimal distance between references
 307   //
 308   //
 309   // Process tracks
 310   //
 311   Int_t npart = fMCinfo->GetNumberOfTracks();
 312   if (npart==0) return;
 313   Double_t vertex[4]={0,0,0,0};
 314
 315   TClonesArray * trefs = new TClonesArray("AliTrackReference");
 316   TParticle * particle= new TParticle;
 317
 318   if(particle && trefs) {
 319     fMCinfo->GetParticleAndTR(0, particle, trefs);
 320   }
 321
 322
 323   if (particle){
 324     vertex[0]=particle->Vx();
 325     vertex[1]=particle->Vy();
 326     vertex[2]=particle->Vz();
 327     vertex[3]=particle->R();
 328   }
 329   //
 330   //
 331   AliTrackReference dummy,*pdummy= &dummy;
 332   AliExternalTrackParam epdummy,*pepdummy= &epdummy;
 333   Int_t nRefITS =0;
 334   Int_t nRefTPC =0;
 335   Int_t nRefTRD =0;
 336   Int_t nRefTOF =0;
 337   AliTrackReference * refITS0, *refITS1;
 338   AliTrackReference * refTPC0, *refTPC1;
 339   AliTrackReference * refTPCIn0, *refTPCIn1;
 340   AliTrackReference * refTRD0, *refTRD1;
 341   AliTrackReference * refTOF0, *refTOF1;
 342   AliTrackReference *refMinR;
 343   //
 344   if(!particle) return;
 345   for (Int_t ipart=0;ipart<npart;ipart++){
 346     //Int_t status = fMCinfo->GetParticleAndTR(ipart, particle, trefs);
 347     AliMCParticle * pp = (AliMCParticle*) fMCinfo->GetTrack(ipart);
 348     if (!pp) continue;
 349     if (particle->P()<kPcut) continue;
 350     Double_t mass = particle->GetMass();
 351
 352     // RESET
 353     nRefITS =0;
 354     nRefTPC =0;
 355     nRefTRD =0;
 356     nRefTOF =0;
 357     refITS0=pdummy;    refITS1=pdummy;
 358     refTPC0=pdummy;    refTPC1=pdummy;
 359     refTPCIn0=pdummy;    refTPCIn1=pdummy;
 360     refTRD0=pdummy;    refTRD1=pdummy;
 361     refTOF0=pdummy;    refTOF1=pdummy;
 362     refMinR = pdummy;
 363     //
 364     Int_t nref = pp->GetNumberOfTrackReferences();
 365     if (nref==0) continue;
 366     for (Int_t iref = 0; iref < nref; iref++) {
 367       AliTrackReference *ref = pp->GetTrackReference(iref);
 368       if (ref->DetectorId()==AliTrackReference::kDisappeared) continue;
 369       //      if (ref.Px()*particle.Px()+ref.Py()*particle.Py()<0) break; // returning track
 370       //
 371       if (ref->DetectorId()==AliTrackReference::kITS){
 372         if (TMath::Abs(ref->R()-refITS1->R())>kMinDrITS) {
 373           refITS1 = ref;
 374           nRefITS++;
 375         }
 376         if (refITS0==pdummy) refITS0=ref;
 377       }
 378       if (ref->DetectorId()==AliTrackReference::kTPC){
 379         nRefTPC++;
 380         refTPC1 = ref;
 381         if (refTPC0==pdummy) refTPC0=ref;
 382       }
 383       if (ref->DetectorId()==AliTrackReference::kTRD){
 384         if (TMath::Abs(ref->R()-refTRD1->R())>kMinDrTRD) {
 385           refTRD1 = ref;
 386           nRefTRD++;
 387         }
 388         if (refTRD0==pdummy) refTRD0=ref;
 389       }
 390       if (ref->DetectorId()==AliTrackReference::kTOF){
 391         if (TMath::Abs(ref->X()-refTOF1->X()) + TMath::Abs(ref->Y()-refTOF1->Y()) + TMath::Abs(ref->Z()-refTOF1->Z())>kMinDrTOF) {
 392           refTOF1 = ref;
 393           nRefTOF++;
 394         }
 395         if (refTOF0==pdummy) refTOF0=ref;
 396       }
 397       //
 398       // "find inner track ref"
 399       if (ref->DetectorId()==AliTrackReference::kTPC){
 400         if (ref->Px()*ref->X()+ref->Py()*ref->Y()<0){
 401           //  track in
 402           if (refTPCIn0 == pdummy) refTPCIn0=ref;
 403           if (refTPCIn0 != pdummy &&  refTPCIn0->R()>ref->R())
 404             refTPCIn0=ref;
 405         }
 406         if (ref->Px()*ref->X()+ref->Py()*ref->Y()>0){
 407           //  track in
 408           if (refTPCIn1 == pdummy) refTPCIn1=ref;
 409           if (refTPCIn1 != pdummy &&  refTPCIn1->R()>ref->R())
 410             refTPCIn1=ref;
 411         }
 412       }
 413
 414
 415       if (refMinR==pdummy && ref->P()>0  ){
 416         refMinR=ref;
 417       }
 418       if (refMinR->R()>ref->R() && ref->P()>0 ) refMinR=ref;
 419
 420     }
 421     //
 422     AliExternalTrackParam * trackMC = pepdummy;
 423     //track0->GetDZ(0,0,0,bz,dvertex0)
 424     Float_t dist[2]={0,0};
 425     AliMagF* field = (AliMagF*) TGeoGlobalMagField::Instance()->GetField();
 426     Double_t esdfield= fESD->GetMagneticField();
 427     Double_t xyz[3]={0,0,0};
 428     Double_t bxyz[3]={0,0,0};
 429     field->Field(xyz,bxyz);
 430     if (refMinR->P()>0) {
 431       trackMC = MakeTrack(refMinR,particle);
 432       trackMC->GetDZ(0,0,0,bxyz[2],dist);
 433     }
 434     Double_t alphaTOF0 = TMath::ATan2(refTOF0->Y(),refTOF0->X());
 435     Double_t alphaTOF1 = TMath::ATan2(refTOF1->Y(),refTOF1->X());
 436     Int_t dsecTOF   = TMath::Nint(180*(alphaTOF0-alphaTOF1)/(TMath::Pi()*20.)-0.5);
 437     //
 438     // make the two different TPC tracks and propagate them to their DCA
 439     //
 440     Double_t dP = 0;
 441     Bool_t statusProp = kFALSE;
 442     Double_t dY = 0;
 443     Double_t dZ = 0;
 444     AliExternalTrackParam * track0  = pepdummy;
 445     AliExternalTrackParam * track1  = pepdummy;
 446     AliExternalTrackParam * otrack0 = pepdummy;
 447     AliExternalTrackParam * otrack1 = pepdummy;
 448     if (refTPCIn0!=pdummy && refTPCIn1!=pdummy) {
 449       track0 = MakeTrack(refTPCIn0,particle);
 450       track1 = MakeTrack(refTPCIn1,particle);
 451       otrack0 = MakeTrack(refTPCIn0,particle);
 452       otrack1 = MakeTrack(refTPCIn1,particle);
 453       dP = track0->P() - track1->P(); // momentum loss
 454       statusProp = AliMaterialBudget::PropagateCosmicToDCA(track0,track1,mass);
 455       if (statusProp) {
 456         dY = track0->GetY() - track1->GetY();
 457         dZ = track0->GetZ() - track1->GetZ();
 458       }
 459     }
 460     //
 461     TTreeSRedirector *pcstream = GetDebugStreamer();
 462     if (pcstream){
 463       char name[100];
 464       for (Int_t id=0; id<3;id++){
 465
 466         //      if (id==0) sprintf(name,"mcAll"); // all tracks: inconvenient to cut if on is only interest in tracks which reach the TPC
 467         if (id==0) continue; // require TPC
 468         if (id==1) snprintf(name,100,"mcITS");
 469         if (id==2) snprintf(name,100,"mcTPC");
 470         if (id==1&& nRefITS==0) continue;
 471         if (id==2&& nRefTPC==0) continue;
 472
 473         (*pcstream)<<name<<
 474           "ipart="<<ipart<<
 475           "p.="<<particle<<
 476           "mass="<<mass<<
 477           "tbfield="<<bxyz[2]<<
 478           "esdbfield="<<esdfield<<
 479           // counter
 480           "nref="<<nref<<
 481           "nRefITS="<<nRefITS<<
 482           "nRefTPC="<<nRefTPC<<
 483           "nRefTRD="<<nRefTRD<<
 484           "nRefTOF="<<nRefTOF<<
 485           //references
 486           "refMinR.="<<refMinR<<
 487           "refITS0.="<<refITS0<<
 488           "refITS1.="<<refITS1<<
 489           "refTPC0.="<<refTPC0<<
 490           "refTPC1.="<<refTPC1<<
 491           "refTPCIn0.="<<refTPCIn0<<
 492           "refTPCIn1.="<<refTPCIn1<<
 493           "refTRD0.="<<refTRD0<<
 494           "refTRD1.="<<refTRD1<<
 495           "refTOF0.="<<refTOF0<<
 496           "refTOF1.="<<refTOF1<<
 497           //trigger variables
 498           "dsecTOF="<<dsecTOF<<    // delta TOF sectors
 499           "aTOF0="<<alphaTOF0<<
 500           "aTOF1="<<alphaTOF1<<
 501           //
 502           // track
 503           "dr="<<dist[0]<<
 504           "dz="<<dist[1]<<
 505           //
 506           // "two" TPC tracks
 507           "status="<<statusProp<<
 508           "dP="<<dP<<
 509           "otrack0.="<<otrack0<<
 510           "otrack1.="<<otrack1<<
 511           "track0.="<<track0<<
 512           "track1.="<<track1<<
 513           "dY="<<dY<<
 514           "dZ="<<dZ<<
 515           "\n";
 516       }
 517     }
 518   }
 519 }
 520
 521
 522
 523 void AliMaterialBudget::ProcessRefTracker(AliTrackReference* refIn,  AliTrackReference* refOut, TParticle*part,Int_t type){
 524   //
 525   // Test propagation from In to out
 526   //
 527   AliExternalTrackParam *param = 0;
 528   AliExternalTrackParam *paramMC = 0;
 529   Double_t xyzIn[3]={refIn->X(),refIn->Y(), refIn->Z()};
 530   Double_t mass = part->GetMass();
 531   Double_t step=1;
 532   //
 533   param=MakeTrack(refOut,part);
 534   paramMC=MakeTrack(refOut,part);
 535   if (!param) return;
 536   if (type<3) PropagateToPoint(param,xyzIn, mass, step);
 537   if (type==3) {
 538     AliTPCseed seed;
 539     seed.Set(param->GetX(),param->GetAlpha(),param->GetParameter(),param->GetCovariance());
 540     Float_t alpha= TMath::ATan2(refIn->Y(),refIn->X());
 541     if(seed.Rotate(alpha-seed.GetAlpha()) == kFALSE) return;
 542     seed.SetMass(mass);
 543     for (Float_t xlayer= seed.GetX(); xlayer>refIn->R(); xlayer-=step){
 544       seed.PropagateTo(xlayer);
 545     }
 546     seed.PropagateTo(refIn->R());
 547     param->Set(seed.GetX(),seed.GetAlpha(),seed.GetParameter(),seed.GetCovariance());
 548   }
 549   TTreeSRedirector *pcstream = GetDebugStreamer();
 550   TVectorD gpos(3);
 551   TVectorD gmom(3);
 552   param->GetXYZ(gpos.GetMatrixArray());
 553   param->GetPxPyPz(gmom.GetMatrixArray());
 554   if (pcstream){
 555     (*pcstream)<<"MC"<<
 556       "type="<<type<<
 557       "step="<<step<<
 558       "refIn.="<<refIn<<
 559       "refOut.="<<refOut<<
 560       "p.="<<part<<
 561       "par.="<<param<<
 562       "parMC.="<<paramMC<<
 563       "gpos.="<<&gpos<<
 564       "gmom.="<<&gmom<<
 565       "\n";
 566   }
 567 }
 568
 569
 570 void  AliMaterialBudget::FitTrackRefs(TParticle * part, TClonesArray * trefs){
 571   //
 572   //
 573   //
 574   //
 575   const Int_t kMinRefs=6;
 576   Int_t nrefs = trefs->GetEntries();
 577   if (nrefs<kMinRefs) return; // we should have enough references
 578   Int_t iref0 =-1;
 579   Int_t iref1 =-1;
 580
 581   for (Int_t iref=0; iref<nrefs; iref++){
 582     AliTrackReference * ref = (AliTrackReference*)trefs->At(iref);
 583     if (!ref) continue;
 584     Float_t dir = ref->X()*ref->Px()+ref->Y()*ref->Py();
 585     if (dir<0) break;
 586     if (ref->DetectorId()!=AliTrackReference::kTPC) continue;
 587     if (iref0<0) iref0 = iref;
 588     iref1 = iref;
 589   }
 590   if (iref1-iref0<kMinRefs) return;
 591   Double_t covar[15];
 592   for (Int_t icov=0; icov<15; icov++) covar[icov]=0;
 593   covar[0]=1;
 594   covar[2]=1;
 595   covar[5]=1;
 596   covar[9]=1;
 597   covar[14]=1;
 598
 599   AliTrackReference * refIn = (AliTrackReference*)trefs->At(iref0);
 600   AliTrackReference * refOut = (AliTrackReference*)trefs->At(iref1);
 601   AliExternalTrackParam *paramPropagate= MakeTrack(refIn,part);
 602   AliExternalTrackParam *paramUpdate   = MakeTrack(refIn,part);
 603   paramUpdate->AddCovariance(covar);
 604   Double_t mass = part->GetMass();
 605   Double_t charge = part->GetPDG()->Charge()/3.;
 606 /*
 607   Float_t alphaIn= TMath::ATan2(refIn->Y(),refIn->X());
 608   Float_t radiusIn= refIn->R();
 609   Float_t alphaOut= TMath::ATan2(refOut->Y(),refOut->X());
 610   Float_t radiusOut= refOut->R();
 611 */
 612   Bool_t isOKP=kTRUE;
 613   Bool_t isOKU=kTRUE;
 614   AliMagF * field = (AliMagF*) TGeoGlobalMagField::Instance()->GetField();
 615   for (Int_t iref = iref0; iref<=iref1; iref++){
 616     AliTrackReference * ref = (AliTrackReference*)trefs->At(iref);
 617     Float_t alphaC= TMath::ATan2(ref->Y(),ref->X());
 618     Double_t pos[3] = {ref->X(), ref->Y(), ref->Z()};
 619     Double_t mag[3];
 620     field->Field(pos,mag);
 621     isOKP&=paramPropagate->Rotate(alphaC);
 622     isOKU&=paramUpdate->Rotate(alphaC);
 623     for (Float_t xref= paramPropagate->GetX(); xref<ref->R(); xref++){
 624       isOKP&=paramPropagate->PropagateTo(xref, mag[2]);
 625       isOKU&=paramUpdate->PropagateTo(xref, mag[2]);
 626     }
 627     isOKP&=paramPropagate->PropagateTo(ref->R(), mag[2]);
 628     isOKU&=paramUpdate->PropagateTo(ref->R(), mag[2]);
 629     Double_t clpos[2] = {0, ref->Z()};
 630     Double_t clcov[3] = { 0.005,0,0.005};
 631     isOKU&= paramUpdate->Update(clpos, clcov);
 632   }
 633   TTreeSRedirector *pcstream = GetDebugStreamer();
 634   if (pcstream){
 635     TVectorD gposU(3);
 636     TVectorD gmomU(3);
 637     TVectorD gposP(3);
 638     TVectorD gmomP(3);
 639     paramUpdate->GetXYZ(gposU.GetMatrixArray());
 640     paramUpdate->GetPxPyPz(gmomU.GetMatrixArray());
 641     paramPropagate->GetXYZ(gposP.GetMatrixArray());
 642     paramPropagate->GetPxPyPz(gmomP.GetMatrixArray());
 643
 644      (*pcstream)<<"MCupdate"<<
 645        "isOKU="<<isOKU<<
 646        "isOKP="<<isOKP<<
 647        "m="<<mass<<
 648        "q="<<charge<<
 649        "part.="<<part<<
 650        "refIn.="<<refIn<<
 651        "refOut.="<<refOut<<
 652        "pP.="<<paramPropagate<<
 653        "pU.="<<paramUpdate<<
 654        "gposU.="<<&gposU<<
 655        "gmomU.="<<&gmomU<<
 656        "gposP.="<<&gposP<<
 657        "gmomP.="<<&gmomP<<
 658        "\n";
 659    }
 660 }
 661
 662
 663
 664 void AliMaterialBudget::FindPairs(AliESDEvent * event) {
 665   //
 666   // This function matches the cosmic tracks and calculates the energy loss in the material.
 667   // If accessible the "true" energy loss is determined with the MC track references.
 668   //
 669   //
 670   // Find cosmic pairs
 671   //
 672   // Track0 is choosen in upper TPC part
 673   // Track1 is choosen in lower TPC part
 674   //
 675   if(fDebugLevel>3)
 676     cout << "AliMaterialBudget::FindPairs()" << endl;
 677
 678
 679   AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
 680   Int_t ntracks=event->GetNumberOfTracks();
 681   TObjArray  tpcSeeds(ntracks);
 682   if (ntracks==0) return;
 683   Double_t vtxx[3]={0,0,0};
 684   Double_t svtxx[3]={0.000001,0.000001,100.};
 685   AliESDVertex vtx(vtxx,svtxx);
 686   //
 687   //track loop
 688   //
 689   for (Int_t i=0;i<ntracks;++i) {
 690    AliESDtrack *track = event->GetTrack(i);
 691      const AliExternalTrackParam * trackIn = track->GetInnerParam();
 692    const AliExternalTrackParam * trackOut = track->GetOuterParam();
 693    if (!trackIn) continue;
 694    if (!trackOut) continue;
 695    AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
 696    if (!friendTrack) continue;
 697    TObject *calibObject;
 698    AliTPCseed *seed = 0;
 699    for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
 700      if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
 701    }
 702    if (seed) tpcSeeds.AddAt(seed,i);
 703   }
 704
 705   if (ntracks<2) return;
 706   //
 707   // Find pairs
 708   //
 709   for (Int_t i=0;i<ntracks;++i) {
 710     AliESDtrack *track0 = event->GetTrack(i);
 711     // track0 - choosen upper part
 712     if (!track0) continue;
 713     if (!track0->GetOuterParam()) continue;
 714     if (track0->GetOuterParam()->GetAlpha()<0) continue;
 715     Double_t dir0[3];
 716     track0->GetDirection(dir0);
 717     for (Int_t j=0;j<ntracks;++j) {
 718       if (i==j) continue;
 719       AliESDtrack *track1 = event->GetTrack(j);
 720       //track 1 lower part
 721       if (!track1) continue;
 722       if (!track1->GetOuterParam()) continue;
 723       if (track1->GetOuterParam()->GetAlpha()>0) continue;
 724       //
 725       Double_t dir1[3];
 726       track1->GetDirection(dir1);
 727
 728       AliTPCseed * seed0 = (AliTPCseed*) tpcSeeds.At(i);
 729       AliTPCseed * seed1 = (AliTPCseed*) tpcSeeds.At(j);
 730       if (! seed0) continue;
 731       if (! seed1) continue;
 732       //
 733       Float_t dir = (dir0[0]*dir1[0] + dir0[1]*dir1[1] + dir0[2]*dir1[2]);
 734       Float_t d0  = track0->GetLinearD(0,0);
 735       Float_t d1  = track1->GetLinearD(0,0);
 736       //
 737       // conservative cuts - convergence to be guarantied
 738       // applying before track propagation
 739       if (TMath::Abs(d0+d1)>fCutMaxD) continue;   // distance to the 0,0
 740       if (dir>fCutMinDir) continue;               // direction vector product
 741       Float_t bz = AliTracker::GetBz();
 742       Float_t dvertex0[2];   //distance to 0,0
 743       Float_t dvertex1[2];   //distance to 0,0
 744       track0->GetDZ(0,0,0,bz,dvertex0);
 745       track1->GetDZ(0,0,0,bz,dvertex1);
 746       if (TMath::Abs(dvertex0[1])>250) continue;
 747       if (TMath::Abs(dvertex1[1])>250) continue;
 748       //
 749       //
 750       //
 751       Float_t dmax = TMath::Max(TMath::Abs(d0),TMath::Abs(d1));
 752       AliExternalTrackParam param0(*track0);
 753       AliExternalTrackParam param1(*track1);
 754       //
 755       // Propagate using Magnetic field and correct fo material budget
 756       //
 757       AliTracker::PropagateTrackToBxByBz(&param0,dmax+1,0.0005,3,kTRUE);
 758       AliTracker::PropagateTrackToBxByBz(&param1,dmax+1,0.0005,3,kTRUE);
 759       //
 760       // Propagate rest to the 0,0 DCA - z should be ignored
 761       //
 762       Bool_t b0 = param0.PropagateToDCA(&vtx,bz,1000);
 763       Bool_t b1 = param1.PropagateToDCA(&vtx,bz,1000);
 764       //
 765       param0.GetDZ(0,0,0,bz,dvertex0);
 766       param1.GetDZ(0,0,0,bz,dvertex1);
 767       if (TMath::Abs(param0.GetZ()-param1.GetZ())>fCutMaxDz) continue;
 768       //
 769       Double_t xyz0[3];//,pxyz0[3];
 770       Double_t xyz1[3];//,pxyz1[3];
 771       param0.GetXYZ(xyz0);
 772       param1.GetXYZ(xyz1);
 773       Bool_t isPair = IsPair(&param0,&param1);
 774       //
 775       // HERE WE WILL PUT THE ACCESS TO THE MC TRACKS AND MATCH THESE !!!!
 776       //
 777       Int_t label0 = TMath::Abs(track0->GetLabel());
 778       AliMCParticle *mcParticle0 = (AliMCParticle*) fMCinfo->GetTrack(label0);
 779       TParticle *particle0 = mcParticle0->Particle();
 780       AliTrackReference *ref0 = GetFirstTPCTrackRef(mcParticle0); // get the first TPC track reference
 781       if (!ref0) continue;
 782       AliExternalTrackParam *paramMC0 = 0;
 783       paramMC0 = MakeTrack(ref0, particle0);
 784       //
 785       Int_t label1 = TMath::Abs(track1->GetLabel());
 786       AliMCParticle *mcParticle1 = (AliMCParticle*) fMCinfo->GetTrack(label1);
 787       TParticle *particle1 = mcParticle1->Particle();
 788       AliTrackReference *ref1 = GetFirstTPCTrackRef(mcParticle1); // get the first TPC track reference
 789       if (!ref1) continue;
 790       AliExternalTrackParam *paramMC1 = 0;
 791       paramMC1 = MakeTrack(ref1, particle1);
 792       //
 793       // ACCESS TOF INFORMATION
 794       Int_t nTrackRefTOF0 = 0;
 795       Int_t nTrackRefITS0 = 0;
 796       AliTrackReference * refLastTOF0 = 0;
 797       AliTrackReference * refFirstTOF0 = GetAllTOFinfo(mcParticle0, nTrackRefTOF0, nTrackRefITS0);
 798       Float_t alphaTOF0 = 0;
 799       if (refFirstTOF0) alphaTOF0 = refFirstTOF0->Alpha();
 800       //
 801       Int_t nTrackRefTOF1 = 0;
 802       Int_t nTrackRefITS1 = 0;
 803       AliTrackReference  * refLastTOF1 = 0;
 804       AliTrackReference  * refFirstTOF1 =GetAllTOFinfo(mcParticle1, nTrackRefTOF1, nTrackRefITS1);
 805       Float_t alphaTOF1 = 0;
 806       if (refFirstTOF1) alphaTOF1 = refFirstTOF1->Alpha();
 807       //cout << " STATUS "<<nTrackRefTOF0<<" "<<refFirstTOF0<<" "<<refLastTOF0<<" " <<refFirstTOF1<<" "<<refLastTOF1<<endl;
 808       //
 809       //
 810       //
 811       if (fStreamLevel>0){
 812         TTreeSRedirector * cstream =  GetDebugStreamer();
 813         AliExternalTrackParam *ip0 = (AliExternalTrackParam *)track0->GetInnerParam();
 814         AliExternalTrackParam *ip1 = (AliExternalTrackParam *)track1->GetInnerParam();
 815         AliExternalTrackParam *op0 = (AliExternalTrackParam *)track0->GetOuterParam();
 816         AliExternalTrackParam *op1 = (AliExternalTrackParam *)track1->GetOuterParam();
 817         //
 818         //
 819         //
 820         if (cstream) {
 821           (*cstream) << "Track0" <<
 822             "dir="<<dir<<               //  direction
 823             "OK="<<isPair<<             //  will be accepted
 824             "b0="<<b0<<                 //  propagate status
 825             "b1="<<b1<<                 //  propagate status
 826             //
 827             "Particle.="<<particle0<<             // TParticle with generated momentum
 828             "NTrackRefTOF0="<<nTrackRefTOF0<<     // Number of TOF track references upper
 829             "NTrackRefTOF1="<<nTrackRefTOF1<<     // Number of TOF track references lower
 830             "NTrackRefITS0="<<nTrackRefITS0<<     // Number of ITS track references upper
 831             "NTrackRefITS1="<<nTrackRefITS1<<     // Number of ITS track references lower
 832             "Alpha0="<<alphaTOF0<<                   // alpha upper
 833             "Alpha1="<<alphaTOF1<<                   // alpha lower
 834             "RefFirstTOF0.="<<refFirstTOF0<<        // first tof reference upper
 835             "RefLastTOF0.="<<refLastTOF0<<          // last tof reference upper
 836             "RefFirstTOF1.="<<refFirstTOF1<<        // first tof reference lower
 837             "RefLastTOF1.="<<refLastTOF1<<          // last tof reference lower
 838             //
 839             "Orig0.=" << track0 <<      //  original track  0
 840             "Orig1.=" << track1 <<      //  original track  1
 841             "Tr0.="<<&param0<<          //  track propagated to the DCA 0,0
 842             "Tr1.="<<&param1<<          //  track propagated to the DCA 0,0
 843             "Ip0.="<<ip0<<              //  inner param - upper
 844             "Ip1.="<<ip1<<              //  inner param - lower
 845             "Op0.="<<op0<<              //  outer param - upper
 846             "Op1.="<<op1<<              //  outer param - lower
 847             //
 848             "paramTrackRef0.="<<paramMC0<< //  "ideal" MC true track parameters from track references - upper
 849             "paramTrackRef1.="<<paramMC1<< //  "ideal" MC true track parameters from track references - lower
 850             //
 851             "v00="<<dvertex0[0]<<       //  distance using kalman
 852             "v01="<<dvertex0[1]<<       //
 853             "v10="<<dvertex1[0]<<       //
 854             "v11="<<dvertex1[1]<<       //
 855             "d0="<<d0<<                 //  linear distance to 0,0
 856             "d1="<<d1<<                 //  linear distance to 0,0
 857             //
 858             //
 859             //
 860             "x00="<<xyz0[0]<<           // global position close to vertex
 861             "x01="<<xyz0[1]<<
 862             "x02="<<xyz0[2]<<
 863             //
 864             "x10="<<xyz1[0]<<           // global position close to vertex
 865             "x11="<<xyz1[1]<<
 866             "x12="<<xyz1[2]<<
 867             //
 868             "dir00="<<dir0[0]<<           // direction upper
 869             "dir01="<<dir0[1]<<
 870             "dir02="<<dir0[2]<<
 871             //
 872             "dir10="<<dir1[0]<<           // direction lower
 873             "dir11="<<dir1[1]<<
 874             "dir12="<<dir1[2]<<
 875             //
 876             //
 877             "Seed0.=" << seed0 <<       //  original seed 0
 878             "Seed1.=" << seed1 <<       //  original seed 1
 879             //
 880             "\n";
 881         }
 882       }
 883       delete paramMC0;
 884       delete paramMC1;
 885     }
 886   }
 887
 888   return;
 889
 890
 891 }
 892
 893 Bool_t  AliMaterialBudget::IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1){
 894   //
 895   //
 896   /*
 897   // 0. Same direction - OPOSITE  - cutDir +cutT
 898   TCut cutDir("cutDir","dir<-0.99")
 899   // 1.
 900   TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03")
 901   //
 902   // 2. The same rphi
 903   TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5")
 904   //
 905   //
 906   //
 907   TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
 908   // 1/Pt diff cut
 909   */
 910   const Double_t *p0 = tr0->GetParameter();
 911   const Double_t *p1 = tr1->GetParameter();
 912   if (TMath::Abs(p0[3]+p1[3])>fCutTheta) return kFALSE;
 913   if (TMath::Abs(p0[1]-p1[1])>fCutMaxDz) return kFALSE;
 914   if (TMath::Abs(p0[0]+p1[0])>fCutMaxD)  return kFALSE;
 915
 916   Double_t d0[3], d1[3];
 917   tr0->GetDirection(d0);
 918   tr1->GetDirection(d1);
 919   if (d0[0]*d1[0] + d0[1]*d1[1] + d0[2]*d1[2] >fCutMinDir) return kFALSE;
 920   //
 921   return kTRUE;
 922 }
 923
 924
 925 AliTrackReference * AliMaterialBudget::GetFirstTPCTrackRef(AliMCParticle *mcParticle)
 926 {
 927   // return first TPC track reference
 928   if(!mcParticle) return 0;
 929
 930   // find first track reference
 931   // check direction to select proper reference point for looping tracks
 932   Int_t nTrackRef = mcParticle->GetNumberOfTrackReferences();
 933   AliTrackReference *ref = 0;
 934   AliTrackReference *refIn = 0;
 935   for (Int_t iref = 0; iref < nTrackRef; iref++) {
 936     ref = mcParticle->GetTrackReference(iref);
 937     if(ref && (ref->DetectorId()==AliTrackReference::kTPC))
 938     {
 939       //Float_t dir = ref->X()*ref->Px()+ref->Y()*ref->Py();
 940       //if(dir < 0.) break;
 941
 942       refIn = ref;
 943       break;
 944     }
 945   }
 946
 947   return refIn;
 948 }
 949
 950
 951 AliTrackReference *  AliMaterialBudget::GetAllTOFinfo(AliMCParticle *mcParticle, Int_t &nTrackRef, Int_t &nTrackRefITS, Int_t retValue) {
 952   //
 953   //
 954   //
 955
 956   if(!mcParticle) return 0;
 957   Int_t counter = 0;
 958   nTrackRef = 0;
 959   nTrackRefITS = 0;
 960   AliTrackReference *ref = 0;
 961   for (Int_t iref = 0; iref < mcParticle->GetNumberOfTrackReferences(); iref++) {
 962     ref = mcParticle->GetTrackReference(iref);
 963     if(ref && (ref->DetectorId()==AliTrackReference::kTOF)) {
 964       counter = iref;
 965       nTrackRef++;
 966     }
 967     if(ref && (ref->DetectorId()==AliTrackReference::kITS)) nTrackRefITS++;
 968   }
 969   if (nTrackRef ==0) return 0;
 970   if (retValue == 0) return mcParticle->GetTrackReference(counter - nTrackRef +1);
 971   return mcParticle->GetTrackReference(counter);
 972
 973 }
 974
 975
 976 void AliMaterialBudget::FinishTaskOutput()
 977 {
 978   //
 979   // According description in AliAnalisysTask this method is call
 980   // on the slaves before sending data
 981   //
 982   Terminate("slave");
 983   gSystem->Exec("pwd");
 984   RegisterDebugOutput();
 985
 986 }
 987
 988
 989 void AliMaterialBudget::RegisterDebugOutput(){
 990   //
 991   //
 992   //
 993   //
 994   // store  - copy debug output to the destination position
 995   // currently ONLY for local copy
 996   TString dsName;
 997   dsName=GetName();
 998   dsName+="Debug.root";
 999   dsName.ReplaceAll(" ","");
1000   TString dsName2=fDebugOutputPath.Data();
1001   gSystem->MakeDirectory(dsName2.Data());
1002   dsName2+=gSystem->HostName();
1003   gSystem->MakeDirectory(dsName2.Data());
1004   dsName2+="/";
1005   dsName2+=gSystem->BaseName(gSystem->pwd());
1006   dsName2+="/";
1007   gSystem->MakeDirectory(dsName2.Data());
1008   dsName2+=dsName;
1009   AliInfo(Form("copy %s\t%s\n",dsName.Data(),dsName2.Data()));
1010   printf("copy %s\t%s\n",dsName.Data(),dsName2.Data());
1011   TFile::Cp(dsName.Data(),dsName2.Data());
1012 }
1013
1014
1015 Bool_t AliMaterialBudget::PropagateCosmicToDCA(AliExternalTrackParam *param0, AliExternalTrackParam *param1, Double_t mass){
1016   //
1017   // param0 - upper part of cosmic track
1018   // param1 - lower part of cosmic track
1019   //
1020   // 0. Propagate both tracks to DCA to (0,0,0)
1021   // 1. After propagation to DCA rotate track param1 to corrdinate system of track1 <-> rotate param0 to coordinate system of param 1 ????
1022   // 2. Propagate track 1 to refernce x from track0
1023   //
1024
1025   // step 0.
1026
1027   Float_t d0  = param0->GetLinearD(0,0);
1028   Float_t d1  = param1->GetLinearD(0,0);
1029   Float_t dmax = TMath::Max(TMath::Abs(d0),TMath::Abs(d1));
1030   //
1031   // propagate in the beginning taking all material into account
1032   //
1033   AliTracker::PropagateTrackToBxByBz(param0,dmax+1.,mass,0.5,kTRUE,0.99,-1.);
1034   AliTracker::PropagateTrackToBxByBz(param1,dmax+1.,mass,0.5,kTRUE,0.99,1.);
1035   //
1036   Double_t vtxx[3]={0,0,0};
1037   Double_t svtxx[3]={0.000001,0.000001,100.};
1038   AliESDVertex vtx(vtxx,svtxx);
1039   //
1040   Bool_t b0 = param0->PropagateToDCA(&vtx,AliTracker::GetBz(),1000);
1041   Bool_t b1 = param1->PropagateToDCA(&vtx,AliTracker::GetBz(),1000);
1042
1043   if (!(b0 && b1)) return 0;
1044
1045   // step 1.
1046
1047   Float_t dAlpha = param0->GetAlpha();
1048   param1->Rotate(dAlpha);
1049
1050   // step 2.
1051
1052   Float_t refX = param0->GetX();
1053   param1->PropagateTo(refX,AliTracker::GetBz());
1054
1055   return kTRUE;
1056
1057
1058 }
1059
1060