]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGHF/hfe/AliHFEsmearDCA.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliHFEsmearDCA.cxx
CommitLineData
959ea9d8 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// Get improved dca info. by ITS upgrade implemented by the ALICE Heavy Flavour Electron Group
17//
18// Authors:
19// MinJung Kweon <minjung@physi.uni-heidelberg.de>
20//
21#include <TBits.h>
22#include <TClass.h>
23#include <TH1F.h>
24#include <TH2F.h>
25#include <TList.h>
26#include <TString.h>
27#include <TMath.h>
28#include <TGraph.h>
29#include <TFile.h>
30
31#include "AliAODTrack.h"
32#include "AliAODPid.h"
33#include "AliESDEvent.h"
34#include "AliESDVertex.h"
35#include "AliESDtrack.h"
36#include "AliLog.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"
45
46#include "AliHFEsmearDCA.h"
47
48ClassImp(AliHFEsmearDCA)
49
50
51//______________________________________________________
52AliHFEsmearDCA::AliHFEsmearDCA(const Char_t */*name*/, const char *resfileCurURI, const char *resfileUpgURI, const Char_t */*title*/):
53 fEvent(NULL),
54 fMCEvent(NULL),
55 fD0ZResPCur (0),
56 fD0ZResKCur (0),
57 fD0ZResPiCur (0),
58 fD0ZResECur (0),
59 fD0RPResPCur (0),
60 fD0RPResKCur (0),
61 fD0RPResPiCur(0),
62 fD0RPResECur (0),
63 fPt1ResPCur (0),
64 fPt1ResKCur (0),
65 fPt1ResPiCur (0),
66 fPt1ResECur (0),
67 fD0ZResPUpg (0),
68 fD0ZResKUpg (0),
69 fD0ZResPiUpg (0),
70 fD0ZResEUpg (0),
71 fD0RPResPUpg (0),
72 fD0RPResKUpg (0),
73 fD0RPResPiUpg(0),
74 fD0RPResEUpg (0),
75 fPt1ResPUpg (0),
76 fPt1ResKUpg (0),
77 fPt1ResPiUpg (0),
78 fPt1ResEUpg (0)
79{
80 //
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.
87 //
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" );
113 delete resfileCur;
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" );
139 delete resfileUpg;
140
141}
142
143//______________________________________________________
144AliHFEsmearDCA::AliHFEsmearDCA(const AliHFEsmearDCA &c):
145 TObject(c),
146 fEvent(c.fEvent),
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)
172{
173 //
174 // Copy constructor
175 // Performs a deep copy
176 //
177}
178
179//______________________________________________________
180AliHFEsmearDCA::~AliHFEsmearDCA(){
181 //
182 // Destructor
183 //
184}
185
186//______________________________________________________
187void AliHFEsmearDCA::SetRecEventInfo(const TObject *event){
188 //
189 // Set Virtual event an make a copy
190 //
191 if (!event) {
192 AliError("Pointer to AliVEvent !");
193 return;
194 }
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 !");
198 return ;
199 }
200 fEvent = (AliVEvent*) event;
201
202}
203
204//______________________________________________________
205void 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){
206 //
207 // Get HFE impact parameter (with recalculated primary vertex)
208 //
209 if(!fEvent){
210 AliDebug(1, "No Input event available\n");
211 return;
212 }
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.};
220
221 // Get reconstructed track parameters
222
223 Double_t xyz[3],pxpypz[3],cv[21];
224 track->GetXYZ(xyz);
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();
230
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());
234
235 const AliMCParticle *mc=static_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(track->GetLabel())));
236 Double_t mcx[3];
237 Double_t mcp[3];
238 Double_t mccv[36]={0.};
239 Short_t mcc;
240 mc->XvYvZv(mcx);
241 mc->PxPyPz(mcp);
242 mcc=mc->Charge();
243 AliExternalTrackParam mct(mcx,mcp,mccv,mcc);
244 const Double_t *parammc=mct.GetParameter();
245 const AliVVertex *tmpvtx = fMCEvent->GetPrimaryVertex();
246
247 AliVertex vtx(mcx,1.,1);
248 AliVertex *tmpvtx2 = (AliVertex *)fMCEvent->GetPrimaryVertex();
249
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());
256
257
258 // Select appropriate smearing functions
259 Double_t ptmc=TMath::Abs(mc->Pt());
260 Double_t sd0rpn=0.;
261 Double_t sd0zn =0.;
262 Double_t spt1n =0.;
263 Double_t sd0rpo=0.;
264 Double_t sd0zo =0.;
265 Double_t spt1o =0.;
266
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);
275 break;
276 case 321: case -321:
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);
283 break;
284 case 211: case -211:
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);
291 break;
292 case 11: case -11:
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);
299 break;
300 default:
301 return;
302 }
303
304 //mj special to check resolution smeared by 10 %
305 sd0rpn = sd0rpo * 1.1;
306 sd0zn = sd0zo * 1.1;
307 spt1n = spt1o;
308
309 // Use the same units (i.e. cm and GeV/c)! TODO: pt!
310 sd0rpo*=1.e-4;
311 sd0zo *=1.e-4;
312 sd0rpn*=1.e-4;
313 sd0zn *=1.e-4;
314
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;
331
332 param[0]=d0rpn;
333 param[1]=d0zn ;
334 param[4]=pt1n ;
335
336 //printf("d0rpo= %lf d0rpmc= %lf dd0rpo= %lf dd0rpn= %lf\n",d0rpo,d0rpmc,dd0rpo,dd0rpn);
337 // Copy the smeared parameters to the AOD track
338 Double_t x[3];
339 Double_t p[3];
340 et.GetXYZ(x);
341 et.GetPxPyPz(p);
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);
345
346 AliESDEvent *esdevent = dynamic_cast<AliESDEvent *>(fEvent);
91e50e2b 347 if(!esdevent) return;
348 //const AliVVertex *vtxESDSkip = esdevent->GetPrimaryVertex();
959ea9d8 349 // recalculate primary vertex for peri. and pp
350 AliVertexerTracks vertexer(fEvent->GetMagneticField());
351 vertexer.SetITSMode();
352 vertexer.SetMinClusters(4);
353 Int_t skipped[2];
354 skipped[0] = track->GetID();
355 vertexer.SetSkipTracks(1,skipped);
356 vertexer.SetConstraintOn();
357 //----diamond constraint-----------------------------
358 Float_t diamondcovxy[3];
359 esdevent->GetDiamondCovXY(diamondcovxy);
360 Double_t pos[3]={esdevent->GetDiamondX(),esdevent->GetDiamondY(),0.};
361 Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
362 AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
363 vertexer.SetVtxStart(diamond);
364 delete diamond; diamond=NULL;
365 //----------------------------------------------------
91e50e2b 366 const AliVVertex *vtxESDSkip = (const AliVVertex *) vertexer.FindPrimaryVertex(fEvent);
367 //vtxESDSkip = vertexer.FindPrimaryVertex(fEvent);
959ea9d8 368
369 // Getting the DCA
370 // Propagation always done on a working copy to not disturb the track params of the original track
371 AliESDtrack *esdtrack = NULL;
372 if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
373 // Case ESD track: take copy constructor
374 AliESDtrack *tmptrack = dynamic_cast<AliESDtrack *>(track);
375 if(tmptrack) esdtrack = new AliESDtrack(*tmptrack);
376 } else {
377 // Case AOD track: take different constructor
378 esdtrack = new AliESDtrack(track);
379 }
380 if(esdtrack){
381 et.PropagateToDCA(tmpvtx, fEvent->GetMagneticField(), kBeampiperadiusnew, dcan, covn);
382 //et.PropagateToDCA(vtxESDSkip, fEvent->GetMagneticField(), kBeampiperadiusnew, dcan, covn);
383 esdtrack->PropagateToDCA(vtxESDSkip, fEvent->GetMagneticField(), kBeampiperadius, dcao, covo);
384 dcaxyn = dcan[0];
385 dcaxyo = dcao[0];
386 dcazn = dcan[1];
387 dcazo = dcao[1];
388 Double_t resfactor=1., resfactorz=1.;
389 if(sd0rpo) resfactor = sd0rpn/sd0rpo;
390 if(sd0zo) resfactorz = sd0zn/sd0zo;
391 if(covo[0]) dcaxySign = dcan[0]/(TMath::Sqrt(covo[0])*resfactor);
392 if(covo[0]) dcaxySigo = dcao[0]/TMath::Sqrt(covo[0]);
393 if(covo[2]) dcazSign = dcan[1]/(TMath::Sqrt(covo[2])*resfactorz);
394 if(covo[2]) dcazSigo = dcao[1]/TMath::Sqrt(covo[2]);
395 //if(dcaxyo) printf("pt= %lf resolusion xy= %lf dcaxyo= %lf dcaxyn= %lf ratio= %lf\n",pt1mc,resfactor,dcaxyo,dcaxyn,dcaxyn/dcaxyo);
396
397 if(!covo[0]){
398 dcaxySigo = -49.;
399 dcaxySign = -49.;
400 }
401 if(!covo[2]){
402 dcazSigo = -49.;
403 dcazSign = -49.;
404 }
405 }
406 delete esdtrack;
407 //delete vtxESDSkip;
408
409}
410
411//______________________________________________________
412Double_t AliHFEsmearDCA::EvalGraph(const TGraph *graph,Double_t x) const {
413 //
414 // Evaluates a TGraph without linear extrapolation. Instead the last
415 // valid point of the graph is used when out of range.
416 // The function assumes an ascending order of X.
417 //
418
419 // TODO: find a pretty solution for this:
420 Int_t n =graph->GetN();
421 Double_t xmin=graph->GetX()[0 ];
422 Double_t xmax=graph->GetX()[n-1];
423 if (x<xmin) return graph->Eval(xmin);
424 if (x>xmax) return graph->Eval(xmax);
425 return graph->Eval(x);
426}