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