]>
Commit | Line | Data |
---|---|---|
f7f33dec | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
54b76c13 | 16 | /* |
108953e9 | 17 | Comments to be written here: |
54b76c13 | 18 | 1. What do we calibrate. |
19 | 2. How to interpret results | |
20 | 3. Simple example | |
21 | 4. Analysis using debug streamers. | |
22 | ||
23 | ||
24 | ||
25 | 3.Simple example | |
26 | // To make cosmic scan the user interaction neccessary | |
27 | // | |
28 | .x ~/UliStyle.C | |
29 | gSystem->Load("libANALYSIS"); | |
30 | gSystem->Load("libTPCcalib"); | |
31 | TFile fcalib("CalibObjects.root"); | |
32 | TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib"); | |
33 | AliTPCcalibCosmic * cosmic = ( AliTPCcalibCosmic *)array->FindObject("cosmicTPC"); | |
34 | ||
35 | ||
3e55050f | 36 | // |
37 | // | |
38 | gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros"); | |
39 | gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+") | |
40 | AliXRDPROOFtoolkit tool; | |
a390f11f | 41 | TChain * chainCosmic = tool.MakeChainRandom("cosmic.txt","Track0",0,10000); |
42 | TChain * chainBudget = tool.MakeChainRandom("cosmic.txt","material",0,1000); | |
43 | TCut cutptV="abs(1/pt0V-1/pt1V)<0.1"; | |
44 | TCut cutptI="abs(1/pt0In-1/pt1In)<0.5"; | |
45 | TCut cutncl="nclmin>120"; | |
46 | TCut cutDz="abs(p0.fP[1])<50"; | |
47 | TCut cutDr="abs(p0.fP[0])<50"; | |
48 | // | |
49 | chainBudget->Draw(">>listB",cutptV+cutptI+cutncl+cutDr+cutDz,"entryList"); | |
50 | TEntryList *elistB = (TEntryList*)gDirectory->Get("listB"); | |
51 | chainBudget->SetEntryList(elistB); | |
52 | ||
53 | chainBudget->SetAlias("dptrel","(pt0V-pt1V)/((pt0V+pt1V)*0.5)"); | |
54 | chainBudget->SetAlias("dptInrel","(pt0In-pt1In)/((pt0In+pt1In)*0.5)"); | |
55 | chainBudget->SetAlias("ptcorr","(pt0In-pt0V)/(pt0V)+(pt1V-pt1In)/(pt1In)"); | |
54b76c13 | 56 | */ |
57 | ||
58 | ||
59 | ||
f7f33dec | 60 | #include "Riostream.h" |
61 | #include "TChain.h" | |
62 | #include "TTree.h" | |
63 | #include "TH1F.h" | |
64 | #include "TH2F.h" | |
65 | #include "TList.h" | |
66 | #include "TMath.h" | |
67 | #include "TCanvas.h" | |
68 | #include "TFile.h" | |
54b76c13 | 69 | #include "TF1.h" |
91fd44c9 | 70 | #include "THnSparse.h" |
9963b5e2 | 71 | #include "TDatabasePDG.h" |
f7f33dec | 72 | |
54b76c13 | 73 | #include "AliTPCclusterMI.h" |
f7f33dec | 74 | #include "AliTPCseed.h" |
75 | #include "AliESDVertex.h" | |
76 | #include "AliESDEvent.h" | |
77 | #include "AliESDfriend.h" | |
78 | #include "AliESDInputHandler.h" | |
54b76c13 | 79 | #include "AliAnalysisManager.h" |
f7f33dec | 80 | |
81 | #include "AliTracker.h" | |
f7a1cc68 | 82 | #include "AliMagF.h" |
54b76c13 | 83 | #include "AliTPCCalROC.h" |
f7f33dec | 84 | |
85 | #include "AliLog.h" | |
86 | ||
87 | #include "AliTPCcalibCosmic.h" | |
f7f33dec | 88 | #include "TTreeStream.h" |
89 | #include "AliTPCTracklet.h" | |
3326b323 | 90 | //#include "AliESDcosmic.h" |
15e48021 | 91 | |
f7f33dec | 92 | |
93 | ClassImp(AliTPCcalibCosmic) | |
94 | ||
95 | ||
96 | AliTPCcalibCosmic::AliTPCcalibCosmic() | |
97 | :AliTPCcalibBase(), | |
9b27d39b | 98 | fHistNTracks(0), |
99 | fClusters(0), | |
100 | fModules(0), | |
101 | fHistPt(0), | |
9b27d39b | 102 | fDeDx(0), |
54b76c13 | 103 | fDeDxMIP(0), |
104 | fMIPvalue(1), | |
9b27d39b | 105 | fCutMaxD(5), // maximal distance in rfi ditection |
a6dc0cf6 | 106 | fCutMaxDz(40), // maximal distance in z ditection |
9b27d39b | 107 | fCutTheta(0.03), // maximal distan theta |
108 | fCutMinDir(-0.99) // direction vector products | |
f7f33dec | 109 | { |
54b76c13 | 110 | AliInfo("Default Constructor"); |
91fd44c9 | 111 | for (Int_t ihis=0; ihis<6;ihis++){ |
112 | fHistoDelta[ihis]=0; | |
113 | fHistoPull[ihis]=0; | |
8a92e133 | 114 | } |
115 | for (Int_t ihis=0; ihis<4;ihis++){ | |
116 | fHistodEdxMax[ihis] =0; | |
117 | fHistodEdxTot[ihis] =0; | |
91fd44c9 | 118 | } |
f7f33dec | 119 | } |
120 | ||
121 | ||
122 | AliTPCcalibCosmic::AliTPCcalibCosmic(const Text_t *name, const Text_t *title) | |
123 | :AliTPCcalibBase(), | |
9b27d39b | 124 | fHistNTracks(0), |
125 | fClusters(0), | |
126 | fModules(0), | |
127 | fHistPt(0), | |
9b27d39b | 128 | fDeDx(0), |
54b76c13 | 129 | fDeDxMIP(0), |
130 | fMIPvalue(1), | |
a6dc0cf6 | 131 | fCutMaxD(5), // maximal distance in rfi ditection |
132 | fCutMaxDz(40), // maximal distance in z ditection | |
9b27d39b | 133 | fCutTheta(0.03), // maximal distan theta |
134 | fCutMinDir(-0.99) // direction vector products | |
f7f33dec | 135 | { |
136 | SetName(name); | |
137 | SetTitle(title); | |
54b76c13 | 138 | |
139 | fHistNTracks = new TH1F("ntracks","Number of Tracks per Event; number of tracks per event; number of tracks",501,-0.5,500.5); | |
140 | fClusters = new TH1F("signal","Number of Clusters per track; number of clusters per track n_{cl}; counts",160,0,160); | |
141 | fModules = new TH2F("sector","Acorde hits; z (cm); x(cm)",1200,-650,650,600,-700,700); | |
142 | fHistPt = new TH1F("Pt","Pt distribution; p_{T} (GeV); counts",2000,0,50); | |
143 | fDeDx = new TH2F("DeDx","dEdx; momentum p (GeV); TPC signal (a.u.)",500,0.01,100.,500,2.,1000); | |
f7f33dec | 144 | BinLogX(fDeDx); |
54b76c13 | 145 | fDeDxMIP = new TH1F("DeDxMIP","MIP region; TPC signal (a.u.);counts ",500,2.,1000); |
91fd44c9 | 146 | Init(); |
9b27d39b | 147 | AliInfo("Non Default Constructor"); |
bca20570 | 148 | // |
f7f33dec | 149 | } |
150 | ||
151 | AliTPCcalibCosmic::~AliTPCcalibCosmic(){ | |
152 | // | |
153 | // | |
154 | // | |
91fd44c9 | 155 | for (Int_t ihis=0; ihis<6;ihis++){ |
156 | delete fHistoDelta[ihis]; | |
157 | delete fHistoPull[ihis]; | |
91fd44c9 | 158 | } |
8a92e133 | 159 | for (Int_t ihis=0; ihis<4;ihis++){ |
160 | delete fHistodEdxTot[ihis]; | |
161 | delete fHistodEdxMax[ihis]; | |
162 | } | |
163 | ||
164 | delete fHistNTracks; // histogram showing number of ESD tracks per event | |
165 | delete fClusters; // histogram showing the number of clusters per track | |
166 | delete fModules; // 2d histogram of tracks which are propagated to the ACORDE scintillator array | |
167 | delete fHistPt; // Pt histogram of reconstructed tracks | |
168 | delete fDeDx; // dEdx spectrum showing the different particle types | |
169 | delete fDeDxMIP; // TPC signal close to the MIP region of muons 0.4 < p < 0.45 GeV | |
f7f33dec | 170 | } |
171 | ||
172 | ||
91fd44c9 | 173 | void AliTPCcalibCosmic::Init(){ |
174 | // | |
175 | // init component | |
176 | // Make performance histograms | |
177 | // | |
178 | ||
179 | // tracking performance bins | |
180 | // 0 - delta of interest | |
181 | // 1 - min (track0, track1) number of clusters | |
182 | // 2 - R - vertex radius | |
183 | // 3 - P1 - mean z | |
184 | // 4 - P2 - snp(phi) at inner wall of TPC | |
185 | // 5 - P3 - tan(theta) at inner wall of TPC | |
186 | // 6 - P4 - 1/pt mean | |
187 | // 7 - pt - pt mean | |
188 | // 8 - alpha | |
189 | ||
190 | Double_t xminTrack[9], xmaxTrack[9]; | |
191 | Int_t binsTrack[9]; | |
192 | TString axisName[9]; | |
193 | // | |
194 | binsTrack[0] =100; | |
195 | axisName[0] ="#Delta"; | |
196 | // | |
197 | binsTrack[1] =8; | |
198 | xminTrack[1] =80; xmaxTrack[1]=160; | |
199 | axisName[1] ="N_{cl}"; | |
200 | // | |
201 | binsTrack[2] =10; | |
202 | xminTrack[2] =0; xmaxTrack[2]=90; // | |
203 | axisName[2] ="dca_{r} (cm)"; | |
204 | // | |
205 | binsTrack[3] =25; | |
206 | xminTrack[3] =-250; xmaxTrack[3]=250; // | |
207 | axisName[3] ="z (cm)"; | |
208 | // | |
209 | binsTrack[4] =10; | |
210 | xminTrack[4] =-0.8; xmaxTrack[4]=0.8; // | |
211 | axisName[4] ="sin(#phi)"; | |
212 | // | |
213 | binsTrack[5] =10; | |
214 | xminTrack[5] =-1; xmaxTrack[5]=1; // | |
8a92e133 | 215 | axisName[5] ="tan(#theta)"; |
91fd44c9 | 216 | // |
a390f11f | 217 | binsTrack[6] =40; |
218 | xminTrack[6] =-2; xmaxTrack[6]=2; // | |
91fd44c9 | 219 | axisName[6] ="1/pt (1/GeV)"; |
220 | // | |
8a92e133 | 221 | binsTrack[7] =40; |
91fd44c9 | 222 | xminTrack[7] =0.2; xmaxTrack[7]=50; // |
8a92e133 | 223 | axisName[7] ="pt (GeV)"; |
91fd44c9 | 224 | // |
a390f11f | 225 | binsTrack[8] =18; |
91fd44c9 | 226 | xminTrack[8] =0; xmaxTrack[8]=TMath::Pi(); // |
227 | axisName[8] ="alpha"; | |
228 | // | |
229 | // delta y | |
230 | xminTrack[0] =-1; xmaxTrack[0]=1; // | |
231 | fHistoDelta[0] = new THnSparseS("#Delta_{Y} (cm)","#Delta_{Y} (cm)", 9, binsTrack,xminTrack, xmaxTrack); | |
232 | xminTrack[0] =-5; xmaxTrack[0]=5; // | |
233 | fHistoPull[0] = new THnSparseS("#Delta_{Y} (unit)","#Delta_{Y} (unit)", 9, binsTrack,xminTrack, xmaxTrack); | |
234 | // | |
235 | // delta z | |
236 | xminTrack[0] =-1; xmaxTrack[0]=1; // | |
237 | fHistoDelta[1] = new THnSparseS("#Delta_{Z} (cm)","#Delta_{Z} (cm)", 9, binsTrack,xminTrack, xmaxTrack); | |
238 | xminTrack[0] =-5; xmaxTrack[0]=5; // | |
239 | fHistoPull[1] = new THnSparseS("#Delta_{Z} (unit)","#Delta_{Z} (unit)", 9, binsTrack,xminTrack, xmaxTrack); | |
240 | // | |
241 | // delta P2 | |
242 | xminTrack[0] =-10; xmaxTrack[0]=10; // | |
243 | fHistoDelta[2] = new THnSparseS("#Delta_{#phi} (mrad)","#Delta_{#phi} (mrad)", 9, binsTrack,xminTrack, xmaxTrack); | |
244 | xminTrack[0] =-5; xmaxTrack[0]=5; // | |
245 | fHistoPull[2] = new THnSparseS("#Delta_{#phi} (unit)","#Delta_{#phi} (unit)", 9, binsTrack,xminTrack, xmaxTrack); | |
246 | // | |
247 | // delta P3 | |
248 | xminTrack[0] =-10; xmaxTrack[0]=10; // | |
249 | fHistoDelta[3] = new THnSparseS("#Delta_{#theta} (mrad)","#Delta_{#theta} (mrad)", 9, binsTrack,xminTrack, xmaxTrack); | |
250 | xminTrack[0] =-5; xmaxTrack[0]=5; // | |
251 | fHistoPull[3] = new THnSparseS("#Delta_{#theta} (unit)","#Delta_{#theta} (unit)", 9, binsTrack,xminTrack, xmaxTrack); | |
252 | // | |
253 | // delta P4 | |
254 | xminTrack[0] =-0.2; xmaxTrack[0]=0.2; // | |
255 | fHistoDelta[4] = new THnSparseS("#Delta_{1/pt} (1/GeV)","#Delta_{1/pt} (1/GeV)", 9, binsTrack,xminTrack, xmaxTrack); | |
256 | xminTrack[0] =-5; xmaxTrack[0]=5; // | |
257 | fHistoPull[4] = new THnSparseS("#Delta_{1/pt} (unit)","#Delta_{1/pt} (unit)", 9, binsTrack,xminTrack, xmaxTrack); | |
258 | ||
259 | // | |
260 | // delta Pt | |
261 | xminTrack[0] =-0.5; xmaxTrack[0]=0.5; // | |
262 | fHistoDelta[5] = new THnSparseS("#Delta_{pt}/p_{t}","#Delta_{pt}/p_{t}", 9, binsTrack,xminTrack, xmaxTrack); | |
263 | xminTrack[0] =-5; xmaxTrack[0]=5; // | |
264 | fHistoPull[5] = new THnSparseS("#Delta_{pt}/p_{t} (unit)","#Delta_{pt}/p_{t} (unit)", 9, binsTrack,xminTrack, xmaxTrack); | |
265 | // | |
8a92e133 | 266 | |
267 | for (Int_t idedx=0;idedx<4;idedx++){ | |
268 | xminTrack[0] =0.5; xmaxTrack[0]=1.5; // | |
269 | binsTrack[1] =40; | |
270 | xminTrack[1] =10; xmaxTrack[1]=160; | |
271 | ||
272 | fHistodEdxMax[idedx] = new THnSparseS(Form("dEdx_{MaxUp}/dEdx_{MaxDown}_Pad%d",idedx), | |
273 | Form("dEdx_{MaxUp}/dEdx_{MaxDown}_Pad%d",idedx), | |
274 | 9, binsTrack,xminTrack, xmaxTrack); | |
275 | fHistodEdxTot[idedx] = new THnSparseS(Form("dEdx_{TotUp}/dEdx_{TotDown}_Pad%d",idedx), | |
276 | Form("dEdx_{TotUp}/dEdx_{TotDown}_Pad%d",idedx), | |
277 | 9, binsTrack,xminTrack, xmaxTrack); | |
278 | } | |
279 | ||
280 | ||
281 | ||
91fd44c9 | 282 | for (Int_t ivar=0;ivar<6;ivar++){ |
8a92e133 | 283 | for (Int_t ivar2=0;ivar2<9;ivar2++){ |
91fd44c9 | 284 | fHistoDelta[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data()); |
285 | fHistoDelta[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data()); | |
286 | fHistoPull[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data()); | |
287 | fHistoPull[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data()); | |
288 | BinLogX(fHistoDelta[ivar],7); | |
289 | BinLogX(fHistoPull[ivar],7); | |
8a92e133 | 290 | if (ivar<4){ |
291 | fHistodEdxMax[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data()); | |
292 | fHistodEdxMax[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data()); | |
293 | fHistodEdxTot[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data()); | |
294 | fHistodEdxTot[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data()); | |
295 | BinLogX(fHistodEdxMax[ivar],7); | |
296 | BinLogX(fHistodEdxTot[ivar],7); | |
297 | } | |
91fd44c9 | 298 | } |
299 | } | |
300 | } | |
301 | ||
302 | ||
303 | void AliTPCcalibCosmic::Add(const AliTPCcalibCosmic* cosmic){ | |
304 | // | |
305 | // | |
306 | // | |
307 | for (Int_t ivar=0; ivar<6;ivar++){ | |
308 | if (fHistoDelta[ivar] && cosmic->fHistoDelta[ivar]){ | |
309 | fHistoDelta[ivar]->Add(cosmic->fHistoDelta[ivar]); | |
310 | } | |
311 | if (fHistoPull[ivar] && cosmic->fHistoPull[ivar]){ | |
312 | fHistoPull[ivar]->Add(cosmic->fHistoPull[ivar]); | |
313 | } | |
314 | } | |
8a92e133 | 315 | for (Int_t ivar=0; ivar<4;ivar++){ |
316 | if (fHistodEdxMax[ivar] && cosmic->fHistodEdxMax[ivar]){ | |
317 | fHistodEdxMax[ivar]->Add(cosmic->fHistodEdxMax[ivar]); | |
318 | } | |
319 | if (fHistodEdxTot[ivar] && cosmic->fHistodEdxTot[ivar]){ | |
320 | fHistodEdxTot[ivar]->Add(cosmic->fHistodEdxTot[ivar]); | |
321 | } | |
322 | } | |
91fd44c9 | 323 | } |
324 | ||
9b27d39b | 325 | |
326 | ||
327 | ||
f7f33dec | 328 | void AliTPCcalibCosmic::Process(AliESDEvent *event) { |
9b27d39b | 329 | // |
330 | // | |
331 | // | |
f7f33dec | 332 | if (!event) { |
333 | Printf("ERROR: ESD not available"); | |
334 | return; | |
9b27d39b | 335 | } |
76273318 | 336 | AliESDfriend *esdFriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend")); |
337 | if (!esdFriend) { | |
338 | Printf("ERROR: esdFriend not available"); | |
f7f33dec | 339 | return; |
340 | } | |
2acad464 | 341 | |
108953e9 | 342 | |
54b76c13 | 343 | FindPairs(event); // nearly everything takes place in find pairs... |
344 | ||
2acad464 | 345 | if (GetDebugLevel()>20) printf("Hallo world: Im here and processing an event\n"); |
9b27d39b | 346 | Int_t ntracks=event->GetNumberOfTracks(); |
347 | fHistNTracks->Fill(ntracks); | |
9b27d39b | 348 | if (ntracks==0) return; |
cc65e4f5 | 349 | // AliESDcosmic cosmicESD; |
350 | // TTreeSRedirector * cstream = GetDebugStreamer(); | |
351 | // cosmicESD.SetDebugStreamer(cstream); | |
352 | // cosmicESD.ProcessEvent(event); | |
353 | // if (cstream) cosmicESD.DumpToTree(); | |
15e48021 | 354 | |
355 | ||
54b76c13 | 356 | } |
9b27d39b | 357 | |
f7f33dec | 358 | |
3326b323 | 359 | void AliTPCcalibCosmic::FillHistoPerformance(const AliExternalTrackParam *par0, const AliExternalTrackParam *par1, const AliExternalTrackParam *inner0, const AliExternalTrackParam */*inner1*/, AliTPCseed *seed0, AliTPCseed *seed1, const AliExternalTrackParam *param0Combined ){ |
91fd44c9 | 360 | // |
a390f11f | 361 | // par0,par1 - parameter of tracks at DCA 0 |
362 | // inner0,inner1 - parameter of tracks at the TPC entrance | |
363 | // seed0, seed1 - detailed track information | |
364 | // param0Combined - Use combined track parameters for binning | |
91fd44c9 | 365 | // |
8a92e133 | 366 | Int_t kMinCldEdx =20; |
367 | Int_t ncl0 = seed0->GetNumberOfClusters(); | |
368 | Int_t ncl1 = seed1->GetNumberOfClusters(); | |
369 | ||
91fd44c9 | 370 | const Double_t kpullCut = 10; |
371 | Double_t x[9]; | |
372 | Double_t xyz0[3]; | |
373 | Double_t xyz1[3]; | |
374 | par0->GetXYZ(xyz0); | |
375 | par1->GetXYZ(xyz1); | |
376 | Double_t radius0 = TMath::Sqrt(xyz0[0]*xyz0[0]+xyz0[1]*xyz0[1]); | |
377 | Double_t radius1 = TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]); | |
378 | inner0->GetXYZ(xyz0); | |
379 | Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]); | |
380 | // bin parameters | |
381 | x[1] = TMath::Min(ncl0,ncl1); | |
382 | x[2] = (radius0+radius1)*0.5; | |
a390f11f | 383 | x[3] = param0Combined->GetZ(); |
384 | x[4] = inner0->GetSnp(); | |
385 | x[5] = param0Combined->GetTgl(); | |
386 | x[6] = param0Combined->GetSigned1Pt(); | |
387 | x[7] = param0Combined->Pt(); | |
91fd44c9 | 388 | x[8] = alpha; |
389 | // deltas | |
390 | Double_t delta[6]; | |
391 | Double_t sigma[6]; | |
392 | delta[0] = (par0->GetY()+par1->GetY()); | |
393 | delta[1] = (par0->GetZ()-par1->GetZ()); | |
394 | delta[2] = (par0->GetAlpha()-par1->GetAlpha()-TMath::Pi()); | |
395 | delta[3] = (par0->GetTgl()+par1->GetTgl()); | |
396 | delta[4] = (par0->GetParameter()[4]+par1->GetParameter()[4]); | |
397 | delta[5] = (par0->Pt()-par1->Pt())/((par0->Pt()+par1->Pt())*0.5); | |
398 | // | |
399 | sigma[0] = TMath::Sqrt(par0->GetSigmaY2()+par1->GetSigmaY2()); | |
400 | sigma[1] = TMath::Sqrt(par0->GetSigmaZ2()+par1->GetSigmaZ2()); | |
401 | sigma[2] = TMath::Sqrt(par0->GetSigmaSnp2()+par1->GetSigmaSnp2()); | |
402 | sigma[3] = TMath::Sqrt(par0->GetSigmaTgl2()+par1->GetSigmaTgl2()); | |
403 | sigma[4] = TMath::Sqrt(par0->GetSigma1Pt2()+par1->GetSigma1Pt2()); | |
404 | sigma[5] = sigma[4]*((par0->Pt()+par1->Pt())*0.5); | |
405 | // | |
406 | Bool_t isOK = kTRUE; | |
407 | for (Int_t ivar=0;ivar<6;ivar++){ | |
408 | if (sigma[ivar]==0) isOK=kFALSE; | |
409 | x[0]= delta[ivar]/sigma[ivar]; | |
410 | if (TMath::Abs(x[0])>kpullCut) isOK = kFALSE; | |
411 | } | |
412 | // | |
413 | ||
414 | if (isOK) for (Int_t ivar=0;ivar<6;ivar++){ | |
415 | x[0]= delta[ivar]/TMath::Sqrt(2); | |
416 | if (ivar==2 || ivar ==3) x[0]*=1000; | |
417 | fHistoDelta[ivar]->Fill(x); | |
418 | if (sigma[ivar]>0){ | |
419 | x[0]= delta[ivar]/sigma[ivar]; | |
420 | fHistoPull[ivar]->Fill(x); | |
421 | } | |
422 | } | |
8a92e133 | 423 | |
424 | // | |
425 | // Fill dedx performance | |
426 | // | |
427 | for (Int_t ipad=0; ipad<4;ipad++){ | |
428 | // | |
429 | // | |
430 | // | |
431 | Int_t row0=0; | |
432 | Int_t row1=160; | |
433 | if (ipad==0) row1=63; | |
434 | if (ipad==1) {row0=63; row1=63+64;} | |
435 | if (ipad==2) {row0=128;} | |
436 | Int_t nclUp = TMath::Nint(seed0->CookdEdxAnalytical(0.01,0.7,0,row0,row1,2)); | |
437 | Int_t nclDown = TMath::Nint(seed1->CookdEdxAnalytical(0.01,0.7,0,row0,row1,2)); | |
438 | Int_t minCl = TMath::Min(nclUp,nclDown); | |
439 | if (minCl<kMinCldEdx) continue; | |
440 | x[1] = minCl; | |
441 | // | |
442 | Float_t dEdxTotUp = seed0->CookdEdxAnalytical(0.01,0.7,0,row0,row1); | |
443 | Float_t dEdxTotDown = seed1->CookdEdxAnalytical(0.01,0.7,0,row0,row1); | |
444 | Float_t dEdxMaxUp = seed0->CookdEdxAnalytical(0.01,0.7,1,row0,row1); | |
445 | Float_t dEdxMaxDown = seed1->CookdEdxAnalytical(0.01,0.7,1,row0,row1); | |
446 | // | |
447 | if (dEdxTotDown<=0) continue; | |
448 | if (dEdxMaxDown<=0) continue; | |
449 | x[0]=dEdxTotUp/dEdxTotDown; | |
450 | fHistodEdxTot[ipad]->Fill(x); | |
451 | x[0]=dEdxMaxUp/dEdxMaxDown; | |
452 | fHistodEdxMax[ipad]->Fill(x); | |
453 | } | |
454 | ||
455 | ||
91fd44c9 | 456 | |
457 | } | |
458 | ||
3326b323 | 459 | void AliTPCcalibCosmic::MaterialBudgetDump(AliExternalTrackParam *const par0, AliExternalTrackParam *const par1, const AliExternalTrackParam *inner0, const AliExternalTrackParam *inner1, AliTPCseed *const seed0, AliTPCseed *const seed1, AliExternalTrackParam *const param0Combined, AliExternalTrackParam *const param1Combined){ |
a390f11f | 460 | // |
461 | // matrial budget AOD dump | |
462 | // | |
463 | // par0,par1 - parameter of tracks at DCA 0 | |
464 | // inner0,inner1 - parameter of tracks at the TPC entrance | |
465 | // seed0, seed1 - detailed track information | |
466 | // param0Combined - Use combined track parameters for binning | |
467 | // param1Combined - | |
468 | Double_t p0In = inner0->GetP(); | |
469 | Double_t p1In = inner1->GetP(); | |
470 | Double_t p0V = par0->GetP(); | |
471 | Double_t p1V = par1->GetP(); | |
472 | // | |
473 | Double_t pt0In = inner0->Pt(); | |
474 | Double_t pt1In = inner1->Pt(); | |
475 | Double_t pt0V = par0->Pt(); | |
476 | Double_t pt1V = par1->Pt(); | |
477 | Int_t ncl0 = seed0->GetNumberOfClusters(); | |
478 | Int_t ncl1 = seed1->GetNumberOfClusters(); | |
479 | Int_t nclmin=TMath::Min(ncl0,ncl1); | |
480 | Double_t sign = (param0Combined->GetSigned1Pt()>0) ? 1:-1.; | |
481 | // | |
482 | TTreeSRedirector * pcstream = GetDebugStreamer(); | |
483 | if (pcstream){ | |
484 | (*pcstream)<<"material"<< | |
485 | "run="<<fRun<< // run number | |
486 | "event="<<fEvent<< // event number | |
487 | "time="<<fTime<< // time stamp of event | |
488 | "trigger="<<fTrigger<< // trigger | |
489 | "triggerClass="<<&fTriggerClass<< // trigger | |
490 | "mag="<<fMagF<< // magnetic field | |
491 | "sign="<<sign<< // sign of the track | |
492 | // | |
493 | "ncl0="<<ncl0<< | |
494 | "ncl1="<<ncl1<< | |
495 | "nclmin="<<nclmin<< | |
496 | // | |
497 | "p0In="<<p0In<< | |
498 | "p1In="<<p1In<< | |
499 | "p0V="<<p0V<< | |
500 | "p1V="<<p1V<< | |
501 | "pt0In="<<pt0In<< | |
502 | "pt1In="<<pt1In<< | |
503 | "pt0V="<<pt0V<< | |
504 | "pt1V="<<pt1V<< | |
505 | "p0.="<<par0<< | |
506 | "p1.="<<par1<< | |
507 | "up0.="<<param0Combined<< | |
508 | "up1.="<<param1Combined<< | |
509 | "\n"; | |
510 | } | |
511 | ||
512 | } | |
f7f33dec | 513 | |
54b76c13 | 514 | void AliTPCcalibCosmic::Analyze() { |
515 | ||
516 | fMIPvalue = CalculateMIPvalue(fDeDxMIP); | |
517 | ||
518 | return; | |
519 | ||
520 | } | |
521 | ||
f7f33dec | 522 | |
f7f33dec | 523 | |
9b27d39b | 524 | void AliTPCcalibCosmic::FindPairs(AliESDEvent *event) { |
525 | // | |
526 | // Find cosmic pairs | |
04087794 | 527 | // |
528 | // Track0 is choosen in upper TPC part | |
529 | // Track1 is choosen in lower TPC part | |
9b27d39b | 530 | // |
2acad464 | 531 | if (GetDebugLevel()>20) printf("Hallo world: Im here\n"); |
76273318 | 532 | AliESDfriend *esdFriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend")); |
9b27d39b | 533 | Int_t ntracks=event->GetNumberOfTracks(); |
9b27d39b | 534 | TObjArray tpcSeeds(ntracks); |
535 | if (ntracks==0) return; | |
536 | Double_t vtxx[3]={0,0,0}; | |
537 | Double_t svtxx[3]={0.000001,0.000001,100.}; | |
538 | AliESDVertex vtx(vtxx,svtxx); | |
539 | // | |
540 | //track loop | |
541 | // | |
54b76c13 | 542 | for (Int_t i=0;i<ntracks;++i) { |
543 | AliESDtrack *track = event->GetTrack(i); | |
544 | fClusters->Fill(track->GetTPCNcls()); | |
545 | ||
546 | const AliExternalTrackParam * trackIn = track->GetInnerParam(); | |
547 | const AliExternalTrackParam * trackOut = track->GetOuterParam(); | |
548 | if (!trackIn) continue; | |
549 | if (!trackOut) continue; | |
860b3d93 | 550 | if (ntracks>4 && TMath::Abs(trackIn->GetTgl())<0.0015) continue; // filter laser |
551 | ||
552 | ||
76273318 | 553 | AliESDfriendTrack *friendTrack = esdFriend->GetTrack(i); |
9b27d39b | 554 | TObject *calibObject; |
555 | AliTPCseed *seed = 0; | |
556 | for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) { | |
557 | if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break; | |
558 | } | |
559 | if (seed) tpcSeeds.AddAt(seed,i); | |
54b76c13 | 560 | |
561 | Double_t meanP = 0.5*(trackIn->GetP() + trackOut->GetP()); | |
562 | if (seed && track->GetTPCNcls() > 80 + 60/(1+TMath::Exp(-meanP+5))) { | |
15e48021 | 563 | fDeDx->Fill(meanP, seed->CookdEdxNorm(0.0,0.45,0,0,159)); |
54b76c13 | 564 | // |
15e48021 | 565 | if (meanP > 0.4 && meanP < 0.45) fDeDxMIP->Fill(seed->CookdEdxNorm(0.0,0.45,0,0,159)); |
54b76c13 | 566 | // |
82628455 | 567 | // if (GetDebugLevel()>0&&meanP>0.2&&seed->CookdEdxNorm(0.0,0.45,0,0,159)>300) { |
568 | // //TFile *curfile = AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile(); | |
569 | // //if (curfile) printf(">>> p+ in file: %s \t event: %i \t Number of ESD tracks: %i \n", curfile->GetName(), (int)event->GetEventNumberInFile(), (int)ntracks); | |
570 | // // if (track->GetOuterParam()->GetAlpha()<0) cout << " Polartiy: " << track->GetSign() << endl; | |
571 | // } | |
54b76c13 | 572 | |
573 | } | |
574 | ||
9b27d39b | 575 | } |
54b76c13 | 576 | |
9b27d39b | 577 | if (ntracks<2) return; |
578 | // | |
579 | // Find pairs | |
580 | // | |
581 | for (Int_t i=0;i<ntracks;++i) { | |
582 | AliESDtrack *track0 = event->GetTrack(i); | |
04087794 | 583 | // track0 - choosen upper part |
584 | if (!track0) continue; | |
585 | if (!track0->GetOuterParam()) continue; | |
586 | if (track0->GetOuterParam()->GetAlpha()<0) continue; | |
e9362f9d | 587 | Double_t dir0[3]; |
588 | track0->GetDirection(dir0); | |
04087794 | 589 | for (Int_t j=0;j<ntracks;++j) { |
590 | if (i==j) continue; | |
591 | AliESDtrack *track1 = event->GetTrack(j); | |
592 | //track 1 lower part | |
593 | if (!track1) continue; | |
594 | if (!track1->GetOuterParam()) continue; | |
595 | if (track1->GetOuterParam()->GetAlpha()>0) continue; | |
596 | // | |
e9362f9d | 597 | Double_t dir1[3]; |
598 | track1->GetDirection(dir1); | |
54b76c13 | 599 | |
9b27d39b | 600 | AliTPCseed * seed0 = (AliTPCseed*) tpcSeeds.At(i); |
601 | AliTPCseed * seed1 = (AliTPCseed*) tpcSeeds.At(j); | |
602 | if (! seed0) continue; | |
603 | if (! seed1) continue; | |
15e48021 | 604 | Float_t dedx0 = seed0->CookdEdxNorm(0.05,0.55,0,0,159); |
605 | Float_t dedx1 = seed1->CookdEdxNorm(0.05,0.55,0,0,159); | |
bca20570 | 606 | // |
15e48021 | 607 | Float_t dedx0I = seed0->CookdEdxNorm(0.05,0.55,0,0,63); |
608 | Float_t dedx1I = seed1->CookdEdxNorm(0.05,0.55,0,0,63); | |
bca20570 | 609 | // |
15e48021 | 610 | Float_t dedx0O = seed0->CookdEdxNorm(0.05,0.55,0,64,159); |
611 | Float_t dedx1O = seed1->CookdEdxNorm(0.05,0.55,0,64,159); | |
bca20570 | 612 | // |
e9362f9d | 613 | Float_t dir = (dir0[0]*dir1[0] + dir0[1]*dir1[1] + dir0[2]*dir1[2]); |
9b27d39b | 614 | Float_t d0 = track0->GetLinearD(0,0); |
615 | Float_t d1 = track1->GetLinearD(0,0); | |
616 | // | |
617 | // conservative cuts - convergence to be guarantied | |
618 | // applying before track propagation | |
619 | if (TMath::Abs(d0+d1)>fCutMaxD) continue; // distance to the 0,0 | |
620 | if (dir>fCutMinDir) continue; // direction vector product | |
621 | Float_t bz = AliTracker::GetBz(); | |
622 | Float_t dvertex0[2]; //distance to 0,0 | |
623 | Float_t dvertex1[2]; //distance to 0,0 | |
624 | track0->GetDZ(0,0,0,bz,dvertex0); | |
625 | track1->GetDZ(0,0,0,bz,dvertex1); | |
626 | if (TMath::Abs(dvertex0[1])>250) continue; | |
627 | if (TMath::Abs(dvertex1[1])>250) continue; | |
628 | // | |
629 | // | |
630 | // | |
631 | Float_t dmax = TMath::Max(TMath::Abs(d0),TMath::Abs(d1)); | |
632 | AliExternalTrackParam param0(*track0); | |
633 | AliExternalTrackParam param1(*track1); | |
634 | // | |
635 | // Propagate using Magnetic field and correct fo material budget | |
636 | // | |
a390f11f | 637 | Double_t sign0=-1; |
638 | Double_t sign1=1; | |
639 | Double_t maxsnp=0.90; | |
9963b5e2 | 640 | AliTracker::PropagateTrackToBxByBz(¶m0,dmax+1,TDatabasePDG::Instance()->GetParticle("e-")->Mass(),3,kTRUE,maxsnp,sign0); |
641 | AliTracker::PropagateTrackToBxByBz(¶m1,dmax+1,TDatabasePDG::Instance()->GetParticle("e-")->Mass(),3,kTRUE,maxsnp,sign1); | |
9b27d39b | 642 | // |
643 | // Propagate rest to the 0,0 DCA - z should be ignored | |
644 | // | |
645 | Bool_t b0 = param0.PropagateToDCA(&vtx,bz,1000); | |
646 | Bool_t b1 = param1.PropagateToDCA(&vtx,bz,1000); | |
647 | // | |
648 | param0.GetDZ(0,0,0,bz,dvertex0); | |
649 | param1.GetDZ(0,0,0,bz,dvertex1); | |
a6dc0cf6 | 650 | if (TMath::Abs(param0.GetZ()-param1.GetZ())>fCutMaxDz) continue; |
9b27d39b | 651 | // |
652 | Double_t xyz0[3];//,pxyz0[3]; | |
653 | Double_t xyz1[3];//,pxyz1[3]; | |
654 | param0.GetXYZ(xyz0); | |
655 | param1.GetXYZ(xyz1); | |
656 | Bool_t isPair = IsPair(¶m0,¶m1); | |
657 | // | |
54b76c13 | 658 | if (isPair) FillAcordeHist(track0); |
15e48021 | 659 | // |
660 | // combined track params | |
661 | // | |
662 | AliExternalTrackParam *par0U=MakeCombinedTrack(¶m0,¶m1); | |
663 | AliExternalTrackParam *par1U=MakeCombinedTrack(¶m1,¶m0); | |
664 | ||
665 | ||
54b76c13 | 666 | // |
9b27d39b | 667 | if (fStreamLevel>0){ |
668 | TTreeSRedirector * cstream = GetDebugStreamer(); | |
108953e9 | 669 | //printf("My stream=%p\n",(void*)cstream); |
e9362f9d | 670 | AliExternalTrackParam *ip0 = (AliExternalTrackParam *)track0->GetInnerParam(); |
671 | AliExternalTrackParam *ip1 = (AliExternalTrackParam *)track1->GetInnerParam(); | |
672 | AliExternalTrackParam *op0 = (AliExternalTrackParam *)track0->GetOuterParam(); | |
673 | AliExternalTrackParam *op1 = (AliExternalTrackParam *)track1->GetOuterParam(); | |
674 | Bool_t isCrossI = ip0->GetZ()*ip1->GetZ()<0; | |
675 | Bool_t isCrossO = op0->GetZ()*op1->GetZ()<0; | |
676 | Double_t alpha0 = TMath::ATan2(dir0[1],dir0[0]); | |
677 | Double_t alpha1 = TMath::ATan2(dir1[1],dir1[0]); | |
91fd44c9 | 678 | // |
679 | // | |
680 | // | |
a390f11f | 681 | FillHistoPerformance(¶m0, ¶m1, ip0, ip1, seed0, seed1,par0U); |
682 | MaterialBudgetDump(¶m0, ¶m1, ip0, ip1, seed0, seed1,par0U,par1U); | |
9b27d39b | 683 | if (cstream) { |
684 | (*cstream) << "Track0" << | |
108953e9 | 685 | "run="<<fRun<< // run number |
686 | "event="<<fEvent<< // event number | |
687 | "time="<<fTime<< // time stamp of event | |
688 | "trigger="<<fTrigger<< // trigger | |
15e48021 | 689 | "triggerClass="<<&fTriggerClass<< // trigger |
108953e9 | 690 | "mag="<<fMagF<< // magnetic field |
9b27d39b | 691 | "dir="<<dir<< // direction |
54b76c13 | 692 | "OK="<<isPair<< // will be accepted |
9b27d39b | 693 | "b0="<<b0<< // propagate status |
694 | "b1="<<b1<< // propagate status | |
e9362f9d | 695 | "crossI="<<isCrossI<< // cross inner |
696 | "crossO="<<isCrossO<< // cross outer | |
697 | // | |
9b27d39b | 698 | "Orig0.=" << track0 << // original track 0 |
699 | "Orig1.=" << track1 << // original track 1 | |
700 | "Tr0.="<<¶m0<< // track propagated to the DCA 0,0 | |
701 | "Tr1.="<<¶m1<< // track propagated to the DCA 0,0 | |
e9362f9d | 702 | "Ip0.="<<ip0<< // inner param - upper |
703 | "Ip1.="<<ip1<< // inner param - lower | |
704 | "Op0.="<<op0<< // outer param - upper | |
705 | "Op1.="<<op1<< // outer param - lower | |
15e48021 | 706 | "Up0.="<<par0U<< // combined track 0 |
707 | "Up1.="<<par1U<< // combined track 1 | |
e9362f9d | 708 | // |
9b27d39b | 709 | "v00="<<dvertex0[0]<< // distance using kalman |
710 | "v01="<<dvertex0[1]<< // | |
711 | "v10="<<dvertex1[0]<< // | |
712 | "v11="<<dvertex1[1]<< // | |
713 | "d0="<<d0<< // linear distance to 0,0 | |
714 | "d1="<<d1<< // linear distance to 0,0 | |
715 | // | |
e9362f9d | 716 | // |
717 | // | |
718 | "x00="<<xyz0[0]<< // global position close to vertex | |
9b27d39b | 719 | "x01="<<xyz0[1]<< |
720 | "x02="<<xyz0[2]<< | |
721 | // | |
e9362f9d | 722 | "x10="<<xyz1[0]<< // global position close to vertex |
9b27d39b | 723 | "x11="<<xyz1[1]<< |
724 | "x12="<<xyz1[2]<< | |
725 | // | |
e9362f9d | 726 | "alpha0="<<alpha0<< |
727 | "alpha1="<<alpha1<< | |
728 | "dir00="<<dir0[0]<< // direction upper | |
729 | "dir01="<<dir0[1]<< | |
730 | "dir02="<<dir0[2]<< | |
731 | // | |
732 | "dir10="<<dir1[0]<< // direction lower | |
733 | "dir11="<<dir1[1]<< | |
734 | "dir12="<<dir1[2]<< | |
735 | // | |
736 | // | |
54b76c13 | 737 | "Seed0.=" << seed0 << // original seed 0 |
738 | "Seed1.=" << seed1 << // original seed 1 | |
bca20570 | 739 | // |
740 | "dedx0="<<dedx0<< // dedx0 - all | |
741 | "dedx1="<<dedx1<< // dedx1 - all | |
742 | // | |
54b76c13 | 743 | "dedx0I="<<dedx0I<< // dedx0 - inner ROC |
744 | "dedx1I="<<dedx1I<< // dedx1 - inner ROC | |
bca20570 | 745 | // |
54b76c13 | 746 | "dedx0O="<<dedx0O<< // dedx0 - outer ROC |
747 | "dedx1O="<<dedx1O<< // dedx1 - outer ROC | |
9b27d39b | 748 | "\n"; |
749 | } | |
15e48021 | 750 | } |
751 | delete par0U; | |
752 | delete par1U; | |
9b27d39b | 753 | } |
754 | } | |
755 | } | |
756 | ||
9b27d39b | 757 | |
bca20570 | 758 | |
759 | ||
54b76c13 | 760 | void AliTPCcalibCosmic::FillAcordeHist(AliESDtrack *upperTrack) { |
761 | ||
762 | // Pt cut to select straight tracks which can be easily propagated to ACORDE which is outside the magnetic field | |
763 | if (upperTrack->Pt() < 10 || upperTrack->GetTPCNcls() < 80) return; | |
764 | ||
76273318 | 765 | const Double_t acordePlane = 850.; // distance of the central Acorde detectors to the beam line at y =0 |
54b76c13 | 766 | const Double_t roof = 210.5; // distance from x =0 to end of magnet roof |
767 | ||
768 | Double_t r[3]; | |
769 | upperTrack->GetXYZ(r); | |
770 | Double_t d[3]; | |
771 | upperTrack->GetDirection(d); | |
772 | Double_t x,z; | |
76273318 | 773 | z = r[2] + (d[2]/d[1])*(acordePlane - r[1]); |
774 | x = r[0] + (d[0]/d[1])*(acordePlane - r[1]); | |
54b76c13 | 775 | |
776 | if (x > roof) { | |
76273318 | 777 | x = r[0] + (d[0]/(d[0]+d[1]))*(acordePlane+roof-r[0]-r[1]); |
778 | z = r[2] + (d[2]/(d[0]+d[1]))*(acordePlane+roof-r[0]-r[1]); | |
54b76c13 | 779 | } |
780 | if (x < -roof) { | |
76273318 | 781 | x = r[0] + (d[0]/(d[1]-d[0]))*(acordePlane+roof+r[0]-r[1]); |
782 | z = r[2] + (d[2]/(d[1]-d[0]))*(acordePlane+roof+r[0]-r[1]); | |
54b76c13 | 783 | } |
784 | ||
785 | fModules->Fill(z, x); | |
786 | ||
787 | } | |
788 | ||
789 | ||
790 | ||
3326b323 | 791 | Long64_t AliTPCcalibCosmic::Merge(TCollection *const li) { |
bca20570 | 792 | |
793 | TIterator* iter = li->MakeIterator(); | |
794 | AliTPCcalibCosmic* cal = 0; | |
795 | ||
796 | while ((cal = (AliTPCcalibCosmic*)iter->Next())) { | |
797 | if (!cal->InheritsFrom(AliTPCcalibCosmic::Class())) { | |
860b3d93 | 798 | //Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName()); |
bca20570 | 799 | return -1; |
800 | } | |
801 | ||
54b76c13 | 802 | fHistNTracks->Add(cal->GetHistNTracks()); |
803 | fClusters->Add(cal-> GetHistClusters()); | |
804 | fModules->Add(cal->GetHistAcorde()); | |
805 | fHistPt->Add(cal->GetHistPt()); | |
806 | fDeDx->Add(cal->GetHistDeDx()); | |
807 | fDeDxMIP->Add(cal->GetHistMIP()); | |
91fd44c9 | 808 | Add(cal); |
bca20570 | 809 | } |
bca20570 | 810 | return 0; |
9b27d39b | 811 | |
f7f33dec | 812 | } |
813 | ||
54b76c13 | 814 | |
9b27d39b | 815 | Bool_t AliTPCcalibCosmic::IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1){ |
816 | // | |
817 | // | |
818 | /* | |
819 | // 0. Same direction - OPOSITE - cutDir +cutT | |
820 | TCut cutDir("cutDir","dir<-0.99") | |
821 | // 1. | |
822 | TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03") | |
823 | // | |
824 | // 2. The same rphi | |
825 | TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5") | |
826 | // | |
827 | // | |
828 | // | |
829 | TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10"); | |
830 | // 1/Pt diff cut | |
831 | */ | |
832 | const Double_t *p0 = tr0->GetParameter(); | |
833 | const Double_t *p1 = tr1->GetParameter(); | |
834 | if (TMath::Abs(p0[3]+p1[3])>fCutTheta) return kFALSE; | |
a6dc0cf6 | 835 | if (TMath::Abs(p0[1]-p1[1])>fCutMaxDz) return kFALSE; |
9b27d39b | 836 | if (TMath::Abs(p0[0]+p1[0])>fCutMaxD) return kFALSE; |
a6dc0cf6 | 837 | |
9b27d39b | 838 | Double_t d0[3], d1[3]; |
839 | tr0->GetDirection(d0); | |
840 | tr1->GetDirection(d1); | |
841 | if (d0[0]*d1[0] + d0[1]*d1[1] + d0[2]*d1[2] >fCutMinDir) return kFALSE; | |
842 | // | |
843 | return kTRUE; | |
844 | } | |
54b76c13 | 845 | |
846 | ||
847 | ||
848 | Double_t AliTPCcalibCosmic::CalculateMIPvalue(TH1F * hist) { | |
849 | ||
850 | TF1 * funcDoubleGaus = new TF1("funcDoubleGaus", "gaus(0)+gaus(3)",0,1000); | |
851 | funcDoubleGaus->SetParameters(hist->GetEntries()*0.75,hist->GetMean()/1.3,hist->GetMean()*0.10, | |
852 | hist->GetEntries()*0.25,hist->GetMean()*1.3,hist->GetMean()*0.10); | |
853 | hist->Fit(funcDoubleGaus); | |
3326b323 | 854 | Double_t aMIPvalue = TMath::Min(funcDoubleGaus->GetParameter(1),funcDoubleGaus->GetParameter(4)); |
54b76c13 | 855 | |
856 | delete funcDoubleGaus; | |
857 | ||
3326b323 | 858 | return aMIPvalue; |
54b76c13 | 859 | |
860 | } | |
9b27d39b | 861 | |
862 | ||
f7f33dec | 863 | |
54b76c13 | 864 | |
57dc06f2 | 865 | void AliTPCcalibCosmic::CalculateBetheParams(TH2F */*hist*/, Double_t * /*initialParam*/) { |
866 | // | |
867 | // Not implemented yet | |
868 | // | |
54b76c13 | 869 | return; |
870 | ||
871 | } | |
872 | ||
873 | ||
3326b323 | 874 | void AliTPCcalibCosmic::BinLogX(THnSparse *const h, Int_t axisDim) { |
91fd44c9 | 875 | |
876 | // Method for the correct logarithmic binning of histograms | |
877 | ||
878 | TAxis *axis = h->GetAxis(axisDim); | |
879 | int bins = axis->GetNbins(); | |
880 | ||
881 | Double_t from = axis->GetXmin(); | |
882 | Double_t to = axis->GetXmax(); | |
3326b323 | 883 | Double_t *newBins = new Double_t[bins + 1]; |
91fd44c9 | 884 | |
3326b323 | 885 | newBins[0] = from; |
91fd44c9 | 886 | Double_t factor = pow(to/from, 1./bins); |
887 | ||
888 | for (int i = 1; i <= bins; i++) { | |
3326b323 | 889 | newBins[i] = factor * newBins[i-1]; |
91fd44c9 | 890 | } |
3326b323 | 891 | axis->Set(bins, newBins); |
892 | delete newBins; | |
91fd44c9 | 893 | |
894 | } | |
895 | ||
896 | ||
3326b323 | 897 | void AliTPCcalibCosmic::BinLogX(TH1 *const h) { |
f7f33dec | 898 | |
899 | // Method for the correct logarithmic binning of histograms | |
900 | ||
901 | TAxis *axis = h->GetXaxis(); | |
902 | int bins = axis->GetNbins(); | |
903 | ||
904 | Double_t from = axis->GetXmin(); | |
905 | Double_t to = axis->GetXmax(); | |
3326b323 | 906 | Double_t *newBins = new Double_t[bins + 1]; |
f7f33dec | 907 | |
3326b323 | 908 | newBins[0] = from; |
f7f33dec | 909 | Double_t factor = pow(to/from, 1./bins); |
910 | ||
911 | for (int i = 1; i <= bins; i++) { | |
3326b323 | 912 | newBins[i] = factor * newBins[i-1]; |
f7f33dec | 913 | } |
3326b323 | 914 | axis->Set(bins, newBins); |
915 | delete newBins; | |
f7f33dec | 916 | |
917 | } | |
918 | ||
7b18d067 | 919 | |
860b3d93 | 920 | AliExternalTrackParam *AliTPCcalibCosmic::MakeTrack(const AliExternalTrackParam *track0, const AliExternalTrackParam *track1){ |
921 | // | |
922 | // | |
923 | // | |
924 | AliExternalTrackParam *par1R= new AliExternalTrackParam(*track1); | |
925 | par1R->Rotate(track0->GetAlpha()); | |
15e48021 | 926 | par1R->PropagateTo(track0->GetX(),AliTracker::GetBz()); |
860b3d93 | 927 | // |
928 | // | |
929 | Double_t * param = (Double_t*)par1R->GetParameter(); | |
930 | Double_t * covar = (Double_t*)par1R->GetCovariance(); | |
15e48021 | 931 | |
860b3d93 | 932 | param[0]*=1; //OK |
933 | param[1]*=1; //OK | |
934 | param[2]*=1; //? | |
935 | param[3]*=-1; //OK | |
936 | param[4]*=-1; //OK | |
937 | // | |
938 | covar[6] *=-1.; covar[7] *=-1.; covar[8] *=-1.; | |
939 | //covar[10]*=-1.; covar[11]*=-1.; covar[12]*=-1.; | |
940 | covar[13]*=-1.; | |
860b3d93 | 941 | return par1R; |
942 | } | |
943 | ||
15e48021 | 944 | AliExternalTrackParam *AliTPCcalibCosmic::MakeCombinedTrack(const AliExternalTrackParam *track0, const AliExternalTrackParam *track1){ |
945 | // | |
946 | // Make combined track | |
947 | // | |
948 | // | |
949 | AliExternalTrackParam * par1T = MakeTrack(track0,track1); | |
950 | AliExternalTrackParam * par0U = new AliExternalTrackParam(*track0); | |
951 | // | |
952 | UpdateTrack(*par0U,*par1T); | |
953 | delete par1T; | |
954 | return par0U; | |
955 | } | |
956 | ||
957 | ||
860b3d93 | 958 | void AliTPCcalibCosmic::UpdateTrack(AliExternalTrackParam &track1, const AliExternalTrackParam &track2){ |
959 | // | |
960 | // Update track 1 with track 2 | |
961 | // | |
962 | // | |
963 | // | |
964 | TMatrixD vecXk(5,1); // X vector | |
965 | TMatrixD covXk(5,5); // X covariance | |
966 | TMatrixD matHk(5,5); // vector to mesurement | |
967 | TMatrixD measR(5,5); // measurement error | |
968 | TMatrixD vecZk(5,1); // measurement | |
969 | // | |
970 | TMatrixD vecYk(5,1); // Innovation or measurement residual | |
971 | TMatrixD matHkT(5,5); | |
972 | TMatrixD matSk(5,5); // Innovation (or residual) covariance | |
973 | TMatrixD matKk(5,5); // Optimal Kalman gain | |
974 | TMatrixD mat1(5,5); // update covariance matrix | |
975 | TMatrixD covXk2(5,5); // | |
976 | TMatrixD covOut(5,5); | |
977 | // | |
978 | Double_t *param1=(Double_t*) track1.GetParameter(); | |
979 | Double_t *covar1=(Double_t*) track1.GetCovariance(); | |
980 | Double_t *param2=(Double_t*) track2.GetParameter(); | |
981 | Double_t *covar2=(Double_t*) track2.GetCovariance(); | |
982 | // | |
983 | // copy data to the matrix | |
984 | for (Int_t ipar=0; ipar<5; ipar++){ | |
860b3d93 | 985 | for (Int_t jpar=0; jpar<5; jpar++){ |
986 | covXk(ipar,jpar) = covar1[track1.GetIndex(ipar, jpar)]; | |
987 | measR(ipar,jpar) = covar2[track2.GetIndex(ipar, jpar)]; | |
15e48021 | 988 | matHk(ipar,jpar)=0; |
989 | mat1(ipar,jpar)=0; | |
860b3d93 | 990 | } |
15e48021 | 991 | vecXk(ipar,0) = param1[ipar]; |
992 | vecZk(ipar,0) = param2[ipar]; | |
993 | matHk(ipar,ipar)=1; | |
994 | mat1(ipar,ipar)=0; | |
860b3d93 | 995 | } |
996 | // | |
997 | // | |
998 | // | |
999 | // | |
860b3d93 | 1000 | // |
1001 | vecYk = vecZk-matHk*vecXk; // Innovation or measurement residual | |
1002 | matHkT=matHk.T(); matHk.T(); | |
1003 | matSk = (matHk*(covXk*matHkT))+measR; // Innovation (or residual) covariance | |
1004 | matSk.Invert(); | |
1005 | matKk = (covXk*matHkT)*matSk; // Optimal Kalman gain | |
1006 | vecXk += matKk*vecYk; // updated vector | |
860b3d93 | 1007 | covXk2 = (mat1-(matKk*matHk)); |
1008 | covOut = covXk2*covXk; | |
1009 | // | |
1010 | // | |
1011 | // | |
1012 | // copy from matrix to parameters | |
1013 | if (0) { | |
1014 | vecXk.Print(); | |
1015 | vecZk.Print(); | |
1016 | // | |
1017 | measR.Print(); | |
1018 | covXk.Print(); | |
1019 | covOut.Print(); | |
1020 | // | |
1021 | track1.Print(); | |
1022 | track2.Print(); | |
1023 | } | |
1024 | ||
1025 | for (Int_t ipar=0; ipar<5; ipar++){ | |
1026 | param1[ipar]= vecXk(ipar,0) ; | |
1027 | for (Int_t jpar=0; jpar<5; jpar++){ | |
1028 | covar1[track1.GetIndex(ipar, jpar)]=covOut(ipar,jpar); | |
1029 | } | |
1030 | } | |
1031 | } | |
1032 | ||
860b3d93 | 1033 | |
1034 | ||
1035 |