#include "TFile.h"
#include "TTree.h"
+#include "TRandom3.h"
#include "AliAODEvent.h"
#include "AliAODTrack.h"
#include "AliAODVertex.h"
#include "AliPID.h"
#include "AliAODpidUtil.h"
+#include "AliAnalysisUtils.h"
+#include "assert.h"
+#include "AliGenHijingEventHeader.h"
ClassImp(AliFemtoEventReaderAOD)
#if !(ST_NO_NAMESPACES)
- using namespace units;
+using namespace units;
#endif
using namespace std;
double fV1[3];
//____________________________
-//constructor with 0 parameters , look at default settings
+//constructor with 0 parameters , look at default settings
AliFemtoEventReaderAOD::AliFemtoEventReaderAOD():
fNumberofEvent(0),
fCurEvent(0),
fEstEventMult(kCentrality),
fAODpidUtil(0),
fAODheader(0),
- fInputFile(" "),
- fFileName(" "),
+ fInputFile(""),
fTree(0x0),
fAodFile(0x0),
- fMagFieldSign(1),
- fisEPVZ(kTRUE)
+ fMagFieldSign(1),
+ fisEPVZ(kTRUE),
+ fpA2013(kFALSE),
+ fisPileUp(kFALSE),
+ fMVPlp(kFALSE),
+ fMinVtxContr(0),
+ fMinPlpContribMV(0),
+ fMinPlpContribSPD(0),
+ fDCAglobalTrack(kFALSE),
+ fFlatCent(kFALSE)
{
// default constructor
fAllTrue.ResetAllBits(kTRUE);
fEstEventMult(kCentrality),
fAODpidUtil(0),
fAODheader(0),
- fInputFile(" "),
- fFileName(" "),
+ fInputFile(""),
fTree(0x0),
fAodFile(0x0),
- fMagFieldSign(1),
- fisEPVZ(kTRUE)
+ fMagFieldSign(1),
+ fisEPVZ(kTRUE),
+ fpA2013(kFALSE),
+ fisPileUp(kFALSE),
+ fMVPlp(kFALSE),
+ fMinVtxContr(0),
+ fMinPlpContribMV(0),
+ fMinPlpContribSPD(0),
+ fDCAglobalTrack(kFALSE),
+ fFlatCent(kFALSE)
{
// copy constructor
fReadMC = aReader.fReadMC;
fReadV0 = aReader.fReadV0;
fInputFile = aReader.fInputFile;
- fFileName = aReader.fFileName;
fNumberofEvent = aReader.fNumberofEvent;
fCurEvent = aReader.fCurEvent;
fEvent = new AliAODEvent();
fCentRange[1] = aReader.fCentRange[1];
fEstEventMult = aReader.fEstEventMult;
fUsePreCent = aReader.fUsePreCent;
+ fpA2013 = aReader.fpA2013;
+ fisPileUp = aReader.fisPileUp;
+ fMVPlp = aReader.fMVPlp;
+ fMinVtxContr = aReader.fMinVtxContr;
+ fMinPlpContribMV = aReader.fMinPlpContribMV;
+ fMinPlpContribSPD = aReader.fMinPlpContribSPD;
+ fDCAglobalTrack = aReader.fDCAglobalTrack;
+
}
//__________________
//Destructor
{
// assignment operator
if (this == &aReader)
- return *this;
+ return *this;
fInputFile = aReader.fInputFile;
- fFileName = aReader.fFileName;
fNumberofEvent = aReader.fNumberofEvent;
fCurEvent = aReader.fCurEvent;
if (fTree) delete fTree;
fCentRange[1] = aReader.fCentRange[1];
fUsePreCent = aReader.fUsePreCent;
fEstEventMult = aReader.fEstEventMult;
+ fpA2013 = aReader.fpA2013;
+ fisPileUp = aReader.fisPileUp;
+ fMVPlp = aReader.fMVPlp;
+ fMinVtxContr = aReader.fMinVtxContr;
+ fMinPlpContribMV = aReader.fMinPlpContribMV;
+ fMinPlpContribSPD = aReader.fMinPlpContribSPD;
+ fDCAglobalTrack = aReader.fDCAglobalTrack;
+ fFlatCent= aReader.fFlatCent;
return *this;
}
//__________________
void AliFemtoEventReaderAOD::SetInputFile(const char* inputFile)
{
- //setting the name of file where names of AOD file are written
+ //setting the name of file where names of AOD file are written
//it takes only this files which have good trees
char buffer[256];
fInputFile=string(inputFile);
fTree = new TChain("aodTree");
if(infile.good()==true)
- {
- //checking if all give files have good tree inside
- while (infile.eof()==false)
- {
- infile.getline(buffer,256);
- TFile *aodFile=TFile::Open(buffer,"READ");
- if (aodFile!=0x0)
- {
+ {
+ //checking if all give files have good tree inside
+ while (infile.eof()==false)
+ {
+ infile.getline(buffer,256);
+ TFile *aodFile=TFile::Open(buffer,"READ");
+ if (aodFile!=0x0)
+ {
TTree* tree = (TTree*) aodFile->Get("aodTree");
if (tree!=0x0)
- {
- // cout<<"putting file "<<string(buffer)<<" into analysis"<<endl;
- fTree->AddFile(buffer);
- delete tree;
- }
- aodFile->Close();
+ {
+ // cout<<"putting file "<<string(buffer)<<" into analysis"<<endl;
+ fTree->AddFile(buffer);
+ delete tree;
+ }
+ aodFile->Close();
}
- delete aodFile;
- }
+ delete aodFile;
}
+ }
}
AliFemtoEvent* AliFemtoEventReaderAOD::ReturnHbtEvent()
// convert it to AliFemtoEvent and return
// for further analysis
AliFemtoEvent *hbtEvent = 0;
- // cout<<"reader"<<endl;
- if (fCurEvent==fNumberofEvent)//open next file
+ // cout<<"reader"<<endl;
+ if (fCurEvent==fNumberofEvent)//open next file
+ {
+ if(fNumberofEvent==0)
{
- if(fNumberofEvent==0)
- {
- fEvent=new AliAODEvent();
- fEvent->ReadFromTree(fTree);
+ fEvent=new AliAODEvent();
+ fEvent->ReadFromTree(fTree);
- // Check for the existence of the additional information
+ // Check for the existence of the additional information
// fPWG2AODTracks = (TClonesArray *) fEvent->GetList()->FindObject("pwg2aodtracks");
// if (fPWG2AODTracks) {
// cout << "Reading only tracks with the additional information" << endl;
// }
- fNumberofEvent=fTree->GetEntries();
- // cout<<"Number of Entries in file "<<fNumberofEvent<<endl;
- fCurEvent=0;
- }
- else //no more data to read
- {
- // cout<<"no more files "<<hbtEvent<<endl;
- fReaderStatus=1;
- return hbtEvent;
- }
- }
+ fNumberofEvent=fTree->GetEntries();
+ // cout<<"Number of Entries in file "<<fNumberofEvent<<endl;
+ fCurEvent=0;
+ }
+ else //no more data to read
+ {
+ // cout<<"no more files "<<hbtEvent<<endl;
+ fReaderStatus=1;
+ return hbtEvent;
+ }
+ }
- // cout<<"starting to read event "<<fCurEvent<<endl;
+ // cout<<"starting to read event "<<fCurEvent<<endl;
fTree->GetEvent(fCurEvent);//getting next event
// cout << "Read event " << fEvent << " from file " << fTree << endl;
-
- hbtEvent = new AliFemtoEvent;
- CopyAODtoFemtoEvent(hbtEvent);
+ //hbtEvent = new AliFemtoEvent;
+
+ hbtEvent = CopyAODtoFemtoEvent();
fCurEvent++;
- return hbtEvent;
+ return hbtEvent;
}
-void AliFemtoEventReaderAOD::CopyAODtoFemtoEvent(AliFemtoEvent *tEvent)
+AliFemtoEvent* AliFemtoEventReaderAOD::CopyAODtoFemtoEvent()
{
// A function that reads in the AOD event
// and transfers the neccessary information into
// the internal AliFemtoEvent
+ AliFemtoEvent *tEvent = new AliFemtoEvent();
+
// setting global event characteristics
tEvent->SetRunNumber(fEvent->GetRunNumber());
tEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
tEvent->SetZDCParticipants(0);
tEvent->SetTriggerMask(fEvent->GetTriggerMask());
tEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
-
+
// Attempt to access MC header
AliAODMCHeader *mcH;
TClonesArray *mcP=0;
}
}
- tEvent->SetReactionPlaneAngle(fEvent->GetHeader()->GetQTheta(0)/2.0);
-
- Int_t *motherids=0;
- if (mcP) {
- motherids = new Int_t[((AliAODMCParticle *) mcP->At(mcP->GetEntries()-1))->GetLabel()];
- for (int ip=0; ip<mcP->GetEntries(); ip++) motherids[ip] = 0;
-
- // Read in mother ids
- AliAODMCParticle *motherpart;
- for (int ip=0; ip<mcP->GetEntries(); ip++) {
- motherpart = (AliAODMCParticle *) mcP->At(ip);
- if (motherpart->GetDaughter(0) > 0)
- motherids[motherpart->GetDaughter(0)] = ip;
- if (motherpart->GetDaughter(1) > 0)
- motherids[motherpart->GetDaughter(1)] = ip;
+ AliAODHeader * header = dynamic_cast<AliAODHeader*>(fEvent->GetHeader());
+ assert(header&&"Not a standard AOD");
+
+ tEvent->SetReactionPlaneAngle(header->GetQTheta(0)/2.0);
+ // Int_t *motherids=0;
+ // if (mcP) {
+ // const int motherTabSize = ((AliAODMCParticle *) mcP->At(mcP->GetEntries()-1))->GetLabel();
+ // motherids = new int[motherTabSize+1];
+ // for (int ip=0; ip<motherTabSize+1; ip++) motherids[ip] = 0;
+
+ // // Read in mother ids
+ // AliAODMCParticle *motherpart;
+ // for (int ip=0; ip<mcP->GetEntries(); ip++) {
+ // motherpart = (AliAODMCParticle *) mcP->At(ip);
+ // if (motherpart->GetDaughter(0) > 0)
+ // motherids[motherpart->GetDaughter(0)] = ip;
+ // if (motherpart->GetDaughter(1) > 0)
+ // motherids[motherpart->GetDaughter(1)] = ip;
+ // }
+ // }
+
+ //AliAnalysisUtils
+ if(fisPileUp||fpA2013)
+ {
+ AliAnalysisUtils *anaUtil=new AliAnalysisUtils();
+ if(fMinVtxContr)
+ anaUtil->SetMinVtxContr(fMinVtxContr);
+ if(fpA2013)
+ if(anaUtil->IsVertexSelected2013pA(fEvent)==kFALSE)
+ {
+ delete tEvent;
+ return NULL; //Vertex rejection for pA analysis.
+ }
+ if(fMVPlp) anaUtil->SetUseMVPlpSelection(kTRUE);
+ else anaUtil->SetUseMVPlpSelection(kFALSE);
+ if(fMinPlpContribMV) anaUtil->SetMinPlpContribMV(fMinPlpContribMV);
+ if(fMinPlpContribSPD) anaUtil->SetMinPlpContribSPD(fMinPlpContribSPD);
+ if(fisPileUp)
+ if(anaUtil->IsPileUpEvent(fEvent)) { delete tEvent;return NULL;} //Pile-up rejection.
+ delete anaUtil;
}
- }
// Primary Vertex position
- // double fV1[3];
- fEvent->GetPrimaryVertex()->GetPosition(fV1);
+ const AliAODVertex* aodvertex = (AliAODVertex*) fEvent->GetPrimaryVertex();
+ if(!aodvertex || aodvertex->GetNContributors() < 1) { delete tEvent;return NULL;} //Bad vertex, skip event.
+ aodvertex->GetPosition(fV1);
AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
tEvent->SetPrimVertPos(vertex);
-
+
//starting to reading tracks
int nofTracks=0; //number of reconstructed tracks in event
- // Check to see whether the additional info exists
- // if (fPWG2AODTracks)
- // nofTracks=fPWG2AODTracks->GetEntries();
- // else
nofTracks=fEvent->GetNumberOfTracks();
-
+
AliEventplane *ep = fEvent->GetEventplane();
if (ep) {
tEvent->SetEP(ep);
- if (fisEPVZ)
- tEvent->SetReactionPlaneAngle(ep->GetEventplane("V0",fEvent,2));
- else
- tEvent->SetReactionPlaneAngle(ep->GetEventplane("Q"));
+ if (fisEPVZ)
+ tEvent->SetReactionPlaneAngle(ep->GetEventplane("V0",fEvent,2));
+ else
+ tEvent->SetReactionPlaneAngle(ep->GetEventplane("Q"));
}
AliCentrality *cent = fEvent->GetCentrality();
-
+
if (!fEstEventMult && cent && fUsePreCent) {
if ((cent->GetCentralityPercentile("V0M")*10 < fCentRange[0]) ||
- (cent->GetCentralityPercentile("V0M")*10 > fCentRange[1]))
- {
- // cout << "Centrality " << cent->GetCentralityPercentile("V0M") << " outside of preselection range " << fCentRange[0] << " - " << fCentRange[1] << endl;
-
- return;
- }
+ (cent->GetCentralityPercentile("V0M")*10 > fCentRange[1]))
+ {
+ // cout << "Centrality " << cent->GetCentralityPercentile("V0M") << " outside of preselection range " << fCentRange[0] << " - " << fCentRange[1] << endl;
+ delete tEvent;
+ return NULL;
+ }
+ }
+
+ float percent = cent->GetCentralityPercentile("V0M");
+//flatten centrality dist.
+ if(percent < 9){
+ if(fFlatCent){
+ if(RejectEventCentFlat(fEvent->GetMagneticField(),percent)) { delete tEvent; return NULL;}
+ }
}
int realnofTracks=0; // number of track which we use in a analysis
- int tracksPrim=0;
+ int tracksPrim=0;
- int labels[20000];
+ int labels[20000];
for (int il=0; il<20000; il++) labels[il] = -1;
// looking for global tracks and saving their numbers to copy from them PID information to TPC-only tracks in the main loop over tracks
for (int i=0;i<nofTracks;i++) {
- const AliAODTrack *aodtrack=fEvent->GetTrack(i);
+ const AliAODTrack *aodtrack=dynamic_cast<const AliAODTrack*>(fEvent->GetTrack(i));
+ assert(aodtrack&&"Not a standard AOD");
if (!aodtrack->TestFilterBit(fFilterBit)) {
if(aodtrack->GetID() < 0) continue;
labels[aodtrack->GetID()] = i;
int tNormMult = 0;
for (int i=0;i<nofTracks;i++)
- {
- AliFemtoTrack* trackCopy = new AliFemtoTrack();
+ {
+ AliFemtoTrack* trackCopy;// = new AliFemtoTrack();
// if (fPWG2AODTracks) {
// // Read tracks from the additional pwg2 specific AOD part
// // if they exist
-// // Note that in that case all the AOD tracks without the
+// // Note that in that case all the AOD tracks without the
// // additional information will be ignored !
// AliPWG2AODTrack *pwg2aodtrack = (AliPWG2AODTrack *) fPWG2AODTracks->At(i);
// // Getting the AOD track through the ref of the additional info
-// AliAODTrack *aodtrack = pwg2aodtrack->GetRefAODTrack();
+// AliAODTrack *aodtrack = pwg2aodtrack->GetRefAODTrack();
// if (!aodtrack->TestFilterBit(fFilterBit)) {
// delete trackCopy;
// continue;
// }
-
+
// if (aodtrack->IsOn(AliESDtrack::kTPCrefit))
-// if (aodtrack->Chi2perNDF() < 6.0)
+// if (aodtrack->Chi2perNDF() < 6.0)
// if (aodtrack->Eta() < 0.9)
// tNormMult++;
// CopyAODtoFemtoTrack(aodtrack, trackCopy, pwg2aodtrack);
-
+
// if (mcP) {
// // Fill the hidden information with the simulated data
// // Int_t pLabel = aodtrack->GetLabel();
// AliAODMCParticle *tPart = GetParticleWithLabel(mcP, (TMath::Abs(aodtrack->GetLabel())));
// // Check the mother information
-
+
// // Using the new way of storing the freeze-out information
// // Final state particle is stored twice on the stack
// // one copy (mother) is stored with original freeze-out information
// // and is not tracked
// // the other one (daughter) is stored with primary vertex position
// // and is tracked
-
+
// // Freeze-out coordinates
// double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
// fpx = tPart->Xv() - fV1[0];
// fpy *= 1e13;
// fpz *= 1e13;
// fpt *= 1e13;
-
+
// // cout << "Looking for mother ids " << endl;
// if (motherids[TMath::Abs(aodtrack->GetLabel())]>0) {
// // cout << "Got mother id" << endl;
// // It is the same particle
// // Read in the original freeze-out information
// // and convert it from to [fm]
-
-// // EPOS style
+
+// // EPOS style
// // fpx = mother->Xv()*1e13*0.197327;
// // fpy = mother->Yv()*1e13*0.197327;
// // fpz = mother->Zv()*1e13*0.197327;
// // fpt = mother->T() *1e13*0.197327*0.5;
-
-
-// // Therminator style
+
+
+// // Therminator style
// fpx = mother->Xv()*1e13;
// fpy = mother->Yv()*1e13;
// fpz = mother->Zv()*1e13;
// fpt = mother->T() *1e13*3e10;
-
+
// }
// }
-
+
// // if (fRotateToEventPlane) {
// // double tPhi = TMath::ATan2(fpy, fpx);
// // double tRad = TMath::Hypot(fpx, fpy);
-
+
// // fpx = tRad*TMath::Cos(tPhi - tReactionPlane);
// // fpy = tRad*TMath::Sin(tPhi - tReactionPlane);
// // }
// // if (fRotateToEventPlane) {
// // double tPhi = TMath::ATan2(tPart->Py(), tPart->Px());
// // double tRad = TMath::Hypot(tPart->Px(), tPart->Py());
-
+
// // tInfo->SetTrueMomentum(tRad*TMath::Cos(tPhi - tReactionPlane),
// // tRad*TMath::Sin(tPhi - tReactionPlane),
// // tPart->Pz());
// tPart->Pz()*tPart->Pz());
// if (mass2>0.0)
// tInfo->SetMass(TMath::Sqrt(mass2));
-// else
+// else
// tInfo->SetMass(0.0);
-
+
// tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
// trackCopy->SetHiddenInfo(tInfo);
// }
// }
// else {
- // No additional information exists
- // Read in the normal AliAODTracks
+ // No additional information exists
+ // Read in the normal AliAODTracks
- // const AliAODTrack *aodtrack=fEvent->GetTrack(i); // getting the AODtrack directly
- AliAODTrack *aodtrack=fEvent->GetTrack(i); // getting the AODtrack directly
-
+ // const AliAODTrack *aodtrack=dynamic_cast<AliAODTrack*>(fEvent->GetTrack(i));
+ AliAODTrack *aodtrack=dynamic_cast<AliAODTrack*>(fEvent->GetTrack(i));
+ assert(aodtrack&&"Not a standard AOD"); // getting the AODtrack directly
- if (aodtrack->IsPrimaryCandidate()) tracksPrim++;
-
- if (fFilterBit && !aodtrack->TestFilterBit(fFilterBit)) {
- delete trackCopy;
- continue;
- }
- if (fFilterMask && !aodtrack->TestFilterBit(fFilterMask)) {
- delete trackCopy;
- continue;
- }
-
- //counting particles to set multiplicity
- double impact[2];
- double covimpact[3];
- if (aodtrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
- if(impact[0]<0.2 && TMath::Abs(impact[1]+fV1[2])<2.0)
- //if (aodtrack->IsPrimaryCandidate()) //? instead of kinks?
- if (aodtrack->Chi2perNDF() < 4.0)
- if (aodtrack->Pt() > 0.15 && aodtrack->Pt() < 20)
- if (aodtrack->GetTPCNcls() > 70)
- if (aodtrack->Eta() < 0.8)
- tNormMult++;
- }
-
- CopyAODtoFemtoTrack(aodtrack, trackCopy);
-
- // copying PID information from the correspondent track
- // const AliAODTrack *aodtrackpid = fEvent->GetTrack(labels[-1-fEvent->GetTrack(i)->GetID()]);
-
-
- AliAODTrack *aodtrackpid;
- if((fFilterBit == (1 << (7))) || fFilterMask==128) //for TPC Only tracks we have to copy PID information from corresponding global tracks
- aodtrackpid = fEvent->GetTrack(labels[-1-fEvent->GetTrack(i)->GetID()]);
- else
- aodtrackpid = fEvent->GetTrack(i);
- CopyPIDtoFemtoTrack(aodtrackpid, trackCopy);
-
- if (mcP) {
- // Fill the hidden information with the simulated data
- // Int_t pLabel = aodtrack->GetLabel();
- AliAODMCParticle *tPart = GetParticleWithLabel(mcP, (TMath::Abs(aodtrack->GetLabel())));
-
- AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
- double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
- if (!tPart) {
- fpx = fV1[0];
- fpy = fV1[1];
- fpz = fV1[2];
- tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
- tInfo->SetPDGPid(0);
- tInfo->SetTrueMomentum(0.0, 0.0, 0.0);
- tInfo->SetEmissionPoint(0.0, 0.0, 0.0, 0.0);
- tInfo->SetMass(0);
- }
- else {
- // Check the mother information
-
- // Using the new way of storing the freeze-out information
- // Final state particle is stored twice on the stack
- // one copy (mother) is stored with original freeze-out information
- // and is not tracked
- // the other one (daughter) is stored with primary vertex position
- // and is tracked
-
- // Freeze-out coordinates
- fpx = tPart->Xv() - fV1[0];
- fpy = tPart->Yv() - fV1[1];
- fpz = tPart->Zv() - fV1[2];
- // fpt = tPart->T();
-
- tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
-
- fpx *= 1e13;
- fpy *= 1e13;
- fpz *= 1e13;
- // fpt *= 1e13;
-
- // cout << "Looking for mother ids " << endl;
- if (motherids[TMath::Abs(aodtrack->GetLabel())]>0) {
- // cout << "Got mother id" << endl;
- AliAODMCParticle *mother = GetParticleWithLabel(mcP, motherids[TMath::Abs(aodtrack->GetLabel())]);
- // Check if this is the same particle stored twice on the stack
- if (mother) {
- if ((mother->GetPdgCode() == tPart->GetPdgCode() || (mother->Px() == tPart->Px()))) {
- // It is the same particle
- // Read in the original freeze-out information
- // and convert it from to [fm]
-
- // EPOS style
- // fpx = mother->Xv()*1e13*0.197327;
- // fpy = mother->Yv()*1e13*0.197327;
- // fpz = mother->Zv()*1e13*0.197327;
- // fpt = mother->T() *1e13*0.197327*0.5;
-
-
- // Therminator style
- fpx = mother->Xv()*1e13;
- fpy = mother->Yv()*1e13;
- fpz = mother->Zv()*1e13;
- // fpt = mother->T() *1e13*3e10;
-
- }
- }
+ if (aodtrack->IsPrimaryCandidate()) tracksPrim++;
+
+ if (fFilterBit && !aodtrack->TestFilterBit(fFilterBit)) {
+ //delete trackCopy;
+ continue;
+ }
+
+ if (fFilterMask && !aodtrack->TestFilterBit(fFilterMask)) {
+ //delete trackCopy;
+ continue;
+ }
+
+ // cout << "Muon? " << aodtrack->IsMuonTrack() << endl;
+ // cout << "Type? " << aodtrack->GetType() << endl;
+
+ // if (aodtrack->IsMuonTrack()) {
+ // cout << "muon" << endl;
+ // delete trackCopy;
+ // continue;
+ // }
+
+ //counting particles to set multiplicity
+ if(fEstEventMult==kGlobalCount){
+ AliAODTrack* trk_clone = (AliAODTrack*)aodtrack->Clone("trk_clone"); //no DCA cut for global count
+ //if (aodtrack->IsPrimaryCandidate()) //? instead of kinks?
+ if (aodtrack->Chi2perNDF() < 4.0)
+ if (aodtrack->Pt() > 0.15 && aodtrack->Pt() < 20)
+ if (aodtrack->GetTPCNcls() > 70)
+ if (aodtrack->Eta() < 0.8)
+ tNormMult++;
+ delete trk_clone;
+ }
+
+ trackCopy = CopyAODtoFemtoTrack(aodtrack);
+
+ // copying PID information from the correspondent track
+ // const AliAODTrack *aodtrackpid = fEvent->GetTrack(labels[-1-fEvent->GetTrack(i)->GetID()]);
+
+
+ AliAODTrack *aodtrackpid;
+ if((fFilterBit == (1 << (7))) || fFilterMask == 128) {//for TPC Only tracks we have to copy PID information from corresponding global tracks
+ aodtrackpid = dynamic_cast<AliAODTrack*>(fEvent->GetTrack(labels[-1-fEvent->GetTrack(i)->GetID()]));
+ }
+ else {
+ aodtrackpid = dynamic_cast<AliAODTrack*>(fEvent->GetTrack(i));
+ }
+ assert(aodtrackpid&&"Not a standard AOD");
+
+ CopyPIDtoFemtoTrack(aodtrackpid, trackCopy);
+
+ if (mcP) {
+ // Fill the hidden information with the simulated data
+ // Int_t pLabel = aodtrack->GetLabel();
+ // AliAODMCParticle *tPart = GetParticleWithLabel(mcP, (TMath::Abs(aodtrack->GetLabel())));
+ AliAODMCParticle *tPart;
+ if(aodtrack->GetLabel() > -1 ) {
+ tPart = (AliAODMCParticle*)mcP->At(aodtrack->GetLabel());
+ }
+ else {
+ tPart = NULL;
+ }
+ AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
+ double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
+ if (!tPart) {
+ fpx = fV1[0];
+ fpy = fV1[1];
+ fpz = fV1[2];
+ tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
+ tInfo->SetPDGPid(0);
+ tInfo->SetTrueMomentum(0.0, 0.0, 0.0);
+ tInfo->SetEmissionPoint(0.0, 0.0, 0.0, 0.0);
+ tInfo->SetMass(0);
+ }
+ else {
+ // Check the mother information
+
+ // Using the new way of storing the freeze-out information
+ // Final state particle is stored twice on the stack
+ // one copy (mother) is stored with original freeze-out information
+ // and is not tracked
+ // the other one (daughter) is stored with primary vertex position
+ // and is tracked
+
+ // Freeze-out coordinates
+ fpx = tPart->Xv() - fV1[0];
+ fpy = tPart->Yv() - fV1[1];
+ fpz = tPart->Zv() - fV1[2];
+ // fpt = tPart->T();
+
+ tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
+
+
+ fpx *= 1e13;
+ fpy *= 1e13;
+ fpz *= 1e13;
+ // fpt *= 1e13;
+
+ // cout << "Looking for mother ids " << endl;
+
+ //if (motherids[TMath::Abs(aodtrack->GetLabel())]>0) {
+ if(tPart->GetMother() > -1) { //MC particle has a mother
+ // cout << "Got mother id" << endl;
+ // AliAODMCParticle *mother = GetParticleWithLabel(mcP, motherids[TMath::Abs(aodtrack->GetLabel())]);
+ AliAODMCParticle *mother = (AliAODMCParticle*)mcP->At(tPart->GetMother());
+ // Check if this is the same particle stored twice on the stack
+ if (mother) {
+ if ((mother->GetPdgCode() == tPart->GetPdgCode() || (mother->Px() == tPart->Px()))) {
+ // It is the same particle
+ // Read in the original freeze-out information
+ // and convert it from to [fm]
+
+ // EPOS style
+ // fpx = mother->Xv()*1e13*0.197327;
+ // fpy = mother->Yv()*1e13*0.197327;
+ // fpz = mother->Zv()*1e13*0.197327;
+ // fpt = mother->T() *1e13*0.197327*0.5;
+
+
+ // Therminator style
+ fpx = mother->Xv()*1e13;
+ fpy = mother->Yv()*1e13;
+ fpz = mother->Zv()*1e13;
+ // fpt = mother->T() *1e13*3e10;
+
+ }
+ else { //particle's mother exists and the information about it can be added to hiddeninfo:
+ tInfo->SetMotherPdgCode(mother->GetPdgCode());
}
-
- // if (fRotateToEventPlane) {
- // double tPhi = TMath::ATan2(fpy, fpx);
- // double tRad = TMath::Hypot(fpx, fpy);
-
- // fpx = tRad*TMath::Cos(tPhi - tReactionPlane);
- // fpy = tRad*TMath::Sin(tPhi - tReactionPlane);
- // }
-
- tInfo->SetPDGPid(tPart->GetPdgCode());
-
- // if (fRotateToEventPlane) {
- // double tPhi = TMath::ATan2(tPart->Py(), tPart->Px());
- // double tRad = TMath::Hypot(tPart->Px(), tPart->Py());
-
- // tInfo->SetTrueMomentum(tRad*TMath::Cos(tPhi - tReactionPlane),
- // tRad*TMath::Sin(tPhi - tReactionPlane),
- // tPart->Pz());
- // }
- // else
- tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
- Double_t mass2 = (tPart->E() *tPart->E() -
- tPart->Px()*tPart->Px() -
- tPart->Py()*tPart->Py() -
- tPart->Pz()*tPart->Pz());
- if (mass2>0.0)
- tInfo->SetMass(TMath::Sqrt(mass2));
- else
- tInfo->SetMass(0.0);
-
- tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
- }
- trackCopy->SetHiddenInfo(tInfo);
- }
+ }
+ }
+
+ // if (fRotateToEventPlane) {
+ // double tPhi = TMath::ATan2(fpy, fpx);
+ // double tRad = TMath::Hypot(fpx, fpy);
+
+ // fpx = tRad*TMath::Cos(tPhi - tReactionPlane);
+ // fpy = tRad*TMath::Sin(tPhi - tReactionPlane);
+ // }
+
+ tInfo->SetPDGPid(tPart->GetPdgCode());
+
+ // if (fRotateToEventPlane) {
+ // double tPhi = TMath::ATan2(tPart->Py(), tPart->Px());
+ // double tRad = TMath::Hypot(tPart->Px(), tPart->Py());
+
+ // tInfo->SetTrueMomentum(tRad*TMath::Cos(tPhi - tReactionPlane),
+ // tRad*TMath::Sin(tPhi - tReactionPlane),
+ // tPart->Pz());
+ // }
+ // else
+ tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
+ Double_t mass2 = (tPart->E() *tPart->E() -
+ tPart->Px()*tPart->Px() -
+ tPart->Py()*tPart->Py() -
+ tPart->Pz()*tPart->Pz());
+ if (mass2>0.0)
+ tInfo->SetMass(TMath::Sqrt(mass2));
+ else
+ tInfo->SetMass(0.0);
+
+ tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
+
+ // // fillDCA
+ // //if (TMath::Abs(impact[0]) > 0.001) {
+ // if (tPart->IsPhysicalPrimary()){
+ // tInfo->SetPartOrigin(0);
+ // // trackCopy->SetImpactDprim(impact[0]);
+ // //cout << "Read prim" << endl;
+ // }
+ // else if (tPart->IsSecondaryFromWeakDecay()) {
+ // tInfo->SetPartOrigin(1);
+ // // trackCopy->SetImpactDweak(impact[0]);
+ // //cout << "Read wea" << endl;
+ // }
+ // else if (tPart->IsSecondaryFromMaterial()) {
+ // tInfo->SetPartOrigin(2);
+ // // trackCopy->SetImpactDmat(impact[0]);
+ // //cout << "Read mat" << endl;
+ // }
+ // //}
+ // // end fillDCA
- double pxyz[3];
- //AliExternalTrackParam *param = new AliExternalTrackParam(*aodtrack->GetInnerParam());
- trackCopy->SetInnerMomentum(aodtrack->GetTPCmomentum());
+ }
+ trackCopy->SetHiddenInfo(tInfo);
+ }
- aodtrack->PxPyPz(pxyz);//reading noconstarined momentum
- const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
- // Check the sanity of the tracks - reject zero momentum tracks
- if (ktP.Mag() == 0) {
- delete trackCopy;
- continue;
- }
- // }
-
-
- tEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
- realnofTracks++;//real number of tracks
+ double pxyz[3];
+
+ //AliExternalTrackParam *param = new AliExternalTrackParam(*aodtrack->GetInnerParam());
+ trackCopy->SetInnerMomentum(aodtrack->GetTPCmomentum());
+
+ aodtrack->PxPyPz(pxyz);//reading noconstarined momentum
+ const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
+ // Check the sanity of the tracks - reject zero momentum tracks
+ if (ktP.Mag() == 0) {
+ delete trackCopy;
+ continue;
}
-
- tEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
+ // }
+
+ tEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
+ realnofTracks++;//real number of tracks
+ }
+
+ tEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
tEvent->SetNormalizedMult(tracksPrim);
if (cent) {
else if (fEstEventMult==kCentralityZNA) {
if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("ZNA")));
}
- else if (fEstEventMult==kCentralityZNC) {
+ else if (fEstEventMult==kCentralityZNC) {
if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("ZNC")));
}
else if (fEstEventMult==kCentralityCL1) {
else if (fEstEventMult==kCentralityNPA) {
if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("NPA")));
}
- else if (fEstEventMult==kCentralityFMD) {
+ else if (fEstEventMult==kCentralityFMD) {
if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("FMD")));
}
else if(fEstEventMult==kGlobalCount){
tEvent->SetNormalizedMult(tNormMult); //particles counted in the loop, trying to reproduce GetReferenceMultiplicity. If better (default) method appears it should be changed
}
else if(fEstEventMult==kReference)
- {
- tEvent->SetNormalizedMult(fAODheader->GetRefMultiplicity());
- }
+ {
+ tEvent->SetNormalizedMult(fAODheader->GetRefMultiplicity());
+ }
else if(fEstEventMult==kTPCOnlyRef)
- {
- tEvent->SetNormalizedMult(fAODheader->GetTPConlyRefMultiplicity());
- }
+ {
+ tEvent->SetNormalizedMult(fAODheader->GetTPConlyRefMultiplicity());
+ }
else if(fEstEventMult == kVZERO)
- {
- Float_t multV0 = 0;
- for (Int_t i=0; i<64; i++)
- multV0 += fEvent->GetVZEROData()->GetMultiplicity(i);
- tEvent->SetNormalizedMult(multV0);
- }
+ {
+ Float_t multV0 = 0;
+ for (Int_t i=0; i<64; i++)
+ multV0 += fEvent->GetVZEROData()->GetMultiplicity(i);
+ tEvent->SetNormalizedMult(multV0);
+ }
- if (mcP) delete [] motherids;
+ // if (mcP) delete [] motherids;
- // cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
+ // cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
if(fReadV0)
- {
- int count_pass = 0;
- for (Int_t i = 0; i < fEvent->GetNumberOfV0s(); i++) {
- AliAODv0* aodv0 = fEvent->GetV0(i);
- if (!aodv0) continue;
- if(aodv0->GetNDaughters()>2) continue;
- if(aodv0->GetNProngs()>2) continue;
- if(aodv0->GetCharge()!=0) continue;
- if(aodv0->ChargeProng(0)==aodv0->ChargeProng(1)) continue;
- if(aodv0->CosPointingAngle(fV1)<0.998) continue;
- AliFemtoV0* trackCopyV0 = new AliFemtoV0();
- count_pass++;
- CopyAODtoFemtoV0(aodv0, trackCopyV0);
- tEvent->V0Collection()->push_back(trackCopyV0);
- //cout<<"Pushback v0 to v0collection"<<endl;
+ {
+ int count_pass = 0;
+ for (Int_t i = 0; i < fEvent->GetNumberOfV0s(); i++) {
+ AliAODv0* aodv0 = fEvent->GetV0(i);
+ if (!aodv0) continue;
+ if(aodv0->GetNDaughters()>2) continue;
+ if(aodv0->GetNProngs()>2) continue;
+ if(aodv0->GetCharge()!=0) continue;
+ if(aodv0->ChargeProng(0)==aodv0->ChargeProng(1)) continue;
+ if(aodv0->CosPointingAngle(fV1)<0.998) continue;
+
+ AliAODTrack* daughterTrackPos = (AliAODTrack*)aodv0->GetDaughter(0); //getting positive daughter track
+ AliAODTrack* daughterTrackNeg = (AliAODTrack*)aodv0->GetDaughter(1); //getting negative daughter track
+ if(!daughterTrackPos) continue; //daughter tracks must exist
+ if(!daughterTrackNeg) continue;
+ if(daughterTrackNeg->Charge() == daughterTrackPos->Charge() ) continue; //and have different charge
+
+ AliFemtoV0* trackCopyV0 = CopyAODtoFemtoV0(aodv0);
+ if(mcP) {
+ daughterTrackPos->SetAODEvent(fEvent);
+ daughterTrackNeg->SetAODEvent(fEvent);
+ if(daughterTrackPos->GetLabel() > 0 && daughterTrackNeg->GetLabel() > 0 ) {
+ AliAODMCParticle* mcParticlePos = (AliAODMCParticle*)mcP->At(daughterTrackPos->GetLabel());
+ AliAODMCParticle* mcParticleNeg = (AliAODMCParticle*)mcP->At(daughterTrackNeg->GetLabel() );
+ if((mcParticlePos!=NULL) && (mcParticleNeg!=NULL)){
+ int motherOfPosID = mcParticlePos->GetMother();
+ int motherOfNegID = mcParticleNeg->GetMother();
+ if ((motherOfPosID > -1) && (motherOfPosID == motherOfNegID)){
+ AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
+ // Both daughter tracks refer to the same mother, we can continue
+ AliAODMCParticle *v0 = (AliAODMCParticle*)mcP->At(motherOfPosID); //our V0 particle
+
+ tInfo->SetPDGPid(v0->GetPdgCode());
+ int v0MotherId = v0->GetMother();
+ if(v0MotherId>-1) { //V0 particle has a mother
+ AliAODMCParticle* motherOfV0 = (AliAODMCParticle*)mcP->At(v0MotherId);
+ tInfo->SetMotherPdgCode(motherOfV0->GetPdgCode());
+ }
+ trackCopyV0->SetHiddenInfo(tInfo);
+ }
+ }
+ }
}
+ tEvent->V0Collection()->push_back(trackCopyV0);
+ count_pass++;
+ //cout<<"Pushback v0 to v0collection"<<endl;
}
+ }
+ return tEvent;
}
-void AliFemtoEventReaderAOD::CopyAODtoFemtoTrack(AliAODTrack *tAodTrack,
- AliFemtoTrack *tFemtoTrack
- // AliPWG2AODTrack *tPWG2AODTrack
- )
+AliFemtoTrack* AliFemtoEventReaderAOD::CopyAODtoFemtoTrack(AliAODTrack *tAodTrack
+ // AliPWG2AODTrack *tPWG2AODTrack
+ )
{
// Copy the track information from the AOD into the internal AliFemtoTrack
// If it exists, use the additional information from the PWG2 AOD
+ AliFemtoTrack *tFemtoTrack = new AliFemtoTrack();
// Primary Vertex position
-
+
fEvent->GetPrimaryVertex()->GetPosition(fV1);
// fEvent->GetPrimaryVertex()->GetXYZ(fV1);
+ tFemtoTrack->SetPrimaryVertex(fV1);
tFemtoTrack->SetCharge(tAodTrack->Charge());
-
+
double pxyz[3];
tAodTrack->PxPyPz(pxyz);//reading noconstrained momentum
+
AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
tFemtoTrack->SetP(v);//setting momentum
tFemtoTrack->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
- //setting track helix
+ //setting track helix
const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
- AliFmPhysicalHelixD helix(ktP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(tFemtoTrack->Charge()));
+ AliFmPhysicalHelixD helix(ktP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(tFemtoTrack->Charge()));
tFemtoTrack->SetHelix(helix);
-
+
// Flags
tFemtoTrack->SetTrackId(tAodTrack->GetID());
tFemtoTrack->SetFlags(tAodTrack->GetFlags());
tFemtoTrack->SetLabel(tAodTrack->GetLabel());
-
- // Track quality information
+
+ // Track quality information
float covmat[6];
- tAodTrack->GetCovMatrix(covmat);
+ tAodTrack->GetCovMatrix(covmat);
+
+ if (!fDCAglobalTrack) {
+
+ tFemtoTrack->SetImpactD(tAodTrack->DCA());
+ tFemtoTrack->SetImpactZ(tAodTrack->ZAtDCA());
- // ! DCA information is done in CopyPIDtoFemtoTrack()
-
- // double impact[2];
- // double covimpact[3];
- // if (!tAodTrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
- // //cout << "sth went wrong with dca propagation" << endl;
- // tFemtoTrack->SetImpactD(-1000.0);
- // tFemtoTrack->SetImpactZ(-1000.0);
+ // double impact[2];
+ // double covimpact[3];
- // }
- // else {
- // tFemtoTrack->SetImpactD(impact[0]);
- // tFemtoTrack->SetImpactZ(impact[1]+fV1[2]);
- // }
+ // AliAODTrack* trk_clone = (AliAODTrack*)tAodTrack->Clone("trk_clone");
+ // // if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
+ // // if (!tAodTrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
+ // if (!trk_clone->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
+ // //cout << "sth went wrong with dca propagation" << endl;
+ // tFemtoTrack->SetImpactD(-1000.0);
+ // tFemtoTrack->SetImpactZ(-1000.0);
+
+ // }
+ // else {
+ // tFemtoTrack->SetImpactD(impact[0]);
+ // tFemtoTrack->SetImpactZ(impact[1]);
+ // }
+ // delete trk_clone;
+
+ }
// if (TMath::Abs(tAodTrack->Xv()) > 0.00000000001)
// tFemtoTrack->SetImpactD(TMath::Hypot(tAodTrack->Xv(), tAodTrack->Yv())*(tAodTrack->Xv()/TMath::Abs(tAodTrack->Xv())));
// else
// tFemtoTrack->SetImpactD(0.0);
// tFemtoTrack->SetImpactD(tAodTrack->DCA());
-
+
// tFemtoTrack->SetImpactZ(tAodTrack->ZAtDCA());
// tFemtoTrack->SetImpactZ(tAodTrack->Zv() - fV1[2]);
- // cout
- // << "dca" << TMath::Hypot(tAodTrack->Xv() - fV1[0], tAodTrack->Yv() - fV1[1])
- // << "xv - fv10 = "<< tAodTrack->Xv() - fV1[0]
- // << tAodTrack->Yv() - fV1[1]
-// << "xv = " << tAodTrack->Xv() << endl
-// << "fv1[0] = " << fV1[0] << endl
-// << "yv = " << tAodTrack->Yv() << endl
-// << "fv1[1] = " << fV1[1] << endl
-// << "zv = " << tAodTrack->Zv() << endl
-// << "fv1[2] = " << fV1[2] << endl
-// << "impact[0] = " << impact[0] << endl
-// << "impact[1] = " << impact[1] << endl
+ // cout
+ // << "dca" << TMath::Hypot(tAodTrack->Xv() - fV1[0], tAodTrack->Yv() - fV1[1])
+ // << "xv - fv10 = "<< tAodTrack->Xv() - fV1[0]
+ // << tAodTrack->Yv() - fV1[1]
+// << "xv = " << tAodTrack->Xv() << endl
+// << "fv1[0] = " << fV1[0] << endl
+// << "yv = " << tAodTrack->Yv() << endl
+// << "fv1[1] = " << fV1[1] << endl
+// << "zv = " << tAodTrack->Zv() << endl
+// << "fv1[2] = " << fV1[2] << endl
+// << "impact[0] = " << impact[0] << endl
+// << "impact[1] = " << impact[1] << endl
// << endl << endl ;
tFemtoTrack->SetCdd(covmat[0]);
tFemtoTrack->SetTPCchi2(tAodTrack->Chi2perNDF());
tFemtoTrack->SetTPCncls(tAodTrack->GetTPCNcls());
tFemtoTrack->SetTPCnclsF(tAodTrack->GetTPCNcls());
- tFemtoTrack->SetTPCsignalN(1);
- tFemtoTrack->SetTPCsignalS(1);
+ tFemtoTrack->SetTPCsignalN(1);
+ tFemtoTrack->SetTPCsignalS(1);
tFemtoTrack->SetTPCsignal(tAodTrack->GetTPCsignal());
// if (tPWG2AODTrack) {
// // Copy the PWG2 specific information if it exists
// tFemtoTrack->SetTPCClusterMap(tPWG2AODTrack->GetTPCClusterMap());
// tFemtoTrack->SetTPCSharedMap(tPWG2AODTrack->GetTPCSharedMap());
-
+
// double xtpc[3] = {0,0,0};
// tPWG2AODTrack->GetTPCNominalEntrancePoint(xtpc);
// tFemtoTrack->SetNominalTPCEntrancePoint(xtpc);
// tFemtoTrack->SetNominalTPCExitPoint(xtpc);
// }
// else {
- // If not use dummy values
+ // If not use dummy values
tFemtoTrack->SetTPCClusterMap(tAodTrack->GetTPCClusterMap());
tFemtoTrack->SetTPCSharedMap(tAodTrack->GetTPCSharedMap());
-
float globalPositionsAtRadii[9][3];
float bfield = 5*fMagFieldSign;
+
GetGlobalPositionAtGlobalRadiiThroughTPC(tAodTrack,bfield,globalPositionsAtRadii);
double tpcEntrance[3]={globalPositionsAtRadii[0][0],globalPositionsAtRadii[0][1],globalPositionsAtRadii[0][2]};
double **tpcPositions;
tpcPositions[i] = new double[3];
double tpcExit[3]={globalPositionsAtRadii[8][0],globalPositionsAtRadii[8][1],globalPositionsAtRadii[8][2]};
for(int i=0;i<9;i++)
- {
- tpcPositions[i][0] = globalPositionsAtRadii[i][0];
- tpcPositions[i][1] = globalPositionsAtRadii[i][1];
- tpcPositions[i][2] = globalPositionsAtRadii[i][2];
- }
+ {
+ tpcPositions[i][0] = globalPositionsAtRadii[i][0];
+ tpcPositions[i][1] = globalPositionsAtRadii[i][1];
+ tpcPositions[i][2] = globalPositionsAtRadii[i][2];
+ }
+
tFemtoTrack->SetNominalTPCEntrancePoint(tpcEntrance);
tFemtoTrack->SetNominalTPCPoints(tpcPositions);
tFemtoTrack->SetNominalTPCExitPoint(tpcExit);
-
+ for(int i=0;i<9;i++)
+ delete [] tpcPositions[i];
+ delete [] tpcPositions;
// }
-
+
// // cout << "Track has " << TMath::Hypot(tAodTrack->Xv(), tAodTrack->Yv()) << " " << tAodTrack->Zv() << " " << tAodTrack->GetTPCNcls() << endl;
-
-
+
+
int indexes[3];
for (int ik=0; ik<3; ik++) {
indexes[ik] = 0;
tFemtoTrack->SetITSHitOnLayer(ii,tAodTrack->HasPointOnITSLayer(ii));
}
-
+ return tFemtoTrack;
}
-void AliFemtoEventReaderAOD::CopyAODtoFemtoV0(AliAODv0 *tAODv0, AliFemtoV0 *tFemtoV0)
+AliFemtoV0* AliFemtoEventReaderAOD::CopyAODtoFemtoV0(AliAODv0 *tAODv0 )
{
+ AliFemtoV0 *tFemtoV0 = new AliFemtoV0();
+
tFemtoV0->SetdecayLengthV0(tAODv0->DecayLength(fV1));
tFemtoV0->SetdecayVertexV0X(tAODv0->DecayVertexV0X());
tFemtoV0->SetdecayVertexV0Y(tAODv0->DecayVertexV0Y());
tFemtoV0->SetmomNeg(momneg);
//jest cos takiego w AliFemtoV0.h czego nie ma w AliAODv0.h
- //void SettpcHitsPos(const int& i);
- //void SettpcHitsNeg(const int& i);
+ //void SettpcHitsPos(const int& i);
+ //void SettpcHitsNeg(const int& i);
//void SetTrackTopologyMapPos(unsigned int word, const unsigned long& m);
//void SetTrackTopologyMapNeg(unsigned int word, const unsigned long& m);
+
tFemtoV0->SetmomV0X(tAODv0->MomV0X());
tFemtoV0->SetmomV0Y(tAODv0->MomV0Y());
tFemtoV0->SetmomV0Z(tAODv0->MomV0Z());
tFemtoV0->SetmassK0Short(tAODv0->MassK0Short());
tFemtoV0->SetrapLambda(tAODv0->RapLambda());
tFemtoV0->SetrapK0Short(tAODv0->RapK0Short());
-
- //void SetcTauLambda( float x);
- //void SetcTauK0Short( float x);
-
+
+ //void SetcTauLambda( float x);
+ //void SetcTauK0Short( float x);
+
//tFemtoV0->SetptV0(::sqrt(tAODv0->Pt2V0())); //!
tFemtoV0->SetptV0(tAODv0->Pt());
tFemtoV0->SetptotV0(::sqrt(tAODv0->Ptot2V0()));
//tFemtoV0->SetptotPos(::sqrt(tAODv0->Ptot2Pos()));
//tFemtoV0->SetptNeg(::sqrt(tAODv0->MomNegX()*tAODv0->MomNegX()+tAODv0->MomNegY()*tAODv0->MomNegY()));
//tFemtoV0->SetptotNeg(::sqrt(tAODv0->Ptot2Neg()));
-
+
tFemtoV0->SetidNeg(tAODv0->GetNegID());
//cout<<"tAODv0->GetNegID(): "<<tAODv0->GetNegID()<<endl;
//cout<<"tFemtoV0->IdNeg(): "<<tFemtoV0->IdNeg()<<endl;
tFemtoV0->SetCosPointingAngle(tAODv0->CosPointingAngle(fV1));
//tFemtoV0->SetYV0(tAODv0->Y());
+
//void SetdedxNeg(float x);
//void SeterrdedxNeg(float x);//Gael 04Fev2002
//void SetlendedxNeg(float x);//Gael 04Fev2002
AliAODTrack *trackpos = (AliAODTrack*)tAODv0->GetDaughter(0);
AliAODTrack *trackneg = (AliAODTrack*)tAODv0->GetDaughter(1);
+
if(trackpos && trackneg)
+ {
+ tFemtoV0->SetEtaPos(trackpos->Eta());
+ tFemtoV0->SetEtaNeg(trackneg->Eta());
+ tFemtoV0->SetptotPos(tAODv0->PProng(0));
+ tFemtoV0->SetptotNeg(tAODv0->PProng(1));
+ tFemtoV0->SetptPos(trackpos->Pt());
+ tFemtoV0->SetptNeg(trackneg->Pt());
+
+ //tFemtoV0->SetEtaPos(trackpos->Eta()); //tAODv0->PseudoRapPos()
+ //tFemtoV0->SetEtaNeg(trackneg->Eta()); //tAODv0->PseudoRapNeg()
+ tFemtoV0->SetTPCNclsPos(trackpos->GetTPCNcls());
+ tFemtoV0->SetTPCNclsNeg(trackneg->GetTPCNcls());
+ tFemtoV0->SetTPCclustersPos(trackpos->GetTPCClusterMap());
+ tFemtoV0->SetTPCclustersNeg(trackneg->GetTPCClusterMap());
+ tFemtoV0->SetTPCsharingPos(trackpos->GetTPCSharedMap());
+ tFemtoV0->SetTPCsharingNeg(trackneg->GetTPCSharedMap());
+ tFemtoV0->SetNdofPos(trackpos->Chi2perNDF());
+ tFemtoV0->SetNdofNeg(trackneg->Chi2perNDF());
+ tFemtoV0->SetStatusPos(trackpos->GetStatus());
+ tFemtoV0->SetStatusNeg(trackneg->GetStatus());
+
+ tFemtoV0->SetPosNSigmaTPCK(fAODpidUtil->NumberOfSigmasTPC(trackpos,AliPID::kKaon));
+ tFemtoV0->SetNegNSigmaTPCK(fAODpidUtil->NumberOfSigmasTPC(trackneg,AliPID::kKaon));
+ tFemtoV0->SetPosNSigmaTPCP(fAODpidUtil->NumberOfSigmasTPC(trackpos,AliPID::kProton));
+ tFemtoV0->SetNegNSigmaTPCP(fAODpidUtil->NumberOfSigmasTPC(trackneg,AliPID::kProton));
+ tFemtoV0->SetPosNSigmaTPCPi(fAODpidUtil->NumberOfSigmasTPC(trackpos,AliPID::kPion));
+ tFemtoV0->SetNegNSigmaTPCPi(fAODpidUtil->NumberOfSigmasTPC(trackneg,AliPID::kPion));
+
+
+ float bfield = 5*fMagFieldSign;
+ float globalPositionsAtRadiiPos[9][3];
+ GetGlobalPositionAtGlobalRadiiThroughTPC(trackpos,bfield,globalPositionsAtRadiiPos);
+ double tpcEntrancePos[3]={globalPositionsAtRadiiPos[0][0],globalPositionsAtRadiiPos[0][1],globalPositionsAtRadiiPos[0][2]};
+ double tpcExitPos[3]={globalPositionsAtRadiiPos[8][0],globalPositionsAtRadiiPos[8][1],globalPositionsAtRadiiPos[8][2]};
+
+ float globalPositionsAtRadiiNeg[9][3];
+ GetGlobalPositionAtGlobalRadiiThroughTPC(trackneg,bfield,globalPositionsAtRadiiNeg);
+ double tpcEntranceNeg[3]={globalPositionsAtRadiiNeg[0][0],globalPositionsAtRadiiNeg[0][1],globalPositionsAtRadiiNeg[0][2]};
+ double tpcExitNeg[3]={globalPositionsAtRadiiNeg[8][0],globalPositionsAtRadiiNeg[8][1],globalPositionsAtRadiiNeg[8][2]};
+
+ AliFemtoThreeVector tmpVec;
+ tmpVec.SetX(tpcEntrancePos[0]); tmpVec.SetY(tpcEntrancePos[1]); tmpVec.SetZ(tpcEntrancePos[2]);
+ tFemtoV0->SetNominalTpcEntrancePointPos(tmpVec);
+
+ tmpVec.SetX(tpcExitPos[0]); tmpVec.SetY(tpcExitPos[1]); tmpVec.SetZ(tpcExitPos[2]);
+ tFemtoV0->SetNominalTpcExitPointPos(tmpVec);
+
+ tmpVec.SetX(tpcEntranceNeg[0]); tmpVec.SetY(tpcEntranceNeg[1]); tmpVec.SetZ(tpcEntranceNeg[2]);
+ tFemtoV0->SetNominalTpcEntrancePointNeg(tmpVec);
+
+ tmpVec.SetX(tpcExitNeg[0]); tmpVec.SetY(tpcExitNeg[1]); tmpVec.SetZ(tpcExitNeg[2]);
+ tFemtoV0->SetNominalTpcExitPointNeg(tmpVec);
+
+
+ AliFemtoThreeVector vecTpcPos[9];
+ AliFemtoThreeVector vecTpcNeg[9];
+ for(int i=0;i<9;i++)
{
- tFemtoV0->SetEtaPos(trackpos->Eta());
- tFemtoV0->SetEtaNeg(trackneg->Eta());
- tFemtoV0->SetptotPos(tAODv0->PProng(0));
- tFemtoV0->SetptotNeg(tAODv0->PProng(1));
- tFemtoV0->SetptPos(trackpos->Pt());
- tFemtoV0->SetptNeg(trackneg->Pt());
-
- //tFemtoV0->SetEtaPos(trackpos->Eta()); //tAODv0->PseudoRapPos()
- //tFemtoV0->SetEtaNeg(trackneg->Eta()); //tAODv0->PseudoRapNeg()
- tFemtoV0->SetTPCNclsPos(trackpos->GetTPCNcls());
- tFemtoV0->SetTPCNclsNeg(trackneg->GetTPCNcls());
- tFemtoV0->SetTPCclustersPos(trackpos->GetTPCClusterMap());
- tFemtoV0->SetTPCclustersNeg(trackneg->GetTPCClusterMap());
- tFemtoV0->SetTPCsharingPos(trackpos->GetTPCSharedMap());
- tFemtoV0->SetTPCsharingNeg(trackneg->GetTPCSharedMap());
- tFemtoV0->SetNdofPos(trackpos->Chi2perNDF());
- tFemtoV0->SetNdofNeg(trackneg->Chi2perNDF());
- tFemtoV0->SetStatusPos(trackpos->GetStatus());
- tFemtoV0->SetStatusNeg(trackneg->GetStatus());
-
- tFemtoV0->SetPosNSigmaTPCK(fAODpidUtil->NumberOfSigmasTPC(trackpos,AliPID::kKaon));
- tFemtoV0->SetNegNSigmaTPCK(fAODpidUtil->NumberOfSigmasTPC(trackneg,AliPID::kKaon));
- tFemtoV0->SetPosNSigmaTPCP(fAODpidUtil->NumberOfSigmasTPC(trackpos,AliPID::kProton));
- tFemtoV0->SetNegNSigmaTPCP(fAODpidUtil->NumberOfSigmasTPC(trackneg,AliPID::kProton));
- tFemtoV0->SetPosNSigmaTPCPi(fAODpidUtil->NumberOfSigmasTPC(trackpos,AliPID::kPion));
- tFemtoV0->SetNegNSigmaTPCPi(fAODpidUtil->NumberOfSigmasTPC(trackneg,AliPID::kPion));
-
-
- float bfield = 5*fMagFieldSign;
- float globalPositionsAtRadiiPos[9][3];
- GetGlobalPositionAtGlobalRadiiThroughTPC(trackpos,bfield,globalPositionsAtRadiiPos);
- double tpcEntrancePos[3]={globalPositionsAtRadiiPos[0][0],globalPositionsAtRadiiPos[0][1],globalPositionsAtRadiiPos[0][2]};
- double tpcExitPos[3]={globalPositionsAtRadiiPos[8][0],globalPositionsAtRadiiPos[8][1],globalPositionsAtRadiiPos[8][2]};
-
- float globalPositionsAtRadiiNeg[9][3];
- GetGlobalPositionAtGlobalRadiiThroughTPC(trackneg,bfield,globalPositionsAtRadiiNeg);
- double tpcEntranceNeg[3]={globalPositionsAtRadiiNeg[0][0],globalPositionsAtRadiiNeg[0][1],globalPositionsAtRadiiNeg[0][2]};
- double tpcExitNeg[3]={globalPositionsAtRadiiNeg[8][0],globalPositionsAtRadiiNeg[8][1],globalPositionsAtRadiiNeg[8][2]};
-
- AliFemtoThreeVector tmpVec;
- tmpVec.SetX(tpcEntrancePos[0]); tmpVec.SetX(tpcEntrancePos[1]); tmpVec.SetX(tpcEntrancePos[2]);
- tFemtoV0->SetNominalTpcEntrancePointPos(tmpVec);
-
- tmpVec.SetX(tpcExitPos[0]); tmpVec.SetX(tpcExitPos[1]); tmpVec.SetX(tpcExitPos[2]);
- tFemtoV0->SetNominalTpcExitPointPos(tmpVec);
-
- tmpVec.SetX(tpcEntranceNeg[0]); tmpVec.SetX(tpcEntranceNeg[1]); tmpVec.SetX(tpcEntranceNeg[2]);
- tFemtoV0->SetNominalTpcEntrancePointNeg(tmpVec);
-
- tmpVec.SetX(tpcExitNeg[0]); tmpVec.SetX(tpcExitNeg[1]); tmpVec.SetX(tpcExitNeg[2]);
- tFemtoV0->SetNominalTpcExitPointNeg(tmpVec);
-
- AliFemtoThreeVector vecTpcPos[9];
- AliFemtoThreeVector vecTpcNeg[9];
- for(int i=0;i<9;i++)
- {
- vecTpcPos[i].SetX(globalPositionsAtRadiiPos[i][0]); vecTpcPos[i].SetY(globalPositionsAtRadiiPos[i][1]); vecTpcPos[i].SetZ(globalPositionsAtRadiiPos[i][2]);
- vecTpcNeg[i].SetX(globalPositionsAtRadiiNeg[i][0]); vecTpcNeg[i].SetY(globalPositionsAtRadiiNeg[i][1]); vecTpcNeg[i].SetZ(globalPositionsAtRadiiNeg[i][2]);
- }
- tFemtoV0->SetNominalTpcPointPos(vecTpcPos);
- tFemtoV0->SetNominalTpcPointNeg(vecTpcNeg);
+ vecTpcPos[i].SetX(globalPositionsAtRadiiPos[i][0]); vecTpcPos[i].SetY(globalPositionsAtRadiiPos[i][1]); vecTpcPos[i].SetZ(globalPositionsAtRadiiPos[i][2]);
+ vecTpcNeg[i].SetX(globalPositionsAtRadiiNeg[i][0]); vecTpcNeg[i].SetY(globalPositionsAtRadiiNeg[i][1]); vecTpcNeg[i].SetZ(globalPositionsAtRadiiNeg[i][2]);
+ }
+ tFemtoV0->SetNominalTpcPointPos(vecTpcPos);
+ tFemtoV0->SetNominalTpcPointNeg(vecTpcNeg);
+
+ tFemtoV0->SetTPCMomentumPos(trackpos->GetTPCmomentum());
+ tFemtoV0->SetTPCMomentumNeg(trackneg->GetTPCmomentum());
+
+ tFemtoV0->SetdedxPos(trackpos->GetTPCsignal());
+ tFemtoV0->SetdedxNeg(trackneg->GetTPCsignal());
- tFemtoV0->SetTPCMomentumPos(trackpos->GetTPCmomentum());
- tFemtoV0->SetTPCMomentumNeg(trackneg->GetTPCmomentum());
- tFemtoV0->SetdedxPos(trackpos->GetTPCsignal());
- tFemtoV0->SetdedxNeg(trackneg->GetTPCsignal());
+ Float_t probMisPos = 1.0;
+ Float_t probMisNeg = 1.0;
- if((tFemtoV0->StatusPos()&AliESDtrack::kTOFpid)==0 || (tFemtoV0->StatusPos()&AliESDtrack::kTIME)==0 || (tFemtoV0->StatusPos()&AliESDtrack::kTOFout)==0 )
- {
- if((tFemtoV0->StatusNeg()&AliESDtrack::kTOFpid)==0 || (tFemtoV0->StatusNeg()&AliESDtrack::kTIME)==0 || (tFemtoV0->StatusNeg()&AliESDtrack::kTOFout)==0 )
+ if (tFemtoV0->StatusPos() & AliESDtrack::kTOFout & AliESDtrack::kTIME) { //AliESDtrack::kTOFpid=0x8000
+ probMisPos = fAODpidUtil->GetTOFMismatchProbability(trackpos);
+ }
+ if (tFemtoV0->StatusNeg() & AliESDtrack::kTOFout & AliESDtrack::kTIME) { //AliESDtrack::kTOFpid=0x8000
+ probMisNeg = fAODpidUtil->GetTOFMismatchProbability(trackneg);
+ }
+
+ if(// (tFemtoV0->StatusPos()& AliESDtrack::kTOFpid)==0 ||
+ (tFemtoV0->StatusPos()&AliESDtrack::kTIME)==0 || (tFemtoV0->StatusPos()&AliESDtrack::kTOFout)==0 || probMisPos > 0.01)
+ {
+ if(// (tFemtoV0->StatusNeg()&AliESDtrack::kTOFpid)==0 ||
+ (tFemtoV0->StatusNeg()&AliESDtrack::kTIME)==0 || (tFemtoV0->StatusNeg()&AliESDtrack::kTOFout)==0 || probMisNeg > 0.01)
{
tFemtoV0->SetPosNSigmaTOFK(-1000);
tFemtoV0->SetNegNSigmaTOFK(-1000);
tFemtoV0->SetTOFPionTimeNeg(-1000);
tFemtoV0->SetTOFKaonTimeNeg(-1000);
}
- }
- else
- {
- tFemtoV0->SetPosNSigmaTOFK(fAODpidUtil->NumberOfSigmasTOF(trackpos,AliPID::kKaon));
- tFemtoV0->SetNegNSigmaTOFK(fAODpidUtil->NumberOfSigmasTOF(trackneg,AliPID::kKaon));
- tFemtoV0->SetPosNSigmaTOFP(fAODpidUtil->NumberOfSigmasTOF(trackpos,AliPID::kProton));
- tFemtoV0->SetNegNSigmaTOFP(fAODpidUtil->NumberOfSigmasTOF(trackneg,AliPID::kProton));
- tFemtoV0->SetPosNSigmaTOFPi(fAODpidUtil->NumberOfSigmasTOF(trackpos,AliPID::kPion));
- tFemtoV0->SetNegNSigmaTOFPi(fAODpidUtil->NumberOfSigmasTOF(trackneg,AliPID::kPion));
-
- double TOFSignalPos = trackpos->GetTOFsignal();
- double TOFSignalNeg = trackneg->GetTOFsignal();
- double pidPos[5];
- double pidNeg[5];
- trackpos->GetIntegratedTimes(pidPos);
- trackneg->GetIntegratedTimes(pidNeg);
-
- tFemtoV0->SetTOFPionTimePos(TOFSignalPos-pidPos[2]);
- tFemtoV0->SetTOFKaonTimePos(TOFSignalPos-pidPos[3]);
- tFemtoV0->SetTOFProtonTimePos(TOFSignalPos-pidPos[4]);
- tFemtoV0->SetTOFPionTimeNeg(TOFSignalNeg-pidNeg[2]);
- tFemtoV0->SetTOFKaonTimeNeg(TOFSignalNeg-pidNeg[3]);
- tFemtoV0->SetTOFProtonTimeNeg(TOFSignalNeg-pidNeg[4]);
- }
}
- else
+ else
{
- tFemtoV0->SetStatusPos(999);
- tFemtoV0->SetStatusNeg(999);
+ if(trackpos->IsOn(AliESDtrack::kTOFout & AliESDtrack::kTIME)) {
+ tFemtoV0->SetPosNSigmaTOFK(fAODpidUtil->NumberOfSigmasTOF(trackpos,AliPID::kKaon));
+ tFemtoV0->SetPosNSigmaTOFP(fAODpidUtil->NumberOfSigmasTOF(trackpos,AliPID::kProton));
+ tFemtoV0->SetPosNSigmaTOFPi(fAODpidUtil->NumberOfSigmasTOF(trackpos,AliPID::kPion));
+ }
+ if(trackneg->IsOn(AliESDtrack::kTOFout & AliESDtrack::kTIME)) {
+ tFemtoV0->SetNegNSigmaTOFK(fAODpidUtil->NumberOfSigmasTOF(trackneg,AliPID::kKaon));
+ tFemtoV0->SetNegNSigmaTOFP(fAODpidUtil->NumberOfSigmasTOF(trackneg,AliPID::kProton));
+ tFemtoV0->SetNegNSigmaTOFPi(fAODpidUtil->NumberOfSigmasTOF(trackneg,AliPID::kPion));
+ }
+ double TOFSignalPos = trackpos->GetTOFsignal();
+ double TOFSignalNeg = trackneg->GetTOFsignal();
+ TOFSignalPos -= fAODpidUtil->GetTOFResponse().GetStartTime(trackpos->P());
+ TOFSignalNeg -= fAODpidUtil->GetTOFResponse().GetStartTime(trackneg->P());
+ double pidPos[5];
+ double pidNeg[5];
+ trackpos->GetIntegratedTimes(pidPos);
+ trackneg->GetIntegratedTimes(pidNeg);
+
+ tFemtoV0->SetTOFPionTimePos(TOFSignalPos-pidPos[2]);
+ tFemtoV0->SetTOFKaonTimePos(TOFSignalPos-pidPos[3]);
+ tFemtoV0->SetTOFProtonTimePos(TOFSignalPos-pidPos[4]);
+ tFemtoV0->SetTOFPionTimeNeg(TOFSignalNeg-pidNeg[2]);
+ tFemtoV0->SetTOFKaonTimeNeg(TOFSignalNeg-pidNeg[3]);
+ tFemtoV0->SetTOFProtonTimeNeg(TOFSignalNeg-pidNeg[4]);
}
+
+ }
+ else
+ {
+
+ tFemtoV0->SetStatusPos(999);
+ tFemtoV0->SetStatusNeg(999);
+ }
+
tFemtoV0->SetOnFlyStatusV0(tAODv0->GetOnFlyStatus());
+ return tFemtoV0;
}
void AliFemtoEventReaderAOD::SetFilterBit(UInt_t ibit)
}
while (posstack < mcP->GetEntries());
}
-
+
return 0;
}
-void AliFemtoEventReaderAOD::CopyPIDtoFemtoTrack(AliAODTrack *tAodTrack,
- AliFemtoTrack *tFemtoTrack)
+void AliFemtoEventReaderAOD::CopyPIDtoFemtoTrack(AliAODTrack *tAodTrack,
+ AliFemtoTrack *tFemtoTrack)
{
- // copying DCA information (taking it from global tracks gives better resolution than from TPC-only)
+ if (fDCAglobalTrack) {
- double impact[2];
- double covimpact[3];
+ // code from Michael and Prabhat from AliAnalysisTaskDptDptCorrelations
+ const AliAODVertex* vertex = (AliAODVertex*) fEvent->GetPrimaryVertex();
+ float vertexX = -999.;
+ float vertexY = -999.;
+ float vertexZ = -999.;
- if (!tAodTrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
- //cout << "sth went wrong with dca propagation" << endl;
- tFemtoTrack->SetImpactD(-1000.0);
- tFemtoTrack->SetImpactZ(-1000.0);
+ if(vertex) {
+ Double32_t fCov[6];
+ vertex->GetCovarianceMatrix(fCov);
+ if(vertex->GetNContributors() > 0) {
+ if(fCov[5] != 0) {
+ vertexX = vertex->GetX();
+ vertexY = vertex->GetY();
+ vertexZ = vertex->GetZ();
- }
- else {
- tFemtoTrack->SetImpactD(impact[0]);
- tFemtoTrack->SetImpactZ(impact[1]);
- }
+ }
+ }
+ }
+
+ Double_t pos[3];
+ tAodTrack->GetXYZ(pos);
+
+ Double_t DCAX = pos[0] - vertexX;
+ Double_t DCAY = pos[1] - vertexY;
+ Double_t DCAZ = pos[2] - vertexZ;
+
+ Double_t DCAXY = TMath::Sqrt((DCAX*DCAX) + (DCAY*DCAY));
+
+ tFemtoTrack->SetImpactD(DCAXY);
+ tFemtoTrack->SetImpactZ(DCAZ);
+
+
+ // // // copying DCA information (taking it from global tracks gives better resolution than from TPC-only)
+
+ // double impact[2];
+ // double covimpact[3];
+
+ // AliAODTrack* trk_clone = (AliAODTrack*)tAodTrack->Clone("trk_clone");
+ // // if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
+ // // if (!tAodTrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
+ // if (!trk_clone->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
+
+ // // if (!tAodTrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
+ // //cout << "sth went wrong with dca propagation" << endl;
+ // tFemtoTrack->SetImpactD(-1000.0);
+ // tFemtoTrack->SetImpactZ(-1000.0);
+
+ // }
+ // else {
+ // tFemtoTrack->SetImpactD(impact[0]);
+ // tFemtoTrack->SetImpactZ(impact[1]);
+ // }
+ // delete trk_clone;
+ }
double aodpid[10];
tAodTrack->GetPID(aodpid);
aodpid[2] = -100000.0;
aodpid[3] = -100000.0;
aodpid[4] = -100000.0;
-
+
double tTOF = 0.0;
+ Float_t probMis = 1.0;
//what is that code? for what do we need it? nsigma values are not enough?
- if (tAodTrack->GetStatus() & AliESDtrack::kTOFpid) { //AliESDtrack::kTOFpid=0x8000
- tTOF = tAodTrack->GetTOFsignal();
- tAodTrack->GetIntegratedTimes(aodpid);
+ if (tAodTrack->GetStatus() & AliESDtrack::kTOFout & AliESDtrack::kTIME) { //AliESDtrack::kTOFpid=0x8000
+ tTOF = tAodTrack->GetTOFsignal();
+ tAodTrack->GetIntegratedTimes(aodpid);
+
+ tTOF -= fAODpidUtil->GetTOFResponse().GetStartTime(tAodTrack->P());
+
+ probMis = fAODpidUtil->GetTOFMismatchProbability(tAodTrack);
+ }
+
- tTOF -= fAODpidUtil->GetTOFResponse().GetStartTime(tAodTrack->P());
- }
- tFemtoTrack->SetTofExpectedTimes(tTOF-aodpid[2], tTOF-aodpid[3], tTOF-aodpid[4]);
-
+ tFemtoTrack->SetTofExpectedTimes(tTOF-aodpid[2], tTOF-aodpid[3], tTOF-aodpid[4]);
+
////// TPC ////////////////////////////////////////////
- float nsigmaTPCK=-1000.;
- float nsigmaTPCPi=-1000.;
- float nsigmaTPCP=-1000.;
-
+ float nsigmaTPCK=-1000.;
+ float nsigmaTPCPi=-1000.;
+ float nsigmaTPCP=-1000.;
+ float nsigmaTPCE=-1000.;
+
// cout<<"in reader fESDpid"<<fESDpid<<endl;
if (tAodTrack->IsOn(AliESDtrack::kTPCpid)){ //AliESDtrack::kTPCpid=0x0080
nsigmaTPCK = fAODpidUtil->NumberOfSigmasTPC(tAodTrack,AliPID::kKaon);
nsigmaTPCPi = fAODpidUtil->NumberOfSigmasTPC(tAodTrack,AliPID::kPion);
nsigmaTPCP = fAODpidUtil->NumberOfSigmasTPC(tAodTrack,AliPID::kProton);
+ nsigmaTPCE = fAODpidUtil->NumberOfSigmasTPC(tAodTrack,AliPID::kElectron);
}
tFemtoTrack->SetNSigmaTPCPi(nsigmaTPCPi);
tFemtoTrack->SetNSigmaTPCK(nsigmaTPCK);
tFemtoTrack->SetNSigmaTPCP(nsigmaTPCP);
+ tFemtoTrack->SetNSigmaTPCE(nsigmaTPCE);
+
+ tFemtoTrack->SetTPCchi2(tAodTrack->Chi2perNDF());
+ tFemtoTrack->SetTPCncls(tAodTrack->GetTPCNcls());
+ tFemtoTrack->SetTPCnclsF(tAodTrack->GetTPCNcls());
- tFemtoTrack->SetTPCchi2(tAodTrack->Chi2perNDF());
- tFemtoTrack->SetTPCncls(tAodTrack->GetTPCNcls());
- tFemtoTrack->SetTPCnclsF(tAodTrack->GetTPCNcls());
-
- tFemtoTrack->SetTPCsignalN(1);
- tFemtoTrack->SetTPCsignalS(1);
+ tFemtoTrack->SetTPCsignalN(1);
+ tFemtoTrack->SetTPCsignalS(1);
tFemtoTrack->SetTPCsignal(tAodTrack->GetTPCsignal());
-
- ///////TOF//////////////////////
- float vp=-1000.;
- float nsigmaTOFPi=-1000.;
- float nsigmaTOFK=-1000.;
- float nsigmaTOFP=-1000.;
+ ///////TOF//////////////////////
- if ((tAodTrack->GetStatus() & AliESDtrack::kTOFpid) && //AliESDtrack::kTOFpid=0x8000
- (tAodTrack->GetStatus() & AliESDtrack::kTOFout) && //AliESDtrack::kTOFout=0x2000
- (tAodTrack->GetStatus() & AliESDtrack::kTIME) ) //AliESDtrack::kTIME=0x80000000
- {
- if(tAodTrack->IsOn(AliESDtrack::kTOFpid)) //AliESDtrack::kTOFpid=0x8000
+ float vp=-1000.;
+ float nsigmaTOFPi=-1000.;
+ float nsigmaTOFK=-1000.;
+ float nsigmaTOFP=-1000.;
+ float nsigmaTOFE=-1000.;
+
+ if (// (tAodTrack->GetStatus() & AliESDtrack::kTOFpid) &&
+ //AliESDtrack::kTOFpid=0x8000
+ (tAodTrack->GetStatus() & AliESDtrack::kTOFout) && //AliESDtrack::kTOFout=0x2000
+ (tAodTrack->GetStatus() & AliESDtrack::kTIME) //AliESDtrack::kTIME=0x80000000
+ && probMis < 0.01) // TOF mismatch probaility
+ {
+ if(tAodTrack->IsOn(AliESDtrack::kTOFout & AliESDtrack::kTIME)) //AliESDtrack::kTOFpid=0x8000
{
nsigmaTOFPi = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kPion);
nsigmaTOFK = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kKaon);
nsigmaTOFP = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kProton);
+ nsigmaTOFE = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kElectron);
Double_t len=200;// esdtrack->GetIntegratedLength(); !!!!!
Double_t tof=tAodTrack->GetTOFsignal();
if(tof > 0.) vp=len/tof/0.03;
}
- }
- tFemtoTrack->SetVTOF(vp);
- tFemtoTrack->SetNSigmaTOFPi(nsigmaTOFPi);
- tFemtoTrack->SetNSigmaTOFK(nsigmaTOFK);
- tFemtoTrack->SetNSigmaTOFP(nsigmaTOFP);
+ }
+ tFemtoTrack->SetVTOF(vp);
+ tFemtoTrack->SetNSigmaTOFPi(nsigmaTOFPi);
+ tFemtoTrack->SetNSigmaTOFK(nsigmaTOFK);
+ tFemtoTrack->SetNSigmaTOFP(nsigmaTOFP);
+ tFemtoTrack->SetNSigmaTOFE(nsigmaTOFE);
+
-
- //////////////////////////////////////
+ //////////////////////////////////////
}
void AliFemtoEventReaderAOD::SetCentralityPreSelection(double min, double max)
{
fCentRange[0] = min; fCentRange[1] = max;
- fUsePreCent = 1;
+ fUsePreCent = 1;
fEstEventMult = kCentrality;
}
void AliFemtoEventReaderAOD::SetEPVZERO(Bool_t iepvz)
{
- fisEPVZ = iepvz;
+ fisEPVZ = iepvz;
}
void AliFemtoEventReaderAOD::GetGlobalPositionAtGlobalRadiiThroughTPC(AliAODTrack *track, Float_t bfield, Float_t globalPositionsAtRadii[9][3])
// Bigger loop has bad precision, we're nearly one centimeter too far, go back in small steps.
while (globalRadius>Rwanted[iR]){
- x-=.1;
- // printf("propagating to x %5.2f\n",x);
- if(!etp.PropagateTo(x,bfield))break;
- etp.GetXYZ(xyz); // GetXYZ returns global coordinates
- globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii
+ x-=.1;
+ // printf("propagating to x %5.2f\n",x);
+ if(!etp.PropagateTo(x,bfield))break;
+ etp.GetXYZ(xyz); // GetXYZ returns global coordinates
+ globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii
}
//printf("At Radius:%05.2f (local x %5.2f). Setting position to x %4.1f y %4.1f z %4.1f\n",globalRadius,x,xyz[0],xyz[1],xyz[2]);
globalPositionsAtRadii[iR][0]=xyz[0];
}
}
+void AliFemtoEventReaderAOD::SetpA2013(Bool_t pa2013)
+{
+ fpA2013 = pa2013;
+}
+
+void AliFemtoEventReaderAOD::SetUseMVPlpSelection(Bool_t mvplp)
+{
+ fMVPlp = mvplp;
+}
+void AliFemtoEventReaderAOD::SetIsPileUpEvent(Bool_t ispileup)
+{
+ fisPileUp = ispileup;
+}
+
+
+
+void AliFemtoEventReaderAOD::SetDCAglobalTrack(Bool_t dcagt)
+{
+ fDCAglobalTrack = dcagt;
+}
+
+
+bool AliFemtoEventReaderAOD::RejectEventCentFlat(float MagField, float CentPercent)
+{ // to flatten centrality distribution
+ bool RejectEvent = kFALSE;
+ int weightBinSign;
+ TRandom3* fRandomNumber = new TRandom3(); //for 3D, random sign switching
+ fRandomNumber->SetSeed(0);
+
+ if(MagField > 0) weightBinSign = 0;
+ else weightBinSign = 1;
+ float kCentWeight[2][9] = {{.878,.876,.860,.859,.859,.88,.873,.879,.894},
+ {.828,.793,.776,.772,.775,.796,.788,.804,.839}};
+ int weightBinCent = (int) CentPercent;
+ if(fRandomNumber->Rndm() > kCentWeight[weightBinSign][weightBinCent]) RejectEvent = kTRUE;
+ delete fRandomNumber;
+ return RejectEvent;
+}
+
+void AliFemtoEventReaderAOD::SetCentralityFlattening(Bool_t dcagt)
+{
+ fFlatCent = dcagt;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff