1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 // Get improved dca info. by ITS upgrade implemented by the ALICE Heavy Flavour Electron Group
19 // MinJung Kweon <minjung@physi.uni-heidelberg.de>
31 #include "AliAODTrack.h"
32 #include "AliAODPid.h"
33 #include "AliESDEvent.h"
34 #include "AliESDVertex.h"
35 #include "AliESDtrack.h"
37 #include "AliMCParticle.h"
38 #include "AliVEvent.h"
39 #include "AliVTrack.h"
40 #include "AliVParticle.h"
41 #include "AliVertexerTracks.h"
42 #include "AliVVertex.h"
43 #include "AliExternalTrackParam.h"
44 #include "AliMCEvent.h"
46 #include "AliHFEsmearDCA.h"
48 ClassImp(AliHFEsmearDCA)
51 //______________________________________________________
52 AliHFEsmearDCA::AliHFEsmearDCA(const Char_t */*name*/, const char *resfileCurURI, const char *resfileUpgURI, const Char_t */*title*/):
81 // Constructor to be used to create the task.
82 // The the URIs specify the resolution files to be used.
83 // They are expected to contain TGraphs with the resolutions
84 // for the current and the upgraded ITS (see code for details).
85 // One may also specify for how many tracks debug information
86 // is written to the output.
88 TFile *resfileCur=TFile::Open(resfileCurURI);
89 fD0RPResPCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResP" )));
90 fD0RPResKCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResK" )));
91 fD0RPResPiCur=new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResPi")));
92 fD0RPResECur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResE")));
93 fD0ZResPCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResP" )));
94 fD0ZResKCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResK" )));
95 fD0ZResPiCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResPi" )));
96 fD0ZResECur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResE")));
97 fPt1ResPCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResP" )));
98 fPt1ResKCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResK" )));
99 fPt1ResPiCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResPi" )));
100 fPt1ResECur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResE" )));
101 fD0RPResPCur ->SetName("D0RPResPCur" );
102 fD0RPResKCur ->SetName("D0RPResKCur" );
103 fD0RPResPiCur->SetName("D0RPResPiCur");
104 fD0RPResECur ->SetName("D0RPResECur");
105 fD0ZResPCur ->SetName("D0ZResPCur" );
106 fD0ZResKCur ->SetName("D0ZResKCur" );
107 fD0ZResPiCur ->SetName("D0ZResPiCur" );
108 fD0ZResECur ->SetName("D0ZResECur" );
109 fPt1ResPCur ->SetName("Pt1ResPCur" );
110 fPt1ResKCur ->SetName("Pt1ResKCur" );
111 fPt1ResPiCur ->SetName("Pt1ResPiCur" );
112 fPt1ResECur ->SetName("Pt1ResECur" );
114 TFile *resfileUpg=TFile::Open(resfileUpgURI );
115 fD0RPResPUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResP" )));
116 fD0RPResKUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResK" )));
117 fD0RPResPiUpg=new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResPi")));
118 fD0RPResEUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResE")));
119 fD0ZResPUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResP" )));
120 fD0ZResKUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResK" )));
121 fD0ZResPiUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResPi" )));
122 fD0ZResEUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResE" )));
123 fPt1ResPUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResP" )));
124 fPt1ResKUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResK" )));
125 fPt1ResPiUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResPi" )));
126 fPt1ResEUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResE" )));
127 fD0RPResPUpg ->SetName("D0RPResPUpg" );
128 fD0RPResKUpg ->SetName("D0RPResKUpg" );
129 fD0RPResPiUpg->SetName("D0RPResPiUpg");
130 fD0RPResEUpg ->SetName("D0RPResEUpg");
131 fD0ZResPUpg ->SetName("D0ZResPUpg" );
132 fD0ZResKUpg ->SetName("D0ZResKUpg" );
133 fD0ZResPiUpg ->SetName("D0ZResPiUpg" );
134 fD0ZResEUpg ->SetName("D0ZResEUpg" );
135 fPt1ResPUpg ->SetName("Pt1ResPUpg" );
136 fPt1ResKUpg ->SetName("Pt1ResKUpg" );
137 fPt1ResPiUpg ->SetName("Pt1ResPiUpg" );
138 fPt1ResEUpg ->SetName("Pt1ResEUpg" );
143 //______________________________________________________
144 AliHFEsmearDCA::AliHFEsmearDCA(const AliHFEsmearDCA &c):
147 fMCEvent(c.fMCEvent),
148 fD0ZResPCur (c.fD0ZResPCur),
149 fD0ZResKCur (c.fD0ZResKCur),
150 fD0ZResPiCur (c.fD0ZResPiCur),
151 fD0ZResECur (c.fD0ZResECur),
152 fD0RPResPCur (c.fD0RPResPCur),
153 fD0RPResKCur (c.fD0RPResKCur),
154 fD0RPResPiCur(c.fD0RPResPiCur),
155 fD0RPResECur (c.fD0RPResECur),
156 fPt1ResPCur (c.fPt1ResPCur),
157 fPt1ResKCur (c.fPt1ResKCur),
158 fPt1ResPiCur (c.fPt1ResPiCur),
159 fPt1ResECur (c.fPt1ResECur),
160 fD0ZResPUpg (c.fD0ZResPUpg),
161 fD0ZResKUpg (c.fD0ZResKUpg),
162 fD0ZResPiUpg (c.fD0ZResPiUpg),
163 fD0ZResEUpg (c.fD0ZResEUpg),
164 fD0RPResPUpg (c.fD0RPResPUpg),
165 fD0RPResKUpg (c.fD0RPResKUpg),
166 fD0RPResPiUpg(c.fD0RPResPiUpg),
167 fD0RPResEUpg (c.fD0RPResEUpg),
168 fPt1ResPUpg (c.fPt1ResPUpg),
169 fPt1ResKUpg (c.fPt1ResKUpg),
170 fPt1ResPiUpg (c.fPt1ResPiUpg),
171 fPt1ResEUpg (c.fPt1ResEUpg)
175 // Performs a deep copy
179 //______________________________________________________
180 AliHFEsmearDCA::~AliHFEsmearDCA(){
186 //______________________________________________________
187 void AliHFEsmearDCA::SetRecEventInfo(const TObject *event){
189 // Set Virtual event an make a copy
192 AliError("Pointer to AliVEvent !");
195 TString className(event->ClassName());
196 if (! (className.CompareTo("AliESDEvent")==0 || className.CompareTo("AliAODEvent")==0)) {
197 AliError("argument must point to an AliESDEvent or AliAODEvent !");
200 fEvent = (AliVEvent*) event;
204 //______________________________________________________
205 void AliHFEsmearDCA::GetImproveITSImpactParameters(AliVTrack *track, Double_t &dcaxyn, Double_t &dcaxyo, Double_t &dcaxySign, Double_t &dcaxySigo, Double_t &dcazn, Double_t &dcazo, Double_t &dcazSign, Double_t &dcazSigo){
207 // Get HFE impact parameter (with recalculated primary vertex)
210 AliDebug(1, "No Input event available\n");
213 const Double_t kBeampiperadius=3.;
214 const Double_t kBeampiperadiusnew=1.9;
215 TString type = track->IsA()->GetName();
216 Double_t dcan[2]={-999.,-999.};
217 Double_t covn[3]={-999.,-999.,-999.};
218 Double_t dcao[2]={-999.,-999.};
219 Double_t covo[3]={-999.,-999.,-999.};
221 // Get reconstructed track parameters
223 Double_t xyz[3],pxpypz[3],cv[21];
225 pxpypz[0]=track->Px();
226 pxpypz[1]=track->Py();
227 pxpypz[2]=track->Pz();
228 track->GetCovarianceXYZPxPyPz(cv);
229 Short_t sign = (Short_t)track->Charge();
231 //AliExternalTrackParam et(track); //it doesn't work, i don't know why yet! so, I use the line below
232 AliExternalTrackParam et(xyz,pxpypz,cv,sign);
233 Double_t *param=const_cast<Double_t*>(et.GetParameter());
235 const AliMCParticle *mc=static_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(track->GetLabel())));
238 Double_t mccv[36]={0.};
243 AliExternalTrackParam mct(mcx,mcp,mccv,mcc);
244 const Double_t *parammc=mct.GetParameter();
245 const AliVVertex *tmpvtx = fMCEvent->GetPrimaryVertex();
247 AliVertex vtx(mcx,1.,1);
248 AliVertex *tmpvtx2 = (AliVertex *)fMCEvent->GetPrimaryVertex();
250 mct.PropagateToDCA(tmpvtx2,track->GetBz(),10.); //mjtest
251 // Correct reference points and frames according to MC
252 // TODO: B-Field correct?
253 // TODO: failing propagation....
254 et.PropagateToDCA(tmpvtx2,track->GetBz(),10.);
255 et.Rotate(mct.GetAlpha());
258 // Select appropriate smearing functions
259 Double_t ptmc=TMath::Abs(mc->Pt());
267 switch (mc->Particle()->GetPdgCode()) {
268 case 2212: case -2212:
269 sd0rpo=EvalGraph(fD0RPResPCur,ptmc);
270 sd0zo =EvalGraph(fD0ZResPCur ,ptmc);
271 spt1o =EvalGraph(fPt1ResPCur ,ptmc);
272 sd0rpn=EvalGraph(fD0RPResPUpg,ptmc);
273 sd0zn =EvalGraph(fD0ZResPUpg ,ptmc);
274 spt1n =EvalGraph(fPt1ResPUpg ,ptmc);
277 sd0rpo=EvalGraph(fD0RPResKCur,ptmc);
278 sd0zo =EvalGraph(fD0ZResKCur ,ptmc);
279 spt1o =EvalGraph(fPt1ResKCur ,ptmc);
280 sd0rpn=EvalGraph(fD0RPResKUpg,ptmc);
281 sd0zn =EvalGraph(fD0ZResKUpg ,ptmc);
282 spt1n =EvalGraph(fPt1ResKUpg ,ptmc);
285 sd0rpo=EvalGraph(fD0RPResPiCur,ptmc);
286 sd0zo =EvalGraph(fD0ZResPiCur ,ptmc);
287 spt1o =EvalGraph(fPt1ResPiCur ,ptmc);
288 sd0rpn=EvalGraph(fD0RPResPiUpg,ptmc);
289 sd0zn =EvalGraph(fD0ZResPiUpg ,ptmc);
290 spt1n =EvalGraph(fPt1ResPiUpg ,ptmc);
293 sd0rpo=EvalGraph(fD0RPResECur,ptmc);
294 sd0zo =EvalGraph(fD0ZResECur ,ptmc);
295 spt1o =EvalGraph(fPt1ResECur ,ptmc);
296 sd0rpn=EvalGraph(fD0RPResEUpg,ptmc);
297 sd0zn =EvalGraph(fD0ZResEUpg ,ptmc);
298 spt1n =EvalGraph(fPt1ResEUpg ,ptmc);
304 //mj special to check resolution smeared by 10 %
305 sd0rpn = sd0rpo * 1.1;
309 // Use the same units (i.e. cm and GeV/c)! TODO: pt!
315 // Apply the smearing
316 Double_t d0zo =param [1];
317 Double_t d0zmc =parammc[1];
318 Double_t d0rpo =param [0];
319 Double_t d0rpmc=parammc[0];
320 Double_t pt1o =param [4];
321 Double_t pt1mc =parammc[4];
322 Double_t dd0zo =d0zo-d0zmc;
323 Double_t dd0zn =dd0zo *(sd0zo >0. ? (sd0zn /sd0zo ) : 1.);
324 Double_t d0zn =d0zmc+dd0zn;
325 Double_t dd0rpo=d0rpo-d0rpmc;
326 Double_t dd0rpn=dd0rpo*(sd0rpo>0. ? (sd0rpn/sd0rpo) : 1.);
327 Double_t d0rpn =d0rpmc+dd0rpn;
328 Double_t dpt1o =pt1o-pt1mc;
329 Double_t dpt1n =dpt1o *(spt1o >0. ? (spt1n /spt1o ) : 1.);
330 Double_t pt1n =pt1mc+dpt1n;
336 //printf("d0rpo= %lf d0rpmc= %lf dd0rpo= %lf dd0rpn= %lf\n",d0rpo,d0rpmc,dd0rpo,dd0rpn);
337 // Copy the smeared parameters to the AOD track
342 // printf("x[0]= %lf x[1]= %lf x[2]= %lf \n",x[0],x[1],x[2]);
343 // track->SetPosition(x,kFALSE);
344 // track->SetP(p,kTRUE);
346 AliESDEvent *esdevent = dynamic_cast<AliESDEvent *>(fEvent);
347 const AliVVertex *vtxESDSkip = esdevent->GetPrimaryVertex();
348 // recalculate primary vertex for peri. and pp
349 AliVertexerTracks vertexer(fEvent->GetMagneticField());
350 vertexer.SetITSMode();
351 vertexer.SetMinClusters(4);
353 skipped[0] = track->GetID();
354 vertexer.SetSkipTracks(1,skipped);
355 vertexer.SetConstraintOn();
356 //----diamond constraint-----------------------------
357 Float_t diamondcovxy[3];
358 esdevent->GetDiamondCovXY(diamondcovxy);
359 Double_t pos[3]={esdevent->GetDiamondX(),esdevent->GetDiamondY(),0.};
360 Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
361 AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
362 vertexer.SetVtxStart(diamond);
363 delete diamond; diamond=NULL;
364 //----------------------------------------------------
365 vtxESDSkip = vertexer.FindPrimaryVertex(fEvent);
368 // Propagation always done on a working copy to not disturb the track params of the original track
369 AliESDtrack *esdtrack = NULL;
370 if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
371 // Case ESD track: take copy constructor
372 AliESDtrack *tmptrack = dynamic_cast<AliESDtrack *>(track);
373 if(tmptrack) esdtrack = new AliESDtrack(*tmptrack);
375 // Case AOD track: take different constructor
376 esdtrack = new AliESDtrack(track);
379 et.PropagateToDCA(tmpvtx, fEvent->GetMagneticField(), kBeampiperadiusnew, dcan, covn);
380 //et.PropagateToDCA(vtxESDSkip, fEvent->GetMagneticField(), kBeampiperadiusnew, dcan, covn);
381 esdtrack->PropagateToDCA(vtxESDSkip, fEvent->GetMagneticField(), kBeampiperadius, dcao, covo);
386 Double_t resfactor=1., resfactorz=1.;
387 if(sd0rpo) resfactor = sd0rpn/sd0rpo;
388 if(sd0zo) resfactorz = sd0zn/sd0zo;
389 if(covo[0]) dcaxySign = dcan[0]/(TMath::Sqrt(covo[0])*resfactor);
390 if(covo[0]) dcaxySigo = dcao[0]/TMath::Sqrt(covo[0]);
391 if(covo[2]) dcazSign = dcan[1]/(TMath::Sqrt(covo[2])*resfactorz);
392 if(covo[2]) dcazSigo = dcao[1]/TMath::Sqrt(covo[2]);
393 //if(dcaxyo) printf("pt= %lf resolusion xy= %lf dcaxyo= %lf dcaxyn= %lf ratio= %lf\n",pt1mc,resfactor,dcaxyo,dcaxyn,dcaxyn/dcaxyo);
409 //______________________________________________________
410 Double_t AliHFEsmearDCA::EvalGraph(const TGraph *graph,Double_t x) const {
412 // Evaluates a TGraph without linear extrapolation. Instead the last
413 // valid point of the graph is used when out of range.
414 // The function assumes an ascending order of X.
417 // TODO: find a pretty solution for this:
418 Int_t n =graph->GetN();
419 Double_t xmin=graph->GetX()[0 ];
420 Double_t xmax=graph->GetX()[n-1];
421 if (x<xmin) return graph->Eval(xmin);
422 if (x>xmax) return graph->Eval(xmax);
423 return graph->Eval(x);