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