]>
Commit | Line | Data |
---|---|---|
65a39007 | 1 | ////////////////////////////////////////////////////////////////////////// |
2 | // // | |
3 | // AliFast TrackMaker class. // | |
4 | // // | |
5 | // // | |
6 | ////////////////////////////////////////////////////////////////////////// | |
7 | // ---------------------------------------------------------------------// | |
8 | // // | |
9 | // origin: "res.f" fortran by Karel Safarik which was used to // | |
10 | // calculate the track resolution for TP. // | |
11 | // Different detectors and material can be selected. // | |
12 | // The basic routines compute information and error matrices // | |
13 | // used for the calculation of momentum resolution. // | |
14 | // see references: ASK KAREL?? // | |
15 | // // | |
16 | // C++ in AliFast framework: Elzbieta Richter-Was and Yiota Foka // | |
17 | // following general structure od Makers in // | |
18 | // ATLFast by R. Brun. // | |
19 | // // | |
20 | // purpose: provide a Maker which by using general basic routines of // | |
21 | // "res.f" computes the necessary elements of covariance matrix// | |
22 | // for the calculation of Track Resolution. // | |
23 | // Those elements are the product of the TrackResolMaker and // | |
24 | // are hold in TrackResol class. They are expected to be used // | |
25 | // together with additional information for the calculation of // | |
26 | // the smeared momenta. // | |
27 | // Additional information necessary for this calculation // | |
28 | // will be provided via classes or functions specific to the // | |
29 | // specific study and/or detectors. // | |
30 | // One can select the detector and/or material for a specific // | |
31 | // study. // | |
32 | // // | |
33 | // starting point: res.f will be initialy partially contained in // | |
34 | // AliFTrackResolMaker and in AliFDet // | |
35 | // It will be reorganised further in the framework of // | |
36 | // AliFast according to the needs. // | |
37 | // Names of variables are kept as in fortran code. // | |
38 | // // | |
39 | ////////////////////////////////////////////////////////////////////////// | |
40 | ||
41 | ||
42 | #ifdef WIN32 | |
43 | // there is a bug in the Microsoft VisualC++ compiler | |
44 | // this class must be compiled with optimization off on Windows | |
45 | # pragma optimize( "", off ) | |
46 | #endif | |
47 | ||
65a39007 | 48 | #include <TFile.h> |
65a39007 | 49 | #include <TH1.h> |
4c155701 | 50 | #include <TMath.h> |
51 | #include <TParticle.h> | |
52 | #include <TRandom.h> | |
65a39007 | 53 | |
65a39007 | 54 | #include "AliFDet.h" |
4c155701 | 55 | #include "AliFTrack.h" |
56 | #include "AliFTrackMaker.h" | |
57 | #include "AliFast.h" | |
5d12ce38 | 58 | #include "AliMC.h" |
65a39007 | 59 | |
4c155701 | 60 | static const Double_t kPi = TMath::Pi(); |
61 | static const Double_t k2Pi = 2*kPi; | |
62 | static const Double_t kPiHalf = kPi/2.; | |
65a39007 | 63 | extern AliFast * gAliFast; |
64 | ClassImp(AliFTrackMaker) | |
65 | ||
66 | //_____________________________________________________________________________ | |
67 | AliFTrackMaker::AliFTrackMaker() | |
68 | { | |
4c155701 | 69 | // |
70 | // Default constructor | |
71 | // | |
72 | fNTracks = 0; | |
73 | fResID1Test = 0; | |
74 | fResID2Test = 0; | |
75 | fResID3Test = 0; | |
76 | fResID4Test = 0; | |
77 | fResID5Test = 0; | |
65a39007 | 78 | } |
79 | ||
80 | //_____________________________________________________________________________ | |
81 | AliFTrackMaker::AliFTrackMaker(const char *name, const char *title) | |
82 | :AliFMaker(name,title) | |
83 | { | |
4c155701 | 84 | // |
85 | // Standard Setters for tracks | |
86 | // | |
65a39007 | 87 | fFruits = new TClonesArray("AliFTrack",100, kFALSE); |
88 | fBranchName = "Tracks"; | |
89 | fNTracks = 0; | |
90 | // Please, how to do this optionally ??!!! | |
91 | Save(); | |
92 | } | |
93 | ||
4c155701 | 94 | //_______________________________________________________________________ |
95 | AliFTrackMaker::AliFTrackMaker(const AliFTrackMaker& aftmk): | |
96 | AliFMaker(aftmk) | |
97 | { | |
98 | // | |
99 | // Copy constructor for AliRun | |
100 | // | |
101 | aftmk.Copy(*this); | |
102 | } | |
103 | ||
65a39007 | 104 | //_____________________________________________________________________________ |
105 | AliFTrackMaker::~AliFTrackMaker() | |
106 | { | |
4c155701 | 107 | // |
108 | // Dummy constructor | |
109 | // | |
65a39007 | 110 | } |
111 | ||
112 | //_____________________________________________________________________________ | |
113 | AliFTrack *AliFTrackMaker::AddTrack(Int_t code, Double_t charge, | |
114 | Double_t pT, Double_t eta,Double_t phi, | |
115 | Double_t v11, Double_t v22, Double_t v33, | |
116 | Double_t v12, Double_t v13, Double_t v23, Int_t iFlag) | |
117 | { | |
118 | // Add a new track to the list of tracks | |
119 | ||
120 | //Note the use of the "new with placement" to create a new track object. | |
121 | //This complex "new" works in the following way: | |
122 | // tracks[i] is the value of the pointer for track number i in the TClonesArray | |
123 | // if it is zero, then a new track must be generated. This typically | |
124 | // will happen only at the first events | |
125 | // If it is not zero, then the already existing object is overwritten | |
126 | // by the new track parameters. | |
127 | // This technique should save a huge amount of time otherwise spent | |
128 | // in the operators new and delete. | |
129 | ||
130 | TClonesArray &tracks = *(TClonesArray*)fFruits; | |
131 | return new(tracks[fNTracks++]) AliFTrack(code,charge,pT,eta,phi, | |
132 | v11, v22, v33, v12, v13, v23, iFlag); | |
133 | } | |
134 | ||
135 | //_____________________________________________________________________________ | |
136 | void AliFTrackMaker::Clear(Option_t *option) | |
137 | { | |
138 | //Reset Track Maker | |
139 | ||
140 | fNTracks = 0; | |
141 | AliFMaker::Clear(option); | |
142 | } | |
143 | ||
144 | //_____________________________________________________________________________ | |
145 | void AliFTrackMaker::Draw(Option_t *) | |
146 | { | |
147 | // Dummy Draw | |
148 | ||
149 | } | |
150 | ||
151 | //_____________________________________________________________________________ | |
152 | void AliFTrackMaker::Init() | |
153 | { | |
154 | //Create control histograms | |
155 | if(gAliFast->TestTrack() == 0){ | |
156 | ||
157 | fResID11 = new TH1D("ResID11","Elec: delta(1/pTotal)*pTotal",1000,-0.5,0.5); | |
158 | fResID12 = new TH1D("ResID12","Elec: delta(lambda)/lambda",1000,-0.01,0.01); | |
159 | fResID13 = new TH1D("ResID13","Elec: delta(phi)/phi",1000,-0.01,0.01); | |
160 | ||
161 | fResID21 = new TH1D("ResID21","Pion: delta(1/pTotal)*pTotal",1000,-1.0,1.0); | |
162 | fResID22 = new TH1D("ResID22","Pion: delta(lambda)/lambda",1000,-1.0,1.0); | |
163 | fResID23 = new TH1D("ResID23","Pion: delta(phi)/phi",1000,-1.0,1.0); | |
164 | ||
165 | fResID31 = new TH1D("ResID31","Kaon: delta(1/pTotal)*pTotal",1000,-1.0,1.0); | |
166 | fResID32 = new TH1D("ResID32","Kaon: delta(lambda)/lambda",1000,-1.0,1.0); | |
167 | fResID33 = new TH1D("ResID33","Kaon: delta(phi)/phi",1000,-1.0,1.0); | |
168 | ||
169 | fResID41 = new TH1D("ResID41","Proton: delta(1/pTotal)*pTotal",1000,-1.0,1.0); | |
170 | fResID42 = new TH1D("ResID42","Proton: delta(lambda)/lambda",1000,-1.0,1.0); | |
171 | fResID43 = new TH1D("ResID43","Proton: delta(phi)/phi",1000,-1.0,1.0); | |
172 | ||
173 | } | |
174 | //Create test histograms for TestJob only | |
175 | if(gAliFast->TestTrack() == 1){ | |
176 | fResID1Test = new TH1D("ResID1Test","histogram21 from res.f",1000,0.075,10.075); | |
177 | fResID2Test = new TH1D("ResID2Test","histogram21 from res.f",1000,0.075,10.075); | |
178 | fResID3Test = new TH1D("ResID3Test","histogram21 from res.f",1000,0.075,10.075); | |
179 | fResID4Test = new TH1D("ResID4Test","histogram21 from res.f",1000,0.075,10.075); | |
180 | fResID5Test = new TH1D("ResID5Test","histogram21 from res.f",1000,0.075,10.075); | |
181 | } | |
182 | ||
183 | //Set particle masses | |
184 | SetPionMass(); | |
185 | SetKaonMass(); | |
186 | SetElectronMass(); | |
187 | SetProtonMass(); | |
188 | ||
189 | //Switch on/off tracks reconstruction | |
190 | SetRecTrack(); | |
191 | ||
192 | } | |
193 | ||
65a39007 | 194 | void AliFTrackMaker::Make() |
195 | { | |
4c155701 | 196 | // |
197 | // Calculate track and its resolution | |
198 | // | |
65a39007 | 199 | Double_t v11, v22, v33, v12, v13, v23; |
200 | Int_t iFlag; | |
201 | ||
202 | fNTracks = 0; | |
203 | ||
204 | // Check if it is a TestJob | |
205 | if(gAliFast->TestTrack() == 1){ | |
206 | // Run test job | |
207 | MakeTest(10); | |
208 | }else{ | |
209 | // Run production job | |
210 | // Get pointers to Particles arrays and TClonesArray | |
211 | ||
65a39007 | 212 | Int_t idPart, idTrack; |
213 | Double_t charge, pT, eta, phi; | |
214 | TParticle *part; | |
5d12ce38 | 215 | Int_t nparticles = gAlice->GetMCApp()->GetNtrack(); |
65a39007 | 216 | printf("%10s%10d\n","nparticles",nparticles); |
217 | for(Int_t ind=0;ind<nparticles;ind++) { | |
5d12ce38 | 218 | part = gAlice->GetMCApp()->Particle(ind); |
65a39007 | 219 | idPart = part->GetPdgCode(); |
220 | charge = part->GetPDG()->Charge(); | |
221 | pT = part->Pt(); | |
222 | eta = part->Eta(); | |
223 | phi = part->Phi(); | |
224 | printf("%10s%10d%20.5e%20.5e%20.5e%20.5e\n","Particle",idPart,charge,pT,eta,phi); | |
225 | // Check convention for tracks reconstruction | |
226 | idTrack = 0; | |
227 | if(TMath::Abs(idPart) == 11) idTrack = 1; | |
228 | if(TMath::Abs(idPart) == 111 || TMath::Abs(idPart) == 211) idTrack = 2; | |
229 | if(TMath::Abs(idPart) == 311 || TMath::Abs(idPart) == 321) idTrack = 3; | |
230 | if(TMath::Abs(idPart) == 2212) idTrack = 4; | |
231 | ||
232 | if(idTrack > 0 && fRecTrack > 0){ | |
233 | // Check if track should be reconstructed | |
234 | if((fRecTrack == 1 && idTrack == 1) || | |
235 | (fRecTrack == 2 && idTrack == 2) || | |
236 | (fRecTrack == 3 && idTrack == 3) || | |
237 | (fRecTrack == 4 && idTrack == 4) || | |
238 | fRecTrack == 100 ) { | |
239 | // Tracks are reconstructed | |
240 | ErrorMatrix(idTrack,pT,eta, v11, v22, v33, v12, v13, v23, iFlag); | |
241 | ||
242 | // Calculate and smear track parameters | |
243 | Double_t lambda, cosLambda, pTotal,pInverse; | |
244 | Double_t pInverseSmea, lambdaSmea, phiSmea; | |
245 | Double_t a1, a2, a3, b2, b3, c3; | |
246 | Double_t rn1, rn2, rn3; | |
247 | ||
248 | lambda = kPiHalf -2.0*TMath::ATan(TMath::Exp(-eta)); | |
249 | cosLambda = TMath::Cos(lambda); | |
250 | pTotal = pT/cosLambda; | |
251 | pInverse = 1.0/pTotal; | |
252 | ||
253 | a1 = TMath::Sqrt(v11); | |
254 | if(a1 == 0.){ | |
255 | a2 = 0; | |
256 | a3 = 0; | |
257 | }else{ | |
258 | a2 = v12/a1; | |
259 | a3 = v13/a1; | |
260 | } | |
261 | b2 = TMath::Sqrt(v22-a2*a2); | |
262 | if(b2 == 0.){ | |
263 | b3 = 0; | |
264 | }else{ | |
265 | b3 = (v23 - a2*a3)/b2; | |
266 | } | |
267 | c3 = TMath::Sqrt(v33 - a3*a3 -b3*b3); | |
268 | rn1 = gRandom->Gaus(0,1); | |
269 | rn2 = gRandom->Gaus(0,1); | |
270 | rn3 = gRandom->Gaus(0,1); | |
271 | ||
272 | pInverseSmea = pInverse + a1*rn1; | |
273 | lambdaSmea = lambda + a2*rn1 + b2*rn2; | |
274 | phiSmea = phi + a3*rn1 + b3*rn2 + c3*rn3; | |
275 | ||
276 | // Fill control histograms | |
277 | if(idTrack == 1){ | |
278 | fResID11->Fill((pInverseSmea-pInverse)/pInverse); | |
279 | fResID12->Fill((lambdaSmea-lambda)/lambda); | |
280 | fResID13->Fill((phiSmea-phi)/phi); | |
281 | } | |
282 | else if(idTrack == 2){ | |
283 | fResID21->Fill((pInverseSmea-pInverse)/pInverse); | |
284 | fResID22->Fill((lambdaSmea-lambda)/lambda); | |
285 | fResID23->Fill((phiSmea-phi)/phi); | |
286 | } | |
287 | else if(idTrack == 3){ | |
288 | fResID31->Fill((pInverseSmea-pInverse)/pInverse); | |
289 | fResID32->Fill((lambdaSmea-lambda)/lambda); | |
290 | fResID33->Fill((phiSmea-phi)/phi); | |
291 | } | |
292 | else if(idTrack == 4){ | |
293 | fResID41->Fill((pInverseSmea-pInverse)/pInverse); | |
294 | fResID42->Fill((lambdaSmea-lambda)/lambda); | |
295 | fResID43->Fill((phiSmea-phi)/phi); | |
296 | } | |
297 | }else{ | |
298 | // Tracks are not reconstructed | |
299 | v11=0.; | |
300 | v12=0.; | |
301 | v13=0.; | |
302 | v22=0.; | |
303 | v23=0.; | |
304 | v33=0.; | |
305 | iFlag=0; | |
306 | } | |
307 | // Store resolution variables to AliFTrack ClonesArray | |
308 | AddTrack(idTrack, charge, pT, eta, phi, v11, v22, v33, v12, v13, v23, iFlag); | |
309 | printf("%10s%10d%20.5e%20.5e%20.5e%20.5e%20.5e%20.5e%20.5e%20.5e%20.5e%20.5e%10d\n", | |
310 | "Track",idTrack,charge,pT,eta,phi,v11, v22, v33, v12, v13, v23, iFlag); | |
311 | } | |
312 | ||
313 | } | |
314 | } | |
315 | } | |
316 | ||
4c155701 | 317 | //_______________________________________________________________________ |
318 | void AliFTrackMaker::Copy(TObject &) const | |
319 | { | |
320 | Fatal("Copy","Not implemented!\n"); | |
321 | } | |
322 | ||
65a39007 | 323 | //_____________________________________________________________________________ |
324 | void AliFTrackMaker::Finish() | |
325 | { | |
326 | // For TestJob only | |
327 | if(gAliFast->TestTrack() == 1){ | |
328 | /* | |
329 | // Draw test histograms | |
330 | TCanvas *c1 = new TCanvas("c1"," ",200,10,600,480); | |
331 | c1->Divide(2,3); | |
332 | c1->cd(1); fResID1Test->Draw(); | |
333 | c1->cd(2); fResID2Test->Draw(); | |
334 | c1->cd(3); fResID3Test->Draw(); | |
335 | c1->cd(4); fResID4Test->Draw(); | |
336 | c1->cd(5); fResID5Test->Draw(); | |
337 | c1->Update(); | |
338 | // Store TestRes.eps file | |
339 | c1->Print("TestRes.eps"); | |
340 | */ | |
341 | // Store histograms on file | |
342 | TFile f2("TestRes.root","RECREATE","Test Res.f"); | |
343 | fResID1Test->Write(); | |
344 | fResID2Test->Write(); | |
345 | fResID3Test->Write(); | |
346 | fResID4Test->Write(); | |
347 | fResID5Test->Write(); | |
348 | f2.Close(); | |
349 | } | |
350 | } | |
351 | //_____________________________________________________________________________ | |
352 | void AliFTrackMaker::ErrorMatrix(Int_t idTrack, Double_t pT, Double_t eta, | |
353 | Double_t &v11, Double_t &v22, Double_t &v33, Double_t &v12, Double_t &v13, Double_t &v23, | |
354 | Int_t &iFlag) | |
355 | { | |
356 | /////////////////////////////////////////////// | |
357 | //idTrack track type input // | |
358 | //pT transverse mom input // | |
359 | //lambda deep angle input // | |
360 | //v11,v22,v23 error matrix output // | |
361 | //v12,v13,v23 output // | |
362 | //iFlag output // | |
363 | /////////////////////////////////////////////// | |
364 | ||
365 | AliFDet *detector = gAliFast->Detector(); | |
366 | Int_t nDet = detector->NDet(); | |
367 | Int_t nDetActive = detector->NDetActive(); | |
368 | Int_t nTwice = nDetActive + nDetActive; | |
369 | ||
370 | Double_t rTrack, rTrackInverse, pTotal, pInverse, diffPInverse; | |
371 | Double_t safety; | |
372 | Double_t cosLambda, tanLambda, diffLambda; | |
373 | Double_t rDet; | |
374 | ||
375 | Double_t hh0[kNMaxDet2][kNMaxDet2], hhi0[kNMaxDet2][kNMaxDet2]; | |
376 | Double_t hh1[kNMaxDet2][kNMaxDet2], hhi1[kNMaxDet2][kNMaxDet2]; | |
377 | Double_t dhhiOverPInverse[kNMaxDet2][kNMaxDet2]; | |
378 | Double_t dhhiOverLambda[kNMaxDet2][kNMaxDet2]; | |
379 | Double_t a1[kNMaxDet2][kNMaxDet2], a2[kNMaxDet2][kNMaxDet2]; | |
380 | Double_t a0PInverse[kNMaxDet2]; | |
381 | Double_t a0Lambda[kNMaxDet2]; | |
382 | Double_t a0Phi[kNMaxDet2]; | |
383 | ||
384 | Double_t vF11, vF12, vF13, vF22, vF23, vF33, d1, d2, d3, det; | |
385 | Int_t idet, icyl, im, in; | |
386 | Double_t phiHalf; | |
387 | Double_t lambda; | |
388 | ||
389 | lambda = kPiHalf -2.0*TMath::ATan(TMath::Exp(-eta)); | |
390 | rTrack = detector->ConstMag()*pT; | |
391 | safety = 10.0; | |
392 | if(2.0*rTrack < (detector->RDet(nDet) + safety)){ | |
393 | iFlag = 0; | |
394 | v11 = 0; | |
395 | v22 = 0; | |
396 | v33 = 0; | |
397 | v12 = 0; | |
398 | v13 = 0; | |
399 | v23 = 0; | |
400 | return; | |
401 | } | |
402 | iFlag = 1; | |
403 | cosLambda = TMath::Cos(lambda); | |
404 | pTotal = pT/cosLambda; | |
405 | pInverse = 1.0/pTotal; | |
406 | diffPInverse = pInverse*1.0e-5; | |
407 | diffLambda = 1.0e-4; | |
408 | ||
409 | // Compute likelihood and derivatives | |
410 | ||
411 | LogLikelyhood(idTrack, pInverse, lambda); | |
412 | for(icyl=1; icyl<nTwice+1; icyl++){ | |
413 | for(im=1; im<nTwice+1; im++){ | |
414 | hh0[icyl][im] = HH(icyl,im); | |
415 | hhi0[icyl][im] = HHI(icyl,im); | |
416 | } | |
417 | } | |
418 | LogLikelyhood(idTrack, pInverse+diffPInverse,lambda); | |
419 | for(icyl=1; icyl<nTwice+1; icyl++){ | |
420 | for(im=1; im<nTwice+1; im++){ | |
421 | hh1[icyl][im] = HH(icyl,im); | |
422 | hhi1[icyl][im] = HHI(icyl,im); | |
423 | } | |
424 | } | |
425 | for(icyl=1; icyl<nTwice+1; icyl++){ | |
426 | for(im=1; im<icyl+1; im++){ | |
427 | dhhiOverPInverse[icyl][im] = (hhi1[icyl][im]-hhi0[icyl][im])/diffPInverse; | |
428 | } | |
429 | } | |
430 | LogLikelyhood(idTrack, pInverse, lambda+diffLambda); | |
431 | for(icyl=1; icyl<nTwice+1; icyl++){ | |
432 | for(im=1; im<nTwice+1; im++){ | |
433 | hh1[icyl][im] = HH(icyl,im); | |
434 | hhi1[icyl][im] = HHI(icyl,im); | |
435 | } | |
436 | } | |
437 | for(icyl=1; icyl<nTwice+1; icyl++){ | |
438 | for(im=1; im<icyl+1; im++){ | |
439 | dhhiOverLambda[icyl][im] = (hhi1[icyl][im]-hhi0[icyl][im])/diffLambda; | |
440 | } | |
441 | } | |
442 | ||
443 | // Compute additional derivatives | |
444 | rTrackInverse = 1.0/rTrack; | |
445 | tanLambda = TMath::Tan(lambda); | |
446 | icyl = 0; | |
447 | for(idet=1; idet<nDet+1;idet++){ | |
448 | if(detector->IFlagDet(idet) > 0){ | |
449 | icyl = icyl + 1; | |
450 | rDet = detector->RDet(idet); | |
451 | phiHalf = TMath::ASin(0.5*rDet*rTrackInverse); | |
452 | Double_t rHelp = rDet / | |
453 | (2.0 * TMath::Sqrt(1.0-(0.5 *rDet*rTrackInverse)* | |
454 | (0.5 *rDet*rTrackInverse))); | |
455 | a0PInverse[icyl] = - rDet* rHelp | |
456 | /(detector->ConstMag()*cosLambda); | |
457 | a0Lambda[icyl] = - rDet* rHelp | |
458 | * tanLambda * rTrackInverse; | |
459 | a0Phi[icyl] = rDet; | |
460 | a0PInverse[nDetActive+icyl] = 2.0 * tanLambda | |
461 | *rTrack*(rHelp-rTrack*phiHalf) | |
462 | /(detector->ConstMag()*cosLambda); | |
463 | a0Lambda[nDetActive+icyl] = 2.0 * ( rHelp*tanLambda*tanLambda | |
464 | + rTrack*phiHalf); | |
465 | a0Phi[nDetActive+icyl] = 0.0 ; | |
466 | } | |
467 | } | |
468 | ||
469 | // Compute information matrix | |
470 | ||
471 | vF11=0.0; | |
472 | vF12=0.0; | |
473 | vF13=0.0; | |
474 | vF22=0.0; | |
475 | vF23=0.0; | |
476 | vF33=0.0; | |
477 | for(icyl=1; icyl<nTwice+1; icyl++){ | |
478 | d1=0.0; | |
479 | d2=0.0; | |
480 | d3=0.0; | |
481 | for(im=1; im < icyl+1; im++){ | |
482 | d1 = d1 + hhi0[icyl][im]*a0PInverse[im]; | |
483 | d2 = d2 + hhi0[icyl][im]*a0Lambda[im]; | |
484 | d3 = d3 + hhi0[icyl][im]*a0Phi[im]; | |
485 | } | |
486 | vF11 =vF11 + d1*d1; | |
487 | vF12 =vF12 + d1*d2; | |
488 | vF13 =vF13 + d1*d3; | |
489 | vF22 =vF22 + d2*d2; | |
490 | vF23 =vF23 + d2*d3; | |
491 | vF33 =vF33 + d3*d3; | |
492 | } | |
493 | for(icyl=1; icyl<nTwice+1; icyl++){ | |
494 | for(im=1; im<icyl+1; im++){ | |
495 | a1[icyl][im] = 0; | |
496 | a2[icyl][im] = 0; | |
497 | for(in=im; in<icyl+1;in++){ | |
498 | a1[icyl][im]=a1[icyl][im]+dhhiOverPInverse[icyl][in]*hh0[im][in]; | |
499 | a2[icyl][im]=a2[icyl][im]+dhhiOverLambda[icyl][in]*hh0[im][in]; | |
500 | } | |
501 | vF11=vF11+a1[icyl][im]*a1[icyl][im]; | |
502 | vF12=vF12+a1[icyl][im]*a2[icyl][im]; | |
503 | vF22=vF22+a2[icyl][im]*a2[icyl][im]; | |
504 | } | |
505 | vF11=vF11+a1[icyl][icyl]*a1[icyl][icyl]; | |
506 | vF12=vF12+a1[icyl][icyl]*a2[icyl][icyl]; | |
507 | vF22=vF22+a2[icyl][icyl]*a2[icyl][icyl]; | |
508 | } | |
509 | ||
510 | // Invert information matrix | |
511 | ||
512 | det=( vF11*vF22 - vF12*vF12 ) *vF33 + (vF12*vF23 - vF13*vF22)*vF13 | |
513 | + (vF12*vF13 - vF11*vF23)*vF23; | |
514 | ||
515 | v11 = (vF22*vF33 - vF23*vF23)/det; | |
516 | v22 = (vF11*vF33 - vF13*vF13)/det; | |
517 | v33 = (vF11*vF22 - vF12*vF12)/det; | |
518 | v12 = (vF13*vF23 - vF12*vF33)/det; | |
519 | v13 = (vF12*vF23 - vF13*vF22)/det; | |
520 | v23 = (vF12*vF13 - vF11*vF23)/det; | |
521 | ||
522 | } | |
523 | //_____________________________________________________________________________// | |
524 | void AliFTrackMaker::LogLikelyhood(Int_t idTrack, Double_t pInverse,Double_t lambda) | |
525 | { | |
526 | /////////////////////////////////////////////// | |
527 | //hh ?? output // | |
528 | //hhi ?? output // | |
529 | //idTrack track type input // | |
530 | //pInverse inverse momentum input // | |
531 | //lambda polar angle of track input // | |
532 | /////////////////////////////////////////////// | |
533 | ||
534 | ||
535 | AliFDet *detector = gAliFast->Detector(); | |
536 | Int_t nDet = detector->NDet(); | |
537 | Int_t nDetActive = detector->NDetActive(); | |
538 | Int_t nTwice = nDetActive + nDetActive; | |
539 | ||
540 | Double_t rDet, rDetSQ; | |
541 | Int_t idet, icyl, im, imc; | |
542 | Double_t cosLambda, tanLambda, pTotal, pT, rTrack, rTrackSQ; | |
543 | Double_t beta, overPBeta, rTrackInv, thickCorr, temp1, temp2; | |
544 | Double_t partMassSQ; | |
545 | Double_t aShelp[kNMaxDet2], dShelp[kNMaxDet2]; | |
546 | Double_t projXVXT[kNMaxDet2],projYVXT[kNMaxDet2], projZVXT[kNMaxDet2]; | |
547 | Double_t proj[kNMaxDet2][kNMaxDet2]; | |
548 | Double_t erroScatt[kNMaxDet2], variance[kNMaxDet2][kNMaxDet2]; | |
549 | Double_t erroSQ[kNMaxDet2]; | |
550 | Double_t hh[kNMaxDet2][kNMaxDet2]; | |
551 | Double_t hhi[kNMaxDet2][kNMaxDet2]; | |
552 | Double_t errorVX, errorVY, errorVZ; | |
553 | ||
554 | cosLambda = TMath::Cos(lambda); | |
555 | tanLambda = TMath::Tan(lambda); | |
556 | pTotal = 1.0/pInverse; | |
557 | pT = pTotal * cosLambda; | |
558 | rTrack = detector->ConstMag() * pTotal * cosLambda; | |
559 | rTrackSQ = rTrack * rTrack; | |
560 | partMassSQ= ParticleMass(idTrack)*ParticleMass(idTrack); | |
561 | beta = pTotal / TMath::Sqrt(partMassSQ+pTotal*pTotal); | |
562 | overPBeta = 1./(pTotal*beta); | |
563 | rTrackInv = 1./rTrack; | |
564 | errorVX = detector->ErrorVertexX(); | |
565 | errorVY = detector->ErrorVertexY(); | |
566 | errorVZ = detector->ErrorVertexZ(); | |
567 | ||
568 | ||
569 | erroScatt[0]=0.0; | |
570 | erroScatt[1]=0.0; | |
571 | for(idet=1; idet < nDet; idet++){ | |
572 | thickCorr = detector->ThickDet(idet)/TMath::Sqrt(cosLambda* | |
573 | TMath::Sqrt(1.0-0.25*(detector->RDetSQ(idet)/rTrackSQ))); | |
574 | if(detector->IFlagGas(idet) == 0){ | |
575 | thickCorr = thickCorr * (1.3266 + 0.076 * TMath::Log(thickCorr));} | |
576 | thickCorr = overPBeta * thickCorr; | |
577 | erroScatt[idet+1]=thickCorr*thickCorr; | |
578 | } | |
579 | ||
580 | ||
581 | icyl = 0; | |
582 | for(idet=1; idet<nDet+1; idet++){ | |
583 | rDet = detector->RDet(idet); | |
584 | rDetSQ = rDet*rDet; | |
585 | dShelp[idet] = TMath::Sqrt(4.0*rTrackSQ-rDetSQ); | |
586 | aShelp[idet] = TMath::ASin(rDet/(rTrack+rTrack)); | |
587 | if(detector->IFlagDet(idet) > 0) { | |
588 | icyl = icyl + 1; | |
589 | projXVXT[icyl] = rDet * rTrackInv; | |
590 | projXVXT[nDetActive+icyl] = -tanLambda; | |
591 | temp1 = (rTrackSQ + rTrackSQ - rDetSQ)/dShelp[idet]; | |
592 | temp2 = rDet/dShelp[idet]; | |
593 | projYVXT[icyl] = temp1*rTrackInv; | |
594 | projYVXT[nDetActive+icyl] = tanLambda * temp2; | |
595 | projZVXT[icyl] = 0.0; | |
596 | projZVXT[nDetActive+icyl] = 1.0; | |
597 | proj[icyl][1] = 0.0; | |
598 | proj[nDetActive+icyl][0] = 0.0; | |
599 | proj[nDetActive+icyl][nDet] = 0.0; | |
600 | proj[icyl][nDet] = 0.0; | |
601 | for(im=2; im<idet+1; im++){ | |
602 | proj[icyl][im]= (( rDet | |
603 | *(rTrackSQ+rTrackSQ-detector->RDetSQ(im-1)) | |
604 | - detector->RDet(im-1)*temp1*dShelp[im-1]) | |
605 | /((rTrackSQ + rTrackSQ)*cosLambda)); | |
606 | proj[nDetActive+icyl][im]= 0.5 * detector->RDet(im-1) | |
607 | * rTrackInv | |
608 | * tanLambda * (detector->RDet(im-1) | |
609 | - dShelp[im-1]*temp2)/cosLambda; | |
610 | proj[nDetActive+icyl][nDet+im]= (rTrack+rTrack) * (aShelp[idet] - aShelp[im-1]) | |
611 | + ( rDet*dShelp[im-1]-detector->RDet(im-1)*dShelp[idet]) | |
612 | * dShelp[im-1] * tanLambda * tanLambda | |
613 | / (dShelp[idet] * (rTrack+rTrack)); | |
614 | proj[icyl][nDet+im]= tanLambda | |
615 | * (rDet*detector->RDet(im-1)*dShelp[im-1] | |
616 | / (rTrackSQ+rTrackSQ) | |
617 | - (rDetSQ + detector->RDetSQ(im-1) | |
618 | - rDetSQ * detector->RDetSQ(im-1) | |
619 | / (rTrackSQ+rTrackSQ)) | |
620 | / dShelp[idet]); | |
621 | } | |
622 | for(im=idet+1; im < nDet+1; im++){ | |
623 | proj[icyl][im] = 0.0; | |
624 | proj[nDetActive+icyl][im] = 0.0; | |
625 | proj[nDetActive+icyl][nDet+im] = 0.0; | |
626 | proj[icyl][nDet+im] = 0.0; | |
627 | } | |
628 | if(detector->IFlagDet(idet) == 1){ | |
629 | erroSQ[icyl] = detector->ErrorRPhi(idet); | |
630 | erroSQ[nDetActive+icyl] = detector->ErrorZ(idet); | |
631 | }else{ | |
632 | TPCResolution(pT, rDet, lambda); | |
633 | erroSQ[icyl] = SigmaRPhiSQ(); | |
634 | erroSQ[nDetActive+icyl] = SigmaZSQ(); | |
635 | } | |
636 | erroSQ[icyl] = erroSQ[icyl] + detector->ErrorR(idet)*temp2*temp2; | |
637 | erroSQ[nDetActive+icyl] = erroSQ[nDetActive+icyl] | |
638 | + detector->ErrorR(idet)*tanLambda*tanLambda; | |
639 | } | |
640 | } | |
641 | for(icyl=1; icyl<nTwice+1; icyl++){ | |
642 | for(im=1; im<icyl+1; im++){ | |
643 | variance[icyl][im]= | |
644 | projXVXT[icyl]*projXVXT[im]*errorVX | |
645 | +projYVXT[icyl]*projYVXT[im]*errorVY | |
646 | +projZVXT[icyl]*projZVXT[im]*errorVZ; | |
647 | for(imc=1; imc<nDet+1; imc++){ | |
648 | variance[icyl][im]=variance[icyl][im] | |
649 | +(proj[icyl][imc]*proj[im][imc] | |
650 | + proj[icyl][nDet+imc]*proj[im][nDet+imc]) | |
651 | * erroScatt[imc]; | |
652 | } | |
653 | } | |
654 | variance[icyl][icyl] = variance[icyl][icyl]+erroSQ[icyl]; | |
655 | } | |
656 | ||
657 | for(icyl=1; icyl<nTwice+1; icyl++){ | |
658 | for(im=icyl; im<nTwice+1; im++){ | |
659 | hh[icyl][im]=variance[im][icyl]; | |
660 | for(imc=1; imc<icyl;imc++){ | |
661 | hh[icyl][im]=hh[icyl][im]-hh[imc][icyl]*hh[imc][im]; | |
662 | } | |
663 | if(im == icyl){ | |
664 | hh[icyl][im] = TMath::Sqrt(hh[icyl][im]); | |
665 | } else { | |
666 | hh[icyl][im] = hh[icyl][im]/hh[icyl][icyl]; | |
667 | } | |
668 | } | |
669 | } | |
670 | ||
671 | for(icyl=1; icyl<nTwice+1; icyl++){ | |
672 | hhi[icyl][icyl] = 1.0 / hh[icyl][icyl]; | |
673 | for(im=1; im<icyl; im++){ | |
674 | hhi[icyl][im] = 0.0; | |
675 | for(imc=im; imc<icyl; imc++){ | |
676 | hhi[icyl][im] = hhi[icyl][im]-hh[imc][icyl]*hhi[imc][im]; | |
677 | } | |
678 | hhi[icyl][im] = hhi[icyl][im]*hhi[icyl][icyl]; | |
679 | } | |
680 | } | |
681 | ||
682 | for(icyl=1; icyl<nTwice+1; icyl++){ | |
683 | for(im=1; im<nTwice+1; im++){ | |
684 | SetHH(icyl,im,hh[icyl][im]); | |
685 | SetHHI(icyl,im,hhi[icyl][im]); | |
686 | } | |
687 | } | |
688 | ||
689 | ||
690 | } | |
691 | //_____________________________________________________________________________ | |
692 | // translation of routine tpc_resolution of res.f | |
693 | //_____________________________________________________________________________ | |
694 | void AliFTrackMaker::TPCResolution(Double_t pTransv, Double_t radiPad, Double_t lambda) | |
695 | { | |
696 | /////////////////////////////////////////////// | |
697 | //sigmaRPhiSQ resolution in r-phi output // | |
698 | //sigmaZSQ resolution in z output // | |
699 | //pTransv transverse momentum input // | |
700 | //radiPad radius of pad row input // | |
701 | //lambda polar angle of track input // | |
702 | // | |
703 | //units: cm, GeV/c, radian // | |
704 | //parametrisation of TPC resolution // | |
705 | //version of 03.07.1995 // | |
706 | //source: Marek Kowalski, Karel Safarik // | |
707 | /////////////////////////////////////////////// | |
708 | ||
709 | Double_t aRCoeff=0.41818e-2; | |
710 | Double_t bRCoeff=0.17460e-4; | |
711 | Double_t cRCoeff=0.30993e-8; | |
712 | Double_t dRCoeff=0.41061e-6; | |
713 | Double_t aZCoeff=0.39610e-2; | |
714 | Double_t bZCoeff=0.22443e-4; | |
715 | Double_t cZCoeff=0.51504e-1; | |
716 | ||
717 | Double_t sigmaRPhiSQ; | |
718 | Double_t sigmaZSQ; | |
719 | ||
720 | sigmaRPhiSQ = aRCoeff - bRCoeff * radiPad * TMath::Tan(lambda)+ | |
721 | (cRCoeff * (radiPad/pTransv) + dRCoeff) * radiPad/pTransv; | |
722 | ||
723 | sigmaZSQ = aZCoeff - bZCoeff * radiPad * TMath::Tan(lambda)+ | |
724 | cZCoeff * TMath::Tan(lambda)*TMath::Tan(lambda); | |
725 | ||
726 | if(sigmaRPhiSQ < 1.0e-6 ) sigmaRPhiSQ = 1.0e-6; | |
727 | if(sigmaZSQ < 1.0e-6 ) sigmaZSQ = 1.0e-6; | |
728 | ||
729 | sigmaRPhiSQ = (TMath::Sqrt(sigmaRPhiSQ) + 0.005) | |
730 | * (TMath::Sqrt(sigmaRPhiSQ) + 0.005); | |
731 | sigmaZSQ = (TMath::Sqrt(sigmaZSQ) + 0.005) | |
732 | * (TMath::Sqrt(sigmaZSQ) + 0.005); | |
733 | ||
734 | SetSigmaRPhiSQ(sigmaRPhiSQ); | |
735 | SetSigmaZSQ(sigmaZSQ); | |
736 | ||
737 | ||
738 | } | |
739 | ||
65a39007 | 740 | //----------------------------------------------------------------------------- |
4c155701 | 741 | Double_t AliFTrackMaker::ParticleMass(Int_t idTrack) const |
65a39007 | 742 | { |
4c155701 | 743 | // |
744 | // returns the mass given particle ID | |
745 | // | |
65a39007 | 746 | Double_t mass = 0.0; |
747 | ||
748 | if(idTrack == 2){ mass = fPionMass;} | |
749 | else if(idTrack == 3){ mass = fKaonMass;} | |
750 | else if(idTrack == 1) {mass = fElectronMass;} | |
751 | else if(idTrack == 4) {mass = fProtonMass;} | |
752 | ||
753 | return mass; | |
754 | ||
755 | } | |
756 | ||
4c155701 | 757 | //____________________________________________________________________________ |
65a39007 | 758 | Double_t AliFTrackMaker::Rapidity(Double_t pt, Double_t pz) |
759 | { | |
4c155701 | 760 | // |
761 | // returns the rapidity given particle pT, pz | |
762 | // Compute rapidity | |
65a39007 | 763 | |
764 | Double_t etalog = TMath::Log((TMath::Sqrt(pt*pt + pz*pz) + TMath::Abs(pz))/pt); | |
765 | if (pz < 0 ) return -TMath::Abs(etalog); | |
766 | else return TMath::Abs(etalog); | |
767 | } | |
768 | ||
769 | //_____________________________________________________________________________ | |
770 | // returns the phi angle given particle px, py | |
771 | //----------------------------------------------------------------------------- | |
772 | Double_t AliFTrackMaker::Angle(Double_t x, Double_t y) | |
773 | { | |
774 | // Compute phi angle of particle | |
775 | // ... this is a copy of function ULANGL | |
776 | // .. sign(a,b) = -abs(a) if b < 0 | |
777 | // = abs(a) if b >= 0 | |
778 | ||
779 | Double_t angle = 0; | |
780 | Double_t r = TMath::Sqrt(x*x + y*y); | |
781 | if (r < 1e-20) return angle; | |
782 | if (TMath::Abs(x)/r < 0.8) { | |
783 | angle = TMath::Sign((Double_t)TMath::Abs(TMath::ACos(x/r)), y); | |
784 | } else { | |
785 | angle = TMath::ASin(y/r); | |
786 | if (x < 0 ) { | |
787 | if(angle >= 0) angle = kPi - angle; | |
788 | else angle = -kPi - angle; | |
789 | } | |
790 | } | |
791 | return angle; | |
792 | } | |
793 | //_____________________________________________________________________________ | |
794 | Int_t AliFTrackMaker::Charge(Int_t kf) | |
795 | { | |
796 | //...this is copy of function LUCHGE | |
797 | //...Purpose: to give three times the charge for a particle/parton. | |
798 | ||
799 | static Int_t kchg[500] = { -1,2,-1,2,-1,2,-1,2,0,0,-3,0,-3,0,-3,0,-3,0, | |
800 | 0,0,0,0,0,3,0,0,0,0,0,0,0,0,0,3,0,0,3,0,-1,0,0,0,0,0,0,0,0,0,0,0, | |
801 | 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, | |
802 | 2,-1,2,-1,2,3,0,0,0,0,0,0,0,0,0,0,0,3,0,3,3,0,3,0,3,0,3,0,0,0,0,0, | |
803 | 0,0,0,0,0,3,0,3,3,0,3,0,3,0,3,0,0,0,0,0,0,0,0,0,0,3,0,3,3,0,3,0,3, | |
804 | 0,3,0,0,0,0,0,0,0,0,0,0,3,0,3,3,0,3,0,3,0,3,0,0,0,0,0,0,0,0,0,0,3, | |
805 | 0,3,3,0,3,0,3,0,3,0,0,0,0,0,0,0,0,0,0,3,0,3,3,0,3,0,3,0,3,0,0,0,0, | |
806 | 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, | |
807 | 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, | |
808 | 3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,3,0, | |
809 | 0,3,0,0,0,0,0,0,0,0,-3,0,0,0,0,0,0,0,0,3,0,-3,0,3,-3,0,0,0,3,6,0, | |
810 | 3,0,0,0,0,0,-3,0,3,-3,0,-3,0,0,0,0,-3,0,3,6,-3,0,3,-3,0,-3,0,3,6, | |
811 | 0,3,0,0,0,0,0,-3,0,3,-3,0,-3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, | |
812 | 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, | |
813 | 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, | |
814 | 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }; | |
815 | ||
816 | // extern integer kfcomp_(integer *); | |
817 | Int_t ipower; | |
818 | Int_t ret = 0; | |
819 | Int_t kfa = TMath::Abs(kf); | |
820 | Int_t kc = Compress(kfa); | |
821 | ||
822 | //...Initial values. Simple case of direct readout. | |
823 | if (kc == 0) { | |
824 | } else if (kfa <= 100 || kc <= 80 || kc > 100) { | |
825 | ret = kchg[kc-1]; | |
826 | ||
827 | // ...Construction from quark content for heavy meson, diquark, baryon. | |
828 | } else if (kfa/1000 % 10 == 0) { | |
829 | ipower = kfa/100 % 10; | |
830 | ret = (kchg[kfa / 100 % 10 - 1] - kchg[kfa/10 % 10 - 1])*Int_t(TMath::Power(-1, ipower)); | |
831 | } else if (kfa / 10 % 10 == 0) { | |
832 | ret = kchg[kfa/1000 % 10 - 1] + kchg[kfa/100 % 10 - 1]; | |
833 | } else { | |
834 | ret = kchg[kfa/1000 % 10 - 1] + kchg[kfa/100 % 10 - 1] + kchg[kfa/10 % 10 - 1]; | |
835 | } | |
836 | ||
837 | // ...Add on correct sign. | |
838 | if (kf > 0) return ret; | |
839 | else return -ret; | |
840 | } | |
841 | //_____________________________________________________________________________ | |
842 | Int_t AliFTrackMaker::Compress(Int_t kf) | |
843 | { | |
844 | //...this is copy of function LUCOMP | |
845 | //...Purpose: to compress the standard KF codes for use in mass and decay | |
846 | //...arrays; also to check whether a given code actually is defined. | |
847 | // from BLOCK LUDATA | |
848 | static Int_t kchg[500] = { 1,1,1,1,1,1,1,1,0,0,1,1,1,1,1,1,1,1,0,0,0,0, | |
849 | 0,1,0,0,0,0,0,0,0,0,0,1,0,0,1,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, | |
850 | 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,1,1,1,1, | |
851 | 1,1,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,1, | |
852 | 1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,0,0,0,0, | |
853 | 0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1, | |
854 | 1,1,1,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0, | |
855 | 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, | |
856 | 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0, | |
857 | 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,1,0,1,1,0,0,0,0, | |
858 | 0,0,1,0,1,1,0,0,0,0,0,0,1,1,0,1,1,1,1,1,0,1,1,1,1,1,1,0,0,0,0,1,1, | |
859 | 1,1,1,1,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,1,1,1,1,1, | |
860 | 1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, | |
861 | 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, | |
862 | 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, | |
863 | 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }; | |
864 | static Int_t kftab[25] = { 211,111,221,311,321,130,310,213,113,223,313, | |
865 | 323,2112,2212,210,2110,2210,110,220,330,440,30443,30553,0,0 }; | |
866 | static Int_t kctab[25] = { 101,111,112,102,103,221,222,121,131,132,122, | |
867 | 123,332,333,281,282,283,284,285,286,287,231,235,0,0 }; | |
868 | ||
869 | Int_t ret = 0; | |
870 | Int_t kfla, kflb, kflc, kflr, kfls, kfa, ikf; | |
871 | ||
872 | kfa = TMath::Abs(kf); | |
873 | //...Simple cases: direct translation or table. | |
874 | if (kfa == 0 || kfa >= 100000) { | |
875 | return ret; | |
876 | } else if (kfa <= 100) { | |
877 | ret = kfa; | |
878 | if (kf < 0 && kchg[kfa - 1] == 0) ret = 0; | |
879 | return ret; | |
880 | } else { | |
881 | for (ikf = 1; ikf <= 23; ++ikf) { | |
882 | if (kfa == kftab[ikf-1]) { | |
883 | ret = kctab[ikf-1]; | |
884 | if (kf < 0 && kchg[ret-1] == 0) ret = 0; | |
885 | return ret; | |
886 | } | |
887 | } | |
888 | } | |
889 | // ...Subdivide KF code into constituent pieces. | |
890 | kfla = kfa / 1000%10; | |
891 | kflb = kfa / 100%10; | |
892 | kflc = kfa / 10%10; | |
893 | kfls = kfa%10; | |
894 | kflr = kfa / 10000%10; | |
895 | // ...Mesons. | |
896 | if (kfa - kflr*10000 < 1000) { | |
897 | if (kflb == 0 || kflb == 9 || kflc == 0 || kflc == 9) { | |
898 | } else if (kflb < kflc) { | |
899 | } else if (kf < 0 && kflb == kflc) { | |
900 | } else if (kflb == kflc) { | |
901 | if (kflr == 0 && kfls == 1) { ret = kflb + 110; | |
902 | } else if (kflr == 0 && kfls == 3) { ret = kflb + 130; | |
903 | } else if (kflr == 1 && kfls == 3) { ret = kflb + 150; | |
904 | } else if (kflr == 1 && kfls == 1) { ret = kflb + 170; | |
905 | } else if (kflr == 2 && kfls == 3) { ret = kflb + 190; | |
906 | } else if (kflr == 0 && kfls == 5) { ret = kflb + 210; | |
907 | } | |
908 | } else if (kflb <= 5) { | |
909 | if (kflr == 0 && kfls == 1) { ret = (kflb-1)*(kflb-2)/2 + 100 + kflc; | |
910 | } else if (kflr == 0 && kfls == 3) { ret = (kflb-1)*(kflb-2)/2 + 120 + kflc; | |
911 | } else if (kflr == 1 && kfls == 3) { ret = (kflb-1)*(kflb-2)/2 + 140 + kflc; | |
912 | } else if (kflr == 1 && kfls == 1) { ret = (kflb-1)*(kflb-2)/2 + 160 + kflc; | |
913 | } else if (kflr == 2 && kfls == 3) { ret = (kflb-1)*(kflb-2)/2 + 180 + kflc; | |
914 | } else if (kflr == 0 && kfls == 5) { ret = (kflb-1)*(kflb-2)/2 + 200 + kflc; | |
915 | } | |
916 | } else if (kfls == 1 && kflr <= 1 || kfls == 3 && kflr <= 2 || kfls == 5 && kflr == 0) { | |
917 | ret = kflb + 80; | |
918 | } | |
919 | // ...Diquarks. | |
920 | } else if ((kflr == 0 || kflr == 1) && kflc == 0) { | |
921 | if (kfls != 1 && kfls != 3) { | |
922 | } else if (kfla == 9 || kflb == 0 || kflb == 9) { | |
923 | } else if (kfla < kflb) { | |
924 | } else if (kfls == 1 && kfla == kflb) { | |
925 | } else { ret = 90; | |
926 | } | |
927 | // ...Spin 1/2 baryons. | |
928 | } else if (kflr == 0 && kfls == 2) { | |
929 | if (kfla == 9 || kflb == 0 || kflb == 9 || kflc == 9) { | |
930 | } else if (kfla <= kflc || kfla < kflb) { | |
931 | } else if (kfla >= 6 || kflb >= 4 || kflc >= 4) { | |
932 | ret = kfla + 80; | |
933 | } else if (kflb < kflc) { | |
934 | ret = (kfla+1)*kfla*(kfla-1)/6 + 300 + kflc*(kflc-1)/2 + kflb; | |
935 | } else { | |
936 | ret = (kfla+1)*kfla*(kfla-1)/6 + 330 + kflb*(kflb-1)/2 + kflc; | |
937 | } | |
938 | // ...Spin 3/2 baryons. | |
939 | } else if (kflr == 0 && kfls == 4) { | |
940 | if (kfla == 9 || kflb == 0 || kflb == 9 || kflc == 9) { | |
941 | } else if (kfla < kflb || kflb < kflc) { | |
942 | } else if (kfla >= 6 || kflb >= 4) { | |
943 | ret = kfla + 80; | |
944 | } else { | |
945 | ret = (kfla+1)*kfla*(kfla-1) / 6 + 360 + kflb*(kflb -1) / 2 + kflc; | |
946 | } | |
947 | } | |
948 | return ret; | |
949 | } | |
950 | ||
65a39007 | 951 | //_____________________________________________________________________________ |
952 | void AliFTrackMaker::MakeTest(Int_t n) | |
953 | { | |
4c155701 | 954 | // |
955 | // TEST JOB: Calculate tracks resolution | |
956 | // | |
65a39007 | 957 | Double_t v11, v22, v33, v12, v13, v23; |
958 | Int_t iFlag; | |
959 | Int_t idTrack; | |
960 | Double_t pTStart, pT, eta; | |
961 | ||
962 | Double_t sumDPop,sumDDip,sumDPhi; | |
963 | Double_t isum,fm; | |
964 | Double_t pTotal,partMassSQ,beta,lambda; | |
965 | Double_t dPop,dLop,dDip,dPhi,rho12,rho13,rho23; | |
29894e14 | 966 | Double_t dPPStrag,dPPTot=0; |
65a39007 | 967 | // Double_t resol1[1001][11],resol2[10001][11],resol3[1001][11], |
968 | // resol4[1001][11],resol5[10001][11] | |
969 | Double_t store1[1001],store2[10001],store3[1001], | |
970 | store4[1001],store5[10001]; | |
971 | ||
972 | ||
973 | idTrack = 2; | |
974 | pTStart = 0.07; | |
975 | for(Int_t istep=1; istep<n; istep++){ | |
976 | if(istep < 100 && istep > 20) istep = istep -1 + 5; | |
977 | if(istep < 500 && istep > 100) istep = istep -1 + 25; | |
978 | if(istep <1000 && istep > 500) istep = istep -1 + 100; | |
979 | pT = pTStart + 0.01*istep; | |
980 | eta = - 0.044; | |
981 | sumDPop = 0; | |
982 | sumDDip = 0; | |
983 | sumDPhi = 0; | |
984 | isum = 0; | |
985 | for(Int_t in=1; in<11; in++){ | |
986 | eta = eta + 0.088; | |
987 | lambda = kPiHalf -2.0*TMath::ATan(TMath::Exp(-eta)); | |
988 | pTotal = pT / TMath::Cos(lambda); | |
989 | if(idTrack == 1){ | |
990 | dPPStrag = 0.055 /pT;} | |
991 | else{ | |
992 | partMassSQ = ParticleMass(idTrack)*ParticleMass(idTrack); | |
993 | beta = pTotal/ TMath::Sqrt(pTotal*pTotal + partMassSQ); | |
994 | dPPStrag = 0.04/(pT*TMath::Power(beta,2.6666666667)); | |
995 | } | |
996 | ErrorMatrix(idTrack,pT,eta, v11, v22, v33, v12, v13, v23, iFlag); | |
997 | if(iFlag == 1){ | |
998 | dLop = TMath::Sqrt(v11); | |
999 | dDip = TMath::Sqrt(v22); | |
1000 | dPhi = TMath::Sqrt(v33); | |
1001 | rho12 = v12/(dLop*dDip); | |
1002 | rho13 = v13/(dLop*dPhi); | |
1003 | rho23 = v23/(dDip*dPhi); | |
1004 | dPop = 100. *dLop * pTotal; | |
1005 | dDip = 1000. * dDip; | |
1006 | dPhi = 1000. * dPhi; | |
1007 | dPPTot = TMath::Sqrt(dPop*dPop + dPPStrag*dPPStrag); | |
1008 | // resol1[istep][in] = dPop; | |
1009 | // resol2[istep][in] = dDip; | |
1010 | // resol3[istep][in] = dPhi; | |
1011 | // resol4[istep][in] = dPPTot; | |
1012 | // resol5[istep][in] = dPPStrag; | |
1013 | sumDPop = sumDPop + dPop; | |
1014 | sumDDip = sumDDip + dDip; | |
1015 | sumDPhi = sumDPhi + dPhi; | |
1016 | isum = isum + 1;} | |
1017 | else{ | |
1018 | printf("%20s %10.5f %10.5f %20s\n","pT,eta",pT,eta,"cannot smear"); | |
1019 | } | |
1020 | } | |
1021 | if(isum > 0){ | |
1022 | dPop = sumDPop/isum; | |
1023 | dDip = sumDDip/isum; | |
1024 | dPhi = sumDPhi/isum; | |
1025 | dPPTot = TMath::Sqrt(dPop*dPop + dPPStrag*dPPStrag);} | |
1026 | else{ | |
1027 | dPop = 0; | |
1028 | dDip = 0; | |
1029 | dPhi = 0; | |
1030 | } | |
1031 | store1[istep] = dPop; | |
1032 | store2[istep] = dDip; | |
1033 | store3[istep] = dPhi; | |
1034 | store4[istep] = dPPTot; | |
1035 | store5[istep] = dPPStrag; | |
1036 | if(istep > 20 ){ | |
1037 | Int_t im = 5; | |
1038 | if(istep > 100) {im = 25;} | |
1039 | if(istep > 500) {im = 100;} | |
1040 | fm = 1./(1.*im); | |
1041 | for(Int_t ist=1; ist<im; ist++){ | |
1042 | // for(Int_t in=1; in < 11; in++){ | |
1043 | // resol1[istep-im+ist][in] = resol1[istep-im][in]+ | |
1044 | // ist*fm*(resol1[istep][in]-resol1[istep-im][in]); | |
1045 | // resol2[istep-im+ist][in] = resol2[istep-im][in]+ | |
1046 | // ist*fm*(resol2[istep][in]-resol2[istep-im][in]); | |
1047 | // resol3[istep-im+ist][in] = resol3[istep-im][in]+ | |
1048 | // ist*fm*(resol3[istep][in]-resol3[istep-im][in]); | |
1049 | // resol4[istep-im+ist][in] = resol4[istep-im][in]+ | |
1050 | // ist*fm*(resol4[istep][in]-resol4[istep-im][in]); | |
1051 | // resol5[istep-im+ist][in] = resol5[istep-im][in]+ | |
1052 | // ist*fm*(resol5[istep][in]-resol5[istep-im][in]); | |
1053 | // } | |
1054 | store1[istep-im+ist]=store1[istep-im]+ | |
1055 | ist*fm*(store1[istep]-store1[istep-im]); | |
1056 | store2[istep-im+ist]=store2[istep-im]+ | |
1057 | ist*fm*(store2[istep]-store2[istep-im]); | |
1058 | store3[istep-im+ist]=store3[istep-im]+ | |
1059 | ist*fm*(store3[istep]-store3[istep-im]); | |
1060 | store4[istep-im+ist]=store4[istep-im]+ | |
1061 | ist*fm*(store4[istep]-store4[istep-im]); | |
1062 | store5[istep-im+ist]=store5[istep-im]+ | |
1063 | ist*fm*(store5[istep]-store5[istep-im]); | |
1064 | // Fill control histograms | |
1065 | fResID1Test->Fill(pTStart + 0.01*(istep-im+ist),store1[istep-im+ist]); | |
1066 | fResID2Test->Fill(pTStart + 0.01*(istep-im+ist),store2[istep-im+ist]); | |
1067 | fResID3Test->Fill(pTStart + 0.01*(istep-im+ist),store3[istep-im+ist]); | |
1068 | fResID4Test->Fill(pTStart + 0.01*(istep-im+ist),store4[istep-im+ist]); | |
1069 | fResID5Test->Fill(pTStart + 0.01*(istep-im+ist),store5[istep-im+ist]); | |
1070 | } | |
1071 | printf("%10s %10d %20.15f %20.15f %20.15f %20.15f %20.15f \n", | |
1072 | "TestTrack:",istep,store1[istep],store2[istep],store3[istep], | |
1073 | store4[istep],store5[istep]); | |
1074 | } else { | |
1075 | printf("%10s %10d %20.15f %20.15f %20.15f %20.15f %20.15f \n", | |
1076 | "TestTrack:",istep,store1[istep],store2[istep],store3[istep], | |
1077 | store4[istep],store5[istep]); | |
1078 | fResID1Test->Fill(pT,store1[istep]); | |
1079 | fResID2Test->Fill(pT,store2[istep]); | |
1080 | fResID3Test->Fill(pT,store3[istep]); | |
1081 | fResID4Test->Fill(pT,store4[istep]); | |
1082 | fResID5Test->Fill(pT,store5[istep]); | |
1083 | } | |
1084 | } | |
1085 | } | |
1086 | //_____________________________________________________________________________ |