[u/mrichter/AliRoot.git] / PWGPP / TPC / AliESDRecKinkInfo.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, 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.                  *
 **************************************************************************/


///////////////////////////////////////////////////////////////////////////////
//                                                                           //
//  Time Projection Chamber                                                  //
//  Comparison macro for reconstructed tracks - ESDs V0s                                     //
//  responsible: 
//  marian.ivanov@cern.ch                                                    //
//
//

 
#include <stdio.h>
#include <string.h>
//ROOT includes
#include "Rtypes.h"
//
//ALIROOT includes
//
#include "AliESDtrack.h"
#include "AliTPCParam.h"
#include "AliTrackReference.h"
#include "AliTPCParamSR.h"
#include "AliESD.h"
#include "AliESDfriend.h"
#include "AliESDtrack.h"
#include "AliTPCseed.h"
#include "AliITStrackMI.h"
#include "AliHelix.h"
#include "AliESDVertex.h"
#include "AliExternalTrackParam.h"
#include "AliESDkink.h"
#include "AliESDv0.h"
#include "AliV0.h"
//
#include "AliTreeDraw.h"
#include "AliMCInfo.h"
#include "AliGenKinkInfo.h"
#include "AliGenV0Info.h"
#include "AliESDRecKinkInfo.h"


ClassImp(AliESDRecKinkInfo)


AliESDRecKinkInfo::AliESDRecKinkInfo():
  fT1(),      //track1
  fT2(),      //track2  
  fKink(),    //kink
  fDist1(0),    //info about closest distance according closest MC - linear DCA
  fDist2(0),    //info about closest distance parabolic DCA
  fInvMass(0),  //reconstructed invariant mass -
  //
  fRr(0),       // rec position of the vertex 
  fMinR(0),     // minimum radius in rphi intersection
  fDistMinR(0), // distance at minimal radius
  fPointAngleFi(0), //point angle fi
  fPointAngleTh(0), //point angle theta
  fPointAngle(0),   //point angle full
  fStatus(0),       //status -tracks 
  fRecStatus(0),    //kink -status- 0 - not found  1-good -  fake
  fMultiple(0),     // how many times was kink reconstructed
  fKinkMultiple(0) // how many times was kink reconstructed
{
  //
  // Default constructor
  //
  for (Int_t i=0; i<3; i++) {
    fPdr[i] = 0;
    fXr[i] = 0;
    fPm[i] = 0;
    fAngle[i] = 0;
  }
  for (Int_t i=0; i<2; i++) { fLab[i] = 0; } 
}

////
void  AliESDRecKinkInfo::Update()
{

  if ( (fT1.fTPCOn)&& (fT2.fTPCOn)){
    //
    // IF BOTH RECONSTRUCTED
    Float_t distance1,distance2;
    Double_t xx[3],pp[3];
    //
    Double_t xd[3],pd[3],signd;
    Double_t xm[3],pm[3],signm;
    for (Int_t i=0;i<3;i++){
      xd[i] = fT2.fTPCinR1[i];
      pd[i] = fT2.fTPCinP1[i];
      xm[i] = fT1.fTPCinR1[i];
      pm[i] = fT1.fTPCinP1[i];
    }
    signd =  fT2.fSign<0 ? -1:1;
    signm =  fT1.fSign<0 ? -1:1;

    AliHelix dhelix1(xd,pd,signd);
    dhelix1.GetMomentum(0,pp,0);
    dhelix1.Evaluate(0,xx);      
    // 
    //  Double_t x2[3],p2[3];
    //            
    AliHelix mhelix(xm,pm,signm);    
    //
    //find intersection linear
    //
    Double_t phase[2][2] = { {0,0}, {0,0}};
    Double_t radius[2] = {0} ;
    Int_t  points = dhelix1.GetRPHIintersections(mhelix, phase, radius,200);
    Double_t delta1=10000,delta2=10000;  

    if (points==1){
      fMinR = TMath::Sqrt(radius[0]);
    }
    if (points==2){
      fMinR =TMath::Min(TMath::Sqrt(radius[0]),TMath::Sqrt(radius[1]));
    }
    
    if (points>0){
      dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
      dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
      dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
    }
    if (points==2){    
      dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
      dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
      dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
    }
    if (points==1){
      fMinR = TMath::Sqrt(radius[0]);
      fDistMinR = delta1;
    }
    if (points==2){
      if (radius[0]<radius[1]){
	fMinR = TMath::Sqrt(radius[0]);
	fDistMinR = delta1;
      }
      else{
	fMinR = TMath::Sqrt(radius[1]);
	fDistMinR = delta2;
      }
    }
    //
    //
    distance1 = TMath::Min(delta1,delta2);
    //
    //find intersection parabolic
    //
    points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
    delta1=10000,delta2=10000;  
    
    if (points>0){
      dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
    }
    if (points==2){    
      dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
    }
    
    distance2 = TMath::Min(delta1,delta2);
    if (delta1<delta2){
      //get V0 info
      dhelix1.Evaluate(phase[0][0],fXr);
      dhelix1.GetMomentum(phase[0][0],fPdr);
      mhelix.GetMomentum(phase[0][1],fPm);
      dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],fAngle);
      fRr = TMath::Sqrt(radius[0]);
    }
    else{
      dhelix1.Evaluate(phase[1][0],fXr);
      dhelix1.GetMomentum(phase[1][0], fPdr);
      mhelix.GetMomentum(phase[1][1], fPm);
      dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],fAngle);
      fRr = TMath::Sqrt(radius[1]);
    }
    fDist1 = TMath::Sqrt(distance1);
    fDist2 = TMath::Sqrt(distance2);      
    
    if (fDist2<10.5){
      Double_t x,alpha,param[5],cov[15];
      //
      fT1.GetESDtrack()->GetInnerExternalParameters(alpha,x,param);
      fT1.GetESDtrack()->GetInnerExternalCovariance(cov);
      AliExternalTrackParam paramm(x,alpha,param,cov);
      //
      fT2.GetESDtrack()->GetInnerExternalParameters(alpha,x,param);
      fT2.GetESDtrack()->GetInnerExternalCovariance(cov);
      AliExternalTrackParam paramd(x,alpha,param,cov);
      /*
      AliESDkink kink;
      kink.Update(&paramm,&paramd);
      //      kink.Dump();
      Double_t diff  = kink.fRr-fRr;
      Double_t diff2 = kink.fDist2-fDist2;
      printf("Diff\t%f\t%f\n",diff,diff2);
      */
    }
    
    //            
    //
  }

}
Commit	Line	Data
7cc34f08	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, 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	///////////////////////////////////////////////////////////////////////////////
	18	// //
	19	// Time Projection Chamber //
	20	// Comparison macro for reconstructed tracks - ESDs V0s //
	21	// responsible:
	22	// marian.ivanov@cern.ch //
	23	//
	24	//
	25
	26
	27
	28
	29
	30	#include <stdio.h>
	31	#include <string.h>
	32	//ROOT includes
	33	#include "Rtypes.h"
	34	//
	35	//ALIROOT includes
	36	//
	37	#include "AliESDtrack.h"
	38	#include "AliTPCParam.h"
	39	#include "AliTrackReference.h"
	40	#include "AliTPCParamSR.h"
	41	#include "AliESD.h"
	42	#include "AliESDfriend.h"
	43	#include "AliESDtrack.h"
	44	#include "AliTPCseed.h"
	45	#include "AliITStrackMI.h"
7cc34f08	46	#include "AliHelix.h"
	47	#include "AliESDVertex.h"
	48	#include "AliExternalTrackParam.h"
	49	#include "AliESDkink.h"
	50	#include "AliESDv0.h"
	51	#include "AliV0.h"
	52	//
	53	#include "AliTreeDraw.h"
	54	#include "AliMCInfo.h"
	55	#include "AliGenKinkInfo.h"
	56	#include "AliGenV0Info.h"
	57	#include "AliESDRecKinkInfo.h"
	58
	59
	60
	61	ClassImp(AliESDRecKinkInfo)
	62
	63
	64
	65	AliESDRecKinkInfo::AliESDRecKinkInfo():
	66	fT1(), //track1
	67	fT2(), //track2
	68	fKink(), //kink
	69	fDist1(0), //info about closest distance according closest MC - linear DCA
	70	fDist2(0), //info about closest distance parabolic DCA
	71	fInvMass(0), //reconstructed invariant mass -
	72	//
	73	fRr(0), // rec position of the vertex
	74	fMinR(0), // minimum radius in rphi intersection
	75	fDistMinR(0), // distance at minimal radius
	76	fPointAngleFi(0), //point angle fi
	77	fPointAngleTh(0), //point angle theta
	78	fPointAngle(0), //point angle full
	79	fStatus(0), //status -tracks
	80	fRecStatus(0), //kink -status- 0 - not found 1-good - fake
	81	fMultiple(0), // how many times was kink reconstructed
	82	fKinkMultiple(0) // how many times was kink reconstructed
	83	{
	84	//
	85	// Default constructor
	86	//
4a9220d4	87	for (Int_t i=0; i<3; i++) {
	88	fPdr[i] = 0;
	89	fXr[i] = 0;
	90	fPm[i] = 0;
	91	fAngle[i] = 0;
	92	}
	93	for (Int_t i=0; i<2; i++) { fLab[i] = 0; }
7cc34f08	94	}
	95
	96	////
	97	void AliESDRecKinkInfo::Update()
	98	{
	99
	100	if ( (fT1.fTPCOn)&& (fT2.fTPCOn)){
	101	//
	102	// IF BOTH RECONSTRUCTED
	103	Float_t distance1,distance2;
	104	Double_t xx[3],pp[3];
	105	//
	106	Double_t xd[3],pd[3],signd;
	107	Double_t xm[3],pm[3],signm;
	108	for (Int_t i=0;i<3;i++){
	109	xd[i] = fT2.fTPCinR1[i];
	110	pd[i] = fT2.fTPCinP1[i];
	111	xm[i] = fT1.fTPCinR1[i];
	112	pm[i] = fT1.fTPCinP1[i];
	113	}
	114	signd = fT2.fSign<0 ? -1:1;
	115	signm = fT1.fSign<0 ? -1:1;
	116
	117	AliHelix dhelix1(xd,pd,signd);
	118	dhelix1.GetMomentum(0,pp,0);
	119	dhelix1.Evaluate(0,xx);
	120	//
	121	// Double_t x2[3],p2[3];
	122	//
	123	AliHelix mhelix(xm,pm,signm);
	124	//
	125	//find intersection linear
	126	//
4a9220d4	127	Double_t phase[2][2] = { {0,0}, {0,0}};
4a9220d4	128	Double_t radius[2] = {0} ;
7cc34f08	129	Int_t points = dhelix1.GetRPHIintersections(mhelix, phase, radius,200);
	130	Double_t delta1=10000,delta2=10000;
	131
	132	if (points==1){
	133	fMinR = TMath::Sqrt(radius[0]);
	134	}
	135	if (points==2){
	136	fMinR =TMath::Min(TMath::Sqrt(radius[0]),TMath::Sqrt(radius[1]));
	137	}
	138
	139	if (points>0){
	140	dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
	141	dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
	142	dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
	143	}
	144	if (points==2){
	145	dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
	146	dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
	147	dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
	148	}
	149	if (points==1){
	150	fMinR = TMath::Sqrt(radius[0]);
	151	fDistMinR = delta1;
	152	}
	153	if (points==2){
	154	if (radius[0]<radius[1]){
	155	fMinR = TMath::Sqrt(radius[0]);
	156	fDistMinR = delta1;
	157	}
	158	else{
	159	fMinR = TMath::Sqrt(radius[1]);
	160	fDistMinR = delta2;
	161	}
	162	}
	163	//
	164	//
	165	distance1 = TMath::Min(delta1,delta2);
	166	//
	167	//find intersection parabolic
	168	//
	169	points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
	170	delta1=10000,delta2=10000;
	171
	172	if (points>0){
	173	dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
	174	}
	175	if (points==2){
	176	dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
	177	}
	178
	179	distance2 = TMath::Min(delta1,delta2);
	180	if (delta1<delta2){
	181	//get V0 info
	182	dhelix1.Evaluate(phase[0][0],fXr);
	183	dhelix1.GetMomentum(phase[0][0],fPdr);
	184	mhelix.GetMomentum(phase[0][1],fPm);
	185	dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],fAngle);
	186	fRr = TMath::Sqrt(radius[0]);
	187	}
	188	else{
	189	dhelix1.Evaluate(phase[1][0],fXr);
	190	dhelix1.GetMomentum(phase[1][0], fPdr);
	191	mhelix.GetMomentum(phase[1][1], fPm);
	192	dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],fAngle);
193	fRr = TMath::Sqrt(radius[1]);
194	}
195	fDist1 = TMath::Sqrt(distance1);
196	fDist2 = TMath::Sqrt(distance2);
197
198	if (fDist2<10.5){
199	Double_t x,alpha,param[5],cov[15];
200	//
201	fT1.GetESDtrack()->GetInnerExternalParameters(alpha,x,param);
202	fT1.GetESDtrack()->GetInnerExternalCovariance(cov);
203	AliExternalTrackParam paramm(x,alpha,param,cov);
204	//
205	fT2.GetESDtrack()->GetInnerExternalParameters(alpha,x,param);
206	fT2.GetESDtrack()->GetInnerExternalCovariance(cov);
207	AliExternalTrackParam paramd(x,alpha,param,cov);
208	/*
209	AliESDkink kink;
210	kink.Update(&paramm,&paramd);
211	// kink.Dump();
212	Double_t diff = kink.fRr-fRr;
213	Double_t diff2 = kink.fDist2-fDist2;
214	printf("Diff\t%f\t%f\n",diff,diff2);
215	*/
216	}
217
218	//
219	//
220	}
221
222	}
223