]>
Commit | Line | Data |
---|---|---|
10757ee9 | 1 | |
2 | #include <TROOT.h> | |
3 | #include <TChain.h> | |
4 | #include <TFile.h> | |
5 | #include <TH3F.h> | |
5b00528f | 6 | #include <TH2F.h> |
10757ee9 | 7 | // |
8 | #include "TParticle.h" | |
9 | #include "TDatabasePDG.h" | |
10 | #include "AliRunLoader.h" | |
11 | #include "AliStack.h" | |
12 | ||
13 | ||
14 | ||
15 | #include <TPDGCode.h> | |
16 | #include <TStyle.h> | |
17 | #include "TLinearFitter.h" | |
18 | #include "TMatrixD.h" | |
19 | #include "TTreeStream.h" | |
20 | #include "TF1.h" | |
21 | ||
22 | ||
23 | ||
24 | #include "AliMagF.h" | |
25 | #include "AliTracker.h" | |
57e4988a | 26 | #include "AliESDEvent.h" |
10757ee9 | 27 | #include "AliESDtrack.h" |
28 | #include "AliESDfriend.h" | |
29 | #include "AliESDfriendTrack.h" | |
b40afa5e | 30 | #include "AliMathBase.h" |
10757ee9 | 31 | #include "AliTPCseed.h" |
32 | #include "AliTPCclusterMI.h" | |
33 | ||
34 | #include "AliKFParticle.h" | |
35 | #include "AliKFVertex.h" | |
36 | ||
37 | #include "AliTrackPointArray.h" | |
38 | #include "TCint.h" | |
39 | #include "AliTPCcalibV0.h" | |
40 | #include "AliV0.h" | |
46e89793 | 41 | #include "TRandom.h" |
42 | #include "TTreeStream.h" | |
43 | #include "AliTPCcalibDB.h" | |
44 | #include "AliTPCCorrection.h" | |
45 | #include "AliGRPObject.h" | |
46 | #include "AliTPCTransform.h" | |
10757ee9 | 47 | |
48 | ||
49 | ||
50 | ||
51 | ||
52 | ClassImp(AliTPCcalibV0) | |
53 | ||
54 | ||
55 | AliTPCcalibV0::AliTPCcalibV0() : | |
57e4988a | 56 | AliTPCcalibBase(), |
46e89793 | 57 | fV0Tree(0), |
58 | fHPTTree(0), | |
10757ee9 | 59 | fStack(0), |
10757ee9 | 60 | fESD(0), |
61 | fPdg(0), | |
62 | fParticles(0), | |
63 | fV0s(0), | |
46e89793 | 64 | fGammas(0) |
65 | { | |
66 | ||
67 | } | |
68 | AliTPCcalibV0::AliTPCcalibV0(const Text_t *name, const Text_t *title): | |
69 | AliTPCcalibBase(), | |
70 | fV0Tree(0), | |
71 | fHPTTree(0), | |
72 | fStack(0), | |
73 | fESD(0), | |
74 | fPdg(0), | |
75 | fParticles(0), | |
76 | fV0s(0), | |
77 | fGammas(0) | |
10757ee9 | 78 | { |
57e4988a | 79 | fPdg = new TDatabasePDG; |
10757ee9 | 80 | // create output histograms |
46e89793 | 81 | SetName(name); |
82 | SetTitle(title); | |
10757ee9 | 83 | } |
84 | ||
85 | AliTPCcalibV0::~AliTPCcalibV0(){ | |
86 | // | |
87 | // | |
88 | // | |
46e89793 | 89 | delete fV0Tree; |
90 | delete fHPTTree; | |
10757ee9 | 91 | } |
92 | ||
93 | ||
94 | ||
10757ee9 | 95 | |
96 | ||
57e4988a | 97 | void AliTPCcalibV0::ProcessESD(AliESDEvent *esd, AliStack *stack){ |
10757ee9 | 98 | // |
99 | // | |
100 | // | |
101 | fESD = esd; | |
102 | AliKFParticle::SetField(esd->GetMagneticField()); | |
46e89793 | 103 | if (TMath::Abs(AliTracker::GetBz())<1) return; |
104 | DumpToTree(esd); | |
105 | DumpToTreeHPT(esd); | |
106 | return; | |
107 | // | |
10757ee9 | 108 | MakeV0s(); |
109 | if (stack) { | |
110 | fStack = stack; | |
111 | MakeMC(); | |
112 | }else{ | |
113 | fStack =0; | |
114 | } | |
115 | } | |
116 | ||
46e89793 | 117 | void AliTPCcalibV0::DumpToTreeHPT(AliESDEvent *esd){ |
118 | // | |
119 | // Dump V0s fith full firend information to the | |
120 | // | |
121 | if (TMath::Abs(AliTracker::GetBz())<1) return; | |
122 | const Int_t kMinCluster=110; | |
123 | const Float_t kMinPt =3.; | |
124 | AliESDfriend *esdFriend=static_cast<AliESDfriend*>(esd->FindListObject("AliESDfriend")); | |
125 | if (!esdFriend) { | |
126 | Printf("ERROR: esdFriend not available"); | |
127 | return; | |
128 | } | |
129 | // | |
130 | Int_t ntracks=esd->GetNumberOfTracks(); | |
131 | for (Int_t i=0;i<ntracks;++i) { | |
132 | Bool_t isOK=kFALSE; | |
133 | AliESDtrack *track = esd->GetTrack(i); | |
134 | if (track->GetTPCncls()<kMinCluster) continue; | |
135 | AliESDfriendTrack *friendTrack = esdFriend->GetTrack(i); | |
136 | if (!friendTrack) continue; | |
137 | if (TMath::Abs(AliTracker::GetBz())>1){ // cut on momenta if measured | |
138 | if (track->Pt()>kMinPt) isOK=kTRUE; | |
139 | } | |
140 | if (TMath::Abs(AliTracker::GetBz())<1){ // require primary track for the B field OFF data | |
141 | Bool_t isAccepted=kTRUE; | |
142 | if (!track->IsOn(AliESDtrack::kITSrefit)) isAccepted=kFALSE; | |
143 | if (!track->IsOn(AliESDtrack::kTPCrefit)) isAccepted=kFALSE; | |
144 | if (!track->IsOn(AliESDtrack::kTOFout)) isAccepted=kFALSE; | |
145 | Float_t dvertex[2],cvertex[3]; | |
146 | track->GetImpactParametersTPC(dvertex,cvertex); | |
147 | if (TMath::Abs(dvertex[0]/TMath::Sqrt(cvertex[0]+0.01))>20) isAccepted=kFALSE; | |
148 | if (TMath::Abs(dvertex[1]/TMath::Sqrt(TMath::Abs(cvertex[2]+0.01)))>20) isAccepted=kFALSE; | |
149 | track->GetImpactParameters(dvertex,cvertex); | |
150 | if (TMath::Abs(dvertex[0]/TMath::Sqrt(cvertex[0]+0.01))>10) isAccepted=kFALSE; | |
151 | if (TMath::Abs(dvertex[1]/TMath::Sqrt(TMath::Abs(cvertex[2]+0.01)))>10) isAccepted=kFALSE; | |
152 | if (!isAccepted) isOK=kFALSE; | |
153 | } | |
154 | ||
155 | if (!isOK) continue; | |
156 | // | |
157 | ||
158 | TObject *calibObject; | |
159 | AliTPCseed *seed = 0; | |
160 | for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) { | |
161 | if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break; | |
162 | } | |
163 | if (!seed) continue; | |
164 | if (!fHPTTree) { | |
165 | fHPTTree = new TTree("HPT","HPT"); | |
166 | fHPTTree->SetDirectory(0); | |
167 | } | |
168 | if (fHPTTree->GetEntries()==0){ | |
169 | // | |
170 | fHPTTree->SetDirectory(0); | |
171 | fHPTTree->Branch("t.",&track); | |
172 | fHPTTree->Branch("ft.",&friendTrack); | |
173 | fHPTTree->Branch("s.",&seed); | |
174 | }else{ | |
175 | fHPTTree->SetBranchAddress("t.",&track); | |
176 | fHPTTree->SetBranchAddress("ft.",&friendTrack); | |
177 | fHPTTree->SetBranchAddress("s.",&seed); | |
178 | } | |
179 | fHPTTree->Fill(); | |
180 | // | |
181 | } | |
182 | } | |
183 | ||
184 | ||
185 | ||
186 | void AliTPCcalibV0::DumpToTree(AliESDEvent *esd){ | |
187 | // | |
188 | // Dump V0s fith full firend information to the | |
189 | // | |
190 | Int_t nV0s = fESD->GetNumberOfV0s(); | |
191 | const Int_t kMinCluster=110; | |
192 | const Double_t kDownscale=0.01; | |
46e89793 | 193 | const Float_t kMinPt =1.0; |
194 | const Float_t kMinMinPt =0.7; | |
195 | AliESDfriend *esdFriend=static_cast<AliESDfriend*>(esd->FindListObject("AliESDfriend")); | |
196 | if (!esdFriend) { | |
197 | Printf("ERROR: esdFriend not available"); | |
198 | return; | |
199 | } | |
200 | // | |
201 | ||
202 | for (Int_t ivertex=0; ivertex<nV0s; ivertex++){ | |
203 | Bool_t isOK=kFALSE; | |
204 | AliESDv0 * v0 = (AliESDv0*) esd->GetV0(ivertex); | |
205 | AliESDtrack * track0 = fESD->GetTrack(v0->GetIndex(0)); // negative track | |
206 | AliESDtrack * track1 = fESD->GetTrack(v0->GetIndex(1)); // positive track | |
207 | if (track0->GetTPCNcls()<kMinCluster) continue; | |
208 | if (track0->GetKinkIndex(0)>0) continue; | |
209 | if (track1->GetTPCNcls()<kMinCluster) continue; | |
210 | if (track1->GetKinkIndex(0)>0) continue; | |
211 | if (v0->GetOnFlyStatus()==kFALSE) continue; | |
212 | // | |
213 | if (TMath::Min(track0->Pt(),track1->Pt())<kMinMinPt) continue; | |
214 | // | |
215 | // | |
216 | if (TMath::Max(track0->Pt(),track1->Pt())>kMinPt) isOK=kTRUE; | |
217 | if (gRandom->Rndm()<kDownscale) isOK=kTRUE; | |
218 | if (!isOK) continue; | |
219 | // | |
220 | AliESDfriendTrack *ftrack0 = esdFriend->GetTrack(v0->GetIndex(0)); | |
221 | if (!ftrack0) continue; | |
222 | AliESDfriendTrack *ftrack1 = esdFriend->GetTrack(v0->GetIndex(1)); | |
223 | if (!ftrack1) continue; | |
224 | // | |
225 | TObject *calibObject; | |
226 | AliTPCseed *seed0 = 0; | |
227 | AliTPCseed *seed1 = 0; | |
228 | for (Int_t l=0;(calibObject=ftrack0->GetCalibObject(l));++l) { | |
229 | if ((seed0=dynamic_cast<AliTPCseed*>(calibObject))) break; | |
230 | } | |
231 | for (Int_t l=0;(calibObject=ftrack1->GetCalibObject(l));++l) { | |
232 | if ((seed1=dynamic_cast<AliTPCseed*>(calibObject))) break; | |
233 | } | |
234 | if (!seed0) continue; | |
235 | if (!seed1) continue; | |
236 | AliExternalTrackParam * paramIn0 = (AliExternalTrackParam *)track0->GetInnerParam(); | |
237 | AliExternalTrackParam * paramIn1 = (AliExternalTrackParam *)track1->GetInnerParam(); | |
238 | if (!paramIn0) continue; | |
239 | if (!paramIn1) continue; | |
240 | // | |
241 | // | |
242 | if (!fV0Tree) { | |
243 | fV0Tree = new TTree("V0s","V0s"); | |
244 | fV0Tree->SetDirectory(0); | |
245 | } | |
246 | if (fV0Tree->GetEntries()==0){ | |
247 | // | |
248 | fV0Tree->SetDirectory(0); | |
249 | fV0Tree->Branch("v0.",&v0); | |
250 | fV0Tree->Branch("t0.",&track0); | |
251 | fV0Tree->Branch("t1.",&track1); | |
252 | fV0Tree->Branch("ft0.",&ftrack0); | |
253 | fV0Tree->Branch("ft1.",&ftrack1); | |
254 | fV0Tree->Branch("s0.",&seed0); | |
255 | fV0Tree->Branch("s1.",&seed1); | |
256 | }else{ | |
257 | fV0Tree->SetBranchAddress("v0.",&v0); | |
258 | fV0Tree->SetBranchAddress("t0.",&track0); | |
259 | fV0Tree->SetBranchAddress("t1.",&track1); | |
260 | fV0Tree->SetBranchAddress("ft0.",&ftrack0); | |
261 | fV0Tree->SetBranchAddress("ft1.",&ftrack1); | |
262 | fV0Tree->SetBranchAddress("s0.",&seed0); | |
263 | fV0Tree->SetBranchAddress("s1.",&seed1); | |
264 | } | |
265 | fV0Tree->Fill(); | |
266 | } | |
267 | } | |
268 | ||
269 | ||
270 | Long64_t AliTPCcalibV0::Merge(TCollection *const li) { | |
271 | ||
272 | TIterator* iter = li->MakeIterator(); | |
273 | AliTPCcalibV0* cal = 0; | |
274 | ||
275 | while ((cal = (AliTPCcalibV0*)iter->Next())) { | |
276 | if (cal->fV0Tree){ | |
277 | if (!fV0Tree) { | |
278 | fV0Tree = new TTree("V0s","V0s"); | |
279 | fV0Tree->SetDirectory(0); | |
280 | } | |
281 | if (cal->fV0Tree->GetEntries()>0) AliTPCcalibV0::AddTree(cal->fV0Tree); | |
282 | if (cal->fHPTTree->GetEntries()>0) AliTPCcalibV0::AddTreeHPT(cal->fHPTTree); | |
283 | } | |
284 | } | |
285 | return 0; | |
286 | } | |
287 | ||
288 | ||
289 | void AliTPCcalibV0::AddTree(TTree * treeInput){ | |
290 | // | |
291 | // Add the content of tree: | |
292 | // Notice automatic copy of tree in ROOT does not work for such complicated tree | |
293 | // | |
294 | AliESDv0 * v0 = new AliESDv0; | |
295 | Double_t kMinPt=0.8; | |
296 | AliESDtrack * track0 = 0; // negative track | |
297 | AliESDtrack * track1 = 0; // positive track | |
298 | AliESDfriendTrack *ftrack0 = 0; | |
299 | AliESDfriendTrack *ftrack1 = 0; | |
300 | AliTPCseed *seed0 = 0; | |
301 | AliTPCseed *seed1 = 0; | |
302 | treeInput->SetBranchStatus("ft0.",kFALSE); | |
303 | treeInput->SetBranchStatus("ft1.",kFALSE); | |
304 | TDatabasePDG pdg; | |
305 | Double_t massK0= pdg.GetParticle("K0")->Mass(); | |
306 | Double_t massLambda= pdg.GetParticle("Lambda0")->Mass(); | |
307 | ||
308 | Int_t entries= treeInput->GetEntries(); | |
309 | for (Int_t i=0; i<entries; i++){ | |
310 | treeInput->SetBranchAddress("v0.",&v0); | |
311 | treeInput->SetBranchAddress("t0.",&track0); | |
312 | treeInput->SetBranchAddress("t1.",&track1); | |
313 | treeInput->SetBranchAddress("ft0.",&ftrack0); | |
314 | treeInput->SetBranchAddress("ft1.",&ftrack1); | |
315 | treeInput->SetBranchAddress("s0.",&seed0); | |
316 | treeInput->SetBranchAddress("s1.",&seed1); | |
317 | if (fV0Tree->GetEntries()==0){ | |
318 | fV0Tree->SetDirectory(0); | |
319 | fV0Tree->Branch("v0.",&v0); | |
320 | fV0Tree->Branch("t0.",&track0); | |
321 | fV0Tree->Branch("t1.",&track1); | |
322 | fV0Tree->Branch("ft0.",&ftrack0); | |
323 | fV0Tree->Branch("ft1.",&ftrack1); | |
324 | fV0Tree->Branch("s0.",&seed0); | |
325 | fV0Tree->Branch("s1.",&seed1); | |
326 | }else{ | |
327 | fV0Tree->SetBranchAddress("v0.",&v0); | |
328 | fV0Tree->SetBranchAddress("t0.",&track0); | |
329 | fV0Tree->SetBranchAddress("t1.",&track1); | |
330 | fV0Tree->SetBranchAddress("ft0.",&ftrack0); | |
331 | fV0Tree->SetBranchAddress("ft1.",&ftrack1); | |
332 | fV0Tree->SetBranchAddress("s0.",&seed0); | |
333 | fV0Tree->SetBranchAddress("s1.",&seed1); | |
334 | } | |
335 | // | |
336 | treeInput->GetEntry(i); | |
f11122f7 | 337 | //ftrack0->GetCalibContainer()->SetOwner(kTRUE); |
338 | //ftrack1->GetCalibContainer()->SetOwner(kTRUE); | |
46e89793 | 339 | Bool_t isOK=kTRUE; |
340 | if (v0->GetOnFlyStatus()==kFALSE) isOK=kFALSE; | |
341 | if (track0->GetTPCncls()<100) isOK=kFALSE; | |
342 | if (track1->GetTPCncls()<100) isOK=kFALSE; | |
343 | if (TMath::Min(seed0->Pt(),seed1->Pt())<kMinPt) isOK=kFALSE; | |
344 | if (TMath::Min(track0->Pt(),track1->Pt())<kMinPt) isOK=kFALSE; | |
345 | Bool_t isV0=kFALSE; | |
346 | if (TMath::Abs(v0->GetEffMass(2,2)-massK0)<0.05) isV0=kTRUE; | |
347 | if (TMath::Abs(v0->GetEffMass(4,2)-massLambda)<0.05) isV0=kTRUE; | |
348 | if (TMath::Abs(v0->GetEffMass(2,4)-massLambda)<0.05) isV0=kTRUE; | |
349 | if (TMath::Abs(v0->GetEffMass(0,0))<0.02) isV0=kFALSE; //reject electrons | |
350 | if (!isV0) isOK=kFALSE; | |
351 | if (isOK) fV0Tree->Fill(); | |
352 | delete v0; | |
353 | delete track0; | |
354 | delete track1; | |
355 | delete ftrack0; | |
356 | delete ftrack1; | |
357 | delete seed0; | |
358 | delete seed1; | |
359 | v0=0; | |
360 | track0=0; | |
361 | track1=0; | |
362 | ftrack0=0; | |
363 | ftrack1=0; | |
364 | seed0=0; | |
365 | seed1=0; | |
366 | } | |
367 | } | |
368 | ||
369 | void AliTPCcalibV0::AddTreeHPT(TTree * treeInput){ | |
370 | // | |
371 | // Add the content of tree: | |
372 | // Notice automatic copy of tree in ROOT does not work for such complicated tree | |
373 | // | |
374 | AliESDtrack *track = 0; | |
375 | AliESDfriendTrack *friendTrack = 0; | |
376 | AliTPCseed *seed = 0; | |
377 | if (!treeInput) return; | |
378 | if (treeInput->GetEntries()==0) return; | |
379 | // | |
380 | Int_t entries= treeInput->GetEntries(); | |
381 | // | |
382 | for (Int_t i=0; i<entries; i++){ | |
383 | track=0; | |
384 | friendTrack=0; | |
385 | seed=0; | |
386 | // | |
387 | treeInput->SetBranchAddress("t.",&track); | |
388 | treeInput->SetBranchAddress("ft.",&friendTrack); | |
389 | treeInput->SetBranchAddress("s.",&seed); | |
390 | treeInput->GetEntry(i); | |
391 | // | |
392 | if (fHPTTree->GetEntries()==0){ | |
393 | fHPTTree->SetDirectory(0); | |
394 | fHPTTree->Branch("t.",&track); | |
395 | fHPTTree->Branch("ft.",&friendTrack); | |
396 | fHPTTree->Branch("s.",&seed); | |
397 | }else{ | |
398 | fHPTTree->SetBranchAddress("t.",&track); | |
399 | fHPTTree->SetBranchAddress("ft.",&friendTrack); | |
400 | fHPTTree->SetBranchAddress("s.",&seed); | |
401 | } | |
402 | Bool_t isOK=kTRUE; | |
403 | if (!track->IsOn(AliESDtrack::kITSrefit)) isOK=kFALSE; | |
404 | if (!track->IsOn(AliESDtrack::kTOFout)) isOK=kFALSE; | |
405 | if (isOK) fHPTTree->Fill(); | |
406 | // | |
407 | delete track; | |
408 | delete friendTrack; | |
409 | delete seed; | |
410 | } | |
411 | } | |
412 | ||
413 | ||
10757ee9 | 414 | void AliTPCcalibV0::MakeMC(){ |
415 | // | |
416 | // MC comparison | |
417 | // 1. Select interesting particles | |
418 | // 2. Assign the recosntructed particles | |
419 | // | |
420 | //1. Select interesting particles | |
421 | const Float_t kMinP = 0.2; | |
422 | const Float_t kMinPt = 0.1; | |
423 | const Float_t kMaxR = 0.5; | |
424 | const Float_t kMaxTan = 1.2; | |
425 | const Float_t kMaxRad = 150; | |
426 | // | |
427 | if (!fParticles) fParticles = new TObjArray; | |
428 | TParticle *part=0; | |
429 | // | |
430 | Int_t entries = fStack->GetNtrack(); | |
431 | for (Int_t ipart=0; ipart<entries; ipart++){ | |
432 | part = fStack->Particle(ipart); | |
433 | if (!part) continue; | |
434 | if (part->P()<kMinP) continue; | |
435 | if (part->R()>kMaxR) continue; | |
436 | if (TMath::Abs(TMath::Tan(part->Theta()-TMath::Pi()*0.5))>kMaxTan) continue; | |
437 | Bool_t isInteresting = kFALSE; | |
438 | if (part->GetPdgCode()==22) isInteresting =kTRUE; | |
439 | if (part->GetPdgCode()==310) isInteresting =kTRUE; | |
440 | if (part->GetPdgCode()==111) isInteresting =kTRUE; | |
441 | if (TMath::Abs(part->GetPdgCode()==3122)) isInteresting =kTRUE; | |
442 | ||
443 | // | |
444 | if (!isInteresting) continue; | |
445 | fParticles->AddLast(new TParticle(*part)); | |
446 | } | |
447 | if (fParticles->GetEntries()<1) { | |
448 | return; | |
449 | } | |
450 | // | |
451 | // | |
452 | // | |
453 | Int_t sentries=fParticles->GetEntries();; | |
454 | for (Int_t ipart=0; ipart<sentries; ipart++){ | |
57dc06f2 | 455 | part = (TParticle*)fParticles->At(ipart); |
10757ee9 | 456 | TParticle *p0 = 0; |
457 | TParticle *p1 = 0; | |
458 | ||
459 | Int_t nold =0; | |
460 | Int_t nnew =0; | |
461 | Int_t id0 = part->GetDaughter(0); | |
462 | Int_t id1 = part->GetDaughter(1); | |
463 | if (id0>=fStack->GetNtrack() ) id0*=-1; | |
464 | if (id1>=fStack->GetNtrack() ) id1*=-1; | |
465 | Bool_t findable = kTRUE; | |
466 | if (id0<0 || id1<0) findable = kFALSE; | |
467 | Int_t charge =0; | |
468 | if (findable){ | |
469 | p0 = fStack->Particle(id0); | |
470 | if (p0->R()>kMaxRad) findable = kFALSE; | |
471 | if (p0->Pt()<kMinPt) findable = kFALSE; | |
472 | if (p0->Vz()>250) findable= kFALSE; | |
473 | if (TMath::Abs(TMath::Tan(p0->Theta()-TMath::Pi()*0.5))>2) findable=kFALSE; | |
474 | if (fPdg->GetParticle(p0->GetPdgCode())==0) findable =kFALSE; | |
475 | else | |
476 | if (fPdg->GetParticle(p0->GetPdgCode())->Charge()==0) charge++; | |
477 | ||
478 | p1 = fStack->Particle(id1); | |
479 | if (p1->R()>kMaxRad) findable = kFALSE; | |
480 | if (p1->Pt()<kMinPt) findable = kFALSE; | |
481 | if (TMath::Abs(p1->Vz())>250) findable= kFALSE; | |
482 | if (TMath::Abs(TMath::Tan(p1->Theta()-TMath::Pi()*0.5))>2) findable=kFALSE; | |
483 | if (fPdg->GetParticle(p1->GetPdgCode())==0) findable = kFALSE; | |
484 | else | |
485 | if (fPdg->GetParticle(p1->GetPdgCode())->Charge()==0) charge++; | |
486 | ||
487 | } | |
488 | // (*fDebugStream)<<"MC0"<< | |
489 | // "P.="<<part<< | |
490 | // "findable="<<findable<< | |
491 | // "id0="<<id0<< | |
492 | // "id1="<<id1<< | |
493 | // "\n"; | |
494 | if (!findable) continue; | |
495 | Float_t minpt = TMath::Min(p0->Pt(), p1->Pt()); | |
496 | Int_t type=-1; | |
497 | ||
498 | // | |
499 | // | |
500 | AliKFVertex primVtx(*(fESD->GetPrimaryVertex())); | |
501 | for (Int_t ivertex=0; ivertex<fESD->GetNumberOfV0s(); ivertex++){ | |
502 | AliESDv0 * v0 = fESD->GetV0(ivertex); | |
503 | // select coresponding track | |
504 | AliESDtrack * trackN = fESD->GetTrack(v0->GetIndex(0)); | |
505 | if (TMath::Abs(trackN->GetLabel())!=id0 && TMath::Abs(trackN->GetLabel())!=id1) continue; | |
506 | AliESDtrack * trackP = fESD->GetTrack(v0->GetIndex(1)); | |
507 | if (TMath::Abs(trackP->GetLabel())!=id0 && TMath::Abs(trackP->GetLabel())!=id1) continue; | |
508 | TParticle *pn = fStack->Particle(TMath::Abs(trackN->GetLabel())); | |
509 | TParticle *pp = fStack->Particle(TMath::Abs(trackP->GetLabel())); | |
510 | // | |
511 | // | |
512 | if ( v0->GetOnFlyStatus()) nnew++; | |
513 | if (!v0->GetOnFlyStatus()) nold++; | |
514 | if (part->GetPdgCode()==22 && TMath::Abs(pn->GetPdgCode())==11 && TMath::Abs(pp->GetPdgCode())==11) | |
515 | type =1; | |
516 | if (part->GetPdgCode()==310 && TMath::Abs(pn->GetPdgCode())==211 && TMath::Abs(pp->GetPdgCode())==211) | |
517 | type =0; | |
518 | if (part->GetPdgCode()==3122){ | |
519 | if (TMath::Abs(pn->GetPdgCode())==210 ) type=2; | |
520 | else type=3; | |
521 | } | |
522 | AliKFParticle *v0kf = Fit(primVtx,v0,pn->GetPdgCode(),pp->GetPdgCode()); | |
523 | v0kf->SetProductionVertex( primVtx ); | |
524 | AliKFParticle *v0kfc = Fit(primVtx,v0,pn->GetPdgCode(),pp->GetPdgCode()); | |
525 | v0kfc->SetProductionVertex( primVtx ); | |
526 | v0kfc->SetMassConstraint(fPdg->GetParticle(part->GetPdgCode())->Mass()); | |
527 | Float_t chi2 = v0kf->GetChi2(); | |
528 | Float_t chi2C = v0kf->GetChi2(); | |
529 | // | |
530 | // | |
57e4988a | 531 | TTreeSRedirector *cstream = GetDebugStreamer(); |
532 | if (cstream){ | |
533 | (*cstream)<<"MCRC"<< | |
534 | "P.="<<part<< | |
535 | "type="<<type<< | |
536 | "chi2="<<chi2<< | |
537 | "chi2C="<<chi2C<< | |
538 | "minpt="<<minpt<< | |
539 | "id0="<<id0<< | |
540 | "id1="<<id1<< | |
541 | "Pn.="<<pn<< | |
542 | "Pp.="<<pp<< | |
543 | "tn.="<<trackN<< | |
544 | "tp.="<<trackP<< | |
545 | "nold.="<<nold<< | |
546 | "nnew.="<<nnew<< | |
547 | "v0.="<<v0<< | |
548 | "v0kf.="<<v0kf<< | |
549 | "v0kfc.="<<v0kfc<< | |
550 | "\n"; | |
551 | delete v0kf; | |
552 | delete v0kfc; | |
553 | // | |
554 | } | |
555 | ||
556 | if (cstream){ | |
557 | (*cstream)<<"MC"<< | |
558 | "P.="<<part<< | |
559 | "charge="<<charge<< | |
560 | "type="<<type<< | |
561 | "minpt="<<minpt<< | |
562 | "id0="<<id0<< | |
563 | "id1="<<id1<< | |
564 | "P0.="<<p0<< | |
565 | "P1.="<<p1<< | |
566 | "nold="<<nold<< | |
567 | "nnew="<<nnew<< | |
568 | "\n"; | |
569 | } | |
10757ee9 | 570 | } |
57e4988a | 571 | fParticles->Delete(); |
10757ee9 | 572 | } |
10757ee9 | 573 | } |
574 | ||
a8ef8a9c | 575 | void AliTPCcalibV0::MakeFitTreeTrack(const TObjArray * corrArray, Double_t ptCut, Int_t /*run*/){ |
46e89793 | 576 | // |
577 | // Make a fit tree | |
578 | // | |
579 | // 0. Loop over selected tracks | |
580 | // 1. Loop over all transformation - refit the track with and without the | |
581 | // transformtation | |
582 | // 2. Dump the matching paramaeters to the debugStremer | |
583 | // | |
584 | ||
585 | //Connect input | |
586 | const Int_t kMinNcl=120; | |
587 | TFile f("TPCV0Objects.root"); | |
588 | AliTPCcalibV0 *v0TPC = (AliTPCcalibV0*) f.Get("v0TPC"); | |
589 | TTree * treeInput = v0TPC->GetHPTTree(); | |
590 | TTreeSRedirector *pcstream = new TTreeSRedirector("fitHPT.root"); | |
591 | AliESDtrack *track = 0; | |
592 | AliESDfriendTrack *friendTrack = 0; | |
593 | AliTPCseed *seed = 0; | |
594 | if (!treeInput) return; | |
595 | if (treeInput->GetEntries()==0) return; | |
596 | // | |
597 | treeInput->SetBranchAddress("t.",&track); | |
598 | treeInput->SetBranchAddress("ft.",&friendTrack); | |
599 | treeInput->SetBranchAddress("s.",&seed); | |
600 | // | |
601 | Int_t ncorr=0; | |
602 | if (corrArray) ncorr = corrArray->GetEntries(); | |
603 | AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ; | |
f2ebe917 | 604 | // AliTPCParam *param = AliTPCcalibDB::Instance()->GetParameters(); |
605 | // AliGRPObject* grp = AliTPCcalibDB::Instance()->GetGRP(run); | |
606 | // Double_t time=0.5*(grp->GetTimeStart() +grp->GetTimeEnd()); | |
46e89793 | 607 | // |
608 | // | |
609 | // | |
610 | Int_t ntracks= treeInput->GetEntries(); | |
611 | for (Int_t itrack=0; itrack<ntracks; itrack++){ | |
612 | treeInput->GetEntry(itrack); | |
613 | if (!track) continue; | |
614 | if (seed->Pt()<ptCut) continue; | |
615 | if (track->Pt()<ptCut) continue; | |
616 | if (track->GetTPCncls()<kMinNcl) continue; | |
617 | // | |
618 | // Reapply transformation | |
619 | // | |
620 | for (Int_t irow=0; irow<159; irow++){ | |
621 | AliTPCclusterMI *cluster=seed->GetClusterPointer(irow); | |
622 | if (cluster &&cluster->GetX()>10){ | |
623 | Double_t x0[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()}; | |
624 | Int_t index0[1]={cluster->GetDetector()}; | |
625 | transform->Transform(x0,index0,0,1); | |
626 | cluster->SetX(x0[0]); | |
627 | cluster->SetY(x0[1]); | |
628 | cluster->SetZ(x0[2]); | |
629 | // | |
630 | } | |
631 | } | |
632 | // | |
46e89793 | 633 | AliExternalTrackParam* paramInner=0; |
634 | AliExternalTrackParam* paramOuter=0; | |
635 | AliExternalTrackParam* paramIO=0; | |
636 | Bool_t isOK=kTRUE; | |
637 | for (Int_t icorr=-1; icorr<ncorr; icorr++){ | |
638 | // | |
639 | AliTPCCorrection *corr = 0; | |
640 | if (icorr>=0) corr = (AliTPCCorrection*)corrArray->At(icorr); | |
641 | AliExternalTrackParam * trackInner = RefitTrack(seed, corr,85,134,0.1); | |
642 | AliExternalTrackParam * trackIO = RefitTrack(seed, corr,245,85,0.1); | |
643 | AliExternalTrackParam * trackOuter = RefitTrack(seed, corr,245,134,0.1 ); | |
644 | if (trackInner&&trackOuter&&trackIO){ | |
645 | trackOuter->Rotate(trackInner->GetAlpha()); | |
646 | trackOuter->PropagateTo(trackInner->GetX(),AliTracker::GetBz()); | |
647 | if (icorr<0) { | |
648 | paramInner=trackInner; | |
649 | paramOuter=trackOuter; | |
650 | paramIO=trackIO; | |
651 | paramIO->Rotate(seed->GetAlpha()); | |
652 | paramIO->PropagateTo(seed->GetX(),AliTracker::GetBz()); | |
653 | } | |
654 | }else{ | |
655 | isOK=kFALSE; | |
656 | } | |
657 | ||
658 | } | |
659 | if (paramOuter&& paramInner) { | |
660 | // Bool_t isOK=kTRUE; | |
661 | if (paramInner->GetSigmaY2()>0.01) isOK&=kFALSE; | |
662 | if (paramOuter->GetSigmaY2()>0.01) isOK&=kFALSE; | |
663 | if (paramInner->GetSigmaZ2()>0.01) isOK&=kFALSE; | |
664 | if (paramOuter->GetSigmaZ2()>0.01) isOK&=kFALSE; | |
665 | (*pcstream)<<"fit"<< | |
666 | "s.="<<seed<< | |
667 | "io.="<<paramIO<< | |
668 | "pIn.="<<paramInner<< | |
669 | "pOut.="<<paramOuter; | |
670 | } | |
671 | // | |
672 | } | |
673 | delete pcstream; | |
674 | /* | |
675 | .x ~/rootlogon.C | |
676 | Int_t run=117112; | |
677 | .x ../ConfigCalibTrain.C(run) | |
678 | .L ../AddTaskTPCCalib.C | |
679 | ConfigOCDB(run) | |
680 | TFile fexb("../../RegisterCorrectionExB.root"); | |
681 | AliTPCComposedCorrection *cc= (AliTPCComposedCorrection*) fexb.Get("ComposedExB"); | |
682 | cc->Init(); | |
683 | cc->Print("DA"); // Print used correction classes | |
684 | TObjArray *array = cc->GetCorrections() | |
685 | AliTPCcalibV0::MakeFitTreeTrack(array,2,run); | |
686 | ||
687 | */ | |
688 | } | |
10757ee9 | 689 | |
46e89793 | 690 | void AliTPCcalibV0::MakeFitTreeV0(const TObjArray * corrArray, Double_t ptCut, Int_t run){ |
691 | // | |
692 | // Make a fit tree | |
693 | // | |
694 | // 0. Loop over selected tracks | |
695 | // 1. Loop over all transformation - refit the track with and without the | |
696 | // transformtation | |
697 | // 2. Dump the matching paramaeters to the debugStremer | |
698 | // | |
699 | ||
700 | //Connect input | |
46e89793 | 701 | TFile f("TPCV0Objects.root"); |
702 | AliTPCcalibV0 *v0TPC = (AliTPCcalibV0*) f.Get("v0TPC"); | |
703 | TTree * treeInput = v0TPC->GetV0Tree(); | |
704 | TTreeSRedirector *pcstream = new TTreeSRedirector("fitV0.root"); | |
705 | AliESDv0 *v0 = 0; | |
706 | AliESDtrack *track0 = 0; | |
707 | AliESDfriendTrack *friendTrack0 = 0; | |
708 | AliTPCseed *seed0 = 0; | |
709 | AliTPCseed *s0 = 0; | |
710 | AliESDtrack *track1 = 0; | |
711 | AliESDfriendTrack *friendTrack1 = 0; | |
712 | AliTPCseed *s1 = 0; | |
713 | AliTPCseed *seed1 = 0; | |
714 | if (!treeInput) return; | |
715 | if (treeInput->GetEntries()==0) return; | |
716 | // | |
717 | treeInput->SetBranchAddress("v0.",&v0); | |
718 | treeInput->SetBranchAddress("t0.",&track0); | |
719 | treeInput->SetBranchAddress("ft0.",&friendTrack0); | |
720 | treeInput->SetBranchAddress("s0.",&s0); | |
721 | treeInput->SetBranchAddress("t1.",&track1); | |
722 | treeInput->SetBranchAddress("ft1.",&friendTrack1); | |
723 | treeInput->SetBranchAddress("s1.",&s1); | |
724 | // | |
725 | TDatabasePDG pdg; | |
726 | Int_t ncorr=0; | |
727 | if (corrArray) ncorr = corrArray->GetEntries(); | |
728 | AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ; | |
46e89793 | 729 | Double_t massK0= pdg.GetParticle("K0")->Mass(); |
730 | Double_t massLambda= pdg.GetParticle("Lambda0")->Mass(); | |
731 | Double_t massPion=pdg.GetParticle("pi+")->Mass(); | |
732 | Double_t massProton=pdg.GetParticle("proton")->Mass(); | |
f2ebe917 | 733 | Int_t pdgPion=pdg.GetParticle("pi+")->PdgCode(); |
734 | Int_t pdgProton=pdg.GetParticle("proton")->PdgCode(); | |
735 | Double_t rmass0=0; | |
736 | Double_t rmass1=0; | |
46e89793 | 737 | Double_t massV0=0; |
738 | Int_t pdg0=0; | |
739 | Int_t pdg1=0; | |
740 | // | |
741 | // | |
742 | // | |
743 | Int_t nv0s= treeInput->GetEntries(); | |
744 | for (Int_t iv0=0; iv0<nv0s; iv0++){ | |
745 | Int_t v0Type=0; | |
746 | Int_t isK0=0; | |
747 | Int_t isLambda=0; | |
748 | Int_t isAntiLambda=0; | |
749 | treeInput->GetEntry(iv0); | |
750 | if (TMath::Abs(v0->GetEffMass(2,2)-massK0)<0.03) {isK0=1; v0Type=1;} //select K0s | |
751 | if (TMath::Abs(v0->GetEffMass(4,2)-massLambda)<0.01) {isLambda=1; v0Type=2;} //select Lambda | |
752 | if (TMath::Abs(v0->GetEffMass(2,4)-massLambda)<0.01) {isAntiLambda=1;v0Type=3;} //select Anti Lambda | |
753 | if (isK0+isLambda+isAntiLambda!=1) continue; | |
f2ebe917 | 754 | rmass0=massPion; |
755 | rmass1=massPion; | |
46e89793 | 756 | pdg0=pdgPion; |
757 | pdg1=pdgPion; | |
f2ebe917 | 758 | if (isLambda) {rmass0=massProton; pdg0=pdgProton;} |
759 | if (isAntiLambda) {rmass1=massProton; pdg1=pdgProton;} | |
46e89793 | 760 | massV0=massK0; |
761 | if (isK0==0) massV0=massLambda; | |
762 | // | |
763 | if (!s0) continue; | |
764 | seed0=(s0->GetSign()>0)?s0:s1; | |
765 | seed1=(s0->GetSign()>0)?s1:s0; | |
766 | if (seed0->GetZ()*seed1->GetZ()<0) continue; //remove membrane crossed tracks | |
767 | if (seed0->Pt()<ptCut) continue; | |
768 | if (seed1->Pt()<ptCut) continue; | |
769 | // | |
770 | // Reapply transformation | |
771 | // | |
772 | for (Int_t itype=0; itype<2; itype++){ | |
773 | AliTPCseed * seed = (itype==0) ? seed0: seed1; | |
774 | for (Int_t irow=0; irow<159; irow++){ | |
775 | AliTPCclusterMI *cluster=seed->GetClusterPointer(irow); | |
776 | if (cluster &&cluster->GetX()>10){ | |
777 | Double_t x0[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()}; | |
778 | Int_t index0[1]={cluster->GetDetector()}; | |
779 | transform->Transform(x0,index0,0,1); | |
780 | cluster->SetX(x0[0]); | |
781 | cluster->SetY(x0[1]); | |
782 | cluster->SetZ(x0[2]); | |
783 | // | |
784 | } | |
785 | } | |
786 | } | |
787 | Bool_t isOK=kTRUE; | |
788 | Double_t radius = v0->GetRr(); | |
789 | Double_t xyz[3]; | |
790 | v0->GetXYZ(xyz[0],xyz[1],xyz[2]); | |
791 | Double_t alpha = TMath::ATan2(xyz[1],xyz[0]); | |
792 | TObjArray arrayV0in(ncorr+1); | |
793 | TObjArray arrayV0io(ncorr+1); | |
794 | TObjArray arrayT0(ncorr+1); | |
795 | TObjArray arrayT1(ncorr+1); | |
796 | arrayV0in.SetOwner(kTRUE); | |
797 | arrayV0io.SetOwner(kTRUE); | |
798 | // | |
799 | for (Int_t icorr=-1; icorr<ncorr; icorr++){ | |
800 | AliTPCCorrection *corr =0; | |
801 | if (icorr>=0) corr = (AliTPCCorrection*)corrArray->At(icorr); | |
802 | // | |
f2ebe917 | 803 | AliExternalTrackParam * trackInner0 = RefitTrack(seed0, corr,160,85,rmass0); |
804 | AliExternalTrackParam * trackIO0 = RefitTrack(seed0, corr,245,85,rmass0); | |
805 | AliExternalTrackParam * trackInner1 = RefitTrack(seed1, corr,160,85,rmass1); | |
806 | AliExternalTrackParam * trackIO1 = RefitTrack(seed1, corr,245,85,rmass1); | |
46e89793 | 807 | if (!trackInner0) isOK=kFALSE; |
808 | if (!trackInner1) isOK=kFALSE; | |
809 | if (!trackIO0) isOK=kFALSE; | |
810 | if (!trackIO1) isOK=kFALSE; | |
811 | if (isOK){ | |
812 | if (!trackInner0->Rotate(alpha)) isOK=kFALSE; | |
813 | if (!trackInner1->Rotate(alpha)) isOK=kFALSE; | |
814 | if (!trackIO0->Rotate(alpha)) isOK=kFALSE; | |
815 | if (!trackIO1->Rotate(alpha)) isOK=kFALSE; | |
816 | // | |
f2ebe917 | 817 | if (!AliTracker::PropagateTrackToBxByBz(trackInner0, radius, rmass0, 1, kFALSE)) isOK=kFALSE; |
818 | if (!AliTracker::PropagateTrackToBxByBz(trackInner1, radius, rmass1, 1, kFALSE)) isOK=kFALSE; | |
819 | if (!AliTracker::PropagateTrackToBxByBz(trackIO0, radius, rmass0, 1, kFALSE)) isOK=kFALSE; | |
820 | if (!AliTracker::PropagateTrackToBxByBz(trackIO1, radius, rmass1, 1, kFALSE)) isOK=kFALSE; | |
46e89793 | 821 | if (!isOK) continue; |
822 | arrayT0.AddAt(trackIO0->Clone(),icorr+1); | |
823 | arrayT1.AddAt(trackIO1->Clone(),icorr+1); | |
f2ebe917 | 824 | Int_t charge=TMath::Nint(trackIO0->GetSign()); |
46e89793 | 825 | AliKFParticle pin0( *trackInner0, pdg0*charge); |
826 | AliKFParticle pin1( *trackInner1, -pdg1*charge); | |
827 | AliKFParticle pio0( *trackIO0, pdg0*charge); | |
828 | AliKFParticle pio1( *trackIO1, -pdg1*charge); | |
829 | AliKFParticle v0in; | |
830 | AliKFParticle v0io; | |
831 | v0in+=pin0; | |
832 | v0in+=pin1; | |
833 | v0io+=pio0; | |
834 | v0io+=pio1; | |
835 | arrayV0in.AddAt(v0in.Clone(),icorr+1); | |
836 | arrayV0io.AddAt(v0io.Clone(),icorr+1); | |
837 | } | |
838 | } | |
839 | if (!isOK) continue; | |
840 | // | |
841 | //AliKFParticle* pin0= (AliKFParticle*)arrayV0in.At(0); | |
842 | AliKFParticle* pio0= (AliKFParticle*)arrayV0io.At(0); | |
843 | AliExternalTrackParam *param0=(AliExternalTrackParam *)arrayT0.At(0); | |
844 | AliExternalTrackParam *param1=(AliExternalTrackParam *)arrayT1.At(0); | |
845 | Double_t mass0=0, mass0E=0; | |
846 | pio0->GetMass( mass0,mass0E); | |
847 | // | |
848 | Double_t mean=mass0-massV0; | |
849 | if (TMath::Abs(mean)>0.05) continue; | |
850 | Double_t mass[10000]; | |
851 | // | |
852 | Int_t dtype=30; // id for V0 | |
853 | Int_t ptype=5; // id for invariant mass | |
854 | // Int_t id=TMath::Nint(100.*(param0->Pt()-param1->Pt())/(param0->Pt()+param1->Pt())); // K0s V0 asymetry | |
f2ebe917 | 855 | Int_t id=Int_t(1000.*(param0->Pt()-param1->Pt())); // K0s V0 asymetry |
46e89793 | 856 | Double_t gx,gy,gz, px,py,pz; |
857 | Double_t pt = v0->Pt(); | |
858 | v0->GetXYZ(gx,gy,gz); | |
859 | v0->GetPxPyPz(px,py,pz); | |
860 | Double_t theta=pz/TMath::Sqrt(px*px+py*py); | |
861 | Double_t phi=TMath::ATan2(py,px); | |
862 | Double_t snp=0.5*(seed0->GetSnp()+seed1->GetSnp()); | |
863 | Double_t sector=9*phi/TMath::Pi(); | |
864 | Double_t dsec=sector-TMath::Nint(sector); | |
865 | Double_t refX=TMath::Sqrt(gx*gx+gy*gy); | |
866 | //Int_t nentries=v0Type; | |
867 | Double_t bz=AliTracker::GetBz(); | |
868 | Double_t dRrec=0; | |
869 | (*pcstream)<<"fitDebug"<< | |
870 | "id="<<id<< | |
871 | "v0.="<<v0<< | |
872 | "mean="<<mean<< | |
873 | "rms="<<mass0E<< | |
874 | "pio0.="<<pio0<< | |
875 | "p0.="<<param0<< | |
876 | "p1.="<<param1; | |
877 | (*pcstream)<<"fit"<< // dump valus for fit | |
878 | "run="<<run<< //run number | |
879 | "bz="<<bz<< // magnetic filed used | |
880 | "dtype="<<dtype<< // detector match type 30 | |
881 | "ptype="<<ptype<< // parameter type | |
882 | "theta="<<theta<< // theta | |
883 | "phi="<<phi<< // phi | |
884 | "snp="<<snp<< // snp | |
885 | "mean="<<mean<< // mean dist value | |
886 | "rms="<<mass0E<< // rms | |
887 | "sector="<<sector<< | |
888 | "dsec="<<dsec<< | |
889 | // | |
890 | "refX="<<refX<< // reference radius | |
891 | "gx="<<gx<< // global position | |
892 | "gy="<<gy<< // global position | |
893 | "gz="<<gz<< // global position | |
894 | "dRrec="<<dRrec<< // delta Radius in reconstruction | |
895 | "pt="<<pt<< // pt of the particle | |
896 | "id="<<id<< //delta of the momenta | |
897 | "entries="<<v0Type;// type of the V0 | |
898 | for (Int_t icorr=0; icorr<ncorr; icorr++){ | |
899 | AliTPCCorrection *corr =0; | |
900 | if (icorr>=0) corr = (AliTPCCorrection*)corrArray->At(icorr); | |
f2ebe917 | 901 | // AliKFParticle* pin= (AliKFParticle*)arrayV0in.At(icorr+1); |
46e89793 | 902 | AliKFParticle* pio= (AliKFParticle*)arrayV0io.At(icorr+1); |
903 | AliExternalTrackParam *par0=(AliExternalTrackParam *)arrayT0.At(icorr+1); | |
904 | AliExternalTrackParam *par1=(AliExternalTrackParam *)arrayT1.At(icorr+1); | |
905 | Double_t massE=0; | |
906 | pio->GetMass( mass[icorr],massE); | |
907 | mass[icorr]-=mass0; | |
908 | (*pcstream)<<"fit"<< | |
909 | Form("%s=",corr->GetName())<<mass[icorr]; | |
910 | (*pcstream)<<"fitDebug"<< | |
911 | Form("%s=",corr->GetName())<<mass[icorr]<< | |
912 | Form("%sp0.=",corr->GetName())<<par0<< | |
913 | Form("%sp1=",corr->GetName())<<par1; | |
914 | } | |
915 | (*pcstream)<<"fit"<< "isOK="<<isOK<<"\n"; | |
916 | (*pcstream)<<"fitDebug"<< "isOK="<<isOK<<"\n"; | |
917 | } | |
918 | delete pcstream; | |
919 | } | |
10757ee9 | 920 | |
46e89793 | 921 | AliExternalTrackParam * AliTPCcalibV0::RefitTrack(AliTPCseed *seed, AliTPCCorrection * corr, Double_t xstart, Double_t xstop, Double_t mass){ |
10757ee9 | 922 | // |
46e89793 | 923 | // Refit the track: |
924 | // seed - tpc track with cluster | |
925 | // corr - distrotion/correction class - apllied to the points | |
926 | // xstart - radius to start propagate/update | |
927 | // xstop - radius to stop propagate/update | |
928 | // | |
929 | const Double_t kResetCov=20.; | |
930 | const Double_t kSigma=5.; | |
931 | Double_t covar[15]; | |
932 | for (Int_t i=0;i<15;i++) covar[i]=0; | |
933 | covar[0]=kSigma*kSigma; | |
934 | covar[2]=kSigma*kSigma; | |
935 | covar[5]=kSigma*kSigma/Float_t(150*150); | |
936 | covar[9]=kSigma*kSigma/Float_t(150*150); | |
937 | covar[14]=1*1; | |
938 | // | |
939 | AliExternalTrackParam *refit = new AliExternalTrackParam(*seed); | |
940 | refit->PropagateTo(xstart, AliTracker::GetBz()); | |
941 | refit->AddCovariance(covar); | |
942 | refit->ResetCovariance(kResetCov); | |
943 | Double_t xmin = TMath::Min(xstart,xstop); | |
944 | Double_t xmax = TMath::Max(xstart,xstop); | |
945 | Int_t ncl=0; | |
10757ee9 | 946 | // |
46e89793 | 947 | Bool_t isOK=kTRUE; |
948 | for (Int_t index0=0; index0<159; index0++){ | |
949 | Int_t irow= (xstart<xstop)? index0:159-index0; | |
950 | AliTPCclusterMI *cluster=seed->GetClusterPointer(irow); //cluster in local system | |
951 | if (!cluster) continue; | |
952 | if (cluster->GetX()<xmin) continue; | |
953 | if (cluster->GetX()>xmax) continue; | |
954 | Double_t alpha = TMath::Pi()*(cluster->GetDetector()+0.5)/9.; | |
955 | if (!refit->Rotate(alpha)) isOK=kFALSE; | |
956 | Double_t x = cluster->GetX(); | |
957 | Double_t y = cluster->GetY(); | |
958 | Double_t z = cluster->GetZ(); | |
959 | if (corr){ | |
960 | Float_t xyz[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()}; // original position | |
961 | corr->DistortPointLocal(xyz,cluster->GetDetector()); | |
962 | x=xyz[0]; | |
963 | y=xyz[1]; | |
964 | z=xyz[2]; | |
965 | } | |
966 | if (!AliTracker::PropagateTrackToBxByBz(refit, x,mass,1.,kFALSE)) isOK=kFALSE; | |
967 | if (!isOK) continue; | |
968 | Double_t cov[3]={0.01,0.,0.01}; | |
969 | Double_t yz[2]={y,z}; | |
970 | if (!refit->Update(yz,cov)) isOK=kFALSE; | |
971 | ncl++; | |
972 | } | |
973 | if (!AliTracker::PropagateTrackToBxByBz(refit, xstop, mass,1.,kTRUE)) isOK=kFALSE; | |
974 | // | |
975 | if (!isOK) { | |
976 | delete refit; | |
977 | return 0; | |
978 | } | |
979 | return refit; | |
980 | } | |
981 | ||
982 | ||
983 | ||
984 | ||
985 | ||
986 | // | |
987 | // Obsolete part | |
988 | // | |
989 | ||
990 | ||
991 | void AliTPCcalibV0::MakeV0s(){ | |
992 | // | |
993 | // | |
10757ee9 | 994 | // |
995 | const Int_t kMinCluster=40; | |
996 | const Float_t kMinR =0; | |
997 | if (! fV0s) fV0s = new TObjArray(10); | |
998 | fV0s->Clear(); | |
999 | // | |
1000 | // Old V0 finder | |
1001 | // | |
1002 | for (Int_t ivertex=0; ivertex<fESD->GetNumberOfV0s(); ivertex++){ | |
1003 | AliESDv0 * v0 = fESD->GetV0(ivertex); | |
1004 | if (v0->GetOnFlyStatus()) continue; | |
1005 | fV0s->AddLast(v0); | |
1006 | } | |
1007 | ProcessV0(0); | |
1008 | fV0s->Clear(0); | |
1009 | // | |
1010 | // MI V0 finder | |
1011 | // | |
1012 | for (Int_t ivertex=0; ivertex<fESD->GetNumberOfV0s(); ivertex++){ | |
1013 | AliESDv0 * v0 = fESD->GetV0(ivertex); | |
1014 | if (!v0->GetOnFlyStatus()) continue; | |
1015 | fV0s->AddLast(v0); | |
1016 | } | |
1017 | ProcessV0(1); | |
1018 | fV0s->Clear(); | |
1019 | return; | |
1020 | // | |
1021 | // combinatorial | |
1022 | // | |
1023 | Int_t ntracks = fESD->GetNumberOfTracks(); | |
1024 | for (Int_t itrack0=0; itrack0<ntracks; itrack0++){ | |
1025 | AliESDtrack * track0 = fESD->GetTrack(itrack0); | |
1026 | if (track0->GetSign()>0) continue; | |
1027 | if ( track0->GetTPCNcls()<kMinCluster) continue; | |
1028 | if (track0->GetKinkIndex(0)>0) continue; | |
1029 | // | |
1030 | for (Int_t itrack1=0; itrack1<ntracks; itrack1++){ | |
1031 | AliESDtrack * track1 = fESD->GetTrack(itrack1); | |
1032 | if (track1->GetSign()<0) continue; | |
1033 | if ( track1->GetTPCNcls()<kMinCluster) continue; | |
1034 | if (track1->GetKinkIndex(0)>0) continue; | |
1035 | // | |
1036 | // AliExternalTrackParam param0(*track0); | |
1037 | // AliExternalTrackParam param1(*track1); | |
1038 | AliV0 vertex; | |
1039 | vertex.SetParamN(*track0); | |
1040 | vertex.SetParamP(*track1); | |
1041 | Float_t xyz[3]; | |
1042 | xyz[0] = fESD->GetPrimaryVertex()->GetXv(); | |
1043 | xyz[1] = fESD->GetPrimaryVertex()->GetYv(); | |
1044 | xyz[2] = fESD->GetPrimaryVertex()->GetZv(); | |
1045 | vertex.Update(xyz); | |
1046 | if (vertex.GetRr()<kMinR) continue; | |
1047 | if (vertex.GetDcaV0Daughters()>1.) continue; | |
1048 | if (vertex.GetDcaV0Daughters()>0.3*vertex.GetRr()) continue; | |
1049 | // if (vertex.GetPointAngle()<0.9) continue; | |
1050 | vertex.SetIndex(0,itrack0); | |
1051 | vertex.SetIndex(1,itrack1); | |
1052 | fV0s->AddLast(new AliV0(vertex)); | |
1053 | } | |
1054 | } | |
1055 | ProcessV0(2); | |
1056 | for (Int_t i=0;i<fV0s->GetEntries(); i++) delete fV0s->At(i); | |
1057 | fV0s->Clear(); | |
1058 | } | |
1059 | ||
1060 | ||
1061 | ||
1062 | ||
1063 | ||
1064 | ||
10757ee9 | 1065 | |
1066 | ||
1067 | void AliTPCcalibV0::ProcessV0(Int_t ftype){ | |
1068 | // | |
46e89793 | 1069 | // Obsolete |
1070 | // | |
10757ee9 | 1071 | if (! fGammas) fGammas = new TObjArray(10); |
1072 | fGammas->Clear(); | |
1073 | Int_t nV0s = fV0s->GetEntries(); | |
1074 | if (nV0s==0) return; | |
1075 | AliKFVertex primVtx(*(fESD->GetPrimaryVertex())); | |
1076 | // | |
1077 | for (Int_t ivertex=0; ivertex<nV0s; ivertex++){ | |
1078 | AliESDv0 * v0 = (AliESDv0*)fV0s->At(ivertex); | |
1079 | AliESDtrack * trackN = fESD->GetTrack(v0->GetIndex(0)); // negative track | |
1080 | AliESDtrack * trackP = fESD->GetTrack(v0->GetIndex(1)); // positive track | |
1081 | ||
1082 | const AliExternalTrackParam * paramInNeg = trackN->GetInnerParam(); | |
1083 | const AliExternalTrackParam * paramInPos = trackP->GetInnerParam(); | |
1084 | ||
1085 | if (!paramInPos) continue; // in case the inner paramters do not exist | |
1086 | if (!paramInNeg) continue; | |
1087 | // | |
1088 | // | |
1089 | // | |
1090 | AliKFParticle *v0K0 = Fit(primVtx,v0,-211,211); | |
1091 | AliKFParticle *v0Gamma = Fit(primVtx,v0,11,-11); | |
1092 | AliKFParticle *v0Lambda42 = Fit(primVtx,v0,-2212,211); | |
1093 | AliKFParticle *v0Lambda24 = Fit(primVtx,v0,-211,2212); | |
1094 | //Set production vertex | |
1095 | v0K0->SetProductionVertex( primVtx ); | |
1096 | v0Gamma->SetProductionVertex( primVtx ); | |
1097 | v0Lambda42->SetProductionVertex( primVtx ); | |
1098 | v0Lambda24->SetProductionVertex( primVtx ); | |
1099 | Double_t massK0, massGamma, massLambda42,massLambda24, massSigma; | |
1100 | v0K0->GetMass( massK0,massSigma); | |
1101 | v0Gamma->GetMass( massGamma,massSigma); | |
1102 | v0Lambda42->GetMass( massLambda42,massSigma); | |
1103 | v0Lambda24->GetMass( massLambda24,massSigma); | |
1104 | Float_t chi2K0 = v0K0->GetChi2()/v0K0->GetNDF(); | |
1105 | Float_t chi2Gamma = v0Gamma->GetChi2()/v0Gamma->GetNDF(); | |
1106 | Float_t chi2Lambda42 = v0Lambda42->GetChi2()/v0Lambda42->GetNDF(); | |
1107 | Float_t chi2Lambda24 = v0Lambda24->GetChi2()/v0Lambda24->GetNDF(); | |
1108 | // | |
1109 | // Mass Contrained params | |
1110 | // | |
1111 | AliKFParticle *v0K0C = Fit(primVtx,v0,-211,211); | |
1112 | AliKFParticle *v0GammaC = Fit(primVtx,v0,11,-11); | |
1113 | AliKFParticle *v0Lambda42C = Fit(primVtx,v0,-2212,211); //lambdaBar | |
1114 | AliKFParticle *v0Lambda24C = Fit(primVtx,v0,-211,2212); //lambda | |
1115 | // | |
1116 | v0K0C->SetProductionVertex( primVtx ); | |
1117 | v0GammaC->SetProductionVertex( primVtx ); | |
1118 | v0Lambda42C->SetProductionVertex( primVtx ); | |
1119 | v0Lambda24C->SetProductionVertex( primVtx ); | |
1120 | ||
1121 | v0K0C->SetMassConstraint(fPdg->GetParticle(310)->Mass()); | |
1122 | v0GammaC->SetMassConstraint(0); | |
1123 | v0Lambda42C->SetMassConstraint(fPdg->GetParticle(-3122)->Mass()); | |
1124 | v0Lambda24C->SetMassConstraint(fPdg->GetParticle(3122)->Mass()); | |
1125 | // | |
1126 | Double_t timeK0, sigmaTimeK0; | |
1127 | Double_t timeLambda42, sigmaTimeLambda42; | |
1128 | Double_t timeLambda24, sigmaTimeLambda24; | |
1129 | v0K0C->GetLifeTime(timeK0, sigmaTimeK0); | |
1130 | //v0K0Gamma->GetLifeTime(timeK0, sigmaTimeK0); | |
1131 | v0Lambda42C->GetLifeTime(timeLambda42, sigmaTimeLambda42); | |
1132 | v0Lambda24C->GetLifeTime(timeLambda24, sigmaTimeLambda24); | |
1133 | ||
1134 | ||
1135 | // | |
1136 | Float_t chi2K0C = v0K0C->GetChi2()/v0K0C->GetNDF(); | |
1137 | if (chi2K0C<0) chi2K0C=100; | |
1138 | Float_t chi2GammaC = v0GammaC->GetChi2()/v0GammaC->GetNDF(); | |
1139 | if (chi2GammaC<0) chi2GammaC=100; | |
1140 | Float_t chi2Lambda42C = v0Lambda42C->GetChi2()/v0Lambda42C->GetNDF(); | |
1141 | if (chi2Lambda42C<0) chi2Lambda42C=100; | |
1142 | Float_t chi2Lambda24C = v0Lambda24C->GetChi2()/v0Lambda24C->GetNDF(); | |
1143 | if (chi2Lambda24C<0) chi2Lambda24C=100; | |
1144 | // | |
1145 | Float_t minChi2C=99; | |
1146 | Int_t type =-1; | |
1147 | if (chi2K0C<minChi2C) { minChi2C= chi2K0C; type=0;} | |
1148 | if (chi2GammaC<minChi2C) { minChi2C= chi2GammaC; type=1;} | |
1149 | if (chi2Lambda42C<minChi2C) { minChi2C= chi2Lambda42C; type=2;} | |
1150 | if (chi2Lambda24C<minChi2C) { minChi2C= chi2Lambda24C; type=3;} | |
1151 | Float_t minChi2=99; | |
1152 | Int_t type0 =-1; | |
1153 | if (chi2K0<minChi2) { minChi2= chi2K0; type0=0;} | |
1154 | if (chi2Gamma<minChi2) { minChi2= chi2Gamma; type0=1;} | |
1155 | if (chi2Lambda42<minChi2) { minChi2= chi2Lambda42; type0=2;} | |
1156 | if (chi2Lambda24<minChi2) { minChi2= chi2Lambda24; type0=3;} | |
1157 | ||
1158 | // 0 is negative particle; 1 is positive particle | |
1159 | Float_t betaGamma0 = 0; | |
1160 | Float_t betaGamma1 = 0; | |
1161 | ||
1162 | switch (type) { | |
1163 | case 0: | |
1164 | betaGamma0 = paramInNeg->GetP()/fPdg->GetParticle(-211)->Mass(); | |
1165 | betaGamma1 = paramInPos->GetP()/fPdg->GetParticle(211)->Mass(); | |
1166 | break; | |
1167 | case 1: | |
1168 | betaGamma0 = paramInNeg->GetP()/fPdg->GetParticle(11)->Mass(); | |
1169 | betaGamma1 = paramInPos->GetP()/fPdg->GetParticle(-11)->Mass(); | |
1170 | break; | |
1171 | case 2: | |
1172 | betaGamma0 = paramInNeg->GetP()/fPdg->GetParticle(-2212)->Mass(); | |
1173 | betaGamma1 = paramInPos->GetP()/fPdg->GetParticle(211)->Mass(); | |
1174 | break; | |
1175 | case 3: | |
1176 | betaGamma0 = paramInNeg->GetP()/fPdg->GetParticle(-211)->Mass(); | |
1177 | betaGamma1 = paramInPos->GetP()/fPdg->GetParticle(2212)->Mass(); | |
1178 | break; | |
1179 | } | |
1180 | ||
1181 | // cuts and histogram filling | |
1182 | Int_t numCand = 0; // number of particle types which have a chi2 < 10*minChi2 | |
1183 | ||
1184 | if (minChi2C < 2 && ftype == 1) { | |
1185 | // | |
1186 | if (chi2K0C < 10*minChi2C) numCand++; | |
1187 | if (chi2GammaC < 10*minChi2C) numCand++; | |
1188 | if (chi2Lambda42C < 10*minChi2C) numCand++; | |
1189 | if (chi2Lambda24C < 10*minChi2C) numCand++; | |
1190 | // | |
10757ee9 | 1191 | } |
1192 | ||
1193 | // | |
1194 | // | |
1195 | // write output tree | |
1196 | if (minChi2>50) continue; | |
57e4988a | 1197 | TTreeSRedirector *cstream = GetDebugStreamer(); |
1198 | if (cstream){ | |
1199 | (*cstream)<<"V0"<< | |
1200 | "ftype="<<ftype<< | |
1201 | "v0.="<<v0<< | |
1202 | "trackN.="<<trackN<< | |
1203 | "trackP.="<<trackP<< | |
1204 | // | |
1205 | "betaGamma0="<<betaGamma0<< | |
1206 | "betaGamma1="<<betaGamma1<< | |
1207 | // | |
1208 | "type="<<type<< | |
1209 | "chi2C="<<minChi2C<< | |
1210 | "v0K0.="<<v0K0<< | |
1211 | "v0Gamma.="<<v0Gamma<< | |
1212 | "v0Lambda42.="<<v0Lambda42<< | |
1213 | "v0Lambda24.="<<v0Lambda24<< | |
1214 | // | |
1215 | "chi20K0.="<<chi2K0<< | |
1216 | "chi2Gamma.="<<chi2Gamma<< | |
1217 | "chi2Lambda42.="<<chi2Lambda42<< | |
1218 | "chi2Lambda24.="<<chi2Lambda24<< | |
1219 | // | |
1220 | "chi20K0c.="<<chi2K0C<< | |
1221 | "chi2Gammac.="<<chi2GammaC<< | |
1222 | "chi2Lambda42c.="<<chi2Lambda42C<< | |
1223 | "chi2Lambda24c.="<<chi2Lambda24C<< | |
1224 | // | |
1225 | "v0K0C.="<<v0K0C<< | |
1226 | "v0GammaC.="<<v0GammaC<< | |
1227 | "v0Lambda42C.="<<v0Lambda42C<< | |
1228 | "v0Lambda24C.="<<v0Lambda24C<< | |
1229 | // | |
1230 | "massK0="<<massK0<< | |
1231 | "massGamma="<<massGamma<< | |
1232 | "massLambda42="<<massLambda42<< | |
1233 | "massLambda24="<<massLambda24<< | |
1234 | // | |
1235 | "timeK0="<<timeK0<< | |
1236 | "timeLambda42="<<timeLambda42<< | |
1237 | "timeLambda24="<<timeLambda24<< | |
1238 | "\n"; | |
1239 | } | |
10757ee9 | 1240 | if (type==1) fGammas->AddLast(v0); |
1241 | // | |
1242 | // | |
1243 | // | |
1244 | delete v0K0; | |
1245 | delete v0Gamma; | |
1246 | delete v0Lambda42; | |
1247 | delete v0Lambda24; | |
1248 | delete v0K0C; | |
1249 | delete v0GammaC; | |
1250 | delete v0Lambda42C; | |
1251 | delete v0Lambda24C; | |
1252 | } | |
1253 | ProcessPI0(); | |
1254 | } | |
1255 | ||
1256 | ||
1257 | ||
1258 | void AliTPCcalibV0::ProcessPI0(){ | |
1259 | // | |
1260 | // | |
1261 | // | |
1262 | Int_t nentries = fGammas->GetEntries(); | |
1263 | if (nentries<2) return; | |
1264 | // | |
1265 | Double_t m0[3], m1[3]; | |
1266 | AliKFVertex primVtx(*(fESD->GetPrimaryVertex())); | |
1267 | for (Int_t i0=0; i0<nentries; i0++){ | |
1268 | AliESDv0 *v00 = (AliESDv0*)fGammas->At(i0); | |
1269 | v00->GetPxPyPz (m0[0], m0[1], m0[2]); | |
1270 | AliKFParticle *p00 = Fit(primVtx, v00, 11,-11); | |
1271 | p00->SetProductionVertex( primVtx ); | |
1272 | p00->SetMassConstraint(0); | |
1273 | // | |
1274 | for (Int_t i1=i0; i1<nentries; i1++){ | |
1275 | AliESDv0 *v01 = (AliESDv0*)fGammas->At(i1); | |
1276 | v01->GetPxPyPz (m1[0], m1[1], m1[2]); | |
1277 | AliKFParticle *p01 = Fit(primVtx, v01, 11,-11); | |
1278 | p01->SetProductionVertex( primVtx ); | |
1279 | p01->SetMassConstraint(0); | |
1280 | if (v00->GetIndex(0) != v01->GetIndex(0) && | |
1281 | v00->GetIndex(1) != v01->GetIndex(1)){ | |
1282 | AliKFParticle pi0( *p00,*p01); | |
1283 | pi0.SetProductionVertex(primVtx); | |
1284 | Double_t n1 = TMath::Sqrt (m0[0]*m0[0] + m0[1]*m0[1] + m0[2]*m0[2]); | |
1285 | Double_t n2 = TMath::Sqrt (m1[0]*m1[0] + m1[1]*m1[1] + m1[2]*m1[2]); | |
1286 | Double_t mass = TMath::Sqrt(2.*(n1*n2 - (m0[0]*m1[0] + m0[1]*m1[1] + m0[2]*m1[2]))); | |
57e4988a | 1287 | TTreeSRedirector *cstream = GetDebugStreamer(); |
1288 | if (cstream){ | |
1289 | (*cstream)<<"PI0"<< | |
1290 | "v00.="<<v00<< | |
1291 | "v01.="<<v01<< | |
1292 | "mass="<<mass<< | |
1293 | "p00.="<<p00<< | |
1294 | "p01.="<<p01<< | |
1295 | "pi0.="<<&pi0<< | |
1296 | "\n"; | |
1297 | } | |
10757ee9 | 1298 | } |
1299 | delete p01; | |
1300 | } | |
1301 | delete p00; | |
1302 | } | |
1303 | } | |
1304 | ||
1305 | ||
1306 | ||
1307 | ||
1308 | ||
57dc06f2 | 1309 | AliKFParticle * AliTPCcalibV0::Fit(AliKFVertex & /*primVtx*/, AliESDv0 *v0, Int_t PDG1, Int_t PDG2){ |
10757ee9 | 1310 | // |
1311 | // Make KF Particle | |
1312 | // | |
1313 | AliKFParticle p1( *(v0->GetParamN()), PDG1 ); | |
1314 | AliKFParticle p2( *(v0->GetParamP()), PDG2 ); | |
1315 | AliKFParticle *V0 = new AliKFParticle; | |
1316 | Double_t x, y, z; | |
1317 | v0->GetXYZ(x,y,z ); | |
1318 | V0->SetVtxGuess(x,y,z); | |
1319 | *(V0)+=p1; | |
1320 | *(V0)+=p2; | |
1321 | return V0; | |
1322 | } | |
1323 | ||
10757ee9 | 1324 | |
1325 | ||
1326 | ||
1327 | void AliTPCcalibV0::BinLogX(TH2F *h) { | |
1328 | // | |
1329 | // | |
1330 | // | |
1331 | TAxis *axis = h->GetXaxis(); | |
1332 | int bins = axis->GetNbins(); | |
1333 | ||
1334 | Double_t from = axis->GetXmin(); | |
1335 | Double_t to = axis->GetXmax(); | |
1336 | Double_t *new_bins = new Double_t[bins + 1]; | |
1337 | ||
1338 | new_bins[0] = from; | |
1339 | Double_t factor = pow(to/from, 1./bins); | |
1340 | ||
1341 | for (int i = 1; i <= bins; i++) { | |
1342 | new_bins[i] = factor * new_bins[i-1]; | |
1343 | } | |
1344 | axis->Set(bins, new_bins); | |
4ce766eb | 1345 | delete [] new_bins; |
10757ee9 | 1346 | |
1347 | } | |
1348 | ||
1349 | ||
1350 | ||
1351 |