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