[u/mrichter/AliRoot.git] / PWG1 / AliMCTrackingTestTask.cxx

//
// This class is the task to check the 
// Propagation method used in the 
//               1. AliExternalTrackParam 
//               2. AliTracker
//
// Principle - Creates AliExternalTrackParam form 1 Track Refernece - 
//             Propagate it to other
// Magnetic filed and the geomtry has to bec 

//
// ROOT includes
#include <TChain.h>
#include <TMath.h>
#include <TVectorD.h>
#include <TSystem.h>
#include <TFile.h>

// ALIROOT includes
#include <TTreeStream.h>
#include <AliAnalysisManager.h>
#include <AliESDInputHandler.h>
#include "AliStack.h"
#include "AliMCEvent.h"
#include "AliMCEventHandler.h"

#include <AliESD.h>
#include "AliMCTrackingTestTask.h"
#include "AliGenInfoMaker.h"
#include "AliHelix.h"

//
#include "AliMCInfo.h"
#include "AliComparisonObject.h"
#include "AliESDRecInfo.h"
#include "AliTPCParamSR.h"
#include "AliTracker.h"
#include "AliTPCseed.h"

// STL includes
#include <iostream>

using namespace std;

ClassImp(AliMCTrackingTestTask)

//________________________________________________________________________
AliMCTrackingTestTask::AliMCTrackingTestTask() : 
  AliAnalysisTask(), 
  fMCinfo(0),     //! MC event handler
  fESD(0),
  fDebugStreamer(0),
  fStreamLevel(0),
  fDebugLevel(0),
  fDebugOutputPath()
{
  //
  // Default constructor (should not be used)
  //
}

AliMCTrackingTestTask::AliMCTrackingTestTask(const AliMCTrackingTestTask& /*info*/) : 
  AliAnalysisTask(), 
  fMCinfo(0),     //! MC event handler
  fESD(0),
  //
  fDebugStreamer(0),
  fStreamLevel(0),
  fDebugLevel(),
  fDebugOutputPath()
{
  //
  // Default constructor 
  //
}


//________________________________________________________________________
AliMCTrackingTestTask::AliMCTrackingTestTask(const char *name) : 
  AliAnalysisTask(name, "AliMCTrackingTestTask"), 
  fMCinfo(0),     //! MC event handler
  fESD(0),
  fDebugStreamer(0),
  fStreamLevel(0),
  fDebugLevel(0),
  fDebugOutputPath()
{
  //
  // Normal constructor
  //
  // Input slot #0 works with a TChain
  DefineInput(0, TChain::Class());
  // Output slot #0 writes into a TList
  DefineOutput(0, TObjArray::Class());
  //
  //
}

AliMCTrackingTestTask::~AliMCTrackingTestTask(){
  //
  //
  //
  if (fDebugLevel>0)  printf("AliMCTrackingTestTask::~AliMCTrackingTestTask\n");
  if (fDebugStreamer) delete fDebugStreamer;
  fDebugStreamer=0;
}


//________________________________________________________________________
void AliMCTrackingTestTask::ConnectInputData(Option_t *) 
{
  //
  // Connect the input data
  //
  if(fDebugLevel>3)
    cout << "AnalysisTaskTPCCluster::ConnectInputData()" << endl;

  TTree* tree=dynamic_cast<TTree*>(GetInputData(0));
  if (!tree) {
    //Printf("ERROR: Could not read chain from input slot 0");
  }
  else {
    AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
    if (!esdH) {
      //Printf("ERROR: Could not get ESDInputHandler");
    }
    else {
      fESD = esdH->GetEvent();
      //Printf("*** CONNECTED NEW EVENT ****");
    }  
  }
  AliMCEventHandler* mcinfo = (AliMCEventHandler*) (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());  
  mcinfo->SetReadTR(kTRUE);
  
  fMCinfo = mcinfo->MCEvent();


}


//________________________________________________________________________
void AliMCTrackingTestTask::CreateOutputObjects() 
{
  //
  // Connect the output objects
  //
  if(fDebugLevel>3)
    cout << "AnalysisTaskTPCCluster::CreateOutputObjects()" << endl;
  
}


//________________________________________________________________________
void AliMCTrackingTestTask::Exec(Option_t *) {
  //
  // Execute analysis for current event 
  //

  if(fDebugLevel>3)
    cout << "AliMCTrackingTestTask::Exec()" << endl;
    

  // If MC has been connected   

  if (!fMCinfo){
    cout << "Not MC info\n" << endl;
  }else{
    ProcessMCInfo();
    //mcinfo->Print();
    //DumpInfo();
  }
  //
}      


