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 if(!esdevent) return;
348 //const AliVVertex *vtxESDSkip = esdevent->GetPrimaryVertex();
349 // recalculate primary vertex for peri. and pp
350 AliVertexerTracks vertexer(fEvent->GetMagneticField());
351 vertexer.SetITSMode();
352 vertexer.SetMinClusters(4);
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 //----------------------------------------------------
366 const AliVVertex *vtxESDSkip = (const AliVVertex *) vertexer.FindPrimaryVertex(fEvent);
367 //vtxESDSkip = vertexer.FindPrimaryVertex(fEvent);
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);
377 // Case AOD track: take different constructor
378 esdtrack = new AliESDtrack(track);
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);
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);
411 //______________________________________________________
412 Double_t AliHFEsmearDCA::EvalGraph(const TGraph *graph,Double_t x) const {
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.
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);