]>
Commit | Line | Data |
---|---|---|
96b3109a | 1 | //------------------------------------------------------------------------------ |
2 | // Implementation of AliComparisonRes class. It keeps information from | |
3 | // comparison of reconstructed and MC particle tracks. In addtion, | |
4 | // it keeps selection cuts used during comparison. The comparison | |
5 | // information is stored in the ROOT histograms. Analysis of these | |
6 | // histograms can be done by using Analyse() class function. The result of | |
3baa4bfd | 7 | // the analysis (histograms/graphs) are stored in the folder which is |
8 | // a data member of AliComparisonRes. | |
9 | // | |
96b3109a | 10 | // Author: J.Otwinowski 04/02/2008 |
11 | //------------------------------------------------------------------------------ | |
12 | ||
13 | /* | |
3baa4bfd | 14 | |
15 | // after running comparison task, read the file, and get component | |
26e10df5 | 16 | gROOT->LoadMacro("$ALICE_ROOT/PWG1/Macros/LoadMyLibs.C"); |
17 | LoadMyLibs(); | |
18 | ||
96b3109a | 19 | TFile f("Output.root"); |
35771050 | 20 | //AliComparisonRes * compObj = (AliComparisonRes*)f.Get("AliComparisonRes"); |
21 | AliComparisonRes * compObj = (AliComparisonRes*)cOutput->FindObject("AliComparisonRes"); | |
22 | ||
3baa4bfd | 23 | // analyse comparison data |
24 | compObj->Analyse(); | |
25 | ||
26 | // the output histograms/graphs will be stored in the folder "folderRes" | |
27 | compObj->GetAnalysisFolder()->ls("*"); | |
28 | ||
29 | // user can save whole comparison object (or only folder with anlysed histograms) | |
30 | // in the seperate output file (e.g.) | |
31 | TFile fout("Analysed_Res.root","recreate"); | |
32 | compObj->Write(); // compObj->GetAnalysisFolder()->Write(); | |
33 | fout.Close(); | |
96b3109a | 34 | |
96b3109a | 35 | */ |
36 | ||
96b3109a | 37 | #include "TCanvas.h" |
7b392278 | 38 | #include "TH1.h" |
39 | #include "TAxis.h" | |
40 | ||
41 | #include "AliComparisonRes.h" | |
42 | #include "AliESDRecInfo.h" | |
96b3109a | 43 | #include "AliESDVertex.h" |
7b392278 | 44 | #include "AliESDtrack.h" |
96b3109a | 45 | #include "AliLog.h" |
7b392278 | 46 | #include "AliMCInfo.h" |
47 | #include "AliMCInfoCuts.h" | |
48 | #include "AliRecInfoCuts.h" | |
96b3109a | 49 | #include "AliTracker.h" |
96b3109a | 50 | #include "AliTreeDraw.h" |
96b3109a | 51 | |
52 | using namespace std; | |
53 | ||
54 | ClassImp(AliComparisonRes) | |
55 | ||
56 | //_____________________________________________________________________________ | |
57 | AliComparisonRes::AliComparisonRes(): | |
3baa4bfd | 58 | AliComparisonObject("AliComparisonRes"), |
71a14197 | 59 | fResolHisto(0), |
60 | fPullHisto(0), | |
96b3109a | 61 | |
71a14197 | 62 | // Cuts |
63 | fCutsRC(0), | |
64 | fCutsMC(0), | |
96b3109a | 65 | |
71a14197 | 66 | // histogram folder |
67 | fAnalysisFolder(0) | |
68 | { | |
69 | //Init(); | |
70 | } | |
96b3109a | 71 | |
71a14197 | 72 | //_____________________________________________________________________________ |
73 | AliComparisonRes::AliComparisonRes(Char_t* name="AliComparisonRes", Char_t* title="AliComparisonRes",Int_t analysisMode=0,Bool_t hptGenerator=kFALSE): | |
74 | //AliComparisonRes::AliComparisonRes(Char_t* name, Char_t* title,Int_t analysisMode,Bool_t hptGenerator): | |
75 | AliComparisonObject(name,title), | |
76 | fResolHisto(0), | |
77 | fPullHisto(0), | |
96b3109a | 78 | |
96b3109a | 79 | // Cuts |
80 | fCutsRC(0), | |
3baa4bfd | 81 | fCutsMC(0), |
82 | ||
83 | // histogram folder | |
84 | fAnalysisFolder(0) | |
96b3109a | 85 | { |
71a14197 | 86 | // named constructor |
87 | // | |
88 | SetAnalysisMode(analysisMode); | |
89 | SetHptGenerator(hptGenerator); | |
90 | ||
3baa4bfd | 91 | Init(); |
96b3109a | 92 | } |
93 | ||
94 | //_____________________________________________________________________________ | |
71a14197 | 95 | AliComparisonRes::~AliComparisonRes() |
96 | { | |
97 | // destructor | |
98 | ||
99 | if(fResolHisto) delete fResolHisto; fResolHisto=0; | |
100 | if(fPullHisto) delete fPullHisto; fPullHisto=0; | |
96b3109a | 101 | |
71a14197 | 102 | if(fCutsRC) delete fCutsRC; fCutsRC=0; |
103 | if(fCutsMC) delete fCutsMC; fCutsMC=0; | |
96b3109a | 104 | |
3baa4bfd | 105 | if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0; |
96b3109a | 106 | } |
107 | ||
108 | //_____________________________________________________________________________ | |
3baa4bfd | 109 | void AliComparisonRes::Init(){ |
96b3109a | 110 | |
35771050 | 111 | // |
71a14197 | 112 | // histogram bining |
35771050 | 113 | // |
96b3109a | 114 | |
71a14197 | 115 | // set pt bins |
116 | Int_t nPtBins = 31; | |
117 | Double_t binsPt[32] = {0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.7,0.8,0.9,1.0,1.2,1.4,1.6,1.8,2.0,2.25,2.5,2.75,3.,3.5,4.,5.,6.,8.,10.}; | |
118 | Double_t ptMin = 0., ptMax = 10.; | |
96b3109a | 119 | |
96b3109a | 120 | |
71a14197 | 121 | if(IsHptGenerator() == kTRUE) { |
122 | nPtBins = 100; | |
123 | ptMin = 0.; ptMax = 100.; | |
124 | } | |
96b3109a | 125 | |
71a14197 | 126 | // res_y:res_z:res_phi,res_lambda:res_1pt:y:z:eta:phi:pt |
127 | Int_t binsResolHisto[10]={100,100,100,100,100,50,100,30,144,nPtBins}; | |
128 | Double_t minResolHisto[10]={-1.,-1.,-0.03,-0.03,-0.2,-0.03,-20.,-1.5,0.,ptMin}; | |
129 | Double_t maxResolHisto[10]={ 1., 1., 0.03, 0.03, 0.2, 0.03, 20., 1.5,2.*TMath::Pi(), ptMax}; | |
130 | ||
131 | ||
132 | fResolHisto = new THnSparseF("fResolHisto","res_y:res_z:res_phi:res_lambda:res_1pt:y:z:eta:phi:pt",10,binsResolHisto,minResolHisto,maxResolHisto); | |
133 | if(!IsHptGenerator()) fResolHisto->SetBinEdges(9,binsPt); | |
134 | ||
135 | fResolHisto->GetAxis(0)->SetTitle("res_y (cm)"); | |
136 | fResolHisto->GetAxis(1)->SetTitle("res_z (cm)"); | |
137 | fResolHisto->GetAxis(2)->SetTitle("res_phi (rad)"); | |
138 | fResolHisto->GetAxis(3)->SetTitle("res_lambda (rad)"); | |
139 | fResolHisto->GetAxis(4)->SetTitle("res_1pt (GeV/c)^{-1}"); | |
140 | fResolHisto->GetAxis(5)->SetTitle("y (cm)"); | |
141 | fResolHisto->GetAxis(6)->SetTitle("z (cm)"); | |
142 | fResolHisto->GetAxis(7)->SetTitle("eta"); | |
143 | fResolHisto->GetAxis(8)->SetTitle("phi (rad)"); | |
144 | fResolHisto->GetAxis(9)->SetTitle("pt (GeV/c)"); | |
145 | fResolHisto->Sumw2(); | |
146 | ||
147 | //pull_y:pull_z:pull_lambda:pull_phi:pull_1pt:y:z:eta:phi:pt | |
148 | Int_t binsPullHisto[10]={100,100,100,100,100,50,100,30,144,nPtBins}; | |
149 | Double_t minPullHisto[10]={-5.,-5.,-5.,-5.,-5.,-0.03, -20.,-1.5, 0., ptMin}; | |
150 | Double_t maxPullHisto[10]={ 5., 5., 5., 5., 5., 0.03, 20., 1.5, 2.*TMath::Pi(),ptMax}; | |
151 | ||
152 | fPullHisto = new THnSparseF("fPullHisto","pull_y:pull_z:pull_phi:pull_lambda:pull_1pt:y:z:eta:phi:pt",10,binsPullHisto,minPullHisto,maxPullHisto); | |
153 | if(!IsHptGenerator()) fPullHisto->SetBinEdges(9,binsPt); | |
154 | ||
155 | fPullHisto->GetAxis(0)->SetTitle("res_y (cm)"); | |
156 | fPullHisto->GetAxis(1)->SetTitle("res_z (cm)"); | |
157 | fPullHisto->GetAxis(2)->SetTitle("res_phi (rad)"); | |
158 | fPullHisto->GetAxis(3)->SetTitle("res_lambda (rad)"); | |
159 | fPullHisto->GetAxis(4)->SetTitle("res_1pt (GeV/c)^{-1}"); | |
160 | fPullHisto->GetAxis(5)->SetTitle("y (rad)"); | |
161 | fPullHisto->GetAxis(6)->SetTitle("z (rad)"); | |
162 | fPullHisto->GetAxis(7)->SetTitle("eta"); | |
163 | fPullHisto->GetAxis(8)->SetTitle("phi (rad)"); | |
164 | fPullHisto->GetAxis(9)->SetTitle("pt (GeV/c)"); | |
165 | fPullHisto->Sumw2(); | |
96b3109a | 166 | |
96b3109a | 167 | // Init cuts |
168 | if(!fCutsMC) | |
169 | AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object"); | |
170 | if(!fCutsRC) | |
171 | AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object"); | |
3baa4bfd | 172 | |
173 | // init folder | |
174 | fAnalysisFolder = CreateFolder("folderRes","Analysis Resolution Folder"); | |
96b3109a | 175 | } |
176 | ||
177 | //_____________________________________________________________________________ | |
71a14197 | 178 | void AliComparisonRes::ProcessTPC(AliMCInfo* const infoMC, AliESDRecInfo *const infoRC) |
96b3109a | 179 | { |
180 | // Fill resolution comparison information | |
181 | AliExternalTrackParam *track = 0; | |
ff0f1307 | 182 | //Double_t field = AliTracker::GetBz(); // nominal Bz field [kG] |
96b3109a | 183 | Double_t kMaxD = 123456.0; // max distance |
35771050 | 184 | AliESDVertex vertexMC; // MC primary vertex |
96b3109a | 185 | |
96b3109a | 186 | Double_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z |
187 | ||
71a14197 | 188 | Float_t delta1PtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC; |
189 | Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC; | |
190 | //Float_t delta1Pt2TPC, deltaY1PtTPC, deltaZ1PtTPC, deltaPhi1PtTPC, deltaTheta1PtTPC; | |
96b3109a | 191 | |
71a14197 | 192 | Float_t mceta = infoMC->GetParticle().Eta(); |
193 | Float_t mcphi = infoMC->GetParticle().Phi(); | |
194 | if(mcphi<0) mcphi += 2.*TMath::Pi(); | |
96b3109a | 195 | Float_t mcpt = infoMC->GetParticle().Pt(); |
196 | ||
197 | // distance to Prim. vertex | |
198 | const Double_t* dv = infoMC->GetVDist(); | |
199 | Bool_t isPrim = TMath::Sqrt(dv[0]*dv[0] + dv[1]*dv[1])<fCutsMC->GetMaxR() && TMath::Abs(dv[2])<fCutsMC->GetMaxVz(); | |
200 | ||
71a14197 | 201 | // Only 5 charged particle species (e,mu,pi,K,p) |
96b3109a | 202 | if (fCutsMC->IsPdgParticle(TMath::Abs(infoMC->GetParticle().GetPdgCode())) == kFALSE) return; |
203 | if (!isPrim) return; | |
71a14197 | 204 | //if (infoRC->GetStatus(1)!=3) return; // TPC refit |
96b3109a | 205 | if (!infoRC->GetESDtrack()) return; |
206 | if (infoRC->GetESDtrack()->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return; | |
207 | ||
208 | // calculate and set prim. vertex | |
35771050 | 209 | vertexMC.SetXv( infoMC->GetParticle().Vx() - dv[0] ); |
210 | vertexMC.SetYv( infoMC->GetParticle().Vy() - dv[1] ); | |
211 | vertexMC.SetZv( infoMC->GetParticle().Vz() - dv[2] ); | |
212 | ||
96b3109a | 213 | // calculate track parameters at vertex |
214 | const AliExternalTrackParam *innerTPC = 0; | |
215 | if ((innerTPC = infoRC->GetESDtrack()->GetTPCInnerParam()) != 0) | |
216 | { | |
217 | if ((track = new AliExternalTrackParam(*infoRC->GetESDtrack()->GetTPCInnerParam())) != 0 ) | |
218 | { | |
ff0f1307 | 219 | //Bool_t bDCAStatus = track->PropagateToDCA(&vertexMC,field,kMaxD,dca,cov); |
220 | Double_t field[3]; track->GetBxByBz(field); | |
221 | Bool_t bDCAStatus = track->PropagateToDCABxByBz(&vertexMC,field,kMaxD,dca,cov); | |
96b3109a | 222 | |
223 | // Fill parametrisation histograms (only TPC track) | |
3baa4bfd | 224 | if(bDCAStatus) |
71a14197 | 225 | { |
226 | // select primaries | |
227 | if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ()) | |
228 | { | |
229 | ||
230 | deltaYTPC= track->GetY()-infoMC->GetParticle().Vy(); | |
231 | deltaZTPC = track->GetZ()-infoMC->GetParticle().Vz(); | |
232 | deltaLambdaTPC = TMath::ATan2(track->Pz(),track->Pt())-TMath::ATan2(infoMC->GetParticle().Pz(),infoMC->GetParticle().Pt()); | |
233 | deltaPhiTPC = TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(infoMC->GetParticle().Py(),infoMC->GetParticle().Px()); | |
234 | delta1PtTPC = (track->OneOverPt()-1./mcpt)*mcpt; | |
235 | ||
236 | pullYTPC= (track->GetY()-infoMC->GetParticle().Vy()) / TMath::Sqrt(track->GetSigmaY2()); | |
237 | pullZTPC = (track->GetZ()-infoMC->GetParticle().Vz()) / TMath::Sqrt(track->GetSigmaZ2()); | |
238 | pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(track->GetSigmaTgl2()); | |
239 | pullPhiTPC = deltaPhiTPC / TMath::Sqrt(track->GetSigmaSnp2()); | |
240 | pull1PtTPC = (track->OneOverPt()-1./mcpt) / TMath::Sqrt(track->GetSigma1Pt2()); | |
241 | ||
242 | Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,delta1PtTPC,infoMC->GetParticle().Vy(),infoMC->GetParticle().Vz(),mceta,mcphi,mcpt}; | |
243 | fResolHisto->Fill(vResolHisto); | |
244 | ||
245 | Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,infoMC->GetParticle().Vy(),infoMC->GetParticle().Vz(),mceta,mcphi,mcpt}; | |
246 | fPullHisto->Fill(vPullHisto); | |
96b3109a | 247 | } |
71a14197 | 248 | |
249 | /* | |
250 | delta1Pt2TPC = (1/mcpt-innerTPC->OneOverPt())/TMath::Power(1+1/mcpt,2); | |
251 | deltaY1PtTPC= (infoMC->GetParticle().Vy()-innerTPC->GetY()) / (0.2+1/mcpt); | |
252 | deltaZ1PtTPC = (infoMC->GetParticle().Vz()-innerTPC->GetZ()) / (0.2+1/mcpt); | |
253 | deltaPhi1PtTPC = deltaPhiTPC / (0.1+1/mcpt); | |
254 | deltaTheta1PtTPC = deltaLambdaTPC / (0.1+1/mcpt); | |
255 | ||
256 | fPhiEtaPtTPC->Fill(TMath::ATan2(innerTPC->Py(),innerTPC->Px()), innerTPC->Eta(), innerTPC->Pt()); | |
257 | ||
258 | f1Pt2ResolS1PtTPC->Fill(s1mcpt,delta1Pt2TPC); | |
259 | fYResolS1PtTPC->Fill(s1mcpt,deltaY1PtTPC); | |
260 | fZResolS1PtTPC->Fill(s1mcpt,deltaZ1PtTPC); | |
261 | fPhiResolS1PtTPC->Fill(s1mcpt,deltaPhi1PtTPC); | |
262 | fThetaResolS1PtTPC->Fill(s1mcpt,deltaTheta1PtTPC); | |
263 | ||
264 | fPtResolPtTPC->Fill(mcpt,deltaPtTPC); | |
265 | fPtPullPtTPC->Fill(mcpt,pullPtTPC); | |
266 | fPhiResolEtaTPC->Fill(eta,deltaPhiTPC); | |
267 | fPhiPullEtaTPC->Fill(eta,pullPhiTPC); | |
268 | fLambdaResolEtaTPC->Fill(eta,deltaLambdaTPC); | |
269 | fLambdaPullEtaTPC->Fill(eta,pullLambdaTPC); | |
270 | */ | |
271 | } | |
272 | delete track; | |
96b3109a | 273 | } |
274 | } | |
71a14197 | 275 | } |
96b3109a | 276 | |
71a14197 | 277 | //_____________________________________________________________________________ |
278 | void AliComparisonRes::ProcessTPCITS(AliMCInfo* const infoMC, AliESDRecInfo *const infoRC) | |
279 | { | |
280 | Int_t clusterITS[200]; | |
281 | Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z | |
282 | ||
283 | Float_t delta1PtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC; | |
284 | Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC; | |
285 | ||
286 | Float_t mceta = infoMC->GetParticle().Eta(); | |
287 | Float_t mcphi = infoMC->GetParticle().Phi(); | |
288 | if(mcphi<0) mcphi += 2.*TMath::Pi(); | |
289 | Float_t mcpt = infoMC->GetParticle().Pt(); | |
290 | ||
291 | // distance to Prim. vertex | |
292 | const Double_t* dv = infoMC->GetVDist(); | |
293 | Bool_t isPrim = TMath::Sqrt(dv[0]*dv[0] + dv[1]*dv[1])<fCutsMC->GetMaxR() && TMath::Abs(dv[2])<fCutsMC->GetMaxVz(); | |
294 | ||
295 | // Only 5 charged particle species (e,mu,pi,K,p) | |
296 | if (fCutsMC->IsPdgParticle(TMath::Abs(infoMC->GetParticle().GetPdgCode())) == kFALSE) return; | |
297 | if (!isPrim) return; | |
298 | if (infoRC->GetStatus(1)!=3) return; // TPC refit | |
299 | if (!infoRC->GetESDtrack()) return; | |
300 | if (infoRC->GetESDtrack()->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return; | |
301 | ||
302 | infoRC->GetESDtrack()->GetImpactParameters(dca,cov); | |
303 | ||
304 | // select primaries | |
305 | if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ()) | |
306 | { | |
307 | deltaYTPC= infoRC->GetESDtrack()->GetY()-infoMC->GetParticle().Vy(); | |
308 | deltaZTPC = infoRC->GetESDtrack()->GetZ()-infoMC->GetParticle().Vz(); | |
309 | deltaLambdaTPC = TMath::ATan2(infoRC->GetESDtrack()->Pz(),infoRC->GetESDtrack()->Pt())-TMath::ATan2(infoMC->GetParticle().Pz(),infoMC->GetParticle().Pt()); | |
310 | deltaPhiTPC = TMath::ATan2(infoRC->GetESDtrack()->Py(),infoRC->GetESDtrack()->Px())-TMath::ATan2(infoMC->GetParticle().Py(),infoMC->GetParticle().Px()); | |
311 | delta1PtTPC = (infoRC->GetESDtrack()->OneOverPt()-1./mcpt)*mcpt; | |
312 | ||
313 | pullYTPC= (infoRC->GetESDtrack()->GetY()-infoMC->GetParticle().Vy()) / TMath::Sqrt(infoRC->GetESDtrack()->GetSigmaY2()); | |
314 | pullZTPC = (infoRC->GetESDtrack()->GetZ()-infoMC->GetParticle().Vz()) / TMath::Sqrt(infoRC->GetESDtrack()->GetSigmaZ2()); | |
315 | pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(infoRC->GetESDtrack()->GetSigmaTgl2()); | |
316 | pullPhiTPC = deltaPhiTPC / TMath::Sqrt(infoRC->GetESDtrack()->GetSigmaSnp2()); | |
317 | pull1PtTPC = (infoRC->GetESDtrack()->OneOverPt()-1./mcpt) / TMath::Sqrt(infoRC->GetESDtrack()->GetSigma1Pt2()); | |
318 | ||
319 | Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,delta1PtTPC,infoMC->GetParticle().Vy(),infoMC->GetParticle().Vz(),mceta,mcphi,mcpt}; | |
320 | Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,infoMC->GetParticle().Vy(),infoMC->GetParticle().Vz(),mceta,mcphi,mcpt}; | |
321 | ||
322 | // TPC and ITS clusters in the system | |
323 | if(infoRC->GetESDtrack()->GetITSclusters(clusterITS)>fCutsRC->GetMinNClustersITS()) | |
324 | { | |
325 | fResolHisto->Fill(vResolHisto); | |
326 | fPullHisto->Fill(vPullHisto); | |
327 | } | |
35771050 | 328 | } |
96b3109a | 329 | } |
330 | ||
331 | //_____________________________________________________________________________ | |
71a14197 | 332 | void AliComparisonRes::ProcessConstrained(AliMCInfo* const infoMC, AliESDRecInfo *const infoRC) |
96b3109a | 333 | { |
334 | // Fill resolution comparison information (constarained parameters) | |
335 | // | |
336 | AliExternalTrackParam *track = 0; | |
ff0f1307 | 337 | //Double_t field = AliTracker::GetBz(); // nominal Bz field [kG] |
96b3109a | 338 | Double_t kMaxD = 123456.0; // max distance |
35771050 | 339 | AliESDVertex vertexMC; // MC primary vertex |
96b3109a | 340 | |
96b3109a | 341 | Double_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z |
96b3109a | 342 | |
71a14197 | 343 | Float_t delta1PtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC; |
344 | Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC; | |
345 | ||
346 | Float_t mceta = infoMC->GetParticle().Eta(); | |
347 | Float_t mcphi = infoMC->GetParticle().Phi(); | |
348 | if(mcphi<0) mcphi += 2.*TMath::Pi(); | |
96b3109a | 349 | Float_t mcpt = infoMC->GetParticle().Pt(); |
96b3109a | 350 | |
351 | // distance to Prim. vertex | |
352 | const Double_t* dv = infoMC->GetVDist(); | |
353 | Bool_t isPrim = TMath::Sqrt(dv[0]*dv[0] + dv[1]*dv[1])<fCutsMC->GetMaxR() && TMath::Abs(dv[2])<fCutsMC->GetMaxVz(); | |
71a14197 | 354 | |
355 | // calculate and set prim. vertex | |
356 | vertexMC.SetXv( infoMC->GetParticle().Vx() - dv[0] ); | |
357 | vertexMC.SetYv( infoMC->GetParticle().Vy() - dv[1] ); | |
358 | vertexMC.SetZv( infoMC->GetParticle().Vz() - dv[2] ); | |
359 | ||
96b3109a | 360 | // Check selection cuts |
361 | if (fCutsMC->IsPdgParticle(TMath::Abs(infoMC->GetParticle().GetPdgCode())) == kFALSE) return; | |
362 | if (!isPrim) return; | |
71a14197 | 363 | if (infoRC->GetStatus(1)!=3) return; // TPC refit |
96b3109a | 364 | if (!infoRC->GetESDtrack()) return; |
365 | if (infoRC->GetESDtrack()->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return; | |
96b3109a | 366 | |
367 | // constrained parameters resolution | |
368 | const AliExternalTrackParam * cparam = infoRC->GetESDtrack()->GetConstrainedParam(); | |
71a14197 | 369 | if(!cparam) return; |
96b3109a | 370 | |
71a14197 | 371 | if ((track = new AliExternalTrackParam(*cparam)) != 0) |
96b3109a | 372 | { |
ff0f1307 | 373 | //Bool_t bDCAStatus = track->PropagateToDCA(&vertexMC,field,kMaxD,dca,cov); |
374 | Double_t field[3]; track->GetBxByBz(field); | |
375 | Bool_t bDCAStatus = track->PropagateToDCABxByBz(&vertexMC,field,kMaxD,dca,cov); | |
71a14197 | 376 | if(bDCAStatus) { |
377 | if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ()) | |
378 | { | |
379 | deltaYTPC= track->GetY()-infoMC->GetParticle().Vy(); | |
380 | deltaZTPC = track->GetZ()-infoMC->GetParticle().Vz(); | |
381 | deltaLambdaTPC = TMath::ATan2(track->Pz(),track->Pt())-TMath::ATan2(infoMC->GetParticle().Pz(),infoMC->GetParticle().Pt()); | |
382 | deltaPhiTPC = TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(infoMC->GetParticle().Py(),infoMC->GetParticle().Px()); | |
383 | delta1PtTPC = (track->OneOverPt()-1./mcpt)*mcpt; | |
384 | ||
385 | pullYTPC= (track->GetY()-infoMC->GetParticle().Vy()) / TMath::Sqrt(track->GetSigmaY2()); | |
386 | pullZTPC = (track->GetZ()-infoMC->GetParticle().Vz()) / TMath::Sqrt(track->GetSigmaZ2()); | |
387 | pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(track->GetSigmaTgl2()); | |
388 | pullPhiTPC = deltaPhiTPC / TMath::Sqrt(track->GetSigmaSnp2()); | |
389 | pull1PtTPC = (track->OneOverPt()-1./mcpt) / TMath::Sqrt(track->GetSigma1Pt2()); | |
390 | ||
391 | Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,delta1PtTPC,infoMC->GetParticle().Vy(),infoMC->GetParticle().Vz(),mceta,mcphi,mcpt}; | |
392 | fResolHisto->Fill(vResolHisto); | |
393 | ||
394 | Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,infoMC->GetParticle().Vy(),infoMC->GetParticle().Vz(),mceta,mcphi,mcpt}; | |
395 | fPullHisto->Fill(vPullHisto); | |
396 | } | |
96b3109a | 397 | } |
71a14197 | 398 | delete track; |
35771050 | 399 | } |
96b3109a | 400 | } |
401 | ||
402 | //_____________________________________________________________________________ | |
71a14197 | 403 | void AliComparisonRes::Exec(AliMCInfo* const infoMC, AliESDRecInfo* const infoRC){ |
96b3109a | 404 | |
405 | // Process comparison information | |
71a14197 | 406 | if(GetAnalysisMode() == 0) ProcessTPC(infoMC,infoRC); |
407 | else if(GetAnalysisMode() == 1) ProcessTPCITS(infoMC,infoRC); | |
408 | else if(GetAnalysisMode() == 2) ProcessConstrained(infoMC,infoRC); | |
409 | else { | |
410 | printf("ERROR: AnalysisMode %d \n",fAnalysisMode); | |
411 | return; | |
412 | } | |
96b3109a | 413 | } |
414 | ||
415 | //_____________________________________________________________________________ | |
416 | TH1F* AliComparisonRes::MakeResol(TH2F * his, Int_t integ, Bool_t type){ | |
417 | // Create resolution histograms | |
418 | ||
419 | TH1F *hisr, *hism; | |
420 | if (!gPad) new TCanvas; | |
421 | hisr = AliTreeDraw::CreateResHistoI(his,&hism,integ); | |
422 | if (type) return hism; | |
423 | else | |
424 | return hisr; | |
425 | } | |
426 | ||
427 | //_____________________________________________________________________________ | |
428 | void AliComparisonRes::Analyse(){ | |
3baa4bfd | 429 | // Analyse comparison information and store output histograms |
430 | // in the folder "folderRes" | |
431 | // | |
3baa4bfd | 432 | TH1::AddDirectory(kFALSE); |
71a14197 | 433 | TH1F *h=0; |
434 | TH2F *h2D=0; | |
b4126c69 | 435 | TObjArray *aFolderObj = new TObjArray; |
3baa4bfd | 436 | |
437 | // write results in the folder | |
96b3109a | 438 | TCanvas * c = new TCanvas("Phi resol Tan","Phi resol Tan"); |
439 | c->cd(); | |
3baa4bfd | 440 | |
71a14197 | 441 | char name[256]; |
442 | char res_title[256] = {"res_y:res_z:res_lambda:res_phi:res_1pt:y:z:eta:phi:pt"} ; | |
443 | char pull_title[256] = {"pull_y:pull_z:pull_lambda:pull_phi:pull_1pt:y:z:eta:phi:pt"}; | |
35771050 | 444 | |
71a14197 | 445 | for(Int_t i=0; i<5; i++) |
446 | { | |
447 | for(Int_t j=5; j<10; j++) | |
448 | { | |
449 | if(j==7) fResolHisto->GetAxis(9)->SetRangeUser(0.2,10.); // pt threshold | |
450 | if(j==9) fResolHisto->GetAxis(7)->SetRangeUser(-0.9,0.9); // eta window | |
96b3109a | 451 | |
71a14197 | 452 | h2D = (TH2F*)fResolHisto->Projection(i,j); |
453 | h = AliComparisonRes::MakeResol(h2D,1,0); | |
454 | sprintf(name,"h_res_%d_vs_%d",i,j); | |
455 | h->SetName(name); | |
456 | h->SetTitle(res_title); | |
96b3109a | 457 | |
71a14197 | 458 | aFolderObj->Add(h); |
96b3109a | 459 | |
71a14197 | 460 | h = AliComparisonRes::MakeResol(h2D,1,1); |
461 | sprintf(name,"h_mean_res_%d_vs_%d",i,j); | |
462 | h->SetName(name); | |
463 | h->SetTitle(res_title); | |
96b3109a | 464 | |
71a14197 | 465 | aFolderObj->Add(h); |
96b3109a | 466 | |
71a14197 | 467 | if(j==7) fResolHisto->GetAxis(9)->SetRangeUser(0.0,10.); |
468 | if(j==9) fResolHisto->GetAxis(7)->SetRangeUser(-1.5,1.5); | |
96b3109a | 469 | |
71a14197 | 470 | // |
471 | if(j==7) fPullHisto->GetAxis(9)->SetRangeUser(0.2,10.); | |
472 | if(j==9) fPullHisto->GetAxis(7)->SetRangeUser(-0.9,0.9); | |
96b3109a | 473 | |
71a14197 | 474 | h2D = (TH2F*)fPullHisto->Projection(i,j); |
475 | h = AliComparisonRes::MakeResol(h2D,1,0); | |
476 | sprintf(name,"h_pull_%d_vs_%d",i,j); | |
477 | h->SetName(name); | |
478 | h->SetTitle(pull_title); | |
96b3109a | 479 | |
71a14197 | 480 | aFolderObj->Add(h); |
96b3109a | 481 | |
71a14197 | 482 | h = AliComparisonRes::MakeResol(h2D,1,1); |
483 | sprintf(name,"h_mean_pull_%d_vs_%d",i,j); | |
484 | h->SetName(name); | |
485 | h->SetTitle(pull_title); | |
96b3109a | 486 | |
71a14197 | 487 | aFolderObj->Add(h); |
96b3109a | 488 | |
71a14197 | 489 | if(j==7) fPullHisto->GetAxis(9)->SetRangeUser(0.0,10.); |
490 | if(j==9) fPullHisto->GetAxis(7)->SetRangeUser(-1.5,1.5); | |
491 | } | |
492 | } | |
b4126c69 | 493 | |
494 | // export objects to analysis folder | |
495 | fAnalysisFolder = ExportToFolder(aFolderObj); | |
496 | ||
497 | // delete only TObjArray | |
498 | if(aFolderObj) delete aFolderObj; | |
499 | } | |
500 | ||
501 | //_____________________________________________________________________________ | |
502 | TFolder* AliComparisonRes::ExportToFolder(TObjArray * array) | |
503 | { | |
504 | // recreate folder avery time and export objects to new one | |
505 | // | |
506 | AliComparisonRes * comp=this; | |
507 | TFolder *folder = comp->GetAnalysisFolder(); | |
508 | ||
509 | TString name, title; | |
510 | TFolder *newFolder = 0; | |
511 | Int_t i = 0; | |
512 | Int_t size = array->GetSize(); | |
513 | ||
514 | if(folder) { | |
515 | // get name and title from old folder | |
516 | name = folder->GetName(); | |
517 | title = folder->GetTitle(); | |
518 | ||
519 | // delete old one | |
520 | delete folder; | |
521 | ||
522 | // create new one | |
523 | newFolder = CreateFolder(name.Data(),title.Data()); | |
524 | newFolder->SetOwner(); | |
525 | ||
526 | // add objects to folder | |
527 | while(i < size) { | |
528 | newFolder->Add(array->At(i)); | |
529 | i++; | |
530 | } | |
531 | } | |
96b3109a | 532 | |
b4126c69 | 533 | return newFolder; |
96b3109a | 534 | } |
535 | ||
536 | //_____________________________________________________________________________ | |
71a14197 | 537 | Long64_t AliComparisonRes::Merge(TCollection* const list) |
96b3109a | 538 | { |
539 | // Merge list of objects (needed by PROOF) | |
540 | ||
541 | if (!list) | |
542 | return 0; | |
543 | ||
544 | if (list->IsEmpty()) | |
545 | return 1; | |
546 | ||
547 | TIterator* iter = list->MakeIterator(); | |
548 | TObject* obj = 0; | |
549 | ||
550 | // collection of generated histograms | |
551 | Int_t count=0; | |
552 | while((obj = iter->Next()) != 0) | |
553 | { | |
554 | AliComparisonRes* entry = dynamic_cast<AliComparisonRes*>(obj); | |
555 | if (entry == 0) continue; | |
71a14197 | 556 | |
557 | fResolHisto->Add(entry->fResolHisto); | |
558 | fPullHisto->Add(entry->fPullHisto); | |
96b3109a | 559 | |
560 | count++; | |
561 | } | |
562 | ||
563 | return count; | |
564 | } | |
3baa4bfd | 565 | |
566 | //_____________________________________________________________________________ | |
567 | TFolder* AliComparisonRes::CreateFolder(TString name,TString title) { | |
568 | // create folder for analysed histograms | |
569 | // | |
570 | TFolder *folder = 0; | |
571 | folder = new TFolder(name.Data(),title.Data()); | |
572 | ||
573 | return folder; | |
574 | } |