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