]>
Commit | Line | Data |
---|---|---|
1ee39b3a | 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-commercialf 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 | /* $Id: AliTRDresolution.cxx 27496 2008-07-22 08:35:45Z cblume $ */ | |
17 | ||
18 | //////////////////////////////////////////////////////////////////////////// | |
19 | // // | |
20 | // TRD tracking resolution // | |
21 | // | |
22 | // The class performs resolution and residual studies | |
23 | // of the TRD tracks for the following quantities : | |
24 | // - spatial position (y, [z]) | |
25 | // - angular (phi) tracklet | |
26 | // - momentum at the track level | |
27 | // | |
28 | // The class has to be used for regular detector performance checks using the official macros: | |
29 | // - $ALICE_ROOT/TRD/qaRec/run.C | |
30 | // - $ALICE_ROOT/TRD/qaRec/makeResults.C | |
31 | // | |
32 | // For stand alone usage please refer to the following example: | |
33 | // { | |
34 | // gSystem->Load("libANALYSIS.so"); | |
35 | // gSystem->Load("libTRDqaRec.so"); | |
36 | // AliTRDresolution *res = new AliTRDresolution(); | |
37 | // //res->SetMCdata(); | |
38 | // //res->SetVerbose(); | |
39 | // //res->SetVisual(); | |
40 | // res->Load(); | |
41 | // if(!res->PostProcess()) return; | |
42 | // res->GetRefFigure(0); | |
43 | // } | |
44 | // | |
45 | // Authors: // | |
46 | // Alexandru Bercuci <A.Bercuci@gsi.de> // | |
47 | // Markus Fasel <M.Fasel@gsi.de> // | |
48 | // // | |
49 | //////////////////////////////////////////////////////////////////////////// | |
50 | ||
92c40e64 | 51 | #include <TSystem.h> |
3ceb45ae | 52 | #include <TStyle.h> |
1ee39b3a | 53 | #include <TROOT.h> |
54 | #include <TObjArray.h> | |
55 | #include <TH3.h> | |
56 | #include <TH2.h> | |
57 | #include <TH1.h> | |
3ceb45ae | 58 | #include <THnSparse.h> |
1ee39b3a | 59 | #include <TF1.h> |
60 | #include <TCanvas.h> | |
61 | #include <TGaxis.h> | |
62 | #include <TBox.h> | |
92c40e64 | 63 | #include <TLegend.h> |
1ee39b3a | 64 | #include <TGraphErrors.h> |
65 | #include <TGraphAsymmErrors.h> | |
08bdd9d1 | 66 | #include <TLinearFitter.h> |
1ee39b3a | 67 | #include <TMath.h> |
68 | #include <TMatrixT.h> | |
69 | #include <TVectorT.h> | |
70 | #include <TTreeStream.h> | |
71 | #include <TGeoManager.h> | |
92c40e64 | 72 | #include <TDatabasePDG.h> |
1ee39b3a | 73 | |
74 | #include "AliPID.h" | |
92d6d80c | 75 | #include "AliLog.h" |
92c40e64 | 76 | #include "AliESDtrack.h" |
068e2c00 | 77 | #include "AliMathBase.h" |
08bdd9d1 | 78 | #include "AliTrackPointArray.h" |
1ee39b3a | 79 | |
80 | #include "AliTRDresolution.h" | |
81 | #include "AliTRDgeometry.h" | |
3ceb45ae | 82 | #include "AliTRDtransform.h" |
1ee39b3a | 83 | #include "AliTRDpadPlane.h" |
84 | #include "AliTRDcluster.h" | |
85 | #include "AliTRDseedV1.h" | |
86 | #include "AliTRDtrackV1.h" | |
87 | #include "AliTRDReconstructor.h" | |
88 | #include "AliTRDrecoParam.h" | |
92c40e64 | 89 | #include "AliTRDpidUtil.h" |
fd7ffd88 | 90 | #include "AliTRDinfoGen.h" |
4860c9db | 91 | #include "AliAnalysisManager.h" |
1ee39b3a | 92 | #include "info/AliTRDclusterInfo.h" |
566c3d46 | 93 | #include "info/AliTRDeventInfo.h" |
1ee39b3a | 94 | |
95 | ClassImp(AliTRDresolution) | |
3ed01fbe | 96 | //ClassImp(AliTRDresolution::AliTRDresolutionProjection) |
1ee39b3a | 97 | |
3ceb45ae | 98 | Int_t const AliTRDresolution::fgkNbins[kNdim] = { |
99 | Int_t(kNbunchCross)/*bc*/, | |
100 | 180/*phi*/, | |
101 | 50/*eta*/, | |
3ceb45ae | 102 | 50/*dy*/, |
3ed01fbe | 103 | 40/*dphi*/, |
3ceb45ae | 104 | 50/*dz*/, |
566c3d46 | 105 | 3/*chg*species*/, |
3ed01fbe | 106 | kNpt/*pt*/ |
3ceb45ae | 107 | }; //! no of bins/projection |
108 | Double_t const AliTRDresolution::fgkMin[kNdim] = { | |
566c3d46 | 109 | -1.5, |
3ceb45ae | 110 | -TMath::Pi(), |
111 | -1., | |
3ceb45ae | 112 | -1.5, |
3ed01fbe | 113 | -10., |
3ceb45ae | 114 | -2.5, |
566c3d46 | 115 | -1.5, |
3ed01fbe | 116 | -0.5 |
3ceb45ae | 117 | }; //! low limits for projections |
118 | Double_t const AliTRDresolution::fgkMax[kNdim] = { | |
566c3d46 | 119 | 1.5, |
3ceb45ae | 120 | TMath::Pi(), |
121 | 1., | |
3ceb45ae | 122 | 1.5, |
3ed01fbe | 123 | 10., |
3ceb45ae | 124 | 2.5, |
566c3d46 | 125 | 1.5, |
3ed01fbe | 126 | kNpt-0.5 |
3ceb45ae | 127 | }; //! high limits for projections |
128 | Char_t const *AliTRDresolution::fgkTitle[kNdim] = { | |
129 | "bunch cross", | |
130 | "#phi [rad]", | |
131 | "#eta", | |
3ceb45ae | 132 | "#Deltay [cm]", |
3ed01fbe | 133 | "#Delta#phi [deg]", |
3ceb45ae | 134 | "#Deltaz [cm]", |
3ed01fbe | 135 | "chg*spec*rc", |
136 | "bin_p_{t}" | |
3ceb45ae | 137 | }; //! title of projection |
138 | ||
3ceb45ae | 139 | Char_t const * AliTRDresolution::fgPerformanceName[kNclasses] = { |
140 | "Cluster2Track" | |
1ee39b3a | 141 | ,"Tracklet2Track" |
a310e49b | 142 | ,"Tracklet2TRDin" |
1ee39b3a | 143 | ,"Cluster2MC" |
144 | ,"Tracklet2MC" | |
92d6d80c | 145 | ,"TRDin2MC" |
1ee39b3a | 146 | ,"TRD2MC" |
f429b017 | 147 | // ,"Tracklet2TRDout" |
148 | // ,"TRDout2MC" | |
1ee39b3a | 149 | }; |
3ed01fbe | 150 | Float_t AliTRDresolution::fgPtBin[kNpt+1]; |
1ee39b3a | 151 | |
152 | //________________________________________________________ | |
153 | AliTRDresolution::AliTRDresolution() | |
f2e89a4c | 154 | :AliTRDrecoTask() |
1ee39b3a | 155 | ,fIdxPlot(0) |
92c40e64 | 156 | ,fIdxFrame(0) |
566c3d46 | 157 | ,fPtThreshold(.3) |
9dcc64cc | 158 | ,fDyRange(0.75) |
566c3d46 | 159 | ,fBCbinTOF(0) |
160 | ,fBCbinFill(0) | |
3ceb45ae | 161 | ,fProj(NULL) |
92c40e64 | 162 | ,fDBPDG(NULL) |
b91fdd71 | 163 | ,fCl(NULL) |
b91fdd71 | 164 | ,fMCcl(NULL) |
f8f46e4d | 165 | { |
166 | // | |
167 | // Default constructor | |
168 | // | |
f2e89a4c | 169 | SetNameTitle("TRDresolution", "TRD spatial and momentum resolution"); |
3ceb45ae | 170 | MakePtSegmentation(); |
f8f46e4d | 171 | } |
172 | ||
705f8b0a | 173 | //________________________________________________________ |
3ed01fbe | 174 | AliTRDresolution::AliTRDresolution(char* name, Bool_t xchange) |
f2e89a4c | 175 | :AliTRDrecoTask(name, "TRD spatial and momentum resolution") |
f8f46e4d | 176 | ,fIdxPlot(0) |
177 | ,fIdxFrame(0) | |
566c3d46 | 178 | ,fPtThreshold(.3) |
9dcc64cc | 179 | ,fDyRange(0.75) |
566c3d46 | 180 | ,fBCbinTOF(0) |
181 | ,fBCbinFill(0) | |
3ceb45ae | 182 | ,fProj(NULL) |
f8f46e4d | 183 | ,fDBPDG(NULL) |
f8f46e4d | 184 | ,fCl(NULL) |
f8f46e4d | 185 | ,fMCcl(NULL) |
1ee39b3a | 186 | { |
187 | // | |
188 | // Default constructor | |
189 | // | |
190 | ||
1ee39b3a | 191 | InitFunctorList(); |
3ceb45ae | 192 | MakePtSegmentation(); |
3ed01fbe | 193 | if(xchange){ |
194 | SetUseExchangeContainers(); | |
195 | DefineOutput(kClToTrk, TObjArray::Class()); // cluster2track | |
196 | DefineOutput(kClToMC, TObjArray::Class()); // cluster2mc | |
197 | } | |
1ee39b3a | 198 | } |
199 | ||
200 | //________________________________________________________ | |
201 | AliTRDresolution::~AliTRDresolution() | |
202 | { | |
203 | // | |
204 | // Destructor | |
205 | // | |
f0857a6a | 206 | if (AliAnalysisManager::GetAnalysisManager() && AliAnalysisManager::GetAnalysisManager()->IsProofMode()) return; |
3ceb45ae | 207 | if(fProj){fProj->Delete(); delete fProj;} |
1ee39b3a | 208 | if(fCl){fCl->Delete(); delete fCl;} |
1ee39b3a | 209 | if(fMCcl){fMCcl->Delete(); delete fMCcl;} |
1ee39b3a | 210 | } |
211 | ||
212 | ||
213 | //________________________________________________________ | |
f8f46e4d | 214 | void AliTRDresolution::UserCreateOutputObjects() |
1ee39b3a | 215 | { |
216 | // spatial resolution | |
5935a6da | 217 | |
068e2c00 | 218 | AliTRDrecoTask::UserCreateOutputObjects(); |
3ed01fbe | 219 | if(UseExchangeContainers()) InitExchangeContainers(); |
83b44483 | 220 | } |
1ee39b3a | 221 | |
83b44483 | 222 | //________________________________________________________ |
223 | void AliTRDresolution::InitExchangeContainers() | |
224 | { | |
61f6b45e | 225 | // Init containers for subsequent tasks (AliTRDclusterResolution) |
226 | ||
3ceb45ae | 227 | fCl = new TObjArray(200); fCl->SetOwner(kTRUE); |
228 | fMCcl = new TObjArray(); fMCcl->SetOwner(kTRUE); | |
e3147647 | 229 | PostData(kClToTrk, fCl); |
230 | PostData(kClToMC, fMCcl); | |
1ee39b3a | 231 | } |
232 | ||
233 | //________________________________________________________ | |
f8f46e4d | 234 | void AliTRDresolution::UserExec(Option_t *opt) |
1ee39b3a | 235 | { |
236 | // | |
237 | // Execution part | |
238 | // | |
239 | ||
3ed01fbe | 240 | if(fCl) fCl->Delete(); |
241 | if(fMCcl) fMCcl->Delete(); | |
b4414720 | 242 | AliTRDrecoTask::UserExec(opt); |
1ee39b3a | 243 | } |
244 | ||
553528eb | 245 | //________________________________________________________ |
6475ec36 | 246 | Bool_t AliTRDresolution::Pulls(Double_t* /*dyz[2]*/, Double_t* /*cov[3]*/, Double_t /*tilt*/) const |
553528eb | 247 | { |
248 | // Helper function to calculate pulls in the yz plane | |
249 | // using proper tilt rotation | |
250 | // Uses functionality defined by AliTRDseedV1. | |
251 | ||
6475ec36 | 252 | return kTRUE; |
253 | /* | |
553528eb | 254 | Double_t t2(tilt*tilt); |
9dcc64cc | 255 | // exit door until a bug fix is found for AliTRDseedV1::GetCovSqrt |
553528eb | 256 | |
257 | // rotate along pad | |
258 | Double_t cc[3]; | |
259 | cc[0] = cov[0] - 2.*tilt*cov[1] + t2*cov[2]; | |
260 | cc[1] = cov[1]*(1.-t2) + tilt*(cov[0] - cov[2]); | |
261 | cc[2] = t2*cov[0] + 2.*tilt*cov[1] + cov[2]; | |
262 | // do sqrt | |
263 | Double_t sqr[3]={0., 0., 0.}; | |
264 | if(AliTRDseedV1::GetCovSqrt(cc, sqr)) return kFALSE; | |
265 | Double_t invsqr[3]={0., 0., 0.}; | |
266 | if(AliTRDseedV1::GetCovInv(sqr, invsqr)<1.e-40) return kFALSE; | |
267 | Double_t tmp(dyz[0]); | |
268 | dyz[0] = invsqr[0]*tmp + invsqr[1]*dyz[1]; | |
269 | dyz[1] = invsqr[1]*tmp + invsqr[2]*dyz[1]; | |
270 | return kTRUE; | |
6475ec36 | 271 | */ |
553528eb | 272 | } |
273 | ||
1ee39b3a | 274 | //________________________________________________________ |
3ceb45ae | 275 | TH1* AliTRDresolution::PlotCluster(const AliTRDtrackV1 *track) |
1ee39b3a | 276 | { |
277 | // | |
3ceb45ae | 278 | // Plot the cluster distributions |
1ee39b3a | 279 | // |
280 | ||
281 | if(track) fkTrack = track; | |
282 | if(!fkTrack){ | |
3d2a3dff | 283 | AliDebug(4, "No Track defined."); |
b91fdd71 | 284 | return NULL; |
1ee39b3a | 285 | } |
3ed01fbe | 286 | if(TMath::Abs(fkESD->GetTOFbc()) > 1){ |
3ceb45ae | 287 | AliDebug(4, Form("Track with BC_index[%d] not used.", fkESD->GetTOFbc())); |
b91fdd71 | 288 | return NULL; |
1ee39b3a | 289 | } |
3ceb45ae | 290 | if(fPt<fPtThreshold){ |
291 | AliDebug(4, Form("Track with pt[%6.4f] under threshold.", fPt)); | |
b91fdd71 | 292 | return NULL; |
1ee39b3a | 293 | } |
3ceb45ae | 294 | THnSparse *H(NULL); |
295 | if(!fContainer || !(H = ((THnSparse*)fContainer->At(kCluster)))){ | |
1ee39b3a | 296 | AliWarning("No output container defined."); |
b91fdd71 | 297 | return NULL; |
1ee39b3a | 298 | } |
3ceb45ae | 299 | |
300 | AliTRDgeometry *geo(AliTRDinfoGen::Geometry()); | |
3ed01fbe | 301 | Double_t val[kNdim]; //Float_t exb, vd, t0, s2, dl, dt; |
3ceb45ae | 302 | TObjArray *clInfoArr(NULL); |
303 | AliTRDseedV1 *fTracklet(NULL); | |
3ed01fbe | 304 | AliTRDcluster *c(NULL), *cc(NULL); |
1ee39b3a | 305 | for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){ |
306 | if(!(fTracklet = fkTrack->GetTracklet(ily))) continue; | |
307 | if(!fTracklet->IsOK()) continue; | |
3ed01fbe | 308 | //fTracklet->GetCalibParam(exb, vd, t0, s2, dl, dt); |
3ceb45ae | 309 | val[kBC] = ily; |
310 | val[kPhi] = fPhi; | |
311 | val[kEta] = fEta; | |
3ed01fbe | 312 | val[kPt] = TMath::ATan(fTracklet->GetYref(1))*TMath::RadToDeg(); |
313 | Float_t corr = 1./TMath::Sqrt(1.+fTracklet->GetYref(1)*fTracklet->GetYref(1)+fTracklet->GetZref(1)*fTracklet->GetZref(1)); | |
314 | Int_t row0(-1); | |
315 | Float_t padCorr(fTracklet->GetTilt()*fTracklet->GetPadLength()); | |
3ceb45ae | 316 | fTracklet->ResetClusterIter(kTRUE); |
317 | while((c = fTracklet->NextCluster())){ | |
3ed01fbe | 318 | Float_t xc(c->GetX()), |
f429b017 | 319 | q(TMath::Abs(c->GetQ())); |
3ed01fbe | 320 | if(row0<0) row0 = c->GetPadRow(); |
321 | ||
322 | val[kYrez] = c->GetY() + padCorr*(c->GetPadRow() - row0) -fTracklet->GetYat(xc); | |
323 | val[kPrez] = fTracklet->GetX0()-xc; | |
f429b017 | 324 | val[kZrez] = 0.; Int_t ic(0), tb(c->GetLocalTimeBin());; |
325 | if((cc = fTracklet->GetClusters(tb-1))) {val[kZrez] += TMath::Abs(cc->GetQ()); ic++;} | |
326 | if((cc = fTracklet->GetClusters(tb-2))) {val[kZrez] += TMath::Abs(cc->GetQ()); ic++;} | |
3ed01fbe | 327 | if(ic) val[kZrez] /= (ic*q); |
328 | val[kSpeciesChgRC]= fTracklet->IsRowCross()?0.:(TMath::Max(q*corr, Float_t(3.))); | |
3ceb45ae | 329 | H->Fill(val); |
330 | /* // tilt rotation of covariance for clusters | |
553528eb | 331 | Double_t sy2(c->GetSigmaY2()), sz2(c->GetSigmaZ2()); |
a47a87c5 | 332 | cov[0] = (sy2+t2*sz2)*corr; |
333 | cov[1] = tilt*(sz2 - sy2)*corr; | |
334 | cov[2] = (t2*sy2+sz2)*corr; | |
553528eb | 335 | // sum with track covariance |
a47a87c5 | 336 | cov[0]+=covR[0]; cov[1]+=covR[1]; cov[2]+=covR[2]; |
553528eb | 337 | Double_t dyz[2]= {dy[1], dz[1]}; |
3ceb45ae | 338 | Pulls(dyz, cov, tilt);*/ |
1ee39b3a | 339 | |
dfd7d48b | 340 | // Get z-position with respect to anode wire |
3ceb45ae | 341 | Float_t yt(fTracklet->GetYref(0)-val[kZrez]*fTracklet->GetYref(1)), |
342 | zt(fTracklet->GetZref(0)-val[kZrez]*fTracklet->GetZref(1)); | |
fd7ffd88 | 343 | Int_t istk = geo->GetStack(c->GetDetector()); |
344 | AliTRDpadPlane *pp = geo->GetPadPlane(ily, istk); | |
3ed01fbe | 345 | Float_t rowZ = pp->GetRow0(); |
346 | Float_t d = rowZ - zt + pp->GetAnodeWireOffset(); | |
dfd7d48b | 347 | d -= ((Int_t)(2 * d)) / 2.0; |
348 | if (d > 0.25) d = 0.5 - d; | |
349 | ||
8ee59659 | 350 | AliTRDclusterInfo *clInfo(NULL); |
351 | clInfo = new AliTRDclusterInfo; | |
dfd7d48b | 352 | clInfo->SetCluster(c); |
3ceb45ae | 353 | //Float_t covF[] = {cov[0], cov[1], cov[2]}; |
354 | clInfo->SetGlobalPosition(yt, zt, fTracklet->GetYref(1), fTracklet->GetZref(1)/*, covF*/); | |
355 | clInfo->SetResolution(val[kYrez]); | |
dfd7d48b | 356 | clInfo->SetAnisochronity(d); |
3ceb45ae | 357 | clInfo->SetDriftLength(val[kZrez]); |
358 | clInfo->SetTilt(fTracklet->GetTilt()); | |
83b44483 | 359 | if(fCl) fCl->Add(clInfo); |
3ed01fbe | 360 | //else AliDebug(1, "Cl exchange container missing. Activate by calling \"InitExchangeContainers()\""); |
83b44483 | 361 | |
3ed01fbe | 362 | if(DebugLevel()>=2){ |
d14a8ac2 | 363 | if(!clInfoArr){ |
364 | clInfoArr=new TObjArray(AliTRDseedV1::kNclusters); | |
365 | clInfoArr->SetOwner(kFALSE); | |
366 | } | |
dfd7d48b | 367 | clInfoArr->Add(clInfo); |
df0514f6 | 368 | } |
1ee39b3a | 369 | } |
3ed01fbe | 370 | if(DebugLevel()>=2 && clInfoArr){ |
3ceb45ae | 371 | ULong_t status = fkESD->GetStatus(); |
dfd7d48b | 372 | (*DebugStream()) << "cluster" |
373 | <<"status=" << status | |
374 | <<"clInfo.=" << clInfoArr | |
375 | << "\n"; | |
d14a8ac2 | 376 | clInfoArr->Clear(); |
dfd7d48b | 377 | } |
1ee39b3a | 378 | } |
d14a8ac2 | 379 | if(clInfoArr) delete clInfoArr; |
00a56108 | 380 | |
3ceb45ae | 381 | return NULL;//H->Projection(kEta, kPhi); |
1ee39b3a | 382 | } |
383 | ||
384 | ||
385 | //________________________________________________________ | |
386 | TH1* AliTRDresolution::PlotTracklet(const AliTRDtrackV1 *track) | |
387 | { | |
388 | // Plot normalized residuals for tracklets to track. | |
389 | // | |
390 | // We start from the result that if X=N(|m|, |Cov|) | |
391 | // BEGIN_LATEX | |
392 | // (Cov^{-1})^{1/2}X = N((Cov^{-1})^{1/2}*|m|, |1|) | |
393 | // END_LATEX | |
394 | // in our case X=(y_trklt - y_trk z_trklt - z_trk) and |Cov| = |Cov_trklt| + |Cov_trk| at the radial | |
395 | // reference position. | |
396 | if(track) fkTrack = track; | |
397 | if(!fkTrack){ | |
3d2a3dff | 398 | AliDebug(4, "No Track defined."); |
b91fdd71 | 399 | return NULL; |
1ee39b3a | 400 | } |
3ed01fbe | 401 | if(TMath::Abs(fkESD->GetTOFbc())>1){ |
3ceb45ae | 402 | AliDebug(4, Form("Track with BC_index[%d] not used.", fkESD->GetTOFbc())); |
403 | return NULL; | |
404 | } | |
405 | THnSparse *H(NULL); | |
406 | if(!fContainer || !(H = (THnSparse*)fContainer->At(kTracklet))){ | |
1ee39b3a | 407 | AliWarning("No output container defined."); |
b91fdd71 | 408 | return NULL; |
1ee39b3a | 409 | } |
dbb2e0a7 | 410 | // return NULL; |
3ceb45ae | 411 | Double_t val[kNdim+1]; |
412 | AliTRDseedV1 *fTracklet(NULL); | |
3ba3e0a4 | 413 | for(Int_t il(0); il<AliTRDgeometry::kNlayer; il++){ |
1ee39b3a | 414 | if(!(fTracklet = fkTrack->GetTracklet(il))) continue; |
566c3d46 | 415 | if(!fTracklet->IsOK() || !fTracklet->IsChmbGood()) continue; |
3ceb45ae | 416 | val [kBC] = il; |
417 | val[kPhi] = fPhi; | |
418 | val[kEta] = fEta; | |
566c3d46 | 419 | val[kSpeciesChgRC]= fTracklet->IsRowCross()?0:fkTrack->Charge();// fSpecies; |
420 | val[kPt] = fPt<0.8?0:(fPt<1.5?1:2);//GetPtBin(fTracklet->GetMomentum()); | |
35983729 | 421 | Double_t dyt(fTracklet->GetYfit(0) - fTracklet->GetYref(0)), |
422 | dzt(fTracklet->GetZfit(0) - fTracklet->GetZref(0)), | |
3ceb45ae | 423 | dydx(fTracklet->GetYfit(1)), |
424 | tilt(fTracklet->GetTilt()); | |
425 | // correct for tilt rotation | |
426 | val[kYrez] = dyt - dzt*tilt; | |
427 | val[kZrez] = dzt + dyt*tilt; | |
428 | dydx+= tilt*fTracklet->GetZref(1); | |
35983729 | 429 | val[kPrez] = TMath::ATan((dydx - fTracklet->GetYref(1))/(1.+ fTracklet->GetYref(1)*dydx)) * TMath::RadToDeg(); |
3ceb45ae | 430 | if(fTracklet->IsRowCross()){ |
431 | val[kSpeciesChgRC]= 0.; | |
3ed01fbe | 432 | // val[kPrez] = fkTrack->Charge(); // may be better defined |
433 | }/* else { | |
3ceb45ae | 434 | Float_t exb, vd, t0, s2, dl, dt; |
435 | fTracklet->GetCalibParam(exb, vd, t0, s2, dl, dt); | |
436 | val[kZrez] = TMath::ATan((fTracklet->GetYref(1) - exb)/(1+fTracklet->GetYref(1)*exb)); | |
3ed01fbe | 437 | }*/ |
3ceb45ae | 438 | val[kNdim] = fTracklet->GetdQdl(); |
566c3d46 | 439 | H->Fill(val); |
3ceb45ae | 440 | |
441 | // // compute covariance matrix | |
442 | // fTracklet->GetCovAt(x, cov); | |
443 | // fTracklet->GetCovRef(covR); | |
444 | // cov[0] += covR[0]; cov[1] += covR[1]; cov[2] += covR[2]; | |
445 | // Double_t dyz[2]= {dy[1], dz[1]}; | |
446 | // Pulls(dyz, cov, tilt); | |
447 | // ((TH3S*)arr->At(1))->Fill(sgm[fSegmentLevel], dyz[0], dyz[1]); | |
448 | // ((TH3S*)arr->At(3))->Fill(tht, dyz[1], rc); | |
dfd7d48b | 449 | |
3ed01fbe | 450 | if(DebugLevel()>=3){ |
3ceb45ae | 451 | Bool_t rc(fTracklet->IsRowCross()); |
dfd7d48b | 452 | UChar_t err(fTracklet->GetErrorMsg()); |
3ceb45ae | 453 | Double_t x(fTracklet->GetX()), |
454 | pt(fTracklet->GetPt()), | |
455 | yt(fTracklet->GetYref(0)), | |
456 | zt(fTracklet->GetZref(0)), | |
457 | phi(fTracklet->GetYref(1)), | |
458 | tht(fTracklet->GetZref(1)); | |
459 | Int_t ncl(fTracklet->GetN()), | |
460 | det(fTracklet->GetDetector()); | |
dfd7d48b | 461 | (*DebugStream()) << "tracklet" |
3ba3e0a4 | 462 | <<"pt=" << pt |
3ceb45ae | 463 | <<"x=" << x |
6475ec36 | 464 | <<"yt=" << yt |
3ceb45ae | 465 | <<"zt=" << zt |
3ba3e0a4 | 466 | <<"phi=" << phi |
467 | <<"tht=" << tht | |
3ceb45ae | 468 | <<"det=" << det |
6475ec36 | 469 | <<"n=" << ncl |
3ceb45ae | 470 | <<"dy0=" << dyt |
471 | <<"dz0=" << dzt | |
472 | <<"dy=" << val[kYrez] | |
473 | <<"dz=" << val[kZrez] | |
474 | <<"dphi="<< val[kPrez] | |
475 | <<"dQ ="<< val[kNdim] | |
dfd7d48b | 476 | <<"rc=" << rc |
477 | <<"err=" << err | |
478 | << "\n"; | |
479 | } | |
1ee39b3a | 480 | } |
3ceb45ae | 481 | return NULL;//H->Projection(kEta, kPhi); |
1ee39b3a | 482 | } |
483 | ||
484 | ||
485 | //________________________________________________________ | |
a310e49b | 486 | TH1* AliTRDresolution::PlotTrackIn(const AliTRDtrackV1 *track) |
1ee39b3a | 487 | { |
488 | // Store resolution/pulls of Kalman before updating with the TRD information | |
489 | // at the radial position of the first tracklet. The following points are used | |
490 | // for comparison | |
491 | // - the (y,z,snp) of the first TRD tracklet | |
492 | // - the (y, z, snp, tgl, pt) of the MC track reference | |
493 | // | |
494 | // Additionally the momentum resolution/pulls are calculated for usage in the | |
495 | // PID calculation. | |
3ed01fbe | 496 | //printf("AliTRDresolution::PlotTrackIn() :: track[%p]\n", (void*)track); |
497 | ||
1ee39b3a | 498 | if(track) fkTrack = track; |
3ceb45ae | 499 | if(!fkTrack){ |
3d2a3dff | 500 | AliDebug(4, "No Track defined."); |
b91fdd71 | 501 | return NULL; |
1ee39b3a | 502 | } |
3ed01fbe | 503 | //fkTrack->Print(); |
3ceb45ae | 504 | // check container |
505 | THnSparseI *H=(THnSparseI*)fContainer->At(kTrackIn); | |
506 | if(!H){ | |
507 | AliError(Form("Missing container @ %d", Int_t(kTrackIn))); | |
83b44483 | 508 | return NULL; |
509 | } | |
3ceb45ae | 510 | // check input track status |
511 | AliExternalTrackParam *tin(NULL); | |
a310e49b | 512 | if(!(tin = fkTrack->GetTrackIn())){ |
3ceb45ae | 513 | AliError("Track did not entered TRD fiducial volume."); |
b91fdd71 | 514 | return NULL; |
1ee39b3a | 515 | } |
3ceb45ae | 516 | // check first tracklet |
517 | AliTRDseedV1 *fTracklet(fkTrack->GetTracklet(0)); | |
518 | if(!fTracklet){ | |
519 | AliDebug(3, "No Tracklet in ly[0]. Skip track."); | |
b91fdd71 | 520 | return NULL; |
1ee39b3a | 521 | } |
566c3d46 | 522 | if(!fTracklet->IsOK() || !fTracklet->IsChmbGood()){ |
523 | AliDebug(3, "Tracklet or Chamber not OK. Skip track."); | |
524 | return NULL; | |
525 | } | |
3ceb45ae | 526 | // check radial position |
527 | Double_t x = tin->GetX(); | |
528 | if(TMath::Abs(x-fTracklet->GetX())>1.e-3){ | |
529 | AliDebug(1, Form("Tracklet did not match Track. dx[cm]=%+4.1f", x-fTracklet->GetX())); | |
530 | return NULL; | |
1ee39b3a | 531 | } |
3ed01fbe | 532 | //printf("USE y[%+f] dydx[%+f]\n", fTracklet->GetYfit(0), fTracklet->GetYfit(1)); |
1ee39b3a | 533 | |
566c3d46 | 534 | Int_t bc(fkESD->GetTOFbc()/2); |
3ceb45ae | 535 | const Double_t *parR(tin->GetParameter()); |
35983729 | 536 | Double_t dyt(fTracklet->GetYfit(0)-parR[0]), dzt(fTracklet->GetZfit(0)-parR[1]), |
3ceb45ae | 537 | phit(fTracklet->GetYfit(1)), |
35983729 | 538 | tilt(fTracklet->GetTilt()), |
539 | norm(1./TMath::Sqrt((1.-parR[2])*(1.+parR[2]))); | |
3ceb45ae | 540 | |
541 | // correct for tilt rotation | |
542 | Double_t dy = dyt - dzt*tilt, | |
35983729 | 543 | dz = dzt + dyt*tilt, |
544 | dx = dy/(parR[2]*norm-parR[3]*norm*tilt); | |
3ceb45ae | 545 | phit += tilt*parR[3]; |
35983729 | 546 | Double_t dphi = TMath::ATan(phit) - TMath::ASin(parR[2]); |
3ceb45ae | 547 | |
566c3d46 | 548 | Double_t val[kNdim+3]; |
549 | val[kBC] = bc==0?0:(bc<0?-1.:1.); | |
3ceb45ae | 550 | val[kPhi] = fPhi; |
551 | val[kEta] = fEta; | |
35983729 | 552 | val[kSpeciesChgRC]= fTracklet->IsRowCross()?0:fkTrack->Charge(); |
553 | val[kPt] = fPt<0.8?0:(fPt<1.5?1:2);//GetPtBin(fPt); | |
3ceb45ae | 554 | val[kYrez] = dy; |
555 | val[kZrez] = dz; | |
556 | val[kPrez] = dphi*TMath::RadToDeg(); | |
35983729 | 557 | val[kNdim] = fTracklet->GetDetector(); |
558 | val[kNdim+1] = dx; | |
566c3d46 | 559 | val[kNdim+2] = fEvent->GetBunchFill(); |
3ceb45ae | 560 | H->Fill(val); |
3ed01fbe | 561 | if(DebugLevel()>=3){ |
562 | (*DebugStream()) << "trackIn" | |
563 | <<"tracklet.=" << fTracklet | |
564 | <<"trackIn.=" << tin | |
565 | << "\n"; | |
566 | } | |
3ceb45ae | 567 | |
f429b017 | 568 | return NULL; // H->Projection(kEta, kPhi); |
1ee39b3a | 569 | } |
570 | ||
f429b017 | 571 | /* |
a310e49b | 572 | //________________________________________________________ |
573 | TH1* AliTRDresolution::PlotTrackOut(const AliTRDtrackV1 *track) | |
574 | { | |
575 | // Store resolution/pulls of Kalman after last update with the TRD information | |
576 | // at the radial position of the first tracklet. The following points are used | |
577 | // for comparison | |
578 | // - the (y,z,snp) of the first TRD tracklet | |
579 | // - the (y, z, snp, tgl, pt) of the MC track reference | |
580 | // | |
581 | // Additionally the momentum resolution/pulls are calculated for usage in the | |
582 | // PID calculation. | |
583 | ||
584 | if(track) fkTrack = track; | |
596f82fb | 585 | return NULL; |
a310e49b | 586 | } |
f429b017 | 587 | */ |
1ee39b3a | 588 | //________________________________________________________ |
589 | TH1* AliTRDresolution::PlotMC(const AliTRDtrackV1 *track) | |
590 | { | |
591 | // | |
592 | // Plot MC distributions | |
593 | // | |
594 | ||
595 | if(!HasMCdata()){ | |
e9d62bd2 | 596 | AliDebug(2, "No MC defined. Results will not be available."); |
b91fdd71 | 597 | return NULL; |
1ee39b3a | 598 | } |
599 | if(track) fkTrack = track; | |
600 | if(!fkTrack){ | |
3d2a3dff | 601 | AliDebug(4, "No Track defined."); |
b91fdd71 | 602 | return NULL; |
1ee39b3a | 603 | } |
f429b017 | 604 | Int_t bc(TMath::Abs(fkESD->GetTOFbc())); |
605 | ||
606 | THnSparse *H(NULL); | |
83b44483 | 607 | if(!fContainer){ |
608 | AliWarning("No output container defined."); | |
609 | return NULL; | |
610 | } | |
92c40e64 | 611 | // retriev track characteristics |
612 | Int_t pdg = fkMC->GetPDG(), | |
613 | sIdx(AliTRDpidUtil::Pdg2Pid(TMath::Abs(pdg))+1), // species index | |
614 | sign(0), | |
f429b017 | 615 | // sgm[3], |
616 | label(fkMC->GetLabel()); | |
617 | // fSegmentLevel(0); | |
92c40e64 | 618 | if(!fDBPDG) fDBPDG=TDatabasePDG::Instance(); |
619 | TParticlePDG *ppdg(fDBPDG->GetParticle(pdg)); | |
620 | if(ppdg) sign = ppdg->Charge() > 0. ? 1 : -1; | |
92c40e64 | 621 | |
f429b017 | 622 | TH1 *h(NULL); |
3ceb45ae | 623 | AliTRDgeometry *geo(AliTRDinfoGen::Geometry()); |
624 | AliTRDseedV1 *fTracklet(NULL); TObjArray *clInfoArr(NULL); | |
1ee39b3a | 625 | UChar_t s; |
f429b017 | 626 | Double_t x, y, z, pt, dydx, dzdx, dzdl; |
1ee39b3a | 627 | Float_t pt0, x0, y0, z0, dx, dy, dz, dydx0, dzdx0; |
628 | Double_t covR[7]/*, cov[3]*/; | |
3ceb45ae | 629 | |
f429b017 | 630 | /* if(DebugLevel()>=3){ |
3ceb45ae | 631 | // get first detector |
632 | Int_t det = -1; | |
633 | for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){ | |
634 | if(!(fTracklet = fkTrack->GetTracklet(ily))) continue; | |
635 | det = fTracklet->GetDetector(); | |
636 | break; | |
637 | } | |
638 | if(det>=0){ | |
639 | TVectorD X(12), Y(12), Z(12), dX(12), dY(12), dZ(12), vPt(12), dPt(12), budget(12), cCOV(12*15); | |
640 | Double_t m(-1.); | |
641 | m = fkTrack->GetMass(); | |
642 | if(fkMC->PropagateKalman(&X, &Y, &Z, &dX, &dY, &dZ, &vPt, &dPt, &budget, &cCOV, m)){ | |
643 | (*DebugStream()) << "MCkalman" | |
644 | << "pdg=" << pdg | |
645 | << "det=" << det | |
646 | << "x=" << &X | |
647 | << "y=" << &Y | |
648 | << "z=" << &Z | |
649 | << "dx=" << &dX | |
650 | << "dy=" << &dY | |
651 | << "dz=" << &dZ | |
652 | << "pt=" << &vPt | |
653 | << "dpt=" << &dPt | |
654 | << "bgt=" << &budget | |
655 | << "cov=" << &cCOV | |
656 | << "\n"; | |
657 | } | |
658 | } | |
f429b017 | 659 | }*/ |
660 | AliTRDcluster *c(NULL); | |
661 | Double_t val[kNdim+1]; | |
1ee39b3a | 662 | for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){ |
663 | if(!(fTracklet = fkTrack->GetTracklet(ily)))/* || | |
664 | !fTracklet->IsOK())*/ continue; | |
665 | ||
f429b017 | 666 | x= x0 = fTracklet->GetX(); |
667 | Bool_t rc(fTracklet->IsRowCross()); Float_t eta, phi; | |
668 | if(!fkMC->GetDirections(x0, y0, z0, dydx0, dzdx0, pt0, eta, phi, s)) continue; | |
1ee39b3a | 669 | |
670 | // MC track position at reference radial position | |
671 | dx = x0 - x; | |
1ee39b3a | 672 | Float_t ymc = y0 - dx*dydx0; |
673 | Float_t zmc = z0 - dx*dzdx0; | |
f429b017 | 674 | phi -= TMath::Pi(); |
675 | ||
676 | val[kBC] = ily; | |
677 | val[kPhi] = phi; | |
678 | val[kEta] = eta; | |
679 | val[kSpeciesChgRC]= rc?0.:sign*sIdx; | |
35983729 | 680 | val[kPt] = pt0<0.8?0:1;//GetPtBin(pt0); |
f429b017 | 681 | Double_t tilt(fTracklet->GetTilt()); |
682 | // ,t2(tilt*tilt) | |
683 | // ,corr(1./(1. + t2)) | |
684 | // ,cost(TMath::Sqrt(corr)); | |
685 | ||
686 | AliExternalTrackParam *tin(fkTrack->GetTrackIn()); | |
687 | if(ily==0 && tin){ // trackIn residuals | |
688 | // check radial position | |
689 | if(TMath::Abs(tin->GetX()-x)>1.e-3) AliDebug(1, Form("TrackIn radial mismatch. dx[cm]=%+4.1f", tin->GetX()-x)); | |
690 | else{ | |
691 | // Float_t phi = TMath::ATan2(y0, x0); | |
692 | val[kBC] = (bc>=kNbunchCross)?(kNbunchCross-1):bc; | |
693 | val[kYrez] = tin->GetY()-ymc; | |
694 | val[kZrez] = tin->GetZ()-zmc; | |
695 | val[kPrez] = (TMath::ASin(tin->GetSnp())-TMath::ATan(dydx0))*TMath::RadToDeg(); | |
696 | if((H = (THnSparseI*)fContainer->At(kMCtrackIn))) H->Fill(val); | |
697 | } | |
698 | } | |
699 | if(bc>1) break; // do nothing for the rest of TRD objects if satellite bunch | |
1ee39b3a | 700 | |
f429b017 | 701 | // track residuals |
1ee39b3a | 702 | dydx = fTracklet->GetYref(1); |
703 | dzdx = fTracklet->GetZref(1); | |
704 | dzdl = fTracklet->GetTgl(); | |
f429b017 | 705 | y = fTracklet->GetYref(0); |
1ee39b3a | 706 | dy = y - ymc; |
f429b017 | 707 | z = fTracklet->GetZref(0); |
1ee39b3a | 708 | dz = z - zmc; |
709 | pt = TMath::Abs(fTracklet->GetPt()); | |
710 | fTracklet->GetCovRef(covR); | |
711 | ||
f429b017 | 712 | val[kYrez] = dy; |
713 | val[kPrez] = TMath::ATan((dydx - dydx0)/(1.+ dydx*dydx0))*TMath::RadToDeg(); | |
714 | val[kZrez] = dz; | |
6465da91 | 715 | val[kNdim] = 1.e2*(pt/pt0-1.); |
f429b017 | 716 | if((H = (THnSparse*)fContainer->At(kMCtrack))) H->Fill(val); |
717 | /* // theta resolution/ tgl pulls | |
718 | Double_t dzdl0 = dzdx0/TMath::Sqrt(1.+dydx0*dydx0), | |
719 | dtgl = (dzdl - dzdl0)/(1.- dzdl*dzdl0); | |
720 | ((TH2I*)arr->At(6))->Fill(dzdl0, TMath::ATan(dtgl)); | |
721 | ((TH2I*)arr->At(7))->Fill(dzdl0, (dzdl - dzdl0)/TMath::Sqrt(covR[4])); | |
722 | // pt resolution \\ 1/pt pulls \\ p resolution for PID | |
723 | Double_t p0 = TMath::Sqrt(1.+ dzdl0*dzdl0)*pt0, | |
724 | p = TMath::Sqrt(1.+ dzdl*dzdl)*pt; | |
725 | ((TH3S*)((TObjArray*)arr->At(8)))->Fill(pt0, pt/pt0-1., sign*sIdx); | |
726 | ((TH3S*)((TObjArray*)arr->At(9)))->Fill(1./pt0, (1./pt-1./pt0)/TMath::Sqrt(covR[6]), sign*sIdx); | |
727 | ((TH3S*)((TObjArray*)arr->At(10)))->Fill(p0, p/p0-1., sign*sIdx);*/ | |
728 | ||
729 | // Fill Debug stream for MC track | |
3ba3e0a4 | 730 | if(DebugLevel()>=4){ |
f429b017 | 731 | Int_t det(fTracklet->GetDetector()); |
732 | (*DebugStream()) << "MC" | |
733 | << "det=" << det | |
734 | << "pdg=" << pdg | |
735 | << "sgn=" << sign | |
736 | << "pt=" << pt0 | |
737 | << "x=" << x0 | |
738 | << "y=" << y0 | |
739 | << "z=" << z0 | |
740 | << "dydx=" << dydx0 | |
741 | << "dzdx=" << dzdx0 | |
742 | << "\n"; | |
743 | ||
744 | // Fill Debug stream for Kalman track | |
1ee39b3a | 745 | (*DebugStream()) << "MCtrack" |
746 | << "pt=" << pt | |
747 | << "x=" << x | |
748 | << "y=" << y | |
749 | << "z=" << z | |
750 | << "dydx=" << dydx | |
751 | << "dzdx=" << dzdx | |
752 | << "s2y=" << covR[0] | |
753 | << "s2z=" << covR[2] | |
754 | << "\n"; | |
755 | } | |
756 | ||
f429b017 | 757 | // tracklet residuals |
758 | dydx = fTracklet->GetYfit(1) + tilt*dzdx0; | |
759 | dzdx = fTracklet->GetZfit(1); | |
760 | y = fTracklet->GetYfit(0); | |
761 | dy = y - ymc; | |
762 | z = fTracklet->GetZfit(0); | |
763 | dz = z - zmc; | |
764 | val[kYrez] = dy - dz*tilt; | |
765 | val[kPrez] = TMath::ATan((dydx - dydx0)/(1.+ dydx*dydx0))*TMath::RadToDeg(); | |
766 | val[kZrez] = dz + dy*tilt; | |
767 | // val[kNdim] = pt/pt0-1.; | |
768 | if((H = (THnSparse*)fContainer->At(kMCtracklet))) H->Fill(val); | |
769 | ||
770 | ||
1ee39b3a | 771 | // Fill Debug stream for tracklet |
3ba3e0a4 | 772 | if(DebugLevel()>=4){ |
f429b017 | 773 | Float_t s2y = fTracklet->GetS2Y(); |
774 | Float_t s2z = fTracklet->GetS2Z(); | |
1ee39b3a | 775 | (*DebugStream()) << "MCtracklet" |
776 | << "rc=" << rc | |
777 | << "x=" << x | |
778 | << "y=" << y | |
779 | << "z=" << z | |
780 | << "dydx=" << dydx | |
781 | << "s2y=" << s2y | |
782 | << "s2z=" << s2z | |
783 | << "\n"; | |
784 | } | |
785 | ||
f429b017 | 786 | AliTRDpadPlane *pp = geo->GetPadPlane(ily, AliTRDgeometry::GetStack(fTracklet->GetDetector())); |
1ee39b3a | 787 | Float_t zr0 = pp->GetRow0() + pp->GetAnodeWireOffset(); |
1ee39b3a | 788 | //Double_t exb = AliTRDCommonParam::Instance()->GetOmegaTau(1.5); |
789 | ||
f429b017 | 790 | H = (THnSparse*)fContainer->At(kMCcluster); |
791 | val[kPt] = TMath::ATan(dydx0)*TMath::RadToDeg(); | |
792 | //Float_t corr = 1./TMath::Sqrt(1.+dydx0*dydx0+dzdx0*dzdx0); | |
793 | Int_t row0(-1); | |
794 | Float_t padCorr(tilt*fTracklet->GetPadLength()); | |
795 | fTracklet->ResetClusterIter(kTRUE); | |
796 | while((c = fTracklet->NextCluster())){ | |
797 | if(row0<0) row0 = c->GetPadRow(); | |
798 | x = c->GetX();//+fXcorr[c->GetDetector()][c->GetLocalTimeBin()]; | |
799 | y = c->GetY() + padCorr*(c->GetPadRow() - row0); | |
800 | z = c->GetZ(); | |
1ee39b3a | 801 | dx = x0 - x; |
802 | ymc= y0 - dx*dydx0; | |
803 | zmc= z0 - dx*dzdx0; | |
f429b017 | 804 | dy = y - ymc; |
805 | dz = z - zmc; | |
806 | val[kYrez] = dy - dz*tilt; | |
807 | val[kPrez] = dx; | |
808 | val[kZrez] = 0.; AliTRDcluster *cc(NULL); Int_t ic(0), tb(c->GetLocalTimeBin()); Float_t q(TMath::Abs(c->GetQ())); | |
809 | if((cc = fTracklet->GetClusters(tb-1))) {val[kZrez] += TMath::Abs(cc->GetQ()); ic++;} | |
810 | if((cc = fTracklet->GetClusters(tb-2))) {val[kZrez] += TMath::Abs(cc->GetQ()); ic++;} | |
811 | if(ic) val[kZrez] /= (ic*q); | |
812 | if(H) H->Fill(val); | |
813 | ||
1ee39b3a | 814 | |
815 | // Fill calibration container | |
816 | Float_t d = zr0 - zmc; | |
817 | d -= ((Int_t)(2 * d)) / 2.0; | |
818 | if (d > 0.25) d = 0.5 - d; | |
819 | AliTRDclusterInfo *clInfo = new AliTRDclusterInfo; | |
1ee39b3a | 820 | clInfo->SetCluster(c); |
821 | clInfo->SetMC(pdg, label); | |
822 | clInfo->SetGlobalPosition(ymc, zmc, dydx0, dzdx0); | |
823 | clInfo->SetResolution(dy); | |
824 | clInfo->SetAnisochronity(d); | |
2589cf64 | 825 | clInfo->SetDriftLength(dx); |
1ee39b3a | 826 | clInfo->SetTilt(tilt); |
83b44483 | 827 | if(fMCcl) fMCcl->Add(clInfo); |
828 | else AliDebug(1, "MCcl exchange container missing. Activate by calling \"InitExchangeContainers()\""); | |
3ba3e0a4 | 829 | if(DebugLevel()>=5){ |
d14a8ac2 | 830 | if(!clInfoArr){ |
831 | clInfoArr=new TObjArray(AliTRDseedV1::kNclusters); | |
832 | clInfoArr->SetOwner(kFALSE); | |
833 | } | |
d43e2ad1 | 834 | clInfoArr->Add(clInfo); |
1ee39b3a | 835 | } |
836 | } | |
d43e2ad1 | 837 | // Fill Debug Tree |
3ba3e0a4 | 838 | if(DebugLevel()>=5 && clInfoArr){ |
d43e2ad1 | 839 | (*DebugStream()) << "MCcluster" |
840 | <<"clInfo.=" << clInfoArr | |
841 | << "\n"; | |
92c40e64 | 842 | clInfoArr->Clear(); |
d43e2ad1 | 843 | } |
1ee39b3a | 844 | } |
92c40e64 | 845 | if(clInfoArr) delete clInfoArr; |
1ee39b3a | 846 | return h; |
847 | } | |
848 | ||
849 | ||
3ceb45ae | 850 | //__________________________________________________________________________ |
851 | Int_t AliTRDresolution::GetPtBin(Float_t pt) | |
852 | { | |
853 | // Find pt bin according to local pt segmentation | |
3ed01fbe | 854 | Int_t ipt(-1); |
3ceb45ae | 855 | while(ipt<AliTRDresolution::kNpt){ |
3ed01fbe | 856 | if(pt<fgPtBin[ipt+1]) break; |
3ceb45ae | 857 | ipt++; |
858 | } | |
3ed01fbe | 859 | return ipt; |
3ceb45ae | 860 | } |
861 | ||
862 | //________________________________________________________ | |
3ed01fbe | 863 | Float_t AliTRDresolution::GetMeanStat(TH1 *h, Float_t cut, Option_t *opt) |
3ceb45ae | 864 | { |
3ed01fbe | 865 | // return mean number of entries/bin of histogram "h" |
866 | // if option "opt" is given the following values are accepted: | |
867 | // "<" : consider only entries less than "cut" | |
868 | // ">" : consider only entries greater than "cut" | |
869 | ||
870 | //Int_t dim(h->GetDimension()); | |
871 | Int_t nbx(h->GetNbinsX()), nby(h->GetNbinsY()), nbz(h->GetNbinsZ()); | |
872 | Double_t sum(0.); Int_t n(0); | |
873 | for(Int_t ix(1); ix<=nbx; ix++) | |
874 | for(Int_t iy(1); iy<=nby; iy++) | |
875 | for(Int_t iz(1); iz<=nbz; iz++){ | |
876 | if(strcmp(opt, "")==0){sum += h->GetBinContent(ix, iy, iz); n++;} | |
877 | else{ | |
878 | if(strcmp(opt, "<")==0) { | |
879 | if(h->GetBinContent(ix, iy, iz)<cut) {sum += h->GetBinContent(ix, iy, iz); n++;} | |
880 | } else if(strcmp(opt, ">")==0){ | |
881 | if(h->GetBinContent(ix, iy, iz)>cut) {sum += h->GetBinContent(ix, iy, iz); n++;} | |
882 | } else {sum += h->GetBinContent(ix, iy, iz); n++;} | |
883 | } | |
884 | } | |
885 | return n>0?sum/n:0.; | |
3ceb45ae | 886 | } |
887 | ||
1ee39b3a | 888 | //________________________________________________________ |
889 | Bool_t AliTRDresolution::GetRefFigure(Int_t ifig) | |
890 | { | |
891 | // | |
892 | // Get the reference figures | |
893 | // | |
894 | ||
1ee39b3a | 895 | if(!gPad){ |
896 | AliWarning("Please provide a canvas to draw results."); | |
897 | return kFALSE; | |
898 | } | |
3ceb45ae | 899 | /* Int_t selection[100], n(0), selStart(0); // |
a310e49b | 900 | Int_t ly0(0), dly(5); |
3ceb45ae | 901 | TList *l(NULL); TVirtualPad *pad(NULL); */ |
1ee39b3a | 902 | switch(ifig){ |
3ceb45ae | 903 | case 0: |
1ee39b3a | 904 | break; |
1ee39b3a | 905 | } |
906 | AliWarning(Form("Reference plot [%d] missing result", ifig)); | |
907 | return kFALSE; | |
908 | } | |
909 | ||
3ceb45ae | 910 | |
911 | //________________________________________________________ | |
912 | void AliTRDresolution::MakePtSegmentation(Float_t pt0, Float_t dpt) | |
913 | { | |
914 | // Build pt segments | |
915 | for(Int_t j(0); j<=kNpt; j++){ | |
916 | pt0+=(TMath::Exp(j*j*dpt)-1.); | |
917 | fgPtBin[j]=pt0; | |
918 | } | |
919 | } | |
920 | ||
00a3f7a4 | 921 | //________________________________________________________ |
922 | void AliTRDresolution::MakeSummary() | |
923 | { | |
924 | // Build summary plots | |
925 | ||
3ceb45ae | 926 | if(!fProj){ |
00a3f7a4 | 927 | AliError("Missing results"); |
928 | return; | |
929 | } | |
3ceb45ae | 930 | TVirtualPad *p(NULL); TCanvas *cOut(NULL); |
3ed01fbe | 931 | TObjArray *arr(NULL); TH2 *h2(NULL); |
068e2c00 | 932 | |
3ceb45ae | 933 | // cluster resolution |
934 | // define palette | |
935 | gStyle->SetPalette(1); | |
f429b017 | 936 | const Int_t nClViews(11); |
937 | const Char_t *vClName[nClViews] = {"ClY", "ClYn", "ClYp", "ClQn", "ClQp", "ClYXTCp", "ClYXTCn", "ClYXPh", "ClYXPh", "ClY", "ClYn"}; | |
938 | const UChar_t vClOpt[nClViews] = {1, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0}; | |
566c3d46 | 939 | const Int_t nTrkltViews(4); |
940 | const Char_t *vTrkltName[nTrkltViews][6] = { | |
941 | {"TrkltYn", "TrkltYp", "TrkltRCZ", "TrkltPhn", "TrkltPhp", "TrkltQ"}, // general view | |
942 | {"TrkltYnl", "TrkltYni", "TrkltYnh", "TrkltYpl", "TrkltYpi", "TrkltYph"}, // alignment view | |
943 | {"TrkltPhnl", "TrkltPhni", "TrkltPhnh", "TrkltPhpl", "TrkltPhpi", "TrkltPhph"}, // calibration view | |
944 | {"TrkltQnl", "TrkltQni", "TrkltQnh", "TrkltQpl", "TrkltQpi", "TrkltQph"} // PID view | |
945 | }; | |
946 | const Int_t nTrkInViews(5); | |
35983729 | 947 | const Char_t *vTrkInName[nTrkInViews][6] = { |
566c3d46 | 948 | {"TrkInY", "TrkInYn", "TrkInYp", "TrkInRCZ", "TrkInPhn", "TrkInPhp"}, |
949 | {"TrkInRCX", "TrkInRCY", "TrkInRCPh", "TrkInRCZl", "TrkInRCZi", "TrkInRCZh"}, | |
35983729 | 950 | {"TrkInYnl", "TrkInYni", "TrkInYnh", "TrkInYpl", "TrkInYpi", "TrkInYph"}, |
566c3d46 | 951 | {"TrkInXnl", "TrkInXpl", "TrkInXl", "TrkInRCXl", "TrkInRCYl", "TrkInYh"}, |
35983729 | 952 | {"TrkInPhnl", "TrkInPhni", "TrkInPhnh", "TrkInPhpl", "TrkInPhpi", "TrkInPhph"}}; |
f429b017 | 953 | const Int_t nTrkViews(10); |
954 | const Char_t *vTrkName[nTrkViews] = {"TrkY", "TrkYn", "TrkYp", "TrkPhn", "TrkPhp", "TrkZ", "TrkQn", "TrkQp", "TrkPn", "TrkPp"}; | |
955 | const Char_t *typName[] = {"", "MC"}; | |
956 | ||
957 | for(Int_t ityp(0); ityp<(HasMCdata()?2:1); ityp++){ | |
958 | if((arr = (TObjArray*)fProj->At(ityp?kMCcluster:kCluster))){ | |
959 | for(Int_t iview(0); iview<nClViews; iview++){ | |
566c3d46 | 960 | cOut = new TCanvas(Form("%s_%sCl%02d", GetName(), typName[ityp], iview), "Cluster Resolution", 1024, 768); |
f429b017 | 961 | cOut->Divide(3,2, 1.e-5, 1.e-5); |
962 | Int_t nplot(0); | |
963 | for(Int_t iplot(0); iplot<6; iplot++){ | |
964 | p=cOut->cd(iplot+1); p->SetRightMargin(0.1572581);p->SetTopMargin(0.08262712); | |
965 | if(!(h2 = (TH2*)arr->FindObject(Form("H%s%s%d_2D", typName[ityp], vClName[iview], iplot)))) continue; | |
966 | nplot++; | |
967 | if(vClOpt[iview]==0) h2->Draw("colz"); | |
968 | else if(vClOpt[iview]==1) DrawSigma(h2, 1.e4, 2.e2, 6.5e2, "#sigma(#Deltay) [#mum]"); | |
969 | } | |
970 | if(nplot) cOut->SaveAs(Form("%s.gif", cOut->GetName())); | |
971 | else delete cOut; | |
972 | } | |
973 | } | |
974 | // tracklet systematic | |
975 | if((arr = (TObjArray*)fProj->At(ityp?kMCtracklet:kTracklet))){ | |
976 | for(Int_t iview(0); iview<nTrkltViews; iview++){ | |
566c3d46 | 977 | for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){ |
978 | cOut = new TCanvas(Form("%s_%sTrklt%02d_%d", GetName(), typName[ityp], iview, ily), "Tracklet Resolution", 1024, 768); | |
979 | cOut->Divide(3,2, 1.e-5, 1.e-5); | |
980 | Int_t nplot(0); | |
981 | for(Int_t iplot(0); iplot<6; iplot++){ | |
982 | p=cOut->cd(iplot+1); p->SetRightMargin(0.1572581); p->SetTopMargin(0.08262712); | |
983 | if(!(h2 = (TH2*)arr->FindObject(Form("H%s%s%d_2D", typName[ityp], vTrkltName[iview][iplot], ily)))){ | |
984 | AliInfo(Form("Missing H%s%s%d_2D", typName[ityp], vTrkltName[iview][iplot], ily)); | |
985 | continue; | |
986 | } | |
987 | h2->Draw("colz"); nplot++; | |
988 | } | |
989 | if(nplot) cOut->SaveAs(Form("%s.gif", cOut->GetName())); | |
990 | else delete cOut; | |
f429b017 | 991 | } |
f429b017 | 992 | } |
993 | } | |
994 | // trackIn systematic | |
995 | if((arr = (TObjArray*)fProj->At(ityp?kMCtrackIn:kTrackIn))){ | |
35983729 | 996 | for(Int_t iview(0); iview<nTrkInViews; iview++){ |
566c3d46 | 997 | cOut = new TCanvas(Form("%s_%sTrkIn%02d", GetName(), typName[ityp], iview), "Track IN Resolution", 1024, 768); |
35983729 | 998 | cOut->Divide(3,2, 1.e-5, 1.e-5); |
999 | Int_t nplot(0); | |
1000 | for(Int_t iplot(0); iplot<6; iplot++){ | |
566c3d46 | 1001 | p=cOut->cd(iplot+1); p->SetRightMargin(0.1572581); p->SetTopMargin(0.08262712); |
1002 | if(!(h2 = (TH2*)arr->FindObject(Form("H%s%s_2D", typName[ityp], vTrkInName[iview][iplot])))){ | |
1003 | AliInfo(Form("Missing H%s%s_2D", typName[ityp], vTrkInName[iview][iplot])); | |
1004 | continue; | |
1005 | } | |
1006 | h2->Draw("colz"); nplot++; | |
35983729 | 1007 | } |
1008 | if(nplot) cOut->SaveAs(Form("%s.gif", cOut->GetName())); | |
1009 | else delete cOut; | |
3ed01fbe | 1010 | } |
3ed01fbe | 1011 | } |
3ceb45ae | 1012 | } |
f429b017 | 1013 | // track MC systematic |
1014 | if((arr = (TObjArray*)fProj->At(kMCtrack))) { | |
1015 | for(Int_t iview(0); iview<nTrkViews; iview++){ | |
566c3d46 | 1016 | cOut = new TCanvas(Form("%s_MCTrk%02d", GetName(), iview), "Track Resolution", 1024, 768); |
0b30040d | 1017 | cOut->Divide(3,2, 1.e-5, 1.e-5); |
09c85be5 | 1018 | Int_t nplot(0); |
3ed01fbe | 1019 | for(Int_t iplot(0); iplot<6; iplot++){ |
1020 | p=cOut->cd(iplot+1); p->SetRightMargin(0.1572581); p->SetTopMargin(0.08262712); | |
f429b017 | 1021 | if(!(h2 = (TH2*)arr->FindObject(Form("H%s%d_2D", vTrkName[iview], iplot)))) continue; |
09c85be5 | 1022 | h2->Draw("colz"); nplot++; |
3ed01fbe | 1023 | } |
09c85be5 | 1024 | if(nplot) cOut->SaveAs(Form("%s.gif", cOut->GetName())); |
1025 | else delete cOut; | |
3ed01fbe | 1026 | } |
1027 | } | |
f429b017 | 1028 | |
1029 | ||
3ceb45ae | 1030 | gStyle->SetPalette(1); |
068e2c00 | 1031 | } |
1032 | ||
0b30040d | 1033 | //________________________________________________________ |
1034 | void AliTRDresolution::DrawSigma(TH2 *h2, Float_t scale, Float_t m, Float_t M, const Char_t *title) | |
1035 | { | |
1036 | // Draw error bars scaled with "scale" instead of content values | |
1037 | //use range [m,M] if limits are specified | |
1038 | ||
1039 | if(!h2) return; | |
1040 | TH2 *h2e = (TH2F*)h2->Clone(Form("%s_E", h2->GetName())); | |
1041 | h2e->SetContour(10); | |
1042 | if(M>m) h2e->GetZaxis()->SetRangeUser(m, M); | |
1043 | if(title) h2e->GetZaxis()->SetTitle(title); | |
1044 | ||
1045 | for(Int_t ix(1); ix<=h2->GetNbinsX(); ix++){ | |
1046 | for(Int_t iy(1); iy<=h2->GetNbinsY(); iy++){ | |
1047 | if(h2->GetBinContent(ix, iy)<-100.) continue; | |
1048 | Float_t v(scale*h2->GetBinError(ix, iy)); | |
1049 | if(M>m && v<m) v=m+TMath::Abs((M-m)*1.e-3); | |
1050 | h2e->SetBinContent(ix, iy, v); | |
1051 | } | |
1052 | } | |
1053 | h2e->Draw("colz"); | |
1054 | } | |
1055 | ||
068e2c00 | 1056 | //________________________________________________________ |
1057 | void AliTRDresolution::GetRange(TH2 *h2, Char_t mod, Float_t *range) | |
1058 | { | |
1059 | // Returns the range of the bulk of data in histogram h2. Removes outliers. | |
1060 | // The "range" vector should be initialized with 2 elements | |
1061 | // Option "mod" can be any of | |
1062 | // - 0 : gaussian like distribution | |
1063 | // - 1 : tailed distribution | |
1064 | ||
1065 | Int_t nx(h2->GetNbinsX()) | |
1066 | , ny(h2->GetNbinsY()) | |
1067 | , n(nx*ny); | |
1068 | Double_t *data=new Double_t[n]; | |
1069 | for(Int_t ix(1), in(0); ix<=nx; ix++){ | |
1070 | for(Int_t iy(1); iy<=ny; iy++) | |
1071 | data[in++] = h2->GetBinContent(ix, iy); | |
00a3f7a4 | 1072 | } |
068e2c00 | 1073 | Double_t mean, sigm; |
1074 | AliMathBase::EvaluateUni(n, data, mean, sigm, Int_t(n*.8)); | |
1075 | ||
1076 | range[0]=mean-3.*sigm; range[1]=mean+3.*sigm; | |
1077 | if(mod==1) range[0]=TMath::Max(Float_t(1.e-3), range[0]); | |
1078 | AliDebug(2, Form("h[%s] range0[%f %f]", h2->GetName(), range[0], range[1])); | |
1079 | TH1S h1("h1SF0", "", 100, range[0], range[1]); | |
1080 | h1.FillN(n,data,0); | |
1081 | delete [] data; | |
1082 | ||
1083 | switch(mod){ | |
1084 | case 0:// gaussian distribution | |
1085 | { | |
1086 | TF1 fg("fg", "gaus", mean-3.*sigm, mean+3.*sigm); | |
1087 | h1.Fit(&fg, "QN"); | |
1088 | mean=fg.GetParameter(1); sigm=fg.GetParameter(2); | |
1089 | range[0] = mean-2.5*sigm;range[1] = mean+2.5*sigm; | |
1090 | AliDebug(2, Form(" rangeG[%f %f]", range[0], range[1])); | |
1091 | break; | |
00a3f7a4 | 1092 | } |
068e2c00 | 1093 | case 1:// tailed distribution |
1094 | { | |
1095 | Int_t bmax(h1.GetMaximumBin()); | |
1096 | Int_t jBinMin(1), jBinMax(100); | |
1097 | for(Int_t ibin(bmax); ibin--;){ | |
61f6b45e | 1098 | if(h1.GetBinContent(ibin)<1.){ |
068e2c00 | 1099 | jBinMin=ibin; break; |
1100 | } | |
1101 | } | |
1102 | for(Int_t ibin(bmax); ibin++;){ | |
61f6b45e | 1103 | if(h1.GetBinContent(ibin)<1.){ |
068e2c00 | 1104 | jBinMax=ibin; break; |
1105 | } | |
00a3f7a4 | 1106 | } |
068e2c00 | 1107 | range[0]=h1.GetBinCenter(jBinMin); range[1]=h1.GetBinCenter(jBinMax); |
1108 | AliDebug(2, Form(" rangeT[%f %f]", range[0], range[1])); | |
1109 | break; | |
1110 | } | |
00a3f7a4 | 1111 | } |
00a3f7a4 | 1112 | |
068e2c00 | 1113 | return; |
1114 | } | |
1115 | ||
3ceb45ae | 1116 | |
068e2c00 | 1117 | //________________________________________________________ |
f429b017 | 1118 | Bool_t AliTRDresolution::MakeProjectionCluster(Bool_t mc) |
068e2c00 | 1119 | { |
3ceb45ae | 1120 | // Analyse cluster |
3ed01fbe | 1121 | const Int_t kNcontours(9); |
1122 | const Int_t kNstat(300); | |
f429b017 | 1123 | Int_t cidx=mc?kMCcluster:kCluster; |
3ceb45ae | 1124 | if(fProj && fProj->At(cidx)) return kTRUE; |
1125 | if(!fContainer){ | |
1126 | AliError("Missing data container."); | |
1127 | return kFALSE; | |
1128 | } | |
1129 | THnSparse *H(NULL); | |
1130 | if(!(H = (THnSparse*)fContainer->At(cidx))){ | |
1131 | AliError(Form("Missing/Wrong data @ %d.", cidx)); | |
1132 | return kFALSE; | |
1133 | } | |
3ed01fbe | 1134 | Int_t ndim(H->GetNdimensions()); |
1135 | Int_t coord[kNdim]; memset(coord, 0, sizeof(Int_t) * kNdim); Double_t v = 0.; | |
1136 | TAxis *aa[kNdim], *as(NULL), *apt(NULL); memset(aa, 0, sizeof(TAxis*) * kNdim); | |
1137 | for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id); | |
1138 | if(ndim > kPt) apt = H->GetAxis(kPt); | |
1139 | if(ndim > kSpeciesChgRC) as = H->GetAxis(kSpeciesChgRC); | |
1140 | // build list of projections | |
1141 | const Int_t nsel(12), npsel(5); | |
1142 | // define rebinning strategy | |
1143 | const Int_t nEtaPhi(4); Int_t rebinEtaPhiX[nEtaPhi] = {1, 2, 5, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 1, 5}; | |
0b366e97 | 1144 | AliTRDresolutionProjection hp[kClNproj], *php[nsel][npsel]; memset(php, 0, nsel*npsel*sizeof(AliTRDresolutionProjection*)); |
3ed01fbe | 1145 | Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t)); |
3ceb45ae | 1146 | for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){ |
3ed01fbe | 1147 | isel++; // new selection |
f429b017 | 1148 | hp[ih].Build(Form("H%sClY%d", mc?"MC":"", ily), Form("Clusters :: r-#phi residuals ly%d", ily), kEta, kPhi, kYrez, aa); |
3ed01fbe | 1149 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); |
1150 | php[isel][np[isel]++] = &hp[ih++]; | |
f429b017 | 1151 | hp[ih].Build(Form("H%sClYn%d", mc?"MC":"", ily), Form("Clusters[-]:: r-#phi residuals ly%d", ily), kEta, kPhi, kYrez, aa); |
3ed01fbe | 1152 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); |
1153 | php[isel][np[isel]++] = &hp[ih++]; | |
f429b017 | 1154 | hp[ih].Build(Form("H%sClQn%d", mc?"MC":"", ily), Form("Clusters[-]:: Charge distribution ly%d", ily), kEta, kPhi, kSpeciesChgRC, aa); |
3ed01fbe | 1155 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); |
0b30040d | 1156 | hp[ih].SetShowRange(20., 40.); |
3ed01fbe | 1157 | php[isel][np[isel]++] = &hp[ih++]; |
f429b017 | 1158 | hp[ih].Build(Form("H%sClYXTCn%d", mc?"MC":"", ily), Form("Clusters[-]:: r-#phi(x,TC) residuals ly%d", ily), kPrez, kZrez, kYrez, aa); |
3ed01fbe | 1159 | // hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); |
1160 | php[isel][np[isel]++] = &hp[ih++]; | |
f429b017 | 1161 | hp[ih].Build(Form("H%sClYXPh%d", mc?"MC":"", ily), Form("Clusters :: r-#phi(x,#Phi) residuals ly%d", ily), kPrez, kPt, kYrez, aa); |
3ed01fbe | 1162 | // hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); |
1163 | php[isel][np[isel]++] = &hp[ih++]; | |
1164 | isel++; // new selection | |
1165 | php[isel][np[isel]++] = &hp[ih-5]; // relink HClY | |
f429b017 | 1166 | hp[ih].Build(Form("H%sClYp%d", mc?"MC":"", ily), Form("Clusters[+]:: r-#phi residuals ly%d", ily), kEta, kPhi, kYrez, aa); |
3ed01fbe | 1167 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); |
1168 | php[isel][np[isel]++] = &hp[ih++]; | |
f429b017 | 1169 | hp[ih].Build(Form("H%sClQp%d", mc?"MC":"", ily), Form("Clusters[+]:: Charge distribution ly%d", ily), kEta, kPhi, kSpeciesChgRC, aa); |
0b30040d | 1170 | hp[ih].SetShowRange(20., 40.); |
3ed01fbe | 1171 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); |
1172 | php[isel][np[isel]++] = &hp[ih++]; | |
f429b017 | 1173 | hp[ih].Build(Form("H%sClYXTCp%d", mc?"MC":"", ily), Form("Clusters[+]:: r-#phi(x,TC) residuals ly%d", ily), kPrez, kZrez, kYrez, aa); |
3ed01fbe | 1174 | // hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); |
1175 | php[isel][np[isel]++] = &hp[ih++]; | |
1176 | php[isel][np[isel]++] = &hp[ih-4]; // relink HClYXPh | |
3ceb45ae | 1177 | } |
3ed01fbe | 1178 | |
1179 | Int_t ly(0), ch(0), rcBin(as?as->FindBin(0.):-1), chBin(apt?apt->FindBin(0.):-1); | |
3ceb45ae | 1180 | for (Long64_t ib(0); ib < H->GetNbins(); ib++) { |
3ed01fbe | 1181 | v = H->GetBinContent(ib, coord); if(v<1.) continue; |
1182 | ly = coord[kBC]-1; | |
1183 | // RC selection | |
1184 | if(rcBin>0 && coord[kSpeciesChgRC] == rcBin) continue; | |
1185 | ||
1186 | // charge selection | |
1187 | ch = 0; // [-] track | |
1188 | if(chBin>0 && coord[kPt] > chBin) ch = 1; // [+] track | |
1189 | ||
1190 | isel = ly*2+ch; | |
1191 | for(Int_t jh(0); jh<np[isel]; jh++) php[isel][jh]->Increment(coord, v); | |
3ceb45ae | 1192 | } |
3ceb45ae | 1193 | TObjArray *arr(NULL); |
0b366e97 | 1194 | fProj->AddAt(arr = new TObjArray(kClNproj), cidx); |
3ed01fbe | 1195 | |
1196 | TH2 *h2(NULL); | |
1197 | for(; ih--; ){ | |
09c85be5 | 1198 | if(!hp[ih].fH) continue; |
3ed01fbe | 1199 | Int_t mid(1), nstat(kNstat); |
0b366e97 | 1200 | if(strchr(hp[ih].fH->GetName(), 'Q')){ mid=2; nstat=300;} |
3ed01fbe | 1201 | if(!(h2 = hp[ih].Projection2D(nstat, kNcontours, mid))) continue; |
1202 | arr->AddAt(h2, ih); | |
00a3f7a4 | 1203 | } |
3ed01fbe | 1204 | |
3ceb45ae | 1205 | return kTRUE; |
00a3f7a4 | 1206 | } |
1207 | ||
3ceb45ae | 1208 | //________________________________________________________ |
f429b017 | 1209 | Bool_t AliTRDresolution::MakeProjectionTracklet(Bool_t mc) |
3ceb45ae | 1210 | { |
1211 | // Analyse tracklet | |
3ed01fbe | 1212 | const Int_t kNcontours(9); |
566c3d46 | 1213 | const Int_t kNstat(30); |
f429b017 | 1214 | Int_t cidx=mc?kMCtracklet:kTracklet; |
3ceb45ae | 1215 | if(fProj && fProj->At(cidx)) return kTRUE; |
1216 | if(!fContainer){ | |
1217 | AliError("Missing data container."); | |
1218 | return kFALSE; | |
1219 | } | |
1220 | THnSparse *H(NULL); | |
1221 | if(!(H = (THnSparse*)fContainer->At(cidx))){ | |
3ed01fbe | 1222 | AliError(Form("Missing/Wrong data @ %d.", cidx)); |
3ceb45ae | 1223 | return kFALSE; |
1224 | } | |
3ed01fbe | 1225 | Int_t ndim(H->GetNdimensions()); |
1226 | Int_t coord[kNdim+1]; memset(coord, 0, sizeof(Int_t) * (kNdim+1)); Double_t v = 0.; | |
566c3d46 | 1227 | TAxis *aa[kNdim+1], *as(NULL), *ap(NULL); memset(aa, 0, sizeof(TAxis*) * (kNdim+1)); |
3ed01fbe | 1228 | for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id); |
1229 | if(ndim > kSpeciesChgRC) as = H->GetAxis(kSpeciesChgRC); | |
566c3d46 | 1230 | if(ndim > kPt) ap = H->GetAxis(kPt); |
3ed01fbe | 1231 | // build list of projections |
566c3d46 | 1232 | const Int_t nsel(54), npsel(4); |
3ed01fbe | 1233 | // define rebinning strategy |
1234 | const Int_t nEtaPhi(4); Int_t rebinEtaPhiX[nEtaPhi] = {1, 2, 5, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 1, 5}; | |
0b366e97 | 1235 | AliTRDresolutionProjection hp[kTrkltNproj], *php[nsel][npsel]; memset(php, 0, nsel*npsel*sizeof(AliTRDresolutionProjection*)); |
3ed01fbe | 1236 | Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t)); |
566c3d46 | 1237 | const Char_t chName[kNcharge] = {'n', 'p'};const Char_t chSgn[kNcharge] = {'-', '+'}; |
1238 | const Char_t ptName[kNpt] = {'l', 'i', 'h'}; | |
1239 | const Char_t *ptCut[kNpt] = {"p_{t}[GeV/c]<0.8", "0.8<=p_{t}[GeV/c]<1.5", "p_{t}[GeV/c]>=1.5"}; | |
3ed01fbe | 1240 | for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){ |
566c3d46 | 1241 | for(Int_t ipt(0); ipt<kNpt; ipt++){ |
1242 | for(Int_t ich(0); ich<kNcharge; ich++){ | |
1243 | isel++; // new selection | |
1244 | hp[ih].Build(Form("H%sTrkltY%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], ily), | |
1245 | Form("Tracklets[%c]:: #Deltay{%s} Ly[%d]", chSgn[ich], ptCut[ipt], ily), | |
1246 | kEta, kPhi, kYrez, aa); | |
1247 | //hp[ih].SetShowRange(-0.1,0.1); | |
1248 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); | |
1249 | php[isel][np[isel]++] = &hp[ih++]; | |
1250 | hp[ih].Build(Form("H%sTrkltPh%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], ily), | |
1251 | Form("Tracklets[%c]:: #Delta#phi{%s} Ly[%d]", chSgn[ich], ptCut[ipt], ily), | |
1252 | kEta, kPhi, kPrez, aa); | |
1253 | //hp[ih].SetShowRange(-0.5,0.5); | |
1254 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); | |
1255 | php[isel][np[isel]++] = &hp[ih++]; | |
1256 | hp[ih].Build(Form("H%sTrkltQ%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], ily), | |
1257 | Form("Tracklets[%c]:: dQdl{%s} Ly[%d]", chSgn[ich], ptCut[ipt], ily), | |
1258 | kEta, kPhi, kNdim, aa); | |
1259 | hp[ih].SetShowRange(700.,1100.); | |
1260 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); | |
1261 | php[isel][np[isel]++] = &hp[ih++]; | |
1262 | } | |
1263 | isel++; // new selection | |
1264 | hp[ih].Build(Form("H%sTrkltRCZ%c%d", mc?"MC":"", ptName[ipt], ily), | |
1265 | Form("Tracklets[RC]:: #Deltaz{%s} Ly[%d]", ptCut[ipt], ily), | |
1266 | kEta, kPhi, kZrez, aa); | |
1267 | // hp[ih].SetShowRange(-0.1,0.1); | |
1268 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); | |
1269 | php[isel][np[isel]++] = &hp[ih++]; | |
1270 | hp[ih].Build(Form("H%sTrkltRCY%c%d", mc?"MC":"", ptName[ipt], ily), | |
1271 | Form("Tracklets[RC]:: #Deltay{%s} Ly[%d]", ptCut[ipt], ily), | |
1272 | kEta, kPhi, kYrez, aa); | |
1273 | //hp[ih].SetShowRange(-0.1,0.1); | |
1274 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); | |
1275 | php[isel][np[isel]++] = &hp[ih++]; | |
1276 | hp[ih].Build(Form("H%sTrkltRCPh%c%d", mc?"MC":"", ptName[ipt], ily), | |
1277 | Form("Tracklets[RC]:: #Delta#phi{%s} Ly[%d]", ptCut[ipt], ily), | |
1278 | kEta, kPhi, kPrez, aa); | |
1279 | //hp[ih].SetShowRange(-0.1,0.1); | |
1280 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); | |
1281 | php[isel][np[isel]++] = &hp[ih++]; | |
1282 | hp[ih].Build(Form("H%sTrkltRCQ%c%d", mc?"MC":"", ptName[ipt], ily), | |
1283 | Form("Tracklets[RC]:: dQdl{%s} Ly[%d]", ptCut[ipt], ily), | |
1284 | kEta, kPhi, kNdim, aa); | |
1285 | //hp[ih].SetShowRange(-0.1,0.1); | |
1286 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); | |
1287 | php[isel][np[isel]++] = &hp[ih++]; | |
1288 | } | |
3ed01fbe | 1289 | } |
1290 | ||
566c3d46 | 1291 | Int_t ly(0), ch(0), rcBin(as?as->FindBin(0.):-1), pt(0); |
3ed01fbe | 1292 | for (Long64_t ib(0); ib < H->GetNbins(); ib++) { |
1293 | v = H->GetBinContent(ib, coord); | |
1294 | if(v<1.) continue; | |
1295 | ly = coord[kBC]-1; // layer selection | |
1296 | // charge selection | |
1297 | ch = 0; // [-] track | |
1298 | if(rcBin>0){ // debug mode in which species are also saved | |
1299 | if(coord[kSpeciesChgRC] > rcBin) ch = 1; // [+] track | |
1300 | else if(coord[kSpeciesChgRC] == rcBin) ch = 2; // [RC] track | |
1301 | } | |
566c3d46 | 1302 | // pt selection |
1303 | pt = 0; // low pt | |
1304 | if(ap) pt = coord[kPt]-1; | |
1305 | // global selection | |
1306 | isel = ly*9+pt*3+ch; | |
3ed01fbe | 1307 | for(Int_t jh(0); jh<np[isel]; jh++) php[isel][jh]->Increment(coord, v); |
1308 | } | |
3ed01fbe | 1309 | TObjArray *arr(NULL); |
0b366e97 | 1310 | fProj->AddAt(arr = new TObjArray(kTrkltNproj), cidx); |
3ed01fbe | 1311 | |
566c3d46 | 1312 | TH2 *h2(NULL); Int_t jh(0); |
3ed01fbe | 1313 | for(; ih--; ){ |
09c85be5 | 1314 | if(!hp[ih].fH) continue; |
3ed01fbe | 1315 | Int_t mid(0), nstat(kNstat); |
0b366e97 | 1316 | if(strchr(hp[ih].fH->GetName(), 'Q')){ mid=2; nstat=200;} |
3ed01fbe | 1317 | if(!(h2 = hp[ih].Projection2D(nstat, kNcontours, mid))) continue; |
566c3d46 | 1318 | arr->AddAt(h2, jh++); |
3ed01fbe | 1319 | } |
566c3d46 | 1320 | // build combined performance plots |
1321 | Int_t iproj(0); | |
1322 | for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){ | |
1323 | /*!dy negative tracks all momenta*/ | |
1324 | iproj = ily*30; | |
1325 | hp[iproj]+=hp[iproj+10]; hp[iproj]+=hp[iproj+20]; | |
1326 | hp[iproj].fH->SetNameTitle(Form("H%sTrkltYn%d", mc?"MC":"", ily), Form("Tracklet[-]:: #Deltay Ly[%d]", ily)); | |
1327 | if((h2 = hp[iproj].Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++); | |
1328 | /*!dy positive tracks all momenta*/ | |
1329 | iproj = ily*30+3; | |
1330 | hp[iproj]+=hp[iproj+10]; hp[iproj]+=hp[iproj+20]; | |
1331 | hp[iproj].fH->SetNameTitle(Form("H%sTrkltYp%d", mc?"MC":"", ily), Form("Tracklet[+]:: #Deltay Ly[%d]", ily)); | |
1332 | if((h2 = hp[iproj].Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++); | |
1333 | /*!dy all tracks all momenta*/ | |
1334 | iproj = ily*30; | |
1335 | hp[iproj]+=hp[iproj+3];hp[iproj]. | |
1336 | fH->SetNameTitle(Form("H%sTrkltY%d", mc?"MC":"", ily), Form("Tracklet :: #Deltay Ly[%d]", ily)); | |
1337 | if((h2 = hp[iproj].Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++); | |
1338 | /*!dy all tracks high momenta*/ | |
1339 | iproj = ily*30+20; | |
1340 | hp[iproj]+=hp[iproj+3]; | |
1341 | hp[iproj].fH->SetNameTitle(Form("H%sTrkltYh%d", mc?"MC":"", ily), Form("Tracklet :: #Deltay{%s} Ly[%d]", ptCut[2], ily)); | |
1342 | if((h2 = hp[iproj].Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++); | |
1343 | /*!dphi negative tracks all momenta*/ | |
1344 | iproj = ily*30+1; | |
1345 | if(hp[iproj].fH){ | |
1346 | hp[iproj]+=hp[iproj+10]; hp[iproj]+=hp[iproj+20]; | |
1347 | hp[iproj].fH->SetNameTitle(Form("H%sTrkltPhn%d", mc?"MC":"", ily), Form("Tracklet[-]:: #Delta#phi Ly[%d]", ily)); | |
1348 | if((h2 = hp[iproj].Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++); | |
1349 | } | |
1350 | /*!dphi positive tracks all momenta*/ | |
1351 | iproj = ily*30+4; | |
1352 | if(hp[iproj].fH){ | |
1353 | hp[iproj]+=hp[iproj+10]; hp[iproj]+=hp[iproj+20]; | |
1354 | hp[iproj].fH->SetNameTitle(Form("H%sTrkltPhp%d", mc?"MC":"", ily), Form("Tracklet[+]:: #Delta#phi Ly[%d]", ily)); | |
1355 | if((h2 = hp[iproj].Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++); | |
1356 | } | |
1357 | /*!dQdl negative tracks all momenta*/ | |
1358 | iproj = ily*30+2; | |
1359 | if(hp[iproj].fH){ | |
1360 | hp[iproj]+=hp[iproj+10]; hp[iproj]+=hp[iproj+20]; | |
1361 | hp[iproj].fH->SetNameTitle(Form("H%sTrkltQn%d", mc?"MC":"", ily), Form("Tracklet[-]:: dQdl Ly[%d]", ily)); | |
1362 | if((h2 = hp[iproj].Projection2D(kNstat, kNcontours, 2))) arr->AddAt(h2, jh++); | |
1363 | } | |
1364 | /*!dQdl positive tracks all momenta*/ | |
1365 | iproj = ily*30+5; | |
1366 | if(hp[iproj].fH){ | |
1367 | hp[iproj]+=hp[iproj+10]; hp[iproj]+=hp[iproj+20]; | |
1368 | hp[iproj].fH->SetNameTitle(Form("H%sTrkltQp%d", mc?"MC":"", ily), Form("Tracklet[+]:: dQdl Ly[%d]", ily)); | |
1369 | if((h2 = hp[iproj].Projection2D(kNstat, kNcontours, 2))) arr->AddAt(h2, jh++); | |
1370 | } | |
1371 | /*!dQdl all tracks all momenta*/ | |
1372 | iproj = ily*30+2; | |
1373 | if(hp[iproj].fH){ | |
1374 | hp[iproj]+=hp[iproj+3];hp[iproj]. | |
1375 | fH->SetNameTitle(Form("H%sTrkltQ%d", mc?"MC":"", ily), Form("Tracklet :: dQdl Ly[%d]", ily)); | |
1376 | if((h2 = hp[iproj].Projection2D(kNstat, kNcontours, 2))) arr->AddAt(h2, jh++); | |
1377 | } | |
1378 | /*!dz[RC] tracks all momenta*/ | |
1379 | iproj = ily*30+6; | |
1380 | if(hp[iproj].fH){ | |
1381 | hp[iproj]+=hp[iproj+10]; hp[iproj]+=hp[iproj+20]; | |
1382 | hp[iproj].fH->SetNameTitle(Form("H%sTrkltRCZ%d", mc?"MC":"", ily), Form("Tracklet[RC]:: #Deltaz Ly[%d]", ily)); | |
1383 | if((h2 = hp[iproj].Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++); | |
1384 | } | |
1385 | /*!dy[RC] tracks all momenta*/ | |
1386 | iproj = ily*30+7; | |
1387 | hp[iproj]+=hp[iproj+10]; hp[iproj]+=hp[iproj+20]; | |
1388 | hp[iproj].fH->SetNameTitle(Form("H%sTrkltRCY%d", mc?"MC":"", ily), Form("Tracklet[RC]:: #Deltay Ly[%d]", ily)); | |
1389 | if((h2 = hp[iproj].Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++); | |
1390 | /*!dphi[RC] tracks all momenta*/ | |
1391 | iproj = ily*30+8; | |
1392 | if(hp[iproj].fH){ | |
1393 | hp[iproj]+=hp[iproj+10]; hp[iproj]+=hp[iproj+20]; | |
1394 | hp[iproj].fH->SetNameTitle(Form("H%sTrkltRCPh%d", mc?"MC":"", ily), Form("Tracklet[RC]:: #Delta#phi Ly[%d]", ily)); | |
1395 | if((h2 = hp[iproj].Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++); | |
1396 | } | |
1397 | /*!dQdl[RC] tracks all momenta*/ | |
1398 | iproj = ily*30+9; | |
1399 | if(hp[iproj].fH){ | |
1400 | hp[iproj]+=hp[iproj+10]; hp[iproj]+=hp[iproj+20]; | |
1401 | hp[iproj].fH->SetNameTitle(Form("H%sTrkltRCQ%d", mc?"MC":"", ily), Form("Tracklet[RC]:: dQdl Ly[%d]", ily)); | |
1402 | if((h2 = hp[iproj].Projection2D(kNstat, kNcontours, 2))) arr->AddAt(h2, jh++); | |
1403 | } | |
1404 | } | |
1405 | ||
3ceb45ae | 1406 | return kTRUE; |
1407 | } | |
068e2c00 | 1408 | |
00a3f7a4 | 1409 | //________________________________________________________ |
f429b017 | 1410 | Bool_t AliTRDresolution::MakeProjectionTrackIn(Bool_t mc) |
1ee39b3a | 1411 | { |
3ceb45ae | 1412 | // Analyse track in |
61f6b45e | 1413 | |
3ed01fbe | 1414 | const Int_t kNcontours(9); |
1415 | const Int_t kNstat(30); | |
f429b017 | 1416 | Int_t cidx=mc?kMCtrackIn:kTrackIn; |
3ceb45ae | 1417 | if(fProj && fProj->At(cidx)) return kTRUE; |
1418 | if(!fContainer){ | |
1419 | AliError("Missing data container."); | |
1420 | return kFALSE; | |
1421 | } | |
1422 | THnSparse *H(NULL); | |
1423 | if(!(H = (THnSparse*)fContainer->At(cidx))){ | |
1424 | AliError(Form("Missing/Wrong data @ %d.", Int_t(cidx))); | |
1ee39b3a | 1425 | return kFALSE; |
1426 | } | |
61f6b45e | 1427 | |
3ed01fbe | 1428 | Int_t coord[kNdim]; memset(coord, 0, sizeof(Int_t) * kNdim); Double_t v = 0.; |
1429 | Int_t ndim(H->GetNdimensions()); | |
566c3d46 | 1430 | TAxis *aa[kNdim+1], *as(NULL), *ap(NULL), *abf(NULL); memset(aa, 0, sizeof(TAxis*) * (kNdim+1)); |
3ed01fbe | 1431 | for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id); |
1432 | if(ndim > kSpeciesChgRC) as = H->GetAxis(kSpeciesChgRC); | |
35983729 | 1433 | if(ndim > kPt) ap = H->GetAxis(kPt); |
566c3d46 | 1434 | if(ndim > (kNdim+2)) abf = H->GetAxis(kNdim+2); |
3ed01fbe | 1435 | // build list of projections |
566c3d46 | 1436 | const Int_t nsel(16), npsel(4); |
3ed01fbe | 1437 | // define rebinning strategy |
1438 | const Int_t nEtaPhi(4); Int_t rebinEtaPhiX[nEtaPhi] = {1, 2, 5, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 1, 5}; | |
35983729 | 1439 | AliTRDresolutionProjection hp[kMCTrkInNproj], *php[nsel][npsel]; memset(php, 0, nsel*npsel*sizeof(AliTRDresolutionProjection*)); |
3ed01fbe | 1440 | Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t)); |
1441 | // define list of projections | |
35983729 | 1442 | isel++; // negative low pt tracks |
1443 | hp[ih].Build(Form("H%sTrkInYnl", mc?"MC":""), "TrackIn[-]:: #Deltay{p_{t}[GeV/c]<0.8}", kEta, kPhi, kYrez, aa); | |
1444 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); | |
1445 | php[isel][np[isel]++] = &hp[ih++]; | |
1446 | hp[ih].Build(Form("H%sTrkInPhnl", mc?"MC":""), "TrackIn[-]:: #Delta#phi{p_{t}[GeV/c]<0.8}", kEta, kPhi, kPrez, aa); | |
3ed01fbe | 1447 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); |
1448 | php[isel][np[isel]++] = &hp[ih++]; | |
35983729 | 1449 | hp[ih].Build(Form("H%sTrkInXnl", mc?"MC":""), "TrackIn[-]:: #Deltax{p_{t}[GeV/c]<0.8}", kEta, kPhi, kNdim+1, aa); |
3ed01fbe | 1450 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); |
1451 | php[isel][np[isel]++] = &hp[ih++]; | |
35983729 | 1452 | isel++; // negative intermediate pt tracks |
1453 | hp[ih].Build(Form("H%sTrkInYni", mc?"MC":""), "TrackIn[-]:: #Deltay{0.8<=p_{t}[GeV/c]<1.5}", kEta, kPhi, kYrez, aa); | |
3ed01fbe | 1454 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); |
1455 | php[isel][np[isel]++] = &hp[ih++]; | |
35983729 | 1456 | hp[ih].Build(Form("H%sTrkInPhni", mc?"MC":""), "TrackIn[-]:: #Delta#phi{0.8<=p_{t}[GeV/c]<1.5}", kEta, kPhi, kPrez, aa); |
1457 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); | |
3ed01fbe | 1458 | php[isel][np[isel]++] = &hp[ih++]; |
35983729 | 1459 | hp[ih].Build(Form("H%sTrkInXni", mc?"MC":""), "TrackIn[-]:: #Deltax{0.8<=p_{t}[GeV/c]<1.5}", kEta, kPhi, kNdim+1, aa); |
3ed01fbe | 1460 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); |
1461 | php[isel][np[isel]++] = &hp[ih++]; | |
35983729 | 1462 | isel++; // negative high pt tracks |
1463 | hp[ih].Build(Form("H%sTrkInYnh", mc?"MC":""), "TrackIn[-]:: #Deltay{p_{t}[GeV/c]>=1.5}", kEta, kPhi, kYrez, aa); | |
3ed01fbe | 1464 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); |
1465 | php[isel][np[isel]++] = &hp[ih++]; | |
35983729 | 1466 | hp[ih].Build(Form("H%sTrkInPhnh", mc?"MC":""), "TrackIn[-]:: #Delta#phi{p_{t}[GeV/c]>=1.5}", kEta, kPhi, kPrez, aa); |
1467 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); | |
1468 | php[isel][np[isel]++] = &hp[ih++]; | |
1469 | hp[ih].Build(Form("H%sTrkInXnh", mc?"MC":""), "TrackIn[-]:: #Deltax{p_{t}[GeV/c]>=1.5}", kEta, kPhi, kNdim+1, aa); | |
1470 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); | |
3ed01fbe | 1471 | php[isel][np[isel]++] = &hp[ih++]; |
35983729 | 1472 | isel++; // positive low pt tracks |
1473 | hp[ih].Build(Form("H%sTrkInYpl", mc?"MC":""), "TrackIn[+]:: #Deltay{p_{t}[GeV/c]<0.8}", kEta, kPhi, kYrez, aa); | |
3ed01fbe | 1474 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); |
1475 | php[isel][np[isel]++] = &hp[ih++]; | |
35983729 | 1476 | hp[ih].Build(Form("H%sTrkInPhpl", mc?"MC":""), "TrackIn[+]:: #Delta#phi{p_{t}[GeV/c]<0.8}", kEta, kPhi, kPrez, aa); |
1477 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); | |
1478 | php[isel][np[isel]++] = &hp[ih++]; | |
1479 | hp[ih].Build(Form("H%sTrkInXpl", mc?"MC":""), "TrackIn[+]:: #Deltax{p_{t}[GeV/c]<0.8}", kEta, kPhi, kNdim+1, aa); | |
1480 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); | |
1481 | php[isel][np[isel]++] = &hp[ih++]; | |
1482 | isel++; // positive intermediate pt tracks | |
1483 | hp[ih].Build(Form("H%sTrkInYpi", mc?"MC":""), "TrackIn[+]:: #Deltay{0.8<=p_{t}[GeV/c]<1.5}", kEta, kPhi, kYrez, aa); | |
1484 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); | |
1485 | php[isel][np[isel]++] = &hp[ih++]; | |
1486 | hp[ih].Build(Form("H%sTrkInPhpi", mc?"MC":""), "TrackIn[+]:: #Delta#phi{0.8<=p_{t}[GeV/c]<1.5}", kEta, kPhi, kPrez, aa); | |
1487 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); | |
1488 | php[isel][np[isel]++] = &hp[ih++]; | |
1489 | hp[ih].Build(Form("H%sTrkInXpi", mc?"MC":""), "TrackIn[+]:: #Deltax{0.8<=p_{t}[GeV/c]<1.5}", kEta, kPhi, kNdim+1, aa); | |
1490 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); | |
1491 | php[isel][np[isel]++] = &hp[ih++]; | |
1492 | isel++; // positive high pt tracks | |
1493 | hp[ih].Build(Form("H%sTrkInYph", mc?"MC":""), "TrackIn[+]:: #Deltay{p_{t}[GeV/c]>=1.5}", kEta, kPhi, kYrez, aa); | |
1494 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); | |
1495 | php[isel][np[isel]++] = &hp[ih++]; | |
1496 | hp[ih].Build(Form("H%sTrkInPhph", mc?"MC":""), "TrackIn[+]:: #Delta#phi{p_{t}[GeV/c]>=1.5}", kEta, kPhi, kPrez, aa); | |
1497 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); | |
1498 | php[isel][np[isel]++] = &hp[ih++]; | |
1499 | hp[ih].Build(Form("H%sTrkInXph", mc?"MC":""), "TrackIn[+]:: #Deltax{p_{t}[GeV/c]>=1.5}", kEta, kPhi, kNdim+1, aa); | |
1500 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); | |
1501 | php[isel][np[isel]++] = &hp[ih++]; | |
566c3d46 | 1502 | isel++; // RC tracks low pt |
1503 | hp[ih].Build(Form("H%sTrkInRCZl", mc?"MC":""), "TrackIn[RC]:: #Deltaz{p_{t}[GeV/c]<0.8}", kEta, kPhi, kZrez, aa); | |
1504 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); | |
1505 | php[isel][np[isel]++] = &hp[ih++]; | |
1506 | hp[ih].Build(Form("H%sTrkInRCYl", mc?"MC":""), "TrackIn[RC]:: #Deltay{p_{t}[GeV/c]<0.8}", kEta, kPhi, kYrez, aa); | |
1507 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); | |
1508 | php[isel][np[isel]++] = &hp[ih++]; | |
1509 | hp[ih].Build(Form("H%sTrkInRCPhl", mc?"MC":""), "TrackIn[RC]:: #Delta#phi{p_{t}[GeV/c]<0.8}", kEta, kPhi, kPrez, aa); | |
1510 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); | |
1511 | php[isel][np[isel]++] = &hp[ih++]; | |
1512 | hp[ih].Build(Form("H%sTrkInRCXl", mc?"MC":""), "TrackIn[RC]:: #Deltax{p_{t}[GeV/c]<0.8}", kEta, kPhi, kNdim+1, aa); | |
1513 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); | |
1514 | php[isel][np[isel]++] = &hp[ih++]; | |
1515 | isel++; // RC tracks intermediate pt | |
1516 | hp[ih].Build(Form("H%sTrkInRCZi", mc?"MC":""), "TrackIn[RC]:: #Deltaz{0.8<=p_{t}[GeV/c]<1.5}", kEta, kPhi, kZrez, aa); | |
1517 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); | |
1518 | php[isel][np[isel]++] = &hp[ih++]; | |
1519 | hp[ih].Build(Form("H%sTrkInRCYi", mc?"MC":""), "TrackIn[RC]:: #Deltay{0.8<=p_{t}[GeV/c]<1.5}", kEta, kPhi, kYrez, aa); | |
1520 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); | |
1521 | php[isel][np[isel]++] = &hp[ih++]; | |
1522 | hp[ih].Build(Form("H%sTrkInRCPhi", mc?"MC":""), "TrackIn[RC]:: #Delta#phi{0.8<=p_{t}[GeV/c]<1.5}", kEta, kPhi, kPrez, aa); | |
1523 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); | |
1524 | php[isel][np[isel]++] = &hp[ih++]; | |
1525 | hp[ih].Build(Form("H%sTrkInRCXi", mc?"MC":""), "TrackIn[RC]:: #Deltax{0.8<=p_{t}[GeV/c]<1.5}", kEta, kPhi, kNdim+1, aa); | |
1526 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); | |
1527 | php[isel][np[isel]++] = &hp[ih++]; | |
1528 | isel++; // RC tracks high pt | |
1529 | hp[ih].Build(Form("H%sTrkInRCZh", mc?"MC":""), "TrackIn[RC]:: #Deltaz{p_{t}[GeV/c]>=1.5}", kEta, kPhi, kZrez, aa); | |
1530 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); | |
1531 | php[isel][np[isel]++] = &hp[ih++]; | |
1532 | hp[ih].Build(Form("H%sTrkInRCYh", mc?"MC":""), "TrackIn[RC]:: #Deltay{p_{t}[GeV/c]>=1.5}", kEta, kPhi, kYrez, aa); | |
35983729 | 1533 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); |
1534 | php[isel][np[isel]++] = &hp[ih++]; | |
566c3d46 | 1535 | hp[ih].Build(Form("H%sTrkInRCPhh", mc?"MC":""), "TrackIn[RC]:: #Delta#phi{p_{t}[GeV/c]>=1.5}", kEta, kPhi, kPrez, aa); |
35983729 | 1536 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); |
1537 | php[isel][np[isel]++] = &hp[ih++]; | |
566c3d46 | 1538 | hp[ih].Build(Form("H%sTrkInRCXh", mc?"MC":""), "TrackIn[RC]:: #Deltax{p_{t}[GeV/c]>=1.5}", kEta, kPhi, kNdim+1, aa); |
1539 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); | |
1540 | php[isel][np[isel]++] = &hp[ih++]; | |
1541 | ||
35983729 | 1542 | if(mc){ |
1543 | for(Int_t is(0); is<AliPID::kSPECIES; is++){ | |
1544 | isel++; // negative MC tracks | |
1545 | hp[ih].Build(Form("HMCTrkInYn%s", AliPID::ParticleShortName(is)), Form("TrackIn[%s-]:: #Deltay", AliPID::ParticleShortName(is)), kEta, kPhi, kYrez, aa); | |
1546 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); | |
1547 | php[isel][np[isel]++] = &hp[ih++]; | |
1548 | hp[ih].Build(Form("HMCTrkInPhn%s", AliPID::ParticleShortName(is)), Form("TrackIn[%s-]:: #Delta#phi", AliPID::ParticleShortName(is)), kEta, kPhi, kPrez, aa); | |
1549 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); | |
1550 | php[isel][np[isel]++] = &hp[ih++]; | |
1551 | } | |
1552 | for(Int_t is(0); is<AliPID::kSPECIES; is++){ | |
1553 | isel++; // positive MC tracks | |
1554 | hp[ih].Build(Form("HMCTrkInYp%s", AliPID::ParticleShortName(is)), Form("TrackIn[%s+]:: #Deltay", AliPID::ParticleShortName(is)), kEta, kPhi, kYrez, aa); | |
1555 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); | |
1556 | php[isel][np[isel]++] = &hp[ih++]; | |
1557 | hp[ih].Build(Form("HMCTrkInPhp%s", AliPID::ParticleShortName(is)), Form("TrackIn[%s+]:: #Delta#phi", AliPID::ParticleShortName(is)), kEta, kPhi, kPrez, aa); | |
1558 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); | |
1559 | php[isel][np[isel]++] = &hp[ih++]; | |
1560 | } | |
1561 | } | |
3ed01fbe | 1562 | |
1563 | // fill projections | |
35983729 | 1564 | Int_t ch(0), pt(0), rcBin(as?as->FindBin(0.):-1); |
3ceb45ae | 1565 | for (Long64_t ib(0); ib < H->GetNbins(); ib++) { |
1566 | v = H->GetBinContent(ib, coord); | |
1567 | if(v<1.) continue; | |
566c3d46 | 1568 | if(fBCbinTOF>0 && coord[kBC]!=fBCbinTOF) continue; // TOF bunch cross cut |
1569 | if(fBCbinFill>0 && abf && coord[kNdim+2]!=fBCbinTOF) continue; // Fill bunch cut | |
3ed01fbe | 1570 | // charge selection |
1571 | ch = 0; // [-] track | |
1572 | if(rcBin>0){ // debug mode in which species are also saved | |
1573 | if(coord[kSpeciesChgRC] > rcBin) ch = 1; // [+] track | |
1574 | else if(coord[kSpeciesChgRC] == rcBin) ch = 2; // [RC] track | |
3ceb45ae | 1575 | } |
35983729 | 1576 | // pt selection |
1577 | pt = 0; // low pt | |
1578 | if(ap) pt = coord[kPt]-1; | |
1579 | // global selection | |
1580 | Int_t selection = ch*3+pt; | |
1581 | for(Int_t jh(0); jh<np[selection]; jh++) php[selection][jh]->Increment(coord, v); | |
1582 | if(!mc || rcBin<0 || ch==2) continue; | |
1583 | Int_t spec = Int_t(TMath::Abs(as->GetBinCenter(coord[kSpeciesChgRC])))-1; | |
1584 | selection+=(ch*AliPID::kSPECIES+spec); | |
1585 | for(Int_t jh(0); jh<np[selection]; jh++) php[selection][jh]->Increment(coord, v); | |
3ceb45ae | 1586 | } |
1587 | TObjArray *arr(NULL); | |
35983729 | 1588 | fProj->AddAt(arr = new TObjArray(mc?kMCTrkInNproj:kTrkInNproj), cidx); |
3ceb45ae | 1589 | |
566c3d46 | 1590 | TH2 *h2(NULL); Int_t jh(0); |
3ed01fbe | 1591 | for(; ih--; ){ |
09c85be5 | 1592 | if(!hp[ih].fH) continue; |
0b366e97 | 1593 | if(!(h2 = hp[ih].Projection2D(kNstat, kNcontours))) continue; |
566c3d46 | 1594 | arr->AddAt(h2, jh++); |
3ed01fbe | 1595 | } |
35983729 | 1596 | // build combined performance plots |
1597 | /*!dy negative tracks all momenta*/ | |
1598 | Int_t iproj(0); | |
566c3d46 | 1599 | hp[iproj]+=hp[iproj+3]; hp[iproj]+=hp[iproj+6]; hp[iproj].fH->SetNameTitle(Form("H%sTrkInYn", mc?"MC":""), "TrackIn[-]:: #Deltay"); |
1600 | if((h2 = hp[iproj].Projection2D(kNstat, kNcontours))) arr->AddAt(h2, jh++); | |
35983729 | 1601 | /*!dy positive tracks all momenta*/ |
1602 | iproj = 9; | |
566c3d46 | 1603 | hp[iproj]+=hp[iproj+3]; hp[iproj]+=hp[iproj+6]; hp[iproj].fH->SetNameTitle(Form("H%sTrkInYp", mc?"MC":""), "TrackIn[+]:: #Deltay"); |
1604 | if((h2 = hp[iproj].Projection2D(kNstat, kNcontours))) arr->AddAt(h2, jh++); | |
35983729 | 1605 | /*!dy all tracks all momenta*/ |
1606 | hp[0]+=hp[9];hp[0].fH->SetNameTitle(Form("H%sTrkInY", mc?"MC":""), "TrackIn :: #Deltay"); | |
566c3d46 | 1607 | if((h2 = hp[0].Projection2D(kNstat, kNcontours))) arr->AddAt(h2, jh++); |
35983729 | 1608 | /*!dy all tracks high momenta*/ |
1609 | iproj = 6; | |
1610 | hp[iproj]+=hp[iproj+9];hp[iproj].fH->SetNameTitle(Form("H%sTrkInYh", mc?"MC":""), "TrackIn :: #Deltay{p_{t}[GeV/c]>=1.5}"); | |
566c3d46 | 1611 | if((h2 = hp[iproj].Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++); |
35983729 | 1612 | /*!dx all tracks low momenta*/ |
1613 | iproj = 2; | |
d222dc79 | 1614 | if(hp[iproj].fH){ |
1615 | hp[iproj]+=hp[iproj+9];hp[iproj].fH->SetNameTitle(Form("H%sTrkInXl", mc?"MC":""), "TrackIn :: #Deltax{p_{t}[GeV/c]<0.8}"); | |
566c3d46 | 1616 | if((h2 = hp[iproj].Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++); |
d222dc79 | 1617 | } |
35983729 | 1618 | /*!dphi negative tracks all momenta*/ |
1619 | iproj =1; | |
566c3d46 | 1620 | hp[iproj]+=hp[iproj+3]; hp[iproj]+=hp[iproj+6]; hp[iproj].fH->SetNameTitle(Form("H%sTrkInPhn", mc?"MC":""), "TrackIn[-]:: #Delta#phi"); |
1621 | if((h2 = hp[iproj].Projection2D(kNstat, kNcontours))) arr->AddAt(h2, jh++); | |
35983729 | 1622 | /*!dphi positive tracks all momenta*/ |
1623 | iproj = 10; | |
566c3d46 | 1624 | hp[iproj]+=hp[iproj+3]; hp[iproj]+=hp[iproj+6]; hp[iproj].fH->SetNameTitle(Form("H%sTrkInPhp", mc?"MC":""), "TrackIn[+]:: #Delta#phi"); |
1625 | if((h2 = hp[iproj].Projection2D(kNstat, kNcontours))) arr->AddAt(h2, jh++); | |
1626 | /*!dz[RC] tracks all momenta*/ | |
35983729 | 1627 | iproj = 18; |
566c3d46 | 1628 | hp[iproj]+=hp[iproj+4];hp[iproj]+=hp[iproj+8]; hp[iproj].fH->SetNameTitle(Form("H%sTrkInRCZ", mc?"MC":""), "TrackIn[RC]:: #Deltaz"); |
1629 | if((h2 = hp[iproj].Projection2D(kNstat, kNcontours))) arr->AddAt(h2, jh++); | |
1630 | /*!dy[RC] tracks all momenta*/ | |
1631 | iproj = 19; | |
1632 | hp[iproj]+=hp[iproj+4];hp[iproj]+=hp[iproj+8]; hp[iproj].fH->SetNameTitle(Form("H%sTrkInRCY", mc?"MC":""), "TrackIn[RC]:: #Deltay"); | |
1633 | if((h2 = hp[iproj].Projection2D(kNstat, kNcontours))) arr->AddAt(h2, jh++); | |
1634 | /*!dphi[RC] tracks all momenta*/ | |
1635 | iproj = 20; | |
1636 | hp[iproj]+=hp[iproj+4];hp[iproj]+=hp[iproj+8]; hp[iproj].fH->SetNameTitle(Form("H%sTrkInRCPh", mc?"MC":""), "TrackIn[RC]:: #Delta#phi"); | |
1637 | if((h2 = hp[iproj].Projection2D(kNstat, kNcontours))) arr->AddAt(h2, jh++); | |
1638 | /*!dx[RC] tracks all momenta*/ | |
1639 | iproj = 21; | |
1640 | if(hp[iproj].fH){ | |
1641 | hp[iproj]+=hp[iproj+4];hp[iproj]+=hp[iproj+8]; hp[iproj].fH->SetNameTitle(Form("H%sTrkInRCX", mc?"MC":""), "TrackIn[RC]:: #Deltax"); | |
1642 | if((h2 = hp[iproj].Projection2D(kNstat, kNcontours))) arr->AddAt(h2, jh++); | |
1643 | } | |
3ceb45ae | 1644 | return kTRUE; |
1645 | } | |
1ee39b3a | 1646 | |
3ceb45ae | 1647 | |
f429b017 | 1648 | //________________________________________________________ |
1649 | Bool_t AliTRDresolution::MakeProjectionTrack() | |
1650 | { | |
1651 | // Analyse tracklet | |
1652 | const Int_t kNcontours(9); | |
1653 | const Int_t kNstat(100); | |
1654 | Int_t cidx(kMCtrack); | |
1655 | if(fProj && fProj->At(cidx)) return kTRUE; | |
1656 | if(!fContainer){ | |
1657 | AliError("Missing data container."); | |
1658 | return kFALSE; | |
1659 | } | |
1660 | THnSparse *H(NULL); | |
1661 | if(!(H = (THnSparse*)fContainer->At(cidx))){ | |
1662 | AliError(Form("Missing/Wrong data @ %d.", cidx)); | |
1663 | return kFALSE; | |
1664 | } | |
1665 | Int_t ndim(H->GetNdimensions()); | |
1666 | Int_t coord[kNdim+1]; memset(coord, 0, sizeof(Int_t) * (kNdim+1)); Double_t v = 0.; | |
1667 | TAxis *aa[kNdim+1], *as(NULL); memset(aa, 0, sizeof(TAxis*) * (kNdim+1)); | |
1668 | for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id); | |
1669 | if(ndim > kSpeciesChgRC) as = H->GetAxis(kSpeciesChgRC); | |
1670 | // build list of projections | |
1671 | const Int_t nsel(18), npsel(6); | |
1672 | // define rebinning strategy | |
1673 | const Int_t nEtaPhi(4); Int_t rebinEtaPhiX[nEtaPhi] = {1, 2, 5, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 1, 5}; | |
0b366e97 | 1674 | AliTRDresolutionProjection hp[kTrkNproj], *php[nsel][npsel]; memset(php, 0, nsel*npsel*sizeof(AliTRDresolutionProjection*)); |
f429b017 | 1675 | Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t)); |
1676 | for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){ | |
1677 | isel++; // new selection | |
1678 | hp[ih].Build(Form("HTrkY%d", ily), Form("Tracks :: r-#phi residuals ly%d", ily), kEta, kPhi, kYrez, aa); | |
1679 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); | |
1680 | php[isel][np[isel]++] = &hp[ih++]; | |
1681 | hp[ih].Build(Form("HTrkYn%d", ily), Form("Tracks[-]:: r-#phi residuals ly%d", ily), kEta, kPhi, kYrez, aa); | |
1682 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); | |
1683 | php[isel][np[isel]++] = &hp[ih++]; | |
1684 | hp[ih].Build(Form("HTrkPhn%d", ily), Form("Tracks[-]:: #Delta#phi residuals ly%d", ily), kEta, kPhi, kPrez, aa); | |
1685 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); | |
1686 | php[isel][np[isel]++] = &hp[ih++]; | |
1687 | hp[ih].Build(Form("HTrkPn%d", ily), Form("Tracks[-]:: Momentum distribution ly%d", ily), kEta, kPhi, kPt, aa); | |
1688 | hp[ih].SetShowRange(6.,12.); | |
1689 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); | |
1690 | php[isel][np[isel]++] = &hp[ih++]; | |
1691 | hp[ih].Build(Form("HTrkYPn%d", ily), Form("Tracks[-]:: r-#phi/p_{t} residuals ly%d", ily), kPt, kPhi, kYrez, aa); | |
1692 | php[isel][np[isel]++] = &hp[ih++]; | |
1693 | hp[ih].Build(Form("HTrkQn%d", ily), Form("Tracks[-]:: #Deltap_{t}/p_{t} ly%d", ily), kEta, kPhi, kNdim, aa); | |
1694 | //hp[ih].SetShowRange(700.,1100.); | |
1695 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); | |
1696 | php[isel][np[isel]++] = &hp[ih++]; | |
1697 | isel++; // new selection | |
1698 | php[isel][np[isel]++] = &hp[ih-6]; // relink first histo | |
1699 | hp[ih].Build(Form("HTrkYp%d", ily), Form("Tracks[+]:: r-#phi residuals ly%d", ily), kEta, kPhi, kYrez, aa); | |
1700 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); | |
1701 | php[isel][np[isel]++] = &hp[ih++]; | |
1702 | hp[ih].Build(Form("HTrkPhp%d", ily), Form("Tracks[+]:: #Delta#phi residuals ly%d", ily), kEta, kPhi, kPrez, aa); | |
1703 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); | |
1704 | php[isel][np[isel]++] = &hp[ih++]; | |
1705 | hp[ih].Build(Form("HTrkPp%d", ily), Form("Tracks[+]:: Momentum distribution ly%d", ily), kEta, kPhi, kPt, aa); | |
1706 | hp[ih].SetShowRange(6.,12.); | |
1707 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); | |
1708 | php[isel][np[isel]++] = &hp[ih++]; | |
1709 | hp[ih].Build(Form("HTrkYPp%d", ily), Form("Tracks[+]:: r-#phi/p_{t} residuals ly%d", ily), kPt, kPhi, kYrez, aa); | |
1710 | php[isel][np[isel]++] = &hp[ih++]; | |
1711 | hp[ih].Build(Form("HTrkQp%d", ily), Form("Tracks[+]:: #Deltap_{t}/p_{t} ly%d", ily), kEta, kPhi, kNdim, aa); | |
1712 | //hp[ih].SetShowRange(700.,1100.); | |
1713 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); | |
1714 | php[isel][np[isel]++] = &hp[ih++]; | |
1715 | isel++; // new selection | |
1716 | hp[ih].Build(Form("HTrkZ%d", ily), Form("Tracks[RC]:: z residuals ly%d", ily), kEta, kPhi, kZrez, aa); | |
1717 | hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY); | |
1718 | php[isel][np[isel]++] = &hp[ih++]; | |
1719 | } | |
1720 | ||
1721 | Int_t ly(0), ch(0), rcBin(as?as->FindBin(0.):-1); | |
1722 | for (Long64_t ib(0); ib < H->GetNbins(); ib++) { | |
1723 | v = H->GetBinContent(ib, coord); | |
1724 | if(v<1.) continue; | |
1725 | ly = coord[kBC]-1; // layer selection | |
1726 | // charge selection | |
1727 | ch = 0; // [-] track | |
1728 | if(rcBin>0){ // debug mode in which species are also saved | |
1729 | if(coord[kSpeciesChgRC] > rcBin) ch = 1; // [+] track | |
1730 | else if(coord[kSpeciesChgRC] == rcBin) ch = 2; // [RC] track | |
1731 | } | |
1732 | isel = ly*3+ch; | |
1733 | for(Int_t jh(0); jh<np[isel]; jh++) php[isel][jh]->Increment(coord, v); | |
1734 | } | |
f429b017 | 1735 | TObjArray *arr(NULL); |
0b366e97 | 1736 | fProj->AddAt(arr = new TObjArray(kTrkNproj), cidx); |
f429b017 | 1737 | |
1738 | TH2 *h2(NULL); | |
1739 | for(; ih--; ){ | |
1740 | if(!hp[ih].fH) continue; | |
0b366e97 | 1741 | if(!(h2 = hp[ih].Projection2D(kNstat, kNcontours))) continue; |
f429b017 | 1742 | arr->AddAt(h2, ih); |
1743 | } | |
1744 | return kTRUE; | |
1745 | } | |
3ceb45ae | 1746 | |
1747 | //________________________________________________________ | |
1748 | Bool_t AliTRDresolution::PostProcess() | |
1749 | { | |
1750 | // Fit, Project, Combine, Extract values from the containers filled during execution | |
1751 | ||
1752 | if (!fContainer) { | |
1753 | AliError("ERROR: list not available"); | |
1754 | return kFALSE; | |
1755 | } | |
35983729 | 1756 | if(!fProj){ |
1757 | AliInfo("Building array of projections ..."); | |
1758 | fProj = new TObjArray(kNclasses); fProj->SetOwner(kTRUE); | |
1759 | } | |
1ee39b3a | 1760 | |
1ee39b3a | 1761 | //PROCESS EXPERIMENTAL DISTRIBUTIONS |
1ee39b3a | 1762 | // Clusters residuals |
566c3d46 | 1763 | // if(!MakeProjectionCluster()) return kFALSE; |
6558fd69 | 1764 | fNRefFigures = 3; |
1ee39b3a | 1765 | // Tracklet residual/pulls |
3ceb45ae | 1766 | if(!MakeProjectionTracklet()) return kFALSE; |
6558fd69 | 1767 | fNRefFigures = 7; |
a310e49b | 1768 | // TRDin residual/pulls |
3ceb45ae | 1769 | if(!MakeProjectionTrackIn()) return kFALSE; |
6558fd69 | 1770 | fNRefFigures = 11; |
1ee39b3a | 1771 | |
1772 | if(!HasMCdata()) return kTRUE; | |
1ee39b3a | 1773 | //PROCESS MC RESIDUAL DISTRIBUTIONS |
1774 | ||
1775 | // CLUSTER Y RESOLUTION/PULLS | |
f429b017 | 1776 | if(!MakeProjectionCluster(kTRUE)) return kFALSE; |
6558fd69 | 1777 | fNRefFigures = 17; |
1ee39b3a | 1778 | |
1779 | // TRACKLET RESOLUTION/PULLS | |
f429b017 | 1780 | if(!MakeProjectionTracklet(kTRUE)) return kFALSE; |
d25116b6 | 1781 | fNRefFigures = 21; |
1ee39b3a | 1782 | |
1783 | // TRACK RESOLUTION/PULLS | |
f429b017 | 1784 | if(!MakeProjectionTrack()) return kFALSE; |
d25116b6 | 1785 | fNRefFigures+=16; |
92d6d80c | 1786 | |
1787 | // TRACK TRDin RESOLUTION/PULLS | |
f429b017 | 1788 | if(!MakeProjectionTrackIn(kTRUE)) return kFALSE; |
d25116b6 | 1789 | fNRefFigures+=8; |
1ee39b3a | 1790 | |
1791 | return kTRUE; | |
1792 | } | |
1793 | ||
1794 | ||
1795 | //________________________________________________________ | |
1796 | void AliTRDresolution::Terminate(Option_t *opt) | |
1797 | { | |
1798 | AliTRDrecoTask::Terminate(opt); | |
1799 | if(HasPostProcess()) PostProcess(); | |
1800 | } | |
1801 | ||
1802 | //________________________________________________________ | |
1803 | void AliTRDresolution::AdjustF1(TH1 *h, TF1 *f) | |
1804 | { | |
1805 | // Helper function to avoid duplication of code | |
1806 | // Make first guesses on the fit parameters | |
1807 | ||
1808 | // find the intial parameters of the fit !! (thanks George) | |
1809 | Int_t nbinsy = Int_t(.5*h->GetNbinsX()); | |
1810 | Double_t sum = 0.; | |
1811 | for(Int_t jbin=nbinsy-4; jbin<=nbinsy+4; jbin++) sum+=h->GetBinContent(jbin); sum/=9.; | |
1812 | f->SetParLimits(0, 0., 3.*sum); | |
1813 | f->SetParameter(0, .9*sum); | |
1814 | Double_t rms = h->GetRMS(); | |
1815 | f->SetParLimits(1, -rms, rms); | |
1816 | f->SetParameter(1, h->GetMean()); | |
1817 | ||
1818 | f->SetParLimits(2, 0., 2.*rms); | |
1819 | f->SetParameter(2, rms); | |
1820 | if(f->GetNpar() <= 4) return; | |
1821 | ||
1822 | f->SetParLimits(3, 0., sum); | |
1823 | f->SetParameter(3, .1*sum); | |
1824 | ||
1825 | f->SetParLimits(4, -.3, .3); | |
1826 | f->SetParameter(4, 0.); | |
1827 | ||
1828 | f->SetParLimits(5, 0., 1.e2); | |
1829 | f->SetParameter(5, 2.e-1); | |
1830 | } | |
1831 | ||
afca20ef | 1832 | //________________________________________________________ |
9dcc64cc | 1833 | TObjArray* AliTRDresolution::BuildMonitorContainerCluster(const char* name, Bool_t expand, Float_t range) |
afca20ef | 1834 | { |
3d2a3dff | 1835 | // Build performance histograms for AliTRDcluster.vs TRD track or MC |
afca20ef | 1836 | // - y reziduals/pulls |
3d2a3dff | 1837 | |
1838 | TObjArray *arr = new TObjArray(2); | |
afca20ef | 1839 | arr->SetName(name); arr->SetOwner(); |
1840 | TH1 *h(NULL); char hname[100], htitle[300]; | |
1841 | ||
1842 | // tracklet resolution/pull in y direction | |
7fe4e88b | 1843 | snprintf(hname, 100, "%s_%s_Y", GetNameId(), name); |
3ceb45ae | 1844 | snprintf(htitle, 300, "Y res for \"%s\" @ %s;tg(#phi);#Delta y [cm];%s", GetNameId(), name, "Detector"); |
9dcc64cc | 1845 | Float_t rr = range<0.?fDyRange:range; |
3d2a3dff | 1846 | if(!(h = (TH3S*)gROOT->FindObject(hname))){ |
3ceb45ae | 1847 | Int_t nybins=50; |
2589cf64 | 1848 | if(expand) nybins*=2; |
3d2a3dff | 1849 | h = new TH3S(hname, htitle, |
6859821f | 1850 | 48, -.48, .48, // phi |
9dcc64cc | 1851 | 60, -rr, rr, // dy |
6859821f | 1852 | nybins, -0.5, nybins-0.5);// segment |
afca20ef | 1853 | } else h->Reset(); |
1854 | arr->AddAt(h, 0); | |
7fe4e88b | 1855 | snprintf(hname, 100, "%s_%s_YZpull", GetNameId(), name); |
3ceb45ae | 1856 | snprintf(htitle, 300, "YZ pull for \"%s\" @ %s;%s;#Delta y / #sigma_{y};#Delta z / #sigma_{z}", GetNameId(), name, "Detector"); |
81979445 | 1857 | if(!(h = (TH3S*)gROOT->FindObject(hname))){ |
3ceb45ae | 1858 | h = new TH3S(hname, htitle, 540, -0.5, 540-0.5, 100, -4.5, 4.5, 100, -4.5, 4.5); |
afca20ef | 1859 | } else h->Reset(); |
1860 | arr->AddAt(h, 1); | |
1861 | ||
3d2a3dff | 1862 | return arr; |
1863 | } | |
1864 | ||
1865 | //________________________________________________________ | |
2589cf64 | 1866 | TObjArray* AliTRDresolution::BuildMonitorContainerTracklet(const char* name, Bool_t expand) |
3d2a3dff | 1867 | { |
1868 | // Build performance histograms for AliExternalTrackParam.vs TRD tracklet | |
1869 | // - y reziduals/pulls | |
1870 | // - z reziduals/pulls | |
1871 | // - phi reziduals | |
9dcc64cc | 1872 | TObjArray *arr = BuildMonitorContainerCluster(name, expand, 0.05); |
3d2a3dff | 1873 | arr->Expand(5); |
1874 | TH1 *h(NULL); char hname[100], htitle[300]; | |
1875 | ||
afca20ef | 1876 | // tracklet resolution/pull in z direction |
7fe4e88b | 1877 | snprintf(hname, 100, "%s_%s_Z", GetNameId(), name); |
9dcc64cc | 1878 | snprintf(htitle, 300, "Z res for \"%s\" @ %s;tg(#theta);#Delta z [cm]", GetNameId(), name); |
1879 | if(!(h = (TH2S*)gROOT->FindObject(hname))){ | |
1880 | h = new TH2S(hname, htitle, 50, -1., 1., 100, -.05, .05); | |
afca20ef | 1881 | } else h->Reset(); |
1882 | arr->AddAt(h, 2); | |
7fe4e88b | 1883 | snprintf(hname, 100, "%s_%s_Zpull", GetNameId(), name); |
1884 | snprintf(htitle, 300, "Z pull for \"%s\" @ %s;tg(#theta);#Delta z / #sigma_{z};row cross", GetNameId(), name); | |
81979445 | 1885 | if(!(h = (TH3S*)gROOT->FindObject(hname))){ |
1886 | h = new TH3S(hname, htitle, 50, -1., 1., 100, -5.5, 5.5, 2, -0.5, 1.5); | |
dfd7d48b | 1887 | h->GetZaxis()->SetBinLabel(1, "no RC"); |
1888 | h->GetZaxis()->SetBinLabel(2, "RC"); | |
afca20ef | 1889 | } else h->Reset(); |
1890 | arr->AddAt(h, 3); | |
1891 | ||
1892 | // tracklet to track phi resolution | |
7fe4e88b | 1893 | snprintf(hname, 100, "%s_%s_PHI", GetNameId(), name); |
3ceb45ae | 1894 | snprintf(htitle, 300, "#Phi res for \"%s\" @ %s;tg(#phi);#Delta #phi [rad];%s", GetNameId(), name, "Detector"); |
1895 | Int_t nsgms=540; | |
9dcc64cc | 1896 | if(!(h = (TH3S*)gROOT->FindObject(hname))){ |
1897 | h = new TH3S(hname, htitle, 48, -.48, .48, 100, -.5, .5, nsgms, -0.5, nsgms-0.5); | |
afca20ef | 1898 | } else h->Reset(); |
1899 | arr->AddAt(h, 4); | |
1900 | ||
1901 | return arr; | |
1902 | } | |
1903 | ||
1904 | //________________________________________________________ | |
1905 | TObjArray* AliTRDresolution::BuildMonitorContainerTrack(const char* name) | |
1906 | { | |
1907 | // Build performance histograms for AliExternalTrackParam.vs MC | |
1908 | // - y resolution/pulls | |
1909 | // - z resolution/pulls | |
1910 | // - phi resolution, snp pulls | |
1911 | // - theta resolution, tgl pulls | |
1912 | // - pt resolution, 1/pt pulls, p resolution | |
1913 | ||
afca20ef | 1914 | TObjArray *arr = BuildMonitorContainerTracklet(name); |
1915 | arr->Expand(11); | |
3d2a3dff | 1916 | TH1 *h(NULL); char hname[100], htitle[300]; |
395d3507 | 1917 | //TAxis *ax(NULL); |
3d2a3dff | 1918 | |
afca20ef | 1919 | // snp pulls |
7fe4e88b | 1920 | snprintf(hname, 100, "%s_%s_SNPpull", GetNameId(), name); |
1921 | snprintf(htitle, 300, "SNP pull for \"%s\" @ %s;tg(#phi);#Delta snp / #sigma_{snp};entries", GetNameId(), name); | |
afca20ef | 1922 | if(!(h = (TH2I*)gROOT->FindObject(hname))){ |
1923 | h = new TH2I(hname, htitle, 60, -.3, .3, 100, -4.5, 4.5); | |
1924 | } else h->Reset(); | |
1925 | arr->AddAt(h, 5); | |
1926 | ||
1927 | // theta resolution | |
7fe4e88b | 1928 | snprintf(hname, 100, "%s_%s_THT", GetNameId(), name); |
1929 | snprintf(htitle, 300, "#Theta res for \"%s\" @ %s;tg(#theta);#Delta #theta [rad];entries", GetNameId(), name); | |
afca20ef | 1930 | if(!(h = (TH2I*)gROOT->FindObject(hname))){ |
1931 | h = new TH2I(hname, htitle, 100, -1., 1., 100, -5e-3, 5e-3); | |
1932 | } else h->Reset(); | |
1933 | arr->AddAt(h, 6); | |
1934 | // tgl pulls | |
7fe4e88b | 1935 | snprintf(hname, 100, "%s_%s_TGLpull", GetNameId(), name); |
1936 | snprintf(htitle, 300, "TGL pull for \"%s\" @ %s;tg(#theta);#Delta tgl / #sigma_{tgl};entries", GetNameId(), name); | |
afca20ef | 1937 | if(!(h = (TH2I*)gROOT->FindObject(hname))){ |
1938 | h = new TH2I(hname, htitle, 100, -1., 1., 100, -4.5, 4.5); | |
1939 | } else h->Reset(); | |
1940 | arr->AddAt(h, 7); | |
1941 | ||
afca20ef | 1942 | const Int_t kNdpt(150); |
1943 | const Int_t kNspc = 2*AliPID::kSPECIES+1; | |
61f6b45e | 1944 | Float_t lPt=0.1, lDPt=-.1, lSpc=-5.5; |
afca20ef | 1945 | Float_t binsPt[kNpt+1], binsSpc[kNspc+1], binsDPt[kNdpt+1]; |
61f6b45e | 1946 | for(Int_t i=0;i<kNpt+1; i++,lPt=TMath::Exp(i*.15)-1.) binsPt[i]=lPt; |
1947 | for(Int_t i=0; i<kNspc+1; i++,lSpc+=1.) binsSpc[i]=lSpc; | |
1948 | for(Int_t i=0; i<kNdpt+1; i++,lDPt+=2.e-3) binsDPt[i]=lDPt; | |
afca20ef | 1949 | |
1950 | // Pt resolution | |
7fe4e88b | 1951 | snprintf(hname, 100, "%s_%s_Pt", GetNameId(), name); |
1952 | snprintf(htitle, 300, "#splitline{P_{t} res for}{\"%s\" @ %s};p_{t} [GeV/c];#Delta p_{t}/p_{t}^{MC};SPECIES", GetNameId(), name); | |
afca20ef | 1953 | if(!(h = (TH3S*)gROOT->FindObject(hname))){ |
1954 | h = new TH3S(hname, htitle, | |
1955 | kNpt, binsPt, kNdpt, binsDPt, kNspc, binsSpc); | |
dbb2e0a7 | 1956 | //ax = h->GetZaxis(); |
3ceb45ae | 1957 | //for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]); |
afca20ef | 1958 | } else h->Reset(); |
1959 | arr->AddAt(h, 8); | |
1960 | // 1/Pt pulls | |
7fe4e88b | 1961 | snprintf(hname, 100, "%s_%s_1Pt", GetNameId(), name); |
1962 | snprintf(htitle, 300, "#splitline{1/P_{t} pull for}{\"%s\" @ %s};1/p_{t}^{MC} [c/GeV];#Delta(1/p_{t})/#sigma(1/p_{t});SPECIES", GetNameId(), name); | |
afca20ef | 1963 | if(!(h = (TH3S*)gROOT->FindObject(hname))){ |
1964 | h = new TH3S(hname, htitle, | |
1965 | kNpt, 0., 2., 100, -4., 4., kNspc, -5.5, 5.5); | |
395d3507 | 1966 | //ax = h->GetZaxis(); |
3ceb45ae | 1967 | //for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]); |
afca20ef | 1968 | } else h->Reset(); |
1969 | arr->AddAt(h, 9); | |
1970 | // P resolution | |
7fe4e88b | 1971 | snprintf(hname, 100, "%s_%s_P", GetNameId(), name); |
1972 | snprintf(htitle, 300, "P res for \"%s\" @ %s;p [GeV/c];#Delta p/p^{MC};SPECIES", GetNameId(), name); | |
afca20ef | 1973 | if(!(h = (TH3S*)gROOT->FindObject(hname))){ |
1974 | h = new TH3S(hname, htitle, | |
1975 | kNpt, binsPt, kNdpt, binsDPt, kNspc, binsSpc); | |
395d3507 | 1976 | //ax = h->GetZaxis(); |
3ceb45ae | 1977 | //for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]); |
afca20ef | 1978 | } else h->Reset(); |
1979 | arr->AddAt(h, 10); | |
1980 | ||
1981 | return arr; | |
1982 | } | |
1983 | ||
1984 | ||
1ee39b3a | 1985 | //________________________________________________________ |
1986 | TObjArray* AliTRDresolution::Histos() | |
1987 | { | |
1988 | // | |
1989 | // Define histograms | |
1990 | // | |
1991 | ||
1992 | if(fContainer) return fContainer; | |
1993 | ||
3ceb45ae | 1994 | fContainer = new TObjArray(kNclasses); fContainer->SetOwner(kTRUE); |
1995 | THnSparse *H(NULL); | |
1996 | const Int_t nhn(100); Char_t hn[nhn]; TString st; | |
1ee39b3a | 1997 | |
3ceb45ae | 1998 | //++++++++++++++++++++++ |
3ed01fbe | 1999 | // cluster to tracklet residuals/pulls |
3ceb45ae | 2000 | snprintf(hn, nhn, "h%s", fgPerformanceName[kCluster]); |
2001 | if(!(H = (THnSparseI*)gROOT->FindObject(hn))){ | |
3ed01fbe | 2002 | const Char_t *clTitle[kNdim] = {"layer", fgkTitle[kPhi], fgkTitle[kEta], fgkTitle[kYrez], "#Deltax [cm]", "Q</Q", "Q/angle", "#Phi [deg]"}; |
f429b017 | 2003 | const Int_t clNbins[kNdim] = {AliTRDgeometry::kNlayer, fgkNbins[kPhi], fgkNbins[kEta], fgkNbins[kYrez], 45, 10, 30, 15}; |
2004 | const Double_t clMin[kNdim] = {-0.5, fgkMin[kPhi], fgkMin[kEta], fgkMin[kYrez]/10., -.5, 0.1, -2., -45}, | |
2005 | clMax[kNdim] = {AliTRDgeometry::kNlayer-0.5, fgkMax[kPhi], fgkMax[kEta], fgkMax[kYrez]/10., 4., 2.1, 118., 45}; | |
3ceb45ae | 2006 | st = "cluster spatial&charge resolution;"; |
3ed01fbe | 2007 | // define minimum info to be saved in non debug mode |
566c3d46 | 2008 | Int_t ndim=DebugLevel()>=1?Int_t(kNdim):Int_t(kNdimCl); |
3ed01fbe | 2009 | for(Int_t idim(0); idim<ndim; idim++){ st += clTitle[idim]; st+=";";} |
2010 | H = new THnSparseI(hn, st.Data(), ndim, clNbins, clMin, clMax); | |
3ceb45ae | 2011 | } else H->Reset(); |
2012 | fContainer->AddAt(H, kCluster); | |
2013 | //++++++++++++++++++++++ | |
afca20ef | 2014 | // tracklet to TRD track |
3ceb45ae | 2015 | snprintf(hn, nhn, "h%s", fgPerformanceName[kTracklet]); |
2016 | if(!(H = (THnSparseI*)gROOT->FindObject(hn))){ | |
3ed01fbe | 2017 | Char_t *trTitle[kNdim+1]; memcpy(trTitle, fgkTitle, kNdim*sizeof(Char_t*)); |
2018 | Int_t trNbins[kNdim+1]; memcpy(trNbins, fgkNbins, kNdim*sizeof(Int_t)); | |
2019 | Double_t trMin[kNdim+1]; memcpy(trMin, fgkMin, kNdim*sizeof(Double_t)); | |
2020 | Double_t trMax[kNdim+1]; memcpy(trMax, fgkMax, kNdim*sizeof(Double_t)); | |
2021 | // set specific fields | |
f429b017 | 2022 | trMin[kYrez] = -0.45; trMax[kYrez] = 0.45; |
2023 | trMin[kPrez] = -4.5; trMax[kPrez] = 4.5; | |
2024 | trMin[kZrez] = -1.5; trMax[kZrez] = 1.5; | |
3ed01fbe | 2025 | trTitle[kBC]=StrDup("layer"); trNbins[kBC] = AliTRDgeometry::kNlayer; trMin[kBC] = -0.5; trMax[kBC] = AliTRDgeometry::kNlayer-0.5; |
2026 | trTitle[kNdim]=StrDup("dq/dl [a.u.]"); trNbins[kNdim] = 30; trMin[kNdim] = 100.; trMax[kNdim] = 3100; | |
2027 | ||
3ceb45ae | 2028 | st = "tracklet spatial&charge resolution;"; |
3ed01fbe | 2029 | // define minimum info to be saved in non debug mode |
566c3d46 | 2030 | Int_t ndim=DebugLevel()>=1?(kNdim+1):kNdimTrklt; |
3ed01fbe | 2031 | for(Int_t idim(0); idim<ndim; idim++){ st += trTitle[idim]; st+=";";} |
2032 | H = new THnSparseI(hn, st.Data(), ndim, trNbins, trMin, trMax); | |
3ceb45ae | 2033 | } else H->Reset(); |
2034 | fContainer->AddAt(H, kTracklet); | |
2035 | //++++++++++++++++++++++ | |
afca20ef | 2036 | // tracklet to TRDin |
3ceb45ae | 2037 | snprintf(hn, nhn, "h%s", fgPerformanceName[kTrackIn]); |
2038 | if(!(H = (THnSparseI*)gROOT->FindObject(hn))){ | |
35983729 | 2039 | // set specific fields |
566c3d46 | 2040 | const Int_t mdim(kNdim+3); |
35983729 | 2041 | Char_t *trinTitle[mdim]; memcpy(trinTitle, fgkTitle, kNdim*sizeof(Char_t*)); |
2042 | Int_t trinNbins[mdim]; memcpy(trinNbins, fgkNbins, kNdim*sizeof(Int_t)); | |
2043 | Double_t trinMin[mdim]; memcpy(trinMin, fgkMin, kNdim*sizeof(Double_t)); | |
2044 | Double_t trinMax[mdim]; memcpy(trinMax, fgkMax, kNdim*sizeof(Double_t)); | |
35983729 | 2045 | trinTitle[kNdim]=StrDup("detector"); trinNbins[kNdim] = 540; trinMin[kNdim] = -0.5; trinMax[kNdim] = 539.5; |
2046 | trinTitle[kNdim+1]=StrDup("dx [cm]"); trinNbins[kNdim+1]=48; trinMin[kNdim+1]=-2.4; trinMax[kNdim+1]=2.4; | |
566c3d46 | 2047 | trinTitle[kNdim+2]=StrDup("Fill Bunch"); trinNbins[kNdim+2]=3500; trinMin[kNdim+2]=-0.5; trinMax[kNdim+2]=3499.5; |
3ceb45ae | 2048 | st = "r-#phi/z/angular residuals @ TRD entry;"; |
3ed01fbe | 2049 | // define minimum info to be saved in non debug mode |
566c3d46 | 2050 | Int_t ndim=DebugLevel()>=1?mdim:kNdimTrkIn; |
35983729 | 2051 | for(Int_t idim(0); idim<ndim; idim++){st+=trinTitle[idim]; st+=";";} |
2052 | H = new THnSparseI(hn, st.Data(), ndim, trinNbins, trinMin, trinMax); | |
3ceb45ae | 2053 | } else H->Reset(); |
2054 | fContainer->AddAt(H, kTrackIn); | |
a310e49b | 2055 | // tracklet to TRDout |
3ed01fbe | 2056 | // fContainer->AddAt(BuildMonitorContainerTracklet("TrkOUT"), kTrackOut); |
1ee39b3a | 2057 | |
2058 | ||
2059 | // Resolution histos | |
2060 | if(!HasMCdata()) return fContainer; | |
2061 | ||
f429b017 | 2062 | //++++++++++++++++++++++ |
2063 | // cluster to TrackRef residuals/pulls | |
2064 | snprintf(hn, nhn, "h%s", fgPerformanceName[kMCcluster]); | |
2065 | if(!(H = (THnSparseI*)gROOT->FindObject(hn))){ | |
2066 | const Char_t *clTitle[kNdim] = {"layer", fgkTitle[kPhi], fgkTitle[kEta], fgkTitle[kYrez], "#Deltax [cm]", "Q</Q", fgkTitle[kSpeciesChgRC], "#Phi [deg]"}; | |
566c3d46 | 2067 | const Int_t clNbins[kNdim] = {AliTRDgeometry::kNlayer, fgkNbins[kPhi], fgkNbins[kEta], fgkNbins[kYrez], 20, 10, Int_t(kNcharge)*AliPID::kSPECIES+1, 15}; |
2068 | const Double_t clMin[kNdim] = {-0.5, fgkMin[kPhi], fgkMin[kEta], fgkMin[kYrez]/10., 0., 0.1, -AliPID::kSPECIES-0.5, -45}, | |
2069 | clMax[kNdim] = {AliTRDgeometry::kNlayer-0.5, fgkMax[kPhi], fgkMax[kEta], fgkMax[kYrez]/10., 4., 2.1, AliPID::kSPECIES+0.5, 45}; | |
f429b017 | 2070 | st = "MC cluster spatial resolution;"; |
2071 | // define minimum info to be saved in non debug mode | |
2072 | Int_t ndim=DebugLevel()>=1?kNdim:4; | |
2073 | for(Int_t idim(0); idim<ndim; idim++){ st += clTitle[idim]; st+=";";} | |
2074 | H = new THnSparseI(hn, st.Data(), ndim, clNbins, clMin, clMax); | |
2075 | } else H->Reset(); | |
2076 | fContainer->AddAt(H, kMCcluster); | |
2077 | //++++++++++++++++++++++ | |
2078 | // tracklet to TrackRef | |
2079 | snprintf(hn, nhn, "h%s", fgPerformanceName[kMCtracklet]); | |
2080 | if(!(H = (THnSparseI*)gROOT->FindObject(hn))){ | |
2081 | Char_t *trTitle[kNdim]; memcpy(trTitle, fgkTitle, kNdim*sizeof(Char_t*)); | |
2082 | Int_t trNbins[kNdim]; memcpy(trNbins, fgkNbins, kNdim*sizeof(Int_t)); | |
2083 | Double_t trMin[kNdim]; memcpy(trMin, fgkMin, kNdim*sizeof(Double_t)); | |
2084 | Double_t trMax[kNdim]; memcpy(trMax, fgkMax, kNdim*sizeof(Double_t)); | |
2085 | // set specific fields | |
2086 | trTitle[kBC]=StrDup("layer"); trNbins[kBC] = AliTRDgeometry::kNlayer; trMin[kBC] = -0.5; trMax[kBC] = AliTRDgeometry::kNlayer-0.5; | |
2087 | trMin[kYrez] = -0.54; trMax[kYrez] = -trMin[kYrez]; | |
2088 | trMin[kPrez] = -4.5; trMax[kPrez] = -trMin[kPrez]; | |
2089 | trMin[kZrez] = -1.5; trMax[kZrez] = -trMin[kZrez]; | |
566c3d46 | 2090 | trNbins[kSpeciesChgRC] = Int_t(kNcharge)*AliPID::kSPECIES+1;trMin[kSpeciesChgRC] = -AliPID::kSPECIES-0.5; trMax[kSpeciesChgRC] = AliPID::kSPECIES+0.5; |
31c8fa8a | 2091 | |
f429b017 | 2092 | st = "MC tracklet spatial resolution;"; |
2093 | // define minimum info to be saved in non debug mode | |
2094 | Int_t ndim=DebugLevel()>=1?kNdim:4; | |
2095 | for(Int_t idim(0); idim<ndim; idim++){ st += trTitle[idim]; st+=";";} | |
2096 | H = new THnSparseI(hn, st.Data(), ndim, trNbins, trMin, trMax); | |
2097 | } else H->Reset(); | |
2098 | fContainer->AddAt(H, kMCtracklet); | |
2099 | //++++++++++++++++++++++ | |
2100 | // TRDin to TrackRef | |
2101 | snprintf(hn, nhn, "h%s", fgPerformanceName[kMCtrackIn]); | |
2102 | if(!(H = (THnSparseI*)gROOT->FindObject(hn))){ | |
2103 | st = "MC r-#phi/z/angular residuals @ TRD entry;"; | |
2104 | // set specific fields | |
566c3d46 | 2105 | Int_t trNbins[kNdim]; memcpy(trNbins, fgkNbins, kNdim*sizeof(Int_t)); |
f429b017 | 2106 | Double_t trMin[kNdim]; memcpy(trMin, fgkMin, kNdim*sizeof(Double_t)); |
2107 | Double_t trMax[kNdim]; memcpy(trMax, fgkMax, kNdim*sizeof(Double_t)); | |
2108 | trMin[kYrez] = -0.54; trMax[kYrez] = -trMin[kYrez]; | |
2109 | trMin[kPrez] = -2.4; trMax[kPrez] = -trMin[kPrez]; | |
2110 | trMin[kZrez] = -0.9; trMax[kZrez] = -trMin[kZrez]; | |
566c3d46 | 2111 | trNbins[kSpeciesChgRC] = Int_t(kNcharge)*AliPID::kSPECIES+1;trMin[kSpeciesChgRC] = -AliPID::kSPECIES-0.5; trMax[kSpeciesChgRC] = AliPID::kSPECIES+0.5; |
f429b017 | 2112 | // define minimum info to be saved in non debug mode |
2113 | Int_t ndim=DebugLevel()>=1?kNdim:7; | |
2114 | for(Int_t idim(0); idim<ndim; idim++){ st += fgkTitle[idim]; st+=";";} | |
2115 | H = new THnSparseI(hn, st.Data(), ndim, fgkNbins, trMin, trMax); | |
2116 | } else H->Reset(); | |
3ceb45ae | 2117 | fContainer->AddAt(H, kMCtrackIn); |
f429b017 | 2118 | //++++++++++++++++++++++ |
2119 | // track to TrackRef | |
2120 | snprintf(hn, nhn, "h%s", fgPerformanceName[kMCtrack]); | |
2121 | if(!(H = (THnSparseI*)gROOT->FindObject(hn))){ | |
2122 | Char_t *trTitle[kNdim+1]; memcpy(trTitle, fgkTitle, kNdim*sizeof(Char_t*)); | |
2123 | Int_t trNbins[kNdim+1]; memcpy(trNbins, fgkNbins, kNdim*sizeof(Int_t)); | |
2124 | Double_t trMin[kNdim+1]; memcpy(trMin, fgkMin, kNdim*sizeof(Double_t)); | |
2125 | Double_t trMax[kNdim+1]; memcpy(trMax, fgkMax, kNdim*sizeof(Double_t)); | |
2126 | // set specific fields | |
2127 | trTitle[kBC]=StrDup("layer"); trNbins[kBC] = AliTRDgeometry::kNlayer; trMin[kBC] = -0.5; trMax[kBC] = AliTRDgeometry::kNlayer-0.5; | |
f429b017 | 2128 | trMin[kYrez] = -0.9; trMax[kYrez] = -trMin[kYrez]; |
2129 | trMin[kPrez] = -1.5; trMax[kPrez] = -trMin[kPrez]; | |
2130 | trMin[kZrez] = -0.9; trMax[kZrez] = -trMin[kZrez]; | |
566c3d46 | 2131 | trNbins[kSpeciesChgRC] = Int_t(kNcharge)*AliPID::kSPECIES+1;trMin[kSpeciesChgRC] = -AliPID::kSPECIES-0.5; trMax[kSpeciesChgRC] = AliPID::kSPECIES+0.5; |
2132 | trTitle[kNdim]=StrDup("#Deltap_{t}/p_{t} [%]"); trNbins[kNdim] = 25; trMin[kNdim] = -4.5; trMax[kNdim] = 20.5; | |
1ee39b3a | 2133 | |
f429b017 | 2134 | st = "MC track spatial&p_{t} resolution;"; |
2135 | // define minimum info to be saved in non debug mode | |
2136 | Int_t ndim=DebugLevel()>=1?(kNdim+1):4; | |
2137 | for(Int_t idim(0); idim<ndim; idim++){ st += trTitle[idim]; st+=";";} | |
2138 | H = new THnSparseI(hn, st.Data(), ndim, trNbins, trMin, trMax); | |
2139 | } else H->Reset(); | |
2140 | fContainer->AddAt(H, kMCtrack); | |
2141 | ||
1ee39b3a | 2142 | return fContainer; |
2143 | } | |
2144 | ||
5468a262 | 2145 | //________________________________________________________ |
2146 | Bool_t AliTRDresolution::Process(TH2* const h2, TGraphErrors **g, Int_t stat) | |
2147 | { | |
b37d601d | 2148 | // Robust function to process sigma/mean for 2D plot dy(x) |
2149 | // For each x bin a gauss fit is performed on the y projection and the range | |
2150 | // with the minimum chi2/ndf is choosen | |
5468a262 | 2151 | |
2152 | if(!h2) { | |
2153 | if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>0) printf("D-AliTRDresolution::Process() : NULL pointer input container.\n"); | |
2154 | return kFALSE; | |
2155 | } | |
2156 | if(!Int_t(h2->GetEntries())){ | |
2157 | if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>0) printf("D-AliTRDresolution::Process() : Empty h[%s - %s].\n", h2->GetName(), h2->GetTitle()); | |
2158 | return kFALSE; | |
2159 | } | |
2160 | if(!g || !g[0]|| !g[1]) { | |
2161 | if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>0) printf("D-AliTRDresolution::Process() : NULL pointer output container.\n"); | |
2162 | return kFALSE; | |
2163 | } | |
2164 | // prepare | |
2165 | TAxis *ax(h2->GetXaxis()), *ay(h2->GetYaxis()); | |
b37d601d | 2166 | Float_t ymin(ay->GetXmin()), ymax(ay->GetXmax()), dy(ay->GetBinWidth(1)), y0(0.), y1(0.); |
2167 | TF1 f("f", "gaus", ymin, ymax); | |
5468a262 | 2168 | Int_t n(0); |
2169 | if((n=g[0]->GetN())) for(;n--;) g[0]->RemovePoint(n); | |
2170 | if((n=g[1]->GetN())) for(;n--;) g[1]->RemovePoint(n); | |
2171 | TH1D *h(NULL); | |
2172 | if((h=(TH1D*)gROOT->FindObject("py"))) delete h; | |
b37d601d | 2173 | Double_t x(0.), y(0.), ex(0.), ey(0.), sy(0.), esy(0.); |
2174 | ||
5468a262 | 2175 | |
2176 | // do actual loop | |
b37d601d | 2177 | Float_t chi2OverNdf(0.); |
5468a262 | 2178 | for(Int_t ix = 1, np=0; ix<=ax->GetNbins(); ix++){ |
b37d601d | 2179 | x = ax->GetBinCenter(ix); ex= ax->GetBinWidth(ix)*0.288; // w/sqrt(12) |
2180 | ymin = ay->GetXmin(); ymax = ay->GetXmax(); | |
2181 | ||
5468a262 | 2182 | h = h2->ProjectionY("py", ix, ix); |
2183 | if((n=(Int_t)h->GetEntries())<stat){ | |
2184 | if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>1) printf("I-AliTRDresolution::Process() : Low statistics @ x[%f] stat[%d]=%d [%d].\n", x, ix, n, stat); | |
2185 | continue; | |
2186 | } | |
b37d601d | 2187 | // looking for a first order mean value |
2188 | f.SetParameter(1, 0.); f.SetParameter(2, 3.e-2); | |
2189 | h->Fit(&f, "QNW"); | |
2190 | chi2OverNdf = f.GetChisquare()/f.GetNDF(); | |
2191 | printf("x[%f] range[%f %f] chi2[%f] ndf[%d] chi2/ndf[%f]\n", x, ymin, ymax, f.GetChisquare(),f.GetNDF(),chi2OverNdf); | |
2192 | y = f.GetParameter(1); y0 = y-4*dy; y1 = y+4*dy; | |
2193 | ey = f.GetParError(1); | |
2194 | sy = f.GetParameter(2); esy = f.GetParError(2); | |
2195 | // // looking for the best chi2/ndf value | |
2196 | // while(ymin<y0 && ymax>y1){ | |
2197 | // f.SetParameter(1, y); | |
2198 | // f.SetParameter(2, sy); | |
2199 | // h->Fit(&f, "QNW", "", y0, y1); | |
2200 | // printf(" range[%f %f] chi2[%f] ndf[%d] chi2/ndf[%f]\n", y0, y1, f.GetChisquare(),f.GetNDF(),f.GetChisquare()/f.GetNDF()); | |
2201 | // if(f.GetChisquare()/f.GetNDF() < Chi2OverNdf){ | |
2202 | // chi2OverNdf = f.GetChisquare()/f.GetNDF(); | |
2203 | // y = f.GetParameter(1); ey = f.GetParError(1); | |
2204 | // sy = f.GetParameter(2); esy = f.GetParError(2); | |
2205 | // printf(" set y[%f] sy[%f] chi2/ndf[%f]\n", y, sy, chi2OverNdf); | |
2206 | // } | |
2207 | // y0-=dy; y1+=dy; | |
2208 | // } | |
2209 | g[0]->SetPoint(np, x, y); | |
2210 | g[0]->SetPointError(np, ex, ey); | |
2211 | g[1]->SetPoint(np, x, sy); | |
2212 | g[1]->SetPointError(np, ex, esy); | |
5468a262 | 2213 | np++; |
2214 | } | |
2215 | return kTRUE; | |
2216 | } | |
2217 | ||
2589cf64 | 2218 | |
1ee39b3a | 2219 | //________________________________________________________ |
2220 | Bool_t AliTRDresolution::Process(TH2 * const h2, TF1 *f, Float_t k, TGraphErrors **g) | |
2221 | { | |
2222 | // | |
2223 | // Do the processing | |
2224 | // | |
2225 | ||
7fe4e88b | 2226 | Char_t pn[10]; snprintf(pn, 10, "p%03d", fIdxPlot); |
1ee39b3a | 2227 | Int_t n = 0; |
2228 | if((n=g[0]->GetN())) for(;n--;) g[0]->RemovePoint(n); | |
2229 | if((n=g[1]->GetN())) for(;n--;) g[1]->RemovePoint(n); | |
2589cf64 | 2230 | if(Int_t(h2->GetEntries())){ |
2231 | AliDebug(4, Form("%s: g[%s %s]", pn, g[0]->GetName(), g[0]->GetTitle())); | |
2232 | } else { | |
2233 | AliDebug(2, Form("%s: g[%s %s]: Missing entries.", pn, g[0]->GetName(), g[0]->GetTitle())); | |
2234 | fIdxPlot++; | |
2235 | return kTRUE; | |
2236 | } | |
1ee39b3a | 2237 | |
dfd7d48b | 2238 | const Int_t kINTEGRAL=1; |
2239 | for(Int_t ibin = 0; ibin < Int_t(h2->GetNbinsX()/kINTEGRAL); ibin++){ | |
2240 | Int_t abin(ibin*kINTEGRAL+1),bbin(abin+kINTEGRAL-1),mbin(abin+Int_t(kINTEGRAL/2)); | |
2241 | Double_t x = h2->GetXaxis()->GetBinCenter(mbin); | |
2242 | TH1D *h = h2->ProjectionY(pn, abin, bbin); | |
068e2c00 | 2243 | if((n=(Int_t)h->GetEntries())<150){ |
2589cf64 | 2244 | AliDebug(4, Form(" x[%f] range[%d %d] stat[%d] low statistics !", x, abin, bbin, n)); |
2245 | continue; | |
2246 | } | |
1ee39b3a | 2247 | h->Fit(f, "QN"); |
1ee39b3a | 2248 | Int_t ip = g[0]->GetN(); |
2589cf64 | 2249 | AliDebug(4, Form(" x_%d[%f] range[%d %d] stat[%d] M[%f] Sgm[%f]", ip, x, abin, bbin, n, f->GetParameter(1), f->GetParameter(2))); |
1ee39b3a | 2250 | g[0]->SetPoint(ip, x, k*f->GetParameter(1)); |
2251 | g[0]->SetPointError(ip, 0., k*f->GetParError(1)); | |
2252 | g[1]->SetPoint(ip, x, k*f->GetParameter(2)); | |
2253 | g[1]->SetPointError(ip, 0., k*f->GetParError(2)); | |
1ee39b3a | 2254 | /* |
2255 | g[0]->SetPoint(ip, x, k*h->GetMean()); | |
2256 | g[0]->SetPointError(ip, 0., k*h->GetMeanError()); | |
2257 | g[1]->SetPoint(ip, x, k*h->GetRMS()); | |
2258 | g[1]->SetPointError(ip, 0., k*h->GetRMSError());*/ | |
2259 | } | |
2260 | fIdxPlot++; | |
2261 | return kTRUE; | |
2262 | } | |
2263 | ||
1ee39b3a | 2264 | |
08bdd9d1 | 2265 | //____________________________________________________________________ |
2266 | Bool_t AliTRDresolution::FitTrack(const Int_t np, AliTrackPoint *points, Float_t param[10]) | |
2267 | { | |
2268 | // | |
2269 | // Fit track with a staight line using the "np" clusters stored in the array "points". | |
2270 | // The following particularities are stored in the clusters from points: | |
2271 | // 1. pad tilt as cluster charge | |
2272 | // 2. pad row cross or vertex constrain as fake cluster with cluster type 1 | |
2273 | // The parameters of the straight line fit are stored in the array "param" in the following order : | |
2274 | // param[0] - x0 reference radial position | |
2275 | // param[1] - y0 reference r-phi position @ x0 | |
2276 | // param[2] - z0 reference z position @ x0 | |
2277 | // param[3] - slope dy/dx | |
2278 | // param[4] - slope dz/dx | |
2279 | // | |
2280 | // Attention : | |
2281 | // Function should be used to refit tracks for B=0T | |
2282 | // | |
2283 | ||
2284 | if(np<40){ | |
b37d601d | 2285 | if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>1) printf("D-AliTRDresolution::FitTrack: Not enough clusters to fit a track [%d].\n", np); |
08bdd9d1 | 2286 | return kFALSE; |
2287 | } | |
2288 | TLinearFitter yfitter(2, "pol1"), zfitter(2, "pol1"); | |
2289 | ||
2290 | Double_t x0(0.); | |
2291 | for(Int_t ip(0); ip<np; ip++) x0+=points[ip].GetX(); | |
2292 | x0/=Float_t(np); | |
2293 | ||
00a56108 | 2294 | Double_t x, y, z, dx, tilt(0.); |
08bdd9d1 | 2295 | for(Int_t ip(0); ip<np; ip++){ |
2296 | x = points[ip].GetX(); z = points[ip].GetZ(); | |
2297 | dx = x - x0; | |
2298 | zfitter.AddPoint(&dx, z, points[ip].GetClusterType()?1.e-3:1.); | |
2299 | } | |
2300 | if(zfitter.Eval() != 0) return kFALSE; | |
2301 | ||
2302 | Double_t z0 = zfitter.GetParameter(0); | |
2303 | Double_t dzdx = zfitter.GetParameter(1); | |
2304 | for(Int_t ip(0); ip<np; ip++){ | |
2305 | if(points[ip].GetClusterType()) continue; | |
2306 | x = points[ip].GetX(); | |
2307 | dx = x - x0; | |
2308 | y = points[ip].GetY(); | |
2309 | z = points[ip].GetZ(); | |
2310 | tilt = points[ip].GetCharge(); | |
2311 | y -= tilt*(-dzdx*dx + z - z0); | |
2312 | Float_t xyz[3] = {x, y, z}; points[ip].SetXYZ(xyz); | |
2313 | yfitter.AddPoint(&dx, y, 1.); | |
2314 | } | |
2315 | if(yfitter.Eval() != 0) return kFALSE; | |
2316 | Double_t y0 = yfitter.GetParameter(0); | |
2317 | Double_t dydx = yfitter.GetParameter(1); | |
2318 | ||
2319 | param[0] = x0; param[1] = y0; param[2] = z0; param[3] = dydx; param[4] = dzdx; | |
b37d601d | 2320 | if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>3) printf("D-AliTRDresolution::FitTrack: x0[%f] y0[%f] z0[%f] dydx[%f] dzdx[%f].\n", x0, y0, z0, dydx, dzdx); |
08bdd9d1 | 2321 | return kTRUE; |
2322 | } | |
2323 | ||
5f7f1c07 | 2324 | //____________________________________________________________________ |
b37d601d | 2325 | Bool_t AliTRDresolution::FitTracklet(const Int_t ly, const Int_t np, const AliTrackPoint *points, const Float_t param[10], Float_t par[3]) |
5f7f1c07 | 2326 | { |
2327 | // | |
2328 | // Fit tracklet with a staight line using the coresponding subset of clusters out of the total "np" clusters stored in the array "points". | |
2329 | // See function FitTrack for the data stored in the "clusters" array | |
2330 | ||
2331 | // The parameters of the straight line fit are stored in the array "param" in the following order : | |
2332 | // par[0] - x0 reference radial position | |
2333 | // par[1] - y0 reference r-phi position @ x0 | |
2334 | // par[2] - slope dy/dx | |
2335 | // | |
2336 | // Attention : | |
2337 | // Function should be used to refit tracks for B=0T | |
2338 | // | |
2339 | ||
2340 | TLinearFitter yfitter(2, "pol1"); | |
2341 | ||
2342 | // grep data for tracklet | |
2343 | Double_t x0(0.), x[60], y[60], dy[60]; | |
2344 | Int_t nly(0); | |
2345 | for(Int_t ip(0); ip<np; ip++){ | |
2346 | if(points[ip].GetClusterType()) continue; | |
2347 | if(points[ip].GetVolumeID() != ly) continue; | |
2348 | Float_t xt(points[ip].GetX()) | |
2349 | ,yt(param[1] + param[3] * (xt - param[0])); | |
2350 | x[nly] = xt; | |
2351 | y[nly] = points[ip].GetY(); | |
2352 | dy[nly]= y[nly]-yt; | |
2353 | x0 += xt; | |
2354 | nly++; | |
2355 | } | |
2356 | if(nly<10){ | |
b37d601d | 2357 | if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>1) printf("D-AliTRDresolution::FitTracklet: Not enough clusters to fit a tracklet [%d].\n", nly); |
5f7f1c07 | 2358 | return kFALSE; |
2359 | } | |
2360 | // set radial reference for fit | |
2361 | x0 /= Float_t(nly); | |
2362 | ||
2363 | // find tracklet core | |
2364 | Double_t mean(0.), sig(1.e3); | |
2365 | AliMathBase::EvaluateUni(nly, dy, mean, sig, 0); | |
2366 | ||
2367 | // simple cluster error parameterization | |
2368 | Float_t kSigCut = TMath::Sqrt(5.e-4 + param[3]*param[3]*0.018); | |
2369 | ||
2370 | // fit tracklet core | |
2371 | for(Int_t jly(0); jly<nly; jly++){ | |
2372 | if(TMath::Abs(dy[jly]-mean)>kSigCut) continue; | |
2373 | Double_t dx(x[jly]-x0); | |
2374 | yfitter.AddPoint(&dx, y[jly], 1.); | |
2375 | } | |
2376 | if(yfitter.Eval() != 0) return kFALSE; | |
2377 | par[0] = x0; | |
2378 | par[1] = yfitter.GetParameter(0); | |
2379 | par[2] = yfitter.GetParameter(1); | |
2380 | return kTRUE; | |
2381 | } | |
2382 | ||
00a56108 | 2383 | //____________________________________________________________________ |
b37d601d | 2384 | Bool_t AliTRDresolution::UseTrack(const Int_t np, const AliTrackPoint *points, Float_t param[10]) |
00a56108 | 2385 | { |
2386 | // | |
2387 | // Global selection mechanism of tracksbased on cluster to fit residuals | |
2388 | // The parameters are the same as used ni function FitTrack(). | |
2389 | ||
2390 | const Float_t kS(0.6), kM(0.2); | |
2391 | TH1S h("h1", "", 100, -5.*kS, 5.*kS); | |
2392 | Float_t dy, dz, s, m; | |
2393 | for(Int_t ip(0); ip<np; ip++){ | |
2394 | if(points[ip].GetClusterType()) continue; | |
2395 | Float_t x0(points[ip].GetX()) | |
2396 | ,y0(param[1] + param[3] * (x0 - param[0])) | |
2397 | ,z0(param[2] + param[4] * (x0 - param[0])); | |
2398 | dy=points[ip].GetY() - y0; h.Fill(dy); | |
2399 | dz=points[ip].GetZ() - z0; | |
2400 | } | |
2401 | TF1 fg("fg", "gaus", -5.*kS, 5.*kS); | |
2402 | fg.SetParameter(1, 0.); | |
2403 | fg.SetParameter(2, 2.e-2); | |
2404 | h.Fit(&fg, "QN"); | |
2405 | m=fg.GetParameter(1); s=fg.GetParameter(2); | |
2406 | if(s>kS || TMath::Abs(m)>kM) return kFALSE; | |
2407 | return kTRUE; | |
2408 | } | |
2409 | ||
1ee39b3a | 2410 | //________________________________________________________ |
2411 | void AliTRDresolution::GetLandauMpvFwhm(TF1 * const f, Float_t &mpv, Float_t &xm, Float_t &xM) | |
2412 | { | |
2413 | // | |
2414 | // Get the most probable value and the full width half mean | |
2415 | // of a Landau distribution | |
2416 | // | |
2417 | ||
2418 | const Float_t dx = 1.; | |
2419 | mpv = f->GetParameter(1); | |
2420 | Float_t fx, max = f->Eval(mpv); | |
2421 | ||
2422 | xm = mpv - dx; | |
2423 | while((fx = f->Eval(xm))>.5*max){ | |
2424 | if(fx>max){ | |
2425 | max = fx; | |
2426 | mpv = xm; | |
2427 | } | |
2428 | xm -= dx; | |
2429 | } | |
2430 | ||
2431 | xM += 2*(mpv - xm); | |
2432 | while((fx = f->Eval(xM))>.5*max) xM += dx; | |
2433 | } | |
2434 | ||
2435 | ||
3ceb45ae | 2436 | // #include "TFile.h" |
2437 | // //________________________________________________________ | |
2438 | // Bool_t AliTRDresolution::LoadCorrection(const Char_t *file) | |
2439 | // { | |
2440 | // if(!file){ | |
2441 | // AliWarning("Use cluster position as in reconstruction."); | |
2442 | // SetLoadCorrection(); | |
2443 | // return kTRUE; | |
2444 | // } | |
2445 | // TDirectory *cwd(gDirectory); | |
2446 | // TString fileList; | |
2447 | // FILE *filePtr = fopen(file, "rt"); | |
2448 | // if(!filePtr){ | |
2449 | // AliWarning(Form("Couldn't open correction list \"%s\". Use cluster position as in reconstruction.", file)); | |
2450 | // SetLoadCorrection(); | |
2451 | // return kFALSE; | |
2452 | // } | |
2453 | // TH2 *h2 = new TH2F("h2", ";time [time bins];detector;dx [#mum]", 30, -0.5, 29.5, AliTRDgeometry::kNdet, -0.5, AliTRDgeometry::kNdet-0.5); | |
2454 | // while(fileList.Gets(filePtr)){ | |
2455 | // if(!TFile::Open(fileList.Data())) { | |
2456 | // AliWarning(Form("Couldn't open \"%s\"", fileList.Data())); | |
2457 | // continue; | |
2458 | // } else AliInfo(Form("\"%s\"", fileList.Data())); | |
2459 | // | |
2460 | // TTree *tSys = (TTree*)gFile->Get("tSys"); | |
2461 | // h2->SetDirectory(gDirectory); h2->Reset("ICE"); | |
2462 | // tSys->Draw("det:t>>h2", "dx", "goff"); | |
2463 | // for(Int_t idet(0); idet<AliTRDgeometry::kNdet; idet++){ | |
2464 | // for(Int_t it(0); it<30; it++) fXcorr[idet][it]+=(1.e-4*h2->GetBinContent(it+1, idet+1)); | |
2465 | // } | |
2466 | // h2->SetDirectory(cwd); | |
2467 | // gFile->Close(); | |
2468 | // } | |
2469 | // cwd->cd(); | |
2470 | // | |
2471 | // if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>=2){ | |
2472 | // for(Int_t il(0); il<184; il++) printf("-"); printf("\n"); | |
2473 | // printf("DET|");for(Int_t it(0); it<30; it++) printf(" tb%02d|", it); printf("\n"); | |
2474 | // for(Int_t il(0); il<184; il++) printf("-"); printf("\n"); | |
2475 | // FILE *fout = fopen("TRD.ClusterCorrection.txt", "at"); | |
2476 | // fprintf(fout, " static const Double_t dx[AliTRDgeometry::kNdet][30]={\n"); | |
2477 | // for(Int_t idet(0); idet<AliTRDgeometry::kNdet; idet++){ | |
2478 | // printf("%03d|", idet); | |
2479 | // fprintf(fout, " {"); | |
2480 | // for(Int_t it(0); it<30; it++){ | |
2481 | // printf("%+5.0f|", 1.e4*fXcorr[idet][it]); | |
2482 | // fprintf(fout, "%+6.4f%c", fXcorr[idet][it], it==29?' ':','); | |
2483 | // } | |
2484 | // printf("\n"); | |
2485 | // fprintf(fout, "}%c\n", idet==AliTRDgeometry::kNdet-1?' ':','); | |
2486 | // } | |
2487 | // fprintf(fout, " };\n"); | |
2488 | // } | |
2489 | // SetLoadCorrection(); | |
2490 | // return kTRUE; | |
2491 | // } | |
3ed01fbe | 2492 | |
2493 | //________________________________________________________ | |
2494 | AliTRDresolution::AliTRDresolutionProjection::AliTRDresolutionProjection() | |
2495 | :fH(NULL) | |
943761d3 | 2496 | ,fNrebin(0) |
3ed01fbe | 2497 | ,fRebinX(NULL) |
2498 | ,fRebinY(NULL) | |
2499 | { | |
2500 | // constructor | |
2501 | memset(fAx, 0, 3*sizeof(Int_t)); | |
0b30040d | 2502 | memset(fRange, 0, 4*sizeof(Float_t)); |
3ed01fbe | 2503 | } |
2504 | ||
2505 | //________________________________________________________ | |
2506 | AliTRDresolution::AliTRDresolutionProjection::~AliTRDresolutionProjection() | |
2507 | { | |
2508 | // destructor | |
2509 | if(fH) delete fH; | |
2510 | } | |
2511 | ||
2512 | //________________________________________________________ | |
2513 | void AliTRDresolution::AliTRDresolutionProjection::Build(const Char_t *n, const Char_t *t, Int_t ix, Int_t iy, Int_t iz, TAxis *aa[]) | |
2514 | { | |
2515 | // check and build (if neccessary) projection determined by axis "ix", "iy" and "iz" | |
09c85be5 | 2516 | if(!aa[ix] || !aa[iy] || !aa[iz]) return; |
3ed01fbe | 2517 | TAxis *ax(aa[ix]), *ay(aa[iy]), *az(aa[iz]); |
2518 | fH = new TH3I(n, Form("%s;%s;%s;%s", t, ax->GetTitle(), ay->GetTitle(), az->GetTitle()), | |
2519 | ax->GetNbins(), ax->GetXmin(), ax->GetXmax(), | |
2520 | ay->GetNbins(), ay->GetXmin(), ay->GetXmax(), | |
2521 | az->GetNbins(), az->GetXmin(), az->GetXmax()); | |
2522 | fAx[0] = ix; fAx[1] = iy; fAx[2] = iz; | |
2523 | fRange[0] = az->GetXmin()/3.; fRange[1] = az->GetXmax()/3.; | |
2524 | } | |
2525 | ||
35983729 | 2526 | //________________________________________________________ |
2527 | AliTRDresolution::AliTRDresolutionProjection& AliTRDresolution::AliTRDresolutionProjection::operator+=(const AliTRDresolutionProjection& other) | |
2528 | { | |
2529 | // increment projections | |
2530 | if(!fH || !other.fH) return *this; | |
2531 | fH->Add(other.fH); | |
2532 | return *this; | |
2533 | } | |
2534 | ||
3ed01fbe | 2535 | //________________________________________________________ |
2536 | void AliTRDresolution::AliTRDresolutionProjection::Increment(Int_t bin[], Double_t v) | |
2537 | { | |
2538 | // increment bin with value "v" pointed by general coord in "bin" | |
2539 | if(!fH) return; | |
35983729 | 2540 | fH->AddBinContent(fH->GetBin(bin[fAx[0]],bin[fAx[1]],bin[fAx[2]]), Int_t(v)); |
3ed01fbe | 2541 | } |
2542 | ||
2543 | //________________________________________________________ | |
2544 | TH2* AliTRDresolution::AliTRDresolutionProjection::Projection2D(const Int_t nstat, const Int_t ncol, const Int_t mid) | |
2545 | { | |
2546 | // build the 2D projection and adjust binning | |
2547 | ||
0b30040d | 2548 | const Char_t *title[] = {"Mean", "#mu", "MPV"}; |
3ed01fbe | 2549 | if(!fH) return NULL; |
2550 | TAxis *ax(fH->GetXaxis()), *ay(fH->GetYaxis()), *az(fH->GetZaxis()); | |
cc0c5c06 | 2551 | TH2 *h2s(NULL); |
2552 | if(!(h2s = (TH2*)fH->Project3D("yx"))) return NULL; | |
3ed01fbe | 2553 | Int_t irebin(0), dxBin(1), dyBin(1); |
2554 | while(irebin<fNrebin && (AliTRDresolution::GetMeanStat(h2s, .5, ">")<nstat)){ | |
2555 | h2s->Rebin2D(fRebinX[irebin], fRebinY[irebin]); | |
2556 | dxBin*=fRebinX[irebin];dyBin*=fRebinY[irebin]; | |
3ed01fbe | 2557 | irebin++; |
2558 | } | |
2559 | Int_t nx(h2s->GetNbinsX()), ny(h2s->GetNbinsY()); | |
2560 | if(h2s) delete h2s; | |
2561 | ||
2562 | // start projection | |
2563 | TH1 *h(NULL); | |
2564 | Float_t dz=(fRange[1]-fRange[1])/ncol; | |
0b30040d | 2565 | TString titlez(az->GetTitle()); TObjArray *tokenTitle(titlez.Tokenize(" ")); |
3ed01fbe | 2566 | TH2 *h2 = new TH2F(Form("%s_2D", fH->GetName()), |
0b30040d | 2567 | Form("%s;%s;%s;%s(%s) %s", fH->GetTitle(), ax->GetTitle(), ay->GetTitle(), title[mid], (*tokenTitle)[0]->GetName(), tokenTitle->GetEntriesFast()>1?(*tokenTitle)[1]->GetName():""), |
3ed01fbe | 2568 | nx, ax->GetXmin(), ax->GetXmax(), ny, ay->GetXmin(), ay->GetXmax()); |
2569 | h2->SetContour(ncol); | |
2570 | h2->GetZaxis()->CenterTitle(); | |
2571 | h2->GetZaxis()->SetRangeUser(fRange[0], fRange[1]); | |
2572 | //printf("%s[%s] nx[%d] ny[%d]\n", h2->GetName(), h2->GetTitle(), nx, ny); | |
2573 | for(Int_t iy(0); iy<ny; iy++){ | |
2574 | for(Int_t ix(0); ix<nx; ix++){ | |
2575 | h = fH->ProjectionZ(Form("%s_z", h2->GetName()), ix*dxBin+1, (ix+1)*dxBin+1, iy*dyBin+1, (iy+1)*dyBin+1); | |
943761d3 | 2576 | Int_t ne((Int_t)h->Integral()); |
3ed01fbe | 2577 | if(ne<nstat/2){ |
2578 | h2->SetBinContent(ix+1, iy+1, -999); | |
2579 | h2->SetBinError(ix+1, iy+1, 1.); | |
2580 | }else{ | |
2581 | Float_t v(h->GetMean()), ve(h->GetRMS()); | |
2582 | if(mid==1){ | |
2583 | TF1 fg("fg", "gaus", az->GetXmin(), az->GetXmax()); | |
2584 | fg.SetParameter(0, Float_t(ne)); fg.SetParameter(1, v); fg.SetParameter(2, ve); | |
2585 | h->Fit(&fg, "WQ"); | |
2586 | v = fg.GetParameter(1); ve = fg.GetParameter(2); | |
2587 | } else if (mid==2) { | |
2588 | TF1 fl("fl", "landau", az->GetXmin(), az->GetXmax()); | |
2589 | fl.SetParameter(0, Float_t(ne)); fl.SetParameter(1, v); fl.SetParameter(2, ve); | |
2590 | h->Fit(&fl, "WQ"); | |
2591 | v = fl.GetMaximumX(); ve = fl.GetParameter(2); | |
2592 | /* TF1 fgle("gle", "[0]*TMath::Landau(x, [1], [2], 1)*TMath::Exp(-[3]*x/[1])", az->GetXmin(), az->GetXmax()); | |
2593 | fgle.SetParameter(0, fl.GetParameter(0)); | |
2594 | fgle.SetParameter(1, fl.GetParameter(1)); | |
2595 | fgle.SetParameter(2, fl.GetParameter(2)); | |
2596 | fgle.SetParameter(3, 1.);fgle.SetParLimits(3, 0., 5.); | |
2597 | h->Fit(&fgle, "WQ"); | |
2598 | v = fgle.GetMaximumX(); ve = fgle.GetParameter(2);*/ | |
2599 | } | |
2600 | if(v<fRange[0]) h2->SetBinContent(ix+1, iy+1, fRange[0]+0.1*dz); | |
2601 | else h2->SetBinContent(ix+1, iy+1, v); | |
2602 | h2->SetBinError(ix+1, iy+1, ve); | |
2603 | } | |
2604 | } | |
2605 | } | |
2606 | if(h) delete h; | |
2607 | return h2; | |
2608 | } | |
2609 | ||
2610 | void AliTRDresolution::AliTRDresolutionProjection::SetRebinStrategy(Int_t n, Int_t rebx[], Int_t reby[]) | |
2611 | { | |
2612 | // define rebinning strategy for this projection | |
2613 | fNrebin = n; | |
2614 | fRebinX = new Int_t[n]; memcpy(fRebinX, rebx, n*sizeof(Int_t)); | |
2615 | fRebinY = new Int_t[n]; memcpy(fRebinY, reby, n*sizeof(Int_t)); | |
2616 | } | |
2617 | ||
2618 |