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