TENDER/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]);
 113     if (phi<0) phi += 2*TMath::Pi();
 114     //
 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     //
 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     //
 128     if (fDebug>1) {
 129       AliInfo("Main Param after corr.");
 130       trc->AliExternalTrackParam::Print();
 131     }
 132     // correct TPCinner param
 133     if ( (extPar=(AliExternalTrackParam*)trc->GetTPCInnerParam()) ) {
 134       cormode = AliOADBTrackFix::kCorModeTPCInner;
 135       xOrig = extPar->GetX();
 136       xIniCor = fParams->GetXIniPtInvCorr(cormode);
 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       //
 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       }
 148       //
 149       if (fDebug>1) {
 150         AliInfo("TPCinner Param after corr.");
 151         extPar->AliExternalTrackParam::Print();
 152       }
 153     }
 154     //
 155   }
 156   //
 157 }
 158
 159 //_____________________________________________________
 160 void 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 //_____________________________________________________
 170 Bool_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 //_____________________________________________________
 190 Bool_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 //_____________________________________________________
 203 Double_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 }