X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=blobdiff_plain;f=PWGCF%2FFEMTOSCOPY%2FAliFemto%2FAliFemtoEventReaderAOD.cxx;h=6a86baf443918e933b9d75308a112fb97da20e98;hp=a017878b14d4edffef0566afaf6d2b1be142d31f;hb=af468deb14520513af31ef6e597a003d930253ec;hpb=2bcd4b963b56f881d53aed71620308d6b93f9b34 diff --git a/PWGCF/FEMTOSCOPY/AliFemto/AliFemtoEventReaderAOD.cxx b/PWGCF/FEMTOSCOPY/AliFemto/AliFemtoEventReaderAOD.cxx index a017878b14d..6a86baf4439 100644 --- a/PWGCF/FEMTOSCOPY/AliFemto/AliFemtoEventReaderAOD.cxx +++ b/PWGCF/FEMTOSCOPY/AliFemto/AliFemtoEventReaderAOD.cxx @@ -11,6 +11,7 @@ #include "TFile.h" #include "TTree.h" +#include "TRandom3.h" #include "AliAODEvent.h" #include "AliAODTrack.h" #include "AliAODVertex.h" @@ -28,11 +29,14 @@ #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; @@ -40,7 +44,7 @@ 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), @@ -56,12 +60,19 @@ AliFemtoEventReaderAOD::AliFemtoEventReaderAOD(): 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); @@ -86,18 +97,24 @@ AliFemtoEventReaderAOD::AliFemtoEventReaderAOD(const AliFemtoEventReaderAOD &aRe 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(); @@ -112,6 +129,14 @@ AliFemtoEventReaderAOD::AliFemtoEventReaderAOD(const AliFemtoEventReaderAOD &aRe 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 @@ -132,10 +157,9 @@ AliFemtoEventReaderAOD& AliFemtoEventReaderAOD::operator=(const AliFemtoEventRea { // 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; @@ -154,6 +178,14 @@ AliFemtoEventReaderAOD& AliFemtoEventReaderAOD::operator=(const AliFemtoEventRea 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; } @@ -168,7 +200,7 @@ AliFemtoString AliFemtoEventReaderAOD::Report() //__________________ 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); @@ -177,26 +209,26 @@ void AliFemtoEventReaderAOD::SetInputFile(const char* 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 "<AddFile(buffer); - delete tree; - } - aodFile->Close(); + { + // cout<<"putting file "<AddFile(buffer); + delete tree; + } + aodFile->Close(); } - delete aodFile; - } + delete aodFile; } + } } AliFemtoEvent* AliFemtoEventReaderAOD::ReturnHbtEvent() @@ -205,15 +237,15 @@ AliFemtoEvent* AliFemtoEventReaderAOD::ReturnHbtEvent() // convert it to AliFemtoEvent and return // for further analysis AliFemtoEvent *hbtEvent = 0; - // cout<<"reader"<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) { @@ -221,38 +253,40 @@ AliFemtoEvent* AliFemtoEventReaderAOD::ReturnHbtEvent() // cout << "Reading only tracks with the additional information" << endl; // } - fNumberofEvent=fTree->GetEntries(); - // cout<<"Number of Entries in file "<GetEntries(); + // cout<<"Number of Entries in file "<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 @@ -264,7 +298,7 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoEvent(AliFemtoEvent *tEvent) tEvent->SetZDCParticipants(0); tEvent->SetTriggerMask(fEvent->GetTriggerMask()); tEvent->SetTriggerCluster(fEvent->GetTriggerCluster()); - + // Attempt to access MC header AliAODMCHeader *mcH; TClonesArray *mcP=0; @@ -280,70 +314,100 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoEvent(AliFemtoEvent *tEvent) } } - 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; ipGetEntries(); ip++) motherids[ip] = 0; - - // Read in mother ids - AliAODMCParticle *motherpart; - for (int ip=0; ipGetEntries(); 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(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; ipGetEntries(); 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;iGetTrack(i); + const AliAODTrack *aodtrack=dynamic_cast(fEvent->GetTrack(i)); + assert(aodtrack&&"Not a standard AOD"); if (!aodtrack->TestFilterBit(fFilterBit)) { if(aodtrack->GetID() < 0) continue; labels[aodtrack->GetID()] = i; @@ -352,46 +416,46 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoEvent(AliFemtoEvent *tEvent) int tNormMult = 0; for (int i=0;iAt(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]; @@ -406,7 +470,7 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoEvent(AliFemtoEvent *tEvent) // fpy *= 1e13; // fpz *= 1e13; // fpt *= 1e13; - + // // cout << "Looking for mother ids " << endl; // if (motherids[TMath::Abs(aodtrack->GetLabel())]>0) { // // cout << "Got mother id" << endl; @@ -416,27 +480,27 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoEvent(AliFemtoEvent *tEvent) // // 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); // // } @@ -446,7 +510,7 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoEvent(AliFemtoEvent *tEvent) // // 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()); @@ -459,9 +523,9 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoEvent(AliFemtoEvent *tEvent) // 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); @@ -477,182 +541,237 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoEvent(AliFemtoEvent *tEvent) // } // } // 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(fEvent->GetTrack(i)); + AliAODTrack *aodtrack=dynamic_cast(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(fEvent->GetTrack(labels[-1-fEvent->GetTrack(i)->GetID()])); + } + else { + aodtrackpid = dynamic_cast(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) { tEvent->SetCentralityV0(cent->GetCentralityPercentile("V0M")); + tEvent->SetCentralityV0A(cent->GetCentralityPercentile("V0A")); + tEvent->SetCentralityV0C(cent->GetCentralityPercentile("V0C")); tEvent->SetCentralityZNA(cent->GetCentralityPercentile("ZNA")); + tEvent->SetCentralityZNC(cent->GetCentralityPercentile("ZNC")); tEvent->SetCentralityCL1(cent->GetCentralityPercentile("CL1")); - // tEvent->SetCentralityFMD(cent->GetCentralityPercentile("FMD")); - // tEvent->SetCentralityTrk(cent->GetCentralityPercentile("TRK")); + tEvent->SetCentralityCL0(cent->GetCentralityPercentile("CL0")); + tEvent->SetCentralityTKL(cent->GetCentralityPercentile("TKL")); + tEvent->SetCentralityFMD(cent->GetCentralityPercentile("FMD")); + tEvent->SetCentralityFMD(cent->GetCentralityPercentile("NPA")); + // tEvent->SetCentralitySPD1(cent->GetCentralityPercentile("CL1")); + tEvent->SetCentralityTrk(cent->GetCentralityPercentile("TRK")); + tEvent->SetCentralityCND(cent->GetCentralityPercentile("CND")); } if (fEstEventMult==kCentrality) { @@ -661,113 +780,184 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoEvent(AliFemtoEvent *tEvent) if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("V0M"))); // if (cent) tEvent->SetNormalizedMult((int) cent->GetCentralityPercentile("V0M")); } + else if (fEstEventMult==kCentralityV0A) { + if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("V0A"))); + } + else if (fEstEventMult==kCentralityV0C) { + if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("V0C"))); + } else if (fEstEventMult==kCentralityZNA) { if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("ZNA"))); } + else if (fEstEventMult==kCentralityZNC) { + if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("ZNC"))); + } else if (fEstEventMult==kCentralityCL1) { if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("CL1"))); } + else if (fEstEventMult==kCentralityCL0) { + if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("CL0"))); + } + else if (fEstEventMult==kCentralityTRK) { + if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("TRK"))); + } + else if (fEstEventMult==kCentralityTKL) { + if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("TKL"))); + } + else if (fEstEventMult==kCentralityCND) { + if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("CND"))); + } + else if (fEstEventMult==kCentralityNPA) { + if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("NPA"))); + } + 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 "<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"<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"<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) { - // ! DCA information is done in CopyPIDtoFemtoTrack() - - // double impact[2]; - // double covimpact[3]; + tFemtoTrack->SetImpactD(tAodTrack->DCA()); + tFemtoTrack->SetImpactZ(tAodTrack->ZAtDCA()); - // 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]+fV1[2]); - // } + // 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)) { + // //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()); @@ -775,18 +965,18 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoTrack(AliAODTrack *tAodTrack, // 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]); @@ -797,15 +987,15 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoTrack(AliAODTrack *tAodTrack, 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); @@ -813,13 +1003,13 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoTrack(AliAODTrack *tAodTrack, // 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; @@ -828,20 +1018,23 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoTrack(AliAODTrack *tAodTrack, 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; @@ -853,11 +1046,13 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoTrack(AliAODTrack *tAodTrack, 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()); @@ -880,12 +1075,13 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoV0(AliAODv0 *tAODv0, AliFemtoV0 *tFem 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()); @@ -902,10 +1098,10 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoV0(AliAODv0 *tAODv0, AliFemtoV0 *tFem 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())); @@ -913,7 +1109,7 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoV0(AliAODv0 *tAODv0, AliFemtoV0 *tFem //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(): "<GetNegID()<IdNeg(): "<IdNeg()<SetCosPointingAngle(tAODv0->CosPointingAngle(fV1)); //tFemtoV0->SetYV0(tAODv0->Y()); + //void SetdedxNeg(float x); //void SeterrdedxNeg(float x);//Gael 04Fev2002 //void SetlendedxNeg(float x);//Gael 04Fev2002 @@ -937,79 +1134,94 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoV0(AliAODv0 *tAODv0, AliFemtoV0 *tFem 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::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 ) - { - if((tFemtoV0->StatusNeg()&AliESDtrack::kTOFpid)==0 || (tFemtoV0->StatusNeg()&AliESDtrack::kTIME)==0 || (tFemtoV0->StatusNeg()&AliESDtrack::kTOFout)==0 ) + 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); @@ -1025,37 +1237,46 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoV0(AliAODv0 *tAODv0, AliFemtoV0 *tFem 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) @@ -1112,29 +1333,70 @@ AliAODMCParticle* AliFemtoEventReaderAOD::GetParticleWithLabel(TClonesArray *mcP } 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); @@ -1149,82 +1411,95 @@ void AliFemtoEventReaderAOD::CopyPIDtoFemtoTrack(AliAODTrack *tAodTrack, 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"<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; } @@ -1258,7 +1533,7 @@ void AliFemtoEventReaderAOD::SetMagneticFieldSign(int s) void AliFemtoEventReaderAOD::SetEPVZERO(Bool_t iepvz) { - fisEPVZ = iepvz; + fisEPVZ = iepvz; } void AliFemtoEventReaderAOD::GetGlobalPositionAtGlobalRadiiThroughTPC(AliAODTrack *track, Float_t bfield, Float_t globalPositionsAtRadii[9][3]) @@ -1309,11 +1584,11 @@ void AliFemtoEventReaderAOD::GetGlobalPositionAtGlobalRadiiThroughTPC(AliAODTrac // 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]; @@ -1329,4 +1604,47 @@ void AliFemtoEventReaderAOD::GetGlobalPositionAtGlobalRadiiThroughTPC(AliAODTrac } } +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; +}