[u/mrichter/AliRoot.git] / HLT / QA / tasks / AliAnalysisTaskD0Trigger.cxx

// $Id$
//**************************************************************************
//* This file is property of and copyright by the ALICE HLT Project        *
//* ALICE Experiment at CERN, All rights reserved.                         *
//*                                                                        *
//* Primary Authors: Gaute Ovrebekk <st05886@alf.uib.no>,                *                  
//*                  for The ALICE HLT Project.                            *
//*                                                                        *
//* 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.                  *
//**************************************************************************


/** @file  AliAnalysisTaskD0Trigger.cxx  
    @author Gaute Ovrebekk
    @date 
    @brief An analysis task for the D0 Trigger.    
*/

class AliAnalysisTask;
class AliAnalysisManager;

#include "TH1F.h"
#include "TCanvas.h"
#include "AliESDEvent.h"
#include "AliESDInputHandler.h"
#include "TObjArray.h"
#include "TDatabasePDG.h"
#include "AliESDVertex.h"
#include "TMath.h"
#include "AliExternalTrackParam.h"
#include "AliAODVertex.h"
#include "TDatabasePDG.h"
#include "AliESDtrack.h"
#include "TVector3.h"
#include "AliVertexerTracks.h"
#include "AliKFVertex.h"
#include "TDatabasePDG.h"
#include "TVector3.h"

#include "AliAnalysisTaskD0Trigger.h"

ClassImp(AliAnalysisTaskD0Trigger)

AliAnalysisTaskD0Trigger::AliAnalysisTaskD0Trigger()
:
AliAnalysisTaskSE()
  , fOutputList(0)
  , fPtMin(0.0)
  , fdca(0.0)
  , finvMass(0.0)     
  , fcosThetaStar(0.0)
  , fd0(0.0) 
  , fd0d0(0.0)
  , fcosPoint(0.0)
  , fD0PDG(TDatabasePDG::Instance()->GetParticle(421)->Mass())
  , fD0massHLT(NULL)
  , fD0ptHLT(NULL)
  , fD0massOFF(NULL)
  , fD0ptOFF(NULL)
  , fPos()
  , fNeg()
  , ftwoTrackArray(NULL)
  , fTotalD0HLT(0)
  , fTotalD0OFF(0)
  , fField(0)
  , fNevents(0)
  , fuseKF(false)
{
  // Constructor
  // Define input and output slots here
  // Input slot #0 works with a TChain
  // DefineInput(0, TChain::Class());
  // Output slot #0 writes into a TH1 container

  //DefineOutput(1, TList::Class());
}

AliAnalysisTaskD0Trigger::AliAnalysisTaskD0Trigger(const char *name,float cuts[7])
:
AliAnalysisTaskSE(name)
  , fOutputList(0)
  , fPtMin(cuts[0])
  , fdca(cuts[1])
  , finvMass(cuts[2])     
  , fcosThetaStar(cuts[3])
  , fd0(cuts[4]) 
  , fd0d0(cuts[5])
  , fcosPoint(cuts[6])
  , fD0PDG(TDatabasePDG::Instance()->GetParticle(421)->Mass())
  , fD0massHLT(NULL)
  , fD0ptHLT(NULL)
  , fD0massOFF(NULL)
  , fD0ptOFF(NULL)
  , fPos()
  , fNeg()
  , ftwoTrackArray(NULL)
  , fTotalD0HLT(0)
  , fTotalD0OFF(0)
  , fField(0)
  , fNevents(0)
  , fuseKF(false)
{
  // Constructor
  // Define input and output slots here
  // Input slot #0 works with a TChain
  // DefineInput(0, TChain::Class());
  // Output slot #0 writes into a TH1 container

  DefineOutput(1, TList::Class());
}
void AliAnalysisTaskD0Trigger::UserCreateOutputObjects(){
  // Create histograms

  OpenFile(1);

  fOutputList = new TList();
  fOutputList->SetName(GetName());

  fD0massHLT = new TH1F("hMassHLT","D^{0} mass plot from HLT reconstruction",100,1.7,2);
  fD0ptHLT = new TH1F("hPtHLT","D^{0} Pt plot HLT reconstruction",20,0,20);

  fD0massOFF = new TH1F("hMassOFF","D^{0} mass plot Offline reconstruction",100,1.7,2);
  fD0ptOFF = new TH1F("hPtOFF","D^{0} Pt plot Offline reconstruction",20,0,20);

  fOutputList->Add(fD0massHLT);
  fOutputList->Add(fD0ptHLT);
  fOutputList->Add(fD0massOFF);
  fOutputList->Add(fD0ptOFF);
  
}

void AliAnalysisTaskD0Trigger::NotifyRun(){
  // see header file of AliAnalysisTask for documentation
}

void AliAnalysisTaskD0Trigger::UserExec(Option_t *){
  // see header file of AliAnalysisTask for documentation

  AliESDEvent *esdOFF = dynamic_cast<AliESDEvent*>(InputEvent());
  
  if(!esdOFF){
    Printf("ERROR: fESD not available");
    return;
  }
  
  AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*>(fInputHandler);
  AliESDEvent *esdHLT = NULL;   
  if(esdH) esdHLT = esdH->GetHLTEvent();
    
  if(!esdHLT){
    Printf("ERROR: HLTesd not available");
    return;
  }
  
  fNevents++;
  
  ftwoTrackArray = new TObjArray(2);
  Int_t nD0=0;

  fField = esdHLT->GetMagneticField();
  const AliESDVertex* pvHLT = esdHLT->GetPrimaryVertexTracks();
  AliESDVertex *pVertexHLT =  new AliESDVertex(*pvHLT);
  if(pVertexHLT->GetNContributors()<2){
    //Printf("ERROR: Contributiors to HLT Primary vertex is to low or not set");
    return;
  }

  for(Int_t it=0;it<esdHLT->GetNumberOfTracks();it++){
    SingleTrackSelect(esdHLT->GetTrack(it),pVertexHLT);
  }
  
  RecD0(nD0,pVertexHLT,true);
  fTotalD0HLT+=nD0;
  
  fPos.clear();
  fNeg.clear();

  nD0=0;

  fField = esdOFF->GetMagneticField();
  const AliESDVertex* pvOFF = esdOFF->GetPrimaryVertexTracks();
  AliESDVertex *pVertexOFF =  new AliESDVertex(*pvOFF);
  if(pVertexOFF->GetNContributors()<2){
    //Printf("ERROR: Contributiors to OFFLINE Primary vertex is to low or not set");
    return;
  }
  for(Int_t it=0;it<esdOFF->GetNumberOfTracks();it++){
    SingleTrackSelect(esdOFF->GetTrack(it),pVertexOFF);
  }
  
  RecD0(nD0,pVertexOFF,false);  
  fTotalD0OFF+=nD0;
  
  fPos.clear();
  fNeg.clear();

  // Post output data.
  PostData(1, fOutputList);
  ftwoTrackArray->Clear();
  delete pVertexHLT;
  delete pVertexOFF;
}

void AliAnalysisTaskD0Trigger::Terminate(Option_t *){
  Printf("Event Number: %d",fNevents);
  Printf("Total Number of D0 foundfor HLT: %d",fTotalD0HLT);
  Printf("Total Number of D0 found for OFFLINE: %d",fTotalD0OFF);
}

void AliAnalysisTaskD0Trigger::SingleTrackSelect(AliExternalTrackParam* t, AliESDVertex *pV){
  // Offline har || på disse kuttene på de to henfallsproduktene 
  Double_t pv[3];
  pV->GetXYZ(pv);
  
  if(t->Pt()<fPtMin){return;}
  if(TMath::Abs(t->GetD(pv[0],pv[1],fField)) > fd0){return;}
  
  if(t->Charge()>0){
    fPos.push_back(t);
  }
  else{
    fNeg.push_back(t);
  }
}

