PWG1/AliMCTrackingTestTask.cxx

   1 //
   2 // This class is the task to check the
   3 // Propagation method used in the
   4 //               1. AliExternalTrackParam
   5 //               2. AliTracker
   6 //
   7 // Principle - Creates AliExternalTrackParam form 1 Track Refernece -
   8 //             Propagate it to other
   9 // Magnetic filed and the geomtry has to bec
  10
  11 //
  12 // ROOT includes
  13 #include <TChain.h>
  14 #include <TMath.h>
  15 #include <TVectorD.h>
  16 #include <TSystem.h>
  17 #include <TFile.h>
  18
  19 // ALIROOT includes
  20 #include <TTreeStream.h>
  21 #include <AliAnalysisManager.h>
  22 #include <AliESDInputHandler.h>
  23 #include "AliStack.h"
  24 #include "AliMCEvent.h"
  25 #include "AliMCEventHandler.h"
  26
  27 #include <AliESD.h>
  28 #include "AliMCTrackingTestTask.h"
  29 #include "AliGenInfoMaker.h"
  30 #include "AliHelix.h"
  31
  32 //
  33 #include "AliMCInfo.h"
  34 #include "AliComparisonObject.h"
  35 #include "AliESDRecInfo.h"
  36 #include "AliTPCParamSR.h"
  37 #include "AliTracker.h"
  38 #include "AliTPCseed.h"
  39
  40 // STL includes
  41 #include <iostream>
  42
  43 using namespace std;
  44
  45 ClassImp(AliMCTrackingTestTask)
  46
  47 //________________________________________________________________________
  48 AliMCTrackingTestTask::AliMCTrackingTestTask() :
  49   AliAnalysisTask(),
  50   fMCinfo(0),     //! MC event handler
  51   fESD(0),
  52   fDebugStreamer(0),
  53   fStreamLevel(0),
  54   fDebugLevel(0),
  55   fDebugOutputPath()
  56 {
  57   //
  58   // Default constructor (should not be used)
  59   //
  60 }
  61
  62 AliMCTrackingTestTask::AliMCTrackingTestTask(const AliMCTrackingTestTask& /*info*/) :
  63   AliAnalysisTask(),
  64   fMCinfo(0),     //! MC event handler
  65   fESD(0),
  66   //
  67   fDebugStreamer(0),
  68   fStreamLevel(0),
  69   fDebugLevel(),
  70   fDebugOutputPath()
  71 {
  72   //
  73   // Default constructor
  74   //
  75 }
  76
  77
  78
  79 //________________________________________________________________________
  80 AliMCTrackingTestTask::AliMCTrackingTestTask(const char *name) :
  81   AliAnalysisTask(name, "AliMCTrackingTestTask"),
  82   fMCinfo(0),     //! MC event handler
  83   fESD(0),
  84   fDebugStreamer(0),
  85   fStreamLevel(0),
  86   fDebugLevel(0),
  87   fDebugOutputPath()
  88 {
  89   //
  90   // Normal constructor
  91   //
  92   // Input slot #0 works with a TChain
  93   DefineInput(0, TChain::Class());
  94   // Output slot #0 writes into a TList
  95   DefineOutput(0, TObjArray::Class());
  96   //
  97   //
  98 }
  99
 100 AliMCTrackingTestTask::~AliMCTrackingTestTask(){
 101   //
 102   //
 103   //
 104   if (fDebugLevel>0)  printf("AliMCTrackingTestTask::~AliMCTrackingTestTask\n");
 105   if (fDebugStreamer) delete fDebugStreamer;
 106   fDebugStreamer=0;
 107 }
 108
 109
 110 //________________________________________________________________________
 111 void AliMCTrackingTestTask::ConnectInputData(Option_t *)
 112 {
 113   //
 114   // Connect the input data
 115   //
 116   if(fDebugLevel>3)
 117     cout << "AnalysisTaskTPCCluster::ConnectInputData()" << endl;
 118
 119   TTree* tree=dynamic_cast<TTree*>(GetInputData(0));
 120   if (!tree) {
 121     //Printf("ERROR: Could not read chain from input slot 0");
 122   }
 123   else {
 124     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 125     if (!esdH) {
 126       //Printf("ERROR: Could not get ESDInputHandler");
 127     }
 128     else {
 129       fESD = esdH->GetEvent();
 130       //Printf("*** CONNECTED NEW EVENT ****");
 131     }
 132   }
 133   AliMCEventHandler* mcinfo = (AliMCEventHandler*) (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 134   mcinfo->SetReadTR(kTRUE);
 135
 136   fMCinfo = mcinfo->MCEvent();
 137
 138
 139 }
 140
 141
 142
 143
 144
 145
 146 //________________________________________________________________________
 147 void AliMCTrackingTestTask::CreateOutputObjects()
 148 {
 149   //
 150   // Connect the output objects
 151   //
 152   if(fDebugLevel>3)
 153     cout << "AnalysisTaskTPCCluster::CreateOutputObjects()" << endl;
 154
 155 }
 156
 157
 158 //________________________________________________________________________
 159 void AliMCTrackingTestTask::Exec(Option_t *) {
 160   //
 161   // Execute analysis for current event
 162   //
 163
 164   if(fDebugLevel>3)
 165     cout << "AliMCTrackingTestTask::Exec()" << endl;
 166
 167
 168   // If MC has been connected
 169
 170   if (!fMCinfo){
 171     cout << "Not MC info\n" << endl;
 172   }else{
 173     ProcessMCInfo();
 174     //mcinfo->Print();
 175     //DumpInfo();
 176   }
 177   //
 178 }
 179
 180
 181
 182
 183 //________________________________________________________________________
 184 void AliMCTrackingTestTask::Terminate(Option_t *) {
 185     //
 186     // Terminate loop
 187     //
 188   if(fDebugLevel>3)
 189     printf("AliMCTrackingTestTask: Terminate() \n");
 190   //
 191   if (fDebugLevel>0) printf("AliMCtrackingTestTask::Terminate\n");
 192   if (fDebugStreamer) delete fDebugStreamer;
 193   fDebugStreamer = 0;
 194   return;
 195 }
 196
 197
 198
 199 TTreeSRedirector *AliMCTrackingTestTask::GetDebugStreamer(){
 200   //
 201   // Get Debug streamer
 202   // In case debug streamer not yet initialized and StreamLevel>0 create new one
 203   //
 204   if (fStreamLevel==0) return 0;
 205   if (fDebugStreamer) return fDebugStreamer;
 206   TString dsName;
 207   dsName=GetName();
 208   dsName+="Debug.root";
 209   dsName.ReplaceAll(" ","");
 210   fDebugStreamer = new TTreeSRedirector(dsName.Data());
 211   return fDebugStreamer;
 212 }
 213
 214
 215
 216
 217
 218
 219 AliExternalTrackParam * AliMCTrackingTestTask::MakeTrack(const AliTrackReference* ref, TParticle*part)
 220 {
 221   //
 222   // Make track out of the track ref
 223   // part - TParticle used to determine chargr
 224   // the covariance matrix - equal 0 - starting from ideal MC position
 225   Double_t xyz[3]={ref->X(),ref->Y(),ref->Z()};
 226   Double_t pxyz[3]={ref->Px(),ref->Py(),ref->Pz()};
 227   Double_t cv[21];
 228   for (Int_t i=0; i<21;i++) cv[i]=0;
 229   if (!part->GetPDG()) return 0;
 230   AliExternalTrackParam * param = new AliExternalTrackParam(xyz,pxyz,cv,TMath::Nint(part->GetPDG()->Charge()/3.));
 231   return param;
 232 }
 233
 234 Bool_t  AliMCTrackingTestTask::PropagateToPoint(AliExternalTrackParam *param, Double_t *xyz, Double_t mass, Float_t step){
 235   //
 236   // Propagate track to point xyz using
 237   // AliTracker::PropagateTo functionality
 238   //
 239   //  param - track parameters
 240   //  xyz   - position to propagate
 241   //  mass  - particle mass
 242   //  step  - step to be used
 243   Double_t radius=TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
 244   AliTracker::PropagateTrackTo(param, radius+step, mass, step, kTRUE,0.99);
 245   AliTracker::PropagateTrackTo(param, radius+0.5, mass, step*0.1, kTRUE,0.99);
 246   Double_t sxyz[3]={0,0,0};
 247   AliESDVertex vertex(xyz,sxyz);
 248   Bool_t isOK = param->PropagateToDCA(&vertex,AliTracker::GetBz(),10);
 249   return isOK;
 250 }
 251
 252
 253 void  AliMCTrackingTestTask::ProcessMCInfo(){
 254   //
 255   //
 256   //
 257    //
 258   TParticle * particle= new TParticle;
 259   TClonesArray * trefs = new TClonesArray("AliTrackReference");
 260   const Double_t kPcut=0.1;
 261   //
 262   //
 263   // Process tracks
 264   //
 265   Int_t npart = fMCinfo->GetNumberOfTracks();
 266   if (npart==0) return;
 267   Double_t vertex[4]={0,0,0,0};
 268   fMCinfo->GetParticleAndTR(0, particle, trefs);
 269   if (particle){
 270     vertex[0]=particle->Vx();
 271     vertex[1]=particle->Vy();
 272     vertex[2]=particle->Vz();
 273     vertex[3]=particle->R();
 274   }
 275   //
 276   //
 277
 278   for (Int_t ipart=0;ipart<npart;ipart++){
 279     Int_t status = fMCinfo->GetParticleAndTR(ipart, particle, trefs);
 280     if (status<0) continue;
 281     if (!particle) continue;
 282     if (!trefs) continue;
 283     Int_t nref = trefs->GetEntries();
 284     if (nref<5) continue;
 285     AliTrackReference * tpcIn=0;
 286     AliTrackReference * tpcOut=0;
 287     AliTrackReference * trdIn=0;
 288     AliTrackReference * trdOut=0;
 289     AliTrackReference * itsIn=0;
 290     AliTrackReference * itsOut=0;
 291     Double_t rmax=0;
 292     Double_t rmin=1000;
 293     for (Int_t iref=0;iref<nref;iref++){
 294       AliTrackReference * ref = (AliTrackReference*)trefs->At(iref);
 295       if (!ref) continue;
 296
 297       Float_t dir = ref->X()*ref->Px()+ref->Y()*ref->Py();
 298
 299       if (dir<0) break; // oposite direction - looping track - return back
 300       if (ref->P()<kPcut) continue;
 301       if (ref->R()<rmax) break;
 302       //if (ref->R()<rmin)  break;
 303       //TPC
 304       if (ref->DetectorId()==AliTrackReference::kTPC){
 305         if (!tpcIn) {
 306           tpcIn  = ref;
 307         }else{
 308           if (ref->R()>tpcIn->R()) tpcOut = ref;
 309         }
 310       }
 311       //ITS
 312       if (ref->DetectorId()==AliTrackReference::kITS){
 313         if (!itsIn) {
 314           itsIn  = ref;
 315         }else{
 316           if (ref->R()>itsIn->R()) itsOut = ref;
 317         }
 318       }
 319       //TRD
 320       if (ref->DetectorId()==AliTrackReference::kTRD){
 321         if (!trdIn) {
 322           trdIn  = ref;
 323         }else{
 324           if (ref->R()>trdIn->R()) trdOut = ref;
 325         }
 326       }
 327       if (ref->R()<rmin) rmin=ref->R();
 328       if (ref->R()>rmax) rmax=ref->R();
 329     }
 330     if (tpcIn && tpcOut) {
 331       ProcessRefTracker(tpcIn,tpcOut,particle,1);
 332       ProcessRefTracker(tpcIn,tpcOut,particle,3);
 333     }
 334     if (itsIn && itsOut) ProcessRefTracker(itsIn,itsOut,particle,0);
 335     if (trdIn && trdOut) ProcessRefTracker(trdIn,trdOut,particle,2);
 336   }
 337 }
 338
 339
 340
 341 void AliMCTrackingTestTask::ProcessRefTracker(AliTrackReference* refIn,  AliTrackReference* refOut, TParticle*part,Int_t type){
 342   //
 343   // Test propagation from In to out
 344   //
 345   AliExternalTrackParam *param = 0;
 346   AliExternalTrackParam *paramMC = 0;
 347   Double_t xyzIn[3]={refIn->X(),refIn->Y(), refIn->Z()};
 348   Double_t mass = part->GetMass();
 349   Double_t step=1;
 350   //
 351   param=MakeTrack(refOut,part);
 352   paramMC=MakeTrack(refOut,part);
 353   if (!param) return;
 354   if (type<3) PropagateToPoint(param,xyzIn, mass, step);
 355   if (type==3) {
 356     AliTPCseed seed;
 357     seed.Set(param->GetX(),param->GetAlpha(),param->GetParameter(),param->GetCovariance());
 358     Float_t alpha= TMath::ATan2(refIn->Y(),refIn->X());
 359     seed.Rotate(alpha-seed.GetAlpha());
 360     seed.SetMass(mass);
 361     for (Float_t xlayer= seed.GetX(); xlayer>refIn->R(); xlayer-=step){
 362       seed.PropagateTo(xlayer);
 363     }
 364     seed.PropagateTo(refIn->R());
 365     param->Set(seed.GetX(),seed.GetAlpha(),seed.GetParameter(),seed.GetCovariance());
 366   }
 367   TTreeSRedirector *pcstream = GetDebugStreamer();
 368   TVectorD gpos(3);
 369   TVectorD gmom(3);
 370   param->GetXYZ(gpos.GetMatrixArray());
 371   param->GetPxPyPz(gmom.GetMatrixArray());
 372   if (pcstream){
 373     (*pcstream)<<"MC"<<
 374       "type="<<type<<
 375       "step="<<step<<
 376       "refIn.="<<refIn<<
 377       "refOut.="<<refOut<<
 378       "p.="<<part<<
 379       "par.="<<param<<
 380       "parMC.="<<paramMC<<
 381       "gpos.="<<&gpos<<
 382       "gmom.="<<&gmom<<
 383       "\n";
 384   }
 385 }
 386
 387
 388
 389 void AliMCTrackingTestTask::FinishTaskOutput()
 390 {
 391   //
 392   // According description in AliAnalisysTask this method is call
 393   // on the slaves before sending data
 394   //
 395   Terminate("slave");
 396   gSystem->Exec("pwd");
 397   RegisterDebugOutput();
 398
 399 }
 400
 401
 402 void AliMCTrackingTestTask::RegisterDebugOutput(){
 403   //
 404   //
 405   //
 406   //
 407   // store  - copy debug output to the destination position
 408   // currently ONLY for local copy
 409   TString dsName;
 410   dsName=GetName();
 411   dsName+="Debug.root";
 412   dsName.ReplaceAll(" ","");
 413   TString dsName2=fDebugOutputPath.Data();
 414   gSystem->MakeDirectory(dsName2.Data());
 415   dsName2+=gSystem->HostName();
 416   gSystem->MakeDirectory(dsName2.Data());
 417   dsName2+="/";
 418   dsName2+=gSystem->BaseName(gSystem->pwd());
 419   dsName2+="/";
 420   gSystem->MakeDirectory(dsName2.Data());
 421   dsName2+=dsName;
 422   AliInfo(Form("copy %s\t%s\n",dsName.Data(),dsName2.Data()));
 423   printf("copy %s\t%s\n",dsName.Data(),dsName2.Data());
 424   TFile::Cp(dsName.Data(),dsName2.Data());
 425 }
 426
 427
 428 /*
 429   gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
 430   gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
 431   AliXRDPROOFtoolkit tool;
 432   TChain * chain = tool.MakeChain("mctracking.txt","MC",0,100);
 433   chain->Lookup();
 434   //
 435   //
 436   chain->SetAlias("pdg","(p.fPdgCode)");
 437   chain->SetAlias("dPRec","(refOut.P()-par.P())/refIn.P()");
 438   chain->SetAlias("dPMC","(refOut.P()-refIn.P())/refIn->P()");
 439   chain->SetAlias("dPtRec","(refOut.Pt()-par.Pt())/refIn.Pt()");
 440   chain->SetAlias("dPtMC","(refOut.Pt()-refIn.Pt())/refIn->Pt()");
 441
 442
 443   // ITS
 444   chain->Draw("-sqrt(dPRec):-sqrt(dPMC)","abs(pdg)!=11&&type==0&&p.Pt()<0.5&&abs(p.R())<1&&abs(refOut.fZ)<220");
 445   htemp->SetYTitle("#sqrt{#DeltaP_{rec}/P}");
 446   htemp->SetXTitle("#sqrt{#DeltaP_{mc}/P}");
 447   gPad->SaveAs("picLoss/dPcorr_ITS_step1.gif");
 448   gPad->SaveAs("picLoss/dPcorr_ITS_step1.eps");
 449   // TPC
 450   chain->Draw("-sqrt(dPRec):-sqrt(dPMC)","abs(pdg)!=11&&type==1&&p.Pt()<0.5&&abs(p.R())<1&&abs(refOut.fZ)<220");
 451   htemp->SetYTitle("#sqrt{#DeltaP_{rec}/P}");
 452   htemp->SetXTitle("#sqrt{#DeltaP_{mc}/P}");
 453   gPad->SaveAs("picLoss/dPcorr_TPC_step1.gif");
 454   gPad->SaveAs("picLoss/dPcorr_TPC_step1.eps");
 455   //
 456    // TPC
 457   chain->Draw("-sqrt(dPRec):-sqrt(dPMC)","abs(pdg)!=11&&type==3&&p.Pt()<0.5&&abs(p.R())<1&&abs(refOut.fZ)<220");
 458   htemp->SetYTitle("#sqrt{#DeltaP_{rec}/P}");
 459   htemp->SetXTitle("#sqrt{#DeltaP_{mc}/P}");
 460   gPad->SaveAs("picLoss/dPcorr_TPCseed_step1.gif");
 461   gPad->SaveAs("picLoss/dPcorr_TPCseed_step1.eps");
 462
 463
 464   // TRD
 465   chain->Draw("-sqrt(dPRec):-sqrt(dPMC)","abs(pdg)!=11&&type==2&&p.Pt()<0.5&&abs(p.R())<1&&abs(refOut.fZ)<220");
 466   htemp->SetYTitle("#sqrt{#DeltaP_{rec}/P}");
 467   htemp->SetXTitle("#sqrt{#DeltaP_{mc}/P}");
 468   gPad->SaveAs("picLoss/dPcorr_TRD_step1.gif");
 469   gPad->SaveAs("picLoss/dPcorr_TRD_step1.eps");
 470
 471   //
 472   //
 473   //
 474   chain->Draw("(par.Pt()-refIn.Pt())/refIn.Pt()>>his(100,-0.02,0.02)","abs(pdg)!=11&&type==3&&p.Pt()<0.5&&abs(p.R())<1&&abs(refOut.fZ)<220");
 475   his->SetXTitle("(P_{trec}-P_{tmc})/P_{tmc}");
 476   gPad->SaveAs("picLoss/dPtcorr_TPCseed_step1_1D.eps");
 477   gPad->SaveAs("picLoss/dPtcorr_TPCseed_step1_1D.gif");
 478
 479   chain->Draw("(par.P()-refIn.P())/refIn.P()>>his(100,-0.02,0.02)","abs(pdg)!=11&&type==3&&p.Pt()<0.5&&abs(p.R())<1&&abs(refOut.fZ)<220");
 480   his->SetXTitle("(P_{rec}-P_{mc})/P_{mc}");
 481   gPad->SaveAs("picLoss/dPcorr_TPCseed_step1_1D.eps");
 482   gPad->SaveAs("picLoss/dPcorr_TPCseed_step1_1D.gif");
 483 */