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