null alignment file
[u/mrichter/AliRoot.git] / ALIFAST / AliFTrackMaker.cxx
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 <TFile.h>
49 #include <TH1.h>
50 #include <TMath.h>
51 #include <TParticle.h>
52 #include <TRandom.h>
53
54 #include "AliFDet.h"
55 #include "AliFTrack.h"
56 #include "AliFTrackMaker.h"
57 #include "AliFast.h"
58 #include "AliMC.h"
59
60 static const Double_t kPi       = TMath::Pi();
61 static const Double_t k2Pi      = 2*kPi;
62 static const Double_t kPiHalf   = kPi/2.;
63 extern  AliFast * gAliFast;
64 ClassImp(AliFTrackMaker)
65
66 //_____________________________________________________________________________
67 AliFTrackMaker::AliFTrackMaker()
68 {
69   //
70   // Default constructor
71   //
72   fNTracks = 0;
73   fResID1Test = 0;
74   fResID2Test = 0;
75   fResID3Test = 0;
76   fResID4Test = 0;
77   fResID5Test = 0;
78 }
79
80 //_____________________________________________________________________________
81 AliFTrackMaker::AliFTrackMaker(const char *name, const char *title)
82                  :AliFMaker(name,title)
83 {
84   //
85   // Standard Setters for tracks
86   //
87    fFruits     = new TClonesArray("AliFTrack",100, kFALSE);
88    fBranchName = "Tracks";
89    fNTracks    = 0;
90 // Please, how to do this optionally ??!!!
91    Save();
92 }
93
94 //_______________________________________________________________________
95 AliFTrackMaker::AliFTrackMaker(const AliFTrackMaker& aftmk):
96   AliFMaker(aftmk)
97 {
98   //
99   // Copy constructor for AliRun
100   //
101   aftmk.Copy(*this);
102 }
103
104 //_____________________________________________________________________________
105 AliFTrackMaker::~AliFTrackMaker()
106 {
107   //
108   // Dummy constructor
109   //
110 }
111
112 //_____________________________________________________________________________
113 AliFTrack *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 //_____________________________________________________________________________
136 void AliFTrackMaker::Clear(Option_t *option)
137 {
138    //Reset Track Maker
139
140    fNTracks = 0;
141    AliFMaker::Clear(option);
142 }
143
144 //_____________________________________________________________________________
145 void AliFTrackMaker::Draw(Option_t *)
146 {
147 //    Dummy Draw
148
149 }
150
151 //_____________________________________________________________________________
152 void 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
194 void AliFTrackMaker::Make()
195 {
196   //
197   // Calculate track and its resolution
198   //
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
212      Int_t idPart, idTrack;
213      Double_t  charge, pT, eta, phi;
214      TParticle *part;
215      Int_t  nparticles = gAlice->GetMCApp()->GetNtrack();
216      printf("%10s%10d\n","nparticles",nparticles);
217      for(Int_t ind=0;ind<nparticles;ind++) {       
218        part = gAlice->GetMCApp()->Particle(ind);
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
317 //_______________________________________________________________________
318 void AliFTrackMaker::Copy(TObject &) const
319 {
320   Fatal("Copy","Not implemented!\n");
321 }
322
323 //_____________________________________________________________________________
324 void 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 //_____________________________________________________________________________
352 void 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 //_____________________________________________________________________________//
524 void 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 //_____________________________________________________________________________
694 void 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
740 //-----------------------------------------------------------------------------
741 Double_t AliFTrackMaker::ParticleMass(Int_t idTrack) const
742 {
743   //
744   // returns the mass given particle ID 
745   //
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
757 //____________________________________________________________________________
758 Double_t AliFTrackMaker::Rapidity(Double_t pt, Double_t pz)
759 {
760   //
761   // returns the rapidity given particle pT, pz 
762   //   Compute rapidity
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 //-----------------------------------------------------------------------------
772 Double_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 //_____________________________________________________________________________
794 Int_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 //_____________________________________________________________________________
842 Int_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
951 //_____________________________________________________________________________
952 void AliFTrackMaker::MakeTest(Int_t n)
953 {
954   //
955   // TEST JOB: Calculate tracks resolution
956   //
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;
966   Double_t dPPStrag,dPPTot=0;
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 //_____________________________________________________________________________