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