.so cleanup: removed from gSystem->Load()
[u/mrichter/AliRoot.git] / PWGPP / TPC / 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");
41   .L $ALICE_ROOT/PWGPP/AliESDresolParams.cxx+
42   .L $ALICE_ROOT/PWGPP/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 }