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