don't lie in the log!
[u/mrichter/AliRoot.git] / PWGPP / TPC / AliESDresolMakerFast.cxx
CommitLineData
7cc34f08 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");
2bfe5463 41 .L $ALICE_ROOT/PWGPP/AliESDresolParams.cxx+
42 .L $ALICE_ROOT/PWGPP/AliESDresolMakerFast.cxx+
7cc34f08 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
74ClassImp(AliESDresolMakerFast)
75
76
77
78
79AliESDresolMakerFast::AliESDresolMakerFast() :
80 TObject()
81{
82 //
83 // Default constructor
84 //
85}
86
87
88TObjArray * 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/*
153void 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
175TObjArray * 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}