PWG1/AliESDresolMakerFast.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 //   ESD track and V0 fast resolution parameterization
  19 //   Fast algorithm to make parameterization
  20 //   Track covariance used
  21 //
  22 //
  23 //    Origin: Marian Ivanov marian.ivanov@cern.ch
  24 //-------------------------------------------------------------------------
  25
  26 /*
  27   EXAMPLE USAGE:
  28   //
  29   // Make esd chain
  30   //
  31   .x ~/rootlogon.C
  32   gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
  33   gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+");
  34   AliXRDPROOFtoolkit tool;
  35   TChain * tree = tool.MakeChain("esd.txt","esdTree",0,1000);
  36   tree->Lookup();
  37   //
  38   // Load macros
  39   //
  40   gSystem->Load("libSTAT.so");
  41   .L $ALICE_ROOT/PWG1/AliESDresolParams.cxx+
  42   .L $ALICE_ROOT/PWG1/AliESDresolMakerFast.cxx+
  43   TCut cutDCA="Tracks[].fCchi2<100&&abs(Tracks[].fP[4])<8&&abs(Tracks[].fP[3])<1&&sqrt(Tracks[].fC[0])/(0.2+abs(Tracks[].fP[4]))<0.02&&abs(Tracks[].fX)<3&&Tracks[].fITSncls>4&&Tracks.fTPCncls>40"
  44   //
  45   // Create resolution
  46   //
  47   TObjArray * array = AliESDresolMakerFast::MakeParamPrimFast(tree,cutDCA,0.95,2000);
  48   AliESDresolParams params;
  49   params.SetResolPrimFast(array);
  50   params.SetInstance(&params);
  51   TFile f("resolParams.root","recreate");
  52   param.Write("resolParams");
  53   //
  54   //
  55   TF2 f2sy("f2sy","AliESDresolParams::SGetResolPrimFast(0,x,y)",0,10,0,2);
  56   TF1 f1sy("f1sy","AliESDresolParams::SGetResolPrimFast(0,x,0)",0,10);
  57   f2sy->Draw("surf2")
  58   f1sy->Draw();
  59
  60 */
  61
  62
  63
  64 #include "TVectorD.h"
  65 #include "../STAT/TStatToolkit.h"
  66 #include "TMath.h"
  67 #include "TCut.h"
  68 #include "TTree.h"
  69
  70 #include "AliESDresolParams.h"
  71 #include "AliESDresolMakerFast.h"
  72
  73
  74 ClassImp(AliESDresolMakerFast)
  75
  76
  77
  78
  79 AliESDresolMakerFast::AliESDresolMakerFast() :
  80   TObject()
  81 {
  82   //
  83   // Default constructor
  84   //
  85 }
  86
  87
  88 TObjArray * AliESDresolMakerFast::MakeParamPrimFast(TTree * tree, TCut & cutDCA, Float_t fraction, Int_t entries){
  89   //
  90   // DCA resolution linear parameterization
  91   // Only valid in ITS acceptance
  92   // Arguments:
  93   // tree     - esdTree or chain
  94   // cutDCA   - track selection criteria
  95   // fraction - robust fitter fraction
  96   // entries  - total number of entries (tracks) used
  97
  98   /* Example
  99      // its adjusted cuts
 100
 101      TCut cut1pt= "sqrt(Tracks[].fC[14])/(1+abs(Tracks[].fP[4]))^2.5<0.01";
 102      TCut cutsy = "sqrt(Tracks[].fC[0])/(0.2+abs(Tracks[].fP[4]))<0.02";
 103
 104      TCut cutDCA="Tracks[].fCchi2<100&&abs(Tracks[].fP[4])<8&&abs(Tracks[].fP[3])<1&&abs(Tracks[].fX)<3&&Tracks[].fITSncls>4&&Tracks.fTPCncls>40"+cut1pt+cutsy;
 105      fraction =0.90;
 106      entries  = 1000;
 107   */
 108
 109   TObjArray *array = new TObjArray();
 110   Double_t chi2=0;
 111   Int_t    npoints=0;
 112   TVectorD fitParam;
 113   TMatrixD covMatrix;
 114   //
 115   // y param
 116   //
 117   TString * dcayParam=  TStatToolkit::FitPlane(tree,"sqrt(Tracks[].fC[0])/(0.2+abs(Tracks[].fP[4]))","(abs(Tracks[].fP[4]))++(abs(Tracks[].fP[4]))^2++(abs(Tracks[].fP[3]))++(abs(Tracks[].fP[3]))^2++(abs(Tracks[].fP[4]))*(abs(Tracks[].fP[3]))", cutDCA, chi2,npoints,fitParam,covMatrix,fraction, 0, entries);
 118   tree->SetAlias("dcayParam",dcayParam->Data());
 119   array->AddAt(new TVectorD(fitParam),0);
 120   printf("Y resol\t%s\n",dcayParam->Data());
 121   //
 122   // z param
 123   //
 124   TString * dcazParam= TStatToolkit::FitPlane(tree, "sqrt(Tracks[].fC[2])/(0.2+abs(Tracks[].fP[4]))", "(abs(Tracks[].fP[4]))++(abs(Tracks[].fP[4]))^2++(abs(Tracks[].fP[3]))++(abs(Tracks[].fP[3]))^2++(abs(Tracks[].fP[4]))*(abs(Tracks[].fP[3]))",  cutDCA, chi2,npoints,fitParam,covMatrix,fraction,0,entries);
 125   printf("Z resol\t%s\n",dcazParam->Data());
 126   tree->SetAlias("dcazParam",dcazParam->Data());
 127   array->AddAt(new TVectorD(fitParam),1);
 128   //
 129   // Phi param
 130   //
 131   TString * dcaphiParam= TStatToolkit::FitPlane(tree, "sqrt(Tracks[].fC[5])/(0.1+abs(Tracks[].fP[4]))", "(abs(Tracks[].fP[4]))++(abs(Tracks[].fP[4]))^2++(abs(Tracks[].fP[3]))++(abs(Tracks[].fP[3]))^2++(abs(Tracks[].fP[4]))*(abs(Tracks[].fP[3]))",  cutDCA, chi2,npoints,fitParam,covMatrix,fraction,0,entries);
 132   printf("Phi resol\t%s\n",dcaphiParam->Data());
 133   tree->SetAlias("dcaphiParam",dcaphiParam->Data());
 134   array->AddAt(new TVectorD(fitParam),2);
 135   //
 136   // theta param
 137   //
 138   TString * dcathParam=  TStatToolkit::FitPlane(tree, "sqrt(Tracks[].fC[9])/((0.1+abs(Tracks[].fP[4])*(1+abs(Tracks[].fParamP.fP[3]^2))))", "(abs(Tracks[].fP[4]))++(abs(Tracks[].fP[4]))^2++(abs(Tracks[].fP[3]))++(abs(Tracks[].fP[3]))^2++(abs(Tracks[].fP[4]))*(abs(Tracks[].fP[3]))",  cutDCA, chi2,npoints,fitParam,covMatrix,fraction,0,entries);
 139   printf("Theta resol\t%s\n",dcathParam->Data());
 140   tree->SetAlias("dcathParam",dcathParam->Data());
 141   array->AddAt(new TVectorD(fitParam),3);
 142   //
 143   // 1pt param
 144   //
 145   TString * dca1ptParam=  TStatToolkit::FitPlane(tree, "sqrt(Tracks[].fC[14])/(1+abs(Tracks[].fP[4]))^2", "(abs(Tracks[].fP[4]))++(abs(Tracks[].fP[4]))^2++(abs(Tracks[].fP[3]))++(abs(Tracks[].fP[3]))^2++(abs(Tracks[].fP[4]))*(abs(Tracks[].fP[3]))",  cutDCA, chi2,npoints,fitParam,covMatrix,fraction,0,entries);
 146   printf("1pt resol\t%s\n",dca1ptParam->Data());
 147   tree->SetAlias("dca1ptParam",dca1ptParam->Data());
 148   array->AddAt(new TVectorD(fitParam),4);
 149   return array;
 150 }
 151
 152 /*
 153 void scalePt(){
 154   //
 155   TString * dca1ptParamS4=  TStatToolkit::FitPlane(tree, "(Tracks[].fC[14])", "abs(Tracks[].fP[4])^4",  cutDCA, chi2,npoints,fitParam,covMatrix,fraction,0,entries);
 156   printf("1pt resol\t%s\n",dca1ptParamS4->Data());
 157   tree->SetAlias("dca1ptParam4",dca1ptParamS4->Data());
 158   //
 159   TString * dca1ptParamS34=  TStatToolkit::FitPlane(tree, "(Tracks[].fC[14])", "abs(Tracks[].fP[4])^4++abs(Tracks[].fP[4])^3",  cutDCA, chi2,npoints,fitParam,covMatrix,fraction,0,entries);
 160   printf("1pt resol\t%s\n",dca1ptParamS34->Data());
 161   tree->SetAlias("dca1ptParam34",dca1ptParamS34->Data());
 162   //
 163   TString * dca1ptParamS234=  TStatToolkit::FitPlane(tree, "(Tracks[].fC[14])", "abs(Tracks[].fP[4])^4++abs(Tracks[].fP[4])^3++abs(Tracks[].fP[4])^2",  cutDCA, chi2,npoints,fitParam,covMatrix,fraction,0,entries);
 164   printf("1pt resol\t%s\n",dca1ptParamS234->Data());
 165   tree->SetAlias("dca1ptParam234",dca1ptParamS234->Data());
 166 }
 167
 168 */
 169
 170
 171
 172
 173
 174
 175 TObjArray * AliESDresolMakerFast::MakeParamRFast(TTree * tree,  TCut &cutV0, Float_t fraction, Int_t entries){
 176   //
 177   //
 178   //
 179   TObjArray *array = new TObjArray;
 180
 181   Double_t chi2=0;
 182   Int_t    npoints=0;
 183   TVectorD fitParam;
 184   TMatrixD covMatrix;
 185   //
 186   //
 187   /*
 188     //TCut cutGood="abs(V0s[].GetEffMass(0,0))<0.05||abs(V0s[].GetEffMass(2,2)-0.5)<0.05"
 189     TCut  cutV0 = "abs((V0s[].fParamP.fC[0]))<3&&abs(V0s[].fParamP.fP[3])<1&&abs(V0s[].fParamP.fP[4])<8";//
 190   */
 191   //
 192   //
 193   //
 194   TString * v0sigmaY= TStatToolkit::FitPlane(tree,"sqrt(sqrt((V0s[].fParamP.fC[0]))/(0.2+abs(V0s[].fParamP.fP[4])))","abs(V0s[].fParamP.fP[4])++V0s[].fParamP.fX++abs(V0s[].fParamP.fP[4])^2++V0s[].fParamP.fX^2++abs(V0s[].fParamP.fP[4])*V0s[].fParamP.fX++abs(V0s[].fParamP.fP[4])*V0s[].fParamP.fX^2", cutV0, chi2,npoints,fitParam,covMatrix,fraction,0,entries);
 195   tree->SetAlias("v0sigmaY",v0sigmaY->Data());
 196   array->AddAt(new TVectorD(fitParam),0);
 197   //
 198   //
 199   //
 200   TString * v0sigmaZ= TStatToolkit::FitPlane(tree,"sqrt(sqrt((V0s[].fParamP.fC[2])))","abs(V0s[].fParamP.fP[4])++V0s[].fParamP.fX++abs(V0s[].fParamP.fP[4])^2++V0s[].fParamP.fX^2++abs(V0s[].fParamP.fP[4])*V0s[].fParamP.fX++abs(V0s[].fParamP.fP[4])*V0s[].fParamP.fX^2", cutV0, chi2,npoints,fitParam,covMatrix,fraction,0,entries);
 201   tree->SetAlias("v0sigmaZ",v0sigmaZ->Data());
 202   array->AddLast(new TVectorD(fitParam));
 203   //
 204   //
 205   //
 206   TString * v0sigmaPhi= TStatToolkit::FitPlane(tree,"sqrt(sqrt((V0s[].fParamP.fC[5]))/(0.1+abs(V0s[].fParamP.fP[4])))","abs(V0s[].fParamP.fP[4])++V0s[].fParamP.fX++abs(V0s[].fParamP.fP[4])^2++abs(V0s[].fParamP.fP[4])*V0s[].fParamP.fX", cutV0, chi2,npoints,fitParam,covMatrix,fraction,0,entries);
 207   tree->SetAlias("v0sigmaPhi",v0sigmaPhi->Data());
 208   array->AddAt(new TVectorD(fitParam),2);
 209   //
 210   //
 211   //
 212   TString * v0sigmaTh= TStatToolkit::FitPlane(tree,"sqrt(sqrt((V0s[].fParamP.fC[9]))/((0.1+abs(V0s[].fParamP.fP[4])*(1+abs(V0s[].fParamP.fP[3])^2))))","abs(V0s[].fParamP.fP[4])++V0s[].fParamP.fX++abs(V0s[].fParamP.fP[4])^2++abs(V0s[].fParamP.fP[4])*V0s[].fParamP.fX", cutV0, chi2,npoints,fitParam,covMatrix,fraction,0,entries);
 213   tree->SetAlias("v0sigmaTh",v0sigmaTh->Data());
 214   array->AddAt(new TVectorD(fitParam),3);
 215   //
 216   //
 217   //
 218   TString * v0sigma1pt= TStatToolkit::FitPlane(tree,"sqrt(sqrt((V0s[].fParamP.fC[14])))","abs(V0s[].fParamP.fP[4])++V0s[].fParamP.fX++abs(V0s[].fParamP.fP[4])^2++V0s[].fParamP.fX^2++abs(V0s[].fParamP.fP[4])*V0s[].fParamP.fX++abs(V0s[].fParamP.fP[4])*V0s[].fParamP.fX^2", cutV0, chi2,npoints,fitParam,covMatrix,fraction,0,entries);
 219   tree->SetAlias("v0sigma1pt",v0sigma1pt->Data());
 220   array->AddAt(new TVectorD(fitParam),4);
 221   return array;
 222
 223 }