9 #include "TDatabasePDG.h"
10 #include "AliRunLoader.h"
17 #include "TLinearFitter.h"
19 #include "TTreeStream.h"
25 #include "AliTracker.h"
26 #include "AliESDEvent.h"
27 #include "AliESDtrack.h"
28 #include "AliESDfriend.h"
29 #include "AliESDfriendTrack.h"
30 #include "AliMathBase.h"
31 #include "AliTPCseed.h"
32 #include "AliTPCclusterMI.h"
34 #include "AliKFParticle.h"
35 #include "AliKFVertex.h"
37 #include "AliTrackPointArray.h"
39 #include "AliTPCcalibV0.h"
42 #include "TTreeStream.h"
43 #include "AliTPCcalibDB.h"
44 #include "AliTPCCorrection.h"
45 #include "AliGRPObject.h"
46 #include "AliTPCTransform.h"
47 #include "AliAnalysisManager.h"
51 ClassImp(AliTPCcalibV0)
54 AliTPCcalibV0::AliTPCcalibV0() :
67 AliTPCcalibV0::AliTPCcalibV0(const Text_t *name, const Text_t *title):
78 fPdg = new TDatabasePDG;
79 // create output histograms
84 AliTPCcalibV0::~AliTPCcalibV0(){
96 void AliTPCcalibV0::ProcessESD(AliESDEvent *esd){
101 AliKFParticle::SetField(esd->GetMagneticField());
102 if (TMath::Abs(AliTracker::GetBz())<1) return;
107 void AliTPCcalibV0::DumpToTreeHPT(AliESDEvent *esd){
109 // Dump V0s fith full firend information to the
111 if (TMath::Abs(AliTracker::GetBz())<1) return;
112 const Int_t kMinCluster=110;
113 const Float_t kMinPt =4.;
114 AliESDfriend *esdFriend=static_cast<AliESDfriend*>(esd->FindListObject("AliESDfriend"));
116 // Printf("ERROR: esdFriend not available");
120 Int_t ntracks=esd->GetNumberOfTracks();
121 for (Int_t i=0;i<ntracks;++i) {
123 AliESDtrack *track = esd->GetTrack(i);
124 if (track->GetTPCncls()<kMinCluster) continue;
125 if (TMath::Abs(AliTracker::GetBz())>1){ // cut on momenta if measured
126 if (track->Pt()>kMinPt) isOK=kTRUE;
128 if (TMath::Abs(AliTracker::GetBz())<1){ // require primary track for the B field OFF data
129 Bool_t isAccepted=kTRUE;
130 if (!track->IsOn(AliESDtrack::kITSrefit)) isAccepted=kFALSE;
131 if (!track->IsOn(AliESDtrack::kTPCrefit)) isAccepted=kFALSE;
132 if (!track->IsOn(AliESDtrack::kTOFout)) isAccepted=kFALSE;
133 Float_t dvertex[2],cvertex[3];
134 track->GetImpactParametersTPC(dvertex,cvertex);
135 if (TMath::Abs(dvertex[0]/TMath::Sqrt(cvertex[0]+0.01))>20) isAccepted=kFALSE;
136 if (TMath::Abs(dvertex[1]/TMath::Sqrt(TMath::Abs(cvertex[2]+0.01)))>20) isAccepted=kFALSE;
137 track->GetImpactParameters(dvertex,cvertex);
138 if (TMath::Abs(dvertex[0]/TMath::Sqrt(cvertex[0]+0.01))>10) isAccepted=kFALSE;
139 if (TMath::Abs(dvertex[1]/TMath::Sqrt(TMath::Abs(cvertex[2]+0.01)))>10) isAccepted=kFALSE;
140 if (!isAccepted) isOK=kFALSE;
142 if ( track->GetTPCsignal()>100 && track->GetInnerParam()->Pt()>1 ){
143 if (track->IsOn(AliESDtrack::kITSin)||track->IsOn(AliESDtrack::kTRDout)||track->IsOn(AliESDtrack::kTOFin))
146 TString filename(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
147 Int_t eventNumber = esd->GetEventNumberInFile();
148 Bool_t hasFriend=(esdFriend) ? (esdFriend->GetTrack(i)!=0):0;
149 Bool_t hasITS=(track->GetNcls(0)>2);
150 printf("DUMPIONTrack:%s|%f|%d|%d|%d\n",filename.Data(),track->GetInnerParam()->Pt()*track->GetTPCsignal()/50., eventNumber,hasFriend,hasITS);
154 TString filename(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
155 Int_t eventNumber = esd->GetEventNumberInFile();
156 Bool_t hasFriend=(esdFriend) ? (esdFriend->GetTrack(i)!=0):0;
157 Bool_t hasITS=(track->GetNcls(0)>2);
158 printf("DUMPHPTTrack:%s|%f|%d|%d|%d\n",filename.Data(),track->Pt(), eventNumber,hasFriend,hasITS);
160 if (!esdFriend) continue;
161 AliESDfriendTrack *friendTrack = esdFriend->GetTrack(i);
162 if (!friendTrack) continue;
167 TObject *calibObject;
168 AliTPCseed *seed = 0;
169 for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
170 if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
174 fHPTTree = new TTree("HPT","HPT");
175 fHPTTree->SetDirectory(0);
177 if (fHPTTree->GetEntries()==0){
179 fHPTTree->SetDirectory(0);
180 fHPTTree->Branch("t.",&track);
181 fHPTTree->Branch("ft.",&friendTrack);
182 fHPTTree->Branch("s.",&seed);
184 fHPTTree->SetBranchAddress("t.",&track);
185 fHPTTree->SetBranchAddress("ft.",&friendTrack);
186 fHPTTree->SetBranchAddress("s.",&seed);
195 void AliTPCcalibV0::DumpToTree(AliESDEvent *esd){
197 // Dump V0s fith full firend information to the
199 Int_t nV0s = fESD->GetNumberOfV0s();
200 const Int_t kMinCluster=110;
201 const Double_t kDownscale=0.01;
202 const Float_t kMinPt =1.0;
203 const Float_t kMinMinPt =0.7;
204 AliESDfriend *esdFriend=static_cast<AliESDfriend*>(esd->FindListObject("AliESDfriend"));
207 for (Int_t ivertex=0; ivertex<nV0s; ivertex++){
209 AliESDv0 * v0 = (AliESDv0*) esd->GetV0(ivertex);
210 AliESDtrack * track0 = fESD->GetTrack(v0->GetIndex(0)); // negative track
211 AliESDtrack * track1 = fESD->GetTrack(v0->GetIndex(1)); // positive track
212 if (track0->GetTPCNcls()<kMinCluster) continue;
213 if (track0->GetKinkIndex(0)>0) continue;
214 if (track1->GetTPCNcls()<kMinCluster) continue;
215 if (track1->GetKinkIndex(0)>0) continue;
216 if (v0->GetOnFlyStatus()==kFALSE) continue;
218 if (TMath::Min(track0->Pt(),track1->Pt())<kMinMinPt) continue;
221 if (TMath::Max(track0->Pt(),track1->Pt())>kMinPt) isOK=kTRUE;
222 if (gRandom->Rndm()<kDownscale) isOK=kTRUE;
225 TString filename(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
226 Int_t eventNumber = esd->GetEventNumberInFile();
227 Bool_t hasITS=(track0->GetNcls(0)+ track1->GetNcls(0)>4);
228 printf("DUMPHPTV0:%s|%f|%d|%d|%d\n",filename.Data(), (TMath::Min(track0->Pt(),track1->Pt())), eventNumber,(esdFriend!=0), hasITS);
230 if (!esdFriend) continue;
234 AliESDfriendTrack *ftrack0 = esdFriend->GetTrack(v0->GetIndex(0));
235 if (!ftrack0) continue;
236 AliESDfriendTrack *ftrack1 = esdFriend->GetTrack(v0->GetIndex(1));
237 if (!ftrack1) continue;
239 TObject *calibObject;
240 AliTPCseed *seed0 = 0;
241 AliTPCseed *seed1 = 0;
242 for (Int_t l=0;(calibObject=ftrack0->GetCalibObject(l));++l) {
243 if ((seed0=dynamic_cast<AliTPCseed*>(calibObject))) break;
245 for (Int_t l=0;(calibObject=ftrack1->GetCalibObject(l));++l) {
246 if ((seed1=dynamic_cast<AliTPCseed*>(calibObject))) break;
248 if (!seed0) continue;
249 if (!seed1) continue;
250 AliExternalTrackParam * paramIn0 = (AliExternalTrackParam *)track0->GetInnerParam();
251 AliExternalTrackParam * paramIn1 = (AliExternalTrackParam *)track1->GetInnerParam();
252 if (!paramIn0) continue;
253 if (!paramIn1) continue;
257 fV0Tree = new TTree("V0s","V0s");
258 fV0Tree->SetDirectory(0);
260 if (fV0Tree->GetEntries()==0){
262 fV0Tree->SetDirectory(0);
263 fV0Tree->Branch("v0.",&v0);
264 fV0Tree->Branch("t0.",&track0);
265 fV0Tree->Branch("t1.",&track1);
266 fV0Tree->Branch("ft0.",&ftrack0);
267 fV0Tree->Branch("ft1.",&ftrack1);
268 fV0Tree->Branch("s0.",&seed0);
269 fV0Tree->Branch("s1.",&seed1);
271 fV0Tree->SetBranchAddress("v0.",&v0);
272 fV0Tree->SetBranchAddress("t0.",&track0);
273 fV0Tree->SetBranchAddress("t1.",&track1);
274 fV0Tree->SetBranchAddress("ft0.",&ftrack0);
275 fV0Tree->SetBranchAddress("ft1.",&ftrack1);
276 fV0Tree->SetBranchAddress("s0.",&seed0);
277 fV0Tree->SetBranchAddress("s1.",&seed1);
284 Long64_t AliTPCcalibV0::Merge(TCollection *const li) {
286 TIterator* iter = li->MakeIterator();
287 AliTPCcalibV0* cal = 0;
289 while ((cal = (AliTPCcalibV0*)iter->Next())) {
292 fV0Tree = new TTree("V0s","V0s");
293 fV0Tree->SetDirectory(0);
295 if (cal->fV0Tree->GetEntries()>0) AliTPCcalibV0::AddTree(cal->fV0Tree);
296 if (cal->fHPTTree->GetEntries()>0) AliTPCcalibV0::AddTreeHPT(cal->fHPTTree);
303 void AliTPCcalibV0::AddTree(TTree * treeInput){
305 // Add the content of tree:
306 // Notice automatic copy of tree in ROOT does not work for such complicated tree
309 AliESDv0 * v0 = new AliESDv0;
311 AliESDtrack * track0 = 0; // negative track
312 AliESDtrack * track1 = 0; // positive track
313 AliESDfriendTrack *ftrack0 = 0;
314 AliESDfriendTrack *ftrack1 = 0;
315 AliTPCseed *seed0 = 0;
316 AliTPCseed *seed1 = 0;
317 treeInput->SetBranchStatus("ft0.",kFALSE);
318 treeInput->SetBranchStatus("ft1.",kFALSE);
320 Double_t massK0= pdg.GetParticle("K0")->Mass();
321 Double_t massLambda= pdg.GetParticle("Lambda0")->Mass();
323 Int_t entries= treeInput->GetEntries();
324 for (Int_t i=0; i<entries; i++){
325 treeInput->SetBranchAddress("v0.",&v0);
326 treeInput->SetBranchAddress("t0.",&track0);
327 treeInput->SetBranchAddress("t1.",&track1);
328 treeInput->SetBranchAddress("ft0.",&ftrack0);
329 treeInput->SetBranchAddress("ft1.",&ftrack1);
330 treeInput->SetBranchAddress("s0.",&seed0);
331 treeInput->SetBranchAddress("s1.",&seed1);
332 if (fV0Tree->GetEntries()==0){
333 fV0Tree->SetDirectory(0);
334 fV0Tree->Branch("v0.",&v0);
335 fV0Tree->Branch("t0.",&track0);
336 fV0Tree->Branch("t1.",&track1);
337 fV0Tree->Branch("ft0.",&ftrack0);
338 fV0Tree->Branch("ft1.",&ftrack1);
339 fV0Tree->Branch("s0.",&seed0);
340 fV0Tree->Branch("s1.",&seed1);
342 fV0Tree->SetBranchAddress("v0.",&v0);
343 fV0Tree->SetBranchAddress("t0.",&track0);
344 fV0Tree->SetBranchAddress("t1.",&track1);
345 fV0Tree->SetBranchAddress("ft0.",&ftrack0);
346 fV0Tree->SetBranchAddress("ft1.",&ftrack1);
347 fV0Tree->SetBranchAddress("s0.",&seed0);
348 fV0Tree->SetBranchAddress("s1.",&seed1);
351 treeInput->GetEntry(i);
352 //ftrack0->GetCalibContainer()->SetOwner(kTRUE);
353 //ftrack1->GetCalibContainer()->SetOwner(kTRUE);
355 if (v0->GetOnFlyStatus()==kFALSE) isOK=kFALSE;
356 if (track0->GetTPCncls()<100) isOK=kFALSE;
357 if (track1->GetTPCncls()<100) isOK=kFALSE;
358 if (TMath::Min(seed0->Pt(),seed1->Pt())<kMinPt) isOK=kFALSE;
359 if (TMath::Min(track0->Pt(),track1->Pt())<kMinPt) isOK=kFALSE;
361 if (TMath::Abs(v0->GetEffMass(2,2)-massK0)<0.05) isV0=kTRUE;
362 if (TMath::Abs(v0->GetEffMass(4,2)-massLambda)<0.05) isV0=kTRUE;
363 if (TMath::Abs(v0->GetEffMass(2,4)-massLambda)<0.05) isV0=kTRUE;
364 if (TMath::Abs(v0->GetEffMass(0,0))<0.02) isV0=kFALSE; //reject electrons
365 if (!isV0) isOK=kFALSE;
366 if (isOK) fV0Tree->Fill();
384 void AliTPCcalibV0::AddTreeHPT(TTree * treeInput){
386 // Add the content of tree:
387 // Notice automatic copy of tree in ROOT does not work for such complicated tree
390 AliESDtrack *track = 0;
391 AliESDfriendTrack *friendTrack = 0;
392 AliTPCseed *seed = 0;
393 if (!treeInput) return;
394 if (treeInput->GetEntries()==0) return;
396 Int_t entries= treeInput->GetEntries();
398 for (Int_t i=0; i<entries; i++){
403 treeInput->SetBranchAddress("t.",&track);
404 treeInput->SetBranchAddress("ft.",&friendTrack);
405 treeInput->SetBranchAddress("s.",&seed);
406 treeInput->GetEntry(i);
408 if (fHPTTree->GetEntries()==0){
409 fHPTTree->SetDirectory(0);
410 fHPTTree->Branch("t.",&track);
411 fHPTTree->Branch("ft.",&friendTrack);
412 fHPTTree->Branch("s.",&seed);
414 fHPTTree->SetBranchAddress("t.",&track);
415 fHPTTree->SetBranchAddress("ft.",&friendTrack);
416 fHPTTree->SetBranchAddress("s.",&seed);
419 if (!track->IsOn(AliESDtrack::kITSrefit)) isOK=kFALSE;
420 if (!track->IsOn(AliESDtrack::kTOFout)) isOK=kFALSE;
421 if (isOK) fHPTTree->Fill();
430 void AliTPCcalibV0::MakeFitTreeTrack(const TObjArray * corrArray, Double_t ptCut, Int_t /*run*/){
434 // 0. Loop over selected tracks
435 // 1. Loop over all transformation - refit the track with and without the
437 // 2. Dump the matching paramaeters to the debugStremer
441 const Int_t kMinNcl=120;
442 TFile f("TPCV0Objects.root");
443 AliTPCcalibV0 *v0TPC = (AliTPCcalibV0*) f.Get("v0TPC");
444 TTree * treeInput = v0TPC->GetHPTTree();
445 TTreeSRedirector *pcstream = new TTreeSRedirector("fitHPT.root");
446 AliESDtrack *track = 0;
447 AliESDfriendTrack *friendTrack = 0;
448 AliTPCseed *seed = 0;
449 if (!treeInput) return;
450 if (treeInput->GetEntries()==0) return;
452 treeInput->SetBranchAddress("t.",&track);
453 treeInput->SetBranchAddress("ft.",&friendTrack);
454 treeInput->SetBranchAddress("s.",&seed);
457 if (corrArray) ncorr = corrArray->GetEntries();
458 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
459 // AliTPCParam *param = AliTPCcalibDB::Instance()->GetParameters();
460 // AliGRPObject* grp = AliTPCcalibDB::Instance()->GetGRP(run);
461 // Double_t time=0.5*(grp->GetTimeStart() +grp->GetTimeEnd());
465 Int_t ntracks= treeInput->GetEntries();
466 for (Int_t itrack=0; itrack<ntracks; itrack++){
467 treeInput->GetEntry(itrack);
468 if (!track) continue;
469 if (seed->Pt()<ptCut) continue;
470 if (track->Pt()<ptCut) continue;
471 if (track->GetTPCncls()<kMinNcl) continue;
473 // Reapply transformation
475 for (Int_t irow=0; irow<159; irow++){
476 AliTPCclusterMI *cluster=seed->GetClusterPointer(irow);
477 if (cluster &&cluster->GetX()>10){
478 Double_t x0[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
479 Int_t index0[1]={cluster->GetDetector()};
480 transform->Transform(x0,index0,0,1);
481 cluster->SetX(x0[0]);
482 cluster->SetY(x0[1]);
483 cluster->SetZ(x0[2]);
488 AliExternalTrackParam* paramInner=0;
489 AliExternalTrackParam* paramOuter=0;
490 AliExternalTrackParam* paramIO=0;
492 for (Int_t icorr=-1; icorr<ncorr; icorr++){
494 AliTPCCorrection *corr = 0;
495 if (icorr>=0) corr = (AliTPCCorrection*)corrArray->At(icorr);
496 AliExternalTrackParam * trackInner = RefitTrack(seed, corr,85,134,0.1);
497 AliExternalTrackParam * trackIO = RefitTrack(seed, corr,245,85,0.1);
498 AliExternalTrackParam * trackOuter = RefitTrack(seed, corr,245,134,0.1 );
499 if (trackInner&&trackOuter&&trackIO){
500 trackOuter->Rotate(trackInner->GetAlpha());
501 trackOuter->PropagateTo(trackInner->GetX(),AliTracker::GetBz());
503 paramInner=trackInner;
504 paramOuter=trackOuter;
506 paramIO->Rotate(seed->GetAlpha());
507 paramIO->PropagateTo(seed->GetX(),AliTracker::GetBz());
514 if (paramOuter&& paramInner) {
515 // Bool_t isOK=kTRUE;
516 if (paramInner->GetSigmaY2()>0.01) isOK&=kFALSE;
517 if (paramOuter->GetSigmaY2()>0.01) isOK&=kFALSE;
518 if (paramInner->GetSigmaZ2()>0.01) isOK&=kFALSE;
519 if (paramOuter->GetSigmaZ2()>0.01) isOK&=kFALSE;
523 "pIn.="<<paramInner<<
524 "pOut.="<<paramOuter;
532 .x ../ConfigCalibTrain.C(run)
533 .L ../AddTaskTPCCalib.C
535 TFile fexb("../../RegisterCorrectionExB.root");
536 AliTPCComposedCorrection *cc= (AliTPCComposedCorrection*) fexb.Get("ComposedExB");
538 cc->Print("DA"); // Print used correction classes
539 TObjArray *array = cc->GetCorrections()
540 AliTPCcalibV0::MakeFitTreeTrack(array,2,run);
545 void AliTPCcalibV0::MakeFitTreeV0(const TObjArray * corrArray, Double_t ptCut, Int_t run){
549 // 0. Loop over selected tracks
550 // 1. Loop over all transformation - refit the track with and without the
552 // 2. Dump the matching paramaeters to the debugStremer
556 TFile f("TPCV0Objects.root");
557 AliTPCcalibV0 *v0TPC = (AliTPCcalibV0*) f.Get("v0TPC");
558 TTree * treeInput = v0TPC->GetV0Tree();
559 TTreeSRedirector *pcstream = new TTreeSRedirector("fitV0.root");
561 AliESDtrack *track0 = 0;
562 AliESDfriendTrack *friendTrack0 = 0;
563 AliTPCseed *seed0 = 0;
565 AliESDtrack *track1 = 0;
566 AliESDfriendTrack *friendTrack1 = 0;
568 AliTPCseed *seed1 = 0;
569 if (!treeInput) return;
570 if (treeInput->GetEntries()==0) return;
572 treeInput->SetBranchAddress("v0.",&v0);
573 treeInput->SetBranchAddress("t0.",&track0);
574 treeInput->SetBranchAddress("ft0.",&friendTrack0);
575 treeInput->SetBranchAddress("s0.",&s0);
576 treeInput->SetBranchAddress("t1.",&track1);
577 treeInput->SetBranchAddress("ft1.",&friendTrack1);
578 treeInput->SetBranchAddress("s1.",&s1);
582 if (corrArray) ncorr = corrArray->GetEntries();
583 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
584 Double_t massK0= pdg.GetParticle("K0")->Mass();
585 Double_t massLambda= pdg.GetParticle("Lambda0")->Mass();
586 Double_t massPion=pdg.GetParticle("pi+")->Mass();
587 Double_t massProton=pdg.GetParticle("proton")->Mass();
588 Int_t pdgPion=pdg.GetParticle("pi+")->PdgCode();
589 Int_t pdgProton=pdg.GetParticle("proton")->PdgCode();
598 Int_t nv0s= treeInput->GetEntries();
599 for (Int_t iv0=0; iv0<nv0s; iv0++){
603 Int_t isAntiLambda=0;
604 treeInput->GetEntry(iv0);
605 if (TMath::Abs(v0->GetEffMass(2,2)-massK0)<0.03) {isK0=1; v0Type=1;} //select K0s
606 if (TMath::Abs(v0->GetEffMass(4,2)-massLambda)<0.01) {isLambda=1; v0Type=2;} //select Lambda
607 if (TMath::Abs(v0->GetEffMass(2,4)-massLambda)<0.01) {isAntiLambda=1;v0Type=3;} //select Anti Lambda
608 if (isK0+isLambda+isAntiLambda!=1) continue;
613 if (isLambda) {rmass0=massProton; pdg0=pdgProton;}
614 if (isAntiLambda) {rmass1=massProton; pdg1=pdgProton;}
616 if (isK0==0) massV0=massLambda;
619 seed0=(s0->GetSign()>0)?s0:s1;
620 seed1=(s0->GetSign()>0)?s1:s0;
621 if (seed0->GetZ()*seed1->GetZ()<0) continue; //remove membrane crossed tracks
622 if (seed0->Pt()<ptCut) continue;
623 if (seed1->Pt()<ptCut) continue;
625 // Reapply transformation
627 for (Int_t itype=0; itype<2; itype++){
628 AliTPCseed * seed = (itype==0) ? seed0: seed1;
629 for (Int_t irow=0; irow<159; irow++){
630 AliTPCclusterMI *cluster=seed->GetClusterPointer(irow);
631 if (cluster &&cluster->GetX()>10){
632 Double_t x0[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
633 Int_t index0[1]={cluster->GetDetector()};
634 transform->Transform(x0,index0,0,1);
635 cluster->SetX(x0[0]);
636 cluster->SetY(x0[1]);
637 cluster->SetZ(x0[2]);
643 Double_t radius = v0->GetRr();
645 v0->GetXYZ(xyz[0],xyz[1],xyz[2]);
646 Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
647 TObjArray arrayV0in(ncorr+1);
648 TObjArray arrayV0io(ncorr+1);
649 TObjArray arrayT0(ncorr+1);
650 TObjArray arrayT1(ncorr+1);
651 arrayV0in.SetOwner(kTRUE);
652 arrayV0io.SetOwner(kTRUE);
654 for (Int_t icorr=-1; icorr<ncorr; icorr++){
655 AliTPCCorrection *corr =0;
656 if (icorr>=0) corr = (AliTPCCorrection*)corrArray->At(icorr);
658 AliExternalTrackParam * trackInner0 = RefitTrack(seed0, corr,160,85,rmass0);
659 AliExternalTrackParam * trackIO0 = RefitTrack(seed0, corr,245,85,rmass0);
660 AliExternalTrackParam * trackInner1 = RefitTrack(seed1, corr,160,85,rmass1);
661 AliExternalTrackParam * trackIO1 = RefitTrack(seed1, corr,245,85,rmass1);
662 if (!trackInner0) isOK=kFALSE;
663 if (!trackInner1) isOK=kFALSE;
664 if (!trackIO0) isOK=kFALSE;
665 if (!trackIO1) isOK=kFALSE;
667 if (!trackInner0->Rotate(alpha)) isOK=kFALSE;
668 if (!trackInner1->Rotate(alpha)) isOK=kFALSE;
669 if (!trackIO0->Rotate(alpha)) isOK=kFALSE;
670 if (!trackIO1->Rotate(alpha)) isOK=kFALSE;
672 if (!AliTracker::PropagateTrackToBxByBz(trackInner0, radius, rmass0, 1, kFALSE)) isOK=kFALSE;
673 if (!AliTracker::PropagateTrackToBxByBz(trackInner1, radius, rmass1, 1, kFALSE)) isOK=kFALSE;
674 if (!AliTracker::PropagateTrackToBxByBz(trackIO0, radius, rmass0, 1, kFALSE)) isOK=kFALSE;
675 if (!AliTracker::PropagateTrackToBxByBz(trackIO1, radius, rmass1, 1, kFALSE)) isOK=kFALSE;
677 arrayT0.AddAt(trackIO0->Clone(),icorr+1);
678 arrayT1.AddAt(trackIO1->Clone(),icorr+1);
679 Int_t charge=TMath::Nint(trackIO0->GetSign());
680 AliKFParticle pin0( *trackInner0, pdg0*charge);
681 AliKFParticle pin1( *trackInner1, -pdg1*charge);
682 AliKFParticle pio0( *trackIO0, pdg0*charge);
683 AliKFParticle pio1( *trackIO1, -pdg1*charge);
690 arrayV0in.AddAt(v0in.Clone(),icorr+1);
691 arrayV0io.AddAt(v0io.Clone(),icorr+1);
696 //AliKFParticle* pin0= (AliKFParticle*)arrayV0in.At(0);
697 AliKFParticle* pio0= (AliKFParticle*)arrayV0io.At(0);
698 AliExternalTrackParam *param0=(AliExternalTrackParam *)arrayT0.At(0);
699 AliExternalTrackParam *param1=(AliExternalTrackParam *)arrayT1.At(0);
700 Double_t mass0=0, mass0E=0;
701 pio0->GetMass( mass0,mass0E);
703 Double_t mean=mass0-massV0;
704 if (TMath::Abs(mean)>0.05) continue;
705 Double_t mass[10000];
707 Int_t dtype=30; // id for V0
708 Int_t ptype=5; // id for invariant mass
709 // Int_t id=TMath::Nint(100.*(param0->Pt()-param1->Pt())/(param0->Pt()+param1->Pt())); // K0s V0 asymetry
710 Int_t id=Int_t(1000.*(param0->Pt()-param1->Pt())); // K0s V0 asymetry
711 Double_t gx,gy,gz, px,py,pz;
712 Double_t pt = v0->Pt();
713 v0->GetXYZ(gx,gy,gz);
714 v0->GetPxPyPz(px,py,pz);
715 Double_t theta=pz/TMath::Sqrt(px*px+py*py);
716 Double_t phi=TMath::ATan2(py,px);
717 Double_t snp=0.5*(seed0->GetSnp()+seed1->GetSnp());
718 Double_t sector=9*phi/TMath::Pi();
719 Double_t dsec=sector-TMath::Nint(sector);
720 Double_t refX=TMath::Sqrt(gx*gx+gy*gy);
721 //Int_t nentries=v0Type;
722 Double_t bz=AliTracker::GetBz();
724 (*pcstream)<<"fitDebug"<<
732 (*pcstream)<<"fit"<< // dump valus for fit
733 "run="<<run<< //run number
734 "bz="<<bz<< // magnetic filed used
735 "dtype="<<dtype<< // detector match type 30
736 "ptype="<<ptype<< // parameter type
737 "theta="<<theta<< // theta
740 "mean="<<mean<< // mean dist value
741 "rms="<<mass0E<< // rms
745 "refX="<<refX<< // reference radius
746 "gx="<<gx<< // global position
747 "gy="<<gy<< // global position
748 "gz="<<gz<< // global position
749 "dRrec="<<dRrec<< // delta Radius in reconstruction
750 "pt="<<pt<< // pt of the particle
751 "id="<<id<< //delta of the momenta
752 "entries="<<v0Type;// type of the V0
753 for (Int_t icorr=0; icorr<ncorr; icorr++){
754 AliTPCCorrection *corr =0;
755 if (icorr>=0) corr = (AliTPCCorrection*)corrArray->At(icorr);
756 // AliKFParticle* pin= (AliKFParticle*)arrayV0in.At(icorr+1);
757 AliKFParticle* pio= (AliKFParticle*)arrayV0io.At(icorr+1);
758 AliExternalTrackParam *par0=(AliExternalTrackParam *)arrayT0.At(icorr+1);
759 AliExternalTrackParam *par1=(AliExternalTrackParam *)arrayT1.At(icorr+1);
761 pio->GetMass( mass[icorr],massE);
764 Form("%s=",corr->GetName())<<mass[icorr];
765 (*pcstream)<<"fitDebug"<<
766 Form("%s=",corr->GetName())<<mass[icorr]<<
767 Form("%sp0.=",corr->GetName())<<par0<<
768 Form("%sp1=",corr->GetName())<<par1;
770 (*pcstream)<<"fit"<< "isOK="<<isOK<<"\n";
771 (*pcstream)<<"fitDebug"<< "isOK="<<isOK<<"\n";
776 AliExternalTrackParam * AliTPCcalibV0::RefitTrack(AliTPCseed *seed, AliTPCCorrection * corr, Double_t xstart, Double_t xstop, Double_t mass){
779 // seed - tpc track with cluster
780 // corr - distrotion/correction class - apllied to the points
781 // xstart - radius to start propagate/update
782 // xstop - radius to stop propagate/update
784 const Double_t kResetCov=20.;
785 const Double_t kSigma=5.;
787 for (Int_t i=0;i<15;i++) covar[i]=0;
788 covar[0]=kSigma*kSigma;
789 covar[2]=kSigma*kSigma;
790 covar[5]=kSigma*kSigma/Float_t(150*150);
791 covar[9]=kSigma*kSigma/Float_t(150*150);
794 AliExternalTrackParam *refit = new AliExternalTrackParam(*seed);
795 refit->PropagateTo(xstart, AliTracker::GetBz());
796 refit->AddCovariance(covar);
797 refit->ResetCovariance(kResetCov);
798 Double_t xmin = TMath::Min(xstart,xstop);
799 Double_t xmax = TMath::Max(xstart,xstop);
803 for (Int_t index0=0; index0<159; index0++){
804 Int_t irow= (xstart<xstop)? index0:159-index0;
805 AliTPCclusterMI *cluster=seed->GetClusterPointer(irow); //cluster in local system
806 if (!cluster) continue;
807 if (cluster->GetX()<xmin) continue;
808 if (cluster->GetX()>xmax) continue;
809 Double_t alpha = TMath::Pi()*(cluster->GetDetector()+0.5)/9.;
810 if (!refit->Rotate(alpha)) isOK=kFALSE;
811 Double_t x = cluster->GetX();
812 Double_t y = cluster->GetY();
813 Double_t z = cluster->GetZ();
815 Float_t xyz[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()}; // original position
816 corr->DistortPointLocal(xyz,cluster->GetDetector());
821 if (!AliTracker::PropagateTrackToBxByBz(refit, x,mass,1.,kFALSE)) isOK=kFALSE;
823 Double_t cov[3]={0.01,0.,0.01};
824 Double_t yz[2]={y,z};
825 if (!refit->Update(yz,cov)) isOK=kFALSE;
828 if (!AliTracker::PropagateTrackToBxByBz(refit, xstop, mass,1.,kTRUE)) isOK=kFALSE;
848 AliKFParticle * AliTPCcalibV0::Fit(AliKFVertex & /*primVtx*/, AliESDv0 *v0, Int_t PDG1, Int_t PDG2){
852 AliKFParticle p1( *(v0->GetParamN()), PDG1 );
853 AliKFParticle p2( *(v0->GetParamP()), PDG2 );
854 AliKFParticle *V0 = new AliKFParticle;
857 V0->SetVtxGuess(x,y,z);
866 void AliTPCcalibV0::BinLogX(TH2F *h) {
870 TAxis *axis = h->GetXaxis();
871 int bins = axis->GetNbins();
873 Double_t from = axis->GetXmin();
874 Double_t to = axis->GetXmax();
875 Double_t *new_bins = new Double_t[bins + 1];
878 Double_t factor = pow(to/from, 1./bins);
880 for (int i = 1; i <= bins; i++) {
881 new_bins[i] = factor * new_bins[i-1];
883 axis->Set(bins, new_bins);
889 void AliTPCcalibV0::FilterV0s(AliESDEvent* event){
893 const Double_t kChi2Cut=20;
894 const Double_t kMinR=2;
895 const Double_t ptCut=0.2;
896 const Int_t kMinNcl=110;
898 Int_t nv0 = event->GetNumberOfV0s();
899 AliESDVertex *vertex= (AliESDVertex *)event->GetPrimaryVertex();
900 AliKFVertex kfvertex=*vertex;
902 for (Int_t iv0=0;iv0<nv0;iv0++){
903 AliESDv0 *v0 = event->GetV0(iv0);
905 if (v0->GetPindex()<0) continue;
906 if (v0->GetNindex()<0) continue;
907 if (TMath::Max(v0->GetPindex(), v0->GetNindex())>event->GetNumberOfTracks()) continue;
910 AliExternalTrackParam pp=(v0->GetParamP()->GetSign()>0) ? (*(v0->GetParamP())):(*(v0->GetParamN()));
911 AliExternalTrackParam pn=(v0->GetParamP()->GetSign()>0) ? (*(v0->GetParamN())):(*(v0->GetParamP()));
912 AliKFParticle kfp1( pp, 211 );
913 AliKFParticle kfp2( pn, -211 );
914 AliKFParticle *v0KFK0 = new AliKFParticle(kfp1,kfp2);
915 AliKFParticle *v0KFK0CV = new AliKFParticle(*v0KFK0);
916 v0KFK0CV->SetProductionVertex(kfvertex);
917 v0KFK0CV->TransportToProductionVertex();
918 AliKFParticle *v0KFK0CVM = new AliKFParticle(*v0KFK0CV);
919 v0KFK0CVM->SetMassConstraint(pdg.GetParticle("K_S0")->Mass());
920 Double_t chi2K0 = v0KFK0CV->GetChi2();
921 Double_t chi2K0M= v0KFK0CVM->GetChi2();
922 Double_t massK0=v0KFK0CV->GetMass();
923 if (chi2K0>kChi2Cut) continue;
924 if (v0->GetRr()>kMinR) continue;
926 Double_t maxPt = TMath::Max(v0->GetParamP()->Pt(), v0->GetParamN()->Pt());
927 Double_t effMass22=v0->GetEffMass(2,2);
928 Double_t effMass42=v0->GetEffMass(4,2);
929 Double_t effMass24=v0->GetEffMass(2,4);
931 isV0|=TMath::Abs(effMass22-pdg.GetParticle("K_S0")->Mass())<0.1;
932 isV0|=TMath::Abs(effMass42-pdg.GetParticle("Lambda0")->Mass())<0.1;
933 isV0|=TMath::Abs(effMass24-pdg.GetParticle("Lambda0")->Mass())<0.1;
935 Double_t sign= v0->GetParamP()->GetSign()* v0->GetParamN()->GetSign();
936 if (sign<0&&v0->GetOnFlyStatus()>0.5&&maxPt>ptCut&&isV0){
937 AliESDtrack * trackP = event->GetTrack(v0->GetPindex());
938 AliESDtrack * trackN = event->GetTrack(v0->GetNindex());
939 if (!trackN) continue;
940 if (!trackP) continue;
941 Int_t nclP= (Int_t)trackP->GetTPCClusterInfo(2,1);
942 Int_t nclN= (Int_t)trackN->GetTPCClusterInfo(2,1);
943 if (TMath::Min(nclP,nclN)<kMinNcl) continue;
944 Double_t eta = TMath::Max(TMath::Abs(trackP->Eta()), TMath::Abs(trackN->Eta()));
945 Double_t ncls = TMath::Min(TMath::Abs(trackP->GetNcls(0)), TMath::Abs(trackN->GetNcls(0)));
946 if (eta<0.8&&ncls>2){
947 // printf("%d\t%f\t%f\t%d\t%d\t%f\t%f\n",i, v0->Pt(), maxPt, v0->GetNindex(),v0->GetPindex(),v0->GetRr(),effMass22);
948 (*fDebugStreamer)<<"v0tree"<<