]>
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; | |
96b3109a | 182 | Double_t field = AliTracker::GetBz(); // nominal Bz field [kG] |
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 | { | |
35771050 | 219 | Bool_t bDCAStatus = track->PropagateToDCA(&vertexMC,field,kMaxD,dca,cov); |
96b3109a | 220 | |
221 | // Fill parametrisation histograms (only TPC track) | |
3baa4bfd | 222 | if(bDCAStatus) |
71a14197 | 223 | { |
224 | // select primaries | |
225 | if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ()) | |
226 | { | |
227 | ||
228 | deltaYTPC= track->GetY()-infoMC->GetParticle().Vy(); | |
229 | deltaZTPC = track->GetZ()-infoMC->GetParticle().Vz(); | |
230 | deltaLambdaTPC = TMath::ATan2(track->Pz(),track->Pt())-TMath::ATan2(infoMC->GetParticle().Pz(),infoMC->GetParticle().Pt()); | |
231 | deltaPhiTPC = TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(infoMC->GetParticle().Py(),infoMC->GetParticle().Px()); | |
232 | delta1PtTPC = (track->OneOverPt()-1./mcpt)*mcpt; | |
233 | ||
234 | pullYTPC= (track->GetY()-infoMC->GetParticle().Vy()) / TMath::Sqrt(track->GetSigmaY2()); | |
235 | pullZTPC = (track->GetZ()-infoMC->GetParticle().Vz()) / TMath::Sqrt(track->GetSigmaZ2()); | |
236 | pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(track->GetSigmaTgl2()); | |
237 | pullPhiTPC = deltaPhiTPC / TMath::Sqrt(track->GetSigmaSnp2()); | |
238 | pull1PtTPC = (track->OneOverPt()-1./mcpt) / TMath::Sqrt(track->GetSigma1Pt2()); | |
239 | ||
240 | Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,delta1PtTPC,infoMC->GetParticle().Vy(),infoMC->GetParticle().Vz(),mceta,mcphi,mcpt}; | |
241 | fResolHisto->Fill(vResolHisto); | |
242 | ||
243 | Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,infoMC->GetParticle().Vy(),infoMC->GetParticle().Vz(),mceta,mcphi,mcpt}; | |
244 | fPullHisto->Fill(vPullHisto); | |
96b3109a | 245 | } |
71a14197 | 246 | |
247 | /* | |
248 | delta1Pt2TPC = (1/mcpt-innerTPC->OneOverPt())/TMath::Power(1+1/mcpt,2); | |
249 | deltaY1PtTPC= (infoMC->GetParticle().Vy()-innerTPC->GetY()) / (0.2+1/mcpt); | |
250 | deltaZ1PtTPC = (infoMC->GetParticle().Vz()-innerTPC->GetZ()) / (0.2+1/mcpt); | |
251 | deltaPhi1PtTPC = deltaPhiTPC / (0.1+1/mcpt); | |
252 | deltaTheta1PtTPC = deltaLambdaTPC / (0.1+1/mcpt); | |
253 | ||
254 | fPhiEtaPtTPC->Fill(TMath::ATan2(innerTPC->Py(),innerTPC->Px()), innerTPC->Eta(), innerTPC->Pt()); | |
255 | ||
256 | f1Pt2ResolS1PtTPC->Fill(s1mcpt,delta1Pt2TPC); | |
257 | fYResolS1PtTPC->Fill(s1mcpt,deltaY1PtTPC); | |
258 | fZResolS1PtTPC->Fill(s1mcpt,deltaZ1PtTPC); | |
259 | fPhiResolS1PtTPC->Fill(s1mcpt,deltaPhi1PtTPC); | |
260 | fThetaResolS1PtTPC->Fill(s1mcpt,deltaTheta1PtTPC); | |
261 | ||
262 | fPtResolPtTPC->Fill(mcpt,deltaPtTPC); | |
263 | fPtPullPtTPC->Fill(mcpt,pullPtTPC); | |
264 | fPhiResolEtaTPC->Fill(eta,deltaPhiTPC); | |
265 | fPhiPullEtaTPC->Fill(eta,pullPhiTPC); | |
266 | fLambdaResolEtaTPC->Fill(eta,deltaLambdaTPC); | |
267 | fLambdaPullEtaTPC->Fill(eta,pullLambdaTPC); | |
268 | */ | |
269 | } | |
270 | delete track; | |
96b3109a | 271 | } |
272 | } | |
71a14197 | 273 | } |
96b3109a | 274 | |
71a14197 | 275 | //_____________________________________________________________________________ |
276 | void AliComparisonRes::ProcessTPCITS(AliMCInfo* const infoMC, AliESDRecInfo *const infoRC) | |
277 | { | |
278 | Int_t clusterITS[200]; | |
279 | Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z | |
280 | ||
281 | Float_t delta1PtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC; | |
282 | Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC; | |
283 | ||
284 | Float_t mceta = infoMC->GetParticle().Eta(); | |
285 | Float_t mcphi = infoMC->GetParticle().Phi(); | |
286 | if(mcphi<0) mcphi += 2.*TMath::Pi(); | |
287 | Float_t mcpt = infoMC->GetParticle().Pt(); | |
288 | ||
289 | // distance to Prim. vertex | |
290 | const Double_t* dv = infoMC->GetVDist(); | |
291 | Bool_t isPrim = TMath::Sqrt(dv[0]*dv[0] + dv[1]*dv[1])<fCutsMC->GetMaxR() && TMath::Abs(dv[2])<fCutsMC->GetMaxVz(); | |
292 | ||
293 | // Only 5 charged particle species (e,mu,pi,K,p) | |
294 | if (fCutsMC->IsPdgParticle(TMath::Abs(infoMC->GetParticle().GetPdgCode())) == kFALSE) return; | |
295 | if (!isPrim) return; | |
296 | if (infoRC->GetStatus(1)!=3) return; // TPC refit | |
297 | if (!infoRC->GetESDtrack()) return; | |
298 | if (infoRC->GetESDtrack()->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return; | |
299 | ||
300 | infoRC->GetESDtrack()->GetImpactParameters(dca,cov); | |
301 | ||
302 | // select primaries | |
303 | if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ()) | |
304 | { | |
305 | deltaYTPC= infoRC->GetESDtrack()->GetY()-infoMC->GetParticle().Vy(); | |
306 | deltaZTPC = infoRC->GetESDtrack()->GetZ()-infoMC->GetParticle().Vz(); | |
307 | deltaLambdaTPC = TMath::ATan2(infoRC->GetESDtrack()->Pz(),infoRC->GetESDtrack()->Pt())-TMath::ATan2(infoMC->GetParticle().Pz(),infoMC->GetParticle().Pt()); | |
308 | deltaPhiTPC = TMath::ATan2(infoRC->GetESDtrack()->Py(),infoRC->GetESDtrack()->Px())-TMath::ATan2(infoMC->GetParticle().Py(),infoMC->GetParticle().Px()); | |
309 | delta1PtTPC = (infoRC->GetESDtrack()->OneOverPt()-1./mcpt)*mcpt; | |
310 | ||
311 | pullYTPC= (infoRC->GetESDtrack()->GetY()-infoMC->GetParticle().Vy()) / TMath::Sqrt(infoRC->GetESDtrack()->GetSigmaY2()); | |
312 | pullZTPC = (infoRC->GetESDtrack()->GetZ()-infoMC->GetParticle().Vz()) / TMath::Sqrt(infoRC->GetESDtrack()->GetSigmaZ2()); | |
313 | pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(infoRC->GetESDtrack()->GetSigmaTgl2()); | |
314 | pullPhiTPC = deltaPhiTPC / TMath::Sqrt(infoRC->GetESDtrack()->GetSigmaSnp2()); | |
315 | pull1PtTPC = (infoRC->GetESDtrack()->OneOverPt()-1./mcpt) / TMath::Sqrt(infoRC->GetESDtrack()->GetSigma1Pt2()); | |
316 | ||
317 | Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,delta1PtTPC,infoMC->GetParticle().Vy(),infoMC->GetParticle().Vz(),mceta,mcphi,mcpt}; | |
318 | Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,infoMC->GetParticle().Vy(),infoMC->GetParticle().Vz(),mceta,mcphi,mcpt}; | |
319 | ||
320 | // TPC and ITS clusters in the system | |
321 | if(infoRC->GetESDtrack()->GetITSclusters(clusterITS)>fCutsRC->GetMinNClustersITS()) | |
322 | { | |
323 | fResolHisto->Fill(vResolHisto); | |
324 | fPullHisto->Fill(vPullHisto); | |
325 | } | |
35771050 | 326 | } |
96b3109a | 327 | } |
328 | ||
329 | //_____________________________________________________________________________ | |
71a14197 | 330 | void AliComparisonRes::ProcessConstrained(AliMCInfo* const infoMC, AliESDRecInfo *const infoRC) |
96b3109a | 331 | { |
332 | // Fill resolution comparison information (constarained parameters) | |
333 | // | |
334 | AliExternalTrackParam *track = 0; | |
96b3109a | 335 | Double_t field = AliTracker::GetBz(); // nominal Bz field [kG] |
336 | Double_t kMaxD = 123456.0; // max distance | |
35771050 | 337 | AliESDVertex vertexMC; // MC primary vertex |
96b3109a | 338 | |
96b3109a | 339 | Double_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z |
96b3109a | 340 | |
71a14197 | 341 | Float_t delta1PtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC; |
342 | Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC; | |
343 | ||
344 | Float_t mceta = infoMC->GetParticle().Eta(); | |
345 | Float_t mcphi = infoMC->GetParticle().Phi(); | |
346 | if(mcphi<0) mcphi += 2.*TMath::Pi(); | |
96b3109a | 347 | Float_t mcpt = infoMC->GetParticle().Pt(); |
96b3109a | 348 | |
349 | // distance to Prim. vertex | |
350 | const Double_t* dv = infoMC->GetVDist(); | |
351 | Bool_t isPrim = TMath::Sqrt(dv[0]*dv[0] + dv[1]*dv[1])<fCutsMC->GetMaxR() && TMath::Abs(dv[2])<fCutsMC->GetMaxVz(); | |
71a14197 | 352 | |
353 | // calculate and set prim. vertex | |
354 | vertexMC.SetXv( infoMC->GetParticle().Vx() - dv[0] ); | |
355 | vertexMC.SetYv( infoMC->GetParticle().Vy() - dv[1] ); | |
356 | vertexMC.SetZv( infoMC->GetParticle().Vz() - dv[2] ); | |
357 | ||
96b3109a | 358 | // Check selection cuts |
359 | if (fCutsMC->IsPdgParticle(TMath::Abs(infoMC->GetParticle().GetPdgCode())) == kFALSE) return; | |
360 | if (!isPrim) return; | |
71a14197 | 361 | if (infoRC->GetStatus(1)!=3) return; // TPC refit |
96b3109a | 362 | if (!infoRC->GetESDtrack()) return; |
363 | if (infoRC->GetESDtrack()->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return; | |
96b3109a | 364 | |
365 | // constrained parameters resolution | |
366 | const AliExternalTrackParam * cparam = infoRC->GetESDtrack()->GetConstrainedParam(); | |
71a14197 | 367 | if(!cparam) return; |
96b3109a | 368 | |
71a14197 | 369 | if ((track = new AliExternalTrackParam(*cparam)) != 0) |
96b3109a | 370 | { |
71a14197 | 371 | Bool_t bDCAStatus = track->PropagateToDCA(&vertexMC,field,kMaxD,dca,cov); |
372 | if(bDCAStatus) { | |
373 | if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ()) | |
374 | { | |
375 | deltaYTPC= track->GetY()-infoMC->GetParticle().Vy(); | |
376 | deltaZTPC = track->GetZ()-infoMC->GetParticle().Vz(); | |
377 | deltaLambdaTPC = TMath::ATan2(track->Pz(),track->Pt())-TMath::ATan2(infoMC->GetParticle().Pz(),infoMC->GetParticle().Pt()); | |
378 | deltaPhiTPC = TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(infoMC->GetParticle().Py(),infoMC->GetParticle().Px()); | |
379 | delta1PtTPC = (track->OneOverPt()-1./mcpt)*mcpt; | |
380 | ||
381 | pullYTPC= (track->GetY()-infoMC->GetParticle().Vy()) / TMath::Sqrt(track->GetSigmaY2()); | |
382 | pullZTPC = (track->GetZ()-infoMC->GetParticle().Vz()) / TMath::Sqrt(track->GetSigmaZ2()); | |
383 | pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(track->GetSigmaTgl2()); | |
384 | pullPhiTPC = deltaPhiTPC / TMath::Sqrt(track->GetSigmaSnp2()); | |
385 | pull1PtTPC = (track->OneOverPt()-1./mcpt) / TMath::Sqrt(track->GetSigma1Pt2()); | |
386 | ||
387 | Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,delta1PtTPC,infoMC->GetParticle().Vy(),infoMC->GetParticle().Vz(),mceta,mcphi,mcpt}; | |
388 | fResolHisto->Fill(vResolHisto); | |
389 | ||
390 | Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,infoMC->GetParticle().Vy(),infoMC->GetParticle().Vz(),mceta,mcphi,mcpt}; | |
391 | fPullHisto->Fill(vPullHisto); | |
392 | } | |
96b3109a | 393 | } |
71a14197 | 394 | delete track; |
35771050 | 395 | } |
96b3109a | 396 | } |
397 | ||
398 | //_____________________________________________________________________________ | |
71a14197 | 399 | void AliComparisonRes::Exec(AliMCInfo* const infoMC, AliESDRecInfo* const infoRC){ |
96b3109a | 400 | |
401 | // Process comparison information | |
71a14197 | 402 | if(GetAnalysisMode() == 0) ProcessTPC(infoMC,infoRC); |
403 | else if(GetAnalysisMode() == 1) ProcessTPCITS(infoMC,infoRC); | |
404 | else if(GetAnalysisMode() == 2) ProcessConstrained(infoMC,infoRC); | |
405 | else { | |
406 | printf("ERROR: AnalysisMode %d \n",fAnalysisMode); | |
407 | return; | |
408 | } | |
96b3109a | 409 | } |
410 | ||
411 | //_____________________________________________________________________________ | |
412 | TH1F* AliComparisonRes::MakeResol(TH2F * his, Int_t integ, Bool_t type){ | |
413 | // Create resolution histograms | |
414 | ||
415 | TH1F *hisr, *hism; | |
416 | if (!gPad) new TCanvas; | |
417 | hisr = AliTreeDraw::CreateResHistoI(his,&hism,integ); | |
418 | if (type) return hism; | |
419 | else | |
420 | return hisr; | |
421 | } | |
422 | ||
423 | //_____________________________________________________________________________ | |
424 | void AliComparisonRes::Analyse(){ | |
3baa4bfd | 425 | // Analyse comparison information and store output histograms |
426 | // in the folder "folderRes" | |
427 | // | |
3baa4bfd | 428 | TH1::AddDirectory(kFALSE); |
71a14197 | 429 | TH1F *h=0; |
430 | TH2F *h2D=0; | |
b4126c69 | 431 | TObjArray *aFolderObj = new TObjArray; |
3baa4bfd | 432 | |
433 | // write results in the folder | |
96b3109a | 434 | TCanvas * c = new TCanvas("Phi resol Tan","Phi resol Tan"); |
435 | c->cd(); | |
3baa4bfd | 436 | |
71a14197 | 437 | char name[256]; |
438 | char res_title[256] = {"res_y:res_z:res_lambda:res_phi:res_1pt:y:z:eta:phi:pt"} ; | |
439 | char pull_title[256] = {"pull_y:pull_z:pull_lambda:pull_phi:pull_1pt:y:z:eta:phi:pt"}; | |
35771050 | 440 | |
71a14197 | 441 | for(Int_t i=0; i<5; i++) |
442 | { | |
443 | for(Int_t j=5; j<10; j++) | |
444 | { | |
445 | if(j==7) fResolHisto->GetAxis(9)->SetRangeUser(0.2,10.); // pt threshold | |
446 | if(j==9) fResolHisto->GetAxis(7)->SetRangeUser(-0.9,0.9); // eta window | |
96b3109a | 447 | |
71a14197 | 448 | h2D = (TH2F*)fResolHisto->Projection(i,j); |
449 | h = AliComparisonRes::MakeResol(h2D,1,0); | |
450 | sprintf(name,"h_res_%d_vs_%d",i,j); | |
451 | h->SetName(name); | |
452 | h->SetTitle(res_title); | |
96b3109a | 453 | |
71a14197 | 454 | aFolderObj->Add(h); |
96b3109a | 455 | |
71a14197 | 456 | h = AliComparisonRes::MakeResol(h2D,1,1); |
457 | sprintf(name,"h_mean_res_%d_vs_%d",i,j); | |
458 | h->SetName(name); | |
459 | h->SetTitle(res_title); | |
96b3109a | 460 | |
71a14197 | 461 | aFolderObj->Add(h); |
96b3109a | 462 | |
71a14197 | 463 | if(j==7) fResolHisto->GetAxis(9)->SetRangeUser(0.0,10.); |
464 | if(j==9) fResolHisto->GetAxis(7)->SetRangeUser(-1.5,1.5); | |
96b3109a | 465 | |
71a14197 | 466 | // |
467 | if(j==7) fPullHisto->GetAxis(9)->SetRangeUser(0.2,10.); | |
468 | if(j==9) fPullHisto->GetAxis(7)->SetRangeUser(-0.9,0.9); | |
96b3109a | 469 | |
71a14197 | 470 | h2D = (TH2F*)fPullHisto->Projection(i,j); |
471 | h = AliComparisonRes::MakeResol(h2D,1,0); | |
472 | sprintf(name,"h_pull_%d_vs_%d",i,j); | |
473 | h->SetName(name); | |
474 | h->SetTitle(pull_title); | |
96b3109a | 475 | |
71a14197 | 476 | aFolderObj->Add(h); |
96b3109a | 477 | |
71a14197 | 478 | h = AliComparisonRes::MakeResol(h2D,1,1); |
479 | sprintf(name,"h_mean_pull_%d_vs_%d",i,j); | |
480 | h->SetName(name); | |
481 | h->SetTitle(pull_title); | |
96b3109a | 482 | |
71a14197 | 483 | aFolderObj->Add(h); |
96b3109a | 484 | |
71a14197 | 485 | if(j==7) fPullHisto->GetAxis(9)->SetRangeUser(0.0,10.); |
486 | if(j==9) fPullHisto->GetAxis(7)->SetRangeUser(-1.5,1.5); | |
487 | } | |
488 | } | |
b4126c69 | 489 | |
490 | // export objects to analysis folder | |
491 | fAnalysisFolder = ExportToFolder(aFolderObj); | |
492 | ||
493 | // delete only TObjArray | |
494 | if(aFolderObj) delete aFolderObj; | |
495 | } | |
496 | ||
497 | //_____________________________________________________________________________ | |
498 | TFolder* AliComparisonRes::ExportToFolder(TObjArray * array) | |
499 | { | |
500 | // recreate folder avery time and export objects to new one | |
501 | // | |
502 | AliComparisonRes * comp=this; | |
503 | TFolder *folder = comp->GetAnalysisFolder(); | |
504 | ||
505 | TString name, title; | |
506 | TFolder *newFolder = 0; | |
507 | Int_t i = 0; | |
508 | Int_t size = array->GetSize(); | |
509 | ||
510 | if(folder) { | |
511 | // get name and title from old folder | |
512 | name = folder->GetName(); | |
513 | title = folder->GetTitle(); | |
514 | ||
515 | // delete old one | |
516 | delete folder; | |
517 | ||
518 | // create new one | |
519 | newFolder = CreateFolder(name.Data(),title.Data()); | |
520 | newFolder->SetOwner(); | |
521 | ||
522 | // add objects to folder | |
523 | while(i < size) { | |
524 | newFolder->Add(array->At(i)); | |
525 | i++; | |
526 | } | |
527 | } | |
96b3109a | 528 | |
b4126c69 | 529 | return newFolder; |
96b3109a | 530 | } |
531 | ||
532 | //_____________________________________________________________________________ | |
71a14197 | 533 | Long64_t AliComparisonRes::Merge(TCollection* const list) |
96b3109a | 534 | { |
535 | // Merge list of objects (needed by PROOF) | |
536 | ||
537 | if (!list) | |
538 | return 0; | |
539 | ||
540 | if (list->IsEmpty()) | |
541 | return 1; | |
542 | ||
543 | TIterator* iter = list->MakeIterator(); | |
544 | TObject* obj = 0; | |
545 | ||
546 | // collection of generated histograms | |
547 | Int_t count=0; | |
548 | while((obj = iter->Next()) != 0) | |
549 | { | |
550 | AliComparisonRes* entry = dynamic_cast<AliComparisonRes*>(obj); | |
551 | if (entry == 0) continue; | |
71a14197 | 552 | |
553 | fResolHisto->Add(entry->fResolHisto); | |
554 | fPullHisto->Add(entry->fPullHisto); | |
96b3109a | 555 | |
556 | count++; | |
557 | } | |
558 | ||
559 | return count; | |
560 | } | |
3baa4bfd | 561 | |
562 | //_____________________________________________________________________________ | |
563 | TFolder* AliComparisonRes::CreateFolder(TString name,TString title) { | |
564 | // create folder for analysed histograms | |
565 | // | |
566 | TFolder *folder = 0; | |
567 | folder = new TFolder(name.Data(),title.Data()); | |
568 | ||
569 | return folder; | |
570 | } |