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