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