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