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