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