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