void AliAnalysisTaskD0Trigger::RecD0(Int_t& nD0, AliESDVertex *pV,bool isHLT){
  // Reconstructing D0
  Double_t starD0,starD0bar,xdummy,ydummy; 
  Double_t d0[2];
  Double_t svpos[3];
  Double_t pvpos[3];
  ftwoTrackArray->Clear();

  if(!pV){
    Printf("No Primary Vertex");
    return;
  }
  pV->GetXYZ(pvpos);
    
  for(UInt_t ip=0;ip<fPos.size();ip++){
    AliExternalTrackParam *tP=fPos[ip];
    for(UInt_t in=0;in<fNeg.size();in++){
      AliExternalTrackParam *tN=fNeg[in];
          
      tP->PropagateToDCA(pV,fField,kVeryBig);  //do I need this??????
      tN->PropagateToDCA(pV,fField,kVeryBig);
      
      Double_t dcatPtN = tP->GetDCA(tN,fField,xdummy,ydummy);
      if(dcatPtN>fdca) { continue; }
      
      ftwoTrackArray->AddAt(tP,0);
      ftwoTrackArray->AddAt(tN,1);
      AliAODVertex *vertexp1n1 = ReconstructSecondaryVertex(ftwoTrackArray,fField,pV,fuseKF);
      if(!vertexp1n1) { 
        ftwoTrackArray->Clear();
        continue; 
      }
      
      vertexp1n1->GetXYZ(svpos);
      
      tP->PropagateToDCA(vertexp1n1,fField,kVeryBig); 
      tN->PropagateToDCA(vertexp1n1,fField,kVeryBig);
      
      if((TMath::Abs(InvMass(tN,tP)-fD0PDG)) > finvMass && TMath::Abs((InvMass(tP,tN))-fD0PDG) > finvMass){continue;}
      CosThetaStar(tN,tP,starD0,starD0bar);
      if(TMath::Abs(starD0) > fcosThetaStar && TMath::Abs(starD0bar) > fcosThetaStar){continue;}
      d0[0] = tP->GetD(pvpos[0],pvpos[1],fField);
      d0[1] = tN->GetD(pvpos[0],pvpos[1],fField);
      if((d0[0]*d0[1]) > fd0d0){continue;}
      if(PointingAngle(tN,tP,pvpos,svpos) < fcosPoint){continue;}
      
      if(isHLT){
	fD0massHLT->Fill(InvMass(tN,tP));
	fD0massHLT->Fill(InvMass(tP,tN));
	fD0ptHLT->Fill(Pt(tP,tN));
      }
      else{
	fD0massOFF->Fill(InvMass(tN,tP));
	fD0massOFF->Fill(InvMass(tP,tN));
	fD0ptOFF->Fill(Pt(tP,tN));
      }
      nD0++;
      delete vertexp1n1;
    }
  }
}

Double_t AliAnalysisTaskD0Trigger::InvMass(AliExternalTrackParam* d1, AliExternalTrackParam* d2)
{
  // Calculating the invariant mass
  Double_t mpi=TDatabasePDG::Instance()->GetParticle(211)->Mass();
  Double_t mK=TDatabasePDG::Instance()->GetParticle(321)->Mass();

  Double_t energy[2]; 
  energy[1] = TMath::Sqrt(mK*mK+d1->GetP()*d1->GetP());
  energy[0] = TMath::Sqrt(mpi*mpi+d2->GetP()*d2->GetP());

  Double_t p1[3],p2[3];
  d1->GetPxPyPz(p1);
  d2->GetPxPyPz(p2);
  
  Double_t momTot2 = (p1[0]+p2[0])*(p1[0]+p2[0])+
                     (p1[1]+p2[1])*(p1[1]+p2[1])+
                     (p1[2]+p2[2])*(p1[2]+p2[2]);

  return TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-momTot2);

}

void AliAnalysisTaskD0Trigger::CosThetaStar(AliExternalTrackParam* d1, AliExternalTrackParam* d2,Double_t &D0,Double_t &D0bar)
{
  //Calculating the decay angle
  Double_t mD0 = TDatabasePDG::Instance()->GetParticle(421)->Mass();
  Double_t mpi=TDatabasePDG::Instance()->GetParticle(211)->Mass();
  Double_t mK=TDatabasePDG::Instance()->GetParticle(321)->Mass();

  Double_t pStar = TMath::Sqrt(TMath::Power(mD0*mD0-mK*mK-mpi*mpi,2.)-4.*mK*mK*mpi*mpi)/(2.*mD0);
 
  Double_t px = d1->Px()+d2->Px();
  Double_t py = d1->Py()+d2->Py();
  Double_t pz = d1->Pz()+d2->Pz();
  Double_t p = TMath::Sqrt(px*px+py*py+pz*pz);
  Double_t energy = TMath::Sqrt(p*p+mD0*mD0);

  Double_t beta = p/energy;
  Double_t gamma = energy/mD0;
  
  Double_t qL;
  TVector3 mom(d1->Px(),d1->Py(),d1->Pz());
  TVector3 momD(px,py,pz);
  qL = mom.Dot(momD)/momD.Mag();

  D0 = (qL/gamma-beta*TMath::Sqrt(pStar*pStar+mK*mK))/pStar;
  
  TVector3 mom2(d2->Px(),d2->Py(),d2->Pz());
  TVector3 momD2(px,py,pz);
  qL = mom2.Dot(momD2)/momD2.Mag();

  D0bar = (qL/gamma-beta*TMath::Sqrt(pStar*pStar+mK*mK))/pStar;

}

Double_t AliAnalysisTaskD0Trigger::PointingAngle(AliExternalTrackParam* n, AliExternalTrackParam* p, Double_t *pv, Double_t *sv)
{
  // Calcutating the pointing angle
  TVector3 mom(n->Px()+p->Px(),n->Py()+p->Py(),n->Pz()+p->Pz());
  TVector3 flight(sv[0]-pv[0],sv[1]-pv[1],sv[2]-pv[2]);
  
  double pta = mom.Angle(flight);

  return TMath::Cos(pta); 
}

AliAODVertex* AliAnalysisTaskD0Trigger::ReconstructSecondaryVertex(TObjArray *trkArray, Double_t b, const AliESDVertex *v, bool useKF)
{
  // Finding the vertex of the two tracks in trkArray
  AliESDVertex *vertexESD = 0;
  AliAODVertex *vertexAOD = 0;
  
  if(!useKF){
    AliVertexerTracks *vertexer = new AliVertexerTracks(b);
    AliESDVertex* vertex =  const_cast<AliESDVertex*>(v);
    vertexer->SetVtxStart(vertex);
    //if(isESD){vertexESD = (AliESDVertex*)vertexer->VertexForSelectedESDTracks(trkArray);}
    UShort_t *id = new UShort_t[2];
    AliExternalTrackParam *t1 = (AliExternalTrackParam*) trkArray->At(0);
    AliExternalTrackParam *t2 = (AliExternalTrackParam*) trkArray->At(1);
    id[0]=(UShort_t) t1->GetID();
    id[1]=(UShort_t) t2->GetID();
    vertexESD = (AliESDVertex*)vertexer->VertexForSelectedTracks(trkArray,id);
    delete [] id;
    delete vertexer; vertexer=NULL;
    
    if(!vertexESD) return vertexAOD;
    
    if(vertexESD->GetNContributors()!=trkArray->GetEntriesFast()) { 
      //AliDebug(2,"vertexing failed"); 
      delete vertexESD; vertexESD=NULL;
      delete vertex;
      return vertexAOD;
    }
  }
  else{
    AliKFParticle::SetField(b);
    
    AliKFVertex vertexKF;
    
    Int_t nTrks = trkArray->GetEntriesFast();
    for(Int_t i=0; i<nTrks; i++) {
      AliESDtrack *esdTrack = (AliESDtrack*)trkArray->At(i);
      AliKFParticle daughterKF(*esdTrack,211);
      vertexKF.AddDaughter(daughterKF);
    }
    vertexESD = new AliESDVertex(vertexKF.Parameters(),
				 vertexKF.CovarianceMatrix(),
				 vertexKF.GetChi2(),
				 vertexKF.GetNContributors());
  }
  // convert to AliAODVertex
  Double_t pos[3],cov[6],chi2perNDF;
  vertexESD->GetXYZ(pos); // position
  vertexESD->GetCovMatrix(cov); //covariance matrix
  chi2perNDF = vertexESD->GetChi2toNDF();
  //dispersion = vertexESD->GetDispersion();
  delete vertexESD; vertexESD=NULL;

  Int_t nprongs= trkArray->GetEntriesFast();
  vertexAOD = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,nprongs);

  return vertexAOD;

}

Double_t AliAnalysisTaskD0Trigger::Pt(AliExternalTrackParam* d1, AliExternalTrackParam* d2)
{
  //Calculating pT of the two tracks
  Double_t p1[3],p2[3];
  d1->GetPxPyPz(p1);
  d2->GetPxPyPz(p2);
  
  Double_t pt2 = (p1[0]+p2[0])*(p1[0]+p2[0]) + (p1[1]+p2[1])*(p1[1]+p2[1]);

  return TMath::Sqrt(pt2);
}
Commit	Line	Data
cd785861	1	// $Id$
	2	//**************************************************************************
	3	//* This file is property of and copyright by the ALICE HLT Project *
	4	//* ALICE Experiment at CERN, All rights reserved. *
	5	//* *
	6	//* Primary Authors: Gaute Ovrebekk <st05886@alf.uib.no>, *
	7	//* for The ALICE HLT Project. *
	8	//* *
	9	//* Permission to use, copy, modify and distribute this software and its *
	10	//* documentation strictly for non-commercial purposes is hereby granted *
	11	//* without fee, provided that the above copyright notice appears in all *
	12	//* copies and that both the copyright notice and this permission notice *
	13	//* appear in the supporting documentation. The authors make no claims *
	14	//* about the suitability of this software for any purpose. It is *
	15	//* provided "as is" without express or implied warranty. *
	16	//**************************************************************************
	17
	18
	19	/** @file AliAnalysisTaskD0Trigger.cxx
	20	@author Gaute Ovrebekk
	21	@date
	22	@brief An analysis task for the D0 Trigger.
	23	*/
	24
	25	class AliAnalysisTask;
	26	class AliAnalysisManager;
	27
	28	#include "TH1F.h"
	29	#include "TCanvas.h"
	30	#include "AliESDEvent.h"
	31	#include "AliESDInputHandler.h"
	32	#include "TObjArray.h"
	33	#include "TDatabasePDG.h"
	34	#include "AliESDVertex.h"
	35	#include "TMath.h"
	36	#include "AliExternalTrackParam.h"
	37	#include "AliAODVertex.h"
	38	#include "TDatabasePDG.h"
	39	#include "AliESDtrack.h"
	40	#include "TVector3.h"
	41	#include "AliVertexerTracks.h"
	42	#include "AliKFVertex.h"
	43	#include "TDatabasePDG.h"
	44	#include "TVector3.h"
	45
	46	#include "AliAnalysisTaskD0Trigger.h"
	47
	48	ClassImp(AliAnalysisTaskD0Trigger)
	49
	50	AliAnalysisTaskD0Trigger::AliAnalysisTaskD0Trigger()
	51	:
	52	AliAnalysisTaskSE()
	53	, fOutputList(0)
	54	, fPtMin(0.0)
	55	, fdca(0.0)
	56	, finvMass(0.0)
	57	, fcosThetaStar(0.0)
	58	, fd0(0.0)
	59	, fd0d0(0.0)
	60	, fcosPoint(0.0)
