1 //////////////////////////////////////////////////////////////////////////
3 // AliFast TrackMaker class. //
6 //////////////////////////////////////////////////////////////////////////
7 // ---------------------------------------------------------------------//
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?? //
16 // C++ in AliFast framework: Elzbieta Richter-Was and Yiota Foka //
17 // following general structure od Makers in //
18 // ATLFast by R. Brun. //
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 //
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. //
39 //////////////////////////////////////////////////////////////////////////
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 )
48 #include <TParticle.h>
59 //#include "AliFMCMaker.h"
60 #include "AliFTrackMaker.h"
61 #include "AliFTrack.h"
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)
70 //_____________________________________________________________________________
71 AliFTrackMaker::AliFTrackMaker()
76 //_____________________________________________________________________________
77 AliFTrackMaker::AliFTrackMaker(const char *name, const char *title)
78 :AliFMaker(name,title)
80 // Default Setters for tracks
82 fFruits = new TClonesArray("AliFTrack",100, kFALSE);
83 fBranchName = "Tracks";
85 // Please, how to do this optionally ??!!!
89 //_____________________________________________________________________________
90 AliFTrackMaker::~AliFTrackMaker()
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)
101 // Add a new track to the list of tracks
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.
113 TClonesArray &tracks = *(TClonesArray*)fFruits;
114 return new(tracks[fNTracks++]) AliFTrack(code,charge,pT,eta,phi,
115 v11, v22, v33, v12, v13, v23, iFlag);
118 //_____________________________________________________________________________
119 void AliFTrackMaker::Clear(Option_t *option)
124 AliFMaker::Clear(option);
127 //_____________________________________________________________________________
128 void AliFTrackMaker::Draw(Option_t *)
134 //_____________________________________________________________________________
135 void AliFTrackMaker::Init()
137 //Create control histograms
138 if(gAliFast->TestTrack() == 0){
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);
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);
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);
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);
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);
166 //Set particle masses
172 //Switch on/off tracks reconstruction
177 //_____________________________________________________________________________
178 // Calculate track and its resolution
179 //_____________________________________________________________________________
180 void AliFTrackMaker::Make()
182 Double_t v11, v22, v33, v12, v13, v23;
187 // Check if it is a TestJob
188 if(gAliFast->TestTrack() == 1){
192 // Run production job
193 // Get pointers to Particles arrays and TClonesArray
195 TClonesArray *particles = gAliFast->Particles();
196 Int_t idPart, idTrack;
197 Double_t charge, pT, eta, phi;
199 Int_t nparticles = particles->GetEntriesFast();
200 printf("%10s%10d\n","nparticles",nparticles);
201 for(Int_t ind=0;ind<nparticles;ind++) {
202 part = (TParticle*)particles->UncheckedAt(ind);
203 idPart = part->GetPdgCode();
204 charge = part->GetPDG()->Charge();
208 printf("%10s%10d%20.5e%20.5e%20.5e%20.5e\n","Particle",idPart,charge,pT,eta,phi);
209 // Check convention for tracks reconstruction
211 if(TMath::Abs(idPart) == 11) idTrack = 1;
212 if(TMath::Abs(idPart) == 111 || TMath::Abs(idPart) == 211) idTrack = 2;
213 if(TMath::Abs(idPart) == 311 || TMath::Abs(idPart) == 321) idTrack = 3;
214 if(TMath::Abs(idPart) == 2212) idTrack = 4;
216 if(idTrack > 0 && fRecTrack > 0){
217 // Check if track should be reconstructed
218 if((fRecTrack == 1 && idTrack == 1) ||
219 (fRecTrack == 2 && idTrack == 2) ||
220 (fRecTrack == 3 && idTrack == 3) ||
221 (fRecTrack == 4 && idTrack == 4) ||
223 // Tracks are reconstructed
224 ErrorMatrix(idTrack,pT,eta, v11, v22, v33, v12, v13, v23, iFlag);
226 // Calculate and smear track parameters
227 Double_t lambda, cosLambda, pTotal,pInverse;
228 Double_t pInverseSmea, lambdaSmea, phiSmea;
229 Double_t a1, a2, a3, b2, b3, c3;
230 Double_t rn1, rn2, rn3;
232 lambda = kPiHalf -2.0*TMath::ATan(TMath::Exp(-eta));
233 cosLambda = TMath::Cos(lambda);
234 pTotal = pT/cosLambda;
235 pInverse = 1.0/pTotal;
237 a1 = TMath::Sqrt(v11);
245 b2 = TMath::Sqrt(v22-a2*a2);
249 b3 = (v23 - a2*a3)/b2;
251 c3 = TMath::Sqrt(v33 - a3*a3 -b3*b3);
252 rn1 = gRandom->Gaus(0,1);
253 rn2 = gRandom->Gaus(0,1);
254 rn3 = gRandom->Gaus(0,1);
256 pInverseSmea = pInverse + a1*rn1;
257 lambdaSmea = lambda + a2*rn1 + b2*rn2;
258 phiSmea = phi + a3*rn1 + b3*rn2 + c3*rn3;
260 // Fill control histograms
262 fResID11->Fill((pInverseSmea-pInverse)/pInverse);
263 fResID12->Fill((lambdaSmea-lambda)/lambda);
264 fResID13->Fill((phiSmea-phi)/phi);
266 else if(idTrack == 2){
267 fResID21->Fill((pInverseSmea-pInverse)/pInverse);
268 fResID22->Fill((lambdaSmea-lambda)/lambda);
269 fResID23->Fill((phiSmea-phi)/phi);
271 else if(idTrack == 3){
272 fResID31->Fill((pInverseSmea-pInverse)/pInverse);
273 fResID32->Fill((lambdaSmea-lambda)/lambda);
274 fResID33->Fill((phiSmea-phi)/phi);
276 else if(idTrack == 4){
277 fResID41->Fill((pInverseSmea-pInverse)/pInverse);
278 fResID42->Fill((lambdaSmea-lambda)/lambda);
279 fResID43->Fill((phiSmea-phi)/phi);
282 // Tracks are not reconstructed
291 // Store resolution variables to AliFTrack ClonesArray
292 AddTrack(idTrack, charge, pT, eta, phi, v11, v22, v33, v12, v13, v23, iFlag);
293 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",
294 "Track",idTrack,charge,pT,eta,phi,v11, v22, v33, v12, v13, v23, iFlag);
301 //_____________________________________________________________________________
302 void AliFTrackMaker::Finish()
305 if(gAliFast->TestTrack() == 1){
307 // Draw test histograms
308 TCanvas *c1 = new TCanvas("c1"," ",200,10,600,480);
310 c1->cd(1); fResID1Test->Draw();
311 c1->cd(2); fResID2Test->Draw();
312 c1->cd(3); fResID3Test->Draw();
313 c1->cd(4); fResID4Test->Draw();
314 c1->cd(5); fResID5Test->Draw();
316 // Store TestRes.eps file
317 c1->Print("TestRes.eps");
319 // Store histograms on file
320 TFile f2("TestRes.root","RECREATE","Test Res.f");
321 fResID1Test->Write();
322 fResID2Test->Write();
323 fResID3Test->Write();
324 fResID4Test->Write();
325 fResID5Test->Write();
329 //_____________________________________________________________________________
330 void AliFTrackMaker::ErrorMatrix(Int_t idTrack, Double_t pT, Double_t eta,
331 Double_t &v11, Double_t &v22, Double_t &v33, Double_t &v12, Double_t &v13, Double_t &v23,
334 ///////////////////////////////////////////////
335 //idTrack track type input //
336 //pT transverse mom input //
337 //lambda deep angle input //
338 //v11,v22,v23 error matrix output //
339 //v12,v13,v23 output //
341 ///////////////////////////////////////////////
343 AliFDet *detector = gAliFast->Detector();
344 Int_t nDet = detector->NDet();
345 Int_t nDetActive = detector->NDetActive();
346 Int_t nTwice = nDetActive + nDetActive;
348 Double_t rTrack, rTrackInverse, pTotal, pInverse, diffPInverse;
350 Double_t cosLambda, tanLambda, diffLambda;
353 Double_t hh0[kNMaxDet2][kNMaxDet2], hhi0[kNMaxDet2][kNMaxDet2];
354 Double_t hh1[kNMaxDet2][kNMaxDet2], hhi1[kNMaxDet2][kNMaxDet2];
355 Double_t dhhiOverPInverse[kNMaxDet2][kNMaxDet2];
356 Double_t dhhiOverLambda[kNMaxDet2][kNMaxDet2];
357 Double_t a1[kNMaxDet2][kNMaxDet2], a2[kNMaxDet2][kNMaxDet2];
358 Double_t a0PInverse[kNMaxDet2];
359 Double_t a0Lambda[kNMaxDet2];
360 Double_t a0Phi[kNMaxDet2];
362 Double_t vF11, vF12, vF13, vF22, vF23, vF33, d1, d2, d3, det;
363 Int_t idet, icyl, im, in;
367 lambda = kPiHalf -2.0*TMath::ATan(TMath::Exp(-eta));
368 rTrack = detector->ConstMag()*pT;
370 if(2.0*rTrack < (detector->RDet(nDet) + safety)){
381 cosLambda = TMath::Cos(lambda);
382 pTotal = pT/cosLambda;
383 pInverse = 1.0/pTotal;
384 diffPInverse = pInverse*1.0e-5;
387 // Compute likelihood and derivatives
389 LogLikelyhood(idTrack, pInverse, lambda);
390 for(icyl=1; icyl<nTwice+1; icyl++){
391 for(im=1; im<nTwice+1; im++){
392 hh0[icyl][im] = HH(icyl,im);
393 hhi0[icyl][im] = HHI(icyl,im);
396 LogLikelyhood(idTrack, pInverse+diffPInverse,lambda);
397 for(icyl=1; icyl<nTwice+1; icyl++){
398 for(im=1; im<nTwice+1; im++){
399 hh1[icyl][im] = HH(icyl,im);
400 hhi1[icyl][im] = HHI(icyl,im);
403 for(icyl=1; icyl<nTwice+1; icyl++){
404 for(im=1; im<icyl+1; im++){
405 dhhiOverPInverse[icyl][im] = (hhi1[icyl][im]-hhi0[icyl][im])/diffPInverse;
408 LogLikelyhood(idTrack, pInverse, lambda+diffLambda);
409 for(icyl=1; icyl<nTwice+1; icyl++){
410 for(im=1; im<nTwice+1; im++){
411 hh1[icyl][im] = HH(icyl,im);
412 hhi1[icyl][im] = HHI(icyl,im);
415 for(icyl=1; icyl<nTwice+1; icyl++){
416 for(im=1; im<icyl+1; im++){
417 dhhiOverLambda[icyl][im] = (hhi1[icyl][im]-hhi0[icyl][im])/diffLambda;
421 // Compute additional derivatives
422 rTrackInverse = 1.0/rTrack;
423 tanLambda = TMath::Tan(lambda);
425 for(idet=1; idet<nDet+1;idet++){
426 if(detector->IFlagDet(idet) > 0){
428 rDet = detector->RDet(idet);
429 phiHalf = TMath::ASin(0.5*rDet*rTrackInverse);
430 Double_t rHelp = rDet /
431 (2.0 * TMath::Sqrt(1.0-(0.5 *rDet*rTrackInverse)*
432 (0.5 *rDet*rTrackInverse)));
433 a0PInverse[icyl] = - rDet* rHelp
434 /(detector->ConstMag()*cosLambda);
435 a0Lambda[icyl] = - rDet* rHelp
436 * tanLambda * rTrackInverse;
438 a0PInverse[nDetActive+icyl] = 2.0 * tanLambda
439 *rTrack*(rHelp-rTrack*phiHalf)
440 /(detector->ConstMag()*cosLambda);
441 a0Lambda[nDetActive+icyl] = 2.0 * ( rHelp*tanLambda*tanLambda
443 a0Phi[nDetActive+icyl] = 0.0 ;
447 // Compute information matrix
455 for(icyl=1; icyl<nTwice+1; icyl++){
459 for(im=1; im < icyl+1; im++){
460 d1 = d1 + hhi0[icyl][im]*a0PInverse[im];
461 d2 = d2 + hhi0[icyl][im]*a0Lambda[im];
462 d3 = d3 + hhi0[icyl][im]*a0Phi[im];
471 for(icyl=1; icyl<nTwice+1; icyl++){
472 for(im=1; im<icyl+1; im++){
475 for(in=im; in<icyl+1;in++){
476 a1[icyl][im]=a1[icyl][im]+dhhiOverPInverse[icyl][in]*hh0[im][in];
477 a2[icyl][im]=a2[icyl][im]+dhhiOverLambda[icyl][in]*hh0[im][in];
479 vF11=vF11+a1[icyl][im]*a1[icyl][im];
480 vF12=vF12+a1[icyl][im]*a2[icyl][im];
481 vF22=vF22+a2[icyl][im]*a2[icyl][im];
483 vF11=vF11+a1[icyl][icyl]*a1[icyl][icyl];
484 vF12=vF12+a1[icyl][icyl]*a2[icyl][icyl];
485 vF22=vF22+a2[icyl][icyl]*a2[icyl][icyl];
488 // Invert information matrix
490 det=( vF11*vF22 - vF12*vF12 ) *vF33 + (vF12*vF23 - vF13*vF22)*vF13
491 + (vF12*vF13 - vF11*vF23)*vF23;
493 v11 = (vF22*vF33 - vF23*vF23)/det;
494 v22 = (vF11*vF33 - vF13*vF13)/det;
495 v33 = (vF11*vF22 - vF12*vF12)/det;
496 v12 = (vF13*vF23 - vF12*vF33)/det;
497 v13 = (vF12*vF23 - vF13*vF22)/det;
498 v23 = (vF12*vF13 - vF11*vF23)/det;
501 //_____________________________________________________________________________//
502 void AliFTrackMaker::LogLikelyhood(Int_t idTrack, Double_t pInverse,Double_t lambda)
504 ///////////////////////////////////////////////
507 //idTrack track type input //
508 //pInverse inverse momentum input //
509 //lambda polar angle of track input //
510 ///////////////////////////////////////////////
513 AliFDet *detector = gAliFast->Detector();
514 Int_t nDet = detector->NDet();
515 Int_t nDetActive = detector->NDetActive();
516 Int_t nTwice = nDetActive + nDetActive;
518 Double_t rDet, rDetSQ;
519 Int_t idet, icyl, im, imc;
520 Double_t cosLambda, tanLambda, pTotal, pT, rTrack, rTrackSQ;
521 Double_t beta, overPBeta, rTrackInv, thickCorr, temp1, temp2;
523 Double_t aShelp[kNMaxDet2], dShelp[kNMaxDet2];
524 Double_t projXVXT[kNMaxDet2],projYVXT[kNMaxDet2], projZVXT[kNMaxDet2];
525 Double_t proj[kNMaxDet2][kNMaxDet2];
526 Double_t erroScatt[kNMaxDet2], variance[kNMaxDet2][kNMaxDet2];
527 Double_t erroSQ[kNMaxDet2];
528 Double_t hh[kNMaxDet2][kNMaxDet2];
529 Double_t hhi[kNMaxDet2][kNMaxDet2];
530 Double_t errorVX, errorVY, errorVZ;
532 cosLambda = TMath::Cos(lambda);
533 tanLambda = TMath::Tan(lambda);
534 pTotal = 1.0/pInverse;
535 pT = pTotal * cosLambda;
536 rTrack = detector->ConstMag() * pTotal * cosLambda;
537 rTrackSQ = rTrack * rTrack;
538 partMassSQ= ParticleMass(idTrack)*ParticleMass(idTrack);
539 beta = pTotal / TMath::Sqrt(partMassSQ+pTotal*pTotal);
540 overPBeta = 1./(pTotal*beta);
541 rTrackInv = 1./rTrack;
542 errorVX = detector->ErrorVertexX();
543 errorVY = detector->ErrorVertexY();
544 errorVZ = detector->ErrorVertexZ();
549 for(idet=1; idet < nDet; idet++){
550 thickCorr = detector->ThickDet(idet)/TMath::Sqrt(cosLambda*
551 TMath::Sqrt(1.0-0.25*(detector->RDetSQ(idet)/rTrackSQ)));
552 if(detector->IFlagGas(idet) == 0){
553 thickCorr = thickCorr * (1.3266 + 0.076 * TMath::Log(thickCorr));}
554 thickCorr = overPBeta * thickCorr;
555 erroScatt[idet+1]=thickCorr*thickCorr;
560 for(idet=1; idet<nDet+1; idet++){
561 rDet = detector->RDet(idet);
563 dShelp[idet] = TMath::Sqrt(4.0*rTrackSQ-rDetSQ);
564 aShelp[idet] = TMath::ASin(rDet/(rTrack+rTrack));
565 if(detector->IFlagDet(idet) > 0) {
567 projXVXT[icyl] = rDet * rTrackInv;
568 projXVXT[nDetActive+icyl] = -tanLambda;
569 temp1 = (rTrackSQ + rTrackSQ - rDetSQ)/dShelp[idet];
570 temp2 = rDet/dShelp[idet];
571 projYVXT[icyl] = temp1*rTrackInv;
572 projYVXT[nDetActive+icyl] = tanLambda * temp2;
573 projZVXT[icyl] = 0.0;
574 projZVXT[nDetActive+icyl] = 1.0;
576 proj[nDetActive+icyl][0] = 0.0;
577 proj[nDetActive+icyl][nDet] = 0.0;
578 proj[icyl][nDet] = 0.0;
579 for(im=2; im<idet+1; im++){
580 proj[icyl][im]= (( rDet
581 *(rTrackSQ+rTrackSQ-detector->RDetSQ(im-1))
582 - detector->RDet(im-1)*temp1*dShelp[im-1])
583 /((rTrackSQ + rTrackSQ)*cosLambda));
584 proj[nDetActive+icyl][im]= 0.5 * detector->RDet(im-1)
586 * tanLambda * (detector->RDet(im-1)
587 - dShelp[im-1]*temp2)/cosLambda;
588 proj[nDetActive+icyl][nDet+im]= (rTrack+rTrack) * (aShelp[idet] - aShelp[im-1])
589 + ( rDet*dShelp[im-1]-detector->RDet(im-1)*dShelp[idet])
590 * dShelp[im-1] * tanLambda * tanLambda
591 / (dShelp[idet] * (rTrack+rTrack));
592 proj[icyl][nDet+im]= tanLambda
593 * (rDet*detector->RDet(im-1)*dShelp[im-1]
594 / (rTrackSQ+rTrackSQ)
595 - (rDetSQ + detector->RDetSQ(im-1)
596 - rDetSQ * detector->RDetSQ(im-1)
597 / (rTrackSQ+rTrackSQ))
600 for(im=idet+1; im < nDet+1; im++){
601 proj[icyl][im] = 0.0;
602 proj[nDetActive+icyl][im] = 0.0;
603 proj[nDetActive+icyl][nDet+im] = 0.0;
604 proj[icyl][nDet+im] = 0.0;
606 if(detector->IFlagDet(idet) == 1){
607 erroSQ[icyl] = detector->ErrorRPhi(idet);
608 erroSQ[nDetActive+icyl] = detector->ErrorZ(idet);
610 TPCResolution(pT, rDet, lambda);
611 erroSQ[icyl] = SigmaRPhiSQ();
612 erroSQ[nDetActive+icyl] = SigmaZSQ();
614 erroSQ[icyl] = erroSQ[icyl] + detector->ErrorR(idet)*temp2*temp2;
615 erroSQ[nDetActive+icyl] = erroSQ[nDetActive+icyl]
616 + detector->ErrorR(idet)*tanLambda*tanLambda;
619 for(icyl=1; icyl<nTwice+1; icyl++){
620 for(im=1; im<icyl+1; im++){
622 projXVXT[icyl]*projXVXT[im]*errorVX
623 +projYVXT[icyl]*projYVXT[im]*errorVY
624 +projZVXT[icyl]*projZVXT[im]*errorVZ;
625 for(imc=1; imc<nDet+1; imc++){
626 variance[icyl][im]=variance[icyl][im]
627 +(proj[icyl][imc]*proj[im][imc]
628 + proj[icyl][nDet+imc]*proj[im][nDet+imc])
632 variance[icyl][icyl] = variance[icyl][icyl]+erroSQ[icyl];
635 for(icyl=1; icyl<nTwice+1; icyl++){
636 for(im=icyl; im<nTwice+1; im++){
637 hh[icyl][im]=variance[im][icyl];
638 for(imc=1; imc<icyl;imc++){
639 hh[icyl][im]=hh[icyl][im]-hh[imc][icyl]*hh[imc][im];
642 hh[icyl][im] = TMath::Sqrt(hh[icyl][im]);
644 hh[icyl][im] = hh[icyl][im]/hh[icyl][icyl];
649 for(icyl=1; icyl<nTwice+1; icyl++){
650 hhi[icyl][icyl] = 1.0 / hh[icyl][icyl];
651 for(im=1; im<icyl; im++){
653 for(imc=im; imc<icyl; imc++){
654 hhi[icyl][im] = hhi[icyl][im]-hh[imc][icyl]*hhi[imc][im];
656 hhi[icyl][im] = hhi[icyl][im]*hhi[icyl][icyl];
660 for(icyl=1; icyl<nTwice+1; icyl++){
661 for(im=1; im<nTwice+1; im++){
662 SetHH(icyl,im,hh[icyl][im]);
663 SetHHI(icyl,im,hhi[icyl][im]);
669 //_____________________________________________________________________________
670 // translation of routine tpc_resolution of res.f
671 //_____________________________________________________________________________
672 void AliFTrackMaker::TPCResolution(Double_t pTransv, Double_t radiPad, Double_t lambda)
674 ///////////////////////////////////////////////
675 //sigmaRPhiSQ resolution in r-phi output //
676 //sigmaZSQ resolution in z output //
677 //pTransv transverse momentum input //
678 //radiPad radius of pad row input //
679 //lambda polar angle of track input //
681 //units: cm, GeV/c, radian //
682 //parametrisation of TPC resolution //
683 //version of 03.07.1995 //
684 //source: Marek Kowalski, Karel Safarik //
685 ///////////////////////////////////////////////
687 Double_t aRCoeff=0.41818e-2;
688 Double_t bRCoeff=0.17460e-4;
689 Double_t cRCoeff=0.30993e-8;
690 Double_t dRCoeff=0.41061e-6;
691 Double_t aZCoeff=0.39610e-2;
692 Double_t bZCoeff=0.22443e-4;
693 Double_t cZCoeff=0.51504e-1;
695 Double_t sigmaRPhiSQ;
698 sigmaRPhiSQ = aRCoeff - bRCoeff * radiPad * TMath::Tan(lambda)+
699 (cRCoeff * (radiPad/pTransv) + dRCoeff) * radiPad/pTransv;
701 sigmaZSQ = aZCoeff - bZCoeff * radiPad * TMath::Tan(lambda)+
702 cZCoeff * TMath::Tan(lambda)*TMath::Tan(lambda);
704 if(sigmaRPhiSQ < 1.0e-6 ) sigmaRPhiSQ = 1.0e-6;
705 if(sigmaZSQ < 1.0e-6 ) sigmaZSQ = 1.0e-6;
707 sigmaRPhiSQ = (TMath::Sqrt(sigmaRPhiSQ) + 0.005)
708 * (TMath::Sqrt(sigmaRPhiSQ) + 0.005);
709 sigmaZSQ = (TMath::Sqrt(sigmaZSQ) + 0.005)
710 * (TMath::Sqrt(sigmaZSQ) + 0.005);
712 SetSigmaRPhiSQ(sigmaRPhiSQ);
713 SetSigmaZSQ(sigmaZSQ);
718 //_____________________________________________________________________________
719 // returns the mass given particle ID
720 //-----------------------------------------------------------------------------
721 Double_t AliFTrackMaker::ParticleMass(Int_t idTrack)
725 if(idTrack == 2){ mass = fPionMass;}
726 else if(idTrack == 3){ mass = fKaonMass;}
727 else if(idTrack == 1) {mass = fElectronMass;}
728 else if(idTrack == 4) {mass = fProtonMass;}
734 //_____________________________________________________________________________
735 // returns the rapidity given particle pT, pz
736 //-----------------------------------------------------------------------------
737 Double_t AliFTrackMaker::Rapidity(Double_t pt, Double_t pz)
741 Double_t etalog = TMath::Log((TMath::Sqrt(pt*pt + pz*pz) + TMath::Abs(pz))/pt);
742 if (pz < 0 ) return -TMath::Abs(etalog);
743 else return TMath::Abs(etalog);
746 //_____________________________________________________________________________
747 // returns the phi angle given particle px, py
748 //-----------------------------------------------------------------------------
749 Double_t AliFTrackMaker::Angle(Double_t x, Double_t y)
751 // Compute phi angle of particle
752 // ... this is a copy of function ULANGL
753 // .. sign(a,b) = -abs(a) if b < 0
754 // = abs(a) if b >= 0
757 Double_t r = TMath::Sqrt(x*x + y*y);
758 if (r < 1e-20) return angle;
759 if (TMath::Abs(x)/r < 0.8) {
760 angle = TMath::Sign((Double_t)TMath::Abs(TMath::ACos(x/r)), y);
762 angle = TMath::ASin(y/r);
764 if(angle >= 0) angle = kPi - angle;
765 else angle = -kPi - angle;
770 //_____________________________________________________________________________
771 Int_t AliFTrackMaker::Charge(Int_t kf)
773 //...this is copy of function LUCHGE
774 //...Purpose: to give three times the charge for a particle/parton.
776 static Int_t kchg[500] = { -1,2,-1,2,-1,2,-1,2,0,0,-3,0,-3,0,-3,0,-3,0,
777 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,
778 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
779 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,
780 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,
781 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,
782 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,
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 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
785 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,
786 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,
787 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,
788 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,
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,0,0,0,
791 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 };
793 // extern integer kfcomp_(integer *);
796 Int_t kfa = TMath::Abs(kf);
797 Int_t kc = Compress(kfa);
799 //...Initial values. Simple case of direct readout.
801 } else if (kfa <= 100 || kc <= 80 || kc > 100) {
804 // ...Construction from quark content for heavy meson, diquark, baryon.
805 } else if (kfa/1000 % 10 == 0) {
806 ipower = kfa/100 % 10;
807 ret = (kchg[kfa / 100 % 10 - 1] - kchg[kfa/10 % 10 - 1])*Int_t(TMath::Power(-1, ipower));
808 } else if (kfa / 10 % 10 == 0) {
809 ret = kchg[kfa/1000 % 10 - 1] + kchg[kfa/100 % 10 - 1];
811 ret = kchg[kfa/1000 % 10 - 1] + kchg[kfa/100 % 10 - 1] + kchg[kfa/10 % 10 - 1];
814 // ...Add on correct sign.
815 if (kf > 0) return ret;
818 //_____________________________________________________________________________
819 Int_t AliFTrackMaker::Compress(Int_t kf)
821 //...this is copy of function LUCOMP
822 //...Purpose: to compress the standard KF codes for use in mass and decay
823 //...arrays; also to check whether a given code actually is defined.
825 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,
826 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,
827 0,0,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,
828 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,
829 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,
830 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,
831 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,
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,0,0,0,0,0,0,
833 0,0,0,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,
834 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,
835 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,
836 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,
837 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,
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,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
840 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 };
841 static Int_t kftab[25] = { 211,111,221,311,321,130,310,213,113,223,313,
842 323,2112,2212,210,2110,2210,110,220,330,440,30443,30553,0,0 };
843 static Int_t kctab[25] = { 101,111,112,102,103,221,222,121,131,132,122,
844 123,332,333,281,282,283,284,285,286,287,231,235,0,0 };
847 Int_t kfla, kflb, kflc, kflr, kfls, kfa, ikf;
849 kfa = TMath::Abs(kf);
850 //...Simple cases: direct translation or table.
851 if (kfa == 0 || kfa >= 100000) {
853 } else if (kfa <= 100) {
855 if (kf < 0 && kchg[kfa - 1] == 0) ret = 0;
858 for (ikf = 1; ikf <= 23; ++ikf) {
859 if (kfa == kftab[ikf-1]) {
861 if (kf < 0 && kchg[ret-1] == 0) ret = 0;
866 // ...Subdivide KF code into constituent pieces.
867 kfla = kfa / 1000%10;
871 kflr = kfa / 10000%10;
873 if (kfa - kflr*10000 < 1000) {
874 if (kflb == 0 || kflb == 9 || kflc == 0 || kflc == 9) {
875 } else if (kflb < kflc) {
876 } else if (kf < 0 && kflb == kflc) {
877 } else if (kflb == kflc) {
878 if (kflr == 0 && kfls == 1) { ret = kflb + 110;
879 } else if (kflr == 0 && kfls == 3) { ret = kflb + 130;
880 } else if (kflr == 1 && kfls == 3) { ret = kflb + 150;
881 } else if (kflr == 1 && kfls == 1) { ret = kflb + 170;
882 } else if (kflr == 2 && kfls == 3) { ret = kflb + 190;
883 } else if (kflr == 0 && kfls == 5) { ret = kflb + 210;
885 } else if (kflb <= 5) {
886 if (kflr == 0 && kfls == 1) { ret = (kflb-1)*(kflb-2)/2 + 100 + kflc;
887 } else if (kflr == 0 && kfls == 3) { ret = (kflb-1)*(kflb-2)/2 + 120 + kflc;
888 } else if (kflr == 1 && kfls == 3) { ret = (kflb-1)*(kflb-2)/2 + 140 + kflc;
889 } else if (kflr == 1 && kfls == 1) { ret = (kflb-1)*(kflb-2)/2 + 160 + kflc;
890 } else if (kflr == 2 && kfls == 3) { ret = (kflb-1)*(kflb-2)/2 + 180 + kflc;
891 } else if (kflr == 0 && kfls == 5) { ret = (kflb-1)*(kflb-2)/2 + 200 + kflc;
893 } else if (kfls == 1 && kflr <= 1 || kfls == 3 && kflr <= 2 || kfls == 5 && kflr == 0) {
897 } else if ((kflr == 0 || kflr == 1) && kflc == 0) {
898 if (kfls != 1 && kfls != 3) {
899 } else if (kfla == 9 || kflb == 0 || kflb == 9) {
900 } else if (kfla < kflb) {
901 } else if (kfls == 1 && kfla == kflb) {
904 // ...Spin 1/2 baryons.
905 } else if (kflr == 0 && kfls == 2) {
906 if (kfla == 9 || kflb == 0 || kflb == 9 || kflc == 9) {
907 } else if (kfla <= kflc || kfla < kflb) {
908 } else if (kfla >= 6 || kflb >= 4 || kflc >= 4) {
910 } else if (kflb < kflc) {
911 ret = (kfla+1)*kfla*(kfla-1)/6 + 300 + kflc*(kflc-1)/2 + kflb;
913 ret = (kfla+1)*kfla*(kfla-1)/6 + 330 + kflb*(kflb-1)/2 + kflc;
915 // ...Spin 3/2 baryons.
916 } else if (kflr == 0 && kfls == 4) {
917 if (kfla == 9 || kflb == 0 || kflb == 9 || kflc == 9) {
918 } else if (kfla < kflb || kflb < kflc) {
919 } else if (kfla >= 6 || kflb >= 4) {
922 ret = (kfla+1)*kfla*(kfla-1) / 6 + 360 + kflb*(kflb -1) / 2 + kflc;
928 //_____________________________________________________________________________
929 // TEST JOB: Calculate tracks resolution
930 //_____________________________________________________________________________
931 void AliFTrackMaker::MakeTest(Int_t n)
933 Double_t v11, v22, v33, v12, v13, v23;
936 Double_t pTStart, pT, eta;
938 Double_t sumDPop,sumDDip,sumDPhi;
940 Double_t pTotal,partMassSQ,beta,lambda;
941 Double_t dPop,dLop,dDip,dPhi,rho12,rho13,rho23;
942 Double_t dPPStrag,dPPTot;
943 // Double_t resol1[1001][11],resol2[10001][11],resol3[1001][11],
944 // resol4[1001][11],resol5[10001][11]
945 Double_t store1[1001],store2[10001],store3[1001],
946 store4[1001],store5[10001];
951 for(Int_t istep=1; istep<n; istep++){
952 if(istep < 100 && istep > 20) istep = istep -1 + 5;
953 if(istep < 500 && istep > 100) istep = istep -1 + 25;
954 if(istep <1000 && istep > 500) istep = istep -1 + 100;
955 pT = pTStart + 0.01*istep;
961 for(Int_t in=1; in<11; in++){
963 lambda = kPiHalf -2.0*TMath::ATan(TMath::Exp(-eta));
964 pTotal = pT / TMath::Cos(lambda);
966 dPPStrag = 0.055 /pT;}
968 partMassSQ = ParticleMass(idTrack)*ParticleMass(idTrack);
969 beta = pTotal/ TMath::Sqrt(pTotal*pTotal + partMassSQ);
970 dPPStrag = 0.04/(pT*TMath::Power(beta,2.6666666667));
972 ErrorMatrix(idTrack,pT,eta, v11, v22, v33, v12, v13, v23, iFlag);
974 dLop = TMath::Sqrt(v11);
975 dDip = TMath::Sqrt(v22);
976 dPhi = TMath::Sqrt(v33);
977 rho12 = v12/(dLop*dDip);
978 rho13 = v13/(dLop*dPhi);
979 rho23 = v23/(dDip*dPhi);
980 dPop = 100. *dLop * pTotal;
983 dPPTot = TMath::Sqrt(dPop*dPop + dPPStrag*dPPStrag);
984 // resol1[istep][in] = dPop;
985 // resol2[istep][in] = dDip;
986 // resol3[istep][in] = dPhi;
987 // resol4[istep][in] = dPPTot;
988 // resol5[istep][in] = dPPStrag;
989 sumDPop = sumDPop + dPop;
990 sumDDip = sumDDip + dDip;
991 sumDPhi = sumDPhi + dPhi;
994 printf("%20s %10.5f %10.5f %20s\n","pT,eta",pT,eta,"cannot smear");
1000 dPhi = sumDPhi/isum;
1001 dPPTot = TMath::Sqrt(dPop*dPop + dPPStrag*dPPStrag);}
1007 store1[istep] = dPop;
1008 store2[istep] = dDip;
1009 store3[istep] = dPhi;
1010 store4[istep] = dPPTot;
1011 store5[istep] = dPPStrag;
1014 if(istep > 100) {im = 25;}
1015 if(istep > 500) {im = 100;}
1017 for(Int_t ist=1; ist<im; ist++){
1018 // for(Int_t in=1; in < 11; in++){
1019 // resol1[istep-im+ist][in] = resol1[istep-im][in]+
1020 // ist*fm*(resol1[istep][in]-resol1[istep-im][in]);
1021 // resol2[istep-im+ist][in] = resol2[istep-im][in]+
1022 // ist*fm*(resol2[istep][in]-resol2[istep-im][in]);
1023 // resol3[istep-im+ist][in] = resol3[istep-im][in]+
1024 // ist*fm*(resol3[istep][in]-resol3[istep-im][in]);
1025 // resol4[istep-im+ist][in] = resol4[istep-im][in]+
1026 // ist*fm*(resol4[istep][in]-resol4[istep-im][in]);
1027 // resol5[istep-im+ist][in] = resol5[istep-im][in]+
1028 // ist*fm*(resol5[istep][in]-resol5[istep-im][in]);
1030 store1[istep-im+ist]=store1[istep-im]+
1031 ist*fm*(store1[istep]-store1[istep-im]);
1032 store2[istep-im+ist]=store2[istep-im]+
1033 ist*fm*(store2[istep]-store2[istep-im]);
1034 store3[istep-im+ist]=store3[istep-im]+
1035 ist*fm*(store3[istep]-store3[istep-im]);
1036 store4[istep-im+ist]=store4[istep-im]+
1037 ist*fm*(store4[istep]-store4[istep-im]);
1038 store5[istep-im+ist]=store5[istep-im]+
1039 ist*fm*(store5[istep]-store5[istep-im]);
1040 // Fill control histograms
1041 fResID1Test->Fill(pTStart + 0.01*(istep-im+ist),store1[istep-im+ist]);
1042 fResID2Test->Fill(pTStart + 0.01*(istep-im+ist),store2[istep-im+ist]);
1043 fResID3Test->Fill(pTStart + 0.01*(istep-im+ist),store3[istep-im+ist]);
1044 fResID4Test->Fill(pTStart + 0.01*(istep-im+ist),store4[istep-im+ist]);
1045 fResID5Test->Fill(pTStart + 0.01*(istep-im+ist),store5[istep-im+ist]);
1047 printf("%10s %10d %20.15f %20.15f %20.15f %20.15f %20.15f \n",
1048 "TestTrack:",istep,store1[istep],store2[istep],store3[istep],
1049 store4[istep],store5[istep]);
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 fResID1Test->Fill(pT,store1[istep]);
1055 fResID2Test->Fill(pT,store2[istep]);
1056 fResID3Test->Fill(pT,store3[istep]);
1057 fResID4Test->Fill(pT,store4[istep]);
1058 fResID5Test->Fill(pT,store5[istep]);
1062 //_____________________________________________________________________________