PWG1/AliMCTrackingTestTask.cxx

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