//________________________________________________________________________
void AliMCTrackingTestTask::Terminate(Option_t *) {
    //
    // Terminate loop
    //
  if(fDebugLevel>3)
    printf("AliMCTrackingTestTask: Terminate() \n");  
  //
  if (fDebugLevel>0) printf("AliMCtrackingTestTask::Terminate\n");
  if (fDebugStreamer) delete fDebugStreamer;
  fDebugStreamer = 0;
  return;
}


TTreeSRedirector *AliMCTrackingTestTask::GetDebugStreamer(){
  //
  // Get Debug streamer
  // In case debug streamer not yet initialized and StreamLevel>0 create new one
  //
  if (fStreamLevel==0) return 0;
  if (fDebugStreamer) return fDebugStreamer;
  TString dsName;
  dsName=GetName();
  dsName+="Debug.root";
  dsName.ReplaceAll(" ","");
  fDebugStreamer = new TTreeSRedirector(dsName.Data());
  return fDebugStreamer;
}


AliExternalTrackParam * AliMCTrackingTestTask::MakeTrack(const AliTrackReference* ref, TParticle*part)
{
  //
  // Make track out of the track ref
  // part - TParticle used to determine chargr
  // the covariance matrix - equal 0 - starting from ideal MC position
  Double_t xyz[3]={ref->X(),ref->Y(),ref->Z()};
  Double_t pxyz[3]={ref->Px(),ref->Py(),ref->Pz()};
  Double_t cv[21];
  for (Int_t i=0; i<21;i++) cv[i]=0;
  if (!part->GetPDG()) return 0;
  AliExternalTrackParam * param = new AliExternalTrackParam(xyz,pxyz,cv,TMath::Nint(part->GetPDG()->Charge()/3.));
  return param;
}

Bool_t  AliMCTrackingTestTask::PropagateToPoint(AliExternalTrackParam *param, Double_t *xyz, Double_t mass, Float_t step){
  // 
  // Propagate track to point xyz using 
  // AliTracker::PropagateTo functionality
  //
  //  param - track parameters
  //  xyz   - position to propagate
  //  mass  - particle mass
  //  step  - step to be used
  Double_t radius=TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
  AliTracker::PropagateTrackTo(param, radius+step, mass, step, kTRUE,0.99);
  AliTracker::PropagateTrackTo(param, radius+0.5, mass, step*0.1, kTRUE,0.99);
  Double_t sxyz[3]={0,0,0};
  AliESDVertex vertex(xyz,sxyz);
  Bool_t isOK = param->PropagateToDCA(&vertex,AliTracker::GetBz(),10);
  return isOK;
}


void  AliMCTrackingTestTask::ProcessMCInfo(){
  //
  //
  //
   //
  TParticle * particle= new TParticle;
  TClonesArray * trefs = new TClonesArray("AliTrackReference");
  const Double_t kPcut=0.1;
  //
  //
  // Process tracks
  //
  Int_t npart = fMCinfo->GetNumberOfTracks();
  if (npart==0) return;
  Double_t vertex[4]={0,0,0,0};
  fMCinfo->GetParticleAndTR(0, particle, trefs);
  if (particle){
    vertex[0]=particle->Vx();
    vertex[1]=particle->Vy();
    vertex[2]=particle->Vz();
    vertex[3]=particle->R();
  }
  //
  //

  for (Int_t ipart=0;ipart<npart;ipart++){
    Int_t status = fMCinfo->GetParticleAndTR(ipart, particle, trefs);
    if (status<0) continue;
    if (!particle) continue;
    if (!trefs) continue;
    Int_t nref = trefs->GetEntries();
    if (nref<5) continue;
    AliTrackReference * tpcIn=0;
    AliTrackReference * tpcOut=0;
    AliTrackReference * trdIn=0;
    AliTrackReference * trdOut=0;
    AliTrackReference * itsIn=0;
    AliTrackReference * itsOut=0;
    Double_t rmax=0;
    Double_t rmin=1000;
    for (Int_t iref=0;iref<nref;iref++){
      AliTrackReference * ref = (AliTrackReference*)trefs->At(iref);
      if (!ref) continue;
      
      Float_t dir = ref->X()*ref->Px()+ref->Y()*ref->Py();
      
      if (dir<0) break; // oposite direction - looping track - return back
      if (ref->P()<kPcut) continue;
      if (ref->R()<rmax) break;
      //if (ref->R()<rmin)  break; 
      //TPC
      if (ref->DetectorId()==AliTrackReference::kTPC){
	if (!tpcIn) {
	  tpcIn  = ref;
	}else{
	  if (ref->R()>tpcIn->R()) tpcOut = ref;
	}	
      }
      //ITS
      if (ref->DetectorId()==AliTrackReference::kITS){
	if (!itsIn) {
	  itsIn  = ref;
	}else{
	  if (ref->R()>itsIn->R()) itsOut = ref;
	}	
      }
      //TRD
      if (ref->DetectorId()==AliTrackReference::kTRD){
	if (!trdIn) {
	  trdIn  = ref;
	}else{
	  if (ref->R()>trdIn->R()) trdOut = ref;
	}	
      }      
      if (ref->R()<rmin) rmin=ref->R();
      if (ref->R()>rmax) rmax=ref->R();
    }
    if (tpcIn && tpcOut) {
      ProcessRefTracker(tpcIn,tpcOut,particle,1);
      ProcessRefTracker(tpcIn,tpcOut,particle,3);
    }
    if (itsIn && itsOut) ProcessRefTracker(itsIn,itsOut,particle,0);
    if (trdIn && trdOut) ProcessRefTracker(trdIn,trdOut,particle,2);
  }
}


