Now the code is blind to the binning used for pulls, efficiencies etc.
[u/mrichter/AliRoot.git] / TPC / AliTPCtrackerParam.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /*
17 $Log$
18 Revision 1.7  2002/04/10 16:30:07  kowal2
19 logs added
20
21 */
22
23
24 /**************************************************************************
25  *                                                                        *
26  * This class builds AliTPCtrack objects from generated tracks to feed    *
27  * ITS tracking (V2). The AliTPCtrack is built from its first hit in      *
28  * the TPC. The track is assigned a Kalman-like covariance matrix         *
29  * depending on its pT and pseudorapidity and track parameters are        *
30  * smeared according to this covariance matrix.                           *
31  * Output file contains sorted tracks, ready for matching with ITS.       *
32  *                                                                        *
33  * For details:                                                           *
34  * http://www.pd.infn.it/alipd/talks/soft/adIII02/TPCtrackingParam.htm    *
35  *                                                                        *
36  * Test macro is: AliBarrelRec_TPCparam.C                                 *   
37  *                                                                        *
38  * 2002/04/28: Major upgrade of the class                                 *
39  * - Covariance matrices and pulls are now separeted for pions, kaons     *
40  *   and electrons.                                                       *
41  * - A parameterization of dE/dx in the TPC has been included and it is   *
42  *   used to assign a mass to each track according to a rough dE/dx PID.  *
43  * - All the "numbers" have been moved to the file with the DataBase,     *
44  *   they are read as objects of the class AliTPCkineGrid, and assigned   *
45  *   to data memebers of the class AliTPCtrackerParam.                    *
46  * - All the code necessary to create a BataBase has been included in     *
47  *   class (see the macro AliTPCtrackingParamDB.C for the details).       *
48  *                                                                        *
49  *  Origin: Andrea Dainese, Padova - e-mail: andrea.dainese@pd.infn.it    *
50  *                                                                        *
51  **************************************************************************/
52 #include <TChain.h>
53 #include <TF1.h>
54 #include <TGraph.h>
55 #include <TGraphErrors.h>
56 #include <TLegend.h>
57 #include <TLine.h>
58 #include <TMatrixD.h>
59 #include <TNtuple.h>
60 #include <TSystem.h>
61 #include "alles.h"
62 #include "AliGausCorr.h"
63 #include "AliKalmanTrack.h"
64 #include "AliMagF.h"
65 #include "AliMagFCM.h"
66 #include "AliTPCkineGrid.h"
67 #include "AliTPCtrack.h"
68 #include "AliTPCtrackerParam.h"
69
70 Double_t RegFunc(Double_t *x,Double_t *par) {
71 // This is the function used to regularize the covariance matrix
72   Double_t value = par[0]+par[1]/TMath::Power(x[0],par[2])/
73                    TMath::Power(x[1],par[3]);
74   return value;
75 }
76 Double_t FitRegFunc(Double_t *x,Double_t *par) {
77 // This is the function used for the fit of the covariance 
78 // matrix as a function of the momentum. 
79 // For cosl the average value 0.87 is used.
80   Double_t value = par[0]+par[1]/TMath::Power(x[0],par[2])/
81                    TMath::Power(0.87,par[3]);
82   return value;
83 }
84
85 // structure for DB building
86 typedef struct {
87   Int_t    pdg;
88   Int_t    bin;
89   Double_t r;
90   Double_t p;
91   Double_t pt;
92   Double_t cosl;
93   Double_t eta;
94   Double_t dpt;
95   Double_t dP0,dP1,dP2,dP3,dP4;
96   Double_t c00,c10,c11,c20,c21,c22,c30,c31,c32,c33,c40,c41,c42,c43,c44;
97   Double_t dEdx;
98 } COMPTRACK;
99 // cov matrix structure 
100 typedef struct {
101   Double_t c00,c10,c11,c20,c21,c22,c30,c31,c32,c33,c40,c41,c42,c43,c44;
102 } COVMATRIX;
103
104 ClassImp(AliTPCtrackerParam)
105
106 //-----------------------------------------------------------------------------
107 AliTPCtrackerParam::AliTPCtrackerParam(const Int_t coll,const Double_t Bz)
108 {
109 //-----------------------------------------------------------------------------
110 // This is the class conctructor 
111 //-----------------------------------------------------------------------------
112
113   fBz = Bz;             // value of the z component of L3 field (Tesla)
114   fColl = coll;         // collision code (0: PbPb6000)
115   fSelAndSmear = kTRUE; // by default selection and smearing are done
116
117   if(fBz!=0.4) {
118     cerr<<"AliTPCtrackerParam::AliTPCtrackerParam:  Invalid field!\n";
119     cerr<<"      Available:  0.4"<<endl;
120   }
121   if(fColl!=0) { 
122     cerr<<"AliTPCtrackerParam::AliTPCtrackerParam:  Invalid collision!\n";
123     cerr<<"      Available:  0   ->   PbPb6000"<<endl; 
124   }
125
126   fDBfileName = gSystem->Getenv("ALICE_ROOT");  
127   fDBfileName.Append("/TPC/CovMatrixDB_");
128   if(fColl==0) fDBfileName.Append("PbPb6000");
129   if(fBz==0.4) fDBfileName.Append("_B0.4T.root");
130 }
131 //-----------------------------------------------------------------------------
132 AliTPCtrackerParam::~AliTPCtrackerParam() 
133 {}
134 //-----------------------------------------------------------------------------
135
136 Int_t AliTPCtrackerParam::BuildTPCtracks(const TFile *inp, TFile *out, Int_t n)
137 {
138 //-----------------------------------------------------------------------------
139 // This function creates the TPC parameterized tracks
140 //-----------------------------------------------------------------------------
141
142   TFile *fileDB=0;
143   TTree *covTreePi[50];
144   TTree *covTreeKa[50];
145   TTree *covTreeEl[50];
146
147   if(fSelAndSmear) {
148     cerr<<"+++\n+++ Reading DataBase from:\n+++ "<<
149       fDBfileName.Data()<<"\n+++\n"; 
150     // Read paramters from DB file
151     if(!ReadAllData(fDBfileName.Data())) {
152       cerr<<"AliTPCtrackerParam::BuildTPCtracks: 
153              Could not read data from DB\n\n"; return 1; 
154     }
155     // Read the trees with regularized cov. matrices from DB
156     TString str;
157     fileDB = TFile::Open(fDBfileName.Data());
158     Int_t nBinsPi = fDBgridPi.GetTotBins();
159     for(Int_t l=0; l<nBinsPi; l++) {
160       str = "/CovMatrices/Pions/CovTreePi_bin";
161       str+=l;
162       covTreePi[l] = (TTree*)fileDB->Get(str.Data());
163     }
164     Int_t nBinsKa = fDBgridKa.GetTotBins();
165     for(Int_t l=0; l<nBinsKa; l++) {
166       str = "/CovMatrices/Kaons/CovTreeKa_bin";
167       str+=l;
168       covTreeKa[l] = (TTree*)fileDB->Get(str.Data());
169     }
170     Int_t nBinsEl = fDBgridEl.GetTotBins();
171     for(Int_t l=0; l<nBinsEl; l++) {
172       str = "/CovMatrices/Electrons/CovTreeEl_bin";
173       str+=l;
174       covTreeEl[l] = (TTree*)fileDB->Get(str.Data());
175     }
176
177   } else cerr<<"\n ! Creating ALL TRUE tracks at TPC 1st hit !\n\n";
178
179   TFile *infile=(TFile*)inp;
180   infile->cd();
181
182   // Get gAlice object from file
183   if(!(gAlice=(AliRun*)infile->Get("gAlice"))) {
184     cerr<<"gAlice has not been found on galice.root !\n";
185     return 1;
186   }
187
188   // Check if value in the galice file is equal to selected one (fBz)
189   AliMagFCM *fiel = (AliMagFCM*)gAlice->Field();
190   Double_t fieval=(Double_t)fiel->SolenoidField()/10.;
191   printf("Magnetic field is %6.2f Tesla\n",fieval);
192   if(fBz!=fieval) {
193     cerr<<"AliTPCtrackerParam::BuildTPCtracks:  Invalid field!"<<endl;
194     cerr<<"Field selected is: "<<fBz<<" T\n";
195     cerr<<"Field found on file is: "<<fieval<<" T\n";
196     return 0;
197   }
198
199   // Get TPC detector 
200   AliTPC *TPC=(AliTPC*)gAlice->GetDetector("TPC");
201   Int_t ver = TPC->IsVersion(); 
202   cerr<<"+++ TPC version "<<ver<<" has been found !\n";
203   AliTPCParam *digp=(AliTPCParam*)infile->Get("75x40_100x60");
204   if(digp){
205     delete digp;
206     digp = new AliTPCParamSR();
207   }
208   else digp=(AliTPCParam*)infile->Get("75x40_100x60_150x60");
209   
210   if(!digp) { cerr<<"TPC parameters have not been found !\n"; return 1; }
211   TPC->SetParam(digp);
212
213   // Set the conversion constant between curvature and Pt
214   AliKalmanTrack::SetConvConst(100/0.299792458/fBz);
215
216   TParticle   *Part=0;
217   AliTPChit   *tpcHit=0;
218   AliTPCtrack *tpctrack=0;
219   Double_t     hPx,hPy,hPz,hPt,hEta,xg,yg,zg,xl,yl,zl;
220   Double_t     alpha;
221   Float_t      cosAlpha,sinAlpha;
222   Int_t        bin,label,pdg,charge;
223   Int_t        tracks=0;
224   Int_t        nParticles,nTracks,arrentr;
225   Char_t       tname[100];
226   //Int_t nSel=0,nAcc=0;
227
228   // loop over first n events in file
229   for(Int_t evt=0; evt<n; evt++){
230     cerr<<"+++\n+++ Processing event "<<evt<<"\n+++\n";
231
232     // tree for TPC tracks
233     sprintf(tname,"TreeT_TPC_%d",evt);
234     TTree *tracktree = new TTree(tname,"Tree with TPC tracks");
235     tracktree->Branch("tracks","AliTPCtrack",&tpctrack,20000,0);
236
237     // array for TPC tracks
238     TObjArray tarray(20000);
239
240     // get the particles stack
241     nParticles = gAlice->GetEvent(evt);
242     Bool_t done[nParticles];
243     Int_t  pdgCodes[nParticles];
244
245     // loop on particles and store pdg codes
246     for(Int_t l=0; l<nParticles; l++) {
247       Part        = gAlice->Particle(l);
248       pdgCodes[l] = Part->GetPdgCode();
249       done[l]     = kFALSE;
250     }
251
252     // Get TreeH with hits
253     TTree *TH=gAlice->TreeH(); 
254     nTracks=(Int_t)TH->GetEntries();
255     cerr<<"+++\n+++ Number of particles in event "<<evt<<":  "<<nParticles<<
256       "\n+++\n+++ Number of \"primary tracks\"(entries in TreeH): "<<nTracks<<
257       "\n+++\n\n";
258
259     // loop over entries in TreeH
260     for(Int_t l=0; l<nTracks; l++) {
261       if(l%1000==0) cerr<<"  --- Processing primary track "
262                         <<l<<" of "<<nTracks<<" ---\r";
263       TPC->ResetHits();
264       TH->GetEvent(l);
265       // Get FirstHit
266       tpcHit=(AliTPChit*)TPC->FirstHit(-1);
267       for( ; tpcHit; tpcHit=(AliTPChit*)TPC->NextHit() ) {
268         if(tpcHit->fQ !=0.) continue;
269         // Get particle momentum at hit
270         hPx=tpcHit->X(); hPy=tpcHit->Y(); hPz=tpcHit->Z();
271         hPt=TMath::Sqrt(hPx*hPx+hPy*hPy);
272         // reject hits with Pt<mag*0.45 GeV/c
273         if(hPt<(fBz*0.45)) continue;
274
275         // Get track label
276         label=tpcHit->Track();
277         // check if this track has already been processed
278         if(done[label]) continue;
279         // PDG code & electric charge
280         pdg = pdgCodes[label];
281         if(pdg>200 || pdg==-11 || pdg==-13) { charge=1; }
282         else if(pdg<-200 || pdg==11 || pdg==13) { charge=-1; }
283         else continue;
284         pdg = TMath::Abs(pdg);
285         if(pdg>3000) pdg=211;
286         if(fSelAndSmear) SetParticle(pdg);
287
288         if((tpcHit=(AliTPChit*)TPC->NextHit())==0) break;
289         if(tpcHit->fQ != 0.) continue;
290         // Get global coordinates of hit
291         xg=tpcHit->X(); yg=tpcHit->Y(); zg=tpcHit->Z();
292         if(TMath::Sqrt(xg*xg+yg*yg)>90.) continue;
293
294         // Get TPC sector, Alpha angle and local coordinates
295         // cerr<<"Sector "<<tpcHit->fSector<<endl;
296         digp->AdjustCosSin(tpcHit->fSector,cosAlpha,sinAlpha);
297         alpha = TMath::ATan2(sinAlpha,cosAlpha);
298         xl = xg*cosAlpha + yg*sinAlpha;
299         yl =-xg*sinAlpha + yg*cosAlpha;
300         zl = zg;
301         //printf("Alpha %f   xl %f  yl %f  zl %f\n",Alpha,xl,yl,zl);
302
303         // reject tracks which are not in the TPC acceptance
304         if(TMath::Abs(zl+(244.-xl)*hPz/hPt)>252.) continue;
305
306         hEta = -TMath::Log(TMath::Tan(0.25*TMath::Pi()-0.5*TMath::ATan(hPz/hPt)));  
307
308         // Apply selection according to TPC efficiency
309         //if(TMath::Abs(pdg)==211) nAcc++;
310         if(fSelAndSmear && !SelectedTrack(hPt,hEta)) continue; 
311         //if(TMath::Abs(pdg)==211) nSel++;
312
313         // create AliTPCtrack object
314         BuildTrack(alpha,xl,yl,zl,hPx,hPy,hPz,hPt,charge);
315         
316         if(fSelAndSmear) {
317           bin = fDBgrid->GetBin(hPt,hEta);
318           switch (pdg) {
319           case 211:
320             fCovTree = covTreePi[bin];
321             break;
322           case 321:
323             fCovTree = covTreeKa[bin];
324             break;
325           case 2212:
326             fCovTree = covTreePi[bin];
327             break;
328           case 11:
329             fCovTree = covTreeEl[bin];
330             break;
331           case 13:
332             fCovTree = covTreePi[bin];
333             break;
334           }
335           // deal with covariance matrix and smearing of parameters
336           CookTrack(hPt,hEta);
337
338           // assign the track a dE/dx and make a rough PID
339           CookdEdx(hPt,hEta);
340         }
341
342         // put track in array
343         AliTPCtrack *iotrack = new AliTPCtrack(fTrack);
344         iotrack->SetLabel(label);
345         tarray.AddLast(iotrack);
346         // Mark track as "done" and register the pdg code
347         done[label]=kTRUE; 
348         tracks++;
349       }
350  
351     } // loop over entries in TreeH
352
353     // sort array with TPC tracks (decreasing pT)
354     tarray.Sort();
355
356     arrentr = tarray.GetEntriesFast();
357     for(Int_t l=0; l<arrentr; l++) {
358       tpctrack=(AliTPCtrack*)tarray.UncheckedAt(l);
359       tracktree->Fill();
360     }
361
362     // write the tree with tracks in the output file
363     out->cd();
364     tracktree->Write();
365
366     delete tracktree;
367
368     printf("\n\n+++\n+++ Number of TPC tracks: %d\n+++\n",tracks);
369     //cerr<<"Average Eff: "<<(Float_t)nSel/nAcc<<endl;
370
371   } // loop on events
372
373   if(fileDB) fileDB->Close();
374
375   return 0;
376 }
377 //-----------------------------------------------------------------------------
378 void AliTPCtrackerParam::AnalyzedEdx(const Char_t *outName,Int_t pdg) {
379 //-----------------------------------------------------------------------------
380 // This function computes the dE/dx for pions, kaons, protons and electrons,
381 // in the [pT,eta] bins.
382 // Input file is CovMatrix_AllEvts.root.  
383 //-----------------------------------------------------------------------------
384
385   gStyle->SetOptStat(0);
386   gStyle->SetOptFit(10001);
387
388   Char_t *part="PIONS";
389   Double_t ymax=500.;
390
391   // create a chain with compared tracks
392   TChain cmptrkchain("cmptrktree");
393   cmptrkchain.Add("CovMatrix_AllEvts_1.root");
394   cmptrkchain.Add("CovMatrix_AllEvts_2.root");
395   cmptrkchain.Add("CovMatrix_AllEvts_3.root");
396
397   COMPTRACK cmptrk; 
398   cmptrkchain.SetBranchAddress("comptracks",&cmptrk);
399   Int_t entries = (Int_t)cmptrkchain.GetEntries(); 
400   cerr<<" Number of entries: "<<entries<<endl;
401
402   InitializeKineGrid("DB","points");
403   InitializeKineGrid("dEdx","bins");
404
405   switch(pdg) {
406   case 211:
407     part = "PIONS";
408     ymax = 100.;
409     break;
410   case 321:
411     part = "KAONS";
412     ymax = 300.;
413     break;
414   case 2212:
415     part = "PROTONS";
416     ymax = 500.;
417     break;
418   case 11:
419     part = "ELECTRONS";
420     ymax = 100.;
421     break;
422   }
423
424   SetParticle(pdg);
425
426   const Int_t nTotBins = fDBgrid->GetTotBins(); 
427
428   cerr<<" Fit bins: "<<nTotBins<<endl;
429
430   Int_t bin=0;
431   Int_t n[nTotBins];
432   Double_t p[nTotBins],ep[nTotBins],mean[nTotBins],sigma[nTotBins];
433
434   for(Int_t l=0; l<nTotBins; l++) {
435     n[l] = 1; // set to 1 to avoid divisions by 0
436     p[l] = mean[l] = sigma[l] = ep[l] = 0.; 
437   }
438
439   // loop on chain entries for the mean
440   for(Int_t l=0; l<entries; l++) {
441     cmptrkchain.GetEvent(l);
442     if(TMath::Abs(cmptrk.pdg)!=pdg) continue;
443     bin = fDBgrid->GetBin(cmptrk.pt,cmptrk.eta);
444     p[bin] += cmptrk.p;
445     mean[bin] += cmptrk.dEdx;
446     n[bin]++;
447   } // loop on chain entries
448
449   for(Int_t l=0; l<nTotBins; l++) {
450     p[l] /= n[l];
451     mean[l] /= n[l];
452     n[l] = 1; // set to 1 to avoid divisions by 0
453   }
454
455   // loop on chain entries for the sigma
456   for(Int_t l=0; l<entries; l++) {
457     cmptrkchain.GetEvent(l);
458     if(TMath::Abs(cmptrk.pdg)!=pdg) continue;
459     bin = fDBgrid->GetBin(cmptrk.pt,cmptrk.eta);
460     if(cmptrk.p<1. && TMath::Abs(cmptrk.p-p[bin])>0.025) continue;
461     n[bin]++;
462     sigma[bin] += (cmptrk.dEdx-mean[bin])*(cmptrk.dEdx-mean[bin]);
463   } // loop on chain entries
464   
465   for(Int_t l=0; l<nTotBins; l++) {
466     sigma[l] = TMath::Sqrt(sigma[l]/n[l]);
467   }
468
469   // create the canvas
470   TCanvas *canv = new TCanvas("canv","dEdx",0,0,900,700); 
471
472   // create the graph for dEdx vs p
473   TGraphErrors *gr = new TGraphErrors(nTotBins,p,mean,ep,sigma);
474   TString title("  : dE/dx vs momentum"); title.Prepend(part);
475   TH2F *frame = new TH2F("frame1",title.Data(),2,0.1,50,2,0,ymax);
476   frame->SetXTitle("p [GeV/c]");
477   frame->SetYTitle("dE/dx [a.u.]");
478   canv->SetLogx();
479   frame->Draw();
480   gr->Draw("P");
481
482   switch(pdg) {
483   case 211:
484     for(Int_t i=0; i<nTotBins; i++) {
485       fdEdxMeanPi.SetParam(i,mean[i]);
486       fdEdxRMSPi.SetParam(i,sigma[i]);
487     }    
488     break;
489   case 321:
490     for(Int_t i=0; i<nTotBins; i++) {
491       fdEdxMeanKa.SetParam(i,mean[i]);
492       fdEdxRMSKa.SetParam(i,sigma[i]);
493     }    
494     break;
495   case 2212:
496     for(Int_t i=0; i<nTotBins; i++) {
497       fdEdxMeanPr.SetParam(i,mean[i]);
498       fdEdxRMSPr.SetParam(i,sigma[i]);
499     }    
500     break;
501   case 11:
502     for(Int_t i=0; i<nTotBins; i++) {
503       fdEdxMeanEl.SetParam(i,mean[i]);
504       fdEdxRMSEl.SetParam(i,sigma[i]);
505     }    
506     break;
507   }
508
509   // write results to file
510   WritedEdx(outName,pdg);
511
512   
513   return;
514 }
515 //-----------------------------------------------------------------------------
516 void AliTPCtrackerParam::AnalyzePulls(const Char_t *outName) {
517 //-----------------------------------------------------------------------------
518 // This function computes the pulls for pions, kaons and electrons,
519 // in the [pT,eta] bins.
520 // Input file is CovMatrix_AllEvts.root.
521 // Output file is pulls.root.  
522 //-----------------------------------------------------------------------------
523
524   // create a chain with compared tracks
525   TChain cmptrkchain("cmptrktree");
526   cmptrkchain.Add("CovMatrix_AllEvts_1.root");
527   cmptrkchain.Add("CovMatrix_AllEvts_2.root");
528   cmptrkchain.Add("CovMatrix_AllEvts_3.root");
529
530   COMPTRACK cmptrk; 
531   cmptrkchain.SetBranchAddress("comptracks",&cmptrk);
532   Int_t entries = (Int_t)cmptrkchain.GetEntries(); 
533   cerr<<" Number of entries: "<<entries<<endl;
534
535   Int_t thisPdg=0;
536   Char_t hname[100], htitle[100];
537   Int_t nTotBins,bin;
538
539   AliTPCkineGrid  pulls[5];
540   TH1F *hDum = new TH1F("name","title",100,-7.,7.);
541   TF1 *g = new TF1("g","gaus");
542
543   InitializeKineGrid("pulls","bins");
544   InitializeKineGrid("DB","points");
545
546
547
548   // loop on the particles Pi,Ka,El
549   for(Int_t part=0; part<3; part++) {
550
551     switch (part) {
552     case 0: // pions
553       thisPdg=211; 
554       cerr<<" Processing pions ...\n";
555       break;   
556     case 1: // kaons
557       thisPdg=321; 
558       cerr<<" Processing kaons ...\n";
559       break;
560       
561     case 2: // electrons
562       thisPdg=11; 
563       cerr<<" Processing electrons ...\n";
564       break;
565     }
566
567     SetParticle(thisPdg);
568
569     for(Int_t i=0;i<5;i++) {
570       pulls[i].~AliTPCkineGrid(); 
571       new(&pulls[i]) AliTPCkineGrid(*(fPulls+i));
572     }
573     nTotBins = fDBgrid->GetTotBins();
574  
575     // create histograms for the all the bins
576     TH1F *hPulls0_=NULL;
577     TH1F *hPulls1_=NULL;
578     TH1F *hPulls2_=NULL;
579     TH1F *hPulls3_=NULL;
580     TH1F *hPulls4_=NULL;
581
582     hPulls0_ = new TH1F[nTotBins]; 
583     hPulls1_ = new TH1F[nTotBins]; 
584     hPulls2_ = new TH1F[nTotBins]; 
585     hPulls3_ = new TH1F[nTotBins]; 
586     hPulls4_ = new TH1F[nTotBins]; 
587
588
589     for(Int_t i=0; i<nTotBins; i++) {
590       sprintf(hname,"hPulls0_%d",i);
591       sprintf(htitle,"P0 pulls for bin %d",i);
592       hDum->SetName(hname); hDum->SetTitle(htitle);
593       hPulls0_[i] = *hDum;
594       sprintf(hname,"hPulls1_%d",i);
595       sprintf(htitle,"P1 pulls for bin %d",i);
596       hDum->SetName(hname); hDum->SetTitle(htitle);
597       hPulls1_[i] = *hDum;
598       sprintf(hname,"hPulls2_%d",i);
599       sprintf(htitle,"P2 pulls for bin %d",i);
600       hDum->SetName(hname); hDum->SetTitle(htitle);
601       hPulls2_[i] = *hDum;
602       sprintf(hname,"hPulls3_%d",i);
603       sprintf(htitle,"P3 pulls for bin %d",i);
604       hDum->SetName(hname); hDum->SetTitle(htitle);
605       hPulls3_[i] = *hDum;
606       sprintf(hname,"hPulls4_%d",i);
607       sprintf(htitle,"P4 pulls for bin %d",i);
608       hDum->SetName(hname); hDum->SetTitle(htitle);
609       hPulls4_[i] = *hDum;
610     }
611
612     // loop on chain entries 
613     for(Int_t i=0; i<entries; i++) {
614       cmptrkchain.GetEvent(i);
615       if(TMath::Abs(cmptrk.pdg)!=thisPdg) continue;
616       // fill histograms with the pulls
617       bin = fDBgrid->GetBin(cmptrk.pt,cmptrk.eta); 
618       hPulls0_[bin].Fill(cmptrk.dP0/TMath::Sqrt(cmptrk.c00));
619       hPulls1_[bin].Fill(cmptrk.dP1/TMath::Sqrt(cmptrk.c11));
620       hPulls2_[bin].Fill(cmptrk.dP2/TMath::Sqrt(cmptrk.c22));
621       hPulls3_[bin].Fill(cmptrk.dP3/TMath::Sqrt(cmptrk.c33));
622       hPulls4_[bin].Fill(cmptrk.dP4/TMath::Sqrt(cmptrk.c44));
623     } // loop on chain entries
624
625     // compute the sigma of the distributions
626     for(Int_t i=0; i<nTotBins; i++) {
627       if(hPulls0_[i].GetEntries()>1000) {
628         g->SetRange(-3.*hPulls0_[i].GetRMS(),3.*hPulls0_[i].GetRMS());
629         hPulls0_[i].Fit("g","R,Q,N");
630         pulls[0].SetParam(i,g->GetParameter(2));
631       } else pulls[0].SetParam(i,-1.);
632       if(hPulls1_[i].GetEntries()>1000) {
633         g->SetRange(-3.*hPulls1_[i].GetRMS(),3.*hPulls1_[i].GetRMS());
634         hPulls1_[i].Fit("g","R,Q,N");
635         pulls[1].SetParam(i,g->GetParameter(2));
636       } else pulls[1].SetParam(i,-1.);
637       if(hPulls2_[i].GetEntries()>1000) {
638         g->SetRange(-3.*hPulls2_[i].GetRMS(),3.*hPulls2_[i].GetRMS());
639         hPulls2_[i].Fit("g","R,Q,N");
640         pulls[2].SetParam(i,g->GetParameter(2));
641       } else pulls[2].SetParam(i,-1.);
642       if(hPulls3_[i].GetEntries()>1000) {
643         g->SetRange(-3.*hPulls3_[i].GetRMS(),3.*hPulls3_[i].GetRMS());
644         hPulls3_[i].Fit("g","R,Q,N");
645         pulls[3].SetParam(i,g->GetParameter(2));
646       } else pulls[3].SetParam(i,-1.);
647       if(hPulls4_[i].GetEntries()>1000) {
648         g->SetRange(-3.*hPulls4_[i].GetRMS(),3.*hPulls4_[i].GetRMS());
649         hPulls4_[i].Fit("g","R,Q,N");
650         pulls[4].SetParam(i,g->GetParameter(2));
651       } else pulls[4].SetParam(i,-1.);
652     } // loop on bins
653
654
655     switch (part) {
656     case 0: // pions
657       for(Int_t i=0;i<5;i++) {
658         fPullsPi[i].~AliTPCkineGrid();
659         new(&fPullsPi[i]) AliTPCkineGrid(pulls[i]);
660       }
661       break;
662     case 1: // kaons
663       for(Int_t i=0;i<5;i++) {
664         fPullsKa[i].~AliTPCkineGrid();
665         new(&fPullsKa[i]) AliTPCkineGrid(pulls[i]);
666       }
667       break;
668     case 2: // electrons
669       for(Int_t i=0;i<5;i++) {
670         fPullsEl[i].~AliTPCkineGrid();
671         new(&fPullsEl[i]) AliTPCkineGrid(pulls[i]);
672       }
673       break;
674     }
675
676     delete [] hPulls0_;
677     delete [] hPulls1_;
678     delete [] hPulls2_;
679     delete [] hPulls3_;
680     delete [] hPulls4_;
681     
682   } // loop on particle species
683
684   // write pulls to file
685   WritePulls(outName);
686
687
688   return;
689 }
690 //-----------------------------------------------------------------------------
691 void AliTPCtrackerParam::BuildTrack(Double_t alpha,
692                                     Double_t x,Double_t y,Double_t z,
693                                     Double_t px,Double_t py,
694                                     Double_t pz,Double_t pt,
695                                     Int_t ch) {  
696 //-----------------------------------------------------------------------------
697 // This function uses GEANT info to set true track parameters
698 //-----------------------------------------------------------------------------
699   Double_t xref = x;
700   Double_t xx[5],cc[15];
701   cc[0]=cc[2]=cc[5]=cc[9]=cc[14]=10.;
702   cc[1]=cc[3]=cc[4]=cc[6]=cc[7]=cc[8]=cc[10]=cc[11]=cc[12]=cc[13]=0.;
703   
704   // Magnetic field
705   TVector3 bfield(0.,0.,fBz);
706   
707   
708   // radius [cm] of track projection in (x,y) 
709   Double_t rho = pt*100./0.299792458/bfield.Z();
710   // center of track projection in local reference frame
711   TVector3 hmom,hpos;
712
713
714   // position (local) and momentum (local) at the hit
715   // in the bending plane (z=0)
716   hpos.SetXYZ(x,y,0.);
717   hmom.SetXYZ(px*TMath::Cos(alpha)+py*TMath::Sin(alpha),-px*TMath::Sin(alpha)+py*TMath::Cos(alpha),0.);
718   TVector3 vrho = hmom.Cross(bfield);
719   vrho *= ch;
720   vrho.SetMag(rho);
721
722   TVector3 vcenter = hpos+vrho;
723
724   Double_t x0 = vcenter.X();
725
726   // fX     = xref         X-coordinate of this track (reference plane)
727   // fAlpha = Alpha        Rotation angle the local (TPC sector) 
728   // fP0    = YL           Y-coordinate of a track
729   // fP1    = ZG           Z-coordinate of a track
730   // fP2    = C*x0         x0 is center x in rotated frame
731   // fP3    = Tgl          tangent of the track momentum dip angle
732   // fP4    = C            track curvature
733   xx[0] = y;
734   xx[1] = z;
735   xx[3] = pz/pt;
736   xx[4] = -ch/rho;
737   xx[2] = xx[4]*x0;
738
739   // create the object AliTPCtrack    
740   AliTPCtrack track(0,xx,cc,xref,alpha);
741   new(&fTrack) AliTPCtrack(track);
742
743   return;
744 }
745 //-----------------------------------------------------------------------------
746 void AliTPCtrackerParam::CompareTPCtracks(
747                            const Char_t* galiceName="galice.root",
748                            const Char_t* trkGeaName="AliTPCtracksGeant.root",
749                            const Char_t* trkKalName="AliTPCtracksSorted.root",
750                            const Char_t* covmatName="CovMatrix.root",
751                            const Char_t* tpceffName="TPCeff.dat") const {
752 //-----------------------------------------------------------------------------
753 // This function compares tracks from TPC Kalman Filter V2 with 
754 // true tracks at TPC 1st hit. It gives:
755 //   - a tree with Kalman cov. matrix elements, resolutions, dEdx
756 //   - the efficiencies as a function of particle type, pT, eta
757 //-----------------------------------------------------------------------------
758
759   TFile *kalFile    = TFile::Open(trkKalName);
760   TFile *geaFile    = TFile::Open(trkGeaName);
761   TFile *galiceFile = TFile::Open(galiceName);
762
763   // particles from TreeK
764   AliRun *gAlice = (AliRun*)galiceFile->Get("gAlice");
765   const Int_t nparticles = gAlice->GetEvent(0);
766
767   Int_t kalLab[nparticles];
768   for(Int_t i=0; i<nparticles; i++) kalLab[i] = -1; 
769  
770
771   // tracks from Kalman
772   TTree *kaltree=(TTree*)kalFile->Get("TreeT_TPC_0");
773   AliTPCtrack *kaltrack=new AliTPCtrack; 
774   kaltree->SetBranchAddress("tracks",&kaltrack);
775   Int_t kalEntries = (Int_t)kaltree->GetEntries();
776
777   // tracks from 1st hit
778   TTree *geatree=(TTree*)geaFile->Get("TreeT_TPC_0");
779   AliTPCtrack *geatrack=new AliTPCtrack; 
780   geatree->SetBranchAddress("tracks",&geatrack);
781   Int_t geaEntries = (Int_t)geatree->GetEntries();
782
783   cerr<<"+++\n+++ Number of tracks:  TPC Kalman  = "<<kalEntries<<endl<<"+++                    TPC 1st hit = "<<geaEntries<<endl<<"+++\n";
784
785   // set pointers for TPC tracks: 
786   // the entry number of the track labelled l is stored in kalLab[l]
787   Int_t fake=0,mult=0;
788   for (Int_t j=0; j<kalEntries; j++) {
789     kaltree->GetEvent(j);
790     if(kaltrack->GetLabel()>=0) { 
791       if(kalLab[kaltrack->GetLabel()]!=-1) mult++; 
792       kalLab[kaltrack->GetLabel()] = j;
793     }
794     else fake++;
795   }
796   
797   cerr<<"+++ Number of fake tracks in TPC Kalman: "<<fake<<endl;
798   cerr<<"+++ Number of multiply found tracks in TPC Kalman: "<<mult<<endl;
799
800   TParticle *Part;
801   Int_t    label;
802   Double_t cc[15],dAlpha;
803   //                 ---  [Pt,eta] binning ---
804   //
805   //               |eta|<0.3    0.3<=|eta|<0.6    0.6<=|eta|
806   //      Pt<0.3       0              1                2 
807   // 0.3<=Pt<0.5       3              4                5  
808   // 0.5<=Pt<1.0       6              7                8 
809   // 1.0<=Pt<1.5       9             10               11 
810   // 1.5<=Pt<3.0      12             13               14 
811   // 3.0<=Pt<5.0      15             16               17 
812   // 5.0<=Pt<7.0      18             19               20 
813   // 7.0<=Pt<15.      21             22               23 
814   // 15.<=Pt          24             25               26
815   // 
816   Int_t    pi=0,ka=0,mu=0,el=0,pr=0;
817   Int_t    geaPi[27],geaKa[27],geaPr[27],geaEl[27],geaMu[27];
818   Int_t    kalPi[27],kalKa[27],kalPr[27],kalEl[27],kalMu[27];
819   Float_t  effPi[27],effKa[27],effPr[27],effEl[27],effMu[27];
820   Int_t bin;
821   for(Int_t j=0; j<27; j++) {
822     geaPi[j]=geaKa[j]=geaPr[j]=geaEl[j]=geaMu[j]=0;
823     kalPi[j]=kalKa[j]=kalPr[j]=kalEl[j]=kalMu[j]=0;
824     effPi[j]=effKa[j]=effPr[j]=effEl[j]=effMu[j]=-1.;
825   }
826
827   COMPTRACK cmptrk;
828   // create the tree for comparison results
829   TTree *cmptrktree = new TTree("cmptrktree","results of track comparison");
830   cmptrktree->Branch("comptracks",&cmptrk,"pdg/I:bin:r/D:p:pt:cosl:eta:dpt:dP0:dP1:dP2:dP3:dP4:c00:c10:c11:c20:c21:c22:c30:c31:c32:c33:c40:c41:c42:c43:c44:dEdx");
831   
832   cerr<<"Doing track comparison...\n";
833   // loop on tracks at 1st hit
834   for(Int_t j=0; j<geaEntries; j++) {
835     geatree->GetEvent(j);
836     
837     label = geatrack->GetLabel();
838     Part = (TParticle*)gAlice->Particle(label);
839     cmptrk.pdg = Part->GetPdgCode();
840     cmptrk.eta = Part->Eta();
841     cmptrk.r = TMath::Sqrt(Part->Vx()*Part->Vx()+Part->Vy()*Part->Vy());
842
843     cmptrk.pt   = 1/TMath::Abs(geatrack->Get1Pt());
844     cmptrk.cosl = TMath::Cos(TMath::ATan(geatrack->GetTgl()));
845     cmptrk.p    = cmptrk.pt/cmptrk.cosl;
846
847
848     bin = GetBin(cmptrk.pt,cmptrk.eta);
849     cmptrk.bin = bin;
850     if(abs(cmptrk.pdg)==211)  geaPi[bin]++;  
851     if(abs(cmptrk.pdg)==321)  geaKa[bin]++;   
852     if(abs(cmptrk.pdg)==2212) geaPr[bin]++;   
853     if(abs(cmptrk.pdg)==11)   geaEl[bin]++;  
854     if(abs(cmptrk.pdg)==13)   geaMu[bin]++; 
855
856
857     // check if there is the corresponding track in the TPC kalman and get it
858     if(kalLab[label]==-1) continue;
859     kaltree->GetEvent(kalLab[label]);
860
861     // and go on only if it has xref = 84.57 cm (inner pad row)
862     if(kaltrack->GetX()>90.) continue;
863
864     if(abs(cmptrk.pdg)==211)  { kalPi[bin]++; pi++; }
865     if(abs(cmptrk.pdg)==321)  { kalKa[bin]++; ka++; }
866     if(abs(cmptrk.pdg)==2212) { kalPr[bin]++; pr++; }
867     if(abs(cmptrk.pdg)==11)   { kalEl[bin]++; el++; }
868     if(abs(cmptrk.pdg)==13)   { kalMu[bin]++; mu++; }
869
870     kaltrack->PropagateTo(geatrack->GetX());
871
872     cmptrk.dEdx = kaltrack->GetdEdx();
873
874     // compute errors on parameters
875     dAlpha = kaltrack->GetAlpha()-geatrack->GetAlpha();
876     if(TMath::Abs(dAlpha)>0.1) { cerr<<" ! WRONG SECTOR !\n"; continue; }
877
878     cmptrk.dP0 = kaltrack->GetY()-geatrack->GetY();
879     cmptrk.dP1 = kaltrack->GetZ()-geatrack->GetZ();
880     cmptrk.dP2 = kaltrack->GetEta()-geatrack->GetEta();
881     cmptrk.dP3 = kaltrack->GetTgl()-geatrack->GetTgl();
882     cmptrk.dP4 = kaltrack->GetC()-geatrack->GetC();
883     cmptrk.dpt = 1/kaltrack->Get1Pt()-1/geatrack->Get1Pt();
884
885     // get covariance matrix
886     // beware: lines 3 and 4 are inverted!
887     kaltrack->GetCovariance(cc);
888
889     cmptrk.c00 = cc[0];
890     cmptrk.c10 = cc[1];
891     cmptrk.c11 = cc[2];
892     cmptrk.c20 = cc[3];
893     cmptrk.c21 = cc[4];
894     cmptrk.c22 = cc[5];
895     cmptrk.c30 = cc[10];
896     cmptrk.c31 = cc[11];
897     cmptrk.c32 = cc[12];
898     cmptrk.c33 = cc[14];
899     cmptrk.c40 = cc[6];
900     cmptrk.c41 = cc[7];
901     cmptrk.c42 = cc[8];
902     cmptrk.c43 = cc[13];
903     cmptrk.c44 = cc[9];
904
905     // fill tree
906     cmptrktree->Fill();
907
908   } // loop on tracks at TPC 1st hit
909
910   cerr<<"+++\n+++ Number of compared tracks: "<<pi+ka+el+mu+pr<<endl;
911   cerr<<"+++ Pions: "<<pi<<", Kaons: "<<ka<<", Protons : "<<pr<<", Electrons: "<<el<<", Muons: "<<mu<<endl;
912   cerr<<"+++"<<endl;
913
914   // Write tree to file
915   TFile *outfile = new TFile(covmatName,"recreate");
916   cmptrktree->Write();
917   outfile->Close();
918
919
920   // Write efficiencies to file
921   FILE *effFile = fopen(tpceffName,"w");
922   //fprintf(effFile,"%d\n",kalEntries);
923   for(Int_t j=0; j<27; j++) {
924     if(geaPi[j]>=5) effPi[j]=(Float_t)kalPi[j]/geaPi[j];
925     if(geaKa[j]>=5) effKa[j]=(Float_t)kalKa[j]/geaKa[j];
926     if(geaPr[j]>=5) effPr[j]=(Float_t)kalPr[j]/geaPr[j];
927     if(geaEl[j]>=5) effEl[j]=(Float_t)kalEl[j]/geaEl[j];
928     if(geaMu[j]>=5) effMu[j]=(Float_t)kalMu[j]/geaMu[j];
929     fprintf(effFile,"%f  %f  %f  %f  %f\n",effPi[j],effKa[j],effPr[j],effEl[j],effMu[j]);
930   }
931
932   for(Int_t j=0; j<27; j++) {
933     fprintf(effFile,"%d  %d  %d  %d  %d\n",geaPi[j],geaKa[j],geaPr[j],geaEl[j],geaMu[j]);
934   }
935   for(Int_t j=0; j<27; j++) {
936     fprintf(effFile,"%d  %d  %d  %d  %d\n",kalPi[j],kalKa[j],kalPr[j],kalEl[j],kalMu[j]);
937   }
938   fclose(effFile);
939
940   // delete AliRun object
941   delete gAlice; gAlice=0;
942   
943   // close all input files
944   kalFile->Close();
945   geaFile->Close();
946   galiceFile->Close();
947
948   return;
949 }
950 //-----------------------------------------------------------------------------
951 void AliTPCtrackerParam::CookdEdx(Double_t pt,Double_t eta) {
952 //-----------------------------------------------------------------------------
953 // This function assigns the track a dE/dx and makes a rough PID
954 //-----------------------------------------------------------------------------
955
956   Double_t mean = fdEdxMean->GetValueAt(pt,eta);
957   Double_t rms  = fdEdxRMS->GetValueAt(pt,eta);
958
959   Double_t dEdx = gRandom->Gaus(mean,rms);
960
961   fTrack.SetdEdx(dEdx);
962
963   AliTPCtrackParam t(fTrack);
964
965   //Very rough PID
966   Double_t p = TMath::Sqrt(1.+t.GetTgl()*t.GetTgl())*pt;
967
968   if (p<0.6) {
969     if (dEdx < 39.+ 12./(p+0.25)/(p+0.25)) { 
970       t.AssignMass(0.13957); new(&fTrack) AliTPCtrack(t); return;
971     }
972     if (dEdx < 39.+ 12./p/p) { 
973       t.AssignMass(0.49368); new(&fTrack) AliTPCtrack(t); return;
974     }
975     t.AssignMass(0.93827); new(&fTrack) AliTPCtrack(t); return;
976   }
977
978   if (p<1.2) {
979     if (dEdx < 39.+ 12./(p+0.25)/(p+0.25)) { 
980       t.AssignMass(0.13957); new(&fTrack) AliTPCtrack(t); return;
981     }
982     t.AssignMass(0.93827); new(&fTrack) AliTPCtrack(t); return;
983   }
984
985   t.AssignMass(0.13957); new(&fTrack) AliTPCtrack(t); return;
986 }
987 //-----------------------------------------------------------------------------
988 void AliTPCtrackerParam::CookTrack(Double_t pt,Double_t eta) {
989 //-----------------------------------------------------------------------------
990 // This function deals with covariance matrix and smearing
991 //-----------------------------------------------------------------------------
992
993   COVMATRIX covmat;
994   Double_t  p,cosl;
995   Double_t  trkKine[2],trkRegPar[4];    
996   Double_t  xref,alpha,xx[5],xxsm[5],cc[15];
997   Int_t     treeEntries;
998
999   fCovTree->SetBranchAddress("matrix",&covmat);
1000
1001   // get random entry from the tree
1002   treeEntries = (Int_t)fCovTree->GetEntries();
1003   fCovTree->GetEvent(gRandom->Integer(treeEntries));
1004
1005   // get P and Cosl from track
1006   cosl = TMath::Cos(TMath::ATan(fTrack.GetTgl()));
1007   p    = 1./TMath::Abs(fTrack.Get1Pt())/cosl;
1008
1009   trkKine[0] = p;
1010   trkKine[1] = cosl; 
1011
1012   // get covariance matrix from regularized matrix    
1013   for(Int_t l=0;l<4;l++) trkRegPar[l] = (*fRegPar)(0,l);
1014   cc[0] = covmat.c00*RegFunc(trkKine,trkRegPar);
1015   cc[1] = covmat.c10;
1016   for(Int_t l=0;l<4;l++) trkRegPar[l] = (*fRegPar)(1,l);
1017   cc[2] = covmat.c11*RegFunc(trkKine,trkRegPar);
1018   for(Int_t l=0;l<4;l++) trkRegPar[l] = (*fRegPar)(2,l);
1019   cc[3] = covmat.c20*RegFunc(trkKine,trkRegPar);
1020   cc[4] = covmat.c21;
1021   for(Int_t l=0;l<4;l++) trkRegPar[l] = (*fRegPar)(3,l);
1022   cc[5] = covmat.c22*RegFunc(trkKine,trkRegPar);
1023   cc[6] = covmat.c30;
1024   for(Int_t l=0;l<4;l++) trkRegPar[l] = (*fRegPar)(4,l);
1025   cc[7] = covmat.c31*RegFunc(trkKine,trkRegPar);
1026   cc[8] = covmat.c32;
1027   for(Int_t l=0;l<4;l++) trkRegPar[l] = (*fRegPar)(5,l);
1028   cc[9] = covmat.c33*RegFunc(trkKine,trkRegPar);
1029   for(Int_t l=0;l<4;l++) trkRegPar[l] = (*fRegPar)(6,l);
1030   cc[10]= covmat.c40*RegFunc(trkKine,trkRegPar);
1031   cc[11]= covmat.c41;
1032   for(Int_t l=0;l<4;l++) trkRegPar[l] = (*fRegPar)(7,l);
1033   cc[12]= covmat.c42*RegFunc(trkKine,trkRegPar);
1034   cc[13]= covmat.c43;
1035   for(Int_t l=0;l<4;l++) trkRegPar[l] = (*fRegPar)(8,l);
1036   cc[14]= covmat.c44*RegFunc(trkKine,trkRegPar);
1037    
1038   TMatrixD covMatSmear(5,5);
1039     
1040   covMatSmear = GetSmearingMatrix(cc,pt,eta);
1041
1042   // get track original parameters
1043   xref =fTrack.GetX();
1044   alpha=fTrack.GetAlpha();
1045   xx[0]=fTrack.GetY();
1046   xx[1]=fTrack.GetZ();
1047   xx[2]=fTrack.GetX()*fTrack.GetC()-fTrack.GetSnp();
1048   xx[3]=fTrack.GetTgl();
1049   xx[4]=fTrack.GetC();
1050     
1051   // use smearing matrix to smear the original parameters
1052   SmearTrack(xx,xxsm,covMatSmear);
1053     
1054   AliTPCtrack track(0,xxsm,cc,xref,alpha);
1055   new(&fTrack) AliTPCtrack(track);
1056   
1057   return; 
1058 }
1059 //-----------------------------------------------------------------------------
1060 void AliTPCtrackerParam::DrawEffs(const Char_t* inName,Int_t pdg=211) {
1061 //-----------------------------------------------------------------------------
1062 // This function draws the TPC efficiencies in the [pT,eta] bins
1063 //-----------------------------------------------------------------------------
1064
1065   ReadEffs(inName);
1066   SetParticle(pdg);
1067
1068   const Int_t n = fEff->GetPointsPt();
1069   Double_t effsA[n],effsB[n],effsC[n],pt[n];
1070   fEff->GetArrayPt(pt);
1071   for(Int_t i=0;i<n;i++) {
1072     effsA[i] = fEff->GetParam(i,0);
1073     effsB[i] = fEff->GetParam(i,1);
1074     effsC[i] = fEff->GetParam(i,2);
1075   }
1076   
1077   TGraph *grA = new TGraph(n,pt,effsA);
1078   TGraph *grB = new TGraph(n,pt,effsB);
1079   TGraph *grC = new TGraph(n,pt,effsC);
1080
1081   grA->SetMarkerStyle(20); grA->SetMarkerColor(2); grA->SetMarkerSize(1.5);
1082   grB->SetMarkerStyle(21); grB->SetMarkerColor(3); grB->SetMarkerSize(1.5);
1083   grC->SetMarkerStyle(22); grC->SetMarkerColor(4); grC->SetMarkerSize(1.5);
1084
1085   TString title("Distribution of the TPC efficiencies");
1086   switch (pdg) {
1087   case 211:
1088     title.Prepend("PIONS - ");
1089     break;
1090   case 321:
1091     title.Prepend("KAONS - ");
1092     break;
1093   case 2212:
1094     title.Prepend("PROTONS - ");
1095     break;
1096   case 11:
1097     title.Prepend("ELECTRONS - ");
1098     break;
1099   case 13:
1100     title.Prepend("MUONS - ");
1101     break;
1102   }
1103
1104
1105   gStyle->SetOptStat(0);
1106   TCanvas *c = new TCanvas("c","effs",0,0,900,700);
1107   c->SetLogx();
1108   c->SetGridx();
1109   c->SetGridy();
1110
1111   TH2F *frame = new TH2F("frame",title.Data(),2,0.1,30,2,0,1);
1112   frame->SetXTitle("p_{T} [GeV/c]");
1113   frame->Draw();
1114   grA->Draw("P");
1115   grB->Draw("P");
1116   grC->Draw("P");
1117
1118   TLegend *leg = new TLegend(0.2,0.2,0.4,0.4);
1119   leg->AddEntry(grA,"|#eta|<0.3","p");
1120   leg->AddEntry(grB,"0.3<|#eta|<0.6","p");
1121   leg->AddEntry(grC,"0.6<|#eta|<0.9","p");
1122   leg->SetFillColor(0);
1123   leg->Draw();
1124
1125
1126   return;
1127 }
1128 //-----------------------------------------------------------------------------
1129 void AliTPCtrackerParam::DrawPulls(const Char_t* inName,Int_t pdg=211,
1130                                    Int_t par=0) {
1131 //-----------------------------------------------------------------------------
1132 // This function draws the pulls in the [pT,eta] bins
1133 //-----------------------------------------------------------------------------
1134
1135   ReadPulls(inName);
1136   SetParticle(pdg);
1137
1138   const Int_t n = (fPulls+par)->GetPointsPt();
1139   Double_t pullsA[n],pullsB[n],pullsC[n],pt[n];
1140   (fPulls+par)->GetArrayPt(pt);  
1141   for(Int_t i=0;i<n;i++) {
1142     pullsA[i] = (fPulls+par)->GetParam(i,0);
1143     pullsB[i] = (fPulls+par)->GetParam(i,1);
1144     pullsC[i] = (fPulls+par)->GetParam(i,2);
1145   }
1146
1147   TGraph *grA = new TGraph(n,pt,pullsA);
1148   TGraph *grB = new TGraph(n,pt,pullsB);
1149   TGraph *grC = new TGraph(n,pt,pullsC);
1150
1151   grA->SetMarkerStyle(20); grA->SetMarkerColor(2); grA->SetMarkerSize(1.5);
1152   grB->SetMarkerStyle(21); grB->SetMarkerColor(3); grB->SetMarkerSize(1.5);
1153   grC->SetMarkerStyle(22); grC->SetMarkerColor(4); grC->SetMarkerSize(1.5);
1154
1155   TString title("Distribution of the pulls: ");
1156   switch (pdg) {
1157   case 211:
1158     title.Prepend("PIONS - ");
1159     break;
1160   case 321:
1161     title.Prepend("KAONS - ");
1162     break;
1163   case 11:
1164     title.Prepend("ELECTRONS - ");
1165     break;
1166   }
1167   switch (par) {
1168   case 0:
1169     title.Append("y");
1170     break;
1171   case 1:
1172     title.Append("z");
1173     break;
1174   case 2:
1175     title.Append("    #eta");
1176     break;
1177   case 3:
1178     title.Append("tg    #lambda");
1179     break;
1180   case 4:
1181     title.Append("C");
1182     break;
1183   }
1184
1185   gStyle->SetOptStat(0);
1186   TCanvas *c = new TCanvas("c","pulls",0,0,900,700);
1187   c->SetLogx();
1188   c->SetGridx();
1189   c->SetGridy();
1190
1191   TH2F *frame = new TH2F("frame",title.Data(),2,0.1,30,2,0,2);
1192   frame->SetXTitle("p_{T} [GeV/c]");
1193   frame->Draw();
1194   grA->Draw("P");
1195   grB->Draw("P");
1196   grC->Draw("P");
1197
1198   TLegend *leg = new TLegend(0.2,0.2,0.4,0.4);
1199   leg->AddEntry(grA,"|#eta|<0.3","p");
1200   leg->AddEntry(grB,"0.3<|#eta|<0.6","p");
1201   leg->AddEntry(grC,"0.6<|#eta|<0.9","p");
1202   leg->SetFillColor(0);
1203   leg->Draw();
1204
1205
1206   return;
1207 }
1208 //-----------------------------------------------------------------------------
1209 Int_t AliTPCtrackerParam::GetBin(Double_t pt,Double_t eta) const {
1210 //-----------------------------------------------------------------------------
1211 // This function tells bin number in a grid (pT,eta) 
1212 //-----------------------------------------------------------------------------
1213   if(TMath::Abs(eta)<0.3) {
1214     if(pt<0.3)            return 0;
1215     if(pt>=0.3 && pt<0.5) return 3;
1216     if(pt>=0.5 && pt<1.)  return 6;
1217     if(pt>=1. && pt<1.5)  return 9;
1218     if(pt>=1.5 && pt<3.)  return 12;
1219     if(pt>=3. && pt<5.)   return 15;
1220     if(pt>=5. && pt<7.)   return 18;
1221     if(pt>=7. && pt<15.)  return 21;
1222     if(pt>=15.)           return 24;
1223   }
1224   if(TMath::Abs(eta)>=0.3 && TMath::Abs(eta)<0.6) {
1225     if(pt<0.3)            return 1;
1226     if(pt>=0.3 && pt<0.5) return 4;
1227     if(pt>=0.5 && pt<1.)  return 7;
1228     if(pt>=1. && pt<1.5)  return 10;
1229     if(pt>=1.5 && pt<3.)  return 13;
1230     if(pt>=3. && pt<5.)   return 16;
1231     if(pt>=5. && pt<7.)   return 19;
1232     if(pt>=7. && pt<15.)  return 22;
1233     if(pt>=15.)           return 25;
1234   }
1235   if(TMath::Abs(eta)>=0.6) {
1236     if(pt<0.3)            return 2;
1237     if(pt>=0.3 && pt<0.5) return 5;
1238     if(pt>=0.5 && pt<1.)  return 8;
1239     if(pt>=1. && pt<1.5)  return 11;
1240     if(pt>=1.5 && pt<3.)  return 14;
1241     if(pt>=3. && pt<5.)   return 17;
1242     if(pt>=5. && pt<7.)   return 20;
1243     if(pt>=7. && pt<15.)  return 23;
1244     if(pt>=15.)           return 26;
1245   }
1246
1247   return -1;
1248
1249 }
1250 //-----------------------------------------------------------------------------
1251 TMatrixD AliTPCtrackerParam::GetSmearingMatrix(Double_t* cc,Double_t pt,
1252                                                Double_t eta) const {
1253 //-----------------------------------------------------------------------------
1254 // This function stretches the covariance matrix according to the pulls
1255 //-----------------------------------------------------------------------------
1256
1257   TMatrixD covMat(5,5);
1258
1259   covMat(0,0)=cc[0];
1260   covMat(1,0)=cc[1];  covMat(0,1)=covMat(1,0);
1261   covMat(1,1)=cc[2];
1262   covMat(2,0)=cc[3];  covMat(0,2)=covMat(2,0);
1263   covMat(2,1)=cc[4];  covMat(1,2)=covMat(2,1);
1264   covMat(2,2)=cc[5];
1265   covMat(3,0)=cc[6];  covMat(0,3)=covMat(3,0);
1266   covMat(3,1)=cc[7];  covMat(1,3)=covMat(3,1);
1267   covMat(3,2)=cc[8];  covMat(2,3)=covMat(3,2);
1268   covMat(3,3)=cc[9];
1269   covMat(4,0)=cc[10]; covMat(0,4)=covMat(4,0);
1270   covMat(4,1)=cc[11]; covMat(1,4)=covMat(4,1);
1271   covMat(4,2)=cc[12]; covMat(2,4)=covMat(4,2);
1272   covMat(4,3)=cc[13]; covMat(3,4)=covMat(4,3);
1273   covMat(4,4)=cc[14];
1274
1275   TMatrixD stretchMat(5,5);
1276   for(Int_t k=0;k<5;k++) {
1277     for(Int_t l=0;l<5;l++) {
1278       stretchMat(k,l)=0.;
1279     }
1280   }
1281
1282   for(Int_t i=0;i<5;i++) stretchMat(i,i) = 
1283                            TMath::Sqrt((fPulls+i)->GetValueAt(pt,eta)); 
1284
1285   TMatrixD mat(stretchMat,TMatrixD::kMult,covMat);
1286   TMatrixD covMatSmear(mat,TMatrixD::kMult,stretchMat);
1287
1288   return covMatSmear;
1289 }
1290 //-----------------------------------------------------------------------------
1291 void AliTPCtrackerParam::InitializeKineGrid(Option_t* which,Option_t* how) {
1292 //-----------------------------------------------------------------------------
1293 // This function initializes ([pt,eta] points) the data members AliTPCkineGrid
1294 // which = "DB"     -> initialize fDBgrid... members
1295 //         "eff"    -> initialize fEff... members
1296 //         "pulls"  -> initialize fPulls... members
1297 //         "dEdx"   -> initialize fdEdx... members
1298 // how =   "points" -> initialize with the points of the grid
1299 //         "bins"   -> initialize with the bins of the grid
1300 //-----------------------------------------------------------------------------
1301
1302   const char *points = strstr(how,"points");
1303   const char *bins   = strstr(how,"bins");
1304   const char *DB     = strstr(which,"DB");
1305   const char *eff    = strstr(which,"eff");
1306   const char *pulls  = strstr(which,"pulls");
1307   const char *dEdx   = strstr(which,"dEdx");
1308
1309
1310   Int_t nEta=0,nPtPi=0,nPtKa=0,nPtPr=0,nPtEl=0,nPtMu=0;
1311
1312   Double_t etaPoints[2] = {0.3,0.6};
1313   Double_t etaBins[3]   = {0.15,0.45,0.75};
1314   Double_t ptPoints8[8] = {0.3,0.5,1.,1.5,3.,5.,7.,15.};
1315   Double_t ptPoints5[5] = {0.3,0.5,1.,1.5,5.};
1316   Double_t ptBins9[9]   = {0.244,0.390,0.676,1.190,2.36,4.,6.,10.,20.};
1317   Double_t ptBins6[6]   = {0.244,0.390,0.676,1.190,2.36,10.};
1318
1319   Double_t *eta=0,*ptPi=0,*ptKa=0,*ptPr=0,*ptEl=0,*ptMu=0;
1320
1321   if(points) {
1322     nEta  = 2;
1323     nPtPi = 8;
1324     nPtKa = 5;
1325     nPtPr = 4;
1326     nPtEl = 8;
1327     nPtMu = 2;
1328     eta  = etaPoints;
1329     ptPi = ptPoints8;
1330     ptKa = ptPoints5;
1331     ptPr = ptPoints8;
1332     ptEl = ptPoints8;
1333     ptMu = ptPoints8;
1334   }
1335   if(bins) {
1336     nEta  = 3;
1337     nPtPi = 9;
1338     nPtKa = 6;
1339     nPtPr = 5;
1340     nPtEl = 9;
1341     nPtMu = 3;
1342     eta  = etaBins;
1343     ptPi = ptBins9;
1344     ptKa = ptBins6;
1345     ptPr = ptBins9;
1346     ptEl = ptBins9;
1347     ptMu = ptBins9;
1348   }
1349
1350   AliTPCkineGrid *dummy=0;
1351
1352   if(DB) {    
1353     dummy = new AliTPCkineGrid(nPtPi,nEta,ptPi,eta);
1354     new(&fDBgridPi) AliTPCkineGrid(*dummy);
1355     delete dummy;
1356     dummy = new AliTPCkineGrid(nPtKa,nEta,ptKa,eta);
1357     new(&fDBgridKa) AliTPCkineGrid(*dummy);
1358     delete dummy;
1359     dummy = new AliTPCkineGrid(nPtEl,nEta,ptEl,eta);
1360     new(&fDBgridEl) AliTPCkineGrid(*dummy);
1361     delete dummy;
1362   }
1363   if(eff) {    
1364     dummy = new AliTPCkineGrid(nPtPi,nEta,ptPi,eta);
1365     new(&fEffPi) AliTPCkineGrid(*dummy);
1366     delete dummy;
1367     dummy = new AliTPCkineGrid(nPtKa,nEta,ptKa,eta);
1368     new(&fEffKa) AliTPCkineGrid(*dummy);
1369     delete dummy;
1370     dummy = new AliTPCkineGrid(nPtPr,nEta,ptPr,eta);
1371     new(&fEffPr) AliTPCkineGrid(*dummy);
1372     delete dummy;
1373     dummy = new AliTPCkineGrid(nPtEl,nEta,ptEl,eta);
1374     new(&fEffEl) AliTPCkineGrid(*dummy);
1375     delete dummy;
1376     dummy = new AliTPCkineGrid(nPtMu,nEta,ptMu,eta);
1377     new(&fEffMu) AliTPCkineGrid(*dummy);
1378     delete dummy;
1379   }
1380   if(pulls) {    
1381     dummy = new AliTPCkineGrid(nPtPi,nEta,ptPi,eta);
1382     for(Int_t i=0;i<5;i++) new(&fPullsPi[i]) AliTPCkineGrid(*dummy);
1383     delete dummy;
1384     dummy = new AliTPCkineGrid(nPtKa,nEta,ptKa,eta);
1385     for(Int_t i=0;i<5;i++) new(&fPullsKa[i]) AliTPCkineGrid(*dummy);
1386     delete dummy;
1387     dummy = new AliTPCkineGrid(nPtEl,nEta,ptEl,eta);
1388     for(Int_t i=0;i<5;i++) new(&fPullsEl[i]) AliTPCkineGrid(*dummy);
1389     delete dummy;
1390   }
1391   if(dEdx) {    
1392     dummy = new AliTPCkineGrid(nPtPi,nEta,ptPi,eta);
1393     new(&fdEdxMeanPi) AliTPCkineGrid(*dummy);
1394     new(&fdEdxRMSPi) AliTPCkineGrid(*dummy);
1395     delete dummy;
1396     dummy = new AliTPCkineGrid(nPtKa,nEta,ptKa,eta);
1397     new(&fdEdxMeanKa) AliTPCkineGrid(*dummy);
1398     new(&fdEdxRMSKa) AliTPCkineGrid(*dummy);
1399     delete dummy;
1400     dummy = new AliTPCkineGrid(nPtPi,nEta,ptPi,eta);
1401     new(&fdEdxMeanPr) AliTPCkineGrid(*dummy);
1402     new(&fdEdxRMSPr) AliTPCkineGrid(*dummy);
1403     delete dummy;
1404     dummy = new AliTPCkineGrid(nPtEl,nEta,ptEl,eta);
1405     new(&fdEdxMeanEl) AliTPCkineGrid(*dummy);
1406     new(&fdEdxRMSEl) AliTPCkineGrid(*dummy);
1407     delete dummy;
1408   }
1409
1410   return;
1411 }
1412 //-----------------------------------------------------------------------------
1413 void AliTPCtrackerParam::MakeDataBase() {
1414 //-----------------------------------------------------------------------------
1415 // This function creates the DB file and store in it:
1416 //  - TPC Efficiencies for pi,ka,pr,el,mu     (AliTPCkineGrid class)
1417 //  - Pulls for pi,ka,el                      (AliTPCkineGrid class)
1418 //  - Regularization parameters for pi,ka,el  (TMatrixD class)
1419 //  - dE/dx parameterization for pi,ka,pr,el  (AliTPCkineGrid class)  
1420 //  - Regularized cov. matrices for pi,ka,el  (COVMATRIX structure)
1421 //-----------------------------------------------------------------------------
1422
1423   // define some file names
1424   Char_t *effFile  ="CovMatrixDB_partial.root";
1425   Char_t *pullsFile="CovMatrixDB_partial.root";
1426   Char_t *regPiFile="CovMatrixDB_partial.root";
1427   Char_t *regKaFile="CovMatrixDB_partial.root";
1428   Char_t *regElFile="CovMatrixDB_partial.root";
1429   Char_t *dEdxPiFile="dEdxPi.root";
1430   Char_t *dEdxKaFile="dEdxKa.root";
1431   Char_t *dEdxPrFile="dEdxPr.root";
1432   Char_t *dEdxElFile="dEdxEl.root";
1433   Char_t *cmFile1  ="CovMatrix_AllEvts_1.root";
1434   Char_t *cmFile2  ="CovMatrix_AllEvts_2.root";
1435   Char_t *cmFile3  ="CovMatrix_AllEvts_3.root";
1436
1437   // store the effieciencies
1438   ReadEffs(effFile);
1439   WriteEffs(fDBfileName.Data());
1440
1441   // store the pulls
1442   ReadPulls(pullsFile);
1443   WritePulls(fDBfileName.Data());
1444
1445   //store the regularization parameters
1446   ReadRegParams(regPiFile,211);
1447   WriteRegParams(fDBfileName.Data(),211);
1448   ReadRegParams(regKaFile,321);
1449   WriteRegParams(fDBfileName.Data(),321);
1450   ReadRegParams(regElFile,11);
1451   WriteRegParams(fDBfileName.Data(),11);
1452
1453   // store the dEdx parameters
1454   ReaddEdx(dEdxPiFile,211);
1455   WritedEdx(fDBfileName.Data(),211);
1456   ReaddEdx(dEdxKaFile,321);
1457   WritedEdx(fDBfileName.Data(),321);
1458   ReaddEdx(dEdxPrFile,2212);
1459   WritedEdx(fDBfileName.Data(),2212);
1460   ReaddEdx(dEdxElFile,11);
1461   WritedEdx(fDBfileName.Data(),11);
1462
1463
1464   //
1465   // store the regularized covariance matrices
1466   //
1467   InitializeKineGrid("DB","points");
1468
1469   const Int_t nBinsPi = fDBgridPi.GetTotBins();
1470   const Int_t nBinsKa = fDBgridKa.GetTotBins();
1471   const Int_t nBinsEl = fDBgridEl.GetTotBins();
1472
1473
1474   // create the trees for cov. matrices
1475   // trees for pions
1476   TTree *CovTreePi_ = NULL;
1477   CovTreePi_ = new TTree[nBinsPi]; 
1478   // trees for kaons
1479   TTree *CovTreeKa_ = NULL;
1480   CovTreeKa_ = new TTree[nBinsKa]; 
1481   // trees for electrons
1482   TTree *CovTreeEl_ = NULL;
1483   CovTreeEl_ = new TTree[nBinsEl]; 
1484
1485   Char_t hname[100], htitle[100];
1486   COVMATRIX covmat;
1487
1488   
1489   for(Int_t i=0; i<nBinsPi; i++) {
1490     sprintf(hname,"CovTreePi_bin%d",i);
1491     sprintf(htitle,"Tree with cov matrix elements for bin %d",i);
1492     CovTreePi_[i].SetName(hname); CovTreePi_[i].SetTitle(htitle);
1493     CovTreePi_[i].Branch("matrix",&covmat,"c00/D:c10:c11:c20:c21:c22:c30:c31:c32:c33:c40:c41:c42:c43:c44",5000000);
1494   }
1495   for(Int_t i=0; i<nBinsKa; i++) {
1496     sprintf(hname,"CovTreeKa_bin%d",i);
1497     sprintf(htitle,"Tree with cov matrix elements for bin %d",i);
1498     CovTreeKa_[i].SetName(hname); CovTreeKa_[i].SetTitle(htitle);
1499     CovTreeKa_[i].Branch("matrix",&covmat,"c00/D:c10:c11:c20:c21:c22:c30:c31:c32:c33:c40:c41:c42:c43:c44",1000000);
1500   }
1501   for(Int_t i=0; i<nBinsEl; i++) {
1502     sprintf(hname,"CovTreeEl_bin%d",i);
1503     sprintf(htitle,"Tree with cov matrix elements for bin %d",i);
1504     CovTreeEl_[i].SetName(hname); CovTreeEl_[i].SetTitle(htitle);
1505     CovTreeEl_[i].Branch("matrix",&covmat,"c00/D:c10:c11:c20:c21:c22:c30:c31:c32:c33:c40:c41:c42:c43:c44",1000000);
1506   }
1507   
1508   // create the chain with the compared tracks
1509   TChain cmptrkchain("cmptrktree");
1510   cmptrkchain.Add(cmFile1);
1511   cmptrkchain.Add(cmFile2);
1512   cmptrkchain.Add(cmFile3);
1513
1514   COMPTRACK cmptrk; 
1515   cmptrkchain.SetBranchAddress("comptracks",&cmptrk);
1516   Int_t entries = (Int_t)cmptrkchain.GetEntries(); 
1517   cerr<<" Number of entries: "<<entries<<endl;
1518
1519   Int_t trkPdg,trkBin;
1520   Double_t trkKine[2],trkRegPar[4]; 
1521   Int_t nPerBinPi[nBinsPi]; for(Int_t k=0;k<nBinsPi;k++) nPerBinPi[k]=0;
1522   Int_t nPerBinKa[nBinsKa]; for(Int_t k=0;k<nBinsKa;k++) nPerBinKa[k]=0;
1523   Int_t nPerBinEl[nBinsEl]; for(Int_t k=0;k<nBinsEl;k++) nPerBinEl[k]=0;
1524
1525   // loop on chain entries 
1526   for(Int_t l=0; l<entries; l++) {
1527     if(l % 10000 == 0) cerr<<"--- Processing track "<<l<<" of "<<entries<<" ---"<<endl;
1528     // get the event
1529     cmptrkchain.GetEvent(l);
1530     // get the pdg code
1531     trkPdg = TMath::Abs(cmptrk.pdg);
1532     // use only pions, kaons, electrons
1533     if(trkPdg!=211 && trkPdg!=321 && trkPdg!=11) continue;
1534     SetParticle(trkPdg);
1535     trkBin = fDBgrid->GetBin(cmptrk.pt,cmptrk.eta);
1536     //cerr<<cmptrk.pt<<"  "<<cmptrk.eta<<"  "<<trkBin<<endl;
1537
1538     if(trkPdg==211 && nPerBinPi[trkBin]>=5000) continue;
1539     if(trkPdg==321 && nPerBinKa[trkBin]>=5000) continue;
1540     if(trkPdg==11  && nPerBinEl[trkBin]>=5000) continue;
1541
1542     trkKine[0] = cmptrk.p;
1543     trkKine[1] = cmptrk.cosl;
1544     
1545     // get regularized covariance matrix
1546     for(Int_t k=0;k<4;k++) trkRegPar[k] = (*fRegPar)(0,k);
1547     covmat.c00 = cmptrk.c00/RegFunc(trkKine,trkRegPar);
1548     covmat.c10 = cmptrk.c10;
1549     for(Int_t k=0;k<4;k++) trkRegPar[k] = (*fRegPar)(1,k);
1550     covmat.c11 = cmptrk.c11/RegFunc(trkKine,trkRegPar);
1551     for(Int_t k=0;k<4;k++) trkRegPar[k] = (*fRegPar)(2,k);
1552     covmat.c20 = cmptrk.c20/RegFunc(trkKine,trkRegPar);
1553     covmat.c21 = cmptrk.c21;
1554     for(Int_t k=0;k<4;k++) trkRegPar[k] = (*fRegPar)(3,k);
1555     covmat.c22 = cmptrk.c22/RegFunc(trkKine,trkRegPar);
1556     covmat.c30 = cmptrk.c30;
1557     for(Int_t k=0;k<4;k++) trkRegPar[k] = (*fRegPar)(4,k);
1558     covmat.c31 = cmptrk.c31/RegFunc(trkKine,trkRegPar);
1559     covmat.c32 = cmptrk.c32;
1560     for(Int_t k=0;k<4;k++) trkRegPar[k] = (*fRegPar)(5,k);
1561     covmat.c33 = cmptrk.c33/RegFunc(trkKine,trkRegPar);
1562     for(Int_t k=0;k<4;k++) trkRegPar[k] = (*fRegPar)(6,k);
1563     covmat.c40 = cmptrk.c40/RegFunc(trkKine,trkRegPar);
1564     covmat.c41 = cmptrk.c41;
1565     for(Int_t k=0;k<4;k++) trkRegPar[k] = (*fRegPar)(7,k);
1566     covmat.c42 = cmptrk.c42/RegFunc(trkKine,trkRegPar);
1567     covmat.c43 = cmptrk.c43;
1568     for(Int_t k=0;k<4;k++) trkRegPar[k] = (*fRegPar)(8,k);
1569     covmat.c44 = cmptrk.c44/RegFunc(trkKine,trkRegPar);
1570
1571     // fill the tree
1572     switch (trkPdg) {
1573     case 211: // pions
1574       CovTreePi_[trkBin].Fill();
1575       nPerBinPi[trkBin]++;
1576       break;
1577     case 321: // kaons
1578       CovTreeKa_[trkBin].Fill();
1579       nPerBinKa[trkBin]++;
1580       break;
1581     case 11: // electrons
1582       CovTreeEl_[trkBin].Fill();
1583       nPerBinEl[trkBin]++;
1584       break;
1585     }
1586   } // loop on chain entries
1587
1588   // store all trees the DB file
1589   TFile *DBfile = new TFile(fDBfileName.Data(),"update");
1590   DBfile->mkdir("CovMatrices");
1591   gDirectory->cd("/CovMatrices");
1592   gDirectory->mkdir("Pions");
1593   gDirectory->mkdir("Kaons");
1594   gDirectory->mkdir("Electrons");
1595   // store pions
1596   gDirectory->cd("/CovMatrices/Pions");
1597   fDBgridPi.SetName("DBgridPi"); fDBgridPi.Write();
1598   for(Int_t i=0;i<nBinsPi;i++) CovTreePi_[i].Write();
1599   // store kaons
1600   gDirectory->cd("/CovMatrices/Kaons");
1601   fDBgridKa.SetName("DBgridKa"); fDBgridKa.Write();
1602   for(Int_t i=0;i<nBinsKa;i++) CovTreeKa_[i].Write();
1603   // store electrons
1604   gDirectory->cd("/CovMatrices/Electrons");
1605   fDBgridEl.SetName("DBgridEl"); fDBgridEl.Write();
1606   for(Int_t i=0;i<nBinsEl;i++) CovTreeEl_[i].Write();
1607
1608   DBfile->Close();
1609   
1610   return;
1611 }
1612 //-----------------------------------------------------------------------------
1613 void AliTPCtrackerParam::MergeEvents(Int_t evFirst=1,Int_t evLast=1) {
1614 //-----------------------------------------------------------------------------
1615 // This function: 1) merges the files from track comparison
1616 //                   (beware: better no more than 100 events per file)
1617 //                2) computes the average TPC efficiencies
1618 //-----------------------------------------------------------------------------
1619
1620   Char_t *outName="TPCeff.root";
1621
1622   Int_t    nCol;
1623   Float_t effPi,effKa,effPr,effEl,effMu;
1624   Float_t avEffPi[27],avEffKa[27],avEffPr[27],avEffEl[27],avEffMu[27];
1625   Int_t    evtsPi[27],evtsKa[27],evtsPr[27],evtsEl[27],evtsMu[27];
1626   for(Int_t j=0; j<27; j++) {
1627     avEffPi[j]=avEffKa[j]=avEffPr[j]=avEffEl[j]=avEffMu[j]=0.;
1628     evtsPi[j]=evtsKa[j]=evtsPr[j]=evtsEl[j]=evtsMu[j]=0;
1629   }
1630
1631   // create the chain for the tree of compared tracks
1632   TChain ch1("cmptrktree");
1633   TChain ch2("cmptrktree");
1634   TChain ch3("cmptrktree");
1635
1636
1637   for(Int_t j=evFirst; j<=evLast; j++) {
1638     cerr<<"Processing event "<<j<<endl;
1639
1640     TString covName("CovMatrix.");
1641     TString effName("TPCeff.");
1642     covName+=j;
1643     effName+=j;
1644     covName.Append(".root");
1645     effName.Append(".dat");
1646
1647     if(gSystem->AccessPathName(covName.Data(),kFileExists) ||
1648        gSystem->AccessPathName(effName.Data(),kFileExists)) continue;
1649
1650     if(j<=100) ch1.Add(covName.Data());
1651     if(j>100 && j<=200) ch2.Add(covName.Data());
1652     if(j>200) ch3.Add(covName.Data());
1653
1654
1655     FILE *effIn = fopen(effName.Data(),"r");
1656     for(Int_t k=0; k<27; k++) {
1657       nCol = fscanf(effIn,"%f  %f  %f  %f  %f",&effPi,&effKa,&effPr,&effEl,&effMu);
1658       if(effPi>=0.) {avEffPi[k]+=effPi; evtsPi[k]++;}
1659       if(effKa>=0.) {avEffKa[k]+=effKa; evtsKa[k]++;}
1660       if(effPr>=0.) {avEffPr[k]+=effPr; evtsPr[k]++;}
1661       if(effEl>=0.) {avEffEl[k]+=effEl; evtsEl[k]++;}
1662       if(effMu>=0.) {avEffMu[k]+=effMu; evtsMu[k]++;}
1663     }
1664     fclose(effIn);
1665
1666   } // loop on events
1667
1668   // merge chain in one file
1669   TFile *covOut=0;
1670   covOut = new TFile("CovMatrix_AllEvts_1.root","recreate");
1671   ch1.Merge(covOut,1000000000);
1672   covOut->Close();
1673   delete covOut;
1674   covOut = new TFile("CovMatrix_AllEvts_2.root","recreate");
1675   ch2.Merge(covOut,1000000000);
1676   covOut->Close();
1677   delete covOut;
1678   covOut = new TFile("CovMatrix_AllEvts_3.root","recreate");
1679   ch3.Merge(covOut,1000000000);
1680   covOut->Close();
1681   delete covOut;
1682
1683
1684   // efficiencies 
1685
1686   InitializeKineGrid("eff","bins");
1687
1688   // pions
1689   for(Int_t j=0; j<27; j++) {
1690     if(evtsPi[j]==0) evtsPi[j]++;
1691     fEffPi.SetParam(j,(Double_t)avEffPi[j]/evtsPi[j]);
1692   }
1693  
1694   // kaons
1695   Int_t j=0;
1696   for(Int_t k=0; k<27; k++) {
1697     if(k>14 && k!=21 && k!=22 && k!=23) continue;
1698     if(evtsKa[k]==0) evtsKa[k]++;
1699     fEffKa.SetParam(j,(Double_t)avEffKa[k]/evtsKa[k]);
1700     j++;
1701   }
1702
1703   // protons
1704   for(Int_t j=0; j<15; j++) {
1705     if(evtsPr[j]==0) evtsPr[j]++;
1706     fEffPr.SetParam(j,(Double_t)avEffPr[j]/evtsPr[j]);
1707   }
1708
1709   // electrons
1710   for(Int_t j=0; j<27; j++) {
1711     if(evtsEl[j]==0) evtsEl[j]++;
1712     fEffEl.SetParam(j,(Double_t)avEffEl[j]/evtsEl[j]);
1713   }
1714
1715   // muons
1716   for(Int_t j=0; j<9; j++) {
1717     if(evtsMu[j]==0) evtsMu[j]++;
1718     fEffMu.SetParam(j,(Double_t)avEffMu[j]/evtsMu[j]);
1719   }
1720
1721   // write efficiencies to a file
1722   WriteEffs(outName);
1723
1724   return;
1725 }
1726 //-----------------------------------------------------------------------------
1727 Int_t AliTPCtrackerParam::ReadAllData(const Char_t* inName) {
1728 //-----------------------------------------------------------------------------
1729 // This function reads all parameters from the DB
1730 //-----------------------------------------------------------------------------
1731
1732   if(!ReadEffs(inName)) return 0;
1733   if(!ReadPulls(inName)) return 0;        
1734   if(!ReadRegParams(inName,211)) return 0; 
1735   if(!ReadRegParams(inName,321)) return 0; 
1736   if(!ReadRegParams(inName,11)) return 0; 
1737   if(!ReaddEdx(inName,211)) return 0;
1738   if(!ReaddEdx(inName,321)) return 0;
1739   if(!ReaddEdx(inName,2212)) return 0;
1740   if(!ReaddEdx(inName,11)) return 0;
1741   if(!ReadDBgrid(inName)) return 0;
1742
1743   return 1;
1744 }
1745 //-----------------------------------------------------------------------------
1746 Int_t AliTPCtrackerParam::ReaddEdx(const Char_t* inName,Int_t pdg) {
1747 //-----------------------------------------------------------------------------
1748 // This function reads the dEdx parameters from the DB
1749 //-----------------------------------------------------------------------------
1750
1751   if(gSystem->AccessPathName(inName,kFileExists)) { 
1752     cerr<<"AliTPCtrackerParam::ReaddEdx: "<<inName<<" not found\n"; 
1753     return 0; }
1754
1755   TFile *inFile = TFile::Open(inName);
1756   switch (pdg) {
1757   case 211:
1758     fdEdxMean = (AliTPCkineGrid*)inFile->Get("/dEdx/Pions/dEdxMeanPi");
1759     fdEdxMeanPi.~AliTPCkineGrid();
1760     new(&fdEdxMeanPi) AliTPCkineGrid(*fdEdxMean);
1761     fdEdxRMS = (AliTPCkineGrid*)inFile->Get("/dEdx/Pions/dEdxRMSPi");
1762     fdEdxRMSPi.~AliTPCkineGrid();
1763     new(&fdEdxRMSPi) AliTPCkineGrid(*fdEdxRMS);
1764     break;
1765   case 321:
1766     fdEdxMean = (AliTPCkineGrid*)inFile->Get("/dEdx/Kaons/dEdxMeanKa");
1767     fdEdxMeanKa.~AliTPCkineGrid();
1768     new(&fdEdxMeanKa) AliTPCkineGrid(*fdEdxMean);
1769     fdEdxRMS = (AliTPCkineGrid*)inFile->Get("/dEdx/Kaons/dEdxRMSKa");
1770     fdEdxRMSKa.~AliTPCkineGrid();
1771     new(&fdEdxRMSKa) AliTPCkineGrid(*fdEdxRMS);
1772     break;
1773   case 2212:
1774     fdEdxMean = (AliTPCkineGrid*)inFile->Get("/dEdx/Protons/dEdxMeanPr");
1775     fdEdxMeanPr.~AliTPCkineGrid();
1776     new(&fdEdxMeanPr) AliTPCkineGrid(*fdEdxMean);
1777     fdEdxRMS = (AliTPCkineGrid*)inFile->Get("/dEdx/Protons/dEdxRMSPr");
1778     fdEdxRMSPr.~AliTPCkineGrid();
1779     new(&fdEdxRMSPr) AliTPCkineGrid(*fdEdxRMS);
1780     break;
1781   case 11:
1782     fdEdxMean = (AliTPCkineGrid*)inFile->Get("/dEdx/Electrons/dEdxMeanEl");
1783     fdEdxMeanEl.~AliTPCkineGrid();
1784     new(&fdEdxMeanEl) AliTPCkineGrid(*fdEdxMean);
1785     fdEdxRMS = (AliTPCkineGrid*)inFile->Get("/dEdx/Electrons/dEdxRMSEl");
1786     fdEdxRMSEl.~AliTPCkineGrid();
1787     new(&fdEdxRMSEl) AliTPCkineGrid(*fdEdxRMS);
1788     break;
1789   }
1790   inFile->Close();
1791
1792   return 1;
1793 }
1794 //-----------------------------------------------------------------------------
1795 Int_t AliTPCtrackerParam::ReadDBgrid(const Char_t* inName) {
1796 //-----------------------------------------------------------------------------
1797 // This function reads the kine grid from the DB
1798 //-----------------------------------------------------------------------------
1799
1800   if(gSystem->AccessPathName(inName,kFileExists)) {
1801     cerr<<"AliTPCtrackerParam::ReadCovMatrices: "<<inName<<" not found\n"; 
1802     return 0; }
1803
1804   TFile *inFile = TFile::Open(inName);
1805
1806   // first read the DB grid for the different particles
1807   fDBgrid = (AliTPCkineGrid*)inFile->Get("/CovMatrices/Pions/DBgridPi");
1808   fDBgridPi.~AliTPCkineGrid();
1809   new(&fDBgridPi) AliTPCkineGrid(*fDBgrid);
1810   fDBgrid = (AliTPCkineGrid*)inFile->Get("/CovMatrices/Kaons/DBgridKa");
1811   fDBgridKa.~AliTPCkineGrid();
1812   new(&fDBgridKa) AliTPCkineGrid(*fDBgrid);
1813   fDBgrid = (AliTPCkineGrid*)inFile->Get("/CovMatrices/Electrons/DBgridEl");
1814   fDBgridEl.~AliTPCkineGrid();
1815   new(&fDBgridEl) AliTPCkineGrid(*fDBgrid);
1816
1817   inFile->Close();
1818
1819   return 1;
1820 }
1821 //-----------------------------------------------------------------------------
1822 Int_t AliTPCtrackerParam::ReadEffs(const Char_t* inName) {
1823 //-----------------------------------------------------------------------------
1824 // This function reads the TPC efficiencies from the DB
1825 //-----------------------------------------------------------------------------
1826
1827   if(gSystem->AccessPathName(inName,kFileExists)) {
1828     cerr<<"AliTPCtrackerParam::ReadEffs: "<<inName<<" not found\n"; 
1829     return 0; }
1830
1831   TFile *inFile = TFile::Open(inName);
1832
1833   fEff = (AliTPCkineGrid*)inFile->Get("/Efficiencies/Pions/EffPi");
1834   fEffPi.~AliTPCkineGrid();
1835   new(&fEffPi) AliTPCkineGrid(*fEff);
1836   fEff = (AliTPCkineGrid*)inFile->Get("/Efficiencies/Kaons/EffKa");
1837   fEffKa.~AliTPCkineGrid();
1838   new(&fEffKa) AliTPCkineGrid(*fEff);
1839   fEff = (AliTPCkineGrid*)inFile->Get("/Efficiencies/Protons/EffPr");
1840   fEffPr.~AliTPCkineGrid();
1841   new(&fEffPr) AliTPCkineGrid(*fEff);
1842   fEff = (AliTPCkineGrid*)inFile->Get("/Efficiencies/Electrons/EffEl");
1843   fEffEl.~AliTPCkineGrid();
1844   new(&fEffEl) AliTPCkineGrid(*fEff);
1845   fEff = (AliTPCkineGrid*)inFile->Get("/Efficiencies/Muons/EffMu");
1846   fEffMu.~AliTPCkineGrid();
1847   new(&fEffMu) AliTPCkineGrid(*fEff);
1848
1849   inFile->Close();
1850
1851   return 1;
1852 }
1853 //-----------------------------------------------------------------------------
1854 Int_t AliTPCtrackerParam::ReadPulls(const Char_t* inName) {
1855 //-----------------------------------------------------------------------------
1856 // This function reads the pulls from the DB
1857 //-----------------------------------------------------------------------------
1858
1859   if(gSystem->AccessPathName(inName,kFileExists)) { 
1860     cerr<<"AliTPCtrackerParam::ReadPulls: "<<inName<<" not found\n"; 
1861     return 0; }
1862
1863   TFile *inFile = TFile::Open(inName);
1864
1865   for(Int_t i=0; i<5; i++) {
1866     TString pi("/Pulls/Pions/PullsPi_"); pi+=i;
1867     TString ka("/Pulls/Kaons/PullsKa_"); ka+=i;
1868     TString el("/Pulls/Electrons/PullsEl_"); el+=i;
1869     fPulls = (AliTPCkineGrid*)inFile->Get(pi.Data());
1870     fPullsPi[i].~AliTPCkineGrid();
1871     new(&fPullsPi[i]) AliTPCkineGrid(*fPulls);
1872     fPulls = (AliTPCkineGrid*)inFile->Get(ka.Data());
1873     fPullsKa[i].~AliTPCkineGrid();
1874     new(&fPullsKa[i]) AliTPCkineGrid(*fPulls);
1875     fPulls = (AliTPCkineGrid*)inFile->Get(el.Data());
1876     fPullsEl[i].~AliTPCkineGrid();
1877     new(&fPullsEl[i]) AliTPCkineGrid(*fPulls);
1878   }
1879
1880   inFile->Close();
1881
1882   return 1;
1883 }
1884 //-----------------------------------------------------------------------------
1885 Int_t AliTPCtrackerParam::ReadRegParams(const Char_t* inName,Int_t pdg) {
1886 //-----------------------------------------------------------------------------
1887 // This function reads the regularization parameters from the DB
1888 //-----------------------------------------------------------------------------
1889
1890   if(gSystem->AccessPathName(inName,kFileExists)) { 
1891     cerr<<"AliTPCtrackerParam::ReadRegParams: "<<inName<<" not found\n"; 
1892     return 0; }
1893
1894   TFile *inFile = TFile::Open(inName);
1895   switch (pdg) {
1896   case 211:
1897     fRegPar = (TMatrixD*)inFile->Get("/RegParams/Pions/RegPions");
1898     new(&fRegParPi) TMatrixD(*fRegPar);
1899     break;
1900   case 321:
1901     fRegPar = (TMatrixD*)inFile->Get("/RegParams/Kaons/RegKaons");
1902     new(&fRegParKa) TMatrixD(*fRegPar);
1903     break;
1904   case 11:
1905     fRegPar = (TMatrixD*)inFile->Get("/RegParams/Electrons/RegElectrons");
1906     new(&fRegParEl) TMatrixD(*fRegPar);
1907     break;
1908   }
1909   inFile->Close();
1910
1911   return 1;
1912 }
1913 //-----------------------------------------------------------------------------
1914 void AliTPCtrackerParam::RegularizeCovMatrix(const Char_t *outName,Int_t pdg) {
1915 //-----------------------------------------------------------------------------
1916 // This function regularizes the elements of the covariance matrix
1917 // that show a momentum depence:
1918 // c00, c11, c22, c33, c44, c20, c24, c40, c31 
1919 // The regularization is done separately for pions, kaons, electrons:
1920 // give "Pion","Kaon" or "Electron" as first argument.
1921 //-----------------------------------------------------------------------------
1922
1923   gStyle->SetOptStat(0);
1924   gStyle->SetOptFit(10001);
1925
1926   Int_t thisPdg=211;
1927   Char_t *part="Pions - ";
1928
1929   InitializeKineGrid("DB","points");
1930   SetParticle(pdg);
1931   const Int_t fitbins = fDBgrid->GetBinsPt();
1932   cerr<<" Fit bins:  "<<fitbins<<endl;
1933
1934   switch (pdg) {
1935   case 211: // pions
1936     thisPdg=211;  
1937     part="Pions - ";
1938     cerr<<" Processing pions ...\n";
1939     break;
1940   case 321: //kaons
1941     thisPdg=321; 
1942     part="Kaons - ";
1943     cerr<<" Processing kaons ...\n";
1944     break;
1945   case 11: // electrons
1946     thisPdg= 11; 
1947     part="Electrons - ";
1948     cerr<<" Processing electrons ...\n";
1949     break;
1950   }
1951
1952   // create the chain with all compared tracks
1953   TChain cmptrkchain("cmptrktree");
1954   cmptrkchain.Add("CovMatrix_AllEvts_1.root");
1955   cmptrkchain.Add("CovMatrix_AllEvts_2.root");
1956   cmptrkchain.Add("CovMatrix_AllEvts_3.root");
1957
1958   COMPTRACK cmptrk; 
1959   cmptrkchain.SetBranchAddress("comptracks",&cmptrk);
1960   Int_t entries = (Int_t)cmptrkchain.GetEntries(); 
1961
1962     
1963   Int_t pbin;
1964   Int_t n[fitbins];
1965   Int_t n00[fitbins],n11[fitbins],n20[fitbins],n22[fitbins],n31[fitbins],n33[fitbins],n40[fitbins],n42[fitbins],n44[fitbins];
1966   Double_t p[fitbins],ep[fitbins];
1967   Double_t mean00[fitbins],mean11[fitbins],mean20[fitbins],mean22[fitbins],mean31[fitbins],mean33[fitbins],mean40[fitbins],mean42[fitbins],mean44[fitbins];
1968   Double_t sigma00[fitbins],sigma11[fitbins],sigma20[fitbins],sigma22[fitbins],sigma31[fitbins],sigma33[fitbins],sigma40[fitbins],sigma42[fitbins],sigma44[fitbins];
1969   Double_t rmean[fitbins],rsigma[fitbins];
1970   Double_t fitpar[4];
1971
1972   for(Int_t l=0; l<fitbins; l++) {
1973     n[l]=1;
1974     n00[l]=n11[l]=n20[l]=n22[l]=n31[l]=n33[l]=n40[l]=n42[l]=n44[l]=1;
1975     p[l ]=ep[l]=0.;
1976     mean00[l]=mean11[l]=mean20[l]=mean22[l]=mean31[l]=mean33[l]=mean40[l]=mean42[l]=mean44[l]=0.;
1977     sigma00[l]=sigma11[l]=sigma20[l]=sigma22[l]=sigma31[l]=sigma33[l]=sigma40[l]=sigma42[l]=sigma44[l]=0.;
1978   }
1979
1980   // loop on chain entries for mean 
1981   for(Int_t l=0; l<entries; l++) {
1982     cmptrkchain.GetEvent(l);
1983     if(TMath::Abs(cmptrk.pdg)!=thisPdg) continue;
1984     pbin = (Int_t)fDBgrid->GetBin(cmptrk.pt,cmptrk.eta)/fDBgrid->GetBinsEta();
1985     n[pbin]++;
1986     p[pbin]+=cmptrk.p;
1987     mean00[pbin]+=cmptrk.c00;
1988     mean11[pbin]+=cmptrk.c11;
1989     mean20[pbin]+=cmptrk.c20;
1990     mean22[pbin]+=cmptrk.c22;
1991     mean31[pbin]+=cmptrk.c31;
1992     mean33[pbin]+=cmptrk.c33;
1993     mean40[pbin]+=cmptrk.c40;
1994     mean42[pbin]+=cmptrk.c42;
1995     mean44[pbin]+=cmptrk.c44;
1996   } // loop on chain entries
1997
1998   for(Int_t l=0; l<fitbins; l++) {
1999     p[l]/=n[l];
2000     mean00[l]/=n[l];
2001     mean11[l]/=n[l];
2002     mean20[l]/=n[l];
2003     mean22[l]/=n[l];
2004     mean31[l]/=n[l];
2005     mean33[l]/=n[l];
2006     mean40[l]/=n[l];
2007     mean42[l]/=n[l];
2008     mean44[l]/=n[l];
2009   }
2010
2011   // loop on chain entries for sigma
2012   for(Int_t l=0; l<entries; l++) {
2013     cmptrkchain.GetEvent(l);
2014     if(TMath::Abs(cmptrk.pdg)!=thisPdg) continue;
2015     pbin = (Int_t)fDBgrid->GetBin(cmptrk.pt,cmptrk.eta)/fDBgrid->GetBinsEta();
2016     if(TMath::Abs(cmptrk.c00-mean00[pbin])<0.4*mean00[pbin]) { n00[pbin]++;
2017       sigma00[pbin]+=(cmptrk.c00-mean00[pbin])*(cmptrk.c00-mean00[pbin]); }
2018     if(TMath::Abs(cmptrk.c11-mean11[pbin])<0.4*mean11[pbin]) { n11[pbin]++;
2019       sigma11[pbin]+=(cmptrk.c11-mean11[pbin])*(cmptrk.c11-mean11[pbin]); }
2020     if(TMath::Abs(cmptrk.c20-mean20[pbin])<0.4*mean20[pbin]) { n20[pbin]++;
2021       sigma20[pbin]+=(cmptrk.c20-mean20[pbin])*(cmptrk.c20-mean20[pbin]); }
2022     if(TMath::Abs(cmptrk.c22-mean22[pbin])<0.4*mean22[pbin]) { n22[pbin]++;
2023       sigma22[pbin]+=(cmptrk.c22-mean22[pbin])*(cmptrk.c22-mean22[pbin]); }
2024     if(TMath::Abs(cmptrk.c31-mean31[pbin])<-0.4*mean31[pbin]) { n31[pbin]++;
2025       sigma31[pbin]+=(cmptrk.c31-mean31[pbin])*(cmptrk.c31-mean31[pbin]); }
2026     if(TMath::Abs(cmptrk.c33-mean33[pbin])<0.4*mean33[pbin]) { n33[pbin]++;
2027       sigma33[pbin]+=(cmptrk.c33-mean33[pbin])*(cmptrk.c33-mean33[pbin]); }
2028     if(TMath::Abs(cmptrk.c40-mean40[pbin])<0.4*mean40[pbin]) { n40[pbin]++;
2029       sigma40[pbin]+=(cmptrk.c40-mean40[pbin])*(cmptrk.c40-mean40[pbin]); }
2030     if(TMath::Abs(cmptrk.c42-mean42[pbin])<0.4*mean42[pbin]) { n42[pbin]++;
2031       sigma42[pbin]+=(cmptrk.c42-mean42[pbin])*(cmptrk.c42-mean42[pbin]); }
2032     if(TMath::Abs(cmptrk.c44-mean44[pbin])<0.4*mean44[pbin]) { n44[pbin]++;
2033       sigma44[pbin]+=(cmptrk.c44-mean44[pbin])*(cmptrk.c44-mean44[pbin]); }
2034   } // loop on chain entries
2035  
2036   for(Int_t l=0; l<fitbins; l++) {
2037     sigma00[l] = TMath::Sqrt(sigma00[l]/n00[l]);
2038     sigma11[l] = TMath::Sqrt(sigma11[l]/n11[l]);
2039     sigma20[l] = TMath::Sqrt(sigma20[l]/n20[l]);
2040     sigma22[l] = TMath::Sqrt(sigma22[l]/n22[l]);
2041     sigma31[l] = TMath::Sqrt(sigma31[l]/n31[l]);
2042     sigma33[l] = TMath::Sqrt(sigma33[l]/n33[l]);
2043     sigma40[l] = TMath::Sqrt(sigma40[l]/n40[l]);
2044     sigma42[l] = TMath::Sqrt(sigma42[l]/n42[l]);
2045     sigma44[l] = TMath::Sqrt(sigma44[l]/n44[l]);
2046   }
2047   
2048
2049   // Fit function
2050   TF1 *func = new TF1("FitRegFunc",FitRegFunc,0.23,50.,4);
2051   func->SetParNames("A_meas","A_scatt","B1","B2");
2052
2053   // line to draw on the plots
2054   TLine *lin = new TLine(-1,1,1.69,1);
2055   lin->SetLineStyle(2);
2056   lin->SetLineWidth(2);
2057
2058   // matrix used to store fit results
2059   TMatrixD fitRes(9,4);
2060
2061   //    --- c00 ---
2062
2063   // create the canvas
2064   TCanvas *canv00 = new TCanvas("canv00","c00",0,0,700,900); 
2065   canv00->Divide(1,2);
2066   // create the graph for cov matrix
2067   TGraphErrors *gr00 = new TGraphErrors(fitbins,p,mean00,ep,sigma00);
2068   TString title00("C(y,y)"); title00.Prepend(part);
2069   TH2F *frame00 = new TH2F("frame00",title00.Data(),2,0.1,50,2,0,5e-3);
2070   frame00->SetXTitle("p [GeV/c]");
2071   canv00->cd(1);  gPad->SetLogx();
2072   frame00->Draw();
2073   gr00->Draw("P");
2074   // Sets initial values for parameters
2075   func->SetParameters(1.6e-3,1.9e-4,1.5,0.);
2076   // Fit points in range defined by function
2077   gr00->Fit("FitRegFunc","R,Q");
2078   func->GetParameters(fitpar);
2079   for(Int_t i=0; i<4; i++) fitRes(0,i)=fitpar[i];
2080   for(Int_t l=0; l<fitbins; l++) {
2081     rmean[l]  = mean00[l]/FitRegFunc(&p[l],fitpar);
2082     rsigma[l] = sigma00[l]/FitRegFunc(&p[l],fitpar);
2083   }
2084   // create the graph the regularized cov. matrix
2085   TGraphErrors *gr00reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
2086   TString regtitle00("C(y,y)/(A_meas+A_scatt/p^{B1}/cos^{B2} #lambda"); 
2087   regtitle00.Prepend(part);
2088   TH2F *frame00reg = new TH2F("frame00reg",regtitle00.Data(),2,0.1,50,2,0,2);
2089   frame00reg->SetXTitle("p [GeV/c]");
2090   canv00->cd(2); gPad->SetLogx();
2091   frame00reg->Draw();
2092   gr00reg->Draw("P");
2093   lin->Draw("same");
2094
2095
2096   //    --- c11 ---
2097
2098   // create the canvas
2099   TCanvas *canv11 = new TCanvas("canv11","c11",0,0,700,900); 
2100   canv11->Divide(1,2);
2101   // create the graph for cov matrix
2102   TGraphErrors *gr11 = new TGraphErrors(fitbins,p,mean11,ep,sigma11);
2103   TString title11("C(z,z)"); title11.Prepend(part);
2104   TH2F *frame11 = new TH2F("frame11",title11.Data(),2,0.1,50,2,0,6e-3);
2105   frame11->SetXTitle("p [GeV/c]");
2106   canv11->cd(1);  gPad->SetLogx();
2107   frame11->Draw();
2108   gr11->Draw("P");
2109   // Sets initial values for parameters
2110   func->SetParameters(1.2e-3,8.1e-4,1.,0.);
2111   // Fit points in range defined by function
2112   gr11->Fit("FitRegFunc","R,Q");
2113   func->GetParameters(fitpar);
2114   for(Int_t i=0; i<4; i++) fitRes(1,i)=fitpar[i];
2115   for(Int_t l=0; l<fitbins; l++) {
2116     rmean[l]  = mean11[l]/FitRegFunc(&p[l],fitpar);
2117     rsigma[l] = sigma11[l]/FitRegFunc(&p[l],fitpar);
2118   }
2119   // create the graph the regularized cov. matrix
2120   TGraphErrors *gr11reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
2121   TString regtitle11("C(z,z)/(A_meas+A_scatt/p^{B1}/cos^{B2} #lambda"); 
2122   regtitle11.Prepend(part);
2123   TH2F *frame11reg = new TH2F("frame11reg",regtitle11.Data(),2,0.1,50,2,0,2);
2124   frame11reg->SetXTitle("p [GeV/c]");
2125   canv11->cd(2); gPad->SetLogx();
2126   frame11reg->Draw();
2127   gr11reg->Draw("P");
2128   lin->Draw("same");
2129
2130
2131   //    --- c20 ---
2132
2133   // create the canvas
2134   TCanvas *canv20 = new TCanvas("canv20","c20",0,0,700,900); 
2135   canv20->Divide(1,2);
2136   // create the graph for cov matrix
2137   TGraphErrors *gr20 = new TGraphErrors(fitbins,p,mean20,ep,sigma20);
2138   TString title20("C(#eta, y)"); title20.Prepend(part);
2139   TH2F *frame20 = new TH2F("frame20",title20.Data(),2,0.1,50,2,0,2.5e-4);
2140   frame20->SetXTitle("p [GeV/c]");
2141   canv20->cd(1);  gPad->SetLogx();
2142   frame20->Draw();
2143   gr20->Draw("P");
2144   // Sets initial values for parameters
2145   func->SetParameters(7.3e-5,1.2e-5,1.5,0.);
2146   // Fit points in range defined by function
2147   gr20->Fit("FitRegFunc","R,Q");
2148   func->GetParameters(fitpar);
2149   for(Int_t i=0; i<4; i++) fitRes(2,i)=fitpar[i];
2150   for(Int_t l=0; l<fitbins; l++) {
2151     rmean[l]  = mean20[l]/FitRegFunc(&p[l],fitpar);
2152     rsigma[l] = sigma20[l]/FitRegFunc(&p[l],fitpar);
2153   }
2154   // create the graph the regularized cov. matrix
2155   TGraphErrors *gr20reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
2156   TString regtitle20("C(#eta, y)/(A_meas+A_scatt/p^{B1}/cos^{B2} #lambda"); 
2157   regtitle20.Prepend(part);
2158   TH2F *frame20reg = new TH2F("frame20reg",regtitle20.Data(),2,0.1,50,2,0,2);
2159   frame20reg->SetXTitle("p [GeV/c]");
2160   canv20->cd(2); gPad->SetLogx();
2161   frame20reg->Draw();
2162   gr20reg->Draw("P");
2163   lin->Draw("same");
2164
2165
2166   //    --- c22 ---
2167
2168   // create the canvas
2169   TCanvas *canv22 = new TCanvas("canv22","c22",0,0,700,900); 
2170   canv22->Divide(1,2);
2171   // create the graph for cov matrix
2172   TGraphErrors *gr22 = new TGraphErrors(fitbins,p,mean22,ep,sigma22);
2173   TString title22("C(#eta, #eta)"); title22.Prepend(part);
2174   TH2F *frame22 = new TH2F("frame22",title22.Data(),2,0.1,50,2,0,3e-5);
2175   frame22->SetXTitle("p [GeV/c]");
2176   canv22->cd(1);  gPad->SetLogx();
2177   frame22->Draw();
2178   gr22->Draw("P");
2179   // Sets initial values for parameters
2180   func->SetParameters(5.2e-6,1.1e-6,2.,1.);
2181   // Fit points in range defined by function
2182   gr22->Fit("FitRegFunc","R,Q");
2183   func->GetParameters(fitpar);
2184   for(Int_t i=0; i<4; i++) fitRes(3,i)=fitpar[i];
2185   for(Int_t l=0; l<fitbins; l++) {
2186     rmean[l]  = mean22[l]/FitRegFunc(&p[l],fitpar);
2187     rsigma[l] = sigma22[l]/FitRegFunc(&p[l],fitpar);
2188   }
2189   // create the graph the regularized cov. matrix
2190   TGraphErrors *gr22reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
2191   TString regtitle22("C(#eta, #eta)/(A_meas+A_scatt/p^{B1}/cos^{B2} #lambda"); 
2192   regtitle22.Prepend(part);
2193   TH2F *frame22reg = new TH2F("frame22reg",regtitle22.Data(),2,0.1,50,2,0,2);
2194   frame22reg->SetXTitle("p [GeV/c]");
2195   canv22->cd(2); gPad->SetLogx();
2196   frame22reg->Draw();
2197   gr22reg->Draw("P");
2198   lin->Draw("same");
2199
2200
2201   //    --- c31 ---
2202
2203   // create the canvas
2204   TCanvas *canv31 = new TCanvas("canv31","c31",0,0,700,900); 
2205   canv31->Divide(1,2);
2206   // create the graph for cov matrix
2207   TGraphErrors *gr31 = new TGraphErrors(fitbins,p,mean31,ep,sigma31);
2208   TString title31("C(tg #lambda,z)"); title31.Prepend(part);
2209   TH2F *frame31 = new TH2F("frame31",title31.Data(),2,0.1,50,2,-2e-4,0);
2210   frame31->SetXTitle("p [GeV/c]");
2211   canv31->cd(1);  gPad->SetLogx();
2212   frame31->Draw();
2213   gr31->Draw("P");
2214   // Sets initial values for parameters
2215   func->SetParameters(-1.2e-5,-1.2e-5,1.5,3.);
2216   // Fit points in range defined by function
2217   gr31->Fit("FitRegFunc","R,Q");
2218   func->GetParameters(fitpar);
2219   for(Int_t i=0; i<4; i++) fitRes(4,i)=fitpar[i];
2220   for(Int_t l=0; l<fitbins; l++) {
2221     rmean[l]  = mean31[l]/FitRegFunc(&p[l],fitpar);
2222     rsigma[l] = -sigma31[l]/FitRegFunc(&p[l],fitpar);
2223   }
2224   // create the graph the regularized cov. matrix
2225   TGraphErrors *gr31reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
2226   TString regtitle31("C(tg #lambda,z)/(A_meas+A_scatt/p^{B1}/cos^{B2} #lambda"); 
2227   regtitle31.Prepend(part);
2228   TH2F *frame31reg = new TH2F("frame31reg",regtitle31.Data(),2,0.1,50,2,0,2);
2229   frame31reg->SetXTitle("p [GeV/c]");
2230   canv31->cd(2); gPad->SetLogx();
2231   frame31reg->Draw();
2232   gr31reg->Draw("P");
2233   lin->Draw("same");
2234
2235
2236   //    --- c33 ---
2237
2238   // create the canvas
2239   TCanvas *canv33 = new TCanvas("canv33","c33",0,0,700,900); 
2240   canv33->Divide(1,2);
2241   // create the graph for cov matrix
2242   TGraphErrors *gr33 = new TGraphErrors(fitbins,p,mean33,ep,sigma33);
2243   TString title33("C(tg #lambda,tg #lambda)"); title33.Prepend(part);
2244   TH2F *frame33 = new TH2F("frame33",title33.Data(),2,0.1,50,2,0,1e-5);
2245   frame33->SetXTitle("p [GeV/c]");
2246   canv33->cd(1);  gPad->SetLogx();
2247   frame33->Draw();
2248   gr33->Draw("P");
2249   // Sets initial values for parameters
2250   func->SetParameters(1.3e-7,4.6e-7,1.7,4.);
2251   // Fit points in range defined by function
2252   gr33->Fit("FitRegFunc","R,Q");
2253   func->GetParameters(fitpar);
2254   for(Int_t i=0; i<4; i++) fitRes(5,i)=fitpar[i];
2255   for(Int_t l=0; l<fitbins; l++) {
2256     rmean[l]  = mean33[l]/FitRegFunc(&p[l],fitpar);
2257     rsigma[l] = sigma33[l]/FitRegFunc(&p[l],fitpar);
2258   }
2259   // create the graph the regularized cov. matrix
2260   TGraphErrors *gr33reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
2261   TString regtitle33("C(tg #lambda,tg #lambda)/(A_meas+A_scatt/p^{B1}/cos^{B2} #lambda"); 
2262   regtitle33.Prepend(part);
2263   TH2F *frame33reg = new TH2F("frame33reg",regtitle33.Data(),2,0.1,50,2,0,2);
2264   frame33reg->SetXTitle("p [GeV/c]");
2265   canv33->cd(2); gPad->SetLogx();
2266   frame33reg->Draw();
2267   gr33reg->Draw("P");
2268   lin->Draw("same");
2269
2270
2271   //    --- c40 ---
2272
2273   // create the canvas
2274   TCanvas *canv40 = new TCanvas("canv40","c40",0,0,700,900); 
2275   canv40->Divide(1,2);
2276   // create the graph for cov matrix
2277   TGraphErrors *gr40 = new TGraphErrors(fitbins,p,mean40,ep,sigma40);
2278   TString title40("C(C,y)"); title40.Prepend(part);
2279   TH2F *frame40 = new TH2F("frame40",title40.Data(),2,0.1,50,2,0,1e-6);
2280   frame40->SetXTitle("p [GeV/c]");
2281   canv40->cd(1);  gPad->SetLogx();
2282   frame40->Draw();
2283   gr40->Draw("P");
2284   // Sets initial values for parameters
2285   func->SetParameters(4.e-7,4.4e-8,1.5,0.);
2286   // Fit points in range defined by function
2287   gr40->Fit("FitRegFunc","R,Q");
2288   func->GetParameters(fitpar);
2289   for(Int_t i=0; i<4; i++) fitRes(6,i)=fitpar[i];
2290   for(Int_t l=0; l<fitbins; l++) {
2291     rmean[l]  = mean40[l]/FitRegFunc(&p[l],fitpar);
2292     rsigma[l] = sigma40[l]/FitRegFunc(&p[l],fitpar);
2293   }
2294   // create the graph the regularized cov. matrix
2295   TGraphErrors *gr40reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
2296   TString regtitle40("C(C,y)/(A_meas+A_scatt/p^{B1}/cos^{B2} #lambda"); 
2297   regtitle40.Prepend(part);
2298   TH2F *frame40reg = new TH2F("frame40reg",regtitle40.Data(),2,0.1,50,2,0,2);
2299   frame40reg->SetXTitle("p [GeV/c]");
2300   canv40->cd(2); gPad->SetLogx();
2301   frame40reg->Draw();
2302   gr40reg->Draw("P");
2303   lin->Draw("same");
2304
2305
2306   //    --- c42 ---
2307
2308   // create the canvas
2309   TCanvas *canv42 = new TCanvas("canv42","c42",0,0,700,900); 
2310   canv42->Divide(1,2);
2311   // create the graph for cov matrix
2312   TGraphErrors *gr42 = new TGraphErrors(fitbins,p,mean42,ep,sigma42);
2313   TString title42("C(C, #eta)"); title42.Prepend(part);
2314   TH2F *frame42 = new TH2F("frame42",title42.Data(),2,0.1,50,2,0,2.2e-7);
2315   frame42->SetXTitle("p [GeV/c]");
2316   canv42->cd(1);  gPad->SetLogx();
2317   frame42->Draw();
2318   gr42->Draw("P");
2319   // Sets initial values for parameters
2320   func->SetParameters(3.e-8,8.2e-9,2.,1.);
2321   // Fit points in range defined by function
2322   gr42->Fit("FitRegFunc","R,Q");
2323   func->GetParameters(fitpar);
2324   for(Int_t i=0; i<4; i++) fitRes(7,i)=fitpar[i];
2325   for(Int_t l=0; l<fitbins; l++) {
2326     rmean[l]  = mean42[l]/FitRegFunc(&p[l],fitpar);
2327     rsigma[l] = sigma42[l]/FitRegFunc(&p[l],fitpar);
2328   }
2329   // create the graph the regularized cov. matrix
2330   TGraphErrors *gr42reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
2331   TString regtitle42("C(C, #eta)/(A_meas+A_scatt/p^{B1}/cos^{B2} #lambda"); 
2332   regtitle42.Prepend(part);
2333   TH2F *frame42reg = new TH2F("frame42reg",regtitle42.Data(),2,0.1,50,2,0,2);
2334   frame42reg->SetXTitle("p [GeV/c]");
2335   canv42->cd(2); gPad->SetLogx();
2336   frame42reg->Draw();
2337   gr42reg->Draw("P");
2338   lin->Draw("same");
2339
2340
2341   //    --- c44 ---
2342
2343   // create the canvas
2344   TCanvas *canv44 = new TCanvas("canv44","c44",0,0,700,900); 
2345   canv44->Divide(1,2);
2346   // create the graph for cov matrix
2347   TGraphErrors *gr44 = new TGraphErrors(fitbins,p,mean44,ep,sigma44);
2348   TString title44("C(C,C)"); title44.Prepend(part);
2349   TH2F *frame44 = new TH2F("frame44",title44.Data(),2,0.1,50,2,0,2e-9);
2350   frame44->SetXTitle("p [GeV/c]");
2351   canv44->cd(1);  gPad->SetLogx();
2352   frame44->Draw();
2353   gr44->Draw("P");
2354   // Sets initial values for parameters
2355   func->SetParameters(1.8e-10,5.8e-11,2.,3.);
2356   // Fit points in range defined by function
2357   gr44->Fit("FitRegFunc","R,Q");
2358   func->GetParameters(fitpar);
2359   for(Int_t i=0; i<4; i++) fitRes(8,i)=fitpar[i];
2360   for(Int_t l=0; l<fitbins; l++) {
2361     rmean[l]  = mean44[l]/FitRegFunc(&p[l],fitpar);
2362     rsigma[l] = sigma44[l]/FitRegFunc(&p[l],fitpar);
2363   }
2364   // create the graph the regularized cov. matrix
2365   TGraphErrors *gr44reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
2366   TString regtitle44("C(C,C)/(A_meas+A_scatt/p^{B1}/cos^{B2} #lambda"); 
2367   regtitle44.Prepend(part);
2368   TH2F *frame44reg = new TH2F("frame44reg",regtitle44.Data(),2,0.1,50,2,0,2);
2369   frame44reg->SetXTitle("p [GeV/c]");
2370   canv44->cd(2); gPad->SetLogx();
2371   frame44reg->Draw();
2372   gr44reg->Draw("P");
2373   lin->Draw("same");
2374
2375
2376
2377
2378   switch (pdg) {
2379   case 211:
2380     new(&fRegParPi) TMatrixD(fitRes);
2381     break;
2382   case 321:
2383     new(&fRegParKa) TMatrixD(fitRes);
2384     break;
2385   case 11:
2386     new(&fRegParEl) TMatrixD(fitRes);
2387     break;
2388   }
2389
2390   // write fit parameters to file
2391   WriteRegParams(outName,pdg);
2392
2393   return;
2394 }
2395 //-----------------------------------------------------------------------------
2396 Bool_t AliTPCtrackerParam::SelectedTrack(Double_t pt,Double_t eta) const {
2397 //-----------------------------------------------------------------------------
2398 // This function makes a selection according to TPC tracking efficiency
2399 //-----------------------------------------------------------------------------
2400
2401   Double_t eff=0.;
2402  
2403   eff = fEff->GetValueAt(pt,eta);
2404
2405   if(gRandom->Rndm() < eff) return kTRUE;
2406
2407   return kFALSE;
2408 }
2409 //-----------------------------------------------------------------------------
2410 void AliTPCtrackerParam::SetParticle(Int_t pdg) {
2411 //-----------------------------------------------------------------------------
2412 // This function set efficiencies, pulls, etc... for the particle
2413 // specie of the current particle
2414 //-----------------------------------------------------------------------------
2415
2416   switch (pdg) {
2417   case 211:
2418     fDBgrid = &fDBgridPi;
2419     fEff    = &fEffPi;
2420     fPulls  =  fPullsPi;
2421     fRegPar = &fRegParPi;
2422     fdEdxMean = &fdEdxMeanPi;
2423     fdEdxRMS  = &fdEdxRMSPi;
2424     break;
2425   case 321:
2426     fDBgrid = &fDBgridKa;
2427     fEff    = &fEffKa;
2428     fPulls  =  fPullsKa;
2429     fRegPar = &fRegParKa;
2430     fdEdxMean = &fdEdxMeanKa;
2431     fdEdxRMS  = &fdEdxRMSKa;
2432     break;
2433   case 2212:
2434     fDBgrid = &fDBgridPi;
2435     fEff    = &fEffPr;
2436     fPulls  =  fPullsPi;
2437     fRegPar = &fRegParPi;
2438     fdEdxMean = &fdEdxMeanPr;
2439     fdEdxRMS  = &fdEdxRMSPr;
2440     break;
2441   case 11:
2442     fDBgrid = &fDBgridEl;
2443     fEff    = &fEffEl;
2444     fPulls  =  fPullsEl;
2445     fRegPar = &fRegParEl;
2446     fdEdxMean = &fdEdxMeanEl;
2447     fdEdxRMS  = &fdEdxRMSEl;
2448     break;
2449   case 13:
2450     fDBgrid = &fDBgridPi;
2451     fEff    = &fEffMu;
2452     fPulls  =  fPullsPi;
2453     fRegPar = &fRegParPi;
2454     fdEdxMean = &fdEdxMeanPi;
2455     fdEdxRMS  = &fdEdxRMSPi;
2456     break;
2457   }
2458
2459   return;
2460 }
2461 //-----------------------------------------------------------------------------
2462 void AliTPCtrackerParam::SmearTrack(Double_t *xx,Double_t *xxsm,TMatrixD cov)
2463   const {
2464 //-----------------------------------------------------------------------------
2465 // This function smears track parameters according to streched cov. matrix
2466 //-----------------------------------------------------------------------------
2467   AliGausCorr *corgen = new AliGausCorr(cov,5);
2468   TArrayD corr(5);
2469   corgen->GetGaussN(corr);
2470   delete corgen;
2471   corgen = 0;
2472
2473   for(Int_t l=0;l<5;l++) {
2474     xxsm[l] = xx[l]+corr[l];
2475   }
2476
2477   return;
2478 }
2479 //-----------------------------------------------------------------------------
2480 Int_t AliTPCtrackerParam::WritedEdx(const Char_t *outName,Int_t pdg) {
2481 //-----------------------------------------------------------------------------
2482 // This function writes the dEdx parameters to the DB
2483 //-----------------------------------------------------------------------------
2484
2485   Option_t *opt;
2486   Char_t *dirName="Pions";
2487   Char_t *meanName="dEdxMeanPi";
2488   Char_t *RMSName="dEdxRMSPi";
2489
2490   SetParticle(pdg);
2491
2492   if(gSystem->AccessPathName(outName,kFileExists))  { opt="recreate"; 
2493   } else { opt="update"; }
2494
2495   switch (pdg) {
2496   case 211:
2497     dirName="Pions";
2498     meanName="dEdxMeanPi";
2499     RMSName="dEdxRMSPi";
2500     break;
2501   case 321:
2502     dirName="Kaons";
2503     meanName="dEdxMeanKa";
2504     RMSName="dEdxRMSKa";
2505     break;
2506   case 2212:
2507     dirName="Protons";
2508     meanName="dEdxMeanPr";
2509     RMSName="dEdxRMSPr";
2510     break;
2511   case 11:
2512     dirName="Electrons";
2513     meanName="dEdxMeanEl";
2514     RMSName="dEdxRMSEl";
2515     break;
2516   }
2517
2518   TFile *outFile = new TFile(outName,opt);
2519   if(!gDirectory->cd("/dEdx")) {
2520     outFile->mkdir("dEdx");
2521     gDirectory->cd("/dEdx"); 
2522   }
2523   TDirectory *dir2 = gDirectory->mkdir(dirName);
2524   dir2->cd();
2525   fdEdxMean->SetName(meanName); fdEdxMean->Write();
2526   fdEdxRMS->SetName(RMSName);  fdEdxRMS->Write();
2527
2528   outFile->Close();
2529   delete outFile;
2530
2531
2532   return 1;
2533 }
2534 //-----------------------------------------------------------------------------
2535 Int_t AliTPCtrackerParam::WriteEffs(const Char_t *outName) {
2536 //-----------------------------------------------------------------------------
2537 // This function writes the TPC efficiencies to the DB
2538 //-----------------------------------------------------------------------------
2539
2540
2541   Option_t *opt;
2542
2543   if(gSystem->AccessPathName(outName,kFileExists))  { opt="recreate"; 
2544   } else { opt="update"; }
2545
2546   TFile *outFile = new TFile(outName,opt);
2547
2548   outFile->mkdir("Efficiencies");
2549   gDirectory->cd("/Efficiencies");
2550   gDirectory->mkdir("Pions");
2551   gDirectory->mkdir("Kaons");
2552   gDirectory->mkdir("Protons");
2553   gDirectory->mkdir("Electrons");
2554   gDirectory->mkdir("Muons");
2555
2556   gDirectory->cd("/Efficiencies/Pions");
2557   fEffPi.SetName("EffPi");
2558   fEffPi.Write();
2559   gDirectory->cd("/Efficiencies/Kaons");
2560   fEffKa.SetName("EffKa");
2561   fEffKa.Write();
2562   gDirectory->cd("/Efficiencies/Protons");
2563   fEffPr.SetName("EffPr");
2564   fEffPr.Write();
2565   gDirectory->cd("/Efficiencies/Electrons");
2566   fEffEl.SetName("EffEl");
2567   fEffEl.Write();
2568   gDirectory->cd("/Efficiencies/Muons");
2569   fEffMu.SetName("EffMu");
2570   fEffMu.Write();
2571
2572   outFile->Close();
2573
2574   delete outFile;
2575
2576   return 1;
2577 }
2578 //-----------------------------------------------------------------------------
2579 Int_t AliTPCtrackerParam::WritePulls(const Char_t *outName) {
2580 //-----------------------------------------------------------------------------
2581 // This function writes the pulls to the DB
2582 //-----------------------------------------------------------------------------
2583
2584   Option_t *opt;
2585
2586   if(gSystem->AccessPathName(outName,kFileExists))  { opt="recreate"; 
2587   } else { opt="update"; }
2588
2589   TFile *outFile = new TFile(outName,opt);
2590
2591   outFile->mkdir("Pulls");
2592   gDirectory->cd("/Pulls");
2593   gDirectory->mkdir("Pions");
2594   gDirectory->mkdir("Kaons");
2595   gDirectory->mkdir("Electrons");
2596
2597   for(Int_t i=0;i<5;i++) {
2598     TString pi("PullsPi_"); pi+=i;
2599     TString ka("PullsKa_"); ka+=i;
2600     TString el("PullsEl_"); el+=i;
2601     fPullsPi[i].SetName(pi.Data());
2602     fPullsKa[i].SetName(ka.Data());
2603     fPullsEl[i].SetName(el.Data());
2604     gDirectory->cd("/Pulls/Pions");
2605     fPullsPi[i].Write();
2606     gDirectory->cd("/Pulls/Kaons");
2607     fPullsKa[i].Write();
2608     gDirectory->cd("/Pulls/Electrons");
2609     fPullsEl[i].Write();
2610   }
2611   outFile->Close();
2612   delete outFile;
2613
2614   return 1;
2615 }
2616 //-----------------------------------------------------------------------------
2617 Int_t AliTPCtrackerParam::WriteRegParams(const Char_t *outName,Int_t pdg) {
2618 //-----------------------------------------------------------------------------
2619 // This function writes the regularization parameters to the DB
2620 //-----------------------------------------------------------------------------
2621
2622   Option_t *opt;
2623   Char_t *dirName="Pions";
2624   Char_t *keyName="RegPions";
2625
2626   SetParticle(pdg);
2627
2628   if(gSystem->AccessPathName(outName,kFileExists))  { opt="recreate"; 
2629   } else { opt="update"; }
2630
2631   switch (pdg) {
2632   case 211:
2633     dirName="Pions";
2634     keyName="RegPions";
2635     break;
2636   case 321:
2637     dirName="Kaons";
2638     keyName="RegKaons";
2639     break;
2640   case 11:
2641     dirName="Electrons";
2642     keyName="RegElectrons";
2643     break;
2644   }
2645
2646   TFile *outFile = new TFile(outName,opt);
2647   if(!gDirectory->cd("/RegParams")) {
2648     outFile->mkdir("RegParams");
2649     gDirectory->cd("/RegParams"); 
2650   }
2651   TDirectory *dir2 = gDirectory->mkdir(dirName);
2652   dir2->cd();
2653   fRegPar->Write(keyName);
2654
2655   outFile->Close();
2656   delete outFile;
2657
2658
2659   return 1;
2660 }
2661
2662
2663
2664
2665
2666
2667
2668
2669