#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"
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<AliESDfriend*>(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; ipart<entries; ipart++){
- part = fStack->Particle(ipart);
- if (!part) continue;
- if (part->P()<kMinP) continue;
- if (part->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;i<ntracks;++i) {
+ Bool_t isOK=kFALSE;
+ AliESDtrack *track = esd->GetTrack(i);
+ if (track->GetTPCncls()<kMinCluster) continue;
+ if (TMath::Abs(AliTracker::GetBz())>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; ipart<sentries; ipart++){
- part = (TParticle*)fParticles->At(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()<kMinPt) findable = kFALSE;
- if (p0->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()<kMinPt) findable = kFALSE;
- if (TMath::Abs(p1->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<AliTPCseed*>(calibObject))) break;
}
- // (*fDebugStream)<<"MC0"<<
- // "P.="<<part<<
- // "findable="<<findable<<
- // "id0="<<id0<<
- // "id1="<<id1<<
- // "\n";
- if (!findable) continue;
- Float_t minpt = TMath::Min(p0->Pt(), p1->Pt());
- Int_t type=-1;
-
- //
- //
- AliKFVertex primVtx(*(fESD->GetPrimaryVertex()));
- for (Int_t ivertex=0; ivertex<fESD->GetNumberOfV0s(); 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.="<<part<<
- "type="<<type<<
- "chi2="<<chi2<<
- "chi2C="<<chi2C<<
- "minpt="<<minpt<<
- "id0="<<id0<<
- "id1="<<id1<<
- "Pn.="<<pn<<
- "Pp.="<<pp<<
- "tn.="<<trackN<<
- "tp.="<<trackP<<
- "nold.="<<nold<<
- "nnew.="<<nnew<<
- "v0.="<<v0<<
- "v0kf.="<<v0kf<<
- "v0kfc.="<<v0kfc<<
- "\n";
- delete v0kf;
- delete v0kfc;
- //
- }
-
- if (cstream){
- (*cstream)<<"MC"<<
- "P.="<<part<<
- "charge="<<charge<<
- "type="<<type<<
- "minpt="<<minpt<<
- "id0="<<id0<<
- "id1="<<id1<<
- "P0.="<<p0<<
- "P1.="<<p1<<
- "nold="<<nold<<
- "nnew="<<nnew<<
- "\n";
- }
+ 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);
}
- 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; ivertex<fESD->GetNumberOfV0s(); 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; ivertex<fESD->GetNumberOfV0s(); 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<AliESDfriend*>(esd->FindListObject("AliESDfriend"));
//
- Int_t ntracks = fESD->GetNumberOfTracks();
- for (Int_t itrack0=0; itrack0<ntracks; itrack0++){
- AliESDtrack * track0 = fESD->GetTrack(itrack0);
- if (track0->GetSign()>0) continue;
- if ( track0->GetTPCNcls()<kMinCluster) continue;
+
+ for (Int_t ivertex=0; ivertex<nV0s; ivertex++){
+ Bool_t isOK=kFALSE;
+ AliESDv0 * v0 = (AliESDv0*) esd->GetV0(ivertex);
+ AliESDtrack * track0 = fESD->GetTrack(v0->GetIndex(0)); // negative track
+ AliESDtrack * track1 = fESD->GetTrack(v0->GetIndex(1)); // positive track
+ if (track0->GetTPCNcls()<kMinCluster) continue;
if (track0->GetKinkIndex(0)>0) continue;
+ if (track1->GetTPCNcls()<kMinCluster) continue;
+ if (track1->GetKinkIndex(0)>0) continue;
+ if (v0->GetOnFlyStatus()==kFALSE) continue;
+ //
+ if (TMath::Min(track0->Pt(),track1->Pt())<kMinMinPt) continue;
+ //
+ //
+ if (TMath::Max(track0->Pt(),track1->Pt())>kMinPt) isOK=kTRUE;
+ if (gRandom->Rndm()<kDownscale) isOK=kTRUE;
+ if (!isOK) continue;
+ //
+ TString filename(AliAnalysisManager::GetAnalysisManager()->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; itrack1<ntracks; itrack1++){
- AliESDtrack * track1 = fESD->GetTrack(itrack1);
- if (track1->GetSign()<0) continue;
- if ( track1->GetTPCNcls()<kMinCluster) continue;
- if (track1->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<AliTPCseed*>(calibObject))) break;
+ }
+ for (Int_t l=0;(calibObject=ftrack1->GetCalibObject(l));++l) {
+ if ((seed1=dynamic_cast<AliTPCseed*>(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()<kMinR) continue;
- if (vertex.GetDcaV0Daughters()>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;i<fV0s->GetEntries(); 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; ivertex<nV0s; ivertex++){
-// AliESDv0 * v0 = (AliESDv0*)fV0s->At(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 (chi2K0C<minChi2C) { minChi2C= chi2K0C; type=0;}
-// if (chi2GammaC<minChi2C) { minChi2C= chi2GammaC; type=1;}
-// if (chi2Lambda42C<minChi2C) { minChi2C= chi2Lambda42C; type=2;}
-// if (chi2Lambda24C<minChi2C) { minChi2C= chi2Lambda24C; type=3;}
-// Float_t minChi2=99;
-// Int_t type0 =-1;
-// if (chi2K0<minChi2) { minChi2= chi2K0; type0=0;}
-// if (chi2Gamma<minChi2) { minChi2= chi2Gamma; type0=1;}
-// if (chi2Lambda42<minChi2) { minChi2= chi2Lambda42; type0=2;}
-// if (chi2Lambda24<minChi2) { minChi2= chi2Lambda24; type0=3;}
-// Float_t betaGammaP = trackN->GetP()/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="<<ftype<<
-// "v0.="<<v0<<
-// "trackN.="<<trackN<<
-// "trackP.="<<trackP<<
-// //
-// "dedxTeorP="<<dedxTeorP<<
-// "dedxTeorPi="<<dedxTeorPi<<
-// "dedxTeorEl="<<dedxTeorEl<<
-// //
-// "type="<<type<<
-// "chi2C="<<minChi2C<<
-// "v0K0.="<<v0K0<<
-// "v0Gamma.="<<v0Gamma<<
-// "v0Lambda42.="<<v0Lambda42<<
-// "v0Lambda24.="<<v0Lambda24<<
-// //
-// "chi20K0.="<<chi2K0<<
-// "chi2Gamma.="<<chi2Gamma<<
-// "chi2Lambda42.="<<chi2Lambda42<<
-// "chi2Lambda24.="<<chi2Lambda24<<
-// //
-// "chi20K0c.="<<chi2K0C<<
-// "chi2Gammac.="<<chi2GammaC<<
-// "chi2Lambda42c.="<<chi2Lambda42C<<
-// "chi2Lambda24c.="<<chi2Lambda24C<<
-// //
-// "v0K0C.="<<v0K0C<<
-// "v0GammaC.="<<v0GammaC<<
-// "v0Lambda42C.="<<v0Lambda42C<<
-// "v0Lambda24C.="<<v0Lambda24C<<
-// //
-// "massK0="<<massK0<<
-// "massGamma="<<massGamma<<
-// "massLambda42="<<massLambda42<<
-// "massLambda24="<<massLambda24<<
-// //
-// "timeK0="<<timeK0<<
-// "timeLambda42="<<timeLambda42<<
-// "timeLambda24="<<timeLambda24<<
-// "\n";
-// if (type==1) fGammas->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; i<entries; i++){
+ treeInput->SetBranchAddress("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())<kMinPt) isOK=kFALSE;
+ if (TMath::Min(track0->Pt(),track1->Pt())<kMinPt) isOK=kFALSE;
+ Bool_t isV0=kFALSE;
+ if (TMath::Abs(v0->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; i<entries; i++){
+ track=0;
+ friendTrack=0;
+ seed=0;
+ //
+ treeInput->SetBranchAddress("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; ivertex<nV0s; ivertex++){
- AliESDv0 * v0 = (AliESDv0*)fV0s->At(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; itrack<ntracks; itrack++){
+ treeInput->GetEntry(itrack);
+ if (!track) continue;
+ if (seed->Pt()<ptCut) continue;
+ if (track->Pt()<ptCut) continue;
+ if (track->GetTPCncls()<kMinNcl) continue;
//
- 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();
+ // 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<ncorr; icorr++){
+ //
+ AliTPCCorrection *corr = 0;
+ if (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.="<<seed<<
+ "io.="<<paramIO<<
+ "pIn.="<<paramInner<<
+ "pOut.="<<paramOuter;
+ }
+ //
+ }
+ delete pcstream;
+ /*
+ .x ~/rootlogon.C
+ Int_t run=117112;
+ .x ../ConfigCalibTrain.C(run)
+ .L ../AddTaskTPCCalib.C
+ ConfigOCDB(run)
+ TFile fexb("../../RegisterCorrectionExB.root");
+ AliTPCComposedCorrection *cc= (AliTPCComposedCorrection*) fexb.Get("ComposedExB");
+ cc->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; iv0<nv0s; iv0++){
+ Int_t v0Type=0;
+ Int_t isK0=0;
+ Int_t isLambda=0;
+ Int_t isAntiLambda=0;
+ treeInput->GetEntry(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()<ptCut) continue;
+ if (seed1->Pt()<ptCut) continue;
//
- Float_t minChi2C=99;
- Int_t type =-1;
- if (chi2K0C<minChi2C) { minChi2C= chi2K0C; type=0;}
- if (chi2GammaC<minChi2C) { minChi2C= chi2GammaC; type=1;}
- if (chi2Lambda42C<minChi2C) { minChi2C= chi2Lambda42C; type=2;}
- if (chi2Lambda24C<minChi2C) { minChi2C= chi2Lambda24C; type=3;}
- Float_t minChi2=99;
- Int_t type0 =-1;
- if (chi2K0<minChi2) { minChi2= chi2K0; type0=0;}
- if (chi2Gamma<minChi2) { minChi2= chi2Gamma; type0=1;}
- if (chi2Lambda42<minChi2) { minChi2= chi2Lambda42; type0=2;}
- if (chi2Lambda24<minChi2) { minChi2= chi2Lambda24; type0=3;}
-
- // 0 is negative particle; 1 is positive particle
- Float_t betaGamma0 = 0;
- Float_t betaGamma1 = 0;
-
- switch (type) {
- case 0:
- betaGamma0 = paramInNeg->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="<<ftype<<
- "v0.="<<v0<<
- "trackN.="<<trackN<<
- "trackP.="<<trackP<<
- //
- "betaGamma0="<<betaGamma0<<
- "betaGamma1="<<betaGamma1<<
- //
- "type="<<type<<
- "chi2C="<<minChi2C<<
- "v0K0.="<<v0K0<<
- "v0Gamma.="<<v0Gamma<<
- "v0Lambda42.="<<v0Lambda42<<
- "v0Lambda24.="<<v0Lambda24<<
- //
- "chi20K0.="<<chi2K0<<
- "chi2Gamma.="<<chi2Gamma<<
- "chi2Lambda42.="<<chi2Lambda42<<
- "chi2Lambda24.="<<chi2Lambda24<<
- //
- "chi20K0c.="<<chi2K0C<<
- "chi2Gammac.="<<chi2GammaC<<
- "chi2Lambda42c.="<<chi2Lambda42C<<
- "chi2Lambda24c.="<<chi2Lambda24C<<
- //
- "v0K0C.="<<v0K0C<<
- "v0GammaC.="<<v0GammaC<<
- "v0Lambda42C.="<<v0Lambda42C<<
- "v0Lambda24C.="<<v0Lambda24C<<
- //
- "massK0="<<massK0<<
- "massGamma="<<massGamma<<
- "massLambda42="<<massLambda42<<
- "massLambda24="<<massLambda24<<
+ for (Int_t icorr=-1; icorr<ncorr; icorr++){
+ AliTPCCorrection *corr =0;
+ if (icorr>=0) corr = (AliTPCCorrection*)corrArray->At(icorr);
+ //
+ AliExternalTrackParam * trackInner0 = RefitTrack(seed0, corr,160,85,rmass0);
+ AliExternalTrackParam * trackIO0 = RefitTrack(seed0, corr,245,85,rmass0);
+ AliExternalTrackParam * trackInner1 = RefitTrack(seed1, corr,160,85,rmass1);
+ AliExternalTrackParam * trackIO1 = RefitTrack(seed1, corr,245,85,rmass1);
+ if (!trackInner0) isOK=kFALSE;
+ if (!trackInner1) isOK=kFALSE;
+ if (!trackIO0) isOK=kFALSE;
+ if (!trackIO1) isOK=kFALSE;
+ if (isOK){
+ if (!trackInner0->Rotate(alpha)) isOK=kFALSE;
+ if (!trackInner1->Rotate(alpha)) isOK=kFALSE;
+ if (!trackIO0->Rotate(alpha)) isOK=kFALSE;
+ if (!trackIO1->Rotate(alpha)) isOK=kFALSE;
//
- "timeK0="<<timeK0<<
- "timeLambda42="<<timeLambda42<<
- "timeLambda24="<<timeLambda24<<
- "\n";
+ if (!AliTracker::PropagateTrackToBxByBz(trackInner0, radius, rmass0, 1, kFALSE)) isOK=kFALSE;
+ if (!AliTracker::PropagateTrackToBxByBz(trackInner1, radius, rmass1, 1, kFALSE)) isOK=kFALSE;
+ if (!AliTracker::PropagateTrackToBxByBz(trackIO0, radius, rmass0, 1, kFALSE)) isOK=kFALSE;
+ if (!AliTracker::PropagateTrackToBxByBz(trackIO1, radius, rmass1, 1, kFALSE)) isOK=kFALSE;
+ if (!isOK) continue;
+ arrayT0.AddAt(trackIO0->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="<<id<<
+ "v0.="<<v0<<
+ "mean="<<mean<<
+ "rms="<<mass0E<<
+ "pio0.="<<pio0<<
+ "p0.="<<param0<<
+ "p1.="<<param1;
+ (*pcstream)<<"fit"<< // dump valus for fit
+ "run="<<run<< //run number
+ "bz="<<bz<< // magnetic filed used
+ "dtype="<<dtype<< // detector match type 30
+ "ptype="<<ptype<< // parameter type
+ "theta="<<theta<< // theta
+ "phi="<<phi<< // phi
+ "snp="<<snp<< // snp
+ "mean="<<mean<< // mean dist value
+ "rms="<<mass0E<< // rms
+ "sector="<<sector<<
+ "dsec="<<dsec<<
+ //
+ "refX="<<refX<< // reference radius
+ "gx="<<gx<< // global position
+ "gy="<<gy<< // global position
+ "gz="<<gz<< // global position
+ "dRrec="<<dRrec<< // delta Radius in reconstruction
+ "pt="<<pt<< // pt of the particle
+ "id="<<id<< //delta of the momenta
+ "entries="<<v0Type;// type of the V0
+ for (Int_t icorr=0; icorr<ncorr; icorr++){
+ AliTPCCorrection *corr =0;
+ if (icorr>=0) corr = (AliTPCCorrection*)corrArray->At(icorr);
+ // AliKFParticle* pin= (AliKFParticle*)arrayV0in.At(icorr+1);
+ AliKFParticle* pio= (AliKFParticle*)arrayV0io.At(icorr+1);
+ AliExternalTrackParam *par0=(AliExternalTrackParam *)arrayT0.At(icorr+1);
+ AliExternalTrackParam *par1=(AliExternalTrackParam *)arrayT1.At(icorr+1);
+ Double_t massE=0;
+ pio->GetMass( mass[icorr],massE);
+ mass[icorr]-=mass0;
+ (*pcstream)<<"fit"<<
+ Form("%s=",corr->GetName())<<mass[icorr];
+ (*pcstream)<<"fitDebug"<<
+ Form("%s=",corr->GetName())<<mass[icorr]<<
+ Form("%sp0.=",corr->GetName())<<par0<<
+ Form("%sp1=",corr->GetName())<<par1;
+ }
+ (*pcstream)<<"fit"<< "isOK="<<isOK<<"\n";
+ (*pcstream)<<"fitDebug"<< "isOK="<<isOK<<"\n";
}
- ProcessPI0();
+ delete pcstream;
}
-
-
-void AliTPCcalibV0::ProcessPI0(){
- //
+AliExternalTrackParam * AliTPCcalibV0::RefitTrack(AliTPCseed *seed, AliTPCCorrection * corr, Double_t xstart, Double_t xstop, Double_t mass){
//
- //
- Int_t nentries = fGammas->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; i0<nentries; i0++){
- AliESDv0 *v00 = (AliESDv0*)fGammas->At(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; i1<nentries; i1++){
- AliESDv0 *v01 = (AliESDv0*)fGammas->At(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.="<<v00<<
- "v01.="<<v01<<
- "mass="<<mass<<
- "p00.="<<p00<<
- "p01.="<<p01<<
- "pi0.="<<&pi0<<
- "\n";
- }
- }
- delete p01;
+ const Double_t kResetCov=20.;
+ const Double_t kSigma=5.;
+ Double_t covar[15];
+ for (Int_t i=0;i<15;i++) covar[i]=0;
+ covar[0]=kSigma*kSigma;
+ covar[2]=kSigma*kSigma;
+ covar[5]=kSigma*kSigma/Float_t(150*150);
+ covar[9]=kSigma*kSigma/Float_t(150*150);
+ covar[14]=1*1;
+ //
+ AliExternalTrackParam *refit = new AliExternalTrackParam(*seed);
+ refit->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= (xstart<xstop)? index0:159-index0;
+ AliTPCclusterMI *cluster=seed->GetClusterPointer(irow); //cluster in local system
+ if (!cluster) continue;
+ if (cluster->GetX()<xmin) continue;
+ if (cluster->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
return V0;
}
-TH2F * AliTPCcalibV0::GetHistograms() {
- //
- //
- //
- return fTPCdEdx;
-}
-
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;iv0<nv0;iv0++){
+ AliESDv0 *v0 = event->GetV0(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)<kMinNcl) continue;
+ Double_t eta = TMath::Max(TMath::Abs(trackP->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.="<<v0<<
+ "tp.="<<trackP<<
+ "tm.="<<trackN<<
+ //
+ "v.="<<vertex<<
+ "ncls="<<ncls<<
+ "maxPt="<<maxPt<<
+ "\n";
+ }
+ }
+ }
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff