Changing to centrality dependent corrections
[u/mrichter/AliRoot.git] / ANALYSIS / TenderSupplies / AliTrackFixTenderSupply.cxx
CommitLineData
4d44dfb6 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// Redo primary vertex on the fly, using the diamond constraint //
20// //
21///////////////////////////////////////////////////////////////////////////////
22
23#include <TSystem.h>
24#include "AliESDEvent.h"
25#include "AliTender.h"
26#include "AliVParticle.h"
27#include "AliLog.h"
28#include "AliOADBContainer.h"
29#include "AliAnalysisManager.h"
30#include "AliTrackFixTenderSupply.h"
31#include "AliOADBTrackFix.h"
32
33ClassImp(AliTrackFixTenderSupply)
34
35AliTrackFixTenderSupply::AliTrackFixTenderSupply() :
1112c329 36fDebug(0),
37 fBz(0),
4d44dfb6 38 fParams(0),
1112c329 39 fOADBObjPath("$OADB/PWGPP/data/CorrPTInv.root"),
4d44dfb6 40 fOADBObjName("CorrPTInv"),
41 fOADBCont(0)
42{
43 // default ctor
44}
45
46//_____________________________________________________
47AliTrackFixTenderSupply::AliTrackFixTenderSupply(const char *name, const AliTender *tender) :
48 AliTenderSupply(name,tender),
1112c329 49 fDebug(0),
4d44dfb6 50 fBz(0),
51 fParams(0),
1112c329 52 fOADBObjPath("$OADB/PWGPP/data/CorrPTInv.root"),
4d44dfb6 53 fOADBObjName("CorrPTInv"),
54 fOADBCont(0)
55{
56 // named ctor
57 //
58}
59
60//_____________________________________________________
61AliTrackFixTenderSupply::~AliTrackFixTenderSupply()
62{
63 // d-tor
64 delete fOADBCont;
65}
66
67
68//_____________________________________________________
69void AliTrackFixTenderSupply::ProcessEvent()
70{
71 //
72 // Fix track kinematics
73 //
74 AliESDEvent *event=fTender->GetEvent();
75 if (!event) return;
76 //
77 if (fTender->RunChanged() && !GetRunCorrections(fTender->GetRun())) return;
78 //
79 int nTracks = event->GetNumberOfTracks();
80 fBz = event->GetMagneticField();
81 if (TMath::Abs(fBz) < kAlmost0Field) return;
82 //
83 const AliESDVertex *vtx = 0, *vtxTPC = 0;
84 //
85 vtx = event->GetPrimaryVertexTracks(); // vertex to be used for update via RelateToVertex
86 if (!vtx || vtx->GetStatus()<1) {
87 vtx = event->GetPrimaryVertexSPD();
88 if (vtx && vtx->GetStatus()<1) vtx = 0;
89 }
90 vtxTPC = event->GetPrimaryVertexTPC(); // vertex to be used for update via RelateToVertexTPC
1112c329 91 if (vtxTPC && vtxTPC->GetStatus()<1) vtxTPC = 0;
4d44dfb6 92 //
93 AliExternalTrackParam* extPar = 0;
94 double xOrig = 0;
51732659 95 double xyzTPCInner[3] = {0,0,0};
4d44dfb6 96 for (int itr=0;itr<nTracks;itr++) {
97 //
98 AliESDtrack* trc = event->GetTrack(itr);
99 if (!trc->IsOn(AliESDtrack::kTPCin)) continue;
100 //
101 double sideAfraction = GetSideAFraction(trc);
102 // correct the main parameterization
103 int cormode = trc->IsOn(AliESDtrack::kITSin) ? AliOADBTrackFix::kCorModeGlob : AliOADBTrackFix::kCorModeTPCInner;
104 xOrig = trc->GetX();
105 double xIniCor = fParams->GetXIniPtInvCorr(cormode);
106 const AliExternalTrackParam* parInner = trc->GetInnerParam();
107 if (!parInner) {
108 AliError("Failed to extract inner param");
109 continue;
110 }
51732659 111 parInner->GetXYZ(xyzTPCInner);
e505d150 112 double phi = TMath::ATan2(xyzTPCInner[1],xyzTPCInner[0]);
113 if (phi<0) phi += 2*TMath::Pi();
4d44dfb6 114 //
1112c329 115 if (fDebug>1) {
116 AliInfo(Form("Tr:%4d kITSin:%d Phi=%+5.2f at X=%+7.2f | SideA fraction: %.3f",itr,trc->IsOn(AliESDtrack::kITSin),phi,parInner->GetX(),sideAfraction));
117 AliInfo(Form("Main Param before corr. in mode %s, xIni:%.1f",cormode== AliOADBTrackFix::kCorModeGlob ? "Glo":"TPC",xIniCor));
118 trc->AliExternalTrackParam::Print();
119 }
120 //
4d44dfb6 121 if (xIniCor>0) trc->PropagateTo(xIniCor,fBz);
122 CorrectTrackPtInv(trc, cormode, sideAfraction, phi);
123 if (xIniCor>0) { // full update is requested
124 if (vtx) trc->RelateToVertex(vtx, fBz, kVeryBig); // redo DCA if vtx is available
125 else trc->PropagateTo(xOrig, fBz); // otherwise bring to original point
126 }
127 //
1112c329 128 if (fDebug>1) {
129 AliInfo("Main Param after corr.");
130 trc->AliExternalTrackParam::Print();
131 }
4d44dfb6 132 // correct TPCinner param
133 if ( (extPar=(AliExternalTrackParam*)trc->GetTPCInnerParam()) ) {
134 cormode = AliOADBTrackFix::kCorModeTPCInner;
135 xOrig = extPar->GetX();
136 xIniCor = fParams->GetXIniPtInvCorr(cormode);
1112c329 137 if (fDebug>1) {
138 AliInfo(Form("TPCinner Param before corr. in mode %s, xIni:%.1f",cormode== AliOADBTrackFix::kCorModeGlob ? "Glo":"TPC",xIniCor));
139 extPar->AliExternalTrackParam::Print();
140 }
141 //
4d44dfb6 142 if (xIniCor>0) extPar->PropagateTo(xIniCor,fBz);
143 CorrectTrackPtInv(extPar,cormode,sideAfraction, phi);
144 if (xIniCor>0) { // full update is requested
145 if (vtxTPC) trc->RelateToVertexTPC(vtxTPC, fBz, kVeryBig); // redo DCA if vtx is available
146 else extPar->PropagateTo(xOrig, fBz); // otherwise bring to original point
147 }
1112c329 148 //
149 if (fDebug>1) {
150 AliInfo("TPCinner Param after corr.");
151 extPar->AliExternalTrackParam::Print();
152 }
4d44dfb6 153 }
154 //
155 }
156 //
157}
158
159//_____________________________________________________
160void AliTrackFixTenderSupply::CorrectTrackPtInv(AliExternalTrackParam* trc, Int_t mode, double sideAfraction, double phi) const
161{
162 // fix track kinematics
163 if (!trc) return;
164 double *param = (double*)trc->GetParameter();
165 param[4] += fParams->GetPtInvCorr(mode,sideAfraction,phi);
166 //
167}
168
169//_____________________________________________________
170Bool_t AliTrackFixTenderSupply::LoadOADBObjects()
171{
172 // Load OADB parameters
173 TString fileName = fOADBObjPath;
174 if (fileName.BeginsWith("$OADB")) fileName.ReplaceAll("$OADB",Form("%s/",AliAnalysisManager::GetOADBPath()));
175 gSystem->ExpandPathName(fileName);
176 AliInfo(Form("Loading correction parameters %s from %s",fOADBObjName.Data(),fileName.Data()));
177 //
178 fOADBCont = new AliOADBContainer("OADB");
179 if (fOADBCont->InitFromFile(fileName.Data(),fOADBObjName.Data())) {
180 AliError("Failed to load OADB Container");
181 delete fOADBCont;
182 fOADBCont = 0;
183 return kFALSE;
184 }
185 //
186 return kTRUE;
187}
188
189//_____________________________________________________
190Bool_t AliTrackFixTenderSupply::GetRunCorrections(int run)
191{
192 // extract corrections for given run
193 fParams = 0;
194 if (!fOADBCont) if (!LoadOADBObjects()) return kFALSE;
195 fParams = dynamic_cast<AliOADBTrackFix*>(fOADBCont->GetObject(run,"default"));
196 if (!fParams) {AliError(Form("No correction parameters for found for run %d",run)); return kFALSE;}
197 AliInfo(Form("Loaded correction parameters for run %d",run));
198 //
199 return kTRUE;
200}
201
202//_____________________________________________________
203Double_t AliTrackFixTenderSupply::GetSideAFraction(const AliESDtrack* track) const
204{
205 // decide A-side fraction of trajectory length in TPC
206 // for the moment just a placeholder...
207 //
208 const AliExternalTrackParam *trIn = track->GetInnerParam();
209 if (!trIn) {AliError("Failed to onbtain InnerParam"); return 0.5;}
210 //
211 // easy cases:
212 double sideAfrac = 0.5;
213 if (trIn->GetZ()>0) { sideAfrac = 1.; if (trIn->GetTgl()>-kAlmost0) return sideAfrac;}
214 else if (trIn->GetZ()<0) { sideAfrac = 0.; if (trIn->GetTgl()< kAlmost0) return sideAfrac;}
215 //
216 // track may cross the membrane
217 // convert to usual helix param:
218 // x(t)=(dca+R)*cos(phi0) -R*cos(t+phi0);
219 // y(t)=(dca+R)*sin(phi0) -R*sin(t+phi0);
220 // z(t)=z0 - R*t*tgl
221 //
222 int sqb = fBz<0 ? -trIn->Charge() : trIn->Charge();
223 double radq = TMath::Abs(trIn->GetC(fBz));
224 radq = radq>kAlmost0 ? 1./radq : 1e4; // track radius
225 if (sqb<0) radq = -radq;
226 double cstCphi0 = trIn->GetSnp(); // tC is current t
227 double sntCphi0 = -TMath::Sqrt((1.-cstCphi0)*(1+cstCphi0));
228 double xc0 = trIn->GetX() + radq*cstCphi0; // circle center
229 double yc0 = trIn->GetY() + radq*sntCphi0;
230 double phi0 = TMath::ATan2( yc0, xc0);
231 if (sqb<0) phi0 += TMath::Pi();
232 if (phi0 > TMath::Pi()) phi0 -= 2.*TMath::Pi();
233 else if (phi0 <-TMath::Pi()) phi0 += 2.*TMath::Pi();
234 //
235 double csphi0 = TMath::Cos(phi0);
236 double snphi0 = TMath::Sin(phi0);
237 double dca = xc0*csphi0 + yc0*snphi0 - radq;
238 double tC = TMath::ATan2(sntCphi0,cstCphi0) - phi0;
239 double tglR = radq*trIn->GetTgl();
240 double z0 = trIn->GetZ() + tglR*tC;
241 Double_t tCross = z0/tglR; // Z=0 is crossed at this t
242 // t is changed withing tCross:2pi range for the track moving forward
243 if (tCross<tC || tCross>TMath::Pi()) return sideAfrac; // on the way forward does not cross the membrane
244 //
245 // find max t reached by the track in the TPC
246 const double kROutTPC2 = 250.*250.;
247 // find t at which the outer radius is crossed
248 // r^2 = (dca+R)^2 + R^2 - 2R(R+dca)*cos(t)
249 double dcaR = dca+radq;
250 double rmax2 = dcaR*dcaR + radq*radq + 2*radq*dcaR; // max rad at t=pi
251 double tmax = TMath::Pi();
252 if (rmax2>kROutTPC2) { // traverses full TPC
253 double ct = 0.5*(dcaR/radq + radq/dcaR - kROutTPC2/radq/dcaR);
254 if (TMath::Abs(ct)<1-kAlmost0) tmax = TMath::ACos(ct);
255 }
256 if (tmax<tCross) return sideAfrac; // does not cross z=0
257 sideAfrac = sideAfrac>0.5 ? (tCross-tC)/(tmax-tC) : 1.-(tCross-tC)/(tmax-tC);
258 return sideAfrac;
259}