80ec6f2e	61	, fD0PDG(TDatabasePDG::Instance()->GetParticle(421)->Mass())
cd785861	62	, fD0massHLT(NULL)
	63	, fD0ptHLT(NULL)
	64	, fD0massOFF(NULL)
	65	, fD0ptOFF(NULL)
	66	, fPos()
	67	, fNeg()
	68	, ftwoTrackArray(NULL)
	69	, fTotalD0HLT(0)
	70	, fTotalD0OFF(0)
	71	, fField(0)
	72	, fNevents(0)
	73	, fuseKF(false)
	74	{
	75	// Constructor
	76	// Define input and output slots here
	77	// Input slot #0 works with a TChain
	78	// DefineInput(0, TChain::Class());
	79	// Output slot #0 writes into a TH1 container
	80
	81	//DefineOutput(1, TList::Class());
	82	}
	83
	84	AliAnalysisTaskD0Trigger::AliAnalysisTaskD0Trigger(const char *name,float cuts[7])
	85	:
	86	AliAnalysisTaskSE(name)
	87	, fOutputList(0)
	88	, fPtMin(cuts[0])
	89	, fdca(cuts[1])
	90	, finvMass(cuts[2])
	91	, fcosThetaStar(cuts[3])
	92	, fd0(cuts[4])
	93	, fd0d0(cuts[5])
	94	, fcosPoint(cuts[6])
80ec6f2e	95	, fD0PDG(TDatabasePDG::Instance()->GetParticle(421)->Mass())
cd785861	96	, fD0massHLT(NULL)
	97	, fD0ptHLT(NULL)
	98	, fD0massOFF(NULL)
	99	, fD0ptOFF(NULL)
	100	, fPos()
	101	, fNeg()
	102	, ftwoTrackArray(NULL)
	103	, fTotalD0HLT(0)
	104	, fTotalD0OFF(0)
	105	, fField(0)
	106	, fNevents(0)
	107	, fuseKF(false)
	108	{
	109	// Constructor
	110	// Define input and output slots here
	111	// Input slot #0 works with a TChain
	112	// DefineInput(0, TChain::Class());
	113	// Output slot #0 writes into a TH1 container
	114
	115	DefineOutput(1, TList::Class());
	116	}
	117	void AliAnalysisTaskD0Trigger::UserCreateOutputObjects(){
	118	// Create histograms
	119
	120	OpenFile(1);
	121
	122	fOutputList = new TList();
	123	fOutputList->SetName(GetName());
	124
	125	fD0massHLT = new TH1F("hMassHLT","D^{0} mass plot from HLT reconstruction",100,1.7,2);
	126	fD0ptHLT = new TH1F("hPtHLT","D^{0} Pt plot HLT reconstruction",20,0,20);
	127
	128	fD0massOFF = new TH1F("hMassOFF","D^{0} mass plot Offline reconstruction",100,1.7,2);
	129	fD0ptOFF = new TH1F("hPtOFF","D^{0} Pt plot Offline reconstruction",20,0,20);
	130
	131	fOutputList->Add(fD0massHLT);
	132	fOutputList->Add(fD0ptHLT);
	133	fOutputList->Add(fD0massOFF);
	134	fOutputList->Add(fD0ptOFF);
	135
	136	}
	137
	138	void AliAnalysisTaskD0Trigger::NotifyRun(){
	139	// see header file of AliAnalysisTask for documentation
	140	}
	141
	142	void AliAnalysisTaskD0Trigger::UserExec(Option_t *){
	143	// see header file of AliAnalysisTask for documentation
	144
	145	AliESDEvent esdOFF = dynamic_cast<AliESDEvent>(InputEvent());
	146
	147	if(!esdOFF){
	148	Printf("ERROR: fESD not available");
	149	return;
	150	}
	151
	152	AliESDInputHandler esdH = dynamic_cast<AliESDInputHandler>(fInputHandler);
	153	AliESDEvent *esdHLT = NULL;
	154	if(esdH) esdHLT = esdH->GetHLTEvent();
	155
	156	if(!esdHLT){
	157	Printf("ERROR: HLTesd not available");
	158	return;
	159	}
160
161	fNevents++;
162
163	ftwoTrackArray = new TObjArray(2);
164	Int_t nD0=0;
165
166	fField = esdHLT->GetMagneticField();
167	const AliESDVertex* pvHLT = esdHLT->GetPrimaryVertexTracks();
168	AliESDVertex pVertexHLT = new AliESDVertex(pvHLT);
169	if(pVertexHLT->GetNContributors()<2){
170	//Printf("ERROR: Contributiors to HLT Primary vertex is to low or not set");
171	return;
172	}
173
174	for(Int_t it=0;it<esdHLT->GetNumberOfTracks();it++){
175	SingleTrackSelect(esdHLT->GetTrack(it),pVertexHLT);
176	}
177
178	RecD0(nD0,pVertexHLT,true);
179	fTotalD0HLT+=nD0;
180
181	fPos.clear();
182	fNeg.clear();
183
184	nD0=0;
185
186	fField = esdOFF->GetMagneticField();
187	const AliESDVertex* pvOFF = esdOFF->GetPrimaryVertexTracks();
188	AliESDVertex pVertexOFF = new AliESDVertex(pvOFF);
189	if(pVertexOFF->GetNContributors()<2){
190	//Printf("ERROR: Contributiors to OFFLINE Primary vertex is to low or not set");
191	return;
192	}
193	for(Int_t it=0;it<esdOFF->GetNumberOfTracks();it++){
194	SingleTrackSelect(esdOFF->GetTrack(it),pVertexOFF);
195	}
196
197	RecD0(nD0,pVertexOFF,false);
198	fTotalD0OFF+=nD0;
199
200	fPos.clear();
201	fNeg.clear();
202
203	// Post output data.
204	PostData(1, fOutputList);
205	ftwoTrackArray->Clear();
206	delete pVertexHLT;
207	delete pVertexOFF;
208	}
209
210	void AliAnalysisTaskD0Trigger::Terminate(Option_t *){
211	Printf("Event Number: %d",fNevents);
212	Printf("Total Number of D0 foundfor HLT: %d",fTotalD0HLT);
213	Printf("Total Number of D0 found for OFFLINE: %d",fTotalD0OFF);
214	}
215
216	void AliAnalysisTaskD0Trigger::SingleTrackSelect(AliExternalTrackParam* t, AliESDVertex *pV){
217	// Offline har \|\| på disse kuttene på de to henfallsproduktene
218	Double_t pv[3];
219	pV->GetXYZ(pv);
220
221	if(t->Pt()<fPtMin){return;}
222	if(TMath::Abs(t->GetD(pv[0],pv[1],fField)) > fd0){return;}
223
224	if(t->Charge()>0){
225	fPos.push_back(t);
226	}
227	else{
228	fNeg.push_back(t);
229	}
230	}
231
232	void AliAnalysisTaskD0Trigger::RecD0(Int_t& nD0, AliESDVertex *pV,bool isHLT){
80ec6f2e	233	// Reconstructing D0
80ec6f2e	234	Double_t starD0,starD0bar,xdummy,ydummy;
cd785861	235	Double_t d0[2];
	236	Double_t svpos[3];
	237	Double_t pvpos[3];
	238	ftwoTrackArray->Clear();
	239
	240	if(!pV){
	241	Printf("No Primary Vertex");
	242	return;
	243	}
	244	pV->GetXYZ(pvpos);
	245
	246	for(UInt_t ip=0;ip<fPos.size();ip++){
	247	AliExternalTrackParam *tP=fPos[ip];
	248	for(UInt_t in=0;in<fNeg.size();in++){
	249	AliExternalTrackParam *tN=fNeg[in];
	250
	251	tP->PropagateToDCA(pV,fField,kVeryBig); //do I need this??????
	252	tN->PropagateToDCA(pV,fField,kVeryBig);
	253
	254	Double_t dcatPtN = tP->GetDCA(tN,fField,xdummy,ydummy);
	255	if(dcatPtN>fdca) { continue; }
	256
	257	ftwoTrackArray->AddAt(tP,0);
	258	ftwoTrackArray->AddAt(tN,1);
	259	AliAODVertex *vertexp1n1 = ReconstructSecondaryVertex(ftwoTrackArray,fField,pV,fuseKF);
	260	if(!vertexp1n1) {
	261	ftwoTrackArray->Clear();
	262	continue;
	263	}
	264
	265	vertexp1n1->GetXYZ(svpos);
	266
	267	tP->PropagateToDCA(vertexp1n1,fField,kVeryBig);
	268	tN->PropagateToDCA(vertexp1n1,fField,kVeryBig);
	269
80ec6f2e	270	if((TMath::Abs(InvMass(tN,tP)-fD0PDG)) > finvMass && TMath::Abs((InvMass(tP,tN))-fD0PDG) > finvMass){continue;}
	271	CosThetaStar(tN,tP,starD0,starD0bar);
	272	if(TMath::Abs(starD0) > fcosThetaStar && TMath::Abs(starD0bar) > fcosThetaStar){continue;}
cd785861	273	d0[0] = tP->GetD(pvpos[0],pvpos[1],fField);
	274	d0[1] = tN->GetD(pvpos[0],pvpos[1],fField);
	275	if((d0[0]*d0[1]) > fd0d0){continue;}
80ec6f2e	276	if(PointingAngle(tN,tP,pvpos,svpos) < fcosPoint){continue;}
cd785861	277
	278	if(isHLT){
	279	fD0massHLT->Fill(InvMass(tN,tP));
	280	fD0massHLT->Fill(InvMass(tP,tN));
	281	fD0ptHLT->Fill(Pt(tP,tN));
	282	}
	283	else{
	284	fD0massOFF->Fill(InvMass(tN,tP));
	285	fD0massOFF->Fill(InvMass(tP,tN));
	286	fD0ptOFF->Fill(Pt(tP,tN));
	287	}
	288	nD0++;
	289	delete vertexp1n1;
	290	}
	291	}
	292	}
	293
	294	Double_t AliAnalysisTaskD0Trigger::InvMass(AliExternalTrackParam* d1, AliExternalTrackParam* d2)
	295	{
80ec6f2e	296	// Calculating the invariant mass
cd785861	297	Double_t mpi=TDatabasePDG::Instance()->GetParticle(211)->Mass();
	298	Double_t mK=TDatabasePDG::Instance()->GetParticle(321)->Mass();
	299
	300	Double_t energy[2];
	301	energy[1] = TMath::Sqrt(mKmK+d1->GetP()d1->GetP());
	302	energy[0] = TMath::Sqrt(mpimpi+d2->GetP()d2->GetP());
	303
	304	Double_t p1[3],p2[3];
	305	d1->GetPxPyPz(p1);
	306	d2->GetPxPyPz(p2);
	307
	308	Double_t momTot2 = (p1[0]+p2[0])*(p1[0]+p2[0])+
	309	(p1[1]+p2[1])*(p1[1]+p2[1])+
	310	(p1[2]+p2[2])*(p1[2]+p2[2]);
	311
	312	return TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-momTot2);
	313
	314	}
	315
80ec6f2e	316	void AliAnalysisTaskD0Trigger::CosThetaStar(AliExternalTrackParam* d1, AliExternalTrackParam* d2,Double_t &D0,Double_t &D0bar)
cd785861	317	{
80ec6f2e	318	//Calculating the decay angle
cd785861	319	Double_t mD0 = TDatabasePDG::Instance()->GetParticle(421)->Mass();
	320	Double_t mpi=TDatabasePDG::Instance()->GetParticle(211)->Mass();
	321	Double_t mK=TDatabasePDG::Instance()->GetParticle(321)->Mass();
	322
	323	Double_t pStar = TMath::Sqrt(TMath::Power(mD0mD0-mKmK-mpimpi,2.)-4.mKmKmpimpi)/(2.mD0);
	324
	325	Double_t px = d1->Px()+d2->Px();
	326	Double_t py = d1->Py()+d2->Py();
	327	Double_t pz = d1->Pz()+d2->Pz();
	328	Double_t p = TMath::Sqrt(pxpx+pypy+pz*pz);
	329	Double_t energy = TMath::Sqrt(pp+mD0mD0);
	330
	331	Double_t beta = p/energy;
	332	Double_t gamma = energy/mD0;
	333
	334	Double_t qL;
	335	TVector3 mom(d1->Px(),d1->Py(),d1->Pz());
	336	TVector3 momD(px,py,pz);
	337	qL = mom.Dot(momD)/momD.Mag();
	338
	339	D0 = (qL/gamma-betaTMath::Sqrt(pStarpStar+mK*mK))/pStar;
	340
	341	TVector3 mom2(d2->Px(),d2->Py(),d2->Pz());
	342	TVector3 momD2(px,py,pz);
	343	qL = mom2.Dot(momD2)/momD2.Mag();
	344
	345	D0bar = (qL/gamma-betaTMath::Sqrt(pStarpStar+mK*mK))/pStar;
	346
	347	}
	348
