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