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