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 | //_____________________________________________________________________________ |