80ec6f2e	349	Double_t AliAnalysisTaskD0Trigger::PointingAngle(AliExternalTrackParam* n, AliExternalTrackParam* p, Double_t pv, Double_t sv)
cd785861	350	{
80ec6f2e	351	// Calcutating the pointing angle
cd785861	352	TVector3 mom(n->Px()+p->Px(),n->Py()+p->Py(),n->Pz()+p->Pz());
	353	TVector3 flight(sv[0]-pv[0],sv[1]-pv[1],sv[2]-pv[2]);
	354
	355	double pta = mom.Angle(flight);
	356
	357	return TMath::Cos(pta);
	358	}
	359
	360	AliAODVertex* AliAnalysisTaskD0Trigger::ReconstructSecondaryVertex(TObjArray trkArray, Double_t b, const AliESDVertex v, bool useKF)
	361	{
80ec6f2e	362	// Finding the vertex of the two tracks in trkArray
cd785861	363	AliESDVertex *vertexESD = 0;
	364	AliAODVertex *vertexAOD = 0;
	365
	366	if(!useKF){
	367	AliVertexerTracks *vertexer = new AliVertexerTracks(b);
80ec6f2e	368	AliESDVertex* vertex = const_cast<AliESDVertex*>(v);
80ec6f2e	369	vertexer->SetVtxStart(vertex);
cd785861	370	//if(isESD){vertexESD = (AliESDVertex*)vertexer->VertexForSelectedESDTracks(trkArray);}
	371	UShort_t *id = new UShort_t[2];
	372	AliExternalTrackParam t1 = (AliExternalTrackParam) trkArray->At(0);
	373	AliExternalTrackParam t2 = (AliExternalTrackParam) trkArray->At(1);
	374	id[0]=(UShort_t) t1->GetID();
	375	id[1]=(UShort_t) t2->GetID();
	376	vertexESD = (AliESDVertex*)vertexer->VertexForSelectedTracks(trkArray,id);
	377	delete [] id;
	378	delete vertexer; vertexer=NULL;
	379
	380	if(!vertexESD) return vertexAOD;
	381
	382	if(vertexESD->GetNContributors()!=trkArray->GetEntriesFast()) {
	383	//AliDebug(2,"vertexing failed");
	384	delete vertexESD; vertexESD=NULL;
80ec6f2e	385	delete vertex;
cd785861	386	return vertexAOD;
	387	}
	388	}
	389	else{
	390	AliKFParticle::SetField(b);
	391
	392	AliKFVertex vertexKF;
	393
	394	Int_t nTrks = trkArray->GetEntriesFast();
	395	for(Int_t i=0; i<nTrks; i++) {
	396	AliESDtrack esdTrack = (AliESDtrack)trkArray->At(i);
	397	AliKFParticle daughterKF(*esdTrack,211);
	398	vertexKF.AddDaughter(daughterKF);
	399	}
	400	vertexESD = new AliESDVertex(vertexKF.Parameters(),
	401	vertexKF.CovarianceMatrix(),
	402	vertexKF.GetChi2(),
	403	vertexKF.GetNContributors());
	404	}
	405	// convert to AliAODVertex
	406	Double_t pos[3],cov[6],chi2perNDF;
	407	vertexESD->GetXYZ(pos); // position
	408	vertexESD->GetCovMatrix(cov); //covariance matrix
	409	chi2perNDF = vertexESD->GetChi2toNDF();
	410	//dispersion = vertexESD->GetDispersion();
	411	delete vertexESD; vertexESD=NULL;
	412
	413	Int_t nprongs= trkArray->GetEntriesFast();
	414	vertexAOD = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,nprongs);
	415
	416	return vertexAOD;
	417
	418	}
	419
	420	Double_t AliAnalysisTaskD0Trigger::Pt(AliExternalTrackParam* d1, AliExternalTrackParam* d2)
	421	{
80ec6f2e	422	//Calculating pT of the two tracks
cd785861	423	Double_t p1[3],p2[3];
	424	d1->GetPxPyPz(p1);
	425	d2->GetPxPyPz(p2);
	426
	427	Double_t pt2 = (p1[0]+p2[0])(p1[0]+p2[0]) + (p1[1]+p2[1])(p1[1]+p2[1]);
	428
	429	return TMath::Sqrt(pt2);
	430	}