Changes needed on Sun with Root v4-03-04
[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 /* $Id$ */
17
18 /**************************************************************************
19  *                                                                        *
20  * This class builds AliTPCtrack objects from generated tracks to feed    *
21  * ITS tracking (V2). The AliTPCtrack is built from its first hit in      *
22  * the TPC. The track is assigned a Kalman-like covariance matrix         *
23  * depending on its pT and pseudorapidity and track parameters are        *
24  * smeared according to this covariance matrix.                           *
25  * Output file contains sorted tracks, ready for matching with ITS.       *
26  *                                                                        *
27  * For details:                                                           *
28  * Alice Internal Note 2003-011                                           *
29  *                                                                        *
30  * Test macro is: AliBarrelRec_TPCparam.C                                 *   
31  *                                                                        *
32  * 2002/10/01: Introduction of the database for pp collisions (B=0.4 T)   *
33  * - Everything done separately for pions, kaons, protons, electrons and  *
34  *   muons.                                                               *
35  * - Now (only for pp) the tracks are built from the AliTrackReferences   *
36  *   which contain the position and momentum of all tracks at R = 83 cm;  *
37  *   This eliminates the loss of tracks due the dead zone of the TPC      *
38  *   where the 1st hit is not created.                                    *
39  * - In AliBarrelRec_TPCparam.C there many possible ways of determining   *
40  *   the z coordinate of the primary vertex in pp events (from pixels,    *
41  *   from ITS tracks, smearing according to resolution given by tracks.   *
42  *                                                                        *
43  * 2002/04/28: Major upgrade of the class                                 *
44  * - Covariance matrices and pulls are now separeted for pions, kaons     *
45  *   and electrons.                                                       *
46  * - A parameterization of dE/dx in the TPC has been included and it is   *
47  *   used to assign a mass to each track according to a rough dE/dx PID.  *
48  * - All the "numbers" have been moved to the file with the DataBase,     *
49  *   they are read as objects of the class AliTPCkineGrid, and assigned   *
50  *   to data memebers of the class AliTPCtrackerParam.                    *
51  * - All the code necessary to create a BataBase has been included in     *
52  *   class (see the macro AliTPCtrackingParamDB.C for the details).       *
53  *                                                                        *
54  *  Origin: Andrea Dainese, Padova - e-mail: andrea.dainese@pd.infn.it    *
55  *                                                                        *
56  **************************************************************************/
57 //  *
58 // This is a dummy comment
59 //
60 //
61 // *
62 //------- Root headers --------
63 #include <TChain.h>
64 #include <TF1.h>
65 #include <TGraph.h>
66 #include <TGraphErrors.h>
67 #include <TLegend.h>
68 #include <TLine.h>
69 #include <TNtuple.h>
70 #include <TSystem.h>
71 #include <TH2F.h>
72 //------ AliRoot headers ------
73 #include "alles.h"
74 #include "AliGausCorr.h"
75 #include "AliKalmanTrack.h"
76 #include "AliMagF.h"
77 #include "AliMagFCM.h"
78 #include "AliRunLoader.h"
79 #include "AliTPCLoader.h"
80 #include "AliTPCkineGrid.h"
81 #include "AliTPCtrack.h"
82 #include "AliTrackReference.h"
83 #include "AliTPCtrackerParam.h"
84 #include "AliMC.h"
85 //-----------------------------
86
87 Double_t RegFunc(Double_t *x,Double_t *par) {
88 // This is the function used to regularize the covariance matrix
89   Double_t value = par[0]+par[1]/TMath::Power(x[0],par[2]);
90   return value;
91 }
92
93 // structure for DB building
94 typedef struct {
95   Int_t    pdg;
96   Int_t    bin;
97   Double_t r;
98   Double_t p;
99   Double_t pt;
100   Double_t cosl;
101   Double_t eta;
102   Double_t dpt;
103   Double_t dP0,dP1,dP2,dP3,dP4;
104   Double_t c00,c10,c11,c20,c21,c22,c30,c31,c32,c33,c40,c41,c42,c43,c44;
105   Double_t dEdx;
106 } COMPTRACK;
107 // cov matrix structure 
108 typedef struct {
109   Double_t c00,c10,c11,c20,c21,c22,c30,c31,c32,c33,c40,c41,c42,c43,c44;
110 } COVMATRIX;
111
112 ClassImp(AliTPCtrackerParam)
113
114 //-----------------------------------------------------------------------------
115 AliTPCtrackerParam::AliTPCtrackerParam(Int_t kcoll, Double_t kBz,
116                                        Int_t kn, const char* evfoldname):
117   fEvFolderName(evfoldname) {
118 //-----------------------------------------------------------------------------
119 // This is the class conctructor 
120 //-----------------------------------------------------------------------------
121
122   fNevents = kn;         // events to be processed
123   fBz = kBz;             // value of the z component of L3 field (Tesla)
124   fColl = kcoll;         // collision code (0: PbPb6000; 1: pp)
125   fSelAndSmear = kTRUE; // by default selection and smearing are done
126
127   if(fBz!=0.4) {
128     cerr<<"AliTPCtrackerParam::AliTPCtrackerParam:  Invalid field!\n";
129     cerr<<"      Available:  0.4"<<endl;
130   }
131   if(fColl!=0 && fColl!=1) { 
132     cerr<<"AliTPCtrackerParam::AliTPCtrackerParam:  Invalid collision!\n";
133     cerr<<"      Available:  0   ->   PbPb6000"<<endl;
134     cerr<<"                  1   ->   pp"<<endl; 
135   }
136
137   fDBfileName = gSystem->Getenv("ALICE_ROOT");  
138   fDBfileName.Append("/TPC/CovMatrixDB_");
139   //fDBfileName = "CovMatrixDB_";
140   if(fColl==0) fDBfileName.Append("PbPb6000");
141   if(fColl==1) fDBfileName.Append("pp");
142   if(fBz==0.4) fDBfileName.Append("_B0.4T.root");
143 }
144 //-----------------------------------------------------------------------------
145 AliTPCtrackerParam::~AliTPCtrackerParam() {}
146 //____________________________________________________________________________
147 AliTPCtrackerParam::AliTPCtrackerParam( const AliTPCtrackerParam& p):TObject(p)
148 {
149   // dummy copy constructor
150 }
151 //----------------------------------------------------------------------------
152 AliTPCtrackerParam::AliTPCseedGeant::AliTPCseedGeant(
153                     Double_t x,Double_t y,Double_t z,
154                     Double_t px,Double_t py,Double_t pz,
155                     Int_t lab) {
156 //----------------------------------------------------------------------------
157 // Constructor of the geant seeds
158 //----------------------------------------------------------------------------
159       fXg = x;
160       fYg = y;
161       fZg = z;
162       fPx = px;
163       fPy = py;
164       fPz = pz;
165       fLabel = lab;
166       Double_t a = TMath::ATan2(y,x)*180./TMath::Pi();
167       if(a<0) a += 360.;
168       fSector = (Int_t)(a/20.);
169       fAlpha = 10.+20.*fSector;
170       fAlpha /= 180.;
171       fAlpha *= TMath::Pi();
172 }
173 //-----------------------------------------------------------------------------
174 Int_t AliTPCtrackerParam::BuildTPCtracks(const TFile *inp, TFile *out) {
175 //-----------------------------------------------------------------------------
176 // This function creates the TPC parameterized tracks
177 //-----------------------------------------------------------------------------
178
179   Error("BuildTPCtracks","in and out parameters ignored. new io");
180   
181 /********************************************/
182   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEvFolderName);
183   if (rl == 0x0)
184    {
185      Error("BuildTPCtracks","Can not get Run Loader from event folder named %s.",
186            fEvFolderName.Data());
187      return 2;
188    }
189   AliLoader* tpcloader = rl->GetLoader("TPCLoader");
190   if (tpcloader == 0x0)
191    {
192      Error("BuildTPCtracks","Can not get TPC Loader from Run Loader.");
193      return 3;
194    }
195   
196 /********************************************/
197
198   TFile *fileDB=0;
199   TTree *covTreePi[50];
200   TTree *covTreeKa[50];
201   TTree *covTreePr[50];
202   TTree *covTreeEl[50];
203   TTree *covTreeMu[50];
204
205   if(fSelAndSmear) {
206     cerr<<"+++\n+++ Reading DataBase from:\n+++ "<<
207       fDBfileName.Data()<<"\n+++\n"; 
208     // Read paramters from DB file
209     if(!ReadAllData(fDBfileName.Data())) {
210       cerr<<"AliTPCtrackerParam::BuildTPCtracks: \
211              Could not read data from DB\n\n"; return 1; 
212     }
213
214     // Read the trees with regularized cov. matrices from DB
215     TString str;
216     fileDB = TFile::Open(fDBfileName.Data());
217     Int_t nBinsPi = fDBgridPi.GetTotBins();
218     for(Int_t l=0; l<nBinsPi; l++) {
219       str = "/CovMatrices/Pions/CovTreePi_bin";
220       str+=l;
221       covTreePi[l] = (TTree*)fileDB->Get(str.Data());
222     }
223     Int_t nBinsKa = fDBgridKa.GetTotBins();
224     for(Int_t l=0; l<nBinsKa; l++) {
225       str = "/CovMatrices/Kaons/CovTreeKa_bin";
226       str+=l;
227       covTreeKa[l] = (TTree*)fileDB->Get(str.Data());
228     }
229     Int_t nBinsPr = fDBgridPr.GetTotBins();
230     for(Int_t l=0; l<nBinsPr; l++) {
231       if(fColl==0) {      
232         str = "/CovMatrices/Pions/CovTreePi_bin";
233       } else {
234         str = "/CovMatrices/Protons/CovTreePr_bin";
235       }
236       str+=l;
237       covTreePr[l] = (TTree*)fileDB->Get(str.Data());
238     }
239     Int_t nBinsEl = fDBgridEl.GetTotBins();
240     for(Int_t l=0; l<nBinsEl; l++) {
241       str = "/CovMatrices/Electrons/CovTreeEl_bin";
242       str+=l;
243       covTreeEl[l] = (TTree*)fileDB->Get(str.Data());
244     }
245     Int_t nBinsMu = fDBgridMu.GetTotBins();
246     for(Int_t l=0; l<nBinsMu; l++) {
247       if(fColl==0) {      
248         str = "/CovMatrices/Pions/CovTreePi_bin";
249       } else {
250         str = "/CovMatrices/Muons/CovTreeMu_bin";
251       }
252       str+=l;
253       covTreeMu[l] = (TTree*)fileDB->Get(str.Data());
254     }
255
256   } else cerr<<"\n ! Creating ALL TRUE tracks at TPC 1st hit !\n\n";
257
258   TFile *infile=(TFile*)inp;
259   infile->cd();
260
261   // Get gAlice object from file
262   
263   rl->LoadgAlice();
264   rl->LoadHeader();
265   tpcloader->LoadHits("read");
266   
267   if(!(gAlice=rl->GetAliRun())) {
268     cerr<<"Can not get gAlice from Run Loader !\n";
269     return 1;
270   }
271
272   // Check if value in the galice file is equal to selected one (fBz)
273   AliMagF *fiel = (AliMagF*)gAlice->Field();
274   Double_t fieval=(Double_t)fiel->SolenoidField()/10.;
275   printf("Magnetic field is %6.2f Tesla\n",fieval);
276   if(fBz!=fieval) {
277     cerr<<"AliTPCtrackerParam::BuildTPCtracks:  Invalid field!"<<endl;
278     cerr<<"Field selected is: "<<fBz<<" T\n";
279     cerr<<"Field found on file is: "<<fieval<<" T\n";
280     return 0;
281   }
282
283   // Get TPC detector 
284   AliTPC *tpc=(AliTPC*)gAlice->GetDetector("TPC");
285   Int_t ver = tpc->IsVersion(); 
286   cerr<<"+++ TPC version "<<ver<<" has been found !\n";
287
288   rl->CdGAFile();
289   AliTPCParam *digp=(AliTPCParam*)infile->Get("75x40_100x60");
290   if(digp){
291     delete digp;
292     digp = new AliTPCParamSR();
293   }
294   else digp=(AliTPCParam*)infile->Get("75x40_100x60_150x60");
295   
296   if(!digp) { cerr<<"TPC parameters have not been found !\n"; return 1; }
297   tpc->SetParam(digp);
298
299   // Set the conversion constant between curvature and Pt
300   AliKalmanTrack::SetConvConst(100/0.299792458/fBz);
301
302   TParticle       *part=0;
303   AliTPCseedGeant *seed=0;
304   AliTPCtrack     *tpctrack=0;
305   Double_t     sPt,sEta;
306   Int_t        cl=0,bin,label,pdg,charge;
307   Int_t        tracks;
308   Int_t        nParticles,nSeeds,arrentr;
309   Char_t       tname[100];
310   //Int_t nSel=0,nAcc=0;
311
312
313   // loop over first n events in file
314   for(Int_t evt=0; evt<fNevents; evt++){
315     cerr<<"+++\n+++ Processing event "<<evt<<"\n+++\n";
316     tracks=0;
317
318     // tree for TPC tracks
319     sprintf(tname,"TreeT_TPC_%d",evt);
320     TTree *tracktree = new TTree(tname,"Tree with TPC tracks");
321     tracktree->Branch("tracks","AliTPCtrack",&tpctrack,20000,0);
322
323     // array for TPC tracks
324     TObjArray tArray(20000);
325
326     // array for TPC seeds with geant info
327     TObjArray sArray(20000);
328
329     // get the particles stack
330     nParticles = (Int_t)gAlice->GetEvent(evt);
331     
332     Bool_t   *done     = new Bool_t[nParticles];
333     Int_t    *pdgCodes = new Int_t[nParticles];
334     Double_t *ptkine   = new Double_t[nParticles];
335     Double_t *pzkine   = new Double_t[nParticles];
336
337     // loop on particles and store pdg codes
338     for(Int_t l=0; l<nParticles; l++) {
339       part        = (TParticle*)gAlice->GetMCApp()->Particle(l);
340       pdgCodes[l] = part->GetPdgCode();
341       ptkine[l]   = part->Pt();
342       pzkine[l]   = part->Pz();
343       done[l]     = kFALSE;
344     }
345     cerr<<"+++\n+++ Number of particles in event "<<evt<<":  "<<nParticles<<
346       "\n+++\n";
347
348     cerr<<"\n ********** MAKING SEEDS *****************"<<endl<<endl;
349
350
351     // Create the seeds for the TPC tracks at the inner radius of TPC
352     if(fColl==0) {
353       // Get TreeH with hits
354       TTree *th = tpcloader->TreeH(); 
355       MakeSeedsFromHits(tpc,th,sArray);
356     } else {
357       // Get TreeTR with track references
358       TTree *ttr = rl->TreeTR();
359       MakeSeedsFromRefs(ttr,sArray);
360     }
361
362
363     nSeeds = sArray.GetEntries();
364     cerr<<"\n\n+++\n+++ Number of seeds: "<<nSeeds<<"\n+++\n";
365
366
367     cerr<<"\n ********** BUILDING TRACKS **************"<<endl<<endl;
368
369     // loop over entries in sArray
370     for(Int_t l=0; l<nSeeds; l++) {
371       //if(l%1000==0) cerr<<"  --- Processing seed "
372       //                <<l<<" of "<<nSeeds<<" ---\r";
373
374       seed = (AliTPCseedGeant*)sArray.At(l);
375
376       // this is TEMPORARY: only for reconstruction of pp production for charm
377       if(fColl==1) cl = CheckLabel(seed,nParticles,ptkine,pzkine);
378       if(cl) continue;
379
380       // Get track label
381       label = seed->GetLabel();
382
383       // check if this track has already been processed
384       if(done[label]) continue;
385       // PDG code & electric charge
386       pdg = pdgCodes[label];
387       if(pdg>200 || pdg==-11 || pdg==-13) { charge=1; }
388       else if(pdg<-200 || pdg==11 || pdg==13) { charge=-1; }
389       else continue;
390       pdg = TMath::Abs(pdg);
391       if(pdg>3000) pdg=211;
392  
393       if(fSelAndSmear) SetParticle(pdg);
394
395
396       sPt  = seed->GetPt();
397       sEta = seed->GetEta();
398       
399       // Apply selection according to TPC efficiency
400       //if(TMath::Abs(pdg)==211) nAcc++;
401       if(fSelAndSmear && !SelectedTrack(sPt,sEta)) continue; 
402       //if(TMath::Abs(pdg)==211) nSel++;
403
404       // create AliTPCtrack object
405       BuildTrack(seed,charge);
406         
407       if(fSelAndSmear) {
408         bin = fDBgrid->GetBin(sPt,sEta);
409         switch (pdg) {
410         case 211:
411           fCovTree = covTreePi[bin];
412           break;
413         case 321:
414           fCovTree = covTreeKa[bin];
415           break;
416         case 2212:
417           fCovTree = covTreePr[bin];
418           break;
419         case 11:
420           fCovTree = covTreeEl[bin];
421           break;
422         case 13:
423           fCovTree = covTreeMu[bin];
424           break;
425         }
426         // deal with covariance matrix and smearing of parameters
427         CookTrack(sPt,sEta);
428
429         // assign the track a dE/dx and make a rough PID
430         CookdEdx(sPt,sEta);
431       }
432
433       // put track in array
434       AliTPCtrack *iotrack = new AliTPCtrack(fTrack);
435       iotrack->SetLabel(label);
436       tArray.AddLast(iotrack);
437       // Mark track as "done" and register the pdg code
438       done[label] = kTRUE; 
439       tracks++;
440  
441     } // loop over entries in sArray
442
443
444     // sort array with TPC tracks (decreasing pT)
445     tArray.Sort();
446
447     arrentr = tArray.GetEntriesFast();
448     for(Int_t l=0; l<arrentr; l++) {
449       tpctrack=(AliTPCtrack*)tArray.UncheckedAt(l);
450       tracktree->Fill();
451     }
452
453
454     // write the tree with tracks in the output file
455     out->cd();
456     tracktree->Write(); 
457     
458     delete tracktree;
459     delete [] done;
460     delete [] pdgCodes;
461     delete [] ptkine;
462     delete [] pzkine;    
463
464     printf("\n\n+++\n+++ Number of TPC tracks: %d\n+++\n",tracks);
465     //cerr<<"Average Eff: "<<(Float_t)nSel/nAcc<<endl;
466
467     sArray.Delete();
468     tArray.Delete();
469   
470     infile->cd();
471   } // loop on events
472
473   if(fileDB) fileDB->Close();
474
475   return 0;
476 }
477 //-----------------------------------------------------------------------------
478 void AliTPCtrackerParam::AnalyzedEdx(const Char_t *outName,Int_t pdg) {
479 //-----------------------------------------------------------------------------
480 // This function computes the dE/dx for pions, kaons, protons and electrons,
481 // in the [pT,eta] bins.
482 // Input file is CovMatrix_AllEvts.root.  
483 //-----------------------------------------------------------------------------
484
485   gStyle->SetOptStat(0);
486   gStyle->SetOptFit(10001);
487
488   const char *part="PIONS";
489   Double_t ymax=500.;
490
491   /*
492   // create a chain with compared tracks
493   TChain *cmptrkchain = new ("cmptrktree");
494   cmptrkchain.Add("CovMatrix_AllEvts.root");
495   //cmptrkchain.Add("CovMatrix_AllEvts_1.root");
496   //cmptrkchain.Add("CovMatrix_AllEvts_2.root");
497   //cmptrkchain.Add("CovMatrix_AllEvts_3.root");
498   */
499
500
501   TFile *infile = TFile::Open("CovMatrix_AllEvts.root");
502   TTree *cmptrktree = (TTree*)infile->Get("cmptrktree");
503
504   COMPTRACK cmptrk; 
505   cmptrktree->SetBranchAddress("comptracks",&cmptrk);
506   Int_t entries = (Int_t)cmptrktree->GetEntries(); 
507   cerr<<" Number of entries: "<<entries<<endl;
508
509   InitializeKineGrid("DB");
510   InitializeKineGrid("dEdx");
511
512   switch(pdg) {
513   case 211:
514     part = "PIONS";
515     ymax = 100.;
516     break;
517   case 321:
518     part = "KAONS";
519     ymax = 300.;
520     break;
521   case 2212:
522     part = "PROTONS";
523     ymax = 500.;
524     break;
525   case 11:
526     part = "ELECTRONS";
527     ymax = 100.;
528     break;
529   case 13:
530     part = "MUONS";
531     ymax = 100.;
532     break;
533   }
534
535   SetParticle(pdg);
536
537   const Int_t knTotBins = fDBgrid->GetTotBins(); 
538
539   cerr<<" Fit bins: "<<knTotBins<<endl;
540
541   Int_t bin=0;
542   Int_t        *n = new Int_t[knTotBins];
543   Double_t     *p = new Double_t[knTotBins];
544   Double_t    *ep = new Double_t[knTotBins];
545   Double_t  *mean = new Double_t[knTotBins];
546   Double_t *sigma = new Double_t[knTotBins];
547
548   for(Int_t l=0; l<knTotBins; l++) {
549     n[l] = 1; // set to 1 to avoid divisions by 0
550     p[l] = mean[l] = sigma[l] = ep[l] = 0.; 
551   }
552
553   // loop on chain entries for the mean
554   for(Int_t l=0; l<entries; l++) {
555     cmptrktree->GetEvent(l);
556     if(TMath::Abs(cmptrk.pdg)!=pdg) continue;
557     bin = fDBgrid->GetBin(cmptrk.pt,cmptrk.eta);
558     p[bin] += cmptrk.p;
559     mean[bin] += cmptrk.dEdx;
560     n[bin]++;
561   } // loop on chain entries
562
563   for(Int_t l=0; l<knTotBins; l++) {
564     p[l] /= n[l];
565     mean[l] /= n[l];
566     n[l] = 1; // set to 1 to avoid divisions by 0
567   }
568
569   // loop on chain entries for the sigma
570   for(Int_t l=0; l<entries; l++) {
571     cmptrktree->GetEvent(l);
572     if(TMath::Abs(cmptrk.pdg)!=pdg) continue;
573     bin = fDBgrid->GetBin(cmptrk.pt,cmptrk.eta);
574     if(cmptrk.p<1. && TMath::Abs(cmptrk.p-p[bin])>0.025) continue;
575     n[bin]++;
576     sigma[bin] += (cmptrk.dEdx-mean[bin])*(cmptrk.dEdx-mean[bin]);
577   } // loop on chain entries
578   
579   for(Int_t l=0; l<knTotBins; l++) {
580     sigma[l] = TMath::Sqrt(sigma[l]/n[l]);
581   }
582
583   // create the canvas
584   TCanvas *canv = new TCanvas("canv","dEdx",0,0,900,700); 
585
586   // create the graph for dEdx vs p
587   TGraphErrors *gr = new TGraphErrors(knTotBins,p,mean,ep,sigma);
588   TString title("  : dE/dx vs momentum"); title.Prepend(part);
589   TH2F *frame = new TH2F("frame1",title.Data(),2,0.1,50,2,0,ymax);
590   frame->SetXTitle("p [GeV/c]");
591   frame->SetYTitle("dE/dx [a.u.]");
592   canv->SetLogx();
593   frame->Draw();
594   gr->Draw("P");
595
596   switch(pdg) {
597   case 211:
598     for(Int_t i=0; i<knTotBins; i++) {
599       fdEdxMeanPi.SetParam(i,mean[i]);
600       fdEdxRMSPi.SetParam(i,sigma[i]);
601     }    
602     break;
603   case 321:
604     for(Int_t i=0; i<knTotBins; i++) {
605       fdEdxMeanKa.SetParam(i,mean[i]);
606       fdEdxRMSKa.SetParam(i,sigma[i]);
607     }    
608     break;
609   case 2212:
610     for(Int_t i=0; i<knTotBins; i++) {
611       fdEdxMeanPr.SetParam(i,mean[i]);
612       fdEdxRMSPr.SetParam(i,sigma[i]);
613     }    
614     break;
615   case 11:
616     for(Int_t i=0; i<knTotBins; i++) {
617       fdEdxMeanEl.SetParam(i,mean[i]);
618       fdEdxRMSEl.SetParam(i,sigma[i]);
619     }    
620     break;
621   case 13:
622     for(Int_t i=0; i<knTotBins; i++) {
623       fdEdxMeanMu.SetParam(i,mean[i]);
624       fdEdxRMSMu.SetParam(i,sigma[i]);
625     }    
626     break;
627   }
628
629   // write results to file
630   WritedEdx(outName,pdg);
631
632   delete [] n;
633   delete [] p;
634   delete [] ep;
635   delete [] mean;
636   delete [] sigma;
637   
638   return;
639 }
640 //-----------------------------------------------------------------------------
641 void AliTPCtrackerParam::AnalyzePulls(const Char_t *outName) {
642 //-----------------------------------------------------------------------------
643 // This function computes the pulls for pions, kaons and electrons,
644 // in the [pT,eta] bins.
645 // Input file is CovMatrix_AllEvts.root.
646 // Output file is pulls.root.  
647 //-----------------------------------------------------------------------------
648
649   /*
650   // create a chain with compared tracks
651   TChain *cmptrkchain = new ("cmptrktree");
652   cmptrkchain.Add("CovMatrix_AllEvts.root");
653   //cmptrkchain.Add("CovMatrix_AllEvts_1.root");
654   //cmptrkchain.Add("CovMatrix_AllEvts_2.root");
655   //cmptrkchain.Add("CovMatrix_AllEvts_3.root");
656   */
657
658
659   TFile *infile = TFile::Open("CovMatrix_AllEvts.root");
660   TTree *cmptrktree = (TTree*)infile->Get("cmptrktree");
661
662   COMPTRACK cmptrk; 
663   cmptrktree->SetBranchAddress("comptracks",&cmptrk);
664   Int_t entries = (Int_t)cmptrktree->GetEntries(); 
665   cerr<<" Number of entries: "<<entries<<endl;
666
667   Int_t thisPdg=0;
668   Char_t hname[100], htitle[100];
669   Int_t nTotBins,bin;
670
671   AliTPCkineGrid  pulls[5];
672   TH1F *hDum = new TH1F("name","title",100,-7.,7.);
673   TF1 *g = new TF1("g","gaus");
674
675   InitializeKineGrid("pulls");
676   InitializeKineGrid("DB");
677
678
679
680   // loop on the particles Pi,Ka,Pr,El,Mu
681   for(Int_t part=0; part<5; part++) {
682
683     switch (part) {
684     case 0: // pions
685       thisPdg=211; 
686       cerr<<" Processing pions ...\n";
687       break;   
688     case 1: // kaons
689       thisPdg=321; 
690       cerr<<" Processing kaons ...\n";
691       break;
692     case 2: // protons
693       thisPdg=2212; 
694       cerr<<" Processing protons ...\n";
695       break;
696     case 3: // electrons
697       thisPdg=11; 
698       cerr<<" Processing electrons ...\n";
699       break;
700     case 4: // muons
701       thisPdg=13; 
702       cerr<<" Processing muons ...\n";
703       break;
704     }
705
706     SetParticle(thisPdg);
707
708     for(Int_t i=0;i<5;i++) {
709       pulls[i].~AliTPCkineGrid(); 
710       new(&pulls[i]) AliTPCkineGrid(*(fPulls+i));
711     }
712     nTotBins = fDBgrid->GetTotBins();
713     cerr<<"nTotBins = "<<nTotBins<<endl; 
714
715     // create histograms for the all the bins
716     TH1F *hPulls0=NULL;
717     TH1F *hPulls1=NULL;
718     TH1F *hPulls2=NULL;
719     TH1F *hPulls3=NULL;
720     TH1F *hPulls4=NULL;
721
722     hPulls0 = new TH1F[nTotBins]; 
723     hPulls1 = new TH1F[nTotBins]; 
724     hPulls2 = new TH1F[nTotBins]; 
725     hPulls3 = new TH1F[nTotBins]; 
726     hPulls4 = new TH1F[nTotBins]; 
727
728
729     for(Int_t i=0; i<nTotBins; i++) {
730       sprintf(hname,"hPulls0%d",i);
731       sprintf(htitle,"P0 pulls for bin %d",i);
732       hDum->SetName(hname); hDum->SetTitle(htitle);
733       hPulls0[i] = *hDum;
734       sprintf(hname,"hPulls1%d",i);
735       sprintf(htitle,"P1 pulls for bin %d",i);
736       hDum->SetName(hname); hDum->SetTitle(htitle);
737       hPulls1[i] = *hDum;
738       sprintf(hname,"hPulls2%d",i);
739       sprintf(htitle,"P2 pulls for bin %d",i);
740       hDum->SetName(hname); hDum->SetTitle(htitle);
741       hPulls2[i] = *hDum;
742       sprintf(hname,"hPulls3%d",i);
743       sprintf(htitle,"P3 pulls for bin %d",i);
744       hDum->SetName(hname); hDum->SetTitle(htitle);
745       hPulls3[i] = *hDum;
746       sprintf(hname,"hPulls4%d",i);
747       sprintf(htitle,"P4 pulls for bin %d",i);
748       hDum->SetName(hname); hDum->SetTitle(htitle);
749       hPulls4[i] = *hDum;
750     }
751
752     // loop on chain entries 
753     for(Int_t i=0; i<entries; i++) {
754       cmptrktree->GetEvent(i);
755       if(TMath::Abs(cmptrk.pdg)!=thisPdg) continue;
756       // fill histograms with the pulls
757       bin = fDBgrid->GetBin(cmptrk.pt,cmptrk.eta);
758       //cerr<<" pt "<<cmptrk.pt<<"   eta "<<cmptrk.eta<<"   bin "<<bin<<endl; 
759       hPulls0[bin].Fill(cmptrk.dP0/TMath::Sqrt(cmptrk.c00));
760       hPulls1[bin].Fill(cmptrk.dP1/TMath::Sqrt(cmptrk.c11));
761       hPulls2[bin].Fill(cmptrk.dP2/TMath::Sqrt(cmptrk.c22));
762       hPulls3[bin].Fill(cmptrk.dP3/TMath::Sqrt(cmptrk.c33));
763       hPulls4[bin].Fill(cmptrk.dP4/TMath::Sqrt(cmptrk.c44));
764     } // loop on chain entries
765
766     // compute the sigma of the distributions
767     for(Int_t i=0; i<nTotBins; i++) {
768       if(hPulls0[i].GetEntries()>10) {
769         g->SetRange(-3.*hPulls0[i].GetRMS(),3.*hPulls0[i].GetRMS());
770         hPulls0[i].Fit("g","R,Q,N");
771         pulls[0].SetParam(i,g->GetParameter(2));
772       } else pulls[0].SetParam(i,-1.);
773       if(hPulls1[i].GetEntries()>10) {
774         g->SetRange(-3.*hPulls1[i].GetRMS(),3.*hPulls1[i].GetRMS());
775         hPulls1[i].Fit("g","R,Q,N");
776         pulls[1].SetParam(i,g->GetParameter(2));
777       } else pulls[1].SetParam(i,-1.);
778       if(hPulls2[i].GetEntries()>10) {
779         g->SetRange(-3.*hPulls2[i].GetRMS(),3.*hPulls2[i].GetRMS());
780         hPulls2[i].Fit("g","R,Q,N");
781         pulls[2].SetParam(i,g->GetParameter(2));
782       } else pulls[2].SetParam(i,-1.);
783       if(hPulls3[i].GetEntries()>10) {
784         g->SetRange(-3.*hPulls3[i].GetRMS(),3.*hPulls3[i].GetRMS());
785         hPulls3[i].Fit("g","R,Q,N");
786         pulls[3].SetParam(i,g->GetParameter(2));
787       } else pulls[3].SetParam(i,-1.);
788       if(hPulls4[i].GetEntries()>10) {
789         g->SetRange(-3.*hPulls4[i].GetRMS(),3.*hPulls4[i].GetRMS());
790         hPulls4[i].Fit("g","R,Q,N");
791         pulls[4].SetParam(i,g->GetParameter(2));
792       } else pulls[4].SetParam(i,-1.);
793     } // loop on bins
794
795
796     switch (part) {
797     case 0: // pions
798       for(Int_t i=0;i<5;i++) {
799         fPullsPi[i].~AliTPCkineGrid();
800         new(&fPullsPi[i]) AliTPCkineGrid(pulls[i]);
801       }
802       break;
803     case 1: // kaons
804       for(Int_t i=0;i<5;i++) {
805         fPullsKa[i].~AliTPCkineGrid();
806         new(&fPullsKa[i]) AliTPCkineGrid(pulls[i]);
807       }
808       break;
809     case 2: // protons
810       for(Int_t i=0;i<5;i++) {
811         fPullsPr[i].~AliTPCkineGrid();
812         new(&fPullsPr[i]) AliTPCkineGrid(pulls[i]);
813       }
814       break;
815     case 3: // electrons
816       for(Int_t i=0;i<5;i++) {
817         fPullsEl[i].~AliTPCkineGrid();
818         new(&fPullsEl[i]) AliTPCkineGrid(pulls[i]);
819       }
820       break;
821     case 4: // muons
822       for(Int_t i=0;i<5;i++) {
823         fPullsMu[i].~AliTPCkineGrid();
824         new(&fPullsMu[i]) AliTPCkineGrid(pulls[i]);
825         //cerr<<" mu  pulls "<<i<<"  "<<fPullsMu[i].GetParam(0)<<endl;
826       }
827       break;
828     }
829
830     delete [] hPulls0;
831     delete [] hPulls1;
832     delete [] hPulls2;
833     delete [] hPulls3;
834     delete [] hPulls4;
835     
836   } // loop on particle species
837
838   // write pulls to file
839   WritePulls(outName);
840
841
842   return;
843 }
844 //-----------------------------------------------------------------------------
845 void AliTPCtrackerParam::AnalyzeResolutions(Int_t pdg) {
846 //-----------------------------------------------------------------------------
847 // This function computes the resolutions:
848 // - dy 
849 // - dC 
850 // - dPt/Pt
851 // as a function of Pt
852 // Input file is CovMatrix_AllEvts.root.  
853 //-----------------------------------------------------------------------------
854
855   /*
856   // create a chain with compared tracks
857   TChain *cmptrkchain = new ("cmptrktree");
858   cmptrkchain.Add("CovMatrix_AllEvts.root");
859   //cmptrkchain.Add("CovMatrix_AllEvts_1.root");
860   //cmptrkchain.Add("CovMatrix_AllEvts_2.root");
861   //cmptrkchain.Add("CovMatrix_AllEvts_3.root");
862   */
863
864
865   TFile *infile = TFile::Open("CovMatrix_AllEvts.root");
866   TTree *cmptrktree = (TTree*)infile->Get("cmptrktree");
867
868   COMPTRACK cmptrk; 
869   cmptrktree->SetBranchAddress("comptracks",&cmptrk);
870   Int_t entries = (Int_t)cmptrktree->GetEntries(); 
871   cerr<<" Number of entries: "<<entries<<endl;
872
873
874   Int_t bin = 0;
875
876   InitializeKineGrid("DB");
877   InitializeKineGrid("eff");
878
879   SetParticle(pdg);
880
881   const Int_t knPtBins = fEff->GetPointsPt();
882   cerr<<"knPtBins = "<<knPtBins<<endl; 
883   Double_t *dP0     = new Double_t[knPtBins];
884   Double_t *dP4     = new Double_t[knPtBins];
885   Double_t *dPtToPt = new Double_t[knPtBins];
886   Double_t *pt      = new Double_t[knPtBins];
887   fEff->GetArrayPt(pt);
888
889
890   TH1F *hDumP0 = new TH1F("nameP0","dy",100,-0.3,0.3);
891   TH1F *hDumP4 = new TH1F("nameP4","dC",100,-0.0005,0.0005);
892   TH1F *hDumPt = new TH1F("namePt","dp_{T}/p_{T}",100,-0.5,0.5);
893
894   TF1 *g = new TF1("g","gaus");
895
896   // create histograms for the all the bins
897   TH1F *hP0=NULL;
898   TH1F *hP4=NULL;
899   TH1F *hPt=NULL;
900
901   hP0 = new TH1F[knPtBins]; 
902   hP4 = new TH1F[knPtBins]; 
903   hPt = new TH1F[knPtBins]; 
904
905   for(Int_t i=0; i<knPtBins; i++) {
906     hP0[i] = *hDumP0;
907     hP4[i] = *hDumP4;
908     hPt[i] = *hDumPt;
909   }
910
911   // loop on chain entries 
912   for(Int_t i=0; i<entries; i++) {
913     cmptrktree->GetEvent(i);
914     if(TMath::Abs(cmptrk.pdg)!=pdg) continue;
915     // fill histograms with the residuals
916     bin = (Int_t)fDBgrid->GetBin(cmptrk.pt,cmptrk.eta)/fDBgrid->GetBinsEta();
917     //cerr<<" pt "<<cmptrk.pt<<"   eta "<<cmptrk.eta<<"   bin "<<bin<<endl; 
918     hP0[bin].Fill(cmptrk.dP0);
919     hP4[bin].Fill(cmptrk.dP4);
920     hPt[bin].Fill(cmptrk.dpt/cmptrk.pt);
921   } // loop on chain entries
922
923
924   TCanvas *cP0res = new TCanvas("cP0res","cP0res",0,0,1200,700);
925   cP0res->Divide(5,2);
926   TCanvas *cP4res = new TCanvas("cP4res","cP4res",0,0,1200,700);
927   cP4res->Divide(5,2);
928   TCanvas *cPtres = new TCanvas("cPtres","cPtres",0,0,1200,700);
929   cPtres->Divide(5,2);
930
931   // Draw histograms
932   for(Int_t i=0; i<knPtBins; i++) {
933     cP0res->cd(i+1); hP0[i].Draw();
934     cP4res->cd(i+1); hP4[i].Draw();
935     cPtres->cd(i+1); hPt[i].Draw();
936   }
937
938
939   // compute the sigma of the distributions
940   for(Int_t i=0; i<knPtBins; i++) {
941     if(hP0[i].GetEntries()>10) {
942       g->SetRange(-3.*hP0[i].GetRMS(),3.*hP0[i].GetRMS());
943       hP0[i].Fit("g","R,Q,N");
944       dP0[i] = g->GetParameter(2);
945     } else dP0[i] = 0.;
946     if(hP4[i].GetEntries()>10) {
947       g->SetRange(-3.*hP4[i].GetRMS(),3.*hP4[i].GetRMS());
948       hP4[i].Fit("g","R,Q,N");
949       dP4[i] = g->GetParameter(2);
950     } else dP4[i] = 0.;
951     if(hPt[i].GetEntries()>10) {
952       g->SetRange(-3.*hPt[i].GetRMS(),3.*hPt[i].GetRMS());
953       hPt[i].Fit("g","R,Q,N");
954       dPtToPt[i] = 100.*g->GetParameter(2);
955     } else dPtToPt[i] = 0.;
956   } // loop on bins
957
958   
959   TGraph *grdP0 = new TGraph(knPtBins,pt,dP0);
960   TGraph *grdP4 = new TGraph(knPtBins,pt,dP4);
961   TGraph *grdPtToPt = new TGraph(knPtBins,pt,dPtToPt);
962
963   grdP0->SetMarkerStyle(20); grdP0->SetMarkerColor(2); grdP0->SetMarkerSize(1.5);
964   grdP4->SetMarkerStyle(21); grdP4->SetMarkerColor(3); grdP4->SetMarkerSize(1.5);
965   grdPtToPt->SetMarkerStyle(22); grdPtToPt->SetMarkerColor(4); grdPtToPt->SetMarkerSize(1.5);
966
967   // Plot Results
968   gStyle->SetOptStat(0);
969   TCanvas *c1 = new TCanvas("c1","dP0",0,0,900,700);
970   c1->SetLogx();
971   c1->SetGridx();
972   c1->SetGridy();
973
974   TH2F *frame1 = new TH2F("frame1","y resolution VS p_{T} in TPC",2,0.1,30,2,0,0.1);
975   frame1->SetXTitle("p_{T} [GeV/c]");
976   frame1->SetYTitle("#sigma(y) [cm]");
977   frame1->Draw();
978   grdP0->Draw("P");
979
980
981   TCanvas *c2 = new TCanvas("c2","dP4",0,0,900,700);
982   c2->SetLogx();
983   c2->SetGridx();
984   c2->SetGridy();
985
986   TH2F *frame2 = new TH2F("frame2","C resolution VS p_{T} in TPC",2,0.1,30,2,0,0.0001);
987   frame2->SetXTitle("p_{T} [GeV/c]");
988   frame2->SetYTitle("#sigma(C) [1/cm]");
989   frame2->Draw();
990   grdP4->Draw("P");
991
992   TCanvas *c3 = new TCanvas("c3","dPtToPt",0,0,900,700);
993   c3->SetLogx();
994   c3->SetLogy();
995   c3->SetGridx();
996   c3->SetGridy();
997
998   TH2F *frame3 = new TH2F("frame3","Relative p_{T} resolution VS p_{T} in TPC",2,0.1,30,2,0.1,30.);
999   frame3->SetXTitle("p_{T} [GeV/c]");
1000   frame3->SetYTitle("dp_{T}/p_{T} (%)");
1001   frame3->Draw();
1002   grdPtToPt->Draw("P");
1003
1004
1005   delete [] dP0;
1006   delete [] dP4;
1007   delete [] dPtToPt;
1008   delete [] pt;
1009
1010   
1011   delete [] hP0;
1012   delete [] hP4;
1013   delete [] hPt;
1014   
1015   return;
1016 }
1017 //-----------------------------------------------------------------------------
1018 void AliTPCtrackerParam::BuildTrack(AliTPCseedGeant *s,Int_t ch) {  
1019 //-----------------------------------------------------------------------------
1020 // This function uses GEANT info to set true track parameters
1021 //-----------------------------------------------------------------------------
1022   Double_t xref = s->GetXL();
1023   Double_t xx[5],cc[15];
1024   cc[0]=cc[2]=cc[5]=cc[9]=cc[14]=10.;
1025   cc[1]=cc[3]=cc[4]=cc[6]=cc[7]=cc[8]=cc[10]=cc[11]=cc[12]=cc[13]=0.;
1026   
1027   // Magnetic field
1028   TVector3 bfield(0.,0.,fBz);
1029   
1030   
1031   // radius [cm] of track projection in (x,y) 
1032   Double_t rho = s->GetPt()*100./0.299792458/bfield.Z();
1033   // center of track projection in local reference frame
1034   TVector3 sMom,sPos;
1035
1036
1037   // position (local) and momentum (local) at the seed
1038   // in the bending plane (z=0)
1039   sPos.SetXYZ(s->GetXL(),s->GetYL(),0.);
1040   sMom.SetXYZ(s->GetPx()*TMath::Cos(s->GetAlpha())+s->GetPy()*TMath::Sin(s->GetAlpha()),-s->GetPx()*TMath::Sin(s->GetAlpha())+s->GetPy()*TMath::Cos(s->GetAlpha()),0.);
1041   TVector3 vrho = sMom.Cross(bfield);
1042   vrho *= ch;
1043   vrho.SetMag(rho);
1044
1045   TVector3 vcenter = sPos+vrho;
1046
1047   Double_t x0 = vcenter.X();
1048
1049   // fX     = xref         X-coordinate of this track (reference plane)
1050   // fAlpha = Alpha        Rotation angle the local (TPC sector) 
1051   // fP0    = YL           Y-coordinate of a track
1052   // fP1    = ZG           Z-coordinate of a track
1053   // fP2    = C*x0         x0 is center x in rotated frame
1054   // fP3    = Tgl          tangent of the track momentum dip angle
1055   // fP4    = C            track curvature
1056   xx[0] = s->GetYL();
1057   xx[1] = s->GetZL();
1058   xx[3] = s->GetPz()/s->GetPt();
1059   xx[4] = -ch/rho;
1060   xx[2] = xx[4]*x0;
1061
1062   // create the object AliTPCtrack    
1063   AliTPCtrack track(0,xx,cc,xref,s->GetAlpha());
1064   new(&fTrack) AliTPCtrack(track);
1065
1066   return;
1067 }
1068 //-----------------------------------------------------------------------------
1069 Int_t AliTPCtrackerParam::CheckLabel(AliTPCseedGeant *s,Int_t nPart,
1070                                      Double_t *ptkine,Double_t *pzkine) const {
1071 //-----------------------------------------------------------------------------
1072 // This function checks if the label of the seed has been correctly 
1073 // assigned (to do only for pp charm production with AliRoot v3-08-02)
1074 //-----------------------------------------------------------------------------
1075
1076   Int_t sLabel = s->GetLabel();
1077   Double_t sPt = s->GetPt();
1078   Double_t sPz = s->GetPz();
1079   
1080   // check if the label is correct (comparing momentum)
1081   if(sLabel<nPart && 
1082      TMath::Abs(sPt-ptkine[sLabel])*
1083      TMath::Abs(sPz-pzkine[sLabel])<0.001) return 0;
1084
1085   if((sLabel-30)>=nPart) return 1;  
1086
1087   Double_t diff=0,mindiff=1000.;
1088   Int_t bestLabel=0;
1089
1090   for(Int_t i=sLabel-30; i<sLabel; i++) {
1091     if(i<0 || i>=nPart) continue;
1092     diff = TMath::Abs(sPt-ptkine[i])*TMath::Abs(sPz-pzkine[i]); 
1093     if(diff<mindiff) { mindiff = diff; bestLabel = i; }
1094   }
1095
1096   if(mindiff>0.001) return 1;
1097   s->SetLabel(bestLabel);
1098
1099   return 0;
1100 }
1101 //-----------------------------------------------------------------------------
1102 void AliTPCtrackerParam::CompareTPCtracks(
1103                            const Char_t* galiceName,
1104                            const Char_t* trkGeaName,
1105                            const Char_t* trkKalName,
1106                            const Char_t* covmatName,
1107                            const Char_t* tpceffasciiName,
1108                            const Char_t* tpceffrootName) {
1109 //-----------------------------------------------------------------------------
1110 // This function compares tracks from TPC Kalman Filter V2 with 
1111 // geant tracks at TPC 1st hit. It gives:
1112 //   - a tree with Kalman cov. matrix elements, resolutions, dEdx
1113 //   - the efficiencies as a function of particle type, pT, eta
1114 //-----------------------------------------------------------------------------
1115
1116   TFile *kalFile    = TFile::Open(trkKalName);
1117   TFile *geaFile    = TFile::Open(trkGeaName);
1118   TFile *galiceFile = TFile::Open(galiceName);
1119
1120   // get the AliRun object
1121   AliRun *gAlice = (AliRun*)galiceFile->Get("gAlice");
1122
1123
1124   // create the tree for comparison results
1125   COMPTRACK cmptrk;
1126   TTree *cmptrktree = new TTree("cmptrktree","results of track comparison");
1127   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");
1128
1129   InitializeKineGrid("eff");
1130   InitializeKineGrid("DB");
1131   Int_t   effBins = fEffPi.GetTotPoints();
1132   Int_t effBinsPt = fEffPi.GetPointsPt();
1133   Double_t    *pt = new Double_t[effBinsPt];
1134   fEffPi.GetArrayPt(pt);
1135
1136   TParticle *part;
1137   Double_t ptgener;
1138   Bool_t   usethis;
1139   Int_t    label;
1140   Double_t cc[15],dAlpha;
1141   Int_t    pi=0,ka=0,mu=0,el=0,pr=0;
1142   Int_t   *geaPi = new Int_t[effBins];
1143   Int_t   *geaKa = new Int_t[effBins];
1144   Int_t   *geaPr = new Int_t[effBins];
1145   Int_t   *geaEl = new Int_t[effBins];
1146   Int_t   *geaMu = new Int_t[effBins];
1147   Int_t   *kalPi = new Int_t[effBins];
1148   Int_t   *kalKa = new Int_t[effBins];
1149   Int_t   *kalPr = new Int_t[effBins];
1150   Int_t   *kalEl = new Int_t[effBins];
1151   Int_t   *kalMu = new Int_t[effBins];
1152   Float_t *effPi = new Float_t[effBins];
1153   Float_t *effKa = new Float_t[effBins];
1154   Float_t *effPr = new Float_t[effBins];
1155   Float_t *effEl = new Float_t[effBins];
1156   Float_t *effMu = new Float_t[effBins];
1157   Int_t bin;
1158   for(Int_t j=0; j<effBins; j++) {
1159     geaPi[j]=geaKa[j]=geaPr[j]=geaEl[j]=geaMu[j]=0;
1160     kalPi[j]=kalKa[j]=kalPr[j]=kalEl[j]=kalMu[j]=0;
1161     effPi[j]=effKa[j]=effPr[j]=effEl[j]=effMu[j]=-1.;
1162   }
1163
1164   Char_t tname[100];
1165
1166   // loop on events in file
1167   for(Int_t evt=0; evt<fNevents; evt++) {  
1168     cerr<<"\n  --- Reading tracks for event "<<evt<<" ---\n\n";
1169     sprintf(tname,"TreeT_TPC_%d",evt);
1170     
1171     // particles from TreeK
1172     const Int_t knparticles = gAlice->GetEvent(evt);
1173
1174     Int_t *kalLab = new Int_t[knparticles];
1175     for(Int_t i=0; i<knparticles; i++) kalLab[i] = -1; 
1176  
1177
1178     // tracks from Kalman
1179     TTree *kaltree=(TTree*)kalFile->Get(tname);
1180     if(!kaltree) continue;
1181     AliTPCtrack *kaltrack=new AliTPCtrack; 
1182     kaltree->SetBranchAddress("tracks",&kaltrack);
1183     Int_t kalEntries = (Int_t)kaltree->GetEntries();
1184
1185     // tracks from 1st hit
1186     TTree *geatree=(TTree*)geaFile->Get(tname);
1187     if(!geatree) continue;
1188     AliTPCtrack *geatrack=new AliTPCtrack; 
1189     geatree->SetBranchAddress("tracks",&geatrack);
1190     Int_t geaEntries = (Int_t)geatree->GetEntries();
1191     
1192     cerr<<"+++\n+++ Number of tracks:  TPC Kalman  = "<<kalEntries<<endl<<"+++                    TPC 1st hit = "<<geaEntries<<endl<<"+++\n";
1193     
1194     // set pointers for TPC tracks: 
1195     // the entry number of the track labelled l is stored in kalLab[l]
1196     Int_t fake=0,mult=0;
1197     for (Int_t j=0; j<kalEntries; j++) {
1198       kaltree->GetEvent(j);
1199       if(kaltrack->GetLabel()>=0) { 
1200         if(kalLab[kaltrack->GetLabel()]!=-1) mult++; 
1201         kalLab[kaltrack->GetLabel()] = j;
1202       }
1203       else fake++;
1204     }
1205     
1206     cerr<<"+++ Number of fake tracks in TPC Kalman: "<<fake<<endl;
1207     cerr<<"+++ Number of multiply found tracks in TPC Kalman: "<<mult<<endl;
1208     
1209     /*
1210     // Read the labels of the seeds
1211     char sname[100];
1212     Int_t sLabel,ncol;
1213     Bool_t *hasSeed = new Bool_t[knparticles];
1214     for(Int_t i=0; i<knparticles; i++) hasSeed[i] = kFALSE; 
1215     sprintf(sname,"seedLabels.%d.dat",evt);
1216     FILE *seedFile = fopen(sname,"r");
1217     while(1) {
1218       ncol = fscanf(seedFile,"%d",&sLabel);
1219       if(ncol<1) break;
1220       if(sLabel>0) hasSeed[sLabel]=kTRUE;
1221     }
1222     fclose(seedFile);
1223     */
1224
1225     cerr<<"Doing track comparison...\n";
1226     // loop on tracks at 1st hit
1227     for(Int_t j=0; j<geaEntries; j++) {
1228       geatree->GetEvent(j);
1229       
1230       label = geatrack->GetLabel();
1231       part = (TParticle*)gAlice->GetMCApp()->Particle(label);
1232       
1233       // use only injected tracks with fixed values of pT
1234       ptgener = part->Pt();
1235       usethis = kFALSE;
1236       for(Int_t l=0; l<fEffPi.GetPointsPt(); l++) {
1237         if(TMath::Abs(ptgener-pt[l])<0.01) usethis = kTRUE;
1238       }
1239       if(!usethis) continue;
1240
1241       // check if it has the seed
1242       //if(!hasSeed[label]) continue;
1243
1244       /*
1245       // check if track is entirely contained in a TPC sector
1246       Bool_t out = kFALSE;
1247       for(Int_t l=0; l<10; l++) {
1248         Double_t x = 85. + (250.-85.)*(Double_t)l/9.;
1249         //cerr<<"x "<<x<<"    X "<<geatrack->GetX()<<endl;
1250         Double_t y = geatrack->GetY() + (
1251         TMath::Sqrt(1-(geatrack->GetC()*geatrack->GetX()-geatrack->GetEta())*
1252                       (geatrack->GetC()*geatrack->GetX()-geatrack->GetEta()))
1253        -TMath::Sqrt(1-(geatrack->GetC()*x-geatrack->GetEta())*
1254                       (geatrack->GetC()*x-geatrack->GetEta()))
1255         )/geatrack->GetC();
1256         y = TMath::Abs(y);
1257         //cerr<<"y "<<y<<"    Y "<<geatrack->GetY()<<endl;
1258         if(y > 0.8*x*TMath::Tan(TMath::Pi()/18)) { out = kTRUE; break; }
1259       }
1260       if(out) continue;      
1261       */
1262
1263       cmptrk.pdg = part->GetPdgCode();
1264       cmptrk.eta = part->Eta();
1265       cmptrk.r = TMath::Sqrt(part->Vx()*part->Vx()+part->Vy()*part->Vy());
1266       
1267       cmptrk.pt   = 1/TMath::Abs(geatrack->Get1Pt());
1268       cmptrk.cosl = TMath::Cos(TMath::ATan(geatrack->GetTgl()));
1269       cmptrk.p    = cmptrk.pt/cmptrk.cosl;
1270     
1271
1272       bin = fDBgridPi.GetBin(cmptrk.pt,cmptrk.eta);
1273       cmptrk.bin = bin;
1274       if(abs(cmptrk.pdg)==211)  geaPi[bin]++;  
1275       if(abs(cmptrk.pdg)==321)  geaKa[bin]++;   
1276       if(abs(cmptrk.pdg)==2212) geaPr[bin]++;   
1277       if(abs(cmptrk.pdg)==11)   geaEl[bin]++;  
1278       if(abs(cmptrk.pdg)==13)   geaMu[bin]++; 
1279       
1280       
1281       // check if there is the corresponding track in the TPC kalman and get it
1282       if(kalLab[label]==-1) continue;
1283       kaltree->GetEvent(kalLab[label]);
1284       
1285       // and go on only if it has xref = 84.57 cm (inner pad row)
1286       if(kaltrack->GetX()>90.) continue;
1287       
1288       if(abs(cmptrk.pdg)==211)  { kalPi[bin]++; pi++; }
1289       if(abs(cmptrk.pdg)==321)  { kalKa[bin]++; ka++; }
1290       if(abs(cmptrk.pdg)==2212) { kalPr[bin]++; pr++; }
1291       if(abs(cmptrk.pdg)==11)   { kalEl[bin]++; el++; }
1292       if(abs(cmptrk.pdg)==13)   { kalMu[bin]++; mu++; }
1293       
1294       kaltrack->PropagateTo(geatrack->GetX(),1.e-9,0.);
1295       
1296       cmptrk.dEdx = kaltrack->GetdEdx();
1297       
1298       // compute errors on parameters
1299       dAlpha = kaltrack->GetAlpha()-geatrack->GetAlpha();
1300       if(TMath::Abs(dAlpha)>0.1) { cerr<<" ! WRONG SECTOR !\n"; continue; }
1301       
1302       cmptrk.dP0 = kaltrack->GetY()-geatrack->GetY();
1303       cmptrk.dP1 = kaltrack->GetZ()-geatrack->GetZ();
1304       cmptrk.dP2 = kaltrack->GetEta()-geatrack->GetEta();
1305       cmptrk.dP3 = kaltrack->GetTgl()-geatrack->GetTgl();
1306       cmptrk.dP4 = kaltrack->GetC()-geatrack->GetC();
1307       cmptrk.dpt = 1/kaltrack->Get1Pt()-1/geatrack->Get1Pt();
1308     
1309       // get covariance matrix
1310       // beware: lines 3 and 4 in the matrix are inverted!
1311       kaltrack->GetCovariance(cc);
1312
1313       cmptrk.c00 = cc[0];
1314       cmptrk.c10 = cc[1];
1315       cmptrk.c11 = cc[2];
1316       cmptrk.c20 = cc[3];
1317       cmptrk.c21 = cc[4];
1318       cmptrk.c22 = cc[5];
1319       cmptrk.c30 = cc[10];
1320       cmptrk.c31 = cc[11];
1321       cmptrk.c32 = cc[12];
1322       cmptrk.c33 = cc[14];
1323       cmptrk.c40 = cc[6];
1324       cmptrk.c41 = cc[7];
1325       cmptrk.c42 = cc[8];
1326       cmptrk.c43 = cc[13];
1327       cmptrk.c44 = cc[9];
1328     
1329       // fill tree
1330       cmptrktree->Fill();
1331
1332     } // loop on tracks at TPC 1st hit
1333
1334     delete [] kalLab;
1335     //delete [] hasSeed;
1336
1337   } // end loop on events in file
1338
1339
1340   cerr<<"+++\n+++ Number of compared tracks: "<<pi+ka+el+mu+pr<<endl;
1341   cerr<<"+++ Pions: "<<pi<<", Kaons: "<<ka<<", Protons : "<<pr<<", Electrons: "<<el<<", Muons: "<<mu<<endl;
1342   cerr<<"+++"<<endl;
1343
1344   // Write tree to file
1345   TFile *outfile = new TFile(covmatName,"recreate");
1346   cmptrktree->Write();
1347   outfile->Close();
1348
1349
1350   // Write efficiencies to ascii file
1351   FILE *effFile = fopen(tpceffasciiName,"w");
1352   //fprintf(effFile,"%d\n",kalEntries);
1353   for(Int_t j=0; j<effBins; j++) {
1354     if(geaPi[j]>=100) effPi[j]=(Float_t)kalPi[j]/geaPi[j];
1355     if(geaKa[j]>=100) effKa[j]=(Float_t)kalKa[j]/geaKa[j];
1356     if(geaPr[j]>=100) effPr[j]=(Float_t)kalPr[j]/geaPr[j];
1357     if(geaEl[j]>=500) effEl[j]=(Float_t)kalEl[j]/geaEl[j];
1358     if(geaMu[j]>=100) effMu[j]=(Float_t)kalMu[j]/geaMu[j];
1359     fprintf(effFile,"%f  %f  %f  %f  %f\n",effPi[j],effKa[j],effPr[j],effEl[j],effMu[j]);
1360   }
1361
1362   for(Int_t j=0; j<effBins; j++) {
1363     fprintf(effFile,"%d  %d  %d  %d  %d\n",geaPi[j],geaKa[j],geaPr[j],geaEl[j],geaMu[j]);
1364   }
1365   for(Int_t j=0; j<effBins; j++) {
1366     fprintf(effFile,"%d  %d  %d  %d  %d\n",kalPi[j],kalKa[j],kalPr[j],kalEl[j],kalMu[j]);
1367   }
1368   fclose(effFile);
1369
1370   // Write efficiencies to root file
1371   for(Int_t j=0; j<effBins; j++) {
1372     fEffPi.SetParam(j,(Double_t)effPi[j]);
1373     fEffKa.SetParam(j,(Double_t)effKa[j]);
1374     fEffPr.SetParam(j,(Double_t)effPr[j]);
1375     fEffEl.SetParam(j,(Double_t)effEl[j]);
1376     fEffMu.SetParam(j,(Double_t)effMu[j]);
1377   }
1378   WriteEffs(tpceffrootName);
1379
1380   // delete AliRun object
1381   delete gAlice; gAlice=0;
1382   
1383   // close all input files
1384   kalFile->Close();
1385   geaFile->Close();
1386   galiceFile->Close();
1387
1388   delete [] pt;
1389   delete [] effPi;
1390   delete [] effKa;
1391   delete [] effPr;
1392   delete [] effEl;
1393   delete [] effMu;
1394   delete [] geaPi;
1395   delete [] geaKa;
1396   delete [] geaPr;
1397   delete [] geaEl;
1398   delete [] geaMu;
1399   delete [] kalPi;
1400   delete [] kalKa;
1401   delete [] kalPr;
1402   delete [] kalEl;
1403   delete [] kalMu;
1404
1405   return;
1406 }
1407 //-----------------------------------------------------------------------------
1408 void AliTPCtrackerParam::CookdEdx(Double_t pt,Double_t eta) {
1409 //-----------------------------------------------------------------------------
1410 // This function assigns the track a dE/dx and makes a rough PID
1411 //-----------------------------------------------------------------------------
1412
1413   Double_t mean = fdEdxMean->GetValueAt(pt,eta);
1414   Double_t rms  = fdEdxRMS->GetValueAt(pt,eta);
1415
1416   Double_t dEdx = gRandom->Gaus(mean,rms);
1417
1418   fTrack.SetdEdx(dEdx);
1419
1420   AliTPCtrackParam t(fTrack);
1421
1422   //Very rough PID
1423   Double_t p = TMath::Sqrt(1.+t.GetTgl()*t.GetTgl())*pt;
1424
1425   if (p<0.6) {
1426     if (dEdx < 39.+ 12./(p+0.25)/(p+0.25)) { 
1427       t.AssignMass(AliPID::ParticleMass(AliPID::kPion)); new(&fTrack) AliTPCtrack(t); return;
1428     }
1429     if (dEdx < 39.+ 12./p/p) { 
1430       t.AssignMass(AliPID::ParticleMass(AliPID::kKaon)); new(&fTrack) AliTPCtrack(t); return;
1431     }
1432     t.AssignMass(AliPID::ParticleMass(AliPID::kProton)); new(&fTrack) AliTPCtrack(t); return;
1433   }
1434
1435   if (p<1.2) {
1436     if (dEdx < 39.+ 12./(p+0.25)/(p+0.25)) { 
1437       t.AssignMass(AliPID::ParticleMass(AliPID::kPion)); new(&fTrack) AliTPCtrack(t); return;
1438     }
1439     t.AssignMass(AliPID::ParticleMass(AliPID::kProton)); new(&fTrack) AliTPCtrack(t); return;
1440   }
1441
1442   t.AssignMass(AliPID::ParticleMass(AliPID::kPion)); new(&fTrack) AliTPCtrack(t); return;
1443 }
1444 //-----------------------------------------------------------------------------
1445 void AliTPCtrackerParam::CookTrack(Double_t pt,Double_t eta) {
1446 //-----------------------------------------------------------------------------
1447 // This function deals with covariance matrix and smearing
1448 //-----------------------------------------------------------------------------
1449
1450   COVMATRIX covmat;
1451   Double_t  p,cosl;
1452   Double_t  trkKine[1],trkRegPar[3];    
1453   Double_t  xref,alpha,xx[5],xxsm[5],cc[15];
1454   Int_t     treeEntries;
1455
1456   fCovTree->SetBranchAddress("matrix",&covmat);
1457
1458   // get random entry from the tree
1459   treeEntries = (Int_t)fCovTree->GetEntries();
1460   fCovTree->GetEvent(gRandom->Integer(treeEntries));
1461
1462   // get P and Cosl from track
1463   cosl = TMath::Cos(TMath::ATan(fTrack.GetTgl()));
1464   p    = 1./TMath::Abs(fTrack.Get1Pt())/cosl;
1465
1466   trkKine[0] = p;
1467
1468   // get covariance matrix from regularized matrix    
1469   for(Int_t l=0;l<3;l++) trkRegPar[l] = (*fRegPar)(0,l);
1470   cc[0] = covmat.c00*RegFunc(trkKine,trkRegPar);
1471   cc[1] = covmat.c10;
1472   for(Int_t l=0;l<3;l++) trkRegPar[l] = (*fRegPar)(1,l);
1473   cc[2] = covmat.c11*RegFunc(trkKine,trkRegPar);
1474   for(Int_t l=0;l<3;l++) trkRegPar[l] = (*fRegPar)(2,l);
1475   cc[3] = covmat.c20*RegFunc(trkKine,trkRegPar);
1476   cc[4] = covmat.c21;
1477   for(Int_t l=0;l<3;l++) trkRegPar[l] = (*fRegPar)(3,l);
1478   cc[5] = covmat.c22*RegFunc(trkKine,trkRegPar);
1479   cc[6] = covmat.c30;
1480   for(Int_t l=0;l<3;l++) trkRegPar[l] = (*fRegPar)(4,l);
1481   cc[7] = covmat.c31*RegFunc(trkKine,trkRegPar);
1482   cc[8] = covmat.c32;
1483   for(Int_t l=0;l<3;l++) trkRegPar[l] = (*fRegPar)(5,l);
1484   cc[9] = covmat.c33*RegFunc(trkKine,trkRegPar);
1485   for(Int_t l=0;l<3;l++) trkRegPar[l] = (*fRegPar)(6,l);
1486   cc[10]= covmat.c40*RegFunc(trkKine,trkRegPar);
1487   cc[11]= covmat.c41;
1488   for(Int_t l=0;l<3;l++) trkRegPar[l] = (*fRegPar)(7,l);
1489   cc[12]= covmat.c42*RegFunc(trkKine,trkRegPar);
1490   cc[13]= covmat.c43;
1491   for(Int_t l=0;l<3;l++) trkRegPar[l] = (*fRegPar)(8,l);
1492   cc[14]= covmat.c44*RegFunc(trkKine,trkRegPar);
1493    
1494   TMatrixD covMatSmear(5,5);
1495     
1496   covMatSmear = GetSmearingMatrix(cc,pt,eta);
1497
1498   // get track original parameters
1499   xref =fTrack.GetX();
1500   alpha=fTrack.GetAlpha();
1501   xx[0]=fTrack.GetY();
1502   xx[1]=fTrack.GetZ();
1503   xx[2]=fTrack.GetX()*fTrack.GetC()-fTrack.GetSnp();
1504   xx[3]=fTrack.GetTgl();
1505   xx[4]=fTrack.GetC();
1506     
1507   // use smearing matrix to smear the original parameters
1508   SmearTrack(xx,xxsm,covMatSmear);
1509     
1510   AliTPCtrack track(0,xxsm,cc,xref,alpha);
1511   new(&fTrack) AliTPCtrack(track);
1512   
1513   return; 
1514 }
1515 //-----------------------------------------------------------------------------
1516 void AliTPCtrackerParam::DrawEffs(const Char_t* inName,Int_t pdg) {
1517 //-----------------------------------------------------------------------------
1518 // This function draws the TPC efficiencies in the [pT,eta] bins
1519 //-----------------------------------------------------------------------------
1520
1521   ReadEffs(inName);
1522   SetParticle(pdg);
1523
1524   const Int_t kn = fEff->GetPointsPt();
1525   Double_t *effsA = new Double_t[kn];
1526   Double_t *effsB = new Double_t[kn];
1527   Double_t *effsC = new Double_t[kn];
1528   Double_t *pt    = new Double_t[kn];
1529
1530   fEff->GetArrayPt(pt);
1531   for(Int_t i=0;i<kn;i++) {
1532     effsA[i] = fEff->GetParam(i,0);
1533     effsB[i] = fEff->GetParam(i,1);
1534     effsC[i] = fEff->GetParam(i,2);
1535   }
1536   
1537   TGraph *grA = new TGraph(kn,pt,effsA);
1538   TGraph *grB = new TGraph(kn,pt,effsB);
1539   TGraph *grC = new TGraph(kn,pt,effsC);
1540
1541   grA->SetMarkerStyle(20); grA->SetMarkerColor(2); grA->SetMarkerSize(1.5);
1542   grB->SetMarkerStyle(21); grB->SetMarkerColor(3); grB->SetMarkerSize(1.5);
1543   grC->SetMarkerStyle(22); grC->SetMarkerColor(4); grC->SetMarkerSize(1.5);
1544
1545   TString title("Distribution of the TPC efficiencies");
1546   switch (pdg) {
1547   case 211:
1548     title.Prepend("PIONS - ");
1549     break;
1550   case 321:
1551     title.Prepend("KAONS - ");
1552     break;
1553   case 2212:
1554     title.Prepend("PROTONS - ");
1555     break;
1556   case 11:
1557     title.Prepend("ELECTRONS - ");
1558     break;
1559   case 13:
1560     title.Prepend("MUONS - ");
1561     break;
1562   }
1563
1564
1565   gStyle->SetOptStat(0);
1566   TCanvas *c = new TCanvas("c","effs",0,0,900,700);
1567   c->SetLogx();
1568   c->SetGridx();
1569   c->SetGridy();
1570
1571   TH2F *frame = new TH2F("frame",title.Data(),2,0.1,30,2,0,1);
1572   frame->SetXTitle("p_{T} [GeV/c]");
1573   frame->Draw();
1574   grA->Draw("P");
1575   grB->Draw("P");
1576   grC->Draw("P");
1577
1578   TLegend *leg = new TLegend(0.2,0.2,0.4,0.4);
1579   leg->AddEntry(grA,"|#eta|<0.3","p");
1580   leg->AddEntry(grB,"0.3<|#eta|<0.6","p");
1581   leg->AddEntry(grC,"0.6<|#eta|<0.9","p");
1582   leg->SetFillColor(0);
1583   leg->Draw();
1584
1585   delete [] effsA;
1586   delete [] effsB;
1587   delete [] effsC;
1588   delete [] pt;
1589
1590   return;
1591 }
1592 //-----------------------------------------------------------------------------
1593 void AliTPCtrackerParam::DrawPulls(const Char_t* inName,Int_t pdg,
1594                                    Int_t par) {
1595 //-----------------------------------------------------------------------------
1596 // This function draws the pulls in the [pT,eta] bins
1597 //-----------------------------------------------------------------------------
1598
1599   ReadPulls(inName);
1600   SetParticle(pdg);
1601
1602   const Int_t kn = (fPulls+par)->GetPointsPt();
1603   Double_t *pullsA = new Double_t[kn];
1604   Double_t *pullsB = new Double_t[kn];
1605   Double_t *pullsC = new Double_t[kn];
1606   Double_t *pt     = new Double_t[kn];
1607   (fPulls+par)->GetArrayPt(pt);  
1608   for(Int_t i=0;i<kn;i++) {
1609     pullsA[i] = (fPulls+par)->GetParam(i,0);
1610     pullsB[i] = (fPulls+par)->GetParam(i,1);
1611     pullsC[i] = (fPulls+par)->GetParam(i,2);
1612   }
1613
1614   TGraph *grA = new TGraph(kn,pt,pullsA);
1615   TGraph *grB = new TGraph(kn,pt,pullsB);
1616   TGraph *grC = new TGraph(kn,pt,pullsC);
1617
1618   grA->SetMarkerStyle(20); grA->SetMarkerColor(2); grA->SetMarkerSize(1.5);
1619   grB->SetMarkerStyle(21); grB->SetMarkerColor(3); grB->SetMarkerSize(1.5);
1620   grC->SetMarkerStyle(22); grC->SetMarkerColor(4); grC->SetMarkerSize(1.5);
1621
1622   TString title("Distribution of the pulls: ");
1623   switch (pdg) {
1624   case 211:
1625     title.Prepend("PIONS - ");
1626     break;
1627   case 321:
1628     title.Prepend("KAONS - ");
1629     break;
1630   case 2212:
1631     title.Prepend("PROTONS - ");
1632     break;
1633   case 11:
1634     title.Prepend("ELECTRONS - ");
1635     break;
1636   case 13:
1637     title.Prepend("MUONS - ");
1638     break;
1639   }
1640   switch (par) {
1641   case 0:
1642     title.Append("y");
1643     break;
1644   case 1:
1645     title.Append("z");
1646     break;
1647   case 2:
1648     title.Append("    #eta");
1649     break;
1650   case 3:
1651     title.Append("tg    #lambda");
1652     break;
1653   case 4:
1654     title.Append("C");
1655     break;
1656   }
1657
1658   gStyle->SetOptStat(0);
1659   TCanvas *c = new TCanvas("c","pulls",0,0,900,700);
1660   c->SetLogx();
1661   c->SetGridx();
1662   c->SetGridy();
1663
1664   TH2F *frame = new TH2F("frame",title.Data(),2,0.1,30,2,0,2);
1665   frame->SetXTitle("p_{T} [GeV/c]");
1666   frame->Draw();
1667   grA->Draw("P");
1668   grB->Draw("P");
1669   grC->Draw("P");
1670
1671   TLegend *leg = new TLegend(0.2,0.2,0.4,0.4);
1672   leg->AddEntry(grA,"|#eta|<0.3","p");
1673   leg->AddEntry(grB,"0.3<|#eta|<0.6","p");
1674   leg->AddEntry(grC,"0.6<|#eta|<0.9","p");
1675   leg->SetFillColor(0);
1676   leg->Draw();
1677
1678   delete [] pullsA;
1679   delete [] pullsB;
1680   delete [] pullsC;
1681   delete [] pt;
1682
1683   return;
1684 }
1685 //-----------------------------------------------------------------------------
1686 TMatrixD AliTPCtrackerParam::GetSmearingMatrix(Double_t* cc,Double_t pt,
1687                                                Double_t eta) const {
1688 //-----------------------------------------------------------------------------
1689 // This function stretches the covariance matrix according to the pulls
1690 //-----------------------------------------------------------------------------
1691
1692   TMatrixD covMat(5,5);
1693
1694   covMat(0,0)=cc[0];
1695   covMat(1,0)=cc[1];  covMat(0,1)=covMat(1,0);
1696   covMat(1,1)=cc[2];
1697   covMat(2,0)=cc[3];  covMat(0,2)=covMat(2,0);
1698   covMat(2,1)=cc[4];  covMat(1,2)=covMat(2,1);
1699   covMat(2,2)=cc[5];
1700   covMat(3,0)=cc[6];  covMat(0,3)=covMat(3,0);
1701   covMat(3,1)=cc[7];  covMat(1,3)=covMat(3,1);
1702   covMat(3,2)=cc[8];  covMat(2,3)=covMat(3,2);
1703   covMat(3,3)=cc[9];
1704   covMat(4,0)=cc[10]; covMat(0,4)=covMat(4,0);
1705   covMat(4,1)=cc[11]; covMat(1,4)=covMat(4,1);
1706   covMat(4,2)=cc[12]; covMat(2,4)=covMat(4,2);
1707   covMat(4,3)=cc[13]; covMat(3,4)=covMat(4,3);
1708   covMat(4,4)=cc[14];
1709
1710
1711   TMatrixD stretchMat(5,5);
1712   for(Int_t k=0;k<5;k++) {
1713     for(Int_t l=0;l<5;l++) {
1714       stretchMat(k,l)=0.;
1715     }
1716   }
1717
1718   for(Int_t i=0;i<5;i++) {
1719     stretchMat(i,i) = (fPulls+i)->GetValueAt(pt,eta); 
1720     if(stretchMat(i,i)==0.) stretchMat(i,i) = 1.;
1721   }
1722
1723   TMatrixD mat(stretchMat,TMatrixD::kMult,covMat);
1724   TMatrixD covMatSmear(mat,TMatrixD::kMult,stretchMat);
1725
1726   return covMatSmear;
1727 }
1728 //-----------------------------------------------------------------------------
1729 void AliTPCtrackerParam::InitializeKineGrid(Option_t* which) {
1730 //-----------------------------------------------------------------------------
1731 // This function initializes ([pt,eta] points) the data members AliTPCkineGrid
1732 // which = "DB"     -> initialize fDBgrid... members
1733 //         "eff"    -> initialize fEff... members
1734 //         "pulls"  -> initialize fPulls... members
1735 //         "dEdx"   -> initialize fdEdx... members
1736 //-----------------------------------------------------------------------------
1737
1738   const char *db     = strstr(which,"DB");
1739   const char *eff    = strstr(which,"eff");
1740   const char *pulls  = strstr(which,"pulls");
1741   const char *dEdx   = strstr(which,"dEdx");
1742
1743
1744   Int_t nEta=0, nPt=0;
1745
1746   Double_t etaPoints[2] = {0.3,0.6};
1747   Double_t etaBins[3]   = {0.15,0.45,0.75};
1748
1749   Double_t ptPoints[9]  = {0.4,0.6,0.8,1.2,1.7,3.,5.,8.,15.};
1750   Double_t ptBins[10]  = {0.3,0.5,0.7,1.,1.5,2.,4.,6.,10.,20.};
1751
1752
1753   Double_t *eta=0,*pt=0;
1754
1755   if(db) {
1756     nEta = 2;
1757     nPt  = 9;
1758     eta  = etaPoints;
1759     pt   = ptPoints;
1760   } else {
1761     nEta = 3;
1762     nPt  = 10;
1763     eta  = etaBins;
1764     pt   = ptBins;
1765   }
1766
1767   AliTPCkineGrid *dummy=0;
1768
1769   if(db) {    
1770     dummy = new AliTPCkineGrid(nPt,nEta,pt,eta);
1771     new(&fDBgridPi) AliTPCkineGrid(*dummy);
1772     new(&fDBgridKa) AliTPCkineGrid(*dummy);
1773     new(&fDBgridPr) AliTPCkineGrid(*dummy);
1774     new(&fDBgridEl) AliTPCkineGrid(*dummy);
1775     new(&fDBgridMu) AliTPCkineGrid(*dummy);
1776     delete dummy;
1777   }
1778   if(eff) {    
1779     dummy = new AliTPCkineGrid(nPt,nEta,pt,eta);
1780     new(&fEffPi) AliTPCkineGrid(*dummy);
1781     new(&fEffKa) AliTPCkineGrid(*dummy);
1782     new(&fEffPr) AliTPCkineGrid(*dummy);
1783     new(&fEffEl) AliTPCkineGrid(*dummy);
1784     new(&fEffMu) AliTPCkineGrid(*dummy);
1785     delete dummy;
1786   }
1787   if(pulls) {    
1788     dummy = new AliTPCkineGrid(nPt,nEta,pt,eta);
1789     for(Int_t i=0;i<5;i++) new(&fPullsPi[i]) AliTPCkineGrid(*dummy);
1790     for(Int_t i=0;i<5;i++) new(&fPullsKa[i]) AliTPCkineGrid(*dummy);
1791     for(Int_t i=0;i<5;i++) new(&fPullsPr[i]) AliTPCkineGrid(*dummy);
1792     for(Int_t i=0;i<5;i++) new(&fPullsEl[i]) AliTPCkineGrid(*dummy);
1793     for(Int_t i=0;i<5;i++) new(&fPullsMu[i]) AliTPCkineGrid(*dummy);
1794     delete dummy;
1795   }
1796   if(dEdx) {    
1797     dummy = new AliTPCkineGrid(nPt,nEta,pt,eta);
1798     new(&fdEdxMeanPi) AliTPCkineGrid(*dummy);
1799     new(&fdEdxRMSPi) AliTPCkineGrid(*dummy);
1800     new(&fdEdxMeanKa) AliTPCkineGrid(*dummy);
1801     new(&fdEdxRMSKa) AliTPCkineGrid(*dummy);
1802     new(&fdEdxMeanPr) AliTPCkineGrid(*dummy);
1803     new(&fdEdxRMSPr) AliTPCkineGrid(*dummy);
1804     new(&fdEdxMeanEl) AliTPCkineGrid(*dummy);
1805     new(&fdEdxRMSEl) AliTPCkineGrid(*dummy);
1806     new(&fdEdxMeanMu) AliTPCkineGrid(*dummy);
1807     new(&fdEdxRMSMu) AliTPCkineGrid(*dummy);
1808     delete dummy; 
1809   }
1810
1811   return;
1812 }
1813 //-----------------------------------------------------------------------------
1814 void AliTPCtrackerParam::MakeDataBase() {
1815 //-----------------------------------------------------------------------------
1816 // This function creates the DB file and store in it:
1817 //  - TPC Efficiencies for pi,ka,pr,el,mu     (AliTPCkineGrid class)
1818 //  - Pulls for pi,ka,el                      (AliTPCkineGrid class)
1819 //  - Regularization parameters for pi,ka,el  (TMatrixD class)
1820 //  - dE/dx parameterization for pi,ka,pr,el  (AliTPCkineGrid class)  
1821 //  - Regularized cov. matrices for pi,ka,el  (COVMATRIX structure)
1822 //-----------------------------------------------------------------------------
1823
1824   // define some file names
1825   const char *effFile   ="TPCeff.root";
1826   const char *pullsFile ="pulls.root";
1827   const char *regPiFile ="regPi.root";
1828   const char *regKaFile ="regKa.root";
1829   const char *regPrFile ="regPr.root";
1830   const char *regElFile ="regEl.root";
1831   const char *regMuFile ="regMu.root";
1832   const char *dEdxPiFile="dEdxPi.root";
1833   const char *dEdxKaFile="dEdxKa.root";
1834   const char *dEdxPrFile="dEdxPr.root";
1835   const char *dEdxElFile="dEdxEl.root";
1836   const char *dEdxMuFile="dEdxMu.root";
1837   const char *cmFile    ="CovMatrix_AllEvts.root";
1838   /*
1839   Char_t *cmFile1  ="CovMatrix_AllEvts_1.root";
1840   Char_t *cmFile2  ="CovMatrix_AllEvts_2.root";
1841   Char_t *cmFile3  ="CovMatrix_AllEvts_3.root";
1842   */
1843
1844   // store the effieciencies
1845   ReadEffs(effFile);
1846   WriteEffs(fDBfileName.Data());
1847
1848   // store the pulls
1849   ReadPulls(pullsFile);
1850   WritePulls(fDBfileName.Data());
1851
1852   //store the regularization parameters
1853   ReadRegParams(regPiFile,211);
1854   WriteRegParams(fDBfileName.Data(),211);
1855   ReadRegParams(regKaFile,321);
1856   WriteRegParams(fDBfileName.Data(),321);
1857   ReadRegParams(regPrFile,2212);
1858   WriteRegParams(fDBfileName.Data(),2212);
1859   ReadRegParams(regElFile,11);
1860   WriteRegParams(fDBfileName.Data(),11);
1861   ReadRegParams(regMuFile,13);
1862   WriteRegParams(fDBfileName.Data(),13);
1863
1864   // store the dEdx parameters
1865   ReaddEdx(dEdxPiFile,211);
1866   WritedEdx(fDBfileName.Data(),211);
1867   ReaddEdx(dEdxKaFile,321);
1868   WritedEdx(fDBfileName.Data(),321);
1869   ReaddEdx(dEdxPrFile,2212);
1870   WritedEdx(fDBfileName.Data(),2212);
1871   ReaddEdx(dEdxElFile,11);
1872   WritedEdx(fDBfileName.Data(),11);
1873   ReaddEdx(dEdxMuFile,13);
1874   WritedEdx(fDBfileName.Data(),13);
1875
1876
1877   //
1878   // store the regularized covariance matrices
1879   //
1880   InitializeKineGrid("DB");
1881
1882   const Int_t knBinsPi = fDBgridPi.GetTotBins();
1883   const Int_t knBinsKa = fDBgridKa.GetTotBins();
1884   const Int_t knBinsPr = fDBgridPr.GetTotBins();
1885   const Int_t knBinsEl = fDBgridEl.GetTotBins();
1886   const Int_t knBinsMu = fDBgridMu.GetTotBins();
1887
1888
1889   // create the trees for cov. matrices
1890   // trees for pions
1891   TTree *covTreePi_ = NULL;
1892   covTreePi_ = new TTree[knBinsPi]; 
1893   // trees for kaons
1894   TTree *covTreeKa_ = NULL;
1895   covTreeKa_ = new TTree[knBinsKa]; 
1896   // trees for protons
1897   TTree *covTreePr_ = NULL;
1898   covTreePr_ = new TTree[knBinsPr]; 
1899   // trees for electrons
1900   TTree *covTreeEl_ = NULL;
1901   covTreeEl_ = new TTree[knBinsEl]; 
1902   // trees for muons
1903   TTree *covTreeMu_ = NULL;
1904   covTreeMu_ = new TTree[knBinsMu]; 
1905
1906   Char_t hname[100], htitle[100];
1907   COVMATRIX covmat;
1908
1909   
1910   for(Int_t i=0; i<knBinsPi; i++) {
1911     sprintf(hname,"CovTreePi_bin%d",i);
1912     sprintf(htitle,"Tree with cov matrix elements for bin %d",i);
1913     covTreePi_[i].SetName(hname); covTreePi_[i].SetTitle(htitle);
1914     covTreePi_[i].Branch("matrix",&covmat,"c00/D:c10:c11:c20:c21:c22:c30:c31:c32:c33:c40:c41:c42:c43:c44",5000000);
1915   }
1916   for(Int_t i=0; i<knBinsKa; i++) {
1917     sprintf(hname,"CovTreeKa_bin%d",i);
1918     sprintf(htitle,"Tree with cov matrix elements for bin %d",i);
1919     covTreeKa_[i].SetName(hname); covTreeKa_[i].SetTitle(htitle);
1920     covTreeKa_[i].Branch("matrix",&covmat,"c00/D:c10:c11:c20:c21:c22:c30:c31:c32:c33:c40:c41:c42:c43:c44",1000000);
1921   }
1922   for(Int_t i=0; i<knBinsPr; i++) {
1923     sprintf(hname,"CovTreePr_bin%d",i);
1924     sprintf(htitle,"Tree with cov matrix elements for bin %d",i);
1925     covTreePr_[i].SetName(hname); covTreePr_[i].SetTitle(htitle);
1926     covTreePr_[i].Branch("matrix",&covmat,"c00/D:c10:c11:c20:c21:c22:c30:c31:c32:c33:c40:c41:c42:c43:c44",1000000);
1927   }
1928   for(Int_t i=0; i<knBinsEl; i++) {
1929     sprintf(hname,"CovTreeEl_bin%d",i);
1930     sprintf(htitle,"Tree with cov matrix elements for bin %d",i);
1931     covTreeEl_[i].SetName(hname); covTreeEl_[i].SetTitle(htitle);
1932     covTreeEl_[i].Branch("matrix",&covmat,"c00/D:c10:c11:c20:c21:c22:c30:c31:c32:c33:c40:c41:c42:c43:c44",1000000);
1933   }
1934   for(Int_t i=0; i<knBinsMu; i++) {
1935     sprintf(hname,"CovTreeMu_bin%d",i);
1936     sprintf(htitle,"Tree with cov matrix elements for bin %d",i);
1937     covTreeMu_[i].SetName(hname); covTreeMu_[i].SetTitle(htitle);
1938     covTreeMu_[i].Branch("matrix",&covmat,"c00/D:c10:c11:c20:c21:c22:c30:c31:c32:c33:c40:c41:c42:c43:c44",1000000);
1939   }
1940
1941   /*  
1942   // create the chain with the compared tracks
1943   TChain *cmptrktree = new TChain("cmptrktree");
1944   cmptrkchain.Add(cmFile1);
1945   cmptrkchain.Add(cmFile2);
1946   cmptrkchain.Add(cmFile3);
1947   */
1948
1949   TFile *infile = TFile::Open(cmFile);
1950   TTree *cmptrktree = (TTree*)infile->Get("cmptrktree");
1951  
1952   COMPTRACK cmptrk; 
1953   cmptrktree->SetBranchAddress("comptracks",&cmptrk);
1954   Int_t entries = (Int_t)cmptrktree->GetEntries(); 
1955   cerr<<" Number of entries: "<<entries<<endl;
1956
1957   Int_t trkPdg,trkBin;
1958   Double_t trkKine[1],trkRegPar[3]; 
1959   Int_t *nPerBinPi = new Int_t[knBinsPi];
1960   for(Int_t k=0;k<knBinsPi;k++) nPerBinPi[k]=0;
1961   Int_t *nPerBinKa = new Int_t[knBinsKa];
1962   for(Int_t k=0;k<knBinsKa;k++) nPerBinKa[k]=0;
1963   Int_t *nPerBinMu = new Int_t[knBinsMu];
1964   for(Int_t k=0;k<knBinsMu;k++) nPerBinMu[k]=0;
1965   Int_t *nPerBinEl = new Int_t[knBinsEl];
1966   for(Int_t k=0;k<knBinsEl;k++) nPerBinEl[k]=0;
1967   Int_t *nPerBinPr = new Int_t[knBinsPr];
1968   for(Int_t k=0;k<knBinsPr;k++) nPerBinPr[k]=0;
1969
1970   // loop on chain entries 
1971   for(Int_t l=0; l<entries; l++) {
1972     if(l % 10000 == 0) cerr<<"--- Processing track "<<l<<" of "<<entries<<" ---"<<endl;
1973     // get the event
1974     cmptrktree->GetEvent(l);
1975     // get the pdg code
1976     trkPdg = TMath::Abs(cmptrk.pdg);
1977     // use only pions, kaons, protons, electrons, muons
1978     if(trkPdg!=211 && trkPdg!=321 && trkPdg!=2212 && trkPdg!=11 && trkPdg!=13) continue;
1979     SetParticle(trkPdg);
1980     trkBin = fDBgrid->GetBin(cmptrk.pt,cmptrk.eta);
1981     //cerr<<cmptrk.pt<<"  "<<cmptrk.eta<<"  "<<trkBin<<endl;
1982
1983     if(trkPdg==211 && nPerBinPi[trkBin]>=5000) continue;
1984     if(trkPdg==321 && nPerBinKa[trkBin]>=5000) continue;
1985     if(trkPdg==2212 && nPerBinPr[trkBin]>=5000) continue;
1986     if(trkPdg==11  && nPerBinEl[trkBin]>=5000) continue;
1987     if(trkPdg==13 && nPerBinMu[trkBin]>=5000) continue;
1988
1989     trkKine[0] = cmptrk.p;
1990     
1991     // get regularized covariance matrix
1992     for(Int_t k=0;k<3;k++) trkRegPar[k] = (*fRegPar)(0,k);
1993     covmat.c00 = cmptrk.c00/RegFunc(trkKine,trkRegPar);
1994     covmat.c10 = cmptrk.c10;
1995     for(Int_t k=0;k<3;k++) trkRegPar[k] = (*fRegPar)(1,k);
1996     covmat.c11 = cmptrk.c11/RegFunc(trkKine,trkRegPar);
1997     for(Int_t k=0;k<3;k++) trkRegPar[k] = (*fRegPar)(2,k);
1998     covmat.c20 = cmptrk.c20/RegFunc(trkKine,trkRegPar);
1999     covmat.c21 = cmptrk.c21;
2000     for(Int_t k=0;k<3;k++) trkRegPar[k] = (*fRegPar)(3,k);
2001     covmat.c22 = cmptrk.c22/RegFunc(trkKine,trkRegPar);
2002     covmat.c30 = cmptrk.c30;
2003     for(Int_t k=0;k<3;k++) trkRegPar[k] = (*fRegPar)(4,k);
2004     covmat.c31 = cmptrk.c31/RegFunc(trkKine,trkRegPar);
2005     covmat.c32 = cmptrk.c32;
2006     for(Int_t k=0;k<3;k++) trkRegPar[k] = (*fRegPar)(5,k);
2007     covmat.c33 = cmptrk.c33/RegFunc(trkKine,trkRegPar);
2008     for(Int_t k=0;k<3;k++) trkRegPar[k] = (*fRegPar)(6,k);
2009     covmat.c40 = cmptrk.c40/RegFunc(trkKine,trkRegPar);
2010     covmat.c41 = cmptrk.c41;
2011     for(Int_t k=0;k<3;k++) trkRegPar[k] = (*fRegPar)(7,k);
2012     covmat.c42 = cmptrk.c42/RegFunc(trkKine,trkRegPar);
2013     covmat.c43 = cmptrk.c43;
2014     for(Int_t k=0;k<3;k++) trkRegPar[k] = (*fRegPar)(8,k);
2015     covmat.c44 = cmptrk.c44/RegFunc(trkKine,trkRegPar);
2016
2017     // fill the tree
2018     switch (trkPdg) {
2019     case 211: // pions
2020       covTreePi_[trkBin].Fill();
2021       nPerBinPi[trkBin]++;
2022       break;
2023     case 321: // kaons
2024       covTreeKa_[trkBin].Fill();
2025       nPerBinKa[trkBin]++;
2026       break;
2027     case 2212: // protons
2028       covTreePr_[trkBin].Fill();
2029       nPerBinPr[trkBin]++;
2030       break;
2031     case 11: // electrons
2032       covTreeEl_[trkBin].Fill();
2033       nPerBinEl[trkBin]++;
2034       break;
2035     case 13: // muons
2036       covTreeMu_[trkBin].Fill();
2037       nPerBinMu[trkBin]++;
2038       break;
2039     }
2040   } // loop on chain entries
2041
2042   // store all trees the DB file
2043   TFile *dbfile = new TFile(fDBfileName.Data(),"update");
2044   dbfile->mkdir("CovMatrices");
2045   gDirectory->cd("/CovMatrices");
2046   gDirectory->mkdir("Pions");
2047   gDirectory->mkdir("Kaons");
2048   gDirectory->mkdir("Protons");
2049   gDirectory->mkdir("Electrons");
2050   gDirectory->mkdir("Muons");
2051   // store pions
2052   gDirectory->cd("/CovMatrices/Pions");
2053   fDBgridPi.SetName("DBgridPi"); fDBgridPi.Write();
2054   for(Int_t i=0;i<knBinsPi;i++) covTreePi_[i].Write();
2055   // store kaons
2056   gDirectory->cd("/CovMatrices/Kaons");
2057   fDBgridKa.SetName("DBgridKa"); fDBgridKa.Write();
2058   for(Int_t i=0;i<knBinsKa;i++) covTreeKa_[i].Write();
2059   // store kaons
2060   gDirectory->cd("/CovMatrices/Protons");
2061   fDBgridPr.SetName("DBgridPr"); fDBgridPr.Write();
2062   for(Int_t i=0;i<knBinsPr;i++) covTreePr_[i].Write();
2063   // store electrons
2064   gDirectory->cd("/CovMatrices/Electrons");
2065   fDBgridEl.SetName("DBgridEl"); fDBgridEl.Write();
2066   for(Int_t i=0;i<knBinsEl;i++) covTreeEl_[i].Write();
2067   // store kaons
2068   gDirectory->cd("/CovMatrices/Muons");
2069   fDBgridMu.SetName("DBgridMu"); fDBgridMu.Write();
2070   for(Int_t i=0;i<knBinsMu;i++) covTreeMu_[i].Write();
2071
2072   dbfile->Close();
2073   delete [] nPerBinPi;
2074   delete [] nPerBinKa;
2075   delete [] nPerBinPr;
2076   delete [] nPerBinEl;
2077   delete [] nPerBinMu; 
2078  
2079   return;
2080 }
2081 //-----------------------------------------------------------------------------
2082 void AliTPCtrackerParam::MakeSeedsFromHits(AliTPC *tpc,TTree *th,
2083                                            TObjArray &seedArray) const {
2084 //-----------------------------------------------------------------------------
2085 // This function makes the seeds for tracks from the 1st hits in the TPC
2086 //-----------------------------------------------------------------------------
2087
2088   Double_t xg,yg,zg,px,py,pz,pt;
2089   Int_t label;
2090   Int_t nTracks=(Int_t)th->GetEntries();
2091
2092   cerr<<"+++\n+++ Number of \"primary tracks\"(entries in TreeH): "<<nTracks<<
2093          "\n+++\n\n";
2094    
2095   AliTPChit *tpcHit=0;
2096
2097   // loop over entries in TreeH
2098   for(Int_t l=0; l<nTracks; l++) {
2099     if(l%1000==0) cerr<<"  --- Processing primary track "
2100                       <<l<<" of "<<nTracks<<" ---\r";
2101     tpc->ResetHits();
2102     th->GetEvent(l);
2103     // Get FirstHit
2104     tpcHit=(AliTPChit*)tpc->FirstHit(-1);
2105     for( ; tpcHit; tpcHit=(AliTPChit*)tpc->NextHit() ) {
2106       if(tpcHit->fQ !=0.) continue;
2107       // Get particle momentum at hit
2108       px=tpcHit->X(); py=tpcHit->Y(); pz=tpcHit->Z();
2109
2110       pt=TMath::Sqrt(px*px+py*py);
2111       // reject hits with Pt<mag*0.45 GeV/c
2112       if(pt<(fBz*0.45)) continue;
2113
2114       // Get track label
2115       label=tpcHit->Track();
2116       
2117       if((tpcHit=(AliTPChit*)tpc->NextHit())==0) break;
2118       if(tpcHit->fQ != 0.) continue;
2119       // Get global coordinates of hit
2120       xg=tpcHit->X(); yg=tpcHit->Y(); zg=tpcHit->Z();
2121       if(TMath::Sqrt(xg*xg+yg*yg)>90.) continue;
2122
2123       // create the seed
2124       AliTPCseedGeant *ioseed = new AliTPCseedGeant(xg,yg,zg,px,py,pz,label);
2125
2126       // reject tracks which are not in the TPC acceptance
2127       if(!ioseed->InTPCAcceptance()) { delete ioseed; continue; }
2128
2129       // put seed in array
2130       seedArray.AddLast(ioseed);
2131     }
2132     
2133   } // loop over entries in TreeH
2134
2135   return;
2136 }
2137 //-----------------------------------------------------------------------------
2138 void AliTPCtrackerParam::MakeSeedsFromRefs(TTree *ttr,
2139                                            TObjArray &seedArray) const {
2140 //-----------------------------------------------------------------------------
2141 // This function makes the seeds for tracks from the track references
2142 //-----------------------------------------------------------------------------
2143
2144   Double_t xg,yg,zg,px,py,pz,pt;
2145   Int_t label;
2146   Int_t nnn,k;
2147
2148   TClonesArray *tkRefArray = new TClonesArray("AliTrackReference");
2149
2150   TBranch *b =(TBranch*)ttr->GetBranch("TPC");
2151   if(!b) {cerr<<"TPC branch of TreeTR not found"<<endl; return; }
2152   b->SetAddress(&tkRefArray);
2153   Int_t nTkRef = (Int_t)b->GetEntries();
2154   cerr<<"+++\n+++ Number of entries in TreeTR(TPC): "<<nTkRef<<
2155          "\n+++\n\n";
2156
2157   // loop on track references
2158   for(Int_t l=0; l<nTkRef; l++){
2159     if(l%1000==0) cerr<<"  --- Processing primary track "
2160                       <<l<<" of "<<nTkRef<<" ---\r";
2161     if(!b->GetEvent(l)) continue;
2162     nnn = tkRefArray->GetEntriesFast();
2163
2164     if(nnn <= 0) continue;   
2165     for(k=0; k<nnn; k++) {
2166       AliTrackReference *tkref = (AliTrackReference*)tkRefArray->UncheckedAt(k);
2167
2168       xg = tkref->X();
2169       yg = tkref->Y();
2170       zg = tkref->Z();
2171       px = tkref->Px();
2172       py = tkref->Py();
2173       pz = tkref->Pz();
2174       label = tkref->GetTrack();
2175
2176       pt=TMath::Sqrt(px*px+py*py);
2177       // reject hits with Pt<mag*0.45 GeV/c
2178       if(pt<(fBz*0.45)) continue;
2179
2180       // create the seed
2181       AliTPCseedGeant *ioseed = new AliTPCseedGeant(xg,yg,zg,px,py,pz,label);
2182
2183       // reject if not at the inner part of TPC
2184       if(TMath::Abs(ioseed->GetXL()-82.9701) > 0.1) { 
2185         delete ioseed; continue; 
2186       }
2187
2188       // reject tracks which are not in the TPC acceptance
2189       if(!ioseed->InTPCAcceptance()) { 
2190         delete ioseed; continue; 
2191       }
2192
2193       // put seed in array
2194       seedArray.AddLast(ioseed);
2195
2196     }
2197     
2198
2199   } // loop on track references
2200
2201   delete tkRefArray;
2202
2203
2204   return;
2205 }
2206 //-----------------------------------------------------------------------------
2207 void AliTPCtrackerParam::MergeEvents(Int_t evFirst,Int_t evLast) {
2208 //-----------------------------------------------------------------------------
2209 // This function: 1) merges the files from track comparison
2210 //                   (beware: better no more than 100 events per file)
2211 //                2) computes the average TPC efficiencies
2212 //-----------------------------------------------------------------------------
2213
2214   const char *outName="TPCeff.root";
2215
2216   // merge files with tracks
2217   cerr<<" ******** MERGING FILES **********\n\n";
2218
2219   // create the chain for the tree of compared tracks
2220   TChain ch1("cmptrktree");
2221   TChain ch2("cmptrktree");
2222   TChain ch3("cmptrktree");
2223
2224   for(Int_t j=evFirst; j<=evLast; j++) {
2225     cerr<<"Processing event "<<j<<endl;
2226
2227     TString covName("CovMatrix.");
2228     covName+=j;
2229     covName.Append(".root");
2230
2231     if(gSystem->AccessPathName(covName.Data(),kFileExists)) continue;
2232        
2233
2234     if(j<=100) ch1.Add(covName.Data());
2235     if(j>100 && j<=200) ch2.Add(covName.Data());
2236     if(j>200) ch3.Add(covName.Data());
2237
2238   } // loop on events
2239
2240   // merge chain in one file
2241   TFile *covOut=0;
2242   covOut = new TFile("CovMatrix_AllEvts_1.root","recreate");
2243   ch1.Merge(covOut,1000000000);
2244   covOut->Close();
2245   delete covOut;
2246   covOut = new TFile("CovMatrix_AllEvts_2.root","recreate");
2247   ch2.Merge(covOut,1000000000);
2248   covOut->Close();
2249   delete covOut;
2250   covOut = new TFile("CovMatrix_AllEvts_3.root","recreate");
2251   ch3.Merge(covOut,1000000000);
2252   covOut->Close();
2253   delete covOut;
2254
2255
2256   // efficiencies 
2257   cerr<<" ***** EFFICIENCIES ******\n\n";
2258
2259   ReadEffs("TPCeff.1.root");
2260
2261   Int_t n = fEffPi.GetTotPoints();
2262   Double_t *avEffPi = new Double_t[n];
2263   Double_t *avEffKa = new Double_t[n];
2264   Double_t *avEffPr = new Double_t[n];
2265   Double_t *avEffEl = new Double_t[n];
2266   Double_t *avEffMu = new Double_t[n];
2267   Int_t    *evtsPi  = new Int_t[n];
2268   Int_t    *evtsKa  = new Int_t[n];
2269   Int_t    *evtsPr  = new Int_t[n];
2270   Int_t    *evtsEl  = new Int_t[n];
2271   Int_t    *evtsMu  = new Int_t[n];
2272
2273   for(Int_t j=0; j<n; j++) {
2274     avEffPi[j]=avEffKa[j]=avEffPr[j]=avEffEl[j]=avEffMu[j]=0.;
2275     evtsPi[j]=evtsKa[j]=evtsPr[j]=evtsEl[j]=evtsMu[j]=0;
2276   }
2277
2278   for(Int_t j=evFirst; j<=evLast; j++) {
2279     cerr<<"Processing event "<<j<<endl;
2280
2281     TString effName("TPCeff.");
2282     effName+=j;
2283     effName.Append(".root");
2284        
2285     if(gSystem->AccessPathName(effName.Data(),kFileExists)) continue;
2286
2287     ReadEffs(effName.Data());
2288
2289     for(Int_t k=0; k<n; k++) {
2290       if(fEffPi.GetParam(k)>=0.) {avEffPi[k]+=fEffPi.GetParam(k); evtsPi[k]++;}
2291       if(fEffKa.GetParam(k)>=0.) {avEffKa[k]+=fEffKa.GetParam(k); evtsKa[k]++;}
2292       if(fEffPr.GetParam(k)>=0.) {avEffPr[k]+=fEffPr.GetParam(k); evtsPr[k]++;}
2293       if(fEffEl.GetParam(k)>=0.) {avEffEl[k]+=fEffEl.GetParam(k); evtsEl[k]++;}
2294       if(fEffMu.GetParam(k)>=0.) {avEffMu[k]+=fEffMu.GetParam(k); evtsMu[k]++;}
2295     }
2296
2297   } // loop on events
2298
2299   // compute average efficiencies
2300   for(Int_t j=0; j<n; j++) {
2301     if(evtsPi[j]==0) evtsPi[j]++;
2302     fEffPi.SetParam(j,(Double_t)avEffPi[j]/evtsPi[j]);
2303     if(evtsKa[j]==0) evtsKa[j]++;
2304     fEffKa.SetParam(j,(Double_t)avEffKa[j]/evtsKa[j]);
2305     if(evtsPr[j]==0) evtsPr[j]++;
2306     fEffPr.SetParam(j,(Double_t)avEffPr[j]/evtsPr[j]);
2307     if(evtsEl[j]==0) evtsEl[j]++;
2308     fEffEl.SetParam(j,(Double_t)avEffEl[j]/evtsEl[j]);
2309     if(evtsMu[j]==0) evtsMu[j]++;
2310     fEffMu.SetParam(j,(Double_t)avEffMu[j]/evtsMu[j]);
2311   }
2312  
2313   // write efficiencies to a file
2314   WriteEffs(outName);
2315
2316   delete [] avEffPi;
2317   delete [] avEffKa;
2318   delete [] avEffPr;
2319   delete [] avEffEl;
2320   delete [] avEffMu;
2321   delete [] evtsPi;
2322   delete [] evtsKa;
2323   delete [] evtsPr;
2324   delete [] evtsEl;
2325   delete [] evtsMu;
2326
2327   return;
2328 }
2329 //-----------------------------------------------------------------------------
2330 Int_t AliTPCtrackerParam::ReadAllData(const Char_t* inName) {
2331 //-----------------------------------------------------------------------------
2332 // This function reads all parameters from the DB
2333 //-----------------------------------------------------------------------------
2334
2335   if(!ReadEffs(inName)) return 0;
2336   if(!ReadPulls(inName)) return 0;        
2337   if(!ReadRegParams(inName,211)) return 0; 
2338   if(!ReadRegParams(inName,321)) return 0; 
2339   if(!ReadRegParams(inName,2212)) return 0; 
2340   if(!ReadRegParams(inName,11)) return 0; 
2341   if(!ReadRegParams(inName,13)) return 0; 
2342   if(!ReaddEdx(inName,211)) return 0;
2343   if(!ReaddEdx(inName,321)) return 0;
2344   if(!ReaddEdx(inName,2212)) return 0;
2345   if(!ReaddEdx(inName,11)) return 0;
2346   if(!ReaddEdx(inName,13)) return 0;
2347   if(!ReadDBgrid(inName)) return 0;
2348
2349   return 1;
2350 }
2351 //-----------------------------------------------------------------------------
2352 Int_t AliTPCtrackerParam::ReaddEdx(const Char_t* inName,Int_t pdg) {
2353 //-----------------------------------------------------------------------------
2354 // This function reads the dEdx parameters from the DB
2355 //-----------------------------------------------------------------------------
2356
2357   if(gSystem->AccessPathName(inName,kFileExists)) { 
2358     cerr<<"AliTPCtrackerParam::ReaddEdx: "<<inName<<" not found\n"; 
2359     return 0; }
2360
2361   TFile *inFile = TFile::Open(inName);
2362   switch (pdg) {
2363   case 211:
2364     fdEdxMean = (AliTPCkineGrid*)inFile->Get("/dEdx/Pions/dEdxMeanPi");
2365     fdEdxMeanPi.~AliTPCkineGrid();
2366     new(&fdEdxMeanPi) AliTPCkineGrid(*fdEdxMean);
2367     fdEdxRMS = (AliTPCkineGrid*)inFile->Get("/dEdx/Pions/dEdxRMSPi");
2368     fdEdxRMSPi.~AliTPCkineGrid();
2369     new(&fdEdxRMSPi) AliTPCkineGrid(*fdEdxRMS);
2370     break;
2371   case 321:
2372     fdEdxMean = (AliTPCkineGrid*)inFile->Get("/dEdx/Kaons/dEdxMeanKa");
2373     fdEdxMeanKa.~AliTPCkineGrid();
2374     new(&fdEdxMeanKa) AliTPCkineGrid(*fdEdxMean);
2375     fdEdxRMS = (AliTPCkineGrid*)inFile->Get("/dEdx/Kaons/dEdxRMSKa");
2376     fdEdxRMSKa.~AliTPCkineGrid();
2377     new(&fdEdxRMSKa) AliTPCkineGrid(*fdEdxRMS);
2378     break;
2379   case 2212:
2380     fdEdxMean = (AliTPCkineGrid*)inFile->Get("/dEdx/Protons/dEdxMeanPr");
2381     fdEdxMeanPr.~AliTPCkineGrid();
2382     new(&fdEdxMeanPr) AliTPCkineGrid(*fdEdxMean);
2383     fdEdxRMS = (AliTPCkineGrid*)inFile->Get("/dEdx/Protons/dEdxRMSPr");
2384     fdEdxRMSPr.~AliTPCkineGrid();
2385     new(&fdEdxRMSPr) AliTPCkineGrid(*fdEdxRMS);
2386     break;
2387   case 11:
2388     fdEdxMean = (AliTPCkineGrid*)inFile->Get("/dEdx/Electrons/dEdxMeanEl");
2389     fdEdxMeanEl.~AliTPCkineGrid();
2390     new(&fdEdxMeanEl) AliTPCkineGrid(*fdEdxMean);
2391     fdEdxRMS = (AliTPCkineGrid*)inFile->Get("/dEdx/Electrons/dEdxRMSEl");
2392     fdEdxRMSEl.~AliTPCkineGrid();
2393     new(&fdEdxRMSEl) AliTPCkineGrid(*fdEdxRMS);
2394     break;
2395   case 13:
2396     if(fColl==0) {
2397     fdEdxMean = (AliTPCkineGrid*)inFile->Get("/dEdx/Pions/dEdxMeanPi");
2398     fdEdxRMS = (AliTPCkineGrid*)inFile->Get("/dEdx/Pions/dEdxRMSPi");
2399     } else {
2400       fdEdxMean = (AliTPCkineGrid*)inFile->Get("/dEdx/Muons/dEdxMeanMu");
2401       fdEdxRMS = (AliTPCkineGrid*)inFile->Get("/dEdx/Muons/dEdxRMSMu");
2402     }
2403     fdEdxMeanMu.~AliTPCkineGrid();
2404     new(&fdEdxMeanMu) AliTPCkineGrid(*fdEdxMean);
2405     fdEdxRMSMu.~AliTPCkineGrid();
2406     new(&fdEdxRMSMu) AliTPCkineGrid(*fdEdxRMS);
2407     break;
2408   }
2409   inFile->Close();
2410
2411   return 1;
2412 }
2413 //-----------------------------------------------------------------------------
2414 Int_t AliTPCtrackerParam::ReadDBgrid(const Char_t* inName) {
2415 //-----------------------------------------------------------------------------
2416 // This function reads the kine grid from the DB
2417 //-----------------------------------------------------------------------------
2418
2419   if(gSystem->AccessPathName(inName,kFileExists)) {
2420     cerr<<"AliTPCtrackerParam::ReadCovMatrices: "<<inName<<" not found\n"; 
2421     return 0; }
2422
2423   TFile *inFile = TFile::Open(inName);
2424
2425   // first read the DB grid for the different particles
2426   fDBgrid = (AliTPCkineGrid*)inFile->Get("/CovMatrices/Pions/DBgridPi");
2427   fDBgridPi.~AliTPCkineGrid();
2428   new(&fDBgridPi) AliTPCkineGrid(*fDBgrid);
2429   fDBgrid = (AliTPCkineGrid*)inFile->Get("/CovMatrices/Kaons/DBgridKa");
2430   fDBgridKa.~AliTPCkineGrid();
2431   new(&fDBgridKa) AliTPCkineGrid(*fDBgrid);
2432   if(fColl==0) {
2433     fDBgrid = (AliTPCkineGrid*)inFile->Get("/CovMatrices/Pions/DBgridPi");
2434   } else {
2435     fDBgrid = (AliTPCkineGrid*)inFile->Get("/CovMatrices/Protons/DBgridPr");
2436   }
2437   fDBgridPr.~AliTPCkineGrid();
2438   new(&fDBgridPr) AliTPCkineGrid(*fDBgrid);
2439   fDBgrid = (AliTPCkineGrid*)inFile->Get("/CovMatrices/Electrons/DBgridEl");
2440   fDBgridEl.~AliTPCkineGrid();      
2441   new(&fDBgridEl) AliTPCkineGrid(*fDBgrid);
2442   if(fColl==0) {
2443     fDBgrid = (AliTPCkineGrid*)inFile->Get("/CovMatrices/Pions/DBgridPi");
2444   } else {
2445     fDBgrid = (AliTPCkineGrid*)inFile->Get("/CovMatrices/Muons/DBgridMu");
2446   }
2447   fDBgridMu.~AliTPCkineGrid();
2448   new(&fDBgridMu) AliTPCkineGrid(*fDBgrid);
2449
2450   inFile->Close();
2451
2452   return 1;
2453 }
2454 //-----------------------------------------------------------------------------
2455 Int_t AliTPCtrackerParam::ReadEffs(const Char_t* inName) {
2456 //-----------------------------------------------------------------------------
2457 // This function reads the TPC efficiencies from the DB
2458 //-----------------------------------------------------------------------------
2459
2460   if(gSystem->AccessPathName(inName,kFileExists)) {
2461     cerr<<"AliTPCtrackerParam::ReadEffs: "<<inName<<" not found\n"; 
2462     return 0; }
2463
2464   TFile *inFile = TFile::Open(inName);
2465
2466   fEff = (AliTPCkineGrid*)inFile->Get("/Efficiencies/Pions/EffPi");
2467   fEffPi.~AliTPCkineGrid();
2468   new(&fEffPi) AliTPCkineGrid(*fEff);
2469   fEff = (AliTPCkineGrid*)inFile->Get("/Efficiencies/Kaons/EffKa");
2470   fEffKa.~AliTPCkineGrid();
2471   new(&fEffKa) AliTPCkineGrid(*fEff);
2472   fEff = (AliTPCkineGrid*)inFile->Get("/Efficiencies/Protons/EffPr");
2473   fEffPr.~AliTPCkineGrid();
2474   new(&fEffPr) AliTPCkineGrid(*fEff);
2475   fEff = (AliTPCkineGrid*)inFile->Get("/Efficiencies/Electrons/EffEl");
2476   fEffEl.~AliTPCkineGrid();
2477   new(&fEffEl) AliTPCkineGrid(*fEff);
2478   fEff = (AliTPCkineGrid*)inFile->Get("/Efficiencies/Muons/EffMu");
2479   fEffMu.~AliTPCkineGrid();
2480   new(&fEffMu) AliTPCkineGrid(*fEff);
2481
2482   inFile->Close();
2483
2484   return 1;
2485 }
2486 //-----------------------------------------------------------------------------
2487 Int_t AliTPCtrackerParam::ReadPulls(const Char_t* inName) {
2488 //-----------------------------------------------------------------------------
2489 // This function reads the pulls from the DB
2490 //-----------------------------------------------------------------------------
2491
2492   if(gSystem->AccessPathName(inName,kFileExists)) { 
2493     cerr<<"AliTPCtrackerParam::ReadPulls: "<<inName<<" not found\n"; 
2494     return 0; }
2495
2496   TFile *inFile = TFile::Open(inName);
2497
2498   for(Int_t i=0; i<5; i++) {
2499     TString pi("/Pulls/Pions/PullsPi_"); pi+=i;
2500     TString ka("/Pulls/Kaons/PullsKa_"); ka+=i;
2501     TString pr("/Pulls/Protons/PullsPr_"); pr+=i;
2502     TString el("/Pulls/Electrons/PullsEl_"); el+=i;
2503     TString mu("/Pulls/Muons/PullsMu_"); mu+=i;
2504
2505     fPulls = (AliTPCkineGrid*)inFile->Get(pi.Data());
2506     fPullsPi[i].~AliTPCkineGrid();
2507     new(&fPullsPi[i]) AliTPCkineGrid(*fPulls);
2508
2509     fPulls = (AliTPCkineGrid*)inFile->Get(ka.Data());
2510     fPullsKa[i].~AliTPCkineGrid();
2511     new(&fPullsKa[i]) AliTPCkineGrid(*fPulls);
2512
2513     if(fColl==0) {
2514       fPulls = (AliTPCkineGrid*)inFile->Get(pi.Data());
2515     } else {
2516       fPulls = (AliTPCkineGrid*)inFile->Get(pr.Data());
2517     }
2518     fPullsPr[i].~AliTPCkineGrid();
2519     new(&fPullsPr[i]) AliTPCkineGrid(*fPulls);
2520
2521     fPulls = (AliTPCkineGrid*)inFile->Get(el.Data());
2522     fPullsEl[i].~AliTPCkineGrid();
2523     new(&fPullsEl[i]) AliTPCkineGrid(*fPulls);
2524
2525     if(fColl==0) {
2526       fPulls = (AliTPCkineGrid*)inFile->Get(pi.Data());
2527     } else {
2528       fPulls = (AliTPCkineGrid*)inFile->Get(mu.Data());
2529     }
2530     fPullsMu[i].~AliTPCkineGrid();
2531     new(&fPullsMu[i]) AliTPCkineGrid(*fPulls);
2532   }
2533
2534   inFile->Close();
2535
2536   return 1;
2537 }
2538 //-----------------------------------------------------------------------------
2539 Int_t AliTPCtrackerParam::ReadRegParams(const Char_t* inName,Int_t pdg) {
2540 //-----------------------------------------------------------------------------
2541 // This function reads the regularization parameters from the DB
2542 //-----------------------------------------------------------------------------
2543
2544   if(gSystem->AccessPathName(inName,kFileExists)) { 
2545     cerr<<"AliTPCtrackerParam::ReadRegParams: "<<inName<<" not found\n"; 
2546     return 0; }
2547
2548   TFile *inFile = TFile::Open(inName);
2549   switch (pdg) {
2550   case 211:
2551     fRegPar = (TMatrixD*)inFile->Get("/RegParams/Pions/RegPions");
2552     new(&fRegParPi) TMatrixD(*fRegPar);
2553     break;
2554   case 321:
2555     fRegPar = (TMatrixD*)inFile->Get("/RegParams/Kaons/RegKaons");
2556     new(&fRegParKa) TMatrixD(*fRegPar);
2557     break;
2558   case 2212:
2559     if(fColl==0) {
2560       fRegPar = (TMatrixD*)inFile->Get("/RegParams/Pions/RegPions");
2561     } else {
2562       fRegPar = (TMatrixD*)inFile->Get("/RegParams/Protons/RegProtons");
2563     }
2564     new(&fRegParPr) TMatrixD(*fRegPar);
2565     break;
2566   case 11:
2567     fRegPar = (TMatrixD*)inFile->Get("/RegParams/Electrons/RegElectrons");
2568     new(&fRegParEl) TMatrixD(*fRegPar);
2569     break;
2570   case 13:
2571     if(fColl==0) {
2572       fRegPar = (TMatrixD*)inFile->Get("/RegParams/Pions/RegPions");
2573     } else {
2574       fRegPar = (TMatrixD*)inFile->Get("/RegParams/Muons/RegMuons");
2575     }
2576     new(&fRegParMu) TMatrixD(*fRegPar);
2577     break;
2578   }
2579   inFile->Close();
2580
2581   return 1;
2582 }
2583 //-----------------------------------------------------------------------------
2584 void AliTPCtrackerParam::RegularizeCovMatrix(const Char_t *outName,Int_t pdg) {
2585 //-----------------------------------------------------------------------------
2586 // This function regularizes the elements of the covariance matrix
2587 // that show a momentum depence:
2588 // c00, c11, c22, c33, c44, c20, c24, c40, c31 
2589 // The regularization is done separately for pions, kaons, electrons:
2590 // give "Pion","Kaon" or "Electron" as first argument.
2591 //-----------------------------------------------------------------------------
2592
2593   gStyle->SetOptStat(0);
2594   gStyle->SetOptFit(10001);
2595
2596   Int_t thisPdg=211;
2597   const char *part="Pions - ";
2598
2599   InitializeKineGrid("DB");
2600   SetParticle(pdg);
2601   const Int_t kfitbins = fDBgrid->GetBinsPt();
2602   cerr<<" Fit bins:  "<<kfitbins<<endl;
2603
2604   switch (pdg) {
2605   case 211: // pions
2606     thisPdg=211;  
2607     part="Pions - ";
2608     cerr<<" Processing pions ...\n";
2609     break;
2610   case 321: //kaons
2611     thisPdg=321; 
2612     part="Kaons - ";
2613     cerr<<" Processing kaons ...\n";
2614     break;
2615   case 2212: //protons
2616     thisPdg=2212; 
2617     part="Protons - ";
2618     cerr<<" Processing protons ...\n";
2619     break;
2620   case 11: // electrons
2621     thisPdg= 11; 
2622     part="Electrons - ";
2623     cerr<<" Processing electrons ...\n";
2624     break;
2625   case 13: // muons
2626     thisPdg= 13; 
2627     part="Muons - ";
2628     cerr<<" Processing muons ...\n";
2629     break;
2630   }
2631
2632
2633   /*
2634   // create a chain with compared tracks
2635   TChain *cmptrkchain = new ("cmptrktree");
2636   cmptrkchain.Add("CovMatrix_AllEvts.root");
2637   //cmptrkchain.Add("CovMatrix_AllEvts_1.root");
2638   //cmptrkchain.Add("CovMatrix_AllEvts_2.root");
2639   //cmptrkchain.Add("CovMatrix_AllEvts_3.root");
2640   */
2641
2642
2643   TFile *infile = TFile::Open("CovMatrix_AllEvts.root");
2644   TTree *cmptrktree = (TTree*)infile->Get("cmptrktree");
2645
2646   COMPTRACK cmptrk; 
2647   cmptrktree->SetBranchAddress("comptracks",&cmptrk);
2648   Int_t entries = (Int_t)cmptrktree->GetEntries(); 
2649
2650     
2651   Int_t pbin;
2652   Int_t    *n       = new Int_t[kfitbins];
2653   Int_t    *n00     = new Int_t[kfitbins];
2654   Int_t    *n11     = new Int_t[kfitbins];
2655   Int_t    *n20     = new Int_t[kfitbins];
2656   Int_t    *n22     = new Int_t[kfitbins];
2657   Int_t    *n31     = new Int_t[kfitbins];
2658   Int_t    *n33     = new Int_t[kfitbins];
2659   Int_t    *n40     = new Int_t[kfitbins];
2660   Int_t    *n42     = new Int_t[kfitbins];
2661   Int_t    *n44     = new Int_t[kfitbins];
2662   Double_t *p       = new Double_t[kfitbins];
2663   Double_t *ep      = new Double_t[kfitbins];
2664   Double_t *mean00  = new Double_t[kfitbins];
2665   Double_t *mean11  = new Double_t[kfitbins];
2666   Double_t *mean20  = new Double_t[kfitbins];
2667   Double_t *mean22  = new Double_t[kfitbins];
2668   Double_t *mean31  = new Double_t[kfitbins];
2669   Double_t *mean33  = new Double_t[kfitbins];
2670   Double_t *mean40  = new Double_t[kfitbins];
2671   Double_t *mean42  = new Double_t[kfitbins];
2672   Double_t *mean44  = new Double_t[kfitbins];
2673   Double_t *sigma00 = new Double_t[kfitbins];
2674   Double_t *sigma11 = new Double_t[kfitbins];
2675   Double_t *sigma20 = new Double_t[kfitbins];
2676   Double_t *sigma22 = new Double_t[kfitbins];
2677   Double_t *sigma31 = new Double_t[kfitbins];
2678   Double_t *sigma33 = new Double_t[kfitbins];
2679   Double_t *sigma40 = new Double_t[kfitbins];
2680   Double_t *sigma42 = new Double_t[kfitbins];
2681   Double_t *sigma44 = new Double_t[kfitbins];
2682   Double_t *rmean   = new Double_t[kfitbins];
2683   Double_t *rsigma  = new Double_t[kfitbins];
2684   Double_t fitpar[3];
2685
2686   for(Int_t l=0; l<kfitbins; l++) {
2687     n[l]=1;
2688     n00[l]=n11[l]=n20[l]=n22[l]=n31[l]=n33[l]=n40[l]=n42[l]=n44[l]=1;
2689     p[l ]=ep[l]=0.;
2690     mean00[l]=mean11[l]=mean20[l]=mean22[l]=mean31[l]=mean33[l]=mean40[l]=mean42[l]=mean44[l]=0.;
2691     sigma00[l]=sigma11[l]=sigma20[l]=sigma22[l]=sigma31[l]=sigma33[l]=sigma40[l]=sigma42[l]=sigma44[l]=0.;
2692   }
2693
2694   // loop on chain entries for mean 
2695   for(Int_t l=0; l<entries; l++) {
2696     cmptrktree->GetEvent(l);
2697     if(TMath::Abs(cmptrk.pdg)!=thisPdg) continue;
2698     pbin = (Int_t)fDBgrid->GetBin(cmptrk.pt,cmptrk.eta)/fDBgrid->GetBinsEta();
2699     n[pbin]++;
2700     p[pbin]+=cmptrk.p;
2701     mean00[pbin]+=cmptrk.c00;
2702     mean11[pbin]+=cmptrk.c11;
2703     mean20[pbin]+=cmptrk.c20;
2704     mean22[pbin]+=cmptrk.c22;
2705     mean31[pbin]+=cmptrk.c31;
2706     mean33[pbin]+=cmptrk.c33;
2707     mean40[pbin]+=cmptrk.c40;
2708     mean42[pbin]+=cmptrk.c42;
2709     mean44[pbin]+=cmptrk.c44;
2710   } // loop on chain entries
2711
2712   for(Int_t l=0; l<kfitbins; l++) {
2713     p[l]/=n[l];
2714     mean00[l]/=n[l];
2715     mean11[l]/=n[l];
2716     mean20[l]/=n[l];
2717     mean22[l]/=n[l];
2718     mean31[l]/=n[l];
2719     mean33[l]/=n[l];
2720     mean40[l]/=n[l];
2721     mean42[l]/=n[l];
2722     mean44[l]/=n[l];
2723   }
2724
2725   // loop on chain entries for sigma
2726   for(Int_t l=0; l<entries; l++) {
2727     cmptrktree->GetEvent(l);
2728     if(TMath::Abs(cmptrk.pdg)!=thisPdg) continue;
2729     pbin = (Int_t)fDBgrid->GetBin(cmptrk.pt,cmptrk.eta)/fDBgrid->GetBinsEta();
2730     if(TMath::Abs(cmptrk.c00-mean00[pbin])<0.4*mean00[pbin]) { n00[pbin]++;
2731       sigma00[pbin]+=(cmptrk.c00-mean00[pbin])*(cmptrk.c00-mean00[pbin]); }
2732     if(TMath::Abs(cmptrk.c11-mean11[pbin])<0.4*mean11[pbin]) { n11[pbin]++;
2733       sigma11[pbin]+=(cmptrk.c11-mean11[pbin])*(cmptrk.c11-mean11[pbin]); }
2734     if(TMath::Abs(cmptrk.c20-mean20[pbin])<0.4*mean20[pbin]) { n20[pbin]++;
2735       sigma20[pbin]+=(cmptrk.c20-mean20[pbin])*(cmptrk.c20-mean20[pbin]); }
2736     if(TMath::Abs(cmptrk.c22-mean22[pbin])<0.4*mean22[pbin]) { n22[pbin]++;
2737       sigma22[pbin]+=(cmptrk.c22-mean22[pbin])*(cmptrk.c22-mean22[pbin]); }
2738     if(TMath::Abs(cmptrk.c31-mean31[pbin])<-0.4*mean31[pbin]) { n31[pbin]++;
2739       sigma31[pbin]+=(cmptrk.c31-mean31[pbin])*(cmptrk.c31-mean31[pbin]); }
2740     if(TMath::Abs(cmptrk.c33-mean33[pbin])<0.4*mean33[pbin]) { n33[pbin]++;
2741       sigma33[pbin]+=(cmptrk.c33-mean33[pbin])*(cmptrk.c33-mean33[pbin]); }
2742     if(TMath::Abs(cmptrk.c40-mean40[pbin])<0.4*mean40[pbin]) { n40[pbin]++;
2743       sigma40[pbin]+=(cmptrk.c40-mean40[pbin])*(cmptrk.c40-mean40[pbin]); }
2744     if(TMath::Abs(cmptrk.c42-mean42[pbin])<0.4*mean42[pbin]) { n42[pbin]++;
2745       sigma42[pbin]+=(cmptrk.c42-mean42[pbin])*(cmptrk.c42-mean42[pbin]); }
2746     if(TMath::Abs(cmptrk.c44-mean44[pbin])<0.4*mean44[pbin]) { n44[pbin]++;
2747       sigma44[pbin]+=(cmptrk.c44-mean44[pbin])*(cmptrk.c44-mean44[pbin]); }
2748   } // loop on chain entries
2749  
2750   for(Int_t l=0; l<kfitbins; l++) {
2751     sigma00[l] = TMath::Sqrt(sigma00[l]/n00[l]);
2752     sigma11[l] = TMath::Sqrt(sigma11[l]/n11[l]);
2753     sigma20[l] = TMath::Sqrt(sigma20[l]/n20[l]);
2754     sigma22[l] = TMath::Sqrt(sigma22[l]/n22[l]);
2755     sigma31[l] = TMath::Sqrt(sigma31[l]/n31[l]);
2756     sigma33[l] = TMath::Sqrt(sigma33[l]/n33[l]);
2757     sigma40[l] = TMath::Sqrt(sigma40[l]/n40[l]);
2758     sigma42[l] = TMath::Sqrt(sigma42[l]/n42[l]);
2759     sigma44[l] = TMath::Sqrt(sigma44[l]/n44[l]);
2760   }
2761   
2762
2763   // Fit function
2764   TF1 *func = new TF1("RegFunc",RegFunc,0.23,50.,3);
2765   func->SetParNames("A_meas","A_scatt","B");
2766
2767   // line to draw on the plots
2768   TLine *lin = new TLine(-1,1,1.69,1);
2769   lin->SetLineStyle(2);
2770   lin->SetLineWidth(2);
2771
2772   // matrix used to store fit results
2773   TMatrixD fitRes(9,3);
2774
2775   //    --- c00 ---
2776
2777   // create the canvas
2778   TCanvas *canv00 = new TCanvas("canv00","c00",0,0,700,900); 
2779   canv00->Divide(1,2);
2780   // create the graph for cov matrix
2781   TGraphErrors *gr00 = new TGraphErrors(kfitbins,p,mean00,ep,sigma00);
2782   TString title00("C(y,y)"); title00.Prepend(part);
2783   TH2F *frame00 = new TH2F("frame00",title00.Data(),2,0.1,50,2,0,5e-3);
2784   frame00->SetXTitle("p [GeV/c]");
2785   canv00->cd(1);  gPad->SetLogx();
2786   frame00->Draw();
2787   gr00->Draw("P");
2788   // Sets initial values for parameters
2789   func->SetParameters(1.6e-3,1.9e-4,1.5);
2790   // Fit points in range defined by function
2791   gr00->Fit("RegFunc","R,Q");
2792   func->GetParameters(fitpar);
2793   for(Int_t i=0; i<3; i++) fitRes(0,i)=fitpar[i];
2794   for(Int_t l=0; l<kfitbins; l++) {
2795     rmean[l]  = mean00[l]/RegFunc(&p[l],fitpar);
2796     rsigma[l] = sigma00[l]/RegFunc(&p[l],fitpar);
2797   }
2798   // create the graph the regularized cov. matrix
2799   TGraphErrors *gr00reg = new TGraphErrors(kfitbins,p,rmean,ep,rsigma);
2800   TString regtitle00("C(y,y)/(A_meas+A_scatt/p^{B})"); 
2801   regtitle00.Prepend(part);
2802   TH2F *frame00reg = new TH2F("frame00reg",regtitle00.Data(),2,0.1,50,2,0,2);
2803   frame00reg->SetXTitle("p [GeV/c]");
2804   canv00->cd(2); gPad->SetLogx();
2805   frame00reg->Draw();
2806   gr00reg->Draw("P");
2807   lin->Draw("same");
2808
2809
2810   //    --- c11 ---
2811  
2812   // create the canvas
2813   TCanvas *canv11 = new TCanvas("canv11","c11",0,0,700,900); 
2814   canv11->Divide(1,2);
2815   // create the graph for cov matrix
2816   TGraphErrors *gr11 = new TGraphErrors(kfitbins,p,mean11,ep,sigma11);
2817   TString title11("C(z,z)"); title11.Prepend(part);
2818   TH2F *frame11 = new TH2F("frame11",title11.Data(),2,0.1,50,2,0,6e-3);
2819   frame11->SetXTitle("p [GeV/c]");
2820   canv11->cd(1);  gPad->SetLogx();
2821   frame11->Draw();
2822   gr11->Draw("P");
2823   // Sets initial values for parameters
2824   func->SetParameters(1.2e-3,8.1e-4,1.);
2825   // Fit points in range defined by function
2826   gr11->Fit("RegFunc","R,Q");
2827   func->GetParameters(fitpar);
2828   for(Int_t i=0; i<3; i++) fitRes(1,i)=fitpar[i];
2829   for(Int_t l=0; l<kfitbins; l++) {
2830     rmean[l]  = mean11[l]/RegFunc(&p[l],fitpar);
2831     rsigma[l] = sigma11[l]/RegFunc(&p[l],fitpar);
2832   }
2833   // create the graph the regularized cov. matrix
2834   TGraphErrors *gr11reg = new TGraphErrors(kfitbins,p,rmean,ep,rsigma);
2835   TString regtitle11("C(z,z)/(A_meas+A_scatt/p^{B})"); 
2836   regtitle11.Prepend(part);
2837   TH2F *frame11reg = new TH2F("frame11reg",regtitle11.Data(),2,0.1,50,2,0,2);
2838   frame11reg->SetXTitle("p [GeV/c]");
2839   canv11->cd(2); gPad->SetLogx();
2840   frame11reg->Draw();
2841   gr11reg->Draw("P");
2842   lin->Draw("same");
2843
2844
2845   //    --- c20 ---
2846
2847   // create the canvas
2848   TCanvas *canv20 = new TCanvas("canv20","c20",0,0,700,900); 
2849   canv20->Divide(1,2);
2850   // create the graph for cov matrix
2851   TGraphErrors *gr20 = new TGraphErrors(kfitbins,p,mean20,ep,sigma20);
2852   TString title20("C(#eta, y)"); title20.Prepend(part);
2853   TH2F *frame20 = new TH2F("frame20",title20.Data(),2,0.1,50,2,0,2.5e-4);
2854   frame20->SetXTitle("p [GeV/c]");
2855   canv20->cd(1);  gPad->SetLogx();
2856   frame20->Draw();
2857   gr20->Draw("P");
2858   // Sets initial values for parameters
2859   func->SetParameters(7.3e-5,1.2e-5,1.5);
2860   // Fit points in range defined by function
2861   gr20->Fit("RegFunc","R,Q");
2862   func->GetParameters(fitpar);
2863   for(Int_t i=0; i<3; i++) fitRes(2,i)=fitpar[i];
2864   for(Int_t l=0; l<kfitbins; l++) {
2865     rmean[l]  = mean20[l]/RegFunc(&p[l],fitpar);
2866     rsigma[l] = sigma20[l]/RegFunc(&p[l],fitpar);
2867   }
2868   // create the graph the regularized cov. matrix
2869   TGraphErrors *gr20reg = new TGraphErrors(kfitbins,p,rmean,ep,rsigma);
2870   TString regtitle20("C(#eta, y)/(A_meas+A_scatt/p^{B})"); 
2871   regtitle20.Prepend(part);
2872   TH2F *frame20reg = new TH2F("frame20reg",regtitle20.Data(),2,0.1,50,2,0,2);
2873   frame20reg->SetXTitle("p [GeV/c]");
2874   canv20->cd(2); gPad->SetLogx();
2875   frame20reg->Draw();
2876   gr20reg->Draw("P");
2877   lin->Draw("same");
2878
2879
2880   //    --- c22 ---
2881
2882   // create the canvas
2883   TCanvas *canv22 = new TCanvas("canv22","c22",0,0,700,900); 
2884   canv22->Divide(1,2);
2885   // create the graph for cov matrix
2886   TGraphErrors *gr22 = new TGraphErrors(kfitbins,p,mean22,ep,sigma22);
2887   TString title22("C(#eta, #eta)"); title22.Prepend(part);
2888   TH2F *frame22 = new TH2F("frame22",title22.Data(),2,0.1,50,2,0,3e-5);
2889   frame22->SetXTitle("p [GeV/c]");
2890   canv22->cd(1);  gPad->SetLogx();
2891   frame22->Draw();
2892   gr22->Draw("P");
2893   // Sets initial values for parameters
2894   func->SetParameters(5.2e-6,1.1e-6,2.);
2895   // Fit points in range defined by function
2896   gr22->Fit("RegFunc","R,Q");
2897   func->GetParameters(fitpar);
2898   for(Int_t i=0; i<3; i++) fitRes(3,i)=fitpar[i];
2899   for(Int_t l=0; l<kfitbins; l++) {
2900     rmean[l]  = mean22[l]/RegFunc(&p[l],fitpar);
2901     rsigma[l] = sigma22[l]/RegFunc(&p[l],fitpar);
2902   }
2903   // create the graph the regularized cov. matrix
2904   TGraphErrors *gr22reg = new TGraphErrors(kfitbins,p,rmean,ep,rsigma);
2905   TString regtitle22("C(#eta, #eta)/(A_meas+A_scatt/p^{B})"); 
2906   regtitle22.Prepend(part);
2907   TH2F *frame22reg = new TH2F("frame22reg",regtitle22.Data(),2,0.1,50,2,0,2);
2908   frame22reg->SetXTitle("p [GeV/c]");
2909   canv22->cd(2); gPad->SetLogx();
2910   frame22reg->Draw();
2911   gr22reg->Draw("P");
2912   lin->Draw("same");
2913
2914
2915   //    --- c31 ---
2916
2917   // create the canvas
2918   TCanvas *canv31 = new TCanvas("canv31","c31",0,0,700,900); 
2919   canv31->Divide(1,2);
2920   // create the graph for cov matrix
2921   TGraphErrors *gr31 = new TGraphErrors(kfitbins,p,mean31,ep,sigma31);
2922   TString title31("C(tg #lambda,z)"); title31.Prepend(part);
2923   TH2F *frame31 = new TH2F("frame31",title31.Data(),2,0.1,50,2,-2e-4,0);
2924   frame31->SetXTitle("p [GeV/c]");
2925   canv31->cd(1);  gPad->SetLogx();
2926   frame31->Draw();
2927   gr31->Draw("P");
2928   // Sets initial values for parameters
2929   func->SetParameters(-1.2e-5,-1.2e-5,1.5);
2930   // Fit points in range defined by function
2931   gr31->Fit("RegFunc","R,Q");
2932   func->GetParameters(fitpar);
2933   for(Int_t i=0; i<3; i++) fitRes(4,i)=fitpar[i];
2934   for(Int_t l=0; l<kfitbins; l++) {
2935     rmean[l]  = mean31[l]/RegFunc(&p[l],fitpar);
2936     rsigma[l] = -sigma31[l]/RegFunc(&p[l],fitpar);
2937   }
2938   // create the graph the regularized cov. matrix
2939   TGraphErrors *gr31reg = new TGraphErrors(kfitbins,p,rmean,ep,rsigma);
2940   TString regtitle31("C(tg #lambda,z)/(A_meas+A_scatt/p^{B})"); 
2941   regtitle31.Prepend(part);
2942   TH2F *frame31reg = new TH2F("frame31reg",regtitle31.Data(),2,0.1,50,2,0,2);
2943   frame31reg->SetXTitle("p [GeV/c]");
2944   canv31->cd(2); gPad->SetLogx();
2945   frame31reg->Draw();
2946   gr31reg->Draw("P");
2947   lin->Draw("same");
2948
2949
2950   //    --- c33 ---
2951
2952   // create the canvas
2953   TCanvas *canv33 = new TCanvas("canv33","c33",0,0,700,900); 
2954   canv33->Divide(1,2);
2955   // create the graph for cov matrix
2956   TGraphErrors *gr33 = new TGraphErrors(kfitbins,p,mean33,ep,sigma33);
2957   TString title33("C(tg #lambda,tg #lambda)"); title33.Prepend(part);
2958   TH2F *frame33 = new TH2F("frame33",title33.Data(),2,0.1,50,2,0,1e-5);
2959   frame33->SetXTitle("p [GeV/c]");
2960   canv33->cd(1);  gPad->SetLogx();
2961   frame33->Draw();
2962   gr33->Draw("P");
2963   // Sets initial values for parameters
2964   func->SetParameters(1.3e-7,4.6e-7,1.7);
2965   // Fit points in range defined by function
2966   gr33->Fit("RegFunc","R,Q");
2967   func->GetParameters(fitpar);
2968   for(Int_t i=0; i<3; i++) fitRes(5,i)=fitpar[i];
2969   for(Int_t l=0; l<kfitbins; l++) {
2970     rmean[l]  = mean33[l]/RegFunc(&p[l],fitpar);
2971     rsigma[l] = sigma33[l]/RegFunc(&p[l],fitpar);
2972   }
2973   // create the graph the regularized cov. matrix
2974   TGraphErrors *gr33reg = new TGraphErrors(kfitbins,p,rmean,ep,rsigma);
2975   TString regtitle33("C(tg #lambda,tg #lambda)/(A_meas+A_scatt/p^{B})"); 
2976   regtitle33.Prepend(part);
2977   TH2F *frame33reg = new TH2F("frame33reg",regtitle33.Data(),2,0.1,50,2,0,2);
2978   frame33reg->SetXTitle("p [GeV/c]");
2979   canv33->cd(2); gPad->SetLogx();
2980   frame33reg->Draw();
2981   gr33reg->Draw("P");
2982   lin->Draw("same");
2983
2984
2985   //    --- c40 ---
2986
2987   // create the canvas
2988   TCanvas *canv40 = new TCanvas("canv40","c40",0,0,700,900); 
2989   canv40->Divide(1,2);
2990   // create the graph for cov matrix
2991   TGraphErrors *gr40 = new TGraphErrors(kfitbins,p,mean40,ep,sigma40);
2992   TString title40("C(C,y)"); title40.Prepend(part);
2993   TH2F *frame40 = new TH2F("frame40",title40.Data(),2,0.1,50,2,0,1e-6);
2994   frame40->SetXTitle("p [GeV/c]");
2995   canv40->cd(1);  gPad->SetLogx();
2996   frame40->Draw();
2997   gr40->Draw("P");
2998   // Sets initial values for parameters
2999   func->SetParameters(4.e-7,4.4e-8,1.5);
3000   // Fit points in range defined by function
3001   gr40->Fit("RegFunc","R,Q");
3002   func->GetParameters(fitpar);
3003   for(Int_t i=0; i<3; i++) fitRes(6,i)=fitpar[i];
3004   for(Int_t l=0; l<kfitbins; l++) {
3005     rmean[l]  = mean40[l]/RegFunc(&p[l],fitpar);
3006     rsigma[l] = sigma40[l]/RegFunc(&p[l],fitpar);
3007   }
3008   // create the graph the regularized cov. matrix
3009   TGraphErrors *gr40reg = new TGraphErrors(kfitbins,p,rmean,ep,rsigma);
3010   TString regtitle40("C(C,y)/(A_meas+A_scatt/p^{B})"); 
3011   regtitle40.Prepend(part);
3012   TH2F *frame40reg = new TH2F("frame40reg",regtitle40.Data(),2,0.1,50,2,0,2);
3013   frame40reg->SetXTitle("p [GeV/c]");
3014   canv40->cd(2); gPad->SetLogx();
3015   frame40reg->Draw();
3016   gr40reg->Draw("P");
3017   lin->Draw("same");
3018
3019
3020   //    --- c42 ---
3021
3022   // create the canvas
3023   TCanvas *canv42 = new TCanvas("canv42","c42",0,0,700,900); 
3024   canv42->Divide(1,2);
3025   // create the graph for cov matrix
3026   TGraphErrors *gr42 = new TGraphErrors(kfitbins,p,mean42,ep,sigma42);
3027   TString title42("C(C, #eta)"); title42.Prepend(part);
3028   TH2F *frame42 = new TH2F("frame42",title42.Data(),2,0.1,50,2,0,2.2e-7);
3029   frame42->SetXTitle("p [GeV/c]");
3030   canv42->cd(1);  gPad->SetLogx();
3031   frame42->Draw();
3032   gr42->Draw("P");
3033   // Sets initial values for parameters
3034   func->SetParameters(3.e-8,8.2e-9,2.);
3035   // Fit points in range defined by function
3036   gr42->Fit("RegFunc","R,Q");
3037   func->GetParameters(fitpar);
3038   for(Int_t i=0; i<3; i++) fitRes(7,i)=fitpar[i];
3039   for(Int_t l=0; l<kfitbins; l++) {
3040     rmean[l]  = mean42[l]/RegFunc(&p[l],fitpar);
3041     rsigma[l] = sigma42[l]/RegFunc(&p[l],fitpar);
3042   }
3043   // create the graph the regularized cov. matrix
3044   TGraphErrors *gr42reg = new TGraphErrors(kfitbins,p,rmean,ep,rsigma);
3045   TString regtitle42("C(C, #eta)/(A_meas+A_scatt/p^{B})"); 
3046   regtitle42.Prepend(part);
3047   TH2F *frame42reg = new TH2F("frame42reg",regtitle42.Data(),2,0.1,50,2,0,2);
3048   frame42reg->SetXTitle("p [GeV/c]");
3049   canv42->cd(2); gPad->SetLogx();
3050   frame42reg->Draw();
3051   gr42reg->Draw("P");
3052   lin->Draw("same");
3053
3054
3055   //    --- c44 ---
3056
3057   // create the canvas
3058   TCanvas *canv44 = new TCanvas("canv44","c44",0,0,700,900); 
3059   canv44->Divide(1,2);
3060   // create the graph for cov matrix
3061   TGraphErrors *gr44 = new TGraphErrors(kfitbins,p,mean44,ep,sigma44);
3062   TString title44("C(C,C)"); title44.Prepend(part);
3063   TH2F *frame44 = new TH2F("frame44",title44.Data(),2,0.1,50,2,0,2e-9);
3064   frame44->SetXTitle("p [GeV/c]");
3065   canv44->cd(1);  gPad->SetLogx();
3066   frame44->Draw();
3067   gr44->Draw("P");
3068   // Sets initial values for parameters
3069   func->SetParameters(1.8e-10,5.8e-11,2.);
3070   // Fit points in range defined by function
3071   gr44->Fit("RegFunc","R,Q");
3072   func->GetParameters(fitpar);
3073   for(Int_t i=0; i<3; i++) fitRes(8,i)=fitpar[i];
3074   for(Int_t l=0; l<kfitbins; l++) {
3075     rmean[l]  = mean44[l]/RegFunc(&p[l],fitpar);
3076     rsigma[l] = sigma44[l]/RegFunc(&p[l],fitpar);
3077   }
3078   // create the graph the regularized cov. matrix
3079   TGraphErrors *gr44reg = new TGraphErrors(kfitbins,p,rmean,ep,rsigma);
3080   TString regtitle44("C(C,C)/(A_meas+A_scatt/p^{B})"); 
3081   regtitle44.Prepend(part);
3082   TH2F *frame44reg = new TH2F("frame44reg",regtitle44.Data(),2,0.1,50,2,0,2);
3083   frame44reg->SetXTitle("p [GeV/c]");
3084   canv44->cd(2); gPad->SetLogx();
3085   frame44reg->Draw();
3086   gr44reg->Draw("P");
3087   lin->Draw("same");
3088
3089
3090
3091
3092   switch (pdg) {
3093   case 211:
3094     new(&fRegParPi) TMatrixD(fitRes);
3095     break;
3096   case 321:
3097     new(&fRegParKa) TMatrixD(fitRes);
3098     break;
3099   case 2212:
3100     new(&fRegParPr) TMatrixD(fitRes);
3101     break;
3102   case 11:
3103     new(&fRegParEl) TMatrixD(fitRes);
3104     break;
3105   case 13:
3106     new(&fRegParMu) TMatrixD(fitRes);
3107     break;
3108   }
3109
3110   // write fit parameters to file
3111   WriteRegParams(outName,pdg);
3112
3113   delete [] n;
3114   delete [] n00;
3115   delete [] n11;
3116   delete [] n20;
3117   delete [] n22;
3118   delete [] n31;
3119   delete [] n33;
3120   delete [] n40;
3121   delete [] n42;
3122   delete [] n44;
3123   delete [] p;
3124   delete [] ep;
3125   delete [] mean00;
3126   delete [] mean11;
3127   delete [] mean20;
3128   delete [] mean22;
3129   delete [] mean31;
3130   delete [] mean33;
3131   delete [] mean40;
3132   delete [] mean42;
3133   delete [] mean44;
3134   delete [] sigma00;
3135   delete [] sigma11;
3136   delete [] sigma20;
3137   delete [] sigma22;
3138   delete [] sigma31;
3139   delete [] sigma33;
3140   delete [] sigma40;
3141   delete [] sigma42;
3142   delete [] sigma44;
3143   delete [] rmean;
3144   delete [] rsigma;
3145
3146   return;
3147 }
3148 //-----------------------------------------------------------------------------
3149 Bool_t AliTPCtrackerParam::SelectedTrack(Double_t pt,Double_t eta) const {
3150 //-----------------------------------------------------------------------------
3151 // This function makes a selection according to TPC tracking efficiency
3152 //-----------------------------------------------------------------------------
3153
3154   Double_t eff=0.;
3155  
3156   eff = fEff->GetValueAt(pt,eta);
3157
3158   if(gRandom->Rndm() < eff) return kTRUE;
3159
3160   return kFALSE;
3161 }
3162 //-----------------------------------------------------------------------------
3163 void AliTPCtrackerParam::SetParticle(Int_t pdg) {
3164 //-----------------------------------------------------------------------------
3165 // This function set efficiencies, pulls, etc... for the particle
3166 // specie of the current particle
3167 //-----------------------------------------------------------------------------
3168
3169   switch (pdg) {
3170   case 211:
3171     fDBgrid = &fDBgridPi;
3172     fEff    = &fEffPi;
3173     fPulls  =  fPullsPi;
3174     fRegPar = &fRegParPi;
3175     fdEdxMean = &fdEdxMeanPi;
3176     fdEdxRMS  = &fdEdxRMSPi;
3177     break;
3178   case 321:
3179     fDBgrid = &fDBgridKa;
3180     fEff    = &fEffKa;
3181     fPulls  =  fPullsKa;
3182     fRegPar = &fRegParKa;
3183     fdEdxMean = &fdEdxMeanKa;
3184     fdEdxRMS  = &fdEdxRMSKa;
3185     break;
3186   case 2212:
3187     fDBgrid = &fDBgridPr;
3188     fEff    = &fEffPr;
3189     fPulls  =  fPullsPr;
3190     fRegPar = &fRegParPr;
3191     fdEdxMean = &fdEdxMeanPr;
3192     fdEdxRMS  = &fdEdxRMSPr;
3193     break;
3194   case 11:
3195     fDBgrid = &fDBgridEl;
3196     fEff    = &fEffEl;
3197     fPulls  =  fPullsEl;
3198     fRegPar = &fRegParEl;
3199     fdEdxMean = &fdEdxMeanEl;
3200     fdEdxRMS  = &fdEdxRMSEl;
3201     break;
3202   case 13:
3203     fDBgrid = &fDBgridMu;
3204     fEff    = &fEffMu;
3205     fPulls  =  fPullsMu;
3206     fRegPar = &fRegParMu;
3207     fdEdxMean = &fdEdxMeanMu;
3208     fdEdxRMS  = &fdEdxRMSMu;
3209     break;
3210   }
3211
3212   return;
3213 }
3214 //-----------------------------------------------------------------------------
3215 void AliTPCtrackerParam::SmearTrack(Double_t *xx,Double_t *xxsm,TMatrixD cov)
3216   const {
3217 //-----------------------------------------------------------------------------
3218 // This function smears track parameters according to streched cov. matrix
3219 //-----------------------------------------------------------------------------
3220   AliGausCorr *corgen = new AliGausCorr(cov,5);
3221   TArrayD corr(5);
3222   corgen->GetGaussN(corr);
3223   delete corgen;
3224   corgen = 0;
3225
3226   for(Int_t l=0;l<5;l++) {
3227     xxsm[l] = xx[l]+corr[l];
3228   }
3229
3230   return;
3231 }
3232 //-----------------------------------------------------------------------------
3233 Int_t AliTPCtrackerParam::WritedEdx(const Char_t *outName,Int_t pdg) {
3234 //-----------------------------------------------------------------------------
3235 // This function writes the dEdx parameters to the DB
3236 //-----------------------------------------------------------------------------
3237
3238   Option_t *opt;
3239   const char *dirName="Pions";
3240   const char *meanName="dEdxMeanPi";
3241   const char *rmsName="dEdxRMSPi";
3242
3243   SetParticle(pdg);
3244
3245   if(gSystem->AccessPathName(outName,kFileExists))  { opt="recreate"; 
3246   } else { opt="update"; }
3247
3248   switch (pdg) {
3249   case 211:
3250     dirName="Pions";
3251     meanName="dEdxMeanPi";
3252     rmsName="dEdxRMSPi";
3253     break;
3254   case 321:
3255     dirName="Kaons";
3256     meanName="dEdxMeanKa";
3257     rmsName="dEdxRMSKa";
3258     break;
3259   case 2212:
3260     dirName="Protons";
3261     meanName="dEdxMeanPr";
3262     rmsName="dEdxRMSPr";
3263     break;
3264   case 11:
3265     dirName="Electrons";
3266     meanName="dEdxMeanEl";
3267     rmsName="dEdxRMSEl";
3268     break;
3269   case 13:
3270     dirName="Muons";
3271     meanName="dEdxMeanMu";
3272     rmsName="dEdxRMSMu";
3273     break;
3274   }
3275
3276   TFile *outFile = new TFile(outName,opt);
3277   if(!gDirectory->cd("/dEdx")) {
3278     outFile->mkdir("dEdx");
3279     gDirectory->cd("/dEdx"); 
3280   }
3281   TDirectory *dir2 = gDirectory->mkdir(dirName);
3282   dir2->cd();
3283   fdEdxMean->SetName(meanName); fdEdxMean->Write();
3284   fdEdxRMS->SetName(rmsName);  fdEdxRMS->Write();
3285
3286   outFile->Close();
3287   delete outFile;
3288
3289
3290   return 1;
3291 }
3292 //-----------------------------------------------------------------------------
3293 Int_t AliTPCtrackerParam::WriteEffs(const Char_t *outName) {
3294 //-----------------------------------------------------------------------------
3295 // This function writes the TPC efficiencies to the DB
3296 //-----------------------------------------------------------------------------
3297
3298
3299   Option_t *opt;
3300
3301   if(gSystem->AccessPathName(outName,kFileExists))  { opt="recreate"; 
3302   } else { opt="update"; }
3303
3304   TFile *outFile = new TFile(outName,opt);
3305
3306   outFile->mkdir("Efficiencies");
3307   gDirectory->cd("/Efficiencies");
3308   gDirectory->mkdir("Pions");
3309   gDirectory->mkdir("Kaons");
3310   gDirectory->mkdir("Protons");
3311   gDirectory->mkdir("Electrons");
3312   gDirectory->mkdir("Muons");
3313
3314   gDirectory->cd("/Efficiencies/Pions");
3315   fEffPi.SetName("EffPi");
3316   fEffPi.Write();
3317   gDirectory->cd("/Efficiencies/Kaons");
3318   fEffKa.SetName("EffKa");
3319   fEffKa.Write();
3320   gDirectory->cd("/Efficiencies/Protons");
3321   fEffPr.SetName("EffPr");
3322   fEffPr.Write();
3323   gDirectory->cd("/Efficiencies/Electrons");
3324   fEffEl.SetName("EffEl");
3325   fEffEl.Write();
3326   gDirectory->cd("/Efficiencies/Muons");
3327   fEffMu.SetName("EffMu");
3328   fEffMu.Write();
3329
3330   outFile->Close();
3331
3332   delete outFile;
3333
3334   return 1;
3335 }
3336 //-----------------------------------------------------------------------------
3337 Int_t AliTPCtrackerParam::WritePulls(const Char_t *outName) {
3338 //-----------------------------------------------------------------------------
3339 // This function writes the pulls to the DB
3340 //-----------------------------------------------------------------------------
3341
3342   Option_t *opt;
3343
3344   if(gSystem->AccessPathName(outName,kFileExists))  { opt="recreate"; 
3345   } else { opt="update"; }
3346
3347   TFile *outFile = new TFile(outName,opt);
3348
3349   outFile->mkdir("Pulls");
3350   gDirectory->cd("/Pulls");
3351   gDirectory->mkdir("Pions");
3352   gDirectory->mkdir("Kaons");
3353   gDirectory->mkdir("Protons");
3354   gDirectory->mkdir("Electrons");
3355   gDirectory->mkdir("Muons");
3356
3357   for(Int_t i=0;i<5;i++) {
3358     TString pi("PullsPi_"); pi+=i;
3359     TString ka("PullsKa_"); ka+=i;
3360     TString pr("PullsPr_"); pr+=i;
3361     TString el("PullsEl_"); el+=i;
3362     TString mu("PullsMu_"); mu+=i;
3363     fPullsPi[i].SetName(pi.Data());
3364     fPullsKa[i].SetName(ka.Data());
3365     fPullsPr[i].SetName(pr.Data());
3366     fPullsEl[i].SetName(el.Data());
3367     fPullsMu[i].SetName(mu.Data());
3368     gDirectory->cd("/Pulls/Pions");
3369     fPullsPi[i].Write();
3370     gDirectory->cd("/Pulls/Kaons");
3371     fPullsKa[i].Write();
3372     gDirectory->cd("/Pulls/Protons");
3373     fPullsPr[i].Write();
3374     gDirectory->cd("/Pulls/Electrons");
3375     fPullsEl[i].Write();
3376     gDirectory->cd("/Pulls/Muons");
3377     fPullsMu[i].Write();
3378   }
3379   outFile->Close();
3380   delete outFile;
3381
3382   return 1;
3383 }
3384 //-----------------------------------------------------------------------------
3385 Int_t AliTPCtrackerParam::WriteRegParams(const Char_t *outName,Int_t pdg) {
3386 //-----------------------------------------------------------------------------
3387 // This function writes the regularization parameters to the DB
3388 //-----------------------------------------------------------------------------
3389
3390   Option_t *opt;
3391   const char *dirName="Pions";
3392   const char *keyName="RegPions";
3393
3394   SetParticle(pdg);
3395
3396   if(gSystem->AccessPathName(outName,kFileExists))  { opt="recreate"; 
3397   } else { opt="update"; }
3398
3399   switch (pdg) {
3400   case 211:
3401     dirName="Pions";
3402     keyName="RegPions";
3403     break;
3404   case 321:
3405     dirName="Kaons";
3406     keyName="RegKaons";
3407     break;
3408   case 2212:
3409     dirName="Protons";
3410     keyName="RegProtons";
3411     break;
3412   case 11:
3413     dirName="Electrons";
3414     keyName="RegElectrons";
3415     break;
3416   case 13:
3417     dirName="Muons";
3418     keyName="RegMuons";
3419     break;
3420   }
3421
3422   TFile *outFile = new TFile(outName,opt);
3423   if(!gDirectory->cd("/RegParams")) {
3424     outFile->mkdir("RegParams");
3425     gDirectory->cd("/RegParams"); 
3426   }
3427   TDirectory *dir2 = gDirectory->mkdir(dirName);
3428   dir2->cd();
3429   fRegPar->Write(keyName);
3430
3431   outFile->Close();
3432   delete outFile;
3433
3434
3435   return 1;
3436 }
3437 //-----------------------------------------------------------------------------
3438 //*****************************************************************************
3439 //-----------------------------------------------------------------------------
3440
3441
3442
3443
3444
3445
3446
3447
3448
3449
3450
3451
3452
3453
3454
3455
3456
3457
3458
3459
3460
3461
3462