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