if (fIsPbPb) fHistCentralityQuality[1]->Fill(CurrentCentralityQuality);
fHistVertexZ[1]->Fill(vtxz);
- cout<<"Event Selected: "<<select<<endl;
+ //cout<<"Event Selected: "<<select<<endl;
return select;
for (Int_t iPtBin = 1; iPtBin < (GetNPtBinsPID() + 1); iPtBin++) {
aout->AddLast((TH2F*)GetHistDataTPCTOFProjection(charge, species, iPtBin));
- cout<<"aout entries: "<<aout->GetEntriesFast()<<endl;
+
//cout<<"aout entries: "<<aout->GetEntriesFast()<<endl;
}
return aout;
if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
TH3F* hout = new TH3F(Form("%s%s%s%s",name,fHistoName[histoclass].Data(),fParticleName[expspecies].Data(),fPtClassName[ptclass].Data()),
- Form("PID %s (Exp: %s);p_{T} (GeV/c);#Delta t (ms);dE/dx (a.u.)",fHistoName[histoclass].Data(),fParticleName[expspecies].Data()),
+ Form("PID %s (Exp: %s);p_{T} (GeV/c);#Delta t (ps);dE/dx (a.u.)",fHistoName[histoclass].Data(),fParticleName[expspecies].Data()),
fNPtBinsPID[ptclass],fPtBoundaryPID[ptclass],fPtBoundaryPID[ptclass+1],
fTOFbins[ptclass][expspecies],fTOFLowerBound[ptclass][expspecies],fTOFUpperBound[ptclass][expspecies],
fTPCbins[ptclass][expspecies],fTPCLowerBound[ptclass][expspecies],fTPCUpperBound[ptclass][expspecies]);
if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
TH2F* hout = new TH2F(Form("%s%s%s%s",name,fHistoName[histoclass].Data(),fParticleName[expspecies].Data(),fPtClassName[ptclass].Data()),
- Form("TOF Mismatch %s (Exp: %s);p_{T} (GeV/c);#Delta t (ms)",fHistoName[histoclass].Data(),fParticleName[expspecies].Data()),
+ Form("TOF Mismatch %s (Exp: %s);p_{T} (GeV/c);#Delta t (ps)",fHistoName[histoclass].Data(),fParticleName[expspecies].Data()),
fNPtBinsPID[ptclass],fPtBoundaryPID[ptclass],fPtBoundaryPID[ptclass+1],
fTOFbins[ptclass][expspecies],fTOFLowerBound[ptclass][expspecies],fTOFUpperBound[ptclass][expspecies]);
fAssociatedTracks(0x0),
fCurrentAODEvent(0x0),
fOutputList(0x0),
- fPtSpectrum(0x0),
- fCorrelations(0x0),
- fMixedEvents(0x0),
+ fPtSpectrumTOFbins(0x0),
+ fCorrelationsTOFbins(0x0),
+ fMixedEventsTOFbins(0x0),
+ fPtSpectrumTOFTPCbins(0x0),
+ fCorrelationsTOFTPCbins(0x0),
+ fMixedEventsTOFTPCbins(0x0),
fTOFhistos(0x0),
+ fTOFmismatch(0x0),
fTOFPtAxis(0x0),
fTOFTPChistos(0x0),
+ fTOFTPCmismatch(0x0),
fTOFTPCPtAxis(0x0),
fNDEtaBins(32),
fNDPhiBins(32),
fPoolSize(1000),
fMixEvents(kTRUE),
fMixTriggers(kFALSE),
- fCalculateTOFmismatch(kTRUE),
+ fCalculateMismatch(kTRUE),
fT0Fill(0x0),
fLvsEta(0x0),
fLvsEtaProjections(0x0),
fMakeTOFcorrelations(kTRUE),
- fMakeTOFTPCcorrelations(kFALSE),
- fDebug(0)
+ fMakeTOFTPCcorrelations(kFALSE)
+ //fDebug(0)
{
fAssociatedTracks(0x0),
fCurrentAODEvent(0x0),
fOutputList(0x0),
- fPtSpectrum(0x0),
- fCorrelations(0x0),
- fMixedEvents(0x0),
+ fPtSpectrumTOFbins(0x0),
+ fCorrelationsTOFbins(0x0),
+ fMixedEventsTOFbins(0x0),
+ fPtSpectrumTOFTPCbins(0x0),
+ fCorrelationsTOFTPCbins(0x0),
+ fMixedEventsTOFTPCbins(0x0),
fTOFhistos(0x0),
+ fTOFmismatch(0x0),
fTOFPtAxis(0x0),
fTOFTPChistos(0x0),
- fTOFTPCPtAxis(0x0),
+ fTOFTPCmismatch(0x0),
+ fTOFTPCPtAxis(0x0),
fNDEtaBins(32),
fNDPhiBins(32),
fMinNEventsForMixing(5),
fPoolSize(1000),
fMixEvents(kTRUE),
fMixTriggers(kFALSE),
- fCalculateTOFmismatch(kTRUE),
+ fCalculateMismatch(kTRUE),
fT0Fill(0x0),
fLvsEta(0x0),
fLvsEtaProjections(0x0),
fMakeTOFcorrelations(kTRUE),
- fMakeTOFTPCcorrelations(kFALSE),
- fDebug(0)
+ fMakeTOFTPCcorrelations(kFALSE)
+ //fDebug(0)
{
// Getting the pointer to the PID response object.
fPIDResponse = inputHandler->GetPIDResponse();
- // Not very neat - only set up for 0-5% analysis.
- Int_t nCentralityBins = 15;
- Double_t centralityBins[] = {0., 1., 2., 3., 4., 5., 10., 20., 30., 40., 50., 60., 70., 80., 90., 100.1 };
+ // For now we don't bin in multiplicity for pp.
+ Int_t nCentralityBins = -1;
+ Double_t* centralityBins = 0x0;
+ if (fEventCuts->GetIsPbPb()) {
+ nCentralityBins = 15;
+ Double_t tmp[] = {0., 1., 2., 3., 4., 5., 10., 20., 30., 40., 50., 60., 70., 80., 90., 100.1 };
+ centralityBins = tmp;
+ } else {
+ nCentralityBins = 1;
+ Double_t tmp[] = {0.,1.};
+ centralityBins = tmp;
+ }
Int_t nZvtxBins = 7;
Double_t vertexBins[] = {-7., -5., -3., -1., 1., 3., 5., 7.};
fTrackCutsAssociated->CreateHistos();
fOutputList->Add(fTrackCutsAssociated);
- // Get the pT axis for the TOF PID correlations.
- Double_t* ptaxis = fTrackCutsAssociated->GetPtAxisPID();
- Int_t nptbins = fTrackCutsAssociated->GetNPtBinsPID();
- fTOFPtAxis = new TAxis(nptbins, ptaxis);
- fTOFPtAxis->SetName("fTOFPtAxis");
- fTOFPtAxis->SetTitle("p_{T} GeV/c");
- fOutputList->Add(fTOFPtAxis);
-
- // Create Pt spectrum histogram.
- fPtSpectrum = new TH1F("fPtSpectrum","p_{T} Spectrum;p_{T} (GeV/c);Count",nptbins,ptaxis);
- fOutputList->Add(fPtSpectrum);
-
- // Create unidentified correlations histogram.
- fCorrelations = AliHistToolsDiHadronPID::MakeHist3D("fCorrelations","Correlations;#Delta#phi;#Delta#eta;p_{T} (GeV/c)",
- fNDPhiBins,-TMath::Pi()/2.,3.*TMath::Pi()/2.,
- fNDEtaBins,-1.6,1.6,
- nptbins, ptaxis);
- fOutputList->Add(fCorrelations);
-
- // Create unidentified mixed events histogram.
- fMixedEvents = AliHistToolsDiHadronPID::MakeHist3D("fMixedEvents","Mixed Events;#Delta#phi;#Delta#eta;p_{T} (GeV/c)",
- fNDPhiBins,-TMath::Pi()/2.,3.*TMath::Pi()/2.,
- fNDEtaBins,-1.6,1.6,
- nptbins, ptaxis);
- fOutputList->Add(fMixedEvents);
-
TString speciesname[] = {"Pion","Kaon","Proton"};
// Create TOF correlations histograms (DPhi,DEta,TOF).
if (fMakeTOFcorrelations) {
+ // Get the pT axis for the TOF PID correlations.
+ Double_t* ptaxis = fTrackCutsAssociated->GetPtAxisPID();
+ Int_t nptbins = fTrackCutsAssociated->GetNPtBinsPID();
+
+ // Create Pt spectrum histogram.
+ fPtSpectrumTOFbins = new TH1F("fPtSpectrumTOFbins","p_{T} Spectrum;p_{T} (GeV/c);Count",nptbins,ptaxis);
+ fOutputList->Add(fPtSpectrumTOFbins);
+
+ // Create unidentified correlations histogram.
+ fCorrelationsTOFbins = AliHistToolsDiHadronPID::MakeHist3D("fCorrelationsTOFbins","Correlations;#Delta#phi;#Delta#eta;p_{T} (GeV/c)",
+ fNDPhiBins,-TMath::Pi()/2.,3.*TMath::Pi()/2.,
+ fNDEtaBins,-1.6,1.6,
+ nptbins, ptaxis);
+ fOutputList->Add(fCorrelationsTOFbins);
+
+ // Create unidentified mixed events histogram.
+ fMixedEventsTOFbins = AliHistToolsDiHadronPID::MakeHist3D("fMixedEventsTOFbins","Mixed Events;#Delta#phi;#Delta#eta;p_{T} (GeV/c)",
+ fNDPhiBins,-TMath::Pi()/2.,3.*TMath::Pi()/2.,
+ fNDEtaBins,-1.6,1.6,
+ nptbins, ptaxis);
+ fOutputList->Add(fMixedEventsTOFbins);
+
+ // Create TOFPtaxis.
+ fTOFPtAxis = new TAxis(nptbins, ptaxis);
+ fTOFPtAxis->SetName("fTOFPtAxis");
+ fTOFPtAxis->SetTitle("p_{T} GeV/c");
+
+ // Create PID histograms.
fTOFhistos = new TObjArray(3);
fTOFhistos->SetOwner(kTRUE);
fTOFhistos->SetName("CorrelationsTOF");
+ if (fCalculateMismatch) {
+ fTOFmismatch = new TObjArray(3);
+ fTOFmismatch->SetOwner(kTRUE);
+ fTOFmismatch->SetName("MismatchTOF");
+ }
+
for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
- TObjArray* atmp = new TObjArray(fTOFPtAxis->GetNbins());
- atmp->SetOwner(kTRUE);
- atmp->SetName(speciesname[iSpecies].Data());
+ TObjArray* TOFhistosTmp = new TObjArray(fTOFPtAxis->GetNbins());
+ TOFhistosTmp->SetOwner(kTRUE);
+ TOFhistosTmp->SetName(speciesname[iSpecies].Data());
+
+ TObjArray* TOFmismatchTmp = 0x0;
+ if (fCalculateMismatch) {
+ TOFmismatchTmp = new TObjArray(fTOFPtAxis->GetNbins());
+ TOFmismatchTmp->SetOwner(kTRUE);
+ TOFmismatchTmp->SetName(speciesname[iSpecies].Data());
+ }
for (Int_t iBinPt = 1; iBinPt < (fTOFPtAxis->GetNbins() + 1); iBinPt++) {
Double_t TOFmin = fTrackCutsAssociated->GetTOFmin(iPtClass,iSpecies);
Double_t TOFmax = fTrackCutsAssociated->GetTOFmax(iPtClass,iSpecies);
- cout << "ptbin: "<< iBinPt << " class: " << iPtClass << " TOFBins: " << NBinsTOF << " min: " << TOFmin << " max: " << TOFmax << endl;
+ //cout << "ptbin: "<< iBinPt << " class: " << iPtClass << " TOFBins: " << NBinsTOF << " min: " << TOFmin << " max: " << TOFmax << endl;
+ // Correlation histogram.
TH3F* htmp = new TH3F(Form("fCorrelationsTOF_%i",iBinPt),
Form("%5.3f < p_{T} < %5.3f; #Delta#phi; #Delta#eta; t_{TOF} (ps)", fTOFPtAxis->GetBinLowEdge(iBinPt), fTOFPtAxis->GetBinUpEdge(iBinPt)),
fNDPhiBins, -TMath::Pi()/2., 3.*TMath::Pi()/2.,
fNDEtaBins, -1.6, 1.6, NBinsTOF, TOFmin, TOFmax);
htmp->SetDirectory(0);
- atmp->Add(htmp);
+ TOFhistosTmp->Add(htmp);
+ if (fCalculateMismatch) {
+ // Mismatch histogram.
+ TH1F* htmp2 = new TH1F(Form("fMismatchTOF_%i",iBinPt),
+ Form("%5.3f < p_{T} < %5.3f; t_{TOF} (ps)", fTOFPtAxis->GetBinLowEdge(iBinPt), fTOFPtAxis->GetBinUpEdge(iBinPt)),
+ NBinsTOF, TOFmin, TOFmax);
+ htmp2->SetDirectory(0);
+
+ TOFmismatchTmp->Add(htmp2);
+ }
}
- fTOFhistos->Add(atmp);
+ fTOFhistos->Add(TOFhistosTmp);
+ if (fCalculateMismatch) {fTOFmismatch->Add(TOFmismatchTmp);}
}
fOutputList->Add(fTOFhistos);
+ if (fCalculateMismatch) {fOutputList->Add(fTOFmismatch);}
}
Double_t ptarrayTOFTPC[16] = {2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6,
2.8, 3.0, 3.2, 3.4, 3.6, 3.8,
4.2, 4.6, 5.0};
- fTOFTPCPtAxis = new TAxis(15, ptarrayTOFTPC);
+ const Int_t nptbins = 15;
+ fTOFTPCPtAxis = new TAxis(nptbins, ptarrayTOFTPC);
fTOFTPCPtAxis->SetName("fTOFTPCPtAxis");
fTOFTPCPtAxis->SetTitle("p_{T} GeV/c");
- fOutputList->Add(fTOFTPCPtAxis);
+
+ // Create Pt spectrum histogram.
+ fPtSpectrumTOFTPCbins = new TH1F("fPtSpectrumTOFTPCbins","p_{T} Spectrum;p_{T} (GeV/c);Count",nptbins,ptarrayTOFTPC);
+ fOutputList->Add(fPtSpectrumTOFTPCbins);
+
+ // Create unidentified correlations histogram.
+ fCorrelationsTOFTPCbins = AliHistToolsDiHadronPID::MakeHist3D("fCorrelationsTOFTPCbins","Correlations;#Delta#phi;#Delta#eta;p_{T} (GeV/c)",
+ fNDPhiBins,-TMath::Pi()/2.,3.*TMath::Pi()/2.,
+ fNDEtaBins,-1.6,1.6,
+ nptbins, ptarrayTOFTPC);
+ fOutputList->Add(fCorrelationsTOFTPCbins);
+
+ // Create unidentified mixed events histogram.
+ fMixedEventsTOFTPCbins = AliHistToolsDiHadronPID::MakeHist3D("fMixedEventsTOFTPCbins","Mixed Events;#Delta#phi;#Delta#eta;p_{T} (GeV/c)",
+ fNDPhiBins,-TMath::Pi()/2.,3.*TMath::Pi()/2.,
+ fNDEtaBins,-1.6,1.6,
+ nptbins, ptarrayTOFTPC);
+ fOutputList->Add(fMixedEventsTOFTPCbins);
Int_t NBinsTOFTPC[4] = {32, 32, 60, 40};
Double_t minTOFTPC[4] = {-TMath::Pi()/2., -1.6, -1., -1.};
fTOFTPChistos = new TObjArray(3);
fTOFTPChistos->SetOwner(kTRUE);
- fTOFTPChistos->SetName("TOFTPChistos");
+ fTOFTPChistos->SetName("CorrelationsTOFTPC");
+
+ if (fCalculateMismatch) {
+ fTOFTPCmismatch = new TObjArray(3);
+ fTOFTPCmismatch->SetOwner(kTRUE);
+ fTOFTPCmismatch->SetName("MismatchTOFTPC");
+ }
for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
- TObjArray* atmp = new TObjArray(fTOFTPCPtAxis->GetNbins());
- atmp->SetOwner(kTRUE);
- atmp->SetName(speciesname[iSpecies].Data());
+ TObjArray* TOFTPChistosTmp = new TObjArray(fTOFTPCPtAxis->GetNbins());
+ TOFTPChistosTmp->SetOwner(kTRUE);
+ TOFTPChistosTmp->SetName(speciesname[iSpecies].Data());
+
+ TObjArray* TOFTPCmismatchTmp = 0x0;
+ if (fCalculateMismatch) {
+ TOFTPCmismatchTmp = new TObjArray(fTOFTPCPtAxis->GetNbins());
+ TOFTPCmismatchTmp->SetOwner(kTRUE);
+ TOFTPCmismatchTmp->SetName(speciesname[iSpecies].Data());
+ }
for (Int_t iBinPt = 1; iBinPt < (fTOFTPCPtAxis->GetNbins() + 1); iBinPt++) {
minTOFTPC[3] = TPCmin;
maxTOFTPC[3] = TPCmax;
- THnF* htmp = new THnF(Form("fCorrelationsTOF_%i",iBinPt),
+ THnF* htmp = new THnF(Form("fCorrelationsTOFTPC_%i",iBinPt),
Form("%5.3f < p_{T} < %5.3f", fTOFTPCPtAxis->GetBinLowEdge(iBinPt), fTOFTPCPtAxis->GetBinUpEdge(iBinPt)),
4, NBinsTOFTPC, minTOFTPC, maxTOFTPC);
(htmp->GetAxis(2))->SetTitle("t_{TOF} (ps)");
(htmp->GetAxis(3))->SetTitle("dE/dx (a.u.)");
- atmp->Add(htmp);
+ TOFTPChistosTmp->Add(htmp);
+
+ if (fCalculateMismatch) {
+ // Mismatch histogram.
+ TH2F* htmp2 = new TH2F(Form("fMismatchTOFTPC_%i",iBinPt),
+ Form("%5.3f < p_{T} < %5.3f; t_{TOF} (ps); dE/dx (a.u.)", fTOFTPCPtAxis->GetBinLowEdge(iBinPt), fTOFTPCPtAxis->GetBinUpEdge(iBinPt)),
+ NBinsTOFTPC[2], TOFmin, TOFmax, NBinsTOFTPC[3], TPCmin, TPCmax);
+ htmp2->SetDirectory(0);
+
+ TOFTPCmismatchTmp->Add(htmp2);
+
+ }
}
- fTOFTPChistos->Add(atmp);
+ fTOFTPChistos->Add(TOFTPChistosTmp);
+ if (fCalculateMismatch) {fTOFTPCmismatch->Add(TOFTPCmismatchTmp);}
}
fOutputList->Add(fTOFTPChistos);
+ if (fCalculateMismatch) {fOutputList->Add(fTOFTPCmismatch);}
}
// Load external TOF histograms if flag is set.
- if (fCalculateTOFmismatch) {LoadExtMismatchHistos();}
+ if (fCalculateMismatch) {LoadExtMismatchHistos();}
PostData(1,fOutputList);
pidtrack->SetDebugLevel(fDebug);
Double_t rndhittime = -1.e21;
- if (fCalculateTOFmismatch) rndhittime = GenerateRandomHit(pidtrack->Eta());
-
- // Fill p_T spectrum.
- fPtSpectrum->Fill(pidtrack->Pt());
+ if (fCalculateMismatch) rndhittime = GenerateRandomHit(pidtrack->Eta());
// Fill the trigger/associated tracks array.
if (fTrackCutsTrigger->IsSelectedData(pidtrack,rndhittime)) {fTriggerTracks->AddLast(pidtrack);}
- else if (fTrackCutsAssociated->IsSelectedData(pidtrack,rndhittime)) {fAssociatedTracks->AddLast(pidtrack);}
+ else if (fTrackCutsAssociated->IsSelectedData(pidtrack,rndhittime)) {
+
+ fAssociatedTracks->AddLast(pidtrack);
+
+ // Fill p_T spectrum.
+ if (fPtSpectrumTOFbins) fPtSpectrumTOFbins->Fill(pidtrack->Pt());
+ if (fPtSpectrumTOFTPCbins) fPtSpectrumTOFTPCbins->Fill(pidtrack->Pt());
+
+ // Fill mismatch histograms with associateds.
+ if (fCalculateMismatch && (rndhittime > -1.e20)) {
+
+ Double_t apt = pidtrack->Pt();
+
+ if (fMakeTOFcorrelations) {
+
+ for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
+
+ TObjArray* atmp = (TObjArray*)fTOFmismatch->At(iSpecies);
+ Int_t ptbin = fTOFPtAxis->FindBin(apt);
+
+ // Only fill if histogram exists in fTOFmismatch.
+ if ( !(ptbin < 1) && !(ptbin > fTOFPtAxis->GetNbins()) ) {
+
+ TH1F* htmp = (TH1F*)atmp->At(ptbin - 1);
+ htmp->Fill(rndhittime - pidtrack->GetTOFsignalExpected(iSpecies));
+
+ }
+ }
+ }
+
+ if (fMakeTOFTPCcorrelations) {
+
+ for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
+
+ TObjArray* atmp = (TObjArray*)fTOFTPCmismatch->At(iSpecies);
+ Int_t ptbin = fTOFTPCPtAxis->FindBin(apt);
+
+ // Only fill if histogram exists in fTOFTPCmismatch.
+ if ( !(ptbin < 1) && !(ptbin > fTOFTPCPtAxis->GetNbins()) ) {
+
+ TH2F* htmp = (TH2F*)atmp->At(ptbin - 1);
+ htmp->Fill(rndhittime - pidtrack->GetTOFsignalExpected(iSpecies), pidtrack->GetTPCsignalMinusExpected(iSpecies));
+
+ }
+ }
+ }
+
+ }
+
+ }
else {delete pidtrack;}
}
if (DPhi > 3.*TMath::Pi()/2.) {DPhi -= 2.*TMath::Pi();}
Double_t DEta = triggertrack->Eta() - associatedtrack->Eta();
- fCorrelations->Fill(DPhi,DEta,associatedtrack->Pt());
+ if (fCorrelationsTOFbins) fCorrelationsTOFbins->Fill(DPhi,DEta,associatedtrack->Pt());
+ if (fCorrelationsTOFTPCbins) fCorrelationsTOFTPCbins->Fill(DPhi,DEta,associatedtrack->Pt());
for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
}
}
- cout<<"Triggers: "<<fTriggerTracks->GetEntriesFast()<<" Associateds: "<<fAssociatedTracks->GetEntriesFast()<<endl;
-
- // Determine centrality of current event.
- TString centralityestimator = fEventCuts->GetCentralityEstimator();
- AliCentrality* currentcentrality = fCurrentAODEvent->GetCentrality();
- Float_t percentile = currentcentrality->GetCentralityPercentile(centralityestimator.Data());
+ //cout<<"Triggers: "<<fTriggerTracks->GetEntriesFast()<<" Associateds: "<<fAssociatedTracks->GetEntriesFast()<<endl;
// Determine vtxz of current event.
AliAODVertex* currentprimaryvertex = fCurrentAODEvent->GetPrimaryVertex();
Double_t vtxz = currentprimaryvertex->GetZ();
- AliEventPool* poolin = fPoolMgr->GetEventPool(percentile, vtxz);
- if (!poolin) {AliFatal(Form("No pool found for centrality = %f, vtxz = %f", percentile, vtxz));}
+ // Determine centrality of current event (for PbPb).
+ AliEventPool* poolin = 0x0;
+ Float_t percentile = -1.;
+ if (fEventCuts->GetIsPbPb()) {
+ TString centralityestimator = fEventCuts->GetCentralityEstimator();
+ AliCentrality* currentcentrality = fCurrentAODEvent->GetCentrality();
+ percentile = currentcentrality->GetCentralityPercentile(centralityestimator.Data());
+
+ poolin = fPoolMgr->GetEventPool(percentile, vtxz);
+ if (!poolin) {AliFatal(Form("No pool found for centrality = %f, vtxz = %f", percentile, vtxz));}
+ } else {
+ poolin = fPoolMgr->GetEventPool(0.5, vtxz); // There are no multiplicity bins for pp yet.
+ if (!poolin) {AliFatal(Form("No pool found for vtxz = %f", vtxz));}
+ }
+
// TObjArray* fGlobalTracksArray;
// Give a print out of the pool manager's contents.
// Mix events if there are enough events in the pool.
if (poolin->GetCurrentNEvents() >= fMinNEventsForMixing) {
- {cout << "Mixing Events." << endl;}
+ //{cout << "Mixing Events." << endl;}
// Loop over all events in the event pool.
for (Int_t iMixEvent = 0; iMixEvent < poolin->GetCurrentNEvents(); iMixEvent++) {
if (DPhi > 3.*TMath::Pi()/2.) {DPhi -= 2.*TMath::Pi();}
Double_t DEta = mixtrack->Eta() - associatedtrack->Eta();
- fMixedEvents->Fill(DPhi,DEta,associatedtrack->Pt());
+ if (fMixedEventsTOFbins) fMixedEventsTOFbins->Fill(DPhi,DEta,associatedtrack->Pt());
+ if (fMixedEventsTOFTPCbins) fMixedEventsTOFTPCbins->Fill(DPhi,DEta,associatedtrack->Pt());
}
}
if (DPhi > 3.*TMath::Pi()/2.) {DPhi -= 2.*TMath::Pi();}
Double_t DEta = triggertrack->Eta() - mixtrack->Eta();
- fMixedEvents->Fill(DPhi,DEta,mixtrack->Pt());
-
+ if (fMixedEventsTOFbins) fMixedEventsTOFbins->Fill(DPhi,DEta,mixtrack->Pt());
+ if (fMixedEventsTOFTPCbins) fMixedEventsTOFTPCbins->Fill(DPhi,DEta,mixtrack->Pt());
}
}
}
// Update the event pool.
- AliEventPool* poolout = fPoolMgr->GetEventPool(percentile, vtxz); // Get the buffer associated with the current centrality and z-vtx
- if (!poolout) AliFatal(Form("No pool found for centrality = %f, vtx_z = %f", percentile, vtxz));
+ AliEventPool* poolout = 0x0;
+ if (fEventCuts->GetIsPbPb()) {
+ poolout = fPoolMgr->GetEventPool(percentile, vtxz); // Get the buffer associated with the current centrality and z-vtx
+ if (!poolout) AliFatal(Form("No pool found for centrality = %f, vtx_z = %f", percentile, vtxz));
+ } else {
+ poolout = fPoolMgr->GetEventPool(0.5, vtxz); // Get the buffer associated with the current centrality and z-vtx
+ if (!poolout) AliFatal(Form("No pool found for vtx_z = %f", vtxz));
+ }
+
// Q: is it a problem that the fAssociatedTracks array can be bigger than the number of tracks inside?
if (fMixTriggers) {
fin = TFile::Open("alien:///alice/cern.ch/user/m/mveldhoe/rootfiles/TOFmismatchHistos.root");
if (!fin) {
AliWarning("Couln't open TOFmismatchHistos, will not calculate mismatches...");
- fCalculateTOFmismatch = kFALSE;
+ fCalculateMismatch = kFALSE;
return kFALSE;
}
TH1F* tmp1 = (TH1F*)fin->Get("hNewT0Fill");
if (!tmp1) {
AliWarning("Couln't find hNewT0Fill, will not calculate mismatches...");
- fCalculateTOFmismatch = kFALSE;
+ fCalculateMismatch = kFALSE;
return kFALSE;
}
TH2F* tmp2 = (TH2F*)fin->Get("hLvsEta");
if (!tmp2) {
AliWarning("Couln't find hLvsEta, will not calculate mismatches...");
- fCalculateTOFmismatch = kFALSE;
+ fCalculateMismatch = kFALSE;
return kFALSE;
}
for (Int_t iEtaBin = 1; iEtaBin < (nbinseta + 1); iEtaBin++) {
TH1F* tmp = (TH1F*)fLvsEta->ProjectionY(Form("LvsEtaProjection_%i",iEtaBin),iEtaBin,iEtaBin);
tmp->SetDirectory(0);
- fLvsEtaProjections->AddAt(tmp,iEtaBin);
+ fLvsEtaProjections->AddAt(tmp,iEtaBin - 1);
}
return kTRUE;
Double_t rndhittime = -1.e21;
// TOF mismatch flag is not turned on.
- if (!fCalculateTOFmismatch) {
- AliFatal("Called GenerateRandomHit() method, but flag fCalculateTOFmismatch not set.");
+ if (!fCalculateMismatch) {
+ AliFatal("Called GenerateRandomHit() method, but flag fCalculateMismatch not set.");
return rndhittime;
}
// Finding the bin of the eta.
TAxis* etaAxis = fLvsEta->GetXaxis();
Int_t etaBin = etaAxis->FindBin(eta);
- //cout<<"Eta: "<<eta<<" bin: "<<etaBin<<endl;
- const TH1F* lengthDistribution = (const TH1F*)fLvsEtaProjections->At(etaBin);
+ if (etaBin == 0 || (etaBin == etaAxis->GetNbins() + 1)) {return rndhittime;}
+
+ const TH1F* lengthDistribution = (const TH1F*)fLvsEtaProjections->At(etaBin - 1);
if (!lengthDistribution) {
AliFatal("length Distribution not found.");
void SetPoolSize(Int_t poolsize) {fPoolSize = poolsize;}
void SetMixEvents(Bool_t mixevents = kTRUE) {fMixEvents = mixevents;}
void SetMixTriggers(Bool_t mixtriggers = kTRUE) {fMixTriggers = mixtriggers;}
- void SetDebugLevel(Int_t debuglevel) {fDebug = debuglevel;}
+ //void SetDebugLevel(Int_t debuglevel) {fDebug = debuglevel;}
+ void SetCalculateMismatch(Bool_t calcmismatch) {fCalculateMismatch = calcmismatch;}
void SetMakeTOFCorrelations(Bool_t makeTOF) {fMakeTOFcorrelations = makeTOF;}
void SetMakeTOFTPCCorrelations(Bool_t makeTOFTPC) {fMakeTOFTPCcorrelations = makeTOFTPC;}
Int_t GetPoolSize() const {return fPoolSize;}
Bool_t GetMixEvents() const {return fMixEvents;}
Bool_t GetMixTriggers() const {return fMixTriggers;}
- Int_t GetDebugLevel() const {return fDebug;}
+ //Int_t GetDebugLevel() const {return fDebug;}
+ Bool_t GetCalculateMismatch() const {return fCalculateMismatch;}
Bool_t GetMakeTOFCorrelations() const {return fMakeTOFcorrelations;}
Bool_t GetMakeTOFTPCCorrelations() const {return fMakeTOFTPCcorrelations;}
TList* fOutputList; //! Output List.
// Histograms.
- TH1F* fPtSpectrum; //! Pt Spectrum.
- TH3F* fCorrelations; //! Correlations Histogram.
- TH3F* fMixedEvents; //! Mixed Events Histogram.
+ TH1F* fPtSpectrumTOFbins; //! Pt Spectrum.
+ TH3F* fCorrelationsTOFbins; //! Correlations Histogram.
+ TH3F* fMixedEventsTOFbins; //! Mixed Events Histogram.
+ TH1F* fPtSpectrumTOFTPCbins; //! Pt Spectrum.
+ TH3F* fCorrelationsTOFTPCbins; //! Correlations Histogram.
+ TH3F* fMixedEventsTOFTPCbins; //! Mixed Events Histogram.
+
TObjArray* fTOFhistos; //! Array holding all correlation functions with TOF information.
+ TObjArray* fTOFmismatch; //! Array holding mismatches, using fTOFPtAxis.
TAxis* fTOFPtAxis; //! P_t axis used for the TOF correlation histograms.
TObjArray* fTOFTPChistos; //! Array holding all correlation functions with TOF and TPC information.
+ TObjArray* fTOFTPCmismatch; //! Array holding mismatches, using fTOFTPCPtAxis.
TAxis* fTOFTPCPtAxis; //! P_t axis used for the TOF/ TPC correlation histograms.
// Settings.
Int_t fPoolSize; //
Bool_t fMixEvents; // NOT YET IMPLEMENTED.
Bool_t fMixTriggers; // If true, triggers are mixed, if not true, associateds are mixed.
- Bool_t fCalculateTOFmismatch; //
+ Bool_t fCalculateMismatch; //
// TOF mismatch stuff.
TH1F* fT0Fill; //
// Flags.
Bool_t fMakeTOFcorrelations; //
Bool_t fMakeTOFTPCcorrelations; //
- Int_t fDebug; // Debug flag.
+ //Int_t fDebug; // Debug flag.
- ClassDef(AliAnalysisTaskDiHadronPID,3);
+ ClassDef(AliAnalysisTaskDiHadronPID, 3);
};
-/************************************************************************* \r
-* Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. * \r
-* * \r
-* Author: The ALICE Off-line Project. * \r
-* Contributors are mentioned in the code where appropriate. * \r
-* * \r
-* Permission to use, copy, modify and distribute this software and its * \r
-* documentation strictly for non-commercial purposes is hereby granted * \r
-* without fee, provided that the above copyright notice appears in all * \r
-* copies and that both the copyright notice and this permission notice * \r
-* appear in the supporting documentation. The authors make no claims * \r
-* about the suitability of this software for any purpose. It is * \r
-* provided "as is" without express or implied warranty. * \r
-**************************************************************************/\r
-\r
-// -----------------------------------------------------------------------\r
-// Definitions the mathematical functions used in the DiHadronPID\r
-// analysis.\r
-// -----------------------------------------------------------------------\r
-// Author: Misha Veldhoen (misha.veldhoen@cern.ch)\r
-\r
-#include <iostream>\r
-\r
-#include "AliExternalTrackParam.h"\r
-#include "TCanvas.h"\r
-#include "TMath.h"\r
-#include "TF1.h"\r
-\r
-#include "AliFunctionsDiHadronPID.h"\r
-\r
-using namespace std;\r
-\r
-// -----------------------------------------------------------------------\r
-AliFunctionsDiHadronPID::AliFunctionsDiHadronPID()\r
-\r
-{\r
-\r
- // Constructor.\r
-\r
-} \r
-\r
-// -----------------------------------------------------------------------\r
-AliFunctionsDiHadronPID::~AliFunctionsDiHadronPID()\r
-\r
-{\r
-\r
- // Destructor.\r
-\r
-} \r
-\r
-// -----------------------------------------------------------------------\r
-Double_t AliFunctionsDiHadronPID::Gaussian1D(const Double_t *x, const Double_t *par) {\r
-\r
- // - Gaussian, I is the integral -\r
- // f(x) = I/(Sqrt(2*pi)*sigma) * exp{-(x-mu)^2/2sigma^2}\r
- // par = {I,mu,sigma}\r
-\r
- Double_t norm = par[0]/(TMath::Sqrt(2.*TMath::Pi())*par[2]);\r
- Double_t gaussian = TMath::Exp(-(x[0]-par[1])*(x[0]-par[1])/(2.*par[2]*par[2]));\r
-\r
- return (norm*gaussian);\r
-\r
-}\r
-\r
-// -----------------------------------------------------------------------\r
-Double_t AliFunctionsDiHadronPID::Gaussian1D(const Double_t xx, const Double_t integral, const Double_t mu, const Double_t sigma, const Double_t binwidth) {\r
-\r
- // The other implementation should make use of this one.\r
- Double_t norm = (binwidth*integral)/(TMath::Sqrt(2.*TMath::Pi())*sigma);\r
- Double_t gaussian = TMath::Exp(-(xx-mu)*(xx-mu)/(2.*sigma*sigma));\r
-\r
- return (norm*gaussian);\r
-\r
-}\r
-\r
-// -----------------------------------------------------------------------\r
-Double_t AliFunctionsDiHadronPID::Gaussian1DTail(const Double_t *x,const Double_t *par) {\r
-\r
- // Gaussian with exponential tail on the right, I is the integral.\r
- // For function definition see: FitFunctions.nb\r
-\r
- Double_t integral = par[0];\r
- Double_t mu = par[1]; // Top of the gaussian.\r
- Double_t kappa = par[1] + par[2]; // Point at which the gaussian becomes an exponential (w.r.t. to mu).\r
- Double_t sigma_x = par[3];\r
-\r
- if (mu >= kappa) return 0.; // Function becomes ill-defined.\r
-\r
- Double_t beta = sigma_x*sigma_x/(kappa-mu);\r
- Double_t BB = TMath::Exp( (kappa*kappa-mu*mu)/(2.*sigma_x*sigma_x) );\r
- Double_t norm1 = beta*TMath::Exp( -(mu-kappa)*(mu-kappa)/(2.*sigma_x*sigma_x) );\r
- Double_t norm2 = TMath::Sqrt(TMath::Pi()/2.)*sigma_x*TMath::Erfc( (mu-kappa)/(TMath::Sqrt2()*sigma_x) );\r
- Double_t norm = norm1 + norm2;\r
-\r
- Double_t funcleft = (integral/norm)*TMath::Exp(-(x[0]-mu)*(x[0]-mu)/(2.*sigma_x*sigma_x));\r
- Double_t funcright = (integral/norm)*BB*TMath::Exp(-x[0]/beta);\r
-\r
- if (x[0] <= kappa) return funcleft;\r
- else return funcright;\r
-\r
-}\r
-\r
-// -----------------------------------------------------------------------\r
-Double_t AliFunctionsDiHadronPID::Gaussian1DTail(const Double_t xx, const Double_t integral, const Double_t mu, const Double_t sigma, const Double_t tail, const Double_t binwidth) {\r
-\r
- // Gaussian with exponential tail on the right, I is the integral.\r
- // For function definition see: FitFunctions.nb\r
-\r
- Double_t kappa = mu + tail;\r
-\r
- if (mu >= kappa) return 0.; // Function becomes ill-defined.\r
-\r
- Double_t beta = sigma*sigma/(kappa-mu);\r
- Double_t BB = TMath::Exp( (kappa*kappa-mu*mu)/(2.*sigma*sigma) );\r
- Double_t norm1 = beta*TMath::Exp( -(mu-kappa)*(mu-kappa)/(2.*sigma*sigma) );\r
- Double_t norm2 = TMath::Sqrt(TMath::Pi()/2.)*sigma*TMath::Erfc( (mu-kappa)/(TMath::Sqrt2()*sigma) );\r
- Double_t norm = norm1 + norm2;\r
-\r
- Double_t funcleft = binwidth * (integral/norm)*TMath::Exp(-(xx-mu)*(xx-mu)/(2.*sigma*sigma));\r
- Double_t funcright = binwidth * (integral/norm)*BB*TMath::Exp(-xx/beta);\r
-\r
- if (xx <= kappa) return funcleft;\r
- else return funcright;\r
-\r
-}\r
-\r
-// -----------------------------------------------------------------------\r
-Double_t AliFunctionsDiHadronPID::Gaussian2D(const Double_t xx, const Double_t yy, const Double_t integral, \r
- const Double_t mux, const Double_t muy, const Double_t sigmax, const Double_t sigmay, \r
- const Double_t binwidthx, const Double_t binwidthy) {\r
-\r
- // 2D Gaussian.\r
- Double_t GaussianX = Gaussian1D(xx, 1., mux, sigmax, binwidthx);\r
- Double_t GaussianY = Gaussian1D(yy, 1., muy, sigmay, binwidthy);\r
-\r
- return integral * GaussianX * GaussianY; \r
-\r
-}\r
-\r
-// -----------------------------------------------------------------------\r
-Double_t AliFunctionsDiHadronPID::Gaussian2DTailX(const Double_t xx, const Double_t yy, const Double_t integral, \r
- const Double_t mux, const Double_t muy, const Double_t sigmax, const Double_t sigmay, \r
- const Double_t tailx, const Double_t binwidthx, const Double_t binwidthy) {\r
-\r
- // 2D Gaussian with exponential tail in X direction.\r
- Double_t GaussianTailX = Gaussian1DTail(xx, 1., mux, sigmax, tailx, binwidthx);\r
- Double_t GaussianY = Gaussian1D(yy, 1., muy, sigmay, binwidthy);\r
-\r
- return integral * GaussianTailX * GaussianY;\r
-\r
-}\r
-\r
-// -----------------------------------------------------------------------\r
-Double_t AliFunctionsDiHadronPID::Gaussian2DTailY(const Double_t xx, const Double_t yy, const Double_t integral, \r
- const Double_t mux, const Double_t muy, const Double_t sigmax, const Double_t sigmay, \r
- const Double_t taily, const Double_t binwidthx, const Double_t binwidthy) {\r
-\r
- // 2D Gaussian with exponential tail in Y direction.\r
- Double_t GaussianX = Gaussian1D(xx, 1., mux, sigmax, binwidthx);\r
- Double_t GaussianTailY = Gaussian1DTail(yy, 1., muy, sigmay, taily, binwidthy);\r
-\r
- return integral * GaussianX * GaussianTailY;\r
-\r
-}\r
-\r
-// -----------------------------------------------------------------------\r
-Double_t AliFunctionsDiHadronPID::Gaussian2DTailXY(const Double_t xx, const Double_t yy, const Double_t integral, \r
- const Double_t mux, const Double_t muy, const Double_t sigmax, const Double_t sigmay, \r
- const Double_t tailx, const Double_t taily, const Double_t binwidthx, const Double_t binwidthy) {\r
-\r
- // 2D Gaussian with exponential tail in X- and Y direction.\r
- Double_t GaussianTailX = Gaussian1DTail(xx, 1., mux, sigmax, tailx, binwidthx);\r
- Double_t GaussianTailY = Gaussian1DTail(yy, 1., muy, sigmay, taily, binwidthy);\r
-\r
- return integral * GaussianTailX * GaussianTailY;\r
-\r
-}\r
-\r
-// -----------------------------------------------------------------------\r
-Double_t AliFunctionsDiHadronPID::Exponent(const Double_t *x, const Double_t *par) {\r
-\r
- // f(x) = A * Exp(bx)\r
- // par = {A,b}\r
-\r
- return par[0]*TMath::Exp(par[1]*x[0]); \r
-\r
-}\r
-\r
-// -----------------------------------------------------------------------\r
-// COMBINED FUNCTIONS\r
-// -----------------------------------------------------------------------\r
-Double_t AliFunctionsDiHadronPID::SimpleTOFfit(const Double_t *x, const Double_t *par) {\r
-\r
- // Signal fitted with a Gaussian, mismatches by an exponent.\r
- return (Gaussian1D(x,&par[0]) + Exponent(x,&par[3]));\r
-\r
-}\r
-\r
-// -----------------------------------------------------------------------\r
-Double_t AliFunctionsDiHadronPID::SimpleTOFfitWithTail(const Double_t *x, const Double_t *par) {\r
-\r
- // Signal fitted with a Gaussian with a tail, mismatches by an exponent.\r
- return (Gaussian1D(x,&par[0]) + Exponent(x,&par[4]));\r
-\r
-}\r
-\r
-// -----------------------------------------------------------------------\r
-// PENALTY FUNCTIONS\r
-// -----------------------------------------------------------------------\r
-Double_t AliFunctionsDiHadronPID::PolyPenalty(const Double_t xx, const Double_t center, const Double_t flatwidth, const Int_t polyorder) {\r
-\r
- // Penalty function for a chi^2 fit. The function is defined as:\r
- // 1 for |xx - center| < flatwidth,\r
- // (|xx - center| - flatwidth) ^ polyorder for |xx - center| > flatwidth.\r
-\r
- Double_t fx = 1.;\r
- if (TMath::Abs(xx - center) > flatwidth) {\r
- fx = TMath::Power( (TMath::Abs(xx - center) - flatwidth), polyorder ) + 1.;\r
- }\r
-\r
- return fx;\r
-\r
-}\r
-\r
-// -----------------------------------------------------------------------\r
-TCanvas* AliFunctionsDiHadronPID::TestPolyPenalty(const Double_t range, const Double_t center, const Double_t flatwidth, const Int_t polyorder) {\r
-\r
- // Creates an example of the TestPolyPenalty function.\r
- TF1* tf = new TF1("tf",Form("AliFunctionsDiHadronPID::PolyPenalty(x,[0],[1],%i)",polyorder),-range,range);\r
- tf->SetParameters(center,flatwidth);\r
- TCanvas* cvs = TCanvas::MakeDefCanvas();\r
- tf->Draw();\r
-\r
- return cvs;\r
-\r
-}\r
-\r
-// -----------------------------------------------------------------------\r
-// PID Expected signal functions.\r
-// -----------------------------------------------------------------------\r
-Double_t AliFunctionsDiHadronPID::TOFExpTime(const Double_t pT, const Double_t eta, const Double_t mass) {\r
-\r
- // For description see ../Documents/TOFtime.tex\r
-\r
- Double_t AA = (2. * pT) / ( Charge() * BTPC() * GeVperkg() );\r
- Double_t BB = TMath::ASin( (Charge() * BTPC() * 0.01 * RTOF() * GeVperkg() ) / (2. * pT * C()) ); \r
- Double_t CC = TMath::Sqrt( mass*mass/(pT*pT) + TMath::CosH(eta)*TMath::CosH(eta) );\r
-\r
- return (1.e12*AA*BB*CC); // Time returned in ps.\r
-\r
-}\r
-\r
-// -----------------------------------------------------------------------\r
-Double_t AliFunctionsDiHadronPID::TPCExpdEdX(const Double_t pT, const Double_t eta, const Double_t mass) {\r
-\r
- // Not so neat solution, however the easiest for now.\r
-\r
- // Prameters taken from the constructor of AliTPCPIDResponse:\r
- Double_t MIP = 50.;\r
- Double_t Kp[5] = {0.0283086, 2.63394e+01, 5.04114e-11, 2.12543, 4.88663};\r
-\r
- Double_t betaGamma = TMath::Abs( (pT * TMath::CosH(eta)) / mass );\r
-\r
- // Implementation as in AliTPCPIDResponse.\r
- return MIP * AliExternalTrackParam::BetheBlochAleph(betaGamma,Kp[0],Kp[1],Kp[2],Kp[3],Kp[4]);\r
-\r
-}\r
+/*************************************************************************
+* Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
+* *
+* Author: The ALICE Off-line Project. *
+* Contributors are mentioned in the code where appropriate. *
+* *
+* Permission to use, copy, modify and distribute this software and its *
+* documentation strictly for non-commercial purposes is hereby granted *
+* without fee, provided that the above copyright notice appears in all *
+* copies and that both the copyright notice and this permission notice *
+* appear in the supporting documentation. The authors make no claims *
+* about the suitability of this software for any purpose. It is *
+* provided "as is" without express or implied warranty. *
+**************************************************************************/
+
+// -----------------------------------------------------------------------
+// Definitions the mathematical functions used in the DiHadronPID
+// analysis.
+// -----------------------------------------------------------------------
+// Author: Misha Veldhoen (misha.veldhoen@cern.ch)
+
+#include <iostream>
+
+#include "AliExternalTrackParam.h"
+#include "TCanvas.h"
+#include "TMath.h"
+#include "TF1.h"
+
+#include "AliFunctionsDiHadronPID.h"
+
+using namespace std;
+
+// -----------------------------------------------------------------------
+AliFunctionsDiHadronPID::AliFunctionsDiHadronPID()
+
+{
+
+ // Constructor.
+
+}
+
+// -----------------------------------------------------------------------
+AliFunctionsDiHadronPID::~AliFunctionsDiHadronPID()
+
+{
+
+ // Destructor.
+
+}
+
+// -----------------------------------------------------------------------
+Double_t AliFunctionsDiHadronPID::Gaussian1D(const Double_t *x, const Double_t *par) {
+
+ // - Gaussian, I is the integral -
+ // f(x) = I/(Sqrt(2*pi)*sigma) * exp{-(x-mu)^2/2sigma^2}
+ // par = {I,mu,sigma}
+
+ Double_t norm = par[0]/(TMath::Sqrt(2.*TMath::Pi())*par[2]);
+ Double_t gaussian = TMath::Exp(-(x[0]-par[1])*(x[0]-par[1])/(2.*par[2]*par[2]));
+
+ return (norm*gaussian);
+
+}
+
+// -----------------------------------------------------------------------
+Double_t AliFunctionsDiHadronPID::Gaussian1D(const Double_t xx, const Double_t integral, const Double_t mu, const Double_t sigma, const Double_t binwidth) {
+
+ // The other implementation should make use of this one.
+ Double_t norm = (binwidth*integral)/(TMath::Sqrt(2.*TMath::Pi())*sigma);
+ Double_t gaussian = TMath::Exp(-(xx-mu)*(xx-mu)/(2.*sigma*sigma));
+
+ return (norm*gaussian);
+
+}
+
+// -----------------------------------------------------------------------
+Double_t AliFunctionsDiHadronPID::Gaussian1DTail(const Double_t *x,const Double_t *par) {
+
+ // Gaussian with exponential tail on the right, I is the integral.
+ // For function definition see: FitFunctions.nb
+
+ Double_t integral = par[0];
+ Double_t mu = par[1]; // Top of the gaussian.
+ Double_t kappa = par[1] + par[2]; // Point at which the gaussian becomes an exponential (w.r.t. to mu).
+ Double_t sigma_x = par[3];
+
+ if (mu >= kappa) return 0.; // Function becomes ill-defined.
+
+ Double_t beta = sigma_x*sigma_x/(kappa-mu);
+ Double_t BB = TMath::Exp( (kappa*kappa-mu*mu)/(2.*sigma_x*sigma_x) );
+ Double_t norm1 = beta*TMath::Exp( -(mu-kappa)*(mu-kappa)/(2.*sigma_x*sigma_x) );
+ Double_t norm2 = TMath::Sqrt(TMath::Pi()/2.)*sigma_x*TMath::Erfc( (mu-kappa)/(TMath::Sqrt2()*sigma_x) );
+ Double_t norm = norm1 + norm2;
+
+ Double_t funcleft = (integral/norm)*TMath::Exp(-(x[0]-mu)*(x[0]-mu)/(2.*sigma_x*sigma_x));
+ Double_t funcright = (integral/norm)*BB*TMath::Exp(-x[0]/beta);
+
+ if (x[0] <= kappa) return funcleft;
+ else return funcright;
+
+}
+
+// -----------------------------------------------------------------------
+Double_t AliFunctionsDiHadronPID::Gaussian1DTail(const Double_t xx, const Double_t integral, const Double_t mu, const Double_t sigma, const Double_t tail, const Double_t binwidth) {
+
+ // Gaussian with exponential tail on the right, I is the integral.
+ // For function definition see: FitFunctions.nb
+
+ Double_t kappa = mu + tail;
+
+ if (mu >= kappa) return 0.; // Function becomes ill-defined.
+
+ Double_t beta = sigma*sigma/(kappa-mu);
+ Double_t BB = TMath::Exp( (kappa*kappa-mu*mu)/(2.*sigma*sigma) );
+ Double_t norm1 = beta*TMath::Exp( -(mu-kappa)*(mu-kappa)/(2.*sigma*sigma) );
+ Double_t norm2 = TMath::Sqrt(TMath::Pi()/2.)*sigma*TMath::Erfc( (mu-kappa)/(TMath::Sqrt2()*sigma) );
+ Double_t norm = norm1 + norm2;
+
+ Double_t funcleft = binwidth * (integral/norm)*TMath::Exp(-(xx-mu)*(xx-mu)/(2.*sigma*sigma));
+ Double_t funcright = binwidth * (integral/norm)*BB*TMath::Exp(-xx/beta);
+
+ if (xx <= kappa) return funcleft;
+ else return funcright;
+
+}
+
+// -----------------------------------------------------------------------
+Double_t AliFunctionsDiHadronPID::Gaussian2D(const Double_t xx, const Double_t yy, const Double_t integral,
+ const Double_t mux, const Double_t muy, const Double_t sigmax, const Double_t sigmay,
+ const Double_t binwidthx, const Double_t binwidthy) {
+
+ // 2D Gaussian.
+ Double_t GaussianX = Gaussian1D(xx, 1., mux, sigmax, binwidthx);
+ Double_t GaussianY = Gaussian1D(yy, 1., muy, sigmay, binwidthy);
+
+ return integral * GaussianX * GaussianY;
+
+}
+
+// -----------------------------------------------------------------------
+Double_t AliFunctionsDiHadronPID::Gaussian2DTailX(const Double_t xx, const Double_t yy, const Double_t integral,
+ const Double_t mux, const Double_t muy, const Double_t sigmax, const Double_t sigmay,
+ const Double_t tailx, const Double_t binwidthx, const Double_t binwidthy) {
+
+ // 2D Gaussian with exponential tail in X direction.
+ Double_t GaussianTailX = Gaussian1DTail(xx, 1., mux, sigmax, tailx, binwidthx);
+ Double_t GaussianY = Gaussian1D(yy, 1., muy, sigmay, binwidthy);
+
+ return integral * GaussianTailX * GaussianY;
+
+}
+
+// -----------------------------------------------------------------------
+Double_t AliFunctionsDiHadronPID::Gaussian2DTailY(const Double_t xx, const Double_t yy, const Double_t integral,
+ const Double_t mux, const Double_t muy, const Double_t sigmax, const Double_t sigmay,
+ const Double_t taily, const Double_t binwidthx, const Double_t binwidthy) {
+
+ // 2D Gaussian with exponential tail in Y direction.
+ Double_t GaussianX = Gaussian1D(xx, 1., mux, sigmax, binwidthx);
+ Double_t GaussianTailY = Gaussian1DTail(yy, 1., muy, sigmay, taily, binwidthy);
+
+ return integral * GaussianX * GaussianTailY;
+
+}
+
+// -----------------------------------------------------------------------
+Double_t AliFunctionsDiHadronPID::Gaussian2DTailXY(const Double_t xx, const Double_t yy, const Double_t integral,
+ const Double_t mux, const Double_t muy, const Double_t sigmax, const Double_t sigmay,
+ const Double_t tailx, const Double_t taily, const Double_t binwidthx, const Double_t binwidthy) {
+
+ // 2D Gaussian with exponential tail in X- and Y direction.
+ Double_t GaussianTailX = Gaussian1DTail(xx, 1., mux, sigmax, tailx, binwidthx);
+ Double_t GaussianTailY = Gaussian1DTail(yy, 1., muy, sigmay, taily, binwidthy);
+
+ return integral * GaussianTailX * GaussianTailY;
+
+}
+
+// -----------------------------------------------------------------------
+Double_t AliFunctionsDiHadronPID::Exponent(const Double_t *x, const Double_t *par) {
+
+ // f(x) = A * Exp(bx)
+ // par = {A,b}
+
+ return par[0]*TMath::Exp(par[1]*x[0]);
+
+}
+
+// -----------------------------------------------------------------------
+// COMBINED FUNCTIONS
+// -----------------------------------------------------------------------
+Double_t AliFunctionsDiHadronPID::SimpleTOFfit(const Double_t *x, const Double_t *par) {
+
+ // Signal fitted with a Gaussian, mismatches by an exponent.
+ return (Gaussian1D(x,&par[0]) + Exponent(x,&par[3]));
+
+}
+
+// -----------------------------------------------------------------------
+Double_t AliFunctionsDiHadronPID::SimpleTOFfitWithTail(const Double_t *x, const Double_t *par) {
+
+ // Signal fitted with a Gaussian with a tail, mismatches by an exponent.
+ return (Gaussian1D(x,&par[0]) + Exponent(x,&par[4]));
+
+}
+
+// -----------------------------------------------------------------------
+// PENALTY FUNCTIONS
+// -----------------------------------------------------------------------
+Double_t AliFunctionsDiHadronPID::PolyPenalty(const Double_t xx, const Double_t center, const Double_t flatwidth, const Int_t polyorder) {
+
+ // Penalty function for a chi^2 fit. The function is defined as:
+ // 1 for |xx - center| < flatwidth,
+ // (|xx - center| - flatwidth) ^ polyorder for |xx - center| > flatwidth.
+
+ Double_t fx = 1.;
+ if (TMath::Abs(xx - center) > flatwidth) {
+ fx = TMath::Power( (TMath::Abs(xx - center) - flatwidth), polyorder ) + 1.;
+ }
+
+ return fx;
+
+}
+
+// -----------------------------------------------------------------------
+TCanvas* AliFunctionsDiHadronPID::TestPolyPenalty(const Double_t range, const Double_t center, const Double_t flatwidth, const Int_t polyorder) {
+
+ // Creates an example of the TestPolyPenalty function.
+ TF1* tf = new TF1("tf",Form("AliFunctionsDiHadronPID::PolyPenalty(x,[0],[1],%i)",polyorder),-range,range);
+ tf->SetParameters(center,flatwidth);
+ TCanvas* cvs = TCanvas::MakeDefCanvas();
+ tf->Draw();
+
+ return cvs;
+
+}
+
+// -----------------------------------------------------------------------
+// PID Expected signal functions.
+// -----------------------------------------------------------------------
+Double_t AliFunctionsDiHadronPID::TOFExpTime(const Double_t pT, const Double_t eta, const Double_t mass) {
+
+ // For description see ../Documents/TOFtime.tex
+
+ Double_t AA = (2. * pT) / ( Charge() * BTPC() * GeVperkg() );
+ Double_t BB = TMath::ASin( (Charge() * BTPC() * 0.01 * RTOF() * GeVperkg() ) / (2. * pT * C()) );
+ Double_t CC = TMath::Sqrt( mass*mass/(pT*pT) + TMath::CosH(eta)*TMath::CosH(eta) );
+
+ return (1.e12*AA*BB*CC); // Time returned in ps.
+
+}
+
+// -----------------------------------------------------------------------
+Double_t AliFunctionsDiHadronPID::TPCExpdEdX(const Double_t pT, const Double_t eta, const Double_t mass) {
+
+ // Not so neat solution, however the easiest for now.
+
+ // Prameters taken from the constructor of AliTPCPIDResponse:
+ Double_t MIP = 50.;
+ Double_t Kp[5] = {0.0283086, 2.63394e+01, 5.04114e-11, 2.12543, 4.88663};
+
+ Double_t betaGamma = TMath::Abs( (pT * TMath::CosH(eta)) / mass );
+
+ // Implementation as in AliTPCPIDResponse.
+ return MIP * AliExternalTrackParam::BetheBlochAleph(betaGamma,Kp[0],Kp[1],Kp[2],Kp[3],Kp[4]);
+
+}
-#ifndef ALIFUNCTIONSDIHADRONPID_H\r
-#define ALIFUNCTIONSDIHADRONPID_H\r
-/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * \r
-* See cxx source for full Copyright notice */ \r
-/* $Id$ */\r
-\r
-class TCanvas;\r
-\r
-class AliFunctionsDiHadronPID {\r
-\r
-public:\r
- AliFunctionsDiHadronPID();\r
-\r
-protected:\r
- ~AliFunctionsDiHadronPID();\r
-\r
-public:\r
- // Natural constants.\r
- static Double_t Charge() {return 1.60217646e-19;} // (C)\r
- static Double_t C() {return 2.99792458e+8;} // (m/s)\r
- static Double_t Mpion() {return 0.13957018;} // (GeV/c^2)\r
- static Double_t Mkaon() {return 0.493667;} // (GeV/c^2)\r
- static Double_t Mproton() {return 0.938272046;} // (GeV/c^2)\r
- static Double_t Mdeuteron() {return 2.01410178*GeVperu();} // (GeV/c^2)\r
- static Double_t M(const Int_t species) {\r
- switch(species) {\r
- case 0:return Mpion();\r
- case 1:return Mkaon();\r
- case 2:return Mproton();\r
- case 3:return Mdeuteron();\r
- default:return -999.; \r
- }\r
- }\r
-\r
- // Conversions.\r
- static Double_t GeVperu() {return 0.931494061;} // (GeV/c^2) per u\r
- static Double_t GeVperkg() {return 5.608524e+26;} // (GeV/c^2) per kg\r
-\r
- // Detector paramters.\r
- static Double_t RTOF() {return 385.;} // Radius of TOF (cm).\r
- static Double_t BTPC() {return 0.5;} // Magnetic field in TPC (T = kg C^-1 s^-1).\r
-\r
- // Fit Functions.\r
- static Double_t Gaussian1D(const Double_t xx, const Double_t integral, const Double_t mu, const Double_t sigma, const Double_t binwidth = 1.);\r
- static Double_t Gaussian1DTail(const Double_t xx, const Double_t integral, const Double_t mu, const Double_t sigma, const Double_t tail, const Double_t binwidth = 1.);\r
-\r
- static Double_t Gaussian2D(const Double_t xx, const Double_t yy, const Double_t integral, \r
- const Double_t mux, const Double_t muy, const Double_t sigmax, const Double_t sigmay, \r
- const Double_t binwidthx = 1., const Double_t binwidthy = 1.);\r
-\r
- static Double_t Gaussian2DTailX(const Double_t xx, const Double_t yy, const Double_t integral, \r
- const Double_t mux, const Double_t muy, const Double_t sigmax, const Double_t sigmay, \r
- const Double_t tailx, const Double_t binwidthx = 1., const Double_t binwidthy = 1.);\r
-\r
- static Double_t Gaussian2DTailY(const Double_t xx, const Double_t yy, const Double_t integral, \r
- const Double_t mux, const Double_t muy, const Double_t sigmax, const Double_t sigmay, \r
- const Double_t taily, const Double_t binwidthx = 1., const Double_t binwidthy = 1.);\r
-\r
- static Double_t Gaussian2DTailXY(const Double_t xx, const Double_t yy, const Double_t integral, \r
- const Double_t mux, const Double_t muy, const Double_t sigmax, const Double_t sigmay, \r
- const Double_t tailx, const Double_t taily, const Double_t binwidthx = 1., const Double_t binwidthy = 1.);\r
-\r
- // FUNCTIONS THAT ARE NOT BEING USED ANYMORE.\r
- static Double_t Gaussian1D(const Double_t *x, const Double_t *par);\r
- static Double_t Gaussian1DTail(const Double_t *x, const Double_t *par); \r
- static Double_t Exponent(const Double_t *x, const Double_t *par);\r
- static Double_t SimpleTOFfit(const Double_t *x, const Double_t *par);\r
- static Double_t SimpleTOFfitWithTail(const Double_t *x, const Double_t *par);\r
-\r
- // Penalty Functions.\r
- static Double_t PolyPenalty(const Double_t xx, const Double_t center, const Double_t flatwidth, const Int_t polyorder);\r
- static TCanvas* TestPolyPenalty(const Double_t range = 3., const Double_t center = 1., const Double_t flatwidth = 1., const Int_t polyorder = 3);\r
-\r
- // PID Expected signal functions.\r
- static Double_t TOFExpTime(const Double_t pT, const Double_t eta, const Double_t mass);\r
- static Double_t TPCExpdEdX(const Double_t pT, const Double_t eta, const Double_t mass);\r
-\r
-};\r
-\r
-#endif\r
+#ifndef ALIFUNCTIONSDIHADRONPID_H
+#define ALIFUNCTIONSDIHADRONPID_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+* See cxx source for full Copyright notice */
+/* $Id$ */
+
+class TCanvas;
+
+class AliFunctionsDiHadronPID {
+
+public:
+ AliFunctionsDiHadronPID();
+
+protected:
+ ~AliFunctionsDiHadronPID();
+
+public:
+ // Natural constants.
+ static Double_t Charge() {return 1.60217646e-19;} // (C)
+ static Double_t C() {return 2.99792458e+8;} // (m/s)
+ static Double_t Mpion() {return 0.13957018;} // (GeV/c^2)
+ static Double_t Mkaon() {return 0.493667;} // (GeV/c^2)
+ static Double_t Mproton() {return 0.938272046;} // (GeV/c^2)
+ static Double_t Mdeuteron() {return 2.01410178*GeVperu();} // (GeV/c^2)
+ static Double_t M(const Int_t species) {
+ switch(species) {
+ case 0:return Mpion();
+ case 1:return Mkaon();
+ case 2:return Mproton();
+ case 3:return Mdeuteron();
+ default:return -999.;
+ }
+ }
+
+ // Conversions.
+ static Double_t GeVperu() {return 0.931494061;} // (GeV/c^2) per u
+ static Double_t GeVperkg() {return 5.608524e+26;} // (GeV/c^2) per kg
+
+ // Detector paramters.
+ static Double_t RTOF() {return 385.;} // Radius of TOF (cm).
+ static Double_t BTPC() {return 0.5;} // Magnetic field in TPC (T = kg C^-1 s^-1).
+
+ // Fit Functions.
+ static Double_t Gaussian1D(const Double_t xx, const Double_t integral, const Double_t mu, const Double_t sigma, const Double_t binwidth = 1.);
+ static Double_t Gaussian1DTail(const Double_t xx, const Double_t integral, const Double_t mu, const Double_t sigma, const Double_t tail, const Double_t binwidth = 1.);
+
+ static Double_t Gaussian2D(const Double_t xx, const Double_t yy, const Double_t integral,
+ const Double_t mux, const Double_t muy, const Double_t sigmax, const Double_t sigmay,
+ const Double_t binwidthx = 1., const Double_t binwidthy = 1.);
+
+ static Double_t Gaussian2DTailX(const Double_t xx, const Double_t yy, const Double_t integral,
+ const Double_t mux, const Double_t muy, const Double_t sigmax, const Double_t sigmay,
+ const Double_t tailx, const Double_t binwidthx = 1., const Double_t binwidthy = 1.);
+
+ static Double_t Gaussian2DTailY(const Double_t xx, const Double_t yy, const Double_t integral,
+ const Double_t mux, const Double_t muy, const Double_t sigmax, const Double_t sigmay,
+ const Double_t taily, const Double_t binwidthx = 1., const Double_t binwidthy = 1.);
+
+ static Double_t Gaussian2DTailXY(const Double_t xx, const Double_t yy, const Double_t integral,
+ const Double_t mux, const Double_t muy, const Double_t sigmax, const Double_t sigmay,
+ const Double_t tailx, const Double_t taily, const Double_t binwidthx = 1., const Double_t binwidthy = 1.);
+
+ // FUNCTIONS THAT ARE NOT BEING USED ANYMORE.
+ static Double_t Gaussian1D(const Double_t *x, const Double_t *par);
+ static Double_t Gaussian1DTail(const Double_t *x, const Double_t *par);
+ static Double_t Exponent(const Double_t *x, const Double_t *par);
+ static Double_t SimpleTOFfit(const Double_t *x, const Double_t *par);
+ static Double_t SimpleTOFfitWithTail(const Double_t *x, const Double_t *par);
+
+ // Penalty Functions.
+ static Double_t PolyPenalty(const Double_t xx, const Double_t center, const Double_t flatwidth, const Int_t polyorder);
+ static TCanvas* TestPolyPenalty(const Double_t range = 3., const Double_t center = 1., const Double_t flatwidth = 1., const Int_t polyorder = 3);
+
+ // PID Expected signal functions.
+ static Double_t TOFExpTime(const Double_t pT, const Double_t eta, const Double_t mass);
+ static Double_t TPCExpdEdX(const Double_t pT, const Double_t eta, const Double_t mass);
+
+};
+
+#endif
}
// -----------------------------------------------------------------------
-TObjArray* AliHistToolsDiHadronPID::CreateSpectraComparison(const char* name,
+TCanvas* AliHistToolsDiHadronPID::CreateSpectraComparison(const char* name,
const char* title, const TH1F* h1, const TH1F* h2, const Int_t markerstyle, const Bool_t logy) {
// - Creates a window comparing two histograms h1, and h2.
p1->cd();
legend->Draw("same");
+/*
// returning the created objects.
TObjArray* returnobjects = new TObjArray(6);
returnobjects->AddLast(c);
returnobjects->AddLast(spectra[0]);
returnobjects->AddLast(spectra[1]);
returnobjects->AddLast(division);
+*/
- return returnobjects;
+ return c;
}
class TF2;
class TH1F;
+class TH2;
class TH2F;
class TH3F;
+class TCanvas;
class AliHistToolsDiHadronPID {
static TH2F* Function2DToHist2D(const TF2* function, const TH2* grid);
// Histogram Visualization.
- static TObjArray* CreateSpectraComparison(const char* name, const char* title, const TH1F* h1, const TH1F* h2, const Int_t markerstyle = 8, const Bool_t logy = kTRUE);
+ static TCanvas* CreateSpectraComparison(const char* name, const char* title, const TH1F* h1, const TH1F* h2, const Int_t markerstyle = 8, const Bool_t logy = kTRUE);
private:
static Double_t* CreateAxis(const Int_t nbins, const Double_t min, const Double_t max);
AliError("No Track Supplied.");
}
- // Find the Global Track.
- if (fID >= 0) fAODGlobalTrack = fAODTrack;
-
- // Copy DCA and PID info.
- if (fAODGlobalTrack) {
- CopyFlags();
- if (fAODEvent) CopyDCAInfo();
- else AliError("Couln't find AOD Event.");
- CopyITSInfo();
- if (fPIDResponse) CopyTPCInfo();
- CopyTOFInfo();
- } else {
- AliError("Couldn't find Global Track.");
- }
-
- // Copy MC info.
- if (fAODMCParticle) {
- CopyMCInfo();
- }
+ // Copy the rest of the track parameters if the filtermap is nonzero.
+ // If fFiltermap == 0, then propagation to the DCA will result in a floating point error.
+ if (fFilterMap) {
+
+ // Find the Global Track.
+ if (fID >= 0) fAODGlobalTrack = fAODTrack;
+
+ // Copy DCA and PID info.
+ if (fAODGlobalTrack) {
+ CopyFlags();
+ if (fAODEvent) CopyDCAInfo();
+ else AliError("Couln't find AOD Event.");
+ CopyITSInfo();
+ if (fPIDResponse) CopyTPCInfo();
+ CopyTOFInfo();
+ } else {
+ AliError("Couldn't find Global Track.");
+ }
+
+ // Copy MC info.
+ if (fAODMCParticle) {
+ CopyMCInfo();
+ }
+
+ }
// Test
/* Double_t sigmaTOFProton = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(fAODTrack, AliPID::kProton));
triggercuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kAllProton);
triggercuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kPosProton);
triggercuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kNegProton);
- } triggercuts->SetDebugLevel(DebugLevel);
+ }
+ triggercuts->SetDebugLevel(DebugLevel);
DiHadronPIDTask->SetTrackCutsTrigger(triggercuts);
// Configure and add track cuts for associateds.
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff