X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;ds=sidebyside;f=TPC%2FAliTPCcalibV0.cxx;h=66f220d78206936182866ef3137daf58280abaf6;hb=063feb044e455efa978d8487f6207ae2cd09838f;hp=96c0ca4f9c5a16a98cc3448c6c428357948e99a5;hpb=4ce766eb4226c070c454ba72807fbf481119a5cb;p=u%2Fmrichter%2FAliRoot.git diff --git a/TPC/AliTPCcalibV0.cxx b/TPC/AliTPCcalibV0.cxx index 96c0ca4f9c5..66f220d7820 100644 --- a/TPC/AliTPCcalibV0.cxx +++ b/TPC/AliTPCcalibV0.cxx @@ -38,8 +38,13 @@ #include "TCint.h" #include "AliTPCcalibV0.h" #include "AliV0.h" - - +#include "TRandom.h" +#include "TTreeStream.h" +#include "AliTPCcalibDB.h" +#include "AliTPCCorrection.h" +#include "AliGRPObject.h" +#include "AliTPCTransform.h" +#include "AliAnalysisManager.h" @@ -48,693 +53,798 @@ ClassImp(AliTPCcalibV0) AliTPCcalibV0::AliTPCcalibV0() : AliTPCcalibBase(), + fV0Tree(0), + fHPTTree(0), fStack(0), fESD(0), fPdg(0), fParticles(0), fV0s(0), - fGammas(0), - fV0type(0), - fTPCdEdx(0), - fTPCdEdxPi(0), - fTPCdEdxEl(0), - fTPCdEdxP(0) + fGammas(0) +{ + +} +AliTPCcalibV0::AliTPCcalibV0(const Text_t *name, const Text_t *title): + AliTPCcalibBase(), + fV0Tree(0), + fHPTTree(0), + fStack(0), + fESD(0), + fPdg(0), + fParticles(0), + fV0s(0), + fGammas(0) { fPdg = new TDatabasePDG; // create output histograms - fTPCdEdx = new TH2F("TPCdEdX", "dE/dX; BetaGamma; TPC signal (a.u.)", 1000, 0.1, 10000, 300, 0, 300); - BinLogX(fTPCdEdx); - fTPCdEdxPi = new TH2F("TPCdEdXPi","dE/dX; BetaGamma; TPC signal (a.u.)", 1000, 0.1, 10000, 300, 0, 300); - fTPCdEdxEl = new TH2F("TPCdEdXEl","dE/dX; BetaGamma; TPC signal (a.u.)", 1000, 0.1, 10000, 300, 0, 300); - fTPCdEdxP = new TH2F("TPCdEdXP", "dE/dX; BetaGamma; TPC signal (a.u.)", 1000, 0.1, 10000, 300, 0, 300); - - + SetName(name); + SetTitle(title); } AliTPCcalibV0::~AliTPCcalibV0(){ // // // + delete fV0Tree; + delete fHPTTree; } -void AliTPCcalibV0::ProcessESD(AliESDEvent *esd, AliStack *stack){ +void AliTPCcalibV0::ProcessESD(AliESDEvent *esd){ // // // fESD = esd; AliKFParticle::SetField(esd->GetMagneticField()); - MakeV0s(); - if (stack) { - fStack = stack; - MakeMC(); - }else{ - fStack =0; - } + if (TMath::Abs(AliTracker::GetBz())<1) return; + DumpToTree(esd); + DumpToTreeHPT(esd); } -void AliTPCcalibV0::MakeMC(){ - // - // MC comparison - // 1. Select interesting particles - // 2. Assign the recosntructed particles +void AliTPCcalibV0::DumpToTreeHPT(AliESDEvent *esd){ // - //1. Select interesting particles - const Float_t kMinP = 0.2; - const Float_t kMinPt = 0.1; - const Float_t kMaxR = 0.5; - const Float_t kMaxTan = 1.2; - const Float_t kMaxRad = 150; + // Dump V0s fith full firend information to the + // + if (TMath::Abs(AliTracker::GetBz())<1) return; + const Int_t kMinCluster=110; + const Float_t kMinPt =4.; + AliESDfriend *esdFriend=static_cast(esd->FindListObject("AliESDfriend")); +// if (!esdFriend) { +// Printf("ERROR: esdFriend not available"); +// return; +// } // - if (!fParticles) fParticles = new TObjArray; - TParticle *part=0; - // - Int_t entries = fStack->GetNtrack(); - for (Int_t ipart=0; ipartParticle(ipart); - if (!part) continue; - if (part->P()R()>kMaxR) continue; - if (TMath::Abs(TMath::Tan(part->Theta()-TMath::Pi()*0.5))>kMaxTan) continue; - Bool_t isInteresting = kFALSE; - if (part->GetPdgCode()==22) isInteresting =kTRUE; - if (part->GetPdgCode()==310) isInteresting =kTRUE; - if (part->GetPdgCode()==111) isInteresting =kTRUE; - if (TMath::Abs(part->GetPdgCode()==3122)) isInteresting =kTRUE; + Int_t ntracks=esd->GetNumberOfTracks(); + for (Int_t i=0;iGetTrack(i); + if (track->GetTPCncls()1){ // cut on momenta if measured + if (track->Pt()>kMinPt) isOK=kTRUE; + } + if (TMath::Abs(AliTracker::GetBz())<1){ // require primary track for the B field OFF data + Bool_t isAccepted=kTRUE; + if (!track->IsOn(AliESDtrack::kITSrefit)) isAccepted=kFALSE; + if (!track->IsOn(AliESDtrack::kTPCrefit)) isAccepted=kFALSE; + if (!track->IsOn(AliESDtrack::kTOFout)) isAccepted=kFALSE; + Float_t dvertex[2],cvertex[3]; + track->GetImpactParametersTPC(dvertex,cvertex); + if (TMath::Abs(dvertex[0]/TMath::Sqrt(cvertex[0]+0.01))>20) isAccepted=kFALSE; + if (TMath::Abs(dvertex[1]/TMath::Sqrt(TMath::Abs(cvertex[2]+0.01)))>20) isAccepted=kFALSE; + track->GetImpactParameters(dvertex,cvertex); + if (TMath::Abs(dvertex[0]/TMath::Sqrt(cvertex[0]+0.01))>10) isAccepted=kFALSE; + if (TMath::Abs(dvertex[1]/TMath::Sqrt(TMath::Abs(cvertex[2]+0.01)))>10) isAccepted=kFALSE; + if (!isAccepted) isOK=kFALSE; + } + if ( track->GetTPCsignal()>100 && track->GetInnerParam()->Pt()>1 ){ + if (track->IsOn(AliESDtrack::kITSin)||track->IsOn(AliESDtrack::kTRDout)||track->IsOn(AliESDtrack::kTOFin)) + isOK=kTRUE; + if (isOK){ + TString filename(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName()); + Int_t eventNumber = esd->GetEventNumberInFile(); + Bool_t hasFriend=(esdFriend) ? (esdFriend->GetTrack(i)!=0):0; + Bool_t hasITS=(track->GetNcls(0)>2); + printf("DUMPIONTrack:%s|%f|%d|%d|%d\n",filename.Data(),track->GetInnerParam()->Pt()*track->GetTPCsignal()/50., eventNumber,hasFriend,hasITS); + } + } + if (!isOK) continue; + TString filename(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName()); + Int_t eventNumber = esd->GetEventNumberInFile(); + Bool_t hasFriend=(esdFriend) ? (esdFriend->GetTrack(i)!=0):0; + Bool_t hasITS=(track->GetNcls(0)>2); + printf("DUMPHPTTrack:%s|%f|%d|%d|%d\n",filename.Data(),track->Pt(), eventNumber,hasFriend,hasITS); + // + if (!esdFriend) continue; + AliESDfriendTrack *friendTrack = esdFriend->GetTrack(i); + if (!friendTrack) continue; + if (!isOK) continue; // - if (!isInteresting) continue; - fParticles->AddLast(new TParticle(*part)); - } - if (fParticles->GetEntries()<1) { - return; - } - // - // - // - Int_t sentries=fParticles->GetEntries();; - for (Int_t ipart=0; ipartAt(ipart); - TParticle *p0 = 0; - TParticle *p1 = 0; - - Int_t nold =0; - Int_t nnew =0; - Int_t id0 = part->GetDaughter(0); - Int_t id1 = part->GetDaughter(1); - if (id0>=fStack->GetNtrack() ) id0*=-1; - if (id1>=fStack->GetNtrack() ) id1*=-1; - Bool_t findable = kTRUE; - if (id0<0 || id1<0) findable = kFALSE; - Int_t charge =0; - if (findable){ - p0 = fStack->Particle(id0); - if (p0->R()>kMaxRad) findable = kFALSE; - if (p0->Pt()Vz()>250) findable= kFALSE; - if (TMath::Abs(TMath::Tan(p0->Theta()-TMath::Pi()*0.5))>2) findable=kFALSE; - if (fPdg->GetParticle(p0->GetPdgCode())==0) findable =kFALSE; - else - if (fPdg->GetParticle(p0->GetPdgCode())->Charge()==0) charge++; - - p1 = fStack->Particle(id1); - if (p1->R()>kMaxRad) findable = kFALSE; - if (p1->Pt()Vz())>250) findable= kFALSE; - if (TMath::Abs(TMath::Tan(p1->Theta()-TMath::Pi()*0.5))>2) findable=kFALSE; - if (fPdg->GetParticle(p1->GetPdgCode())==0) findable = kFALSE; - else - if (fPdg->GetParticle(p1->GetPdgCode())->Charge()==0) charge++; - + + TObject *calibObject; + AliTPCseed *seed = 0; + for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) { + if ((seed=dynamic_cast(calibObject))) break; } - // (*fDebugStream)<<"MC0"<< - // "P.="<Pt(), p1->Pt()); - Int_t type=-1; - - // - // - AliKFVertex primVtx(*(fESD->GetPrimaryVertex())); - for (Int_t ivertex=0; ivertexGetNumberOfV0s(); ivertex++){ - AliESDv0 * v0 = fESD->GetV0(ivertex); - // select coresponding track - AliESDtrack * trackN = fESD->GetTrack(v0->GetIndex(0)); - if (TMath::Abs(trackN->GetLabel())!=id0 && TMath::Abs(trackN->GetLabel())!=id1) continue; - AliESDtrack * trackP = fESD->GetTrack(v0->GetIndex(1)); - if (TMath::Abs(trackP->GetLabel())!=id0 && TMath::Abs(trackP->GetLabel())!=id1) continue; - TParticle *pn = fStack->Particle(TMath::Abs(trackN->GetLabel())); - TParticle *pp = fStack->Particle(TMath::Abs(trackP->GetLabel())); - // - // - if ( v0->GetOnFlyStatus()) nnew++; - if (!v0->GetOnFlyStatus()) nold++; - if (part->GetPdgCode()==22 && TMath::Abs(pn->GetPdgCode())==11 && TMath::Abs(pp->GetPdgCode())==11) - type =1; - if (part->GetPdgCode()==310 && TMath::Abs(pn->GetPdgCode())==211 && TMath::Abs(pp->GetPdgCode())==211) - type =0; - if (part->GetPdgCode()==3122){ - if (TMath::Abs(pn->GetPdgCode())==210 ) type=2; - else type=3; - } - AliKFParticle *v0kf = Fit(primVtx,v0,pn->GetPdgCode(),pp->GetPdgCode()); - v0kf->SetProductionVertex( primVtx ); - AliKFParticle *v0kfc = Fit(primVtx,v0,pn->GetPdgCode(),pp->GetPdgCode()); - v0kfc->SetProductionVertex( primVtx ); - v0kfc->SetMassConstraint(fPdg->GetParticle(part->GetPdgCode())->Mass()); - Float_t chi2 = v0kf->GetChi2(); - Float_t chi2C = v0kf->GetChi2(); - // + if (!seed) continue; + if (!fHPTTree) { + fHPTTree = new TTree("HPT","HPT"); + fHPTTree->SetDirectory(0); + } + if (fHPTTree->GetEntries()==0){ // - TTreeSRedirector *cstream = GetDebugStreamer(); - if (cstream){ - (*cstream)<<"MCRC"<< - "P.="<SetDirectory(0); + fHPTTree->Branch("t.",&track); + fHPTTree->Branch("ft.",&friendTrack); + fHPTTree->Branch("s.",&seed); + }else{ + fHPTTree->SetBranchAddress("t.",&track); + fHPTTree->SetBranchAddress("ft.",&friendTrack); + fHPTTree->SetBranchAddress("s.",&seed); } - fParticles->Delete(); + fHPTTree->Fill(); + // } } -void AliTPCcalibV0::MakeV0s(){ - // - // +void AliTPCcalibV0::DumpToTree(AliESDEvent *esd){ // - const Int_t kMinCluster=40; - const Float_t kMinR =0; - if (! fV0s) fV0s = new TObjArray(10); - fV0s->Clear(); - // - // Old V0 finder - // - for (Int_t ivertex=0; ivertexGetNumberOfV0s(); ivertex++){ - AliESDv0 * v0 = fESD->GetV0(ivertex); - if (v0->GetOnFlyStatus()) continue; - fV0s->AddLast(v0); - } - ProcessV0(0); - fV0s->Clear(0); - // - // MI V0 finder - // - for (Int_t ivertex=0; ivertexGetNumberOfV0s(); ivertex++){ - AliESDv0 * v0 = fESD->GetV0(ivertex); - if (!v0->GetOnFlyStatus()) continue; - fV0s->AddLast(v0); - } - ProcessV0(1); - fV0s->Clear(); - return; - // - // combinatorial + // Dump V0s fith full firend information to the + // + Int_t nV0s = fESD->GetNumberOfV0s(); + const Int_t kMinCluster=110; + const Double_t kDownscale=0.01; + const Float_t kMinPt =1.0; + const Float_t kMinMinPt =0.7; + AliESDfriend *esdFriend=static_cast(esd->FindListObject("AliESDfriend")); // - Int_t ntracks = fESD->GetNumberOfTracks(); - for (Int_t itrack0=0; itrack0GetTrack(itrack0); - if (track0->GetSign()>0) continue; - if ( track0->GetTPCNcls()GetV0(ivertex); + AliESDtrack * track0 = fESD->GetTrack(v0->GetIndex(0)); // negative track + AliESDtrack * track1 = fESD->GetTrack(v0->GetIndex(1)); // positive track + if (track0->GetTPCNcls()GetKinkIndex(0)>0) continue; + if (track1->GetTPCNcls()GetKinkIndex(0)>0) continue; + if (v0->GetOnFlyStatus()==kFALSE) continue; + // + if (TMath::Min(track0->Pt(),track1->Pt())Pt(),track1->Pt())>kMinPt) isOK=kTRUE; + if (gRandom->Rndm()GetTree()->GetCurrentFile()->GetName()); + Int_t eventNumber = esd->GetEventNumberInFile(); + Bool_t hasITS=(track0->GetNcls(0)+ track1->GetNcls(0)>4); + printf("DUMPHPTV0:%s|%f|%d|%d|%d\n",filename.Data(), (TMath::Min(track0->Pt(),track1->Pt())), eventNumber,(esdFriend!=0), hasITS); + // + if (!esdFriend) continue; + // + + // + AliESDfriendTrack *ftrack0 = esdFriend->GetTrack(v0->GetIndex(0)); + if (!ftrack0) continue; + AliESDfriendTrack *ftrack1 = esdFriend->GetTrack(v0->GetIndex(1)); + if (!ftrack1) continue; // - for (Int_t itrack1=0; itrack1GetTrack(itrack1); - if (track1->GetSign()<0) continue; - if ( track1->GetTPCNcls()GetKinkIndex(0)>0) continue; + TObject *calibObject; + AliTPCseed *seed0 = 0; + AliTPCseed *seed1 = 0; + for (Int_t l=0;(calibObject=ftrack0->GetCalibObject(l));++l) { + if ((seed0=dynamic_cast(calibObject))) break; + } + for (Int_t l=0;(calibObject=ftrack1->GetCalibObject(l));++l) { + if ((seed1=dynamic_cast(calibObject))) break; + } + if (!seed0) continue; + if (!seed1) continue; + AliExternalTrackParam * paramIn0 = (AliExternalTrackParam *)track0->GetInnerParam(); + AliExternalTrackParam * paramIn1 = (AliExternalTrackParam *)track1->GetInnerParam(); + if (!paramIn0) continue; + if (!paramIn1) continue; + // + // + if (!fV0Tree) { + fV0Tree = new TTree("V0s","V0s"); + fV0Tree->SetDirectory(0); + } + if (fV0Tree->GetEntries()==0){ // - // AliExternalTrackParam param0(*track0); - // AliExternalTrackParam param1(*track1); - AliV0 vertex; - vertex.SetParamN(*track0); - vertex.SetParamP(*track1); - Float_t xyz[3]; - xyz[0] = fESD->GetPrimaryVertex()->GetXv(); - xyz[1] = fESD->GetPrimaryVertex()->GetYv(); - xyz[2] = fESD->GetPrimaryVertex()->GetZv(); - vertex.Update(xyz); - if (vertex.GetRr()1.) continue; - if (vertex.GetDcaV0Daughters()>0.3*vertex.GetRr()) continue; - // if (vertex.GetPointAngle()<0.9) continue; - vertex.SetIndex(0,itrack0); - vertex.SetIndex(1,itrack1); - fV0s->AddLast(new AliV0(vertex)); + fV0Tree->SetDirectory(0); + fV0Tree->Branch("v0.",&v0); + fV0Tree->Branch("t0.",&track0); + fV0Tree->Branch("t1.",&track1); + fV0Tree->Branch("ft0.",&ftrack0); + fV0Tree->Branch("ft1.",&ftrack1); + fV0Tree->Branch("s0.",&seed0); + fV0Tree->Branch("s1.",&seed1); + }else{ + fV0Tree->SetBranchAddress("v0.",&v0); + fV0Tree->SetBranchAddress("t0.",&track0); + fV0Tree->SetBranchAddress("t1.",&track1); + fV0Tree->SetBranchAddress("ft0.",&ftrack0); + fV0Tree->SetBranchAddress("ft1.",&ftrack1); + fV0Tree->SetBranchAddress("s0.",&seed0); + fV0Tree->SetBranchAddress("s1.",&seed1); } + fV0Tree->Fill(); } - ProcessV0(2); - for (Int_t i=0;iGetEntries(); i++) delete fV0s->At(i); - fV0s->Clear(); } +Long64_t AliTPCcalibV0::Merge(TCollection *const li) { + TIterator* iter = li->MakeIterator(); + AliTPCcalibV0* cal = 0; + while ((cal = (AliTPCcalibV0*)iter->Next())) { + if (cal->fV0Tree){ + if (!fV0Tree) { + fV0Tree = new TTree("V0s","V0s"); + fV0Tree->SetDirectory(0); + } + if (cal->fV0Tree->GetEntries()>0) AliTPCcalibV0::AddTree(cal->fV0Tree); + if (cal->fHPTTree->GetEntries()>0) AliTPCcalibV0::AddTreeHPT(cal->fHPTTree); + } + } + return 0; +} -// void AliTPCcalibV0::ProcessV0(Int_t ftype){ -// // -// // -// const Double_t ktimeK0 = 2.684; -// const Double_t ktimeLambda = 7.89; - - -// if (! fGammas) fGammas = new TObjArray(10); -// fGammas->Clear(); -// Int_t nV0s = fV0s->GetEntries(); -// if (nV0s==0) return; -// AliKFVertex primVtx(*(fESD->GetPrimaryVertex())); -// // -// for (Int_t ivertex=0; ivertexAt(ivertex); -// AliESDtrack * trackN = fESD->GetTrack(v0->GetIndex(0)); -// AliESDtrack * trackP = fESD->GetTrack(v0->GetIndex(1)); -// // -// // -// // -// AliKFParticle *v0K0 = Fit(primVtx,v0,211,211); -// AliKFParticle *v0Gamma = Fit(primVtx,v0,11,-11); -// AliKFParticle *v0Lambda42 = Fit(primVtx,v0,2212,211); -// AliKFParticle *v0Lambda24 = Fit(primVtx,v0,211,2212); -// //Set production vertex -// v0K0->SetProductionVertex( primVtx ); -// v0Gamma->SetProductionVertex( primVtx ); -// v0Lambda42->SetProductionVertex( primVtx ); -// v0Lambda24->SetProductionVertex( primVtx ); -// Double_t massK0, massGamma, massLambda42,massLambda24, massSigma; -// v0K0->GetMass( massK0,massSigma); -// v0Gamma->GetMass( massGamma,massSigma); -// v0Lambda42->GetMass( massLambda42,massSigma); -// v0Lambda24->GetMass( massLambda24,massSigma); -// Float_t chi2K0 = v0K0->GetChi2()/v0K0->GetNDF(); -// Float_t chi2Gamma = v0Gamma->GetChi2()/v0Gamma->GetNDF(); -// Float_t chi2Lambda42 = v0Lambda42->GetChi2()/v0Lambda42->GetNDF(); -// Float_t chi2Lambda24 = v0Lambda24->GetChi2()/v0Lambda24->GetNDF(); -// // -// // Mass Contrained params -// // -// AliKFParticle *v0K0C = Fit(primVtx,v0,211,211); -// AliKFParticle *v0GammaC = Fit(primVtx,v0,11,-11); -// AliKFParticle *v0Lambda42C = Fit(primVtx,v0,2212,211); -// AliKFParticle *v0Lambda24C = Fit(primVtx,v0,211,2212); -// // -// v0K0C->SetProductionVertex( primVtx ); -// v0GammaC->SetProductionVertex( primVtx ); -// v0Lambda42C->SetProductionVertex( primVtx ); -// v0Lambda24C->SetProductionVertex( primVtx ); - -// v0K0C->SetMassConstraint(fPdg->GetParticle(310)->Mass()); -// v0GammaC->SetMassConstraint(0); -// v0Lambda42C->SetMassConstraint(fPdg->GetParticle(3122)->Mass()); -// v0Lambda24C->SetMassConstraint(fPdg->GetParticle(-3122)->Mass()); -// // -// Double_t timeK0, sigmaTimeK0; -// Double_t timeLambda42, sigmaTimeLambda42; -// Double_t timeLambda24, sigmaTimeLambda24; -// v0K0C->GetLifeTime(timeK0, sigmaTimeK0); -// //v0K0Gamma->GetLifeTime(timeK0, sigmaTimeK0); -// v0Lambda42C->GetLifeTime(timeLambda42, sigmaTimeLambda42); -// v0Lambda24C->GetLifeTime(timeLambda24, sigmaTimeLambda24); - - -// // -// Float_t chi2K0C = v0K0C->GetChi2()/v0K0C->GetNDF(); -// if (chi2K0C<0) chi2K0C=100; -// Float_t chi2GammaC = v0GammaC->GetChi2()/v0GammaC->GetNDF(); -// if (chi2GammaC<0) chi2GammaC=100; -// Float_t chi2Lambda42C = v0Lambda42C->GetChi2()/v0Lambda42C->GetNDF(); -// if (chi2Lambda42C<0) chi2Lambda42C=100; -// Float_t chi2Lambda24C = v0Lambda24C->GetChi2()/v0Lambda24C->GetNDF(); -// if (chi2Lambda24C<0) chi2Lambda24C=100; -// // -// Float_t minChi2C=99; -// Int_t type =-1; -// if (chi2K0CGetP()/fPdg->GetParticle(-2212)->Mass(); -// Float_t betaGammaPi = trackN->GetP()/fPdg->GetParticle(-211)->Mass(); -// Float_t betaGammaEl = trackN->GetP()/fPdg->GetParticle(11)->Mass(); -// Float_t dedxTeorP = BetheBlochAleph(betaGammaP); -// Float_t dedxTeorPi = BetheBlochAleph(betaGammaPi);; -// Float_t dedxTeorEl = BetheBlochAleph(betaGammaEl);; -// // -// // -// if (minChi2>50) continue; -// (*fDebugStream)<<"V0"<< -// "ftype="<AddLast(v0); -// // -// // -// // -// delete v0K0; -// delete v0Gamma; -// delete v0Lambda42; -// delete v0Lambda24; -// delete v0K0C; -// delete v0GammaC; -// delete v0Lambda42C; -// delete v0Lambda24C; -// } -// ProcessPI0(); -// } - +void AliTPCcalibV0::AddTree(TTree * treeInput){ + // + // Add the content of tree: + // Notice automatic copy of tree in ROOT does not work for such complicated tree + // + return ; + AliESDv0 * v0 = new AliESDv0; + Double_t kMinPt=0.8; + AliESDtrack * track0 = 0; // negative track + AliESDtrack * track1 = 0; // positive track + AliESDfriendTrack *ftrack0 = 0; + AliESDfriendTrack *ftrack1 = 0; + AliTPCseed *seed0 = 0; + AliTPCseed *seed1 = 0; + treeInput->SetBranchStatus("ft0.",kFALSE); + treeInput->SetBranchStatus("ft1.",kFALSE); + TDatabasePDG pdg; + Double_t massK0= pdg.GetParticle("K0")->Mass(); + Double_t massLambda= pdg.GetParticle("Lambda0")->Mass(); + + Int_t entries= treeInput->GetEntries(); + for (Int_t i=0; iSetBranchAddress("v0.",&v0); + treeInput->SetBranchAddress("t0.",&track0); + treeInput->SetBranchAddress("t1.",&track1); + treeInput->SetBranchAddress("ft0.",&ftrack0); + treeInput->SetBranchAddress("ft1.",&ftrack1); + treeInput->SetBranchAddress("s0.",&seed0); + treeInput->SetBranchAddress("s1.",&seed1); + if (fV0Tree->GetEntries()==0){ + fV0Tree->SetDirectory(0); + fV0Tree->Branch("v0.",&v0); + fV0Tree->Branch("t0.",&track0); + fV0Tree->Branch("t1.",&track1); + fV0Tree->Branch("ft0.",&ftrack0); + fV0Tree->Branch("ft1.",&ftrack1); + fV0Tree->Branch("s0.",&seed0); + fV0Tree->Branch("s1.",&seed1); + }else{ + fV0Tree->SetBranchAddress("v0.",&v0); + fV0Tree->SetBranchAddress("t0.",&track0); + fV0Tree->SetBranchAddress("t1.",&track1); + fV0Tree->SetBranchAddress("ft0.",&ftrack0); + fV0Tree->SetBranchAddress("ft1.",&ftrack1); + fV0Tree->SetBranchAddress("s0.",&seed0); + fV0Tree->SetBranchAddress("s1.",&seed1); + } + // + treeInput->GetEntry(i); + //ftrack0->GetCalibContainer()->SetOwner(kTRUE); + //ftrack1->GetCalibContainer()->SetOwner(kTRUE); + Bool_t isOK=kTRUE; + if (v0->GetOnFlyStatus()==kFALSE) isOK=kFALSE; + if (track0->GetTPCncls()<100) isOK=kFALSE; + if (track1->GetTPCncls()<100) isOK=kFALSE; + if (TMath::Min(seed0->Pt(),seed1->Pt())Pt(),track1->Pt())GetEffMass(2,2)-massK0)<0.05) isV0=kTRUE; + if (TMath::Abs(v0->GetEffMass(4,2)-massLambda)<0.05) isV0=kTRUE; + if (TMath::Abs(v0->GetEffMass(2,4)-massLambda)<0.05) isV0=kTRUE; + if (TMath::Abs(v0->GetEffMass(0,0))<0.02) isV0=kFALSE; //reject electrons + if (!isV0) isOK=kFALSE; + if (isOK) fV0Tree->Fill(); + delete v0; + delete track0; + delete track1; + delete ftrack0; + delete ftrack1; + delete seed0; + delete seed1; + v0=0; + track0=0; + track1=0; + ftrack0=0; + ftrack1=0; + seed0=0; + seed1=0; + } +} +void AliTPCcalibV0::AddTreeHPT(TTree * treeInput){ + // + // Add the content of tree: + // Notice automatic copy of tree in ROOT does not work for such complicated tree + // + return ; + AliESDtrack *track = 0; + AliESDfriendTrack *friendTrack = 0; + AliTPCseed *seed = 0; + if (!treeInput) return; + if (treeInput->GetEntries()==0) return; + // + Int_t entries= treeInput->GetEntries(); + // + for (Int_t i=0; iSetBranchAddress("t.",&track); + treeInput->SetBranchAddress("ft.",&friendTrack); + treeInput->SetBranchAddress("s.",&seed); + treeInput->GetEntry(i); + // + if (fHPTTree->GetEntries()==0){ + fHPTTree->SetDirectory(0); + fHPTTree->Branch("t.",&track); + fHPTTree->Branch("ft.",&friendTrack); + fHPTTree->Branch("s.",&seed); + }else{ + fHPTTree->SetBranchAddress("t.",&track); + fHPTTree->SetBranchAddress("ft.",&friendTrack); + fHPTTree->SetBranchAddress("s.",&seed); + } + Bool_t isOK=kTRUE; + if (!track->IsOn(AliESDtrack::kITSrefit)) isOK=kFALSE; + if (!track->IsOn(AliESDtrack::kTOFout)) isOK=kFALSE; + if (isOK) fHPTTree->Fill(); + // + delete track; + delete friendTrack; + delete seed; + } +} -void AliTPCcalibV0::ProcessV0(Int_t ftype){ +void AliTPCcalibV0::MakeFitTreeTrack(const TObjArray * corrArray, Double_t ptCut, Int_t /*run*/){ // + // Make a fit tree // - - if (! fGammas) fGammas = new TObjArray(10); - fGammas->Clear(); - Int_t nV0s = fV0s->GetEntries(); - if (nV0s==0) return; - AliKFVertex primVtx(*(fESD->GetPrimaryVertex())); + // 0. Loop over selected tracks + // 1. Loop over all transformation - refit the track with and without the + // transformtation + // 2. Dump the matching paramaeters to the debugStremer // - for (Int_t ivertex=0; ivertexAt(ivertex); - AliESDtrack * trackN = fESD->GetTrack(v0->GetIndex(0)); // negative track - AliESDtrack * trackP = fESD->GetTrack(v0->GetIndex(1)); // positive track - - const AliExternalTrackParam * paramInNeg = trackN->GetInnerParam(); - const AliExternalTrackParam * paramInPos = trackP->GetInnerParam(); - if (!paramInPos) continue; // in case the inner paramters do not exist - if (!paramInNeg) continue; - // - // + //Connect input + const Int_t kMinNcl=120; + TFile f("TPCV0Objects.root"); + AliTPCcalibV0 *v0TPC = (AliTPCcalibV0*) f.Get("v0TPC"); + TTree * treeInput = v0TPC->GetHPTTree(); + TTreeSRedirector *pcstream = new TTreeSRedirector("fitHPT.root"); + AliESDtrack *track = 0; + AliESDfriendTrack *friendTrack = 0; + AliTPCseed *seed = 0; + if (!treeInput) return; + if (treeInput->GetEntries()==0) return; + // + treeInput->SetBranchAddress("t.",&track); + treeInput->SetBranchAddress("ft.",&friendTrack); + treeInput->SetBranchAddress("s.",&seed); + // + Int_t ncorr=0; + if (corrArray) ncorr = corrArray->GetEntries(); + AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ; + // AliTPCParam *param = AliTPCcalibDB::Instance()->GetParameters(); +// AliGRPObject* grp = AliTPCcalibDB::Instance()->GetGRP(run); +// Double_t time=0.5*(grp->GetTimeStart() +grp->GetTimeEnd()); + // + // + // + Int_t ntracks= treeInput->GetEntries(); + for (Int_t itrack=0; itrackGetEntry(itrack); + if (!track) continue; + if (seed->Pt()Pt()GetTPCncls()SetProductionVertex( primVtx ); - v0Gamma->SetProductionVertex( primVtx ); - v0Lambda42->SetProductionVertex( primVtx ); - v0Lambda24->SetProductionVertex( primVtx ); - Double_t massK0, massGamma, massLambda42,massLambda24, massSigma; - v0K0->GetMass( massK0,massSigma); - v0Gamma->GetMass( massGamma,massSigma); - v0Lambda42->GetMass( massLambda42,massSigma); - v0Lambda24->GetMass( massLambda24,massSigma); - Float_t chi2K0 = v0K0->GetChi2()/v0K0->GetNDF(); - Float_t chi2Gamma = v0Gamma->GetChi2()/v0Gamma->GetNDF(); - Float_t chi2Lambda42 = v0Lambda42->GetChi2()/v0Lambda42->GetNDF(); - Float_t chi2Lambda24 = v0Lambda24->GetChi2()/v0Lambda24->GetNDF(); + // Reapply transformation // - // Mass Contrained params + for (Int_t irow=0; irow<159; irow++){ + AliTPCclusterMI *cluster=seed->GetClusterPointer(irow); + if (cluster &&cluster->GetX()>10){ + Double_t x0[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()}; + Int_t index0[1]={cluster->GetDetector()}; + transform->Transform(x0,index0,0,1); + cluster->SetX(x0[0]); + cluster->SetY(x0[1]); + cluster->SetZ(x0[2]); + // + } + } // - AliKFParticle *v0K0C = Fit(primVtx,v0,-211,211); - AliKFParticle *v0GammaC = Fit(primVtx,v0,11,-11); - AliKFParticle *v0Lambda42C = Fit(primVtx,v0,-2212,211); //lambdaBar - AliKFParticle *v0Lambda24C = Fit(primVtx,v0,-211,2212); //lambda - // - v0K0C->SetProductionVertex( primVtx ); - v0GammaC->SetProductionVertex( primVtx ); - v0Lambda42C->SetProductionVertex( primVtx ); - v0Lambda24C->SetProductionVertex( primVtx ); - - v0K0C->SetMassConstraint(fPdg->GetParticle(310)->Mass()); - v0GammaC->SetMassConstraint(0); - v0Lambda42C->SetMassConstraint(fPdg->GetParticle(-3122)->Mass()); - v0Lambda24C->SetMassConstraint(fPdg->GetParticle(3122)->Mass()); - // - Double_t timeK0, sigmaTimeK0; - Double_t timeLambda42, sigmaTimeLambda42; - Double_t timeLambda24, sigmaTimeLambda24; - v0K0C->GetLifeTime(timeK0, sigmaTimeK0); - //v0K0Gamma->GetLifeTime(timeK0, sigmaTimeK0); - v0Lambda42C->GetLifeTime(timeLambda42, sigmaTimeLambda42); - v0Lambda24C->GetLifeTime(timeLambda24, sigmaTimeLambda24); - + AliExternalTrackParam* paramInner=0; + AliExternalTrackParam* paramOuter=0; + AliExternalTrackParam* paramIO=0; + Bool_t isOK=kTRUE; + for (Int_t icorr=-1; icorr=0) corr = (AliTPCCorrection*)corrArray->At(icorr); + AliExternalTrackParam * trackInner = RefitTrack(seed, corr,85,134,0.1); + AliExternalTrackParam * trackIO = RefitTrack(seed, corr,245,85,0.1); + AliExternalTrackParam * trackOuter = RefitTrack(seed, corr,245,134,0.1 ); + if (trackInner&&trackOuter&&trackIO){ + trackOuter->Rotate(trackInner->GetAlpha()); + trackOuter->PropagateTo(trackInner->GetX(),AliTracker::GetBz()); + if (icorr<0) { + paramInner=trackInner; + paramOuter=trackOuter; + paramIO=trackIO; + paramIO->Rotate(seed->GetAlpha()); + paramIO->PropagateTo(seed->GetX(),AliTracker::GetBz()); + } + }else{ + isOK=kFALSE; + } + + } + if (paramOuter&& paramInner) { + // Bool_t isOK=kTRUE; + if (paramInner->GetSigmaY2()>0.01) isOK&=kFALSE; + if (paramOuter->GetSigmaY2()>0.01) isOK&=kFALSE; + if (paramInner->GetSigmaZ2()>0.01) isOK&=kFALSE; + if (paramOuter->GetSigmaZ2()>0.01) isOK&=kFALSE; + (*pcstream)<<"fit"<< + "s.="<Init(); + cc->Print("DA"); // Print used correction classes + TObjArray *array = cc->GetCorrections() + AliTPCcalibV0::MakeFitTreeTrack(array,2,run); + + */ +} +void AliTPCcalibV0::MakeFitTreeV0(const TObjArray * corrArray, Double_t ptCut, Int_t run){ + // + // Make a fit tree + // + // 0. Loop over selected tracks + // 1. Loop over all transformation - refit the track with and without the + // transformtation + // 2. Dump the matching paramaeters to the debugStremer + // + + //Connect input + TFile f("TPCV0Objects.root"); + AliTPCcalibV0 *v0TPC = (AliTPCcalibV0*) f.Get("v0TPC"); + TTree * treeInput = v0TPC->GetV0Tree(); + TTreeSRedirector *pcstream = new TTreeSRedirector("fitV0.root"); + AliESDv0 *v0 = 0; + AliESDtrack *track0 = 0; + AliESDfriendTrack *friendTrack0 = 0; + AliTPCseed *seed0 = 0; + AliTPCseed *s0 = 0; + AliESDtrack *track1 = 0; + AliESDfriendTrack *friendTrack1 = 0; + AliTPCseed *s1 = 0; + AliTPCseed *seed1 = 0; + if (!treeInput) return; + if (treeInput->GetEntries()==0) return; + // + treeInput->SetBranchAddress("v0.",&v0); + treeInput->SetBranchAddress("t0.",&track0); + treeInput->SetBranchAddress("ft0.",&friendTrack0); + treeInput->SetBranchAddress("s0.",&s0); + treeInput->SetBranchAddress("t1.",&track1); + treeInput->SetBranchAddress("ft1.",&friendTrack1); + treeInput->SetBranchAddress("s1.",&s1); + // + TDatabasePDG pdg; + Int_t ncorr=0; + if (corrArray) ncorr = corrArray->GetEntries(); + AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ; + Double_t massK0= pdg.GetParticle("K0")->Mass(); + Double_t massLambda= pdg.GetParticle("Lambda0")->Mass(); + Double_t massPion=pdg.GetParticle("pi+")->Mass(); + Double_t massProton=pdg.GetParticle("proton")->Mass(); + Int_t pdgPion=pdg.GetParticle("pi+")->PdgCode(); + Int_t pdgProton=pdg.GetParticle("proton")->PdgCode(); + Double_t rmass0=0; + Double_t rmass1=0; + Double_t massV0=0; + Int_t pdg0=0; + Int_t pdg1=0; + // + // + // + Int_t nv0s= treeInput->GetEntries(); + for (Int_t iv0=0; iv0GetEntry(iv0); + if (TMath::Abs(v0->GetEffMass(2,2)-massK0)<0.03) {isK0=1; v0Type=1;} //select K0s + if (TMath::Abs(v0->GetEffMass(4,2)-massLambda)<0.01) {isLambda=1; v0Type=2;} //select Lambda + if (TMath::Abs(v0->GetEffMass(2,4)-massLambda)<0.01) {isAntiLambda=1;v0Type=3;} //select Anti Lambda + if (isK0+isLambda+isAntiLambda!=1) continue; + rmass0=massPion; + rmass1=massPion; + pdg0=pdgPion; + pdg1=pdgPion; + if (isLambda) {rmass0=massProton; pdg0=pdgProton;} + if (isAntiLambda) {rmass1=massProton; pdg1=pdgProton;} + massV0=massK0; + if (isK0==0) massV0=massLambda; // - Float_t chi2K0C = v0K0C->GetChi2()/v0K0C->GetNDF(); - if (chi2K0C<0) chi2K0C=100; - Float_t chi2GammaC = v0GammaC->GetChi2()/v0GammaC->GetNDF(); - if (chi2GammaC<0) chi2GammaC=100; - Float_t chi2Lambda42C = v0Lambda42C->GetChi2()/v0Lambda42C->GetNDF(); - if (chi2Lambda42C<0) chi2Lambda42C=100; - Float_t chi2Lambda24C = v0Lambda24C->GetChi2()/v0Lambda24C->GetNDF(); - if (chi2Lambda24C<0) chi2Lambda24C=100; + if (!s0) continue; + seed0=(s0->GetSign()>0)?s0:s1; + seed1=(s0->GetSign()>0)?s1:s0; + if (seed0->GetZ()*seed1->GetZ()<0) continue; //remove membrane crossed tracks + if (seed0->Pt()Pt()GetP()/fPdg->GetParticle(-211)->Mass(); - betaGamma1 = paramInPos->GetP()/fPdg->GetParticle(211)->Mass(); - break; - case 1: - betaGamma0 = paramInNeg->GetP()/fPdg->GetParticle(11)->Mass(); - betaGamma1 = paramInPos->GetP()/fPdg->GetParticle(-11)->Mass(); - break; - case 2: - betaGamma0 = paramInNeg->GetP()/fPdg->GetParticle(-2212)->Mass(); - betaGamma1 = paramInPos->GetP()/fPdg->GetParticle(211)->Mass(); - break; - case 3: - betaGamma0 = paramInNeg->GetP()/fPdg->GetParticle(-211)->Mass(); - betaGamma1 = paramInPos->GetP()/fPdg->GetParticle(2212)->Mass(); - break; - } - - // cuts and histogram filling - Int_t numCand = 0; // number of particle types which have a chi2 < 10*minChi2 - - if (minChi2C < 2 && ftype == 1) { - // - if (chi2K0C < 10*minChi2C) numCand++; - if (chi2GammaC < 10*minChi2C) numCand++; - if (chi2Lambda42C < 10*minChi2C) numCand++; - if (chi2Lambda24C < 10*minChi2C) numCand++; - // - if (numCand < 2) { - if (paramInNeg->GetP() > 0.4) fTPCdEdx->Fill(betaGamma0, trackN->GetTPCsignal()); - if (paramInPos->GetP() > 0.4) fTPCdEdx->Fill(betaGamma1, trackP->GetTPCsignal()); - } - } - + // Reapply transformation // + for (Int_t itype=0; itype<2; itype++){ + AliTPCseed * seed = (itype==0) ? seed0: seed1; + for (Int_t irow=0; irow<159; irow++){ + AliTPCclusterMI *cluster=seed->GetClusterPointer(irow); + if (cluster &&cluster->GetX()>10){ + Double_t x0[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()}; + Int_t index0[1]={cluster->GetDetector()}; + transform->Transform(x0,index0,0,1); + cluster->SetX(x0[0]); + cluster->SetY(x0[1]); + cluster->SetZ(x0[2]); + // + } + } + } + Bool_t isOK=kTRUE; + Double_t radius = v0->GetRr(); + Double_t xyz[3]; + v0->GetXYZ(xyz[0],xyz[1],xyz[2]); + Double_t alpha = TMath::ATan2(xyz[1],xyz[0]); + TObjArray arrayV0in(ncorr+1); + TObjArray arrayV0io(ncorr+1); + TObjArray arrayT0(ncorr+1); + TObjArray arrayT1(ncorr+1); + arrayV0in.SetOwner(kTRUE); + arrayV0io.SetOwner(kTRUE); // - // write output tree - if (minChi2>50) continue; - TTreeSRedirector *cstream = GetDebugStreamer(); - if (cstream){ - (*cstream)<<"V0"<< - "ftype="<Clone(),icorr+1); + arrayT1.AddAt(trackIO1->Clone(),icorr+1); + Int_t charge=TMath::Nint(trackIO0->GetSign()); + AliKFParticle pin0( *trackInner0, pdg0*charge); + AliKFParticle pin1( *trackInner1, -pdg1*charge); + AliKFParticle pio0( *trackIO0, pdg0*charge); + AliKFParticle pio1( *trackIO1, -pdg1*charge); + AliKFParticle v0in; + AliKFParticle v0io; + v0in+=pin0; + v0in+=pin1; + v0io+=pio0; + v0io+=pio1; + arrayV0in.AddAt(v0in.Clone(),icorr+1); + arrayV0io.AddAt(v0io.Clone(),icorr+1); + } } - if (type==1) fGammas->AddLast(v0); + if (!isOK) continue; // + //AliKFParticle* pin0= (AliKFParticle*)arrayV0in.At(0); + AliKFParticle* pio0= (AliKFParticle*)arrayV0io.At(0); + AliExternalTrackParam *param0=(AliExternalTrackParam *)arrayT0.At(0); + AliExternalTrackParam *param1=(AliExternalTrackParam *)arrayT1.At(0); + Double_t mass0=0, mass0E=0; + pio0->GetMass( mass0,mass0E); // + Double_t mean=mass0-massV0; + if (TMath::Abs(mean)>0.05) continue; + Double_t mass[10000]; // - delete v0K0; - delete v0Gamma; - delete v0Lambda42; - delete v0Lambda24; - delete v0K0C; - delete v0GammaC; - delete v0Lambda42C; - delete v0Lambda24C; + Int_t dtype=30; // id for V0 + Int_t ptype=5; // id for invariant mass + // Int_t id=TMath::Nint(100.*(param0->Pt()-param1->Pt())/(param0->Pt()+param1->Pt())); // K0s V0 asymetry + Int_t id=Int_t(1000.*(param0->Pt()-param1->Pt())); // K0s V0 asymetry + Double_t gx,gy,gz, px,py,pz; + Double_t pt = v0->Pt(); + v0->GetXYZ(gx,gy,gz); + v0->GetPxPyPz(px,py,pz); + Double_t theta=pz/TMath::Sqrt(px*px+py*py); + Double_t phi=TMath::ATan2(py,px); + Double_t snp=0.5*(seed0->GetSnp()+seed1->GetSnp()); + Double_t sector=9*phi/TMath::Pi(); + Double_t dsec=sector-TMath::Nint(sector); + Double_t refX=TMath::Sqrt(gx*gx+gy*gy); + //Int_t nentries=v0Type; + Double_t bz=AliTracker::GetBz(); + Double_t dRrec=0; + (*pcstream)<<"fitDebug"<< + "id="<GetEntries(); - if (nentries<2) return; + // Refit the track: + // seed - tpc track with cluster + // corr - distrotion/correction class - apllied to the points + // xstart - radius to start propagate/update + // xstop - radius to stop propagate/update // - Double_t m0[3], m1[3]; - AliKFVertex primVtx(*(fESD->GetPrimaryVertex())); - for (Int_t i0=0; i0At(i0); - v00->GetPxPyPz (m0[0], m0[1], m0[2]); - AliKFParticle *p00 = Fit(primVtx, v00, 11,-11); - p00->SetProductionVertex( primVtx ); - p00->SetMassConstraint(0); - // - for (Int_t i1=i0; i1At(i1); - v01->GetPxPyPz (m1[0], m1[1], m1[2]); - AliKFParticle *p01 = Fit(primVtx, v01, 11,-11); - p01->SetProductionVertex( primVtx ); - p01->SetMassConstraint(0); - if (v00->GetIndex(0) != v01->GetIndex(0) && - v00->GetIndex(1) != v01->GetIndex(1)){ - AliKFParticle pi0( *p00,*p01); - pi0.SetProductionVertex(primVtx); - Double_t n1 = TMath::Sqrt (m0[0]*m0[0] + m0[1]*m0[1] + m0[2]*m0[2]); - Double_t n2 = TMath::Sqrt (m1[0]*m1[0] + m1[1]*m1[1] + m1[2]*m1[2]); - Double_t mass = TMath::Sqrt(2.*(n1*n2 - (m0[0]*m1[0] + m0[1]*m1[1] + m0[2]*m1[2]))); - TTreeSRedirector *cstream = GetDebugStreamer(); - if (cstream){ - (*cstream)<<"PI0"<< - "v00.="<PropagateTo(xstart, AliTracker::GetBz()); + refit->AddCovariance(covar); + refit->ResetCovariance(kResetCov); + Double_t xmin = TMath::Min(xstart,xstop); + Double_t xmax = TMath::Max(xstart,xstop); + Int_t ncl=0; + // + Bool_t isOK=kTRUE; + for (Int_t index0=0; index0<159; index0++){ + Int_t irow= (xstartGetClusterPointer(irow); //cluster in local system + if (!cluster) continue; + if (cluster->GetX()GetX()>xmax) continue; + Double_t alpha = TMath::Pi()*(cluster->GetDetector()+0.5)/9.; + if (!refit->Rotate(alpha)) isOK=kFALSE; + Double_t x = cluster->GetX(); + Double_t y = cluster->GetY(); + Double_t z = cluster->GetZ(); + if (corr){ + Float_t xyz[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()}; // original position + corr->DistortPointLocal(xyz,cluster->GetDetector()); + x=xyz[0]; + y=xyz[1]; + z=xyz[2]; } - delete p00; + if (!AliTracker::PropagateTrackToBxByBz(refit, x,mass,1.,kFALSE)) isOK=kFALSE; + if (!isOK) continue; + Double_t cov[3]={0.01,0.,0.01}; + Double_t yz[2]={y,z}; + if (!refit->Update(yz,cov)) isOK=kFALSE; + ncl++; } + if (!AliTracker::PropagateTrackToBxByBz(refit, xstop, mass,1.,kTRUE)) isOK=kFALSE; + // + if (!isOK) { + delete refit; + return 0; + } + return refit; } +// +// Obsolete part +// + + + + AliKFParticle * AliTPCcalibV0::Fit(AliKFVertex & /*primVtx*/, AliESDv0 *v0, Int_t PDG1, Int_t PDG2){ // // Make KF Particle @@ -750,13 +860,6 @@ AliKFParticle * AliTPCcalibV0::Fit(AliKFVertex & /*primVtx*/, AliESDv0 *v0, Int_ return V0; } -TH2F * AliTPCcalibV0::GetHistograms() { - // - // - // - return fTPCdEdx; -} - @@ -778,10 +881,80 @@ void AliTPCcalibV0::BinLogX(TH2F *h) { new_bins[i] = factor * new_bins[i-1]; } axis->Set(bins, new_bins); - delete [] new_bins; - + delete [] new_bins; } - +void AliTPCcalibV0::FilterV0s(AliESDEvent* event){ + // + // + TDatabasePDG pdg; + const Double_t kChi2Cut=20; + const Double_t kMinR=2; + const Double_t ptCut=0.2; + const Int_t kMinNcl=110; + // + Int_t nv0 = event->GetNumberOfV0s(); + AliESDVertex *vertex= (AliESDVertex *)event->GetPrimaryVertex(); + AliKFVertex kfvertex=*vertex; + // + for (Int_t iv0=0;iv0GetV0(iv0); + if (!v0) continue; + if (v0->GetPindex()<0) continue; + if (v0->GetNindex()<0) continue; + if (TMath::Max(v0->GetPindex(), v0->GetNindex())>event->GetNumberOfTracks()) continue; + // + // + AliExternalTrackParam pp=(v0->GetParamP()->GetSign()>0) ? (*(v0->GetParamP())):(*(v0->GetParamN())); + AliExternalTrackParam pn=(v0->GetParamP()->GetSign()>0) ? (*(v0->GetParamN())):(*(v0->GetParamP())); + AliKFParticle kfp1( pp, 211 ); + AliKFParticle kfp2( pn, -211 ); + AliKFParticle *v0KFK0 = new AliKFParticle(kfp1,kfp2); + AliKFParticle *v0KFK0CV = new AliKFParticle(*v0KFK0); + v0KFK0CV->SetProductionVertex(kfvertex); + v0KFK0CV->TransportToProductionVertex(); + AliKFParticle *v0KFK0CVM = new AliKFParticle(*v0KFK0CV); + v0KFK0CVM->SetMassConstraint(pdg.GetParticle("K_S0")->Mass()); + Double_t chi2K0 = v0KFK0CV->GetChi2(); + Double_t chi2K0M= v0KFK0CVM->GetChi2(); + Double_t massK0=v0KFK0CV->GetMass(); + if (chi2K0>kChi2Cut) continue; + if (v0->GetRr()>kMinR) continue; + // + Double_t maxPt = TMath::Max(v0->GetParamP()->Pt(), v0->GetParamN()->Pt()); + Double_t effMass22=v0->GetEffMass(2,2); + Double_t effMass42=v0->GetEffMass(4,2); + Double_t effMass24=v0->GetEffMass(2,4); + Bool_t isV0= kFALSE; + isV0|=TMath::Abs(effMass22-pdg.GetParticle("K_S0")->Mass())<0.1; + isV0|=TMath::Abs(effMass42-pdg.GetParticle("Lambda0")->Mass())<0.1; + isV0|=TMath::Abs(effMass24-pdg.GetParticle("Lambda0")->Mass())<0.1; + + Double_t sign= v0->GetParamP()->GetSign()* v0->GetParamN()->GetSign(); + if (sign<0&&v0->GetOnFlyStatus()>0.5&&maxPt>ptCut&&isV0){ + AliESDtrack * trackP = event->GetTrack(v0->GetPindex()); + AliESDtrack * trackN = event->GetTrack(v0->GetNindex()); + if (!trackN) continue; + if (!trackP) continue; + Int_t nclP= (Int_t)trackP->GetTPCClusterInfo(2,1); + Int_t nclN= (Int_t)trackN->GetTPCClusterInfo(2,1); + if (TMath::Min(nclP,nclN)Eta()), TMath::Abs(trackN->Eta())); + Double_t ncls = TMath::Min(TMath::Abs(trackP->GetNcls(0)), TMath::Abs(trackN->GetNcls(0))); + if (eta<0.8&&ncls>2){ + // 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); + (*fDebugStreamer)<<"v0tree"<< + "v0.="<