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