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