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