void AliMCTrackingTestTask::ProcessRefTracker(AliTrackReference* refIn,  AliTrackReference* refOut, TParticle*part,Int_t type){
  //
  // Test propagation from In to out
  //
  AliExternalTrackParam *param = 0;
  AliExternalTrackParam *paramMC = 0;
  Double_t xyzIn[3]={refIn->X(),refIn->Y(), refIn->Z()};
  Double_t mass = part->GetMass();
  Double_t step=1;
  //
  param=MakeTrack(refOut,part);
  paramMC=MakeTrack(refOut,part);
  if (!param) return;
  if (type<3) PropagateToPoint(param,xyzIn, mass, step);
  if (type==3) {
    AliTPCseed seed;
    seed.Set(param->GetX(),param->GetAlpha(),param->GetParameter(),param->GetCovariance());
    Float_t alpha= TMath::ATan2(refIn->Y(),refIn->X());
    seed.Rotate(alpha-seed.GetAlpha());
    seed.SetMass(mass);
    for (Float_t xlayer= seed.GetX(); xlayer>refIn->R(); xlayer-=step){
      seed.PropagateTo(xlayer);
    }
    seed.PropagateTo(refIn->R());
    param->Set(seed.GetX(),seed.GetAlpha(),seed.GetParameter(),seed.GetCovariance());
  }
  TTreeSRedirector *pcstream = GetDebugStreamer();
  TVectorD gpos(3);
  TVectorD gmom(3);
  param->GetXYZ(gpos.GetMatrixArray());
  param->GetPxPyPz(gmom.GetMatrixArray());
  if (pcstream){
    (*pcstream)<<"MC"<<
      "type="<<type<<
      "step="<<step<<
      "refIn.="<<refIn<<
      "refOut.="<<refOut<<
      "p.="<<part<<
      "par.="<<param<<   
      "parMC.="<<paramMC<<   
      "gpos.="<<&gpos<<
      "gmom.="<<&gmom<<
      "\n";
  }
}


void AliMCTrackingTestTask::FinishTaskOutput()
{
  //
  // According description in AliAnalisysTask this method is call
  // on the slaves before sending data
  //
  Terminate("slave");
  gSystem->Exec("pwd");
  RegisterDebugOutput();

}


void AliMCTrackingTestTask::RegisterDebugOutput(){
  //
  //
  //
  //
  // store  - copy debug output to the destination position
  // currently ONLY for local copy
  TString dsName;
  dsName=GetName();
  dsName+="Debug.root";
  dsName.ReplaceAll(" ","");
  TString dsName2=fDebugOutputPath.Data();
  gSystem->MakeDirectory(dsName2.Data());
  dsName2+=gSystem->HostName();
  gSystem->MakeDirectory(dsName2.Data());
  dsName2+="/";
  dsName2+=gSystem->BaseName(gSystem->pwd());
  dsName2+="/";
  gSystem->MakeDirectory(dsName2.Data());
  dsName2+=dsName;
  AliInfo(Form("copy %s\t%s\n",dsName.Data(),dsName2.Data()));
  printf("copy %s\t%s\n",dsName.Data(),dsName2.Data());
  TFile::Cp(dsName.Data(),dsName2.Data());
}


