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