[u/mrichter/AliRoot.git] / PWGPP / AliAnalysisTaskIPInfo.cxx

/**************************************************************************
 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

//*************************************************************************
// Class AliAnalysisTaskIPInfo
// AliAnalysisTask to extract from the ESD the IP position and sigma
// as well as to estimate the primary vertex and tracks DCA resolution.
// Uses external class AliIntSpotEstimator
//
// Author: ruben.shahoyan@cern.ch
//*************************************************************************

#include <TChain.h>
#include <TTree.h>
#include <TNtuple.h>
#include <TBranch.h>
#include <TClonesArray.h>
#include <TObjArray.h>
#include <TH1F.h>
#include <TH2F.h>  
#include <TNtuple.h>  
#include <TCanvas.h>

#include "AliAnalysisTask.h"
#include "AliAnalysisManager.h"

#include "AliESDtrack.h"
#include "AliExternalTrackParam.h"
#include "AliESDVertex.h"
#include "AliESDEvent.h"
#include "AliESDfriend.h"
#include "AliVertexerTracks.h"
#include "AliESDInputHandler.h"
#include "AliAnalysisTaskIPInfo.h"
#include "AliIntSpotEstimator.h"
#include "AliMultiplicity.h"


ClassImp(AliAnalysisTaskIPInfo)

const Char_t* AliAnalysisTaskIPInfo::fEstNames[kNEst] = {"ITSTPC","TPC","SPD"};

//________________________________________________________________________
AliAnalysisTaskIPInfo::AliAnalysisTaskIPInfo(const char *name) 
: AliAnalysisTask(name, "IP analysis"),
  fESD(0),fESDfriend(0),fOutput(0),fTracks(50)
{

  // Define input and output slots here
  // Input slot #0 works with a TChain
  DefineInput(0, TChain::Class());
  // Output slot #0 writes into a TList container
  DefineOutput(0, TList::Class());  //My private output
  //
  for (int i=0;i<kNEst;i++) fIPEst[i] = 0;
  //
}

//________________________________________________________________________
AliAnalysisTaskIPInfo::~AliAnalysisTaskIPInfo()
{
  // Destructor
  if (fOutput) {
    delete fOutput;
    fOutput = 0;
  }
}
//________________________________________________________________________
void AliAnalysisTaskIPInfo::SetIPCenIni(Int_t estID, Double_t x,Double_t y,Double_t z)
{
  // set initial estimate of the IP center
  if (estID<0 || estID>= kNEst) return;
  fIPCenIni[estID][0] = x;
  fIPCenIni[estID][1] = y;
  fIPCenIni[estID][2] = z;
}

//________________________________________________________________________
void AliAnalysisTaskIPInfo::SetOptions(Int_t estID, Bool_t recoVtx,
				       Double_t outcut,Int_t ntrIP,Int_t nPhiBins,Int_t nestb,
				       Double_t estmin,Double_t estmax,
				       Int_t ntrBins,Int_t ntMn,Int_t ntMx,
				       Int_t nPBins,Double_t pmn,Double_t pmx,Bool_t fillNt)
{
  // set options for estimators
  if (estID<0 || estID>= kNEst) return;
  fNTrMinIP[estID] = ntrIP;
  fRecoVtx[estID]  = recoVtx;
  fNPhiBins[estID] = nPhiBins;
  fNEstb[estID]    = nestb;
  fNTrBins[estID]  = ntrBins;
  fNPBins[estID]   = nPBins;
  fNTrMin[estID]   = ntMn;
  fNTrMax[estID]   = ntMx;
  fOutCut[estID]   = outcut;
  fEstMin[estID]   = estmin;
  fEstMax[estID]   = estmax;
  fPMin[estID]     = pmn;
  fPMax[estID]     = pmx;
  fFillNt[estID]   = fillNt;
}

//________________________________________________________________________
void AliAnalysisTaskIPInfo::ConnectInputData(Option_t *) 
{
  // Connect ESD or AOD here
  // Called once
  //
  AliInfo("HERE");
  TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
  if (!tree) {
    Printf("ERROR: Could not read chain from input slot 0");
  } 
  else {
    tree->SetBranchAddress("ESDfriend.",&fESDfriend);
    AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());    
    if (!esdH) Printf("ERROR: Could not get ESDInputHandler");
    else fESD = esdH->GetEvent();
  }
  //
  return;
}

//________________________________________________________________________
void AliAnalysisTaskIPInfo::CreateOutputObjects()
{
  // Create estimators
  // Several histograms are more conveniently managed in a TList
  fOutput = new TList;
  fOutput->SetOwner();
  //
  Double_t err[3]={0.0400,0.0400,7.5};
  //
  for (int i=0;i<kNEst;i++) if (fNEstb[i]>1) {
      TString nm   = GetName();
      TString nmes = fEstNames[i];
      nm += nmes;
      fIPEst[i] = new AliIntSpotEstimator(nm.Data(),fOutCut[i],fNTrMinIP[i],fNPhiBins[i],
					  fNEstb[i],fEstMin[i],fEstMax[i],
					 fNTrBins[i],fNTrMin[i],fNTrMax[i],
					 fNPBins[i],fPMin[i],fPMax[i],fFillNt[i]);
      AliESDVertex *initVertex = new AliESDVertex(fIPCenIni[i],err);
      fIPEst[i]->GetVertexer()->SetVtxStart(initVertex);
      delete initVertex;
      fIPEst[i]->GetVertexer()->SetConstraintOff();
      fIPEst[i]->GetVertexer()->SetMinClusters(2);
      fIPEst[i]->SetIPCenIni(fIPCenIni[i]);
      if (nmes == "TPC") fIPEst[i]->GetVertexer()->SetTPCMode();
      else               fIPEst[i]->GetVertexer()->SetITSMode();
      //
      fOutput->Add(fIPEst[i]);
      if (fIPEst[i]->GetNtuple()) fOutput->Add(fIPEst[i]->GetNtuple());
    }
  //
  return;
}

//________________________________________________________________________
void AliAnalysisTaskIPInfo::Exec(Option_t *) 
{
  static TClonesArray tracks("AliExternalTrackParam",50);
  //
  // Main loop
  // Called for each event
  if (!fESD) {
    Printf("ERROR: fESD not available");
    return;
  }
  fESD->SetESDfriend(fESDfriend);
  //
  const AliESDVertex *vtx;
  UShort_t *trackID;
  Int_t ntracks;
  //
  for (int ie=0;ie<kNEst;ie++) {
    //
    if (!fIPEst[ie]) continue;
    if (ie==kTPC) {
      fIPEst[kTPC]->GetVertexer()->SetFieldkG( fESD->GetMagneticField() );
      vtx = fRecoVtx[kTPC] ? fIPEst[kTPC]->GetVertexer()->FindPrimaryVertex(fESD) : fESD->GetPrimaryVertexTPC();
      if (vtx) {
	ntracks = vtx->GetNIndices();
	trackID = (UShort_t*)vtx->GetIndices();
	for (int i=ntracks;i--;) fTracks.Add((TObject*)fESD->GetTrack(trackID[i])->GetTPCInnerParam());
	fIPEst[kTPC]->ProcessEvent(&fTracks);
	fTracks.Clear();
      }
    }
    else if (ie==kITSTPC) {
      fIPEst[kITSTPC]->GetVertexer()->SetFieldkG( fESD->GetMagneticField() );
      vtx = fRecoVtx[kITSTPC] ? fIPEst[kITSTPC]->GetVertexer()->FindPrimaryVertex(fESD) : fESD->GetPrimaryVertex();
      if (vtx) {
	ntracks = vtx->GetNIndices();
	trackID = (UShort_t*)vtx->GetIndices();
	for (int i=ntracks;i--;) fTracks.Add((TObject*)fESD->GetTrack(trackID[i]));
	fIPEst[kITSTPC]->ProcessEvent(&fTracks);
	fTracks.Clear();
      }
    }
    else if (ie==kSPD) {
      fIPEst[kSPD]->GetVertexer()->SetFieldkG( fESD->GetMagneticField() );
      ntracks = CreateSPDTracklets(tracks);
      for (int i=ntracks;i--;) fTracks.Add((TObject*)tracks[i]);
      fIPEst[kSPD]->ProcessEvent(&fTracks);
      fTracks.Clear();
    }
  }
  //
  PostData(0, fOutput);
  //
  return;
}      

//________________________________________________________________________
Int_t AliAnalysisTaskIPInfo::CreateSPDTracklets(TClonesArray& tracks)
{
  // create traclets from multiplicity class
  double cv[21] = {
    25e-4,
    0   ,  25e-4,
    0   ,     0,  40e-2,
    0   ,     0,      0,   1e-2,
    0   ,     0,      0,      0,   1e-2,
    0   ,     0,      0,      0,      0,   1e-2
  };
  //
  double xyzt[3],pxyz[3];
  int nSPDtracks = 0;
  tracks.Delete();
  const AliMultiplicity *mult = fESD->GetMultiplicity();
  const AliESDVertex *spdVtx  = fESD->GetPrimaryVertexSPD();
  int nTracklets = 0;
  if (mult && spdVtx && (nTracklets=mult->GetNumberOfTracklets())>2 ) {
    const Double_t kRLay1=3.9, kRLay2=7.6;
    double xv = spdVtx->GetX();
    double yv = spdVtx->GetY();
    double zv = spdVtx->GetZ();
    for (int i=0;i<nTracklets;i++) { // get cluster coordinates from tracklet
      double phi1 = mult->GetPhi(i);
      double tht1 = mult->GetTheta(i);
      double phi2 = phi1 - mult->GetDeltaPhi(i);
      double tht2 = tht1 - mult->GetDeltaTheta(i);
      double cs = TMath::Cos(phi1);
      double sn = TMath::Sin(tht1);
      double det = xv*sn+yv*sn;
      det = det*det + kRLay1*kRLay1;
      double t  = -(xv*cs+yv*sn) + TMath::Sqrt(det);
      double x1 = cs*t;
      double y1 = sn*t;
      double z1 = zv + TMath::Sqrt(x1*x1+y1*y1)/TMath::Tan(tht1);
      x1 += xv;
      y1 += yv;
      //
      cs = TMath::Cos(phi2);
      sn = TMath::Sin(tht2);
      det = xv*sn+yv*sn;
      det = det*det + kRLay2*kRLay2;
      t  = -(xv*cs+yv*sn) + TMath::Sqrt(det);
      double dx = cs*t;
      double dy = sn*t;
      double dz = zv + TMath::Sqrt(dx*dx+dy*dy)/TMath::Tan(tht2);
      dx += xv-x1;
      dy += yv-y1;
      dz += -z1;
      double dr = TMath::Sqrt(dx*dx+dy*dy+dz*dz);
      pxyz[0] = dx/dr; // direction cosines
      pxyz[1] = dy/dr;
      pxyz[2] = dz/dr;
      t = (xv-x1)*pxyz[0] + (yv-y1)*pxyz[1] + (zv-z1)*pxyz[2];
      xyzt[0] = x1 + t*pxyz[0];  // PCA to vertex
      xyzt[1] = y1 + t*pxyz[1];
      xyzt[2] = z1 + t*pxyz[2];
      //
      new(tracks[nSPDtracks++]) AliExternalTrackParam(xyzt,pxyz,cv,0);
    }
  }
  return nSPDtracks;
}

//________________________________________________________________________
void AliAnalysisTaskIPInfo::Terminate(Option_t *) 
{
  // Draw result to the screen
  // Called once at the end of the query
  fOutput = dynamic_cast<TList*> (GetOutputData(0));
  if (!fOutput) {     
    Printf("ERROR: fOutput not available");
    return;
  }
  //
  return;
}
Commit	Line	Data
5b09c01f	1	/**************************************************************************
	2	* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	//*************************************************************************
	17	// Class AliAnalysisTaskIPInfo
	18	// AliAnalysisTask to extract from the ESD the IP position and sigma
	19	// as well as to estimate the primary vertex and tracks DCA resolution.
	20	// Uses external class AliIntSpotEstimator
	21	//
	22	// Author: ruben.shahoyan@cern.ch
	23	//*************************************************************************
	24
	25	#include <TChain.h>
	26	#include <TTree.h>
	27	#include <TNtuple.h>
	28	#include <TBranch.h>
	29	#include <TClonesArray.h>
	30	#include <TObjArray.h>
	31	#include <TH1F.h>
	32	#include <TH2F.h>
d631e5e7	33	#include <TNtuple.h>
5b09c01f	34	#include <TCanvas.h>
	35
	36	#include "AliAnalysisTask.h"
	37	#include "AliAnalysisManager.h"
	38
	39	#include "AliESDtrack.h"
	40	#include "AliExternalTrackParam.h"
	41	#include "AliESDVertex.h"
	42	#include "AliESDEvent.h"
d631e5e7	43	#include "AliESDfriend.h"
5b09c01f	44	#include "AliVertexerTracks.h"
	45	#include "AliESDInputHandler.h"
	46	#include "AliAnalysisTaskIPInfo.h"
	47	#include "AliIntSpotEstimator.h"
d631e5e7	48	#include "AliMultiplicity.h"
5b09c01f	49
	50
	51	ClassImp(AliAnalysisTaskIPInfo)
	52
d631e5e7	53	const Char_t* AliAnalysisTaskIPInfo::fEstNames[kNEst] = {"ITSTPC","TPC","SPD"};
d631e5e7	54
5b09c01f	55	//________________________________________________________________________
	56	AliAnalysisTaskIPInfo::AliAnalysisTaskIPInfo(const char *name)
	57	: AliAnalysisTask(name, "IP analysis"),
d631e5e7	58	fESD(0),fESDfriend(0),fOutput(0),fTracks(50)
5b09c01f	59	{
	60
	61	// Define input and output slots here
	62	// Input slot #0 works with a TChain
	63	DefineInput(0, TChain::Class());
	64	// Output slot #0 writes into a TList container
	65	DefineOutput(0, TList::Class()); //My private output
	66	//
d631e5e7	67	for (int i=0;i<kNEst;i++) fIPEst[i] = 0;
d631e5e7	68	//
5b09c01f	69	}
	70
	71	//________________________________________________________________________
	72	AliAnalysisTaskIPInfo::~AliAnalysisTaskIPInfo()
	73	{
	74	// Destructor
	75	if (fOutput) {
	76	delete fOutput;
	77	fOutput = 0;
	78	}
	79	}
d631e5e7	80	//________________________________________________________________________
	81	void AliAnalysisTaskIPInfo::SetIPCenIni(Int_t estID, Double_t x,Double_t y,Double_t z)
	82	{
	83	// set initial estimate of the IP center
	84	if (estID<0 \|\| estID>= kNEst) return;
	85	fIPCenIni[estID][0] = x;
	86	fIPCenIni[estID][1] = y;
	87	fIPCenIni[estID][2] = z;
	88	}
5b09c01f	89
	90	//________________________________________________________________________
	91	void AliAnalysisTaskIPInfo::SetOptions(Int_t estID, Bool_t recoVtx,
d631e5e7	92	Double_t outcut,Int_t ntrIP,Int_t nPhiBins,Int_t nestb,
5b09c01f	93	Double_t estmin,Double_t estmax,
5b09c01f	94	Int_t ntrBins,Int_t ntMn,Int_t ntMx,
d631e5e7	95	Int_t nPBins,Double_t pmn,Double_t pmx,Bool_t fillNt)
5b09c01f	96	{
	97	// set options for estimators
	98	if (estID<0 \|\| estID>= kNEst) return;
d631e5e7	99	fNTrMinIP[estID] = ntrIP;
d631e5e7	100	fRecoVtx[estID] = recoVtx;
5b09c01f	101	fNPhiBins[estID] = nPhiBins;
	102	fNEstb[estID] = nestb;
	103	fNTrBins[estID] = ntrBins;
	104	fNPBins[estID] = nPBins;
	105	fNTrMin[estID] = ntMn;
	106	fNTrMax[estID] = ntMx;
	107	fOutCut[estID] = outcut;
	108	fEstMin[estID] = estmin;
	109	fEstMax[estID] = estmax;
	110	fPMin[estID] = pmn;
	111	fPMax[estID] = pmx;
d631e5e7	112	fFillNt[estID] = fillNt;
5b09c01f	113	}
	114
	115	//________________________________________________________________________
	116	void AliAnalysisTaskIPInfo::ConnectInputData(Option_t *)
	117	{
	118	// Connect ESD or AOD here
	119	// Called once
	120	//
d631e5e7	121	AliInfo("HERE");
5b09c01f	122	TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
	123	if (!tree) {
	124	Printf("ERROR: Could not read chain from input slot 0");
	125	}
	126	else {
d631e5e7	127	tree->SetBranchAddress("ESDfriend.",&fESDfriend);
5b09c01f	128	AliESDInputHandler esdH = dynamic_cast<AliESDInputHandler> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
	129	if (!esdH) Printf("ERROR: Could not get ESDInputHandler");
	130	else fESD = esdH->GetEvent();
	131	}
	132	//
	133	return;
	134	}
	135
	136	//________________________________________________________________________
	137	void AliAnalysisTaskIPInfo::CreateOutputObjects()
	138	{
	139	// Create estimators
	140	// Several histograms are more conveniently managed in a TList
	141	fOutput = new TList;
	142	fOutput->SetOwner();
	143	//
d631e5e7	144	Double_t err[3]={0.0400,0.0400,7.5};
5b09c01f	145	//
d631e5e7	146	for (int i=0;i<kNEst;i++) if (fNEstb[i]>1) {
	147	TString nm = GetName();
	148	TString nmes = fEstNames[i];
	149	nm += nmes;
	150	fIPEst[i] = new AliIntSpotEstimator(nm.Data(),fOutCut[i],fNTrMinIP[i],fNPhiBins[i],
	151	fNEstb[i],fEstMin[i],fEstMax[i],
	152	fNTrBins[i],fNTrMin[i],fNTrMax[i],
	153	fNPBins[i],fPMin[i],fPMax[i],fFillNt[i]);
	154	AliESDVertex *initVertex = new AliESDVertex(fIPCenIni[i],err);
	155	fIPEst[i]->GetVertexer()->SetVtxStart(initVertex);
	156	delete initVertex;
	157	fIPEst[i]->GetVertexer()->SetConstraintOff();
	158	fIPEst[i]->GetVertexer()->SetMinClusters(2);
	159	fIPEst[i]->SetIPCenIni(fIPCenIni[i]);
	160	if (nmes == "TPC") fIPEst[i]->GetVertexer()->SetTPCMode();
	161	else fIPEst[i]->GetVertexer()->SetITSMode();
	162	//
	163	fOutput->Add(fIPEst[i]);
	164	if (fIPEst[i]->GetNtuple()) fOutput->Add(fIPEst[i]->GetNtuple());
	165	}
5b09c01f	166	//
	167	return;
	168	}
	169
	170	//________________________________________________________________________
	171	void AliAnalysisTaskIPInfo::Exec(Option_t *)
	172	{
d631e5e7	173	static TClonesArray tracks("AliExternalTrackParam",50);
d631e5e7	174	//
5b09c01f	175	// Main loop
5b09c01f	176	// Called for each event
5b09c01f	177	if (!fESD) {
	178	Printf("ERROR: fESD not available");
	179	return;
	180	}
d631e5e7	181	fESD->SetESDfriend(fESDfriend);
5b09c01f	182	//
5b09c01f	183	const AliESDVertex *vtx;
	184	UShort_t *trackID;
	185	Int_t ntracks;
	186	//
d631e5e7	187	for (int ie=0;ie<kNEst;ie++) {
	188	//
	189	if (!fIPEst[ie]) continue;
	190	if (ie==kTPC) {
	191	fIPEst[kTPC]->GetVertexer()->SetFieldkG( fESD->GetMagneticField() );
	192	vtx = fRecoVtx[kTPC] ? fIPEst[kTPC]->GetVertexer()->FindPrimaryVertex(fESD) : fESD->GetPrimaryVertexTPC();
	193	if (vtx) {
	194	ntracks = vtx->GetNIndices();
	195	trackID = (UShort_t*)vtx->GetIndices();
	196	for (int i=ntracks;i--;) fTracks.Add((TObject*)fESD->GetTrack(trackID[i])->GetTPCInnerParam());
	197	fIPEst[kTPC]->ProcessEvent(&fTracks);
	198	fTracks.Clear();
	199	}
	200	}
	201	else if (ie==kITSTPC) {
	202	fIPEst[kITSTPC]->GetVertexer()->SetFieldkG( fESD->GetMagneticField() );
	203	vtx = fRecoVtx[kITSTPC] ? fIPEst[kITSTPC]->GetVertexer()->FindPrimaryVertex(fESD) : fESD->GetPrimaryVertex();
	204	if (vtx) {
	205	ntracks = vtx->GetNIndices();
	206	trackID = (UShort_t*)vtx->GetIndices();
	207	for (int i=ntracks;i--;) fTracks.Add((TObject*)fESD->GetTrack(trackID[i]));
	208	fIPEst[kITSTPC]->ProcessEvent(&fTracks);
	209	fTracks.Clear();
	210	}
	211	}
	212	else if (ie==kSPD) {
	213	fIPEst[kSPD]->GetVertexer()->SetFieldkG( fESD->GetMagneticField() );
	214	ntracks = CreateSPDTracklets(tracks);
	215	for (int i=ntracks;i--;) fTracks.Add((TObject*)tracks[i]);
	216	fIPEst[kSPD]->ProcessEvent(&fTracks);
	217	fTracks.Clear();
	218	}
5b09c01f	219	}
5b09c01f	220	//
5b09c01f	221	PostData(0, fOutput);
	222	//
	223	return;
	224	}
	225
d631e5e7	226	//________________________________________________________________________
	227	Int_t AliAnalysisTaskIPInfo::CreateSPDTracklets(TClonesArray& tracks)
	228	{
	229	// create traclets from multiplicity class
	230	double cv[21] = {
	231	25e-4,
	232	0 , 25e-4,
	233	0 , 0, 40e-2,
	234	0 , 0, 0, 1e-2,
	235	0 , 0, 0, 0, 1e-2,
	236	0 , 0, 0, 0, 0, 1e-2
	237	};
	238	//
	239	double xyzt[3],pxyz[3];
	240	int nSPDtracks = 0;
	241	tracks.Delete();
	242	const AliMultiplicity *mult = fESD->GetMultiplicity();
	243	const AliESDVertex *spdVtx = fESD->GetPrimaryVertexSPD();
	244	int nTracklets = 0;
	245	if (mult && spdVtx && (nTracklets=mult->GetNumberOfTracklets())>2 ) {
	246	const Double_t kRLay1=3.9, kRLay2=7.6;
e690d4d0	247	double xv = spdVtx->GetX();
	248	double yv = spdVtx->GetY();
	249	double zv = spdVtx->GetZ();
d631e5e7	250	for (int i=0;i<nTracklets;i++) { // get cluster coordinates from tracklet
	251	double phi1 = mult->GetPhi(i);
	252	double tht1 = mult->GetTheta(i);
	253	double phi2 = phi1 - mult->GetDeltaPhi(i);
	254	double tht2 = tht1 - mult->GetDeltaTheta(i);
	255	double cs = TMath::Cos(phi1);
	256	double sn = TMath::Sin(tht1);
	257	double det = xvsn+yvsn;
	258	det = detdet + kRLay1kRLay1;
	259	double t = -(xvcs+yvsn) + TMath::Sqrt(det);
	260	double x1 = cs*t;
	261	double y1 = sn*t;
	262	double z1 = zv + TMath::Sqrt(x1x1+y1y1)/TMath::Tan(tht1);
	263	x1 += xv;
	264	y1 += yv;
	265	//
	266	cs = TMath::Cos(phi2);
	267	sn = TMath::Sin(tht2);
	268	det = xvsn+yvsn;
	269	det = detdet + kRLay2kRLay2;
	270	t = -(xvcs+yvsn) + TMath::Sqrt(det);
	271	double dx = cs*t;
	272	double dy = sn*t;
	273	double dz = zv + TMath::Sqrt(dxdx+dydy)/TMath::Tan(tht2);
	274	dx += xv-x1;
	275	dy += yv-y1;
	276	dz += -z1;
	277	double dr = TMath::Sqrt(dxdx+dydy+dz*dz);
	278	pxyz[0] = dx/dr; // direction cosines
	279	pxyz[1] = dy/dr;
	280	pxyz[2] = dz/dr;
	281	t = (xv-x1)pxyz[0] + (yv-y1)pxyz[1] + (zv-z1)*pxyz[2];
	282	xyzt[0] = x1 + t*pxyz[0]; // PCA to vertex
	283	xyzt[1] = y1 + t*pxyz[1];
	284	xyzt[2] = z1 + t*pxyz[2];
	285	//
	286	new(tracks[nSPDtracks++]) AliExternalTrackParam(xyzt,pxyz,cv,0);
	287	}
	288	}
	289	return nSPDtracks;
	290	}
	291
5b09c01f	292	//________________________________________________________________________
	293	void AliAnalysisTaskIPInfo::Terminate(Option_t *)
	294	{
	295	// Draw result to the screen
	296	// Called once at the end of the query
	297	fOutput = dynamic_cast<TList*> (GetOutputData(0));
	298	if (!fOutput) {
	299	Printf("ERROR: fOutput not available");
	300	return;
	301	}
	302	//
	303	return;
	304	}