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