/*
  gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
  gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
  AliXRDPROOFtoolkit tool; 
  TChain * chain = tool.MakeChain("mctracking.txt","MC",0,100);
  chain->Lookup();
  //
  //
  chain->SetAlias("pdg","(p.fPdgCode)");
  chain->SetAlias("dPRec","(refOut.P()-par.P())/refIn.P()");
  chain->SetAlias("dPMC","(refOut.P()-refIn.P())/refIn->P()");
  chain->SetAlias("dPtRec","(refOut.Pt()-par.Pt())/refIn.Pt()");
  chain->SetAlias("dPtMC","(refOut.Pt()-refIn.Pt())/refIn->Pt()");


  // ITS
  chain->Draw("-sqrt(dPRec):-sqrt(dPMC)","abs(pdg)!=11&&type==0&&p.Pt()<0.5&&abs(p.R())<1&&abs(refOut.fZ)<220");
  htemp->SetYTitle("#sqrt{#DeltaP_{rec}/P}");
  htemp->SetXTitle("#sqrt{#DeltaP_{mc}/P}");
  gPad->SaveAs("picLoss/dPcorr_ITS_step1.gif");
  gPad->SaveAs("picLoss/dPcorr_ITS_step1.eps");
  // TPC
  chain->Draw("-sqrt(dPRec):-sqrt(dPMC)","abs(pdg)!=11&&type==1&&p.Pt()<0.5&&abs(p.R())<1&&abs(refOut.fZ)<220");
  htemp->SetYTitle("#sqrt{#DeltaP_{rec}/P}");
  htemp->SetXTitle("#sqrt{#DeltaP_{mc}/P}");
  gPad->SaveAs("picLoss/dPcorr_TPC_step1.gif");
  gPad->SaveAs("picLoss/dPcorr_TPC_step1.eps");
  //
   // TPC
  chain->Draw("-sqrt(dPRec):-sqrt(dPMC)","abs(pdg)!=11&&type==3&&p.Pt()<0.5&&abs(p.R())<1&&abs(refOut.fZ)<220");
  htemp->SetYTitle("#sqrt{#DeltaP_{rec}/P}");
  htemp->SetXTitle("#sqrt{#DeltaP_{mc}/P}");
  gPad->SaveAs("picLoss/dPcorr_TPCseed_step1.gif");
  gPad->SaveAs("picLoss/dPcorr_TPCseed_step1.eps");


  // TRD
  chain->Draw("-sqrt(dPRec):-sqrt(dPMC)","abs(pdg)!=11&&type==2&&p.Pt()<0.5&&abs(p.R())<1&&abs(refOut.fZ)<220");
  htemp->SetYTitle("#sqrt{#DeltaP_{rec}/P}");
  htemp->SetXTitle("#sqrt{#DeltaP_{mc}/P}");
  gPad->SaveAs("picLoss/dPcorr_TRD_step1.gif");
  gPad->SaveAs("picLoss/dPcorr_TRD_step1.eps");

  //
  //
  //
  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");
  his->SetXTitle("(P_{trec}-P_{tmc})/P_{tmc}");
  gPad->SaveAs("picLoss/dPtcorr_TPCseed_step1_1D.eps");
  gPad->SaveAs("picLoss/dPtcorr_TPCseed_step1_1D.gif");

  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");
  his->SetXTitle("(P_{rec}-P_{mc})/P_{mc}");
  gPad->SaveAs("picLoss/dPcorr_TPCseed_step1_1D.eps");
  gPad->SaveAs("picLoss/dPcorr_TPCseed_step1_1D.gif");
*/
Commit	Line	Data
3880629c	1	//
2a4194b9	2	// This class is the task to check the
	3	// Propagation method used in the
	4	// 1. AliExternalTrackParam
	5	// 2. AliTracker
3880629c	6	//
2a4194b9	7	// Principle - Creates AliExternalTrackParam form 1 Track Refernece -
	8	// Propagate it to other
	9	// Magnetic filed and the geomtry has to bec
3880629c	10
2a4194b9	11	//
3880629c	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"
2a4194b9	38	#include "AliTPCseed.h"
3880629c	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]);
2a4194b9	244	AliTracker::PropagateTrackTo(param, radius+step, mass, step, kTRUE,0.99);
2a4194b9	245	AliTracker::PropagateTrackTo(param, radius+0.5, mass, step*0.1, kTRUE,0.99);
3880629c	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	}
2a4194b9	330	if (tpcIn && tpcOut) {
	331	ProcessRefTracker(tpcIn,tpcOut,particle,1);
	332	ProcessRefTracker(tpcIn,tpcOut,particle,3);
	333	}
3880629c	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;
2a4194b9	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	}
3880629c	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()");
2a4194b9	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");
3880629c	483	*/