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