//________________________________________________________________________\r
AliAnalysisTaskBFPsi::AliAnalysisTaskBFPsi(const char *name) \r
: AliAnalysisTaskSE(name),\r
+ fDebugLevel(kFALSE),\r
fArrayMC(0), //+++++++++++++\r
fBalance(0),\r
fRunShuffling(kFALSE),\r
fHistProbTPCTOFvsPtbeforePID(NULL),\r
fHistNSigmaTPCvsPtbeforePID(NULL), \r
fHistNSigmaTOFvsPtbeforePID(NULL), \r
+ fHistBetaVsdEdXbeforePID(NULL), //+++++++ \r
+ fHistNSigmaTPCTOFvsPtbeforePID(NULL), //++++++\r
fHistdEdxVsPTPCafterPID(NULL),\r
fHistBetavsPTOFafterPID(NULL), \r
fHistProbTPCvsPtafterPID(NULL), \r
fHistProbTPCTOFvsPtafterPID(NULL),\r
fHistNSigmaTPCvsPtafterPID(NULL), \r
fHistNSigmaTOFvsPtafterPID(NULL), \r
+ fHistBetaVsdEdXafterPID(NULL), //+++++++ \r
+ fHistNSigmaTPCTOFvsPtafterPID(NULL), //+++++++\r
+ fHistdEdxVsPTPCbeforePIDelectron(NULL), //+++++++\r
+ fHistNSigmaTPCvsPtbeforePIDelectron(NULL), //+++++++\r
+ fHistdEdxVsPTPCafterPIDelectron(NULL), //+++++++\r
+ fHistNSigmaTPCvsPtafterPIDelectron(NULL), //+++++++\r
fCentralityArrayBinsForCorrections(kCENTRALITY),\r
fPIDResponse(0x0),\r
fPIDCombined(0x0),\r
fCentralityPercentileMax(5.),\r
fImpactParameterMin(0.),\r
fImpactParameterMax(20.),\r
- fMultiplicityEstimator("V0M"),\r
+ fMultiplicityEstimator("V0A"),\r
fUseMultiplicity(kFALSE),\r
fNumberOfAcceptedTracksMin(0),\r
fNumberOfAcceptedTracksMax(10000),\r
fExcludeElectronsInMC(kFALSE),\r
fUseMCPdgCode(kFALSE),\r
fPDGCodeToBeAnalyzed(-1),\r
- fEventClass("EventPlane") \r
-{\r
+ fEventClass("EventPlane"), \r
+ fCustomBinning(""),\r
+ fHistVZEROAGainEqualizationMap(0),\r
+ fHistVZEROCGainEqualizationMap(0),\r
+ fHistVZEROChannelGainEqualizationMap(0) {\r
// Constructor\r
// Define input and output slots here\r
// Input slot #0 works with a TChain\r
13,-0.5,12.5,220,-5,105);\r
for(Int_t i = 1; i <= 13; i++){\r
fHistCentStats->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());\r
- //fHistCentStatsUsed->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data()); //++++++++++++++++++++++\r
+ //fHistCentStatsUsed->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());\r
}\r
fList->Add(fHistCentStats);\r
\r
- fHistCentStatsUsed = new TH2F("fHistCentStatsUsed","Centrality statistics;;Cent percentile", 1,-0.5,0.5,220,-5,105); //++++++++++++++++++++++\r
- fHistCentStatsUsed->GetXaxis()->SetBinLabel(1,fCentralityEstimator.Data()); //++++++++++++++++++++++\r
- fList->Add(fHistCentStatsUsed); //++++++++++++++++++++++\r
+ fHistCentStatsUsed = new TH2F("fHistCentStatsUsed","Centrality statistics;;Cent percentile", 1,-0.5,0.5,220,-5,105);\r
+ fHistCentStatsUsed->GetXaxis()->SetBinLabel(1,fCentralityEstimator.Data());\r
+ fList->Add(fHistCentStatsUsed);\r
\r
fHistTriggerStats = new TH1F("fHistTriggerStats","Trigger statistics;TriggerBit;N_{events}",1025,0,1025);\r
fList->Add(fHistTriggerStats);\r
fHistNumberOfAcceptedTracks = new TH2F("fHistNumberOfAcceptedTracks",";N_{acc.};Centrality percentile;Entries",4001,-0.5,4000.5,220,-5,105);\r
fList->Add(fHistNumberOfAcceptedTracks);\r
\r
- fHistMultiplicity = new TH1F("fHistMultiplicity",";N_{ch.};Entries",4001,-0.5,4000);\r
+ fHistMultiplicity = new TH1F("fHistMultiplicity",";N_{ch.};Entries",30001,-0.5,30000.5);\r
fList->Add(fHistMultiplicity);\r
\r
// Vertex distributions\r
//TPC vs VZERO multiplicity\r
fHistTPCvsVZEROMultiplicity = new TH2F("fHistTPCvsVZEROMultiplicity","VZERO vs TPC multiplicity",3001,-0.5,30000.5,4001,-0.5,4000.5);\r
if(fMultiplicityEstimator == "V0A") \r
- fHistTPCvsVZEROMultiplicity->GetXaxis()->SetTitle("VZERO multiplicity (a.u.)");\r
+ fHistTPCvsVZEROMultiplicity->GetXaxis()->SetTitle("VZERO-A multiplicity (a.u.)");\r
else if(fMultiplicityEstimator == "V0C") \r
- fHistTPCvsVZEROMultiplicity->GetXaxis()->SetTitle("VZERO multiplicity (a.u.)");\r
+ fHistTPCvsVZEROMultiplicity->GetXaxis()->SetTitle("VZERO-C multiplicity (a.u.)");\r
else \r
fHistTPCvsVZEROMultiplicity->GetXaxis()->SetTitle("VZERO multiplicity (a.u.)");\r
fList->Add(fHistTPCvsVZEROMultiplicity);\r
fHistRefTracks->GetXaxis()->SetBinLabel(i,gRefTrackName[i-1].Data());\r
fList->Add(fHistRefTracks);\r
\r
+ // Balance function histograms\r
+ // Initialize histograms if not done yet (including the custom binning)\r
+ if(!fBalance->GetHistNp()){\r
+ AliInfo("Histograms not yet initialized! --> Will be done now");\r
+ fBalance->SetCustomBinning(fCustomBinning);\r
+ fBalance->InitHistograms();\r
+ }\r
+\r
+ if(fRunShuffling) {\r
+ if(!fShuffledBalance->GetHistNp()) {\r
+ AliInfo("Histograms (shuffling) not yet initialized! --> Will be done now");\r
+ fShuffledBalance->SetCustomBinning(fCustomBinning);\r
+ fShuffledBalance->InitHistograms();\r
+ }\r
+ }\r
+\r
+ if(fRunMixing) {\r
+ if(!fMixedBalance->GetHistNp()) {\r
+ AliInfo("Histograms (mixing) not yet initialized! --> Will be done now");\r
+ fMixedBalance->SetCustomBinning(fCustomBinning);\r
+ fMixedBalance->InitHistograms();\r
+ }\r
+ }\r
+\r
// QA histograms for different cuts\r
fList->Add(fBalance->GetQAHistHBTbefore());\r
fList->Add(fBalance->GetQAHistHBTafter());\r
fList->Add(fBalance->GetQAHistQbefore());\r
fList->Add(fBalance->GetQAHistQafter());\r
\r
- // Balance function histograms\r
- // Initialize histograms if not done yet\r
- if(!fBalance->GetHistNp()){\r
- AliWarning("Histograms not yet initialized! --> Will be done now");\r
- AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");\r
- fBalance->InitHistograms();\r
- }\r
-\r
- if(fRunShuffling) {\r
- if(!fShuffledBalance->GetHistNp()) {\r
- AliWarning("Histograms (shuffling) not yet initialized! --> Will be done now");\r
- AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");\r
- fShuffledBalance->InitHistograms();\r
- }\r
- }\r
-\r
//for(Int_t a = 0; a < ANALYSIS_TYPES; a++){\r
fListBF->Add(fBalance->GetHistNp());\r
fListBF->Add(fBalance->GetHistNn());\r
Int_t poolsize = 1000; // Maximum number of events, ignored in the present implemented of AliEventPoolManager\r
\r
// centrality bins\r
- //Double_t centralityBins[] = {0.,1.,2.,3.,4.,5.,7.,10.,20.,30.,40.,50.,60.,70.,80.,100.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!\r
- Double_t centralityBins[] = {0.,1.,2.,3.,4.,5.,7.,10.,15.,20.,25.,30.,35.,40.,45.,50.,55.,60.,65.,70.,75.,80.};\r
- Double_t* centbins = centralityBins;\r
- Int_t nCentralityBins = sizeof(centralityBins) / sizeof(Double_t) - 1;\r
-\r
+ Double_t* centbins;\r
+ Int_t nCentralityBins;\r
+ if(fBalance->IsUseVertexBinning()){\r
+ centbins = fBalance->GetBinning(fBalance->GetBinningString(), "centralityVertex", nCentralityBins);\r
+ }\r
+ else{\r
+ centbins = fBalance->GetBinning(fBalance->GetBinningString(), "centrality", nCentralityBins);\r
+ }\r
+ \r
// multiplicity bins\r
- Double_t multiplicityBins[] = {0,10,20,30,40,50,60,70,80,100,100000}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!\r
- Double_t* multbins = multiplicityBins;\r
- Int_t nMultiplicityBins = sizeof(multiplicityBins) / sizeof(Double_t) - 1;\r
+ Double_t* multbins;\r
+ Int_t nMultiplicityBins;\r
+ multbins = fBalance->GetBinning(fBalance->GetBinningString(), "multiplicity", nMultiplicityBins);\r
\r
// Zvtx bins\r
- Double_t vertexBins[] = {-10., -7., -5., -3., -1., 1., 3., 5., 7., 10.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!\r
- Double_t* vtxbins = vertexBins;\r
- Int_t nVertexBins = sizeof(vertexBins) / sizeof(Double_t) - 1;\r
- \r
+ Double_t* vtxbins; \r
+ Int_t nVertexBins;\r
+ if(fBalance->IsUseVertexBinning()){\r
+ vtxbins = fBalance->GetBinning(fBalance->GetBinningString(), "vertexVertex", nVertexBins);\r
+ }\r
+ else{\r
+ vtxbins = fBalance->GetBinning(fBalance->GetBinningString(), "vertex", nVertexBins);\r
+ }\r
+\r
// Event plane angle (Psi) bins\r
- Double_t psiBins[] = {0.,45.,135.,215.,305.,360.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!\r
- Double_t* psibins = psiBins;\r
- Int_t nPsiBins = sizeof(psiBins) / sizeof(Double_t) - 1;\r
+ Double_t* psibins;\r
+ Int_t nPsiBins; \r
+ psibins = fBalance->GetBinning(fBalance->GetBinningString(), "eventPlane", nPsiBins);\r
+\r
\r
// run the event mixing also in bins of event plane (statistics!)\r
if(fRunMixingEventPlane){\r
fPIDCombined->SetDefaultTPCPriors();\r
\r
fHistdEdxVsPTPCbeforePID = new TH2D ("dEdxVsPTPCbefore","dEdxVsPTPCbefore", 1000, -10.0, 10.0, 1000, 0, 1000); \r
- fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePID); //addition \r
+ fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePID);\r
\r
fHistBetavsPTOFbeforePID = new TH2D ("BetavsPTOFbefore","BetavsPTOFbefore", 1000, -10.0, 10., 1000, 0, 1.2); \r
- fHistListPIDQA->Add(fHistBetavsPTOFbeforePID); //addition\r
+ fHistListPIDQA->Add(fHistBetavsPTOFbeforePID); \r
\r
fHistProbTPCvsPtbeforePID = new TH2D ("ProbTPCvsPtbefore","ProbTPCvsPtbefore", 1000, -10.0,10.0, 1000, 0, 2.0); \r
- fHistListPIDQA->Add(fHistProbTPCvsPtbeforePID); //addition \r
+ fHistListPIDQA->Add(fHistProbTPCvsPtbeforePID); \r
\r
fHistProbTOFvsPtbeforePID = new TH2D ("ProbTOFvsPtbefore","ProbTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0); \r
- fHistListPIDQA->Add(fHistProbTOFvsPtbeforePID); //addition \r
+ fHistListPIDQA->Add(fHistProbTOFvsPtbeforePID);\r
\r
fHistProbTPCTOFvsPtbeforePID =new TH2D ("ProbTPCTOFvsPtbefore","ProbTPCTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0); \r
- fHistListPIDQA->Add(fHistProbTPCTOFvsPtbeforePID); //addition \r
+ fHistListPIDQA->Add(fHistProbTPCTOFvsPtbeforePID);\r
\r
fHistNSigmaTPCvsPtbeforePID = new TH2D ("NSigmaTPCvsPtbefore","NSigmaTPCvsPtbefore", 1000, -10, 10, 1000, 0, 500); \r
- fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePID); //addition \r
+ fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePID);\r
\r
fHistNSigmaTOFvsPtbeforePID = new TH2D ("NSigmaTOFvsPtbefore","NSigmaTOFvsPtbefore", 1000, -10, 10, 1000, 0, 500); \r
- fHistListPIDQA->Add(fHistNSigmaTOFvsPtbeforePID); //addition \r
+ fHistListPIDQA->Add(fHistNSigmaTOFvsPtbeforePID); \r
+\r
+ fHistBetaVsdEdXbeforePID = new TH2D ("BetaVsdEdXbefore","BetaVsdEdXbefore", 1000, 0., 1000, 1000, 0, 1.2); \r
+ fHistListPIDQA->Add(fHistBetaVsdEdXbeforePID); //++++++++\r
+ \r
+ fHistNSigmaTPCTOFvsPtbeforePID = new TH2D ("NSigmaTPCTOFvsPtbefore","NSigmaTPCTOFvsPtbefore", 1000, -10., 10., 1000, 0, 500); \r
+ fHistListPIDQA->Add(fHistNSigmaTPCTOFvsPtbeforePID); //++++++++\r
\r
fHistdEdxVsPTPCafterPID = new TH2D ("dEdxVsPTPCafter","dEdxVsPTPCafter", 1000, -10, 10, 1000, 0, 1000); \r
- fHistListPIDQA->Add(fHistdEdxVsPTPCafterPID); //addition \r
+ fHistListPIDQA->Add(fHistdEdxVsPTPCafterPID);\r
\r
fHistBetavsPTOFafterPID = new TH2D ("BetavsPTOFafter","BetavsPTOFafter", 1000, -10, 10, 1000, 0, 1.2); \r
- fHistListPIDQA->Add(fHistBetavsPTOFafterPID); //addition \r
+ fHistListPIDQA->Add(fHistBetavsPTOFafterPID); \r
\r
fHistProbTPCvsPtafterPID = new TH2D ("ProbTPCvsPtafter","ProbTPCvsPtafter", 1000, -10, 10, 1000, 0, 2); \r
- fHistListPIDQA->Add(fHistProbTPCvsPtafterPID); //addition \r
+ fHistListPIDQA->Add(fHistProbTPCvsPtafterPID);\r
\r
fHistProbTOFvsPtafterPID = new TH2D ("ProbTOFvsPtafter","ProbTOFvsPtafter", 1000, -10, 10, 1000, 0, 2); \r
- fHistListPIDQA->Add(fHistProbTOFvsPtafterPID); //addition \r
+ fHistListPIDQA->Add(fHistProbTOFvsPtafterPID); \r
\r
fHistProbTPCTOFvsPtafterPID =new TH2D ("ProbTPCTOFvsPtafter","ProbTPCTOFvsPtafter", 1000, -50, 50, 1000, 0, 2.0); \r
- fHistListPIDQA->Add(fHistProbTPCTOFvsPtafterPID); //addition \r
+ fHistListPIDQA->Add(fHistProbTPCTOFvsPtafterPID);\r
\r
fHistNSigmaTPCvsPtafterPID = new TH2D ("NSigmaTPCvsPtafter","NSigmaTPCvsPtafter", 1000, -10, 10, 1000, 0, 500); \r
- fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPID); //addition \r
+ fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPID);\r
\r
fHistNSigmaTOFvsPtafterPID = new TH2D ("NSigmaTOFvsPtafter","NSigmaTOFvsPtafter", 1000, -10, 10, 1000, 0, 500); \r
- fHistListPIDQA->Add(fHistNSigmaTOFvsPtafterPID); //addition \r
+ fHistListPIDQA->Add(fHistNSigmaTOFvsPtafterPID);\r
+\r
+ fHistBetaVsdEdXafterPID = new TH2D ("BetaVsdEdXafter","BetaVsdEdXafter", 1000, 0., 1000, 1000, 0, 1.2); \r
+ fHistListPIDQA->Add(fHistBetaVsdEdXafterPID); //++++++++\r
+\r
+ fHistNSigmaTPCTOFvsPtafterPID = new TH2D ("NSigmaTPCTOFvsPtafter","NSigmaTPCTOFvsPtafter", 1000, -10., 10., 1000, 0, 500); \r
+ fHistListPIDQA->Add(fHistNSigmaTPCTOFvsPtafterPID); //++++++++\r
}\r
\r
// for electron rejection only TPC nsigma histograms\r
- if(!fUsePID && fElectronRejection) {\r
+ if(fElectronRejection) {\r
\r
- fHistdEdxVsPTPCbeforePID = new TH2D ("dEdxVsPTPCbefore","dEdxVsPTPCbefore", 1000, -10.0, 10.0, 1000, 0, 1000); \r
- fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePID); //addition \r
+ fHistdEdxVsPTPCbeforePIDelectron = new TH2D ("dEdxVsPTPCbeforeelectron","dEdxVsPTPCbeforeelectron", 1000, -10.0, 10.0, 1000, 0, 1000); \r
+ fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePIDelectron);\r
\r
- fHistNSigmaTPCvsPtbeforePID = new TH2D ("NSigmaTPCvsPtbefore","NSigmaTPCvsPtbefore", 1000, -10, 10, 1000, 0, 500); \r
- fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePID); //addition \r
+ fHistNSigmaTPCvsPtbeforePIDelectron = new TH2D ("NSigmaTPCvsPtbeforeelectron","NSigmaTPCvsPtbeforeelectron", 1000, -10, 10, 1000, 0, 500); \r
+ fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePIDelectron);\r
\r
- fHistdEdxVsPTPCafterPID = new TH2D ("dEdxVsPTPCafter","dEdxVsPTPCafter", 1000, -10, 10, 1000, 0, 1000); \r
- fHistListPIDQA->Add(fHistdEdxVsPTPCafterPID); //addition \r
+ fHistdEdxVsPTPCafterPIDelectron = new TH2D ("dEdxVsPTPCafterelectron","dEdxVsPTPCafterelectron", 1000, -10, 10, 1000, 0, 1000); \r
+ fHistListPIDQA->Add(fHistdEdxVsPTPCafterPIDelectron);\r
\r
- fHistNSigmaTPCvsPtafterPID = new TH2D ("NSigmaTPCvsPtafter","NSigmaTPCvsPtafter", 1000, -10, 10, 1000, 0, 500); \r
- fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPID); //addition \r
+ fHistNSigmaTPCvsPtafterPIDelectron = new TH2D ("NSigmaTPCvsPtafterelectron","NSigmaTPCvsPtafterelectron", 1000, -10, 10, 1000, 0, 500); \r
+ fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPIDelectron); \r
}\r
//====================PID========================//\r
\r
if(fRunMixing) PostData(4, fListBFM);\r
if(fUsePID || fElectronRejection) PostData(5, fHistListPIDQA); //PID\r
\r
+ AliInfo("Finished setting up the Output");\r
+\r
TH1::AddDirectory(oldStatus);\r
}\r
\r
\r
TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
Int_t gNumberOfAcceptedTracks = 0;\r
- Double_t gCentrality = -1.;\r
+ Double_t lMultiplicityVar = -999.; //-1\r
Double_t gReactionPlane = -1.; \r
Float_t bSign = 0.;\r
\r
eventMain = dynamic_cast<AliVEvent*>(MCEvent()); \r
}\r
else{\r
- eventMain = dynamic_cast<AliVEvent*>(InputEvent()); \r
- \r
+ eventMain = dynamic_cast<AliVEvent*>(InputEvent()); \r
// for HBT like cuts need magnetic field sign\r
bSign = (eventMain->GetMagneticField() > 0) ? 1 : -1;\r
}\r
fPIDResponse = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetPIDResponse();\r
if (!fPIDResponse) AliFatal("This Task needs the PID response attached to the inputHandler");\r
}\r
- \r
+ \r
// check event cuts and fill event histograms\r
- if((gCentrality = IsEventAccepted(eventMain)) < 0){\r
+ if((lMultiplicityVar = IsEventAccepted(eventMain)) < 0){ \r
return;\r
}\r
- \r
- //Compute Multiplicity or Centrality variable\r
- Double_t lMultiplicityVar = GetRefMultiOrCentrality( eventMain );\r
-\r
// get the reaction plane\r
if(fEventClass != "Multiplicity") {\r
gReactionPlane = GetEventPlane(eventMain);\r
// Fills Event statistics histograms\r
\r
Bool_t isSelectedMain = kTRUE;\r
- Float_t gCentrality = -1.;\r
Float_t gRefMultiplicity = -1.;\r
TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
\r
AliMCEvent *mcevent = dynamic_cast<AliMCEvent*>(event);\r
\r
- fHistEventStats->Fill(1,gCentrality); //all events\r
+ fHistEventStats->Fill(1,gRefMultiplicity); //all events\r
\r
// check first event in chunk (is not needed for new reconstructions)\r
if(fCheckFirstEventInChunk){\r
AliAnalysisUtils ut;\r
if(ut.IsFirstEventInChunk(event)) \r
return -1.;\r
- fHistEventStats->Fill(6,gCentrality); \r
+ fHistEventStats->Fill(6,gRefMultiplicity); \r
}\r
-\r
// check for pile-up event\r
if(fCheckPileUp){\r
AliAnalysisUtils ut;\r
ut.SetUseOutOfBunchPileUp(kTRUE);\r
if(ut.IsPileUpEvent(event))\r
return -1.;\r
- fHistEventStats->Fill(7,gCentrality); \r
+ fHistEventStats->Fill(7,gRefMultiplicity); \r
}\r
\r
// Event trigger bits\r
isSelectedMain = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r
\r
if(isSelectedMain) {\r
- fHistEventStats->Fill(2,gCentrality); //triggered events\r
-\r
- //Centrality stuff \r
- if(fUseCentrality) {\r
- if(gAnalysisLevel == "AOD"|| gAnalysisLevel == "MCAOD" || gAnalysisLevel == "MCAODrec") { //centrality in AOD header //+++++++++++\r
- AliAODHeader *header = (AliAODHeader*) event->GetHeader();\r
- if(header){\r
- gCentrality = header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());\r
-\r
- // QA for centrality estimators\r
- fHistCentStats->Fill(0.,header->GetCentralityP()->GetCentralityPercentile("V0M"));\r
- fHistCentStats->Fill(1.,header->GetCentralityP()->GetCentralityPercentile("V0A"));\r
- fHistCentStats->Fill(2.,header->GetCentralityP()->GetCentralityPercentile("V0C"));\r
- fHistCentStats->Fill(3.,header->GetCentralityP()->GetCentralityPercentile("FMD"));\r
- fHistCentStats->Fill(4.,header->GetCentralityP()->GetCentralityPercentile("TRK"));\r
- fHistCentStats->Fill(5.,header->GetCentralityP()->GetCentralityPercentile("TKL")); \r
- fHistCentStats->Fill(6.,header->GetCentralityP()->GetCentralityPercentile("CL0"));\r
- fHistCentStats->Fill(7.,header->GetCentralityP()->GetCentralityPercentile("CL1"));\r
- fHistCentStats->Fill(8.,header->GetCentralityP()->GetCentralityPercentile("ZNA"));\r
- fHistCentStats->Fill(9.,header->GetCentralityP()->GetCentralityPercentile("ZPA"));\r
- fHistCentStats->Fill(10.,header->GetCentralityP()->GetCentralityPercentile("V0MvsFMD"));\r
- fHistCentStats->Fill(11.,header->GetCentralityP()->GetCentralityPercentile("TKLvsV0M"));\r
- fHistCentStats->Fill(12.,header->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC"));\r
-\r
- // Centrality estimator USED ++++++++++++++++++++++++++++++\r
- fHistCentStatsUsed->Fill(0.,header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data()));\r
- \r
- // centrality QA (V0M)\r
- fHistV0M->Fill(event->GetVZEROData()->GetMTotV0A(), event->GetVZEROData()->GetMTotV0C());\r
- \r
- // centrality QA (reference tracks)\r
- fHistRefTracks->Fill(0.,header->GetRefMultiplicity());\r
- fHistRefTracks->Fill(1.,header->GetRefMultiplicityPos());\r
- fHistRefTracks->Fill(2.,header->GetRefMultiplicityNeg());\r
- fHistRefTracks->Fill(3.,header->GetTPConlyRefMultiplicity());\r
- fHistRefTracks->Fill(4.,header->GetNumberOfITSClusters(0));\r
- fHistRefTracks->Fill(5.,header->GetNumberOfITSClusters(1));\r
- fHistRefTracks->Fill(6.,header->GetNumberOfITSClusters(2));\r
- fHistRefTracks->Fill(7.,header->GetNumberOfITSClusters(3));\r
- fHistRefTracks->Fill(8.,header->GetNumberOfITSClusters(4));\r
- }//AOD header\r
- }//AOD\r
-\r
- else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD"){ // centrality class for ESDs or MC-ESDs\r
- AliCentrality *centrality = event->GetCentrality();\r
- gCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
-\r
- // QA for centrality estimators\r
- fHistCentStats->Fill(0.,centrality->GetCentralityPercentile("V0M"));\r
- fHistCentStats->Fill(1.,centrality->GetCentralityPercentile("V0A"));\r
- fHistCentStats->Fill(2.,centrality->GetCentralityPercentile("V0C"));\r
- fHistCentStats->Fill(3.,centrality->GetCentralityPercentile("FMD"));\r
- fHistCentStats->Fill(4.,centrality->GetCentralityPercentile("TRK"));\r
- fHistCentStats->Fill(5.,centrality->GetCentralityPercentile("TKL"));\r
- fHistCentStats->Fill(6.,centrality->GetCentralityPercentile("CL0"));\r
- fHistCentStats->Fill(7.,centrality->GetCentralityPercentile("CL1"));\r
- fHistCentStats->Fill(8.,centrality->GetCentralityPercentile("ZNA"));\r
- fHistCentStats->Fill(9.,centrality->GetCentralityPercentile("ZPA"));\r
- fHistCentStats->Fill(10.,centrality->GetCentralityPercentile("V0MvsFMD"));\r
- fHistCentStats->Fill(11.,centrality->GetCentralityPercentile("TKLvsV0M"));\r
- fHistCentStats->Fill(12.,centrality->GetCentralityPercentile("ZEMvsZDC"));\r
-\r
- // Centrality estimator USED ++++++++++++++++++++++++++++++\r
- fHistCentStatsUsed->Fill(0.,centrality->GetCentralityPercentile(fCentralityEstimator.Data()));\r
-\r
- // centrality QA (V0M)\r
- fHistV0M->Fill(event->GetVZEROData()->GetMTotV0A(), event->GetVZEROData()->GetMTotV0C());\r
- }//ESD\r
- else if(gAnalysisLevel == "MC"){\r
- Double_t gImpactParameter = 0.;\r
- if(mcevent) {\r
- AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(dynamic_cast<AliMCEvent*>(mcevent)->GenEventHeader());\r
- if(headerH){\r
- gImpactParameter = headerH->ImpactParameter();\r
- gCentrality = gImpactParameter;\r
- }//MC header\r
- }//MC event cast\r
- }//MC\r
- else{\r
- gCentrality = -1.;\r
- }\r
- }//centrality\r
-\r
- //Multiplicity stuff \r
- if(fUseMultiplicity)\r
- gRefMultiplicity = GetRefMultiOrCentrality(event);\r
-\r
+ fHistEventStats->Fill(2,gRefMultiplicity); //triggered events\r
+ \r
// Event Vertex MC\r
if(gAnalysisLevel == "MC") {\r
if(!event) {\r
//gVertexArray.At(0),\r
//gVertexArray.At(1),\r
//gVertexArray.At(2));\r
- if(fUseMultiplicity) \r
- fHistEventStats->Fill(3,gRefMultiplicity); //events with a proper vertex\r
- else \r
- fHistEventStats->Fill(3,gCentrality); //events with a proper vertex\r
+ fHistEventStats->Fill(3,gRefMultiplicity); //events with a proper vertex\r
if(TMath::Abs(gVertexArray.At(0)) < fVxMax) {\r
if(TMath::Abs(gVertexArray.At(1)) < fVyMax) {\r
if(TMath::Abs(gVertexArray.At(2)) < fVzMax) {\r
- if(fUseMultiplicity) \r
- fHistEventStats->Fill(4,gRefMultiplicity);//analyzed events\r
- else \r
- fHistEventStats->Fill(4,gCentrality); //analyzed events\r
+ fHistEventStats->Fill(4,gRefMultiplicity);//analyzed events\r
+\r
+ // get the reference multiplicty or centrality\r
+ gRefMultiplicity = GetRefMultiOrCentrality(event);\r
+\r
fHistVx->Fill(gVertexArray.At(0));\r
fHistVy->Fill(gVertexArray.At(1));\r
- fHistVz->Fill(gVertexArray.At(2),gCentrality);\r
+ fHistVz->Fill(gVertexArray.At(2),gRefMultiplicity);\r
\r
// take only events inside centrality class\r
if(fUseCentrality) {\r
- if((fImpactParameterMin < gCentrality) && (fImpactParameterMax > gCentrality)){\r
- fHistEventStats->Fill(5,gCentrality); //events with correct centrality\r
- return gCentrality; \r
+ if((fImpactParameterMin < gRefMultiplicity) && (fImpactParameterMax > gRefMultiplicity)){\r
+ fHistEventStats->Fill(5,gRefMultiplicity); //events with correct centrality\r
+ return gRefMultiplicity; \r
}//centrality class\r
}\r
// take events only within the same multiplicity class\r
else if(fUseMultiplicity) {\r
- if((gRefMultiplicity > fNumberOfAcceptedTracksMin)||(gRefMultiplicity < fNumberOfAcceptedTracksMax)) {\r
+ if((gRefMultiplicity > fNumberOfAcceptedTracksMin) && (gRefMultiplicity < fNumberOfAcceptedTracksMax)) {\r
fHistEventStats->Fill(5,gRefMultiplicity); //events with correct multiplicity\r
return gRefMultiplicity;\r
}\r
vertex->GetCovarianceMatrix(fCov);\r
if(vertex->GetNContributors() > 0) {\r
if(fCov[5] != 0) {\r
- if(fUseMultiplicity) \r
- fHistEventStats->Fill(3,gRefMultiplicity); //proper vertex\r
- else if(fUseCentrality)\r
- fHistEventStats->Fill(3,gCentrality); //proper vertex\r
+ fHistEventStats->Fill(3,gRefMultiplicity); //proper vertex\r
if(TMath::Abs(vertex->GetX()) < fVxMax) {\r
if(TMath::Abs(vertex->GetY()) < fVyMax) {\r
if(TMath::Abs(vertex->GetZ()) < fVzMax) {\r
- if(fUseMultiplicity) {\r
- //cout<<"Filling 4 for multiplicity..."<<endl;\r
- fHistEventStats->Fill(4,gRefMultiplicity);//analyzed events\r
- }\r
- else if(fUseCentrality) {\r
- //cout<<"Filling 4 for centrality..."<<endl;\r
- fHistEventStats->Fill(4,gCentrality); //analyzed events\r
- }\r
+ fHistEventStats->Fill(4,gRefMultiplicity);//analyzed events\r
+\r
+ // get the reference multiplicty or centrality\r
+ gRefMultiplicity = GetRefMultiOrCentrality(event);\r
+ \r
fHistVx->Fill(vertex->GetX());\r
fHistVy->Fill(vertex->GetY());\r
- fHistVz->Fill(vertex->GetZ(),gCentrality);\r
+ fHistVz->Fill(vertex->GetZ(),gRefMultiplicity);\r
\r
// take only events inside centrality class\r
if(fUseCentrality) {\r
- //cout<<"Centrality..."<<endl;\r
- if((gCentrality > fCentralityPercentileMin) && (gCentrality < fCentralityPercentileMax)){\r
- fHistEventStats->Fill(5,gCentrality); //events with correct centrality\r
- return gCentrality; \r
+ if((gRefMultiplicity > fCentralityPercentileMin) && (gRefMultiplicity < fCentralityPercentileMax)){\r
+ fHistEventStats->Fill(5,gRefMultiplicity); //events with correct centrality\r
+ return gRefMultiplicity; \r
}//centrality class\r
}\r
// take events only within the same multiplicity class\r
else if(fUseMultiplicity) {\r
- //cout<<"N(min): "<<fNumberOfAcceptedTracksMin<<" - N(max): "<<fNumberOfAcceptedTracksMax<<" - Nref: "<<gRefMultiplicity<<endl;\r
+ //if(fDebugLevel) \r
+ //Printf("N(min): %.0f, N(max): %.0f - N(ref): %.0f",fNumberOfAcceptedTracksMin,\r
+ //fNumberOfAcceptedTracksMax,gRefMultiplicity);\r
\r
- if((gRefMultiplicity > fNumberOfAcceptedTracksMin)||(gRefMultiplicity < fNumberOfAcceptedTracksMax)) {\r
+ if((gRefMultiplicity > fNumberOfAcceptedTracksMin) && (gRefMultiplicity < fNumberOfAcceptedTracksMax)) {\r
fHistEventStats->Fill(5,gRefMultiplicity); //events with correct multiplicity\r
return gRefMultiplicity;\r
}\r
// Fills Event statistics histograms\r
\r
Float_t gCentrality = -1.;\r
- Double_t gMultiplicity = 0.;\r
+ Double_t gMultiplicity = -1.;\r
TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
\r
- if(fEventClass == "Centrality"){\r
- if(gAnalysisLevel == "AOD"|| gAnalysisLevel == "MCAOD" || gAnalysisLevel == "MCAODrec" ) { //centrality in AOD header //++++++++++++++\r
- AliAODHeader *header = (AliAODHeader*) event->GetHeader();\r
- if(header){\r
- gCentrality = header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());\r
- }//AOD header\r
- }//AOD\r
+\r
+ // calculate centrality always (not only in centrality mode)\r
+ if(gAnalysisLevel == "AOD"|| gAnalysisLevel == "MCAOD" || gAnalysisLevel == "MCAODrec" ) { //centrality in AOD header //++++++++++++++\r
+ AliAODHeader *header = (AliAODHeader*) event->GetHeader();\r
+ if(header){\r
+ gCentrality = header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());\r
+\r
+ // QA for centrality estimators\r
+ fHistCentStats->Fill(0.,header->GetCentralityP()->GetCentralityPercentile("V0M"));\r
+ fHistCentStats->Fill(1.,header->GetCentralityP()->GetCentralityPercentile("V0A"));\r
+ fHistCentStats->Fill(2.,header->GetCentralityP()->GetCentralityPercentile("V0C"));\r
+ fHistCentStats->Fill(3.,header->GetCentralityP()->GetCentralityPercentile("FMD"));\r
+ fHistCentStats->Fill(4.,header->GetCentralityP()->GetCentralityPercentile("TRK"));\r
+ fHistCentStats->Fill(5.,header->GetCentralityP()->GetCentralityPercentile("TKL")); \r
+ fHistCentStats->Fill(6.,header->GetCentralityP()->GetCentralityPercentile("CL0"));\r
+ fHistCentStats->Fill(7.,header->GetCentralityP()->GetCentralityPercentile("CL1"));\r
+ fHistCentStats->Fill(8.,header->GetCentralityP()->GetCentralityPercentile("ZNA"));\r
+ fHistCentStats->Fill(9.,header->GetCentralityP()->GetCentralityPercentile("ZPA"));\r
+ fHistCentStats->Fill(10.,header->GetCentralityP()->GetCentralityPercentile("V0MvsFMD"));\r
+ fHistCentStats->Fill(11.,header->GetCentralityP()->GetCentralityPercentile("TKLvsV0M"));\r
+ fHistCentStats->Fill(12.,header->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC"));\r
+ \r
+ // Centrality estimator USED ++++++++++++++++++++++++++++++\r
+ fHistCentStatsUsed->Fill(0.,header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data()));\r
+ \r
+ // centrality QA (V0M)\r
+ fHistV0M->Fill(event->GetVZEROData()->GetMTotV0A(), event->GetVZEROData()->GetMTotV0C());\r
+ \r
+ // centrality QA (reference tracks)\r
+ fHistRefTracks->Fill(0.,header->GetRefMultiplicity());\r
+ fHistRefTracks->Fill(1.,header->GetRefMultiplicityPos());\r
+ fHistRefTracks->Fill(2.,header->GetRefMultiplicityNeg());\r
+ fHistRefTracks->Fill(3.,header->GetTPConlyRefMultiplicity());\r
+ fHistRefTracks->Fill(4.,header->GetNumberOfITSClusters(0));\r
+ fHistRefTracks->Fill(5.,header->GetNumberOfITSClusters(1));\r
+ fHistRefTracks->Fill(6.,header->GetNumberOfITSClusters(2));\r
+ fHistRefTracks->Fill(7.,header->GetNumberOfITSClusters(3));\r
+ fHistRefTracks->Fill(8.,header->GetNumberOfITSClusters(4));\r
+\r
+ }//AOD header\r
+ }//AOD\r
+ \r
+ else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD"){ // centrality class for ESDs or MC-ESDs\r
+ AliCentrality *centrality = event->GetCentrality();\r
+ gCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
+\r
+ // QA for centrality estimators\r
+ fHistCentStats->Fill(0.,centrality->GetCentralityPercentile("V0M"));\r
+ fHistCentStats->Fill(1.,centrality->GetCentralityPercentile("V0A"));\r
+ fHistCentStats->Fill(2.,centrality->GetCentralityPercentile("V0C"));\r
+ fHistCentStats->Fill(3.,centrality->GetCentralityPercentile("FMD"));\r
+ fHistCentStats->Fill(4.,centrality->GetCentralityPercentile("TRK"));\r
+ fHistCentStats->Fill(5.,centrality->GetCentralityPercentile("TKL"));\r
+ fHistCentStats->Fill(6.,centrality->GetCentralityPercentile("CL0"));\r
+ fHistCentStats->Fill(7.,centrality->GetCentralityPercentile("CL1"));\r
+ fHistCentStats->Fill(8.,centrality->GetCentralityPercentile("ZNA"));\r
+ fHistCentStats->Fill(9.,centrality->GetCentralityPercentile("ZPA"));\r
+ fHistCentStats->Fill(10.,centrality->GetCentralityPercentile("V0MvsFMD"));\r
+ fHistCentStats->Fill(11.,centrality->GetCentralityPercentile("TKLvsV0M"));\r
+ fHistCentStats->Fill(12.,centrality->GetCentralityPercentile("ZEMvsZDC"));\r
\r
- else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD"){ // centrality class for ESDs or MC-ESDs\r
- AliCentrality *centrality = event->GetCentrality();\r
- gCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
- }//ESD\r
- else if(gAnalysisLevel == "MC"){\r
- Double_t gImpactParameter = 0.;\r
- AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);\r
- if(gMCEvent){\r
- AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(gMCEvent->GenEventHeader()); \r
- if(headerH){\r
- gImpactParameter = headerH->ImpactParameter();\r
- gCentrality = gImpactParameter;\r
- }//MC header\r
- }//MC event cast\r
- }//MC\r
- else{\r
- gCentrality = -1.;\r
- }\r
- }//End if "Centrality"\r
- if(fEventClass=="Multiplicity"&&gAnalysisLevel == "ESD"){\r
+ // Centrality estimator USED ++++++++++++++++++++++++++++++\r
+ fHistCentStatsUsed->Fill(0.,centrality->GetCentralityPercentile(fCentralityEstimator.Data()));\r
+ \r
+ // centrality QA (V0M)\r
+ fHistV0M->Fill(event->GetVZEROData()->GetMTotV0A(), event->GetVZEROData()->GetMTotV0C());\r
+ }//ESD\r
+\r
+ else if(gAnalysisLevel == "MC"){\r
+ Double_t gImpactParameter = 0.;\r
+ AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);\r
+ if(gMCEvent){\r
+ AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(gMCEvent->GenEventHeader()); \r
+ if(headerH){\r
+ gImpactParameter = headerH->ImpactParameter();\r
+ gCentrality = gImpactParameter;\r
+ }//MC header\r
+ }//MC event cast\r
+ }//MC\r
+\r
+ else{\r
+ gCentrality = -1.;\r
+ }\r
+ \r
+ // calculate reference multiplicity always (not only in multiplicity mode)\r
+ if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD"){\r
AliESDEvent* gESDEvent = dynamic_cast<AliESDEvent*>(event);\r
if(gESDEvent){\r
gMultiplicity = fESDtrackCuts->GetReferenceMultiplicity(gESDEvent, AliESDtrackCuts::kTrackletsITSTPC,0.5);\r
+ fHistMultiplicity->Fill(gMultiplicity);\r
}//AliESDevent cast\r
- }\r
- else if(fEventClass=="Multiplicity"&&gAnalysisLevel == "AOD"){\r
+ }//ESD mode\r
+\r
+ else if(gAnalysisLevel == "AOD"|| gAnalysisLevel == "MCAOD" || gAnalysisLevel == "MCAODrec" ){\r
AliAODHeader *header = (AliAODHeader*) event->GetHeader();\r
if ((fMultiplicityEstimator == "V0M")||\r
(fMultiplicityEstimator == "V0A")||\r
(fMultiplicityEstimator == "V0C") ||\r
- (fMultiplicityEstimator == "TPC"))\r
+ (fMultiplicityEstimator == "TPC")) {\r
gMultiplicity = GetReferenceMultiplicityFromAOD(event);\r
+ if(fDebugLevel) Printf("Reference multiplicity (calculated): %.0f",gMultiplicity);\r
+ }\r
else {\r
if(header)\r
gMultiplicity = header->GetRefMultiplicity();\r
+ if(fDebugLevel) Printf("Reference multiplicity (AOD header): %.0f",gMultiplicity);\r
}\r
- }\r
- else if((fEventClass=="Multiplicity")&&(gAnalysisLevel == "MC")) {\r
+ fHistMultiplicity->Fill(gMultiplicity);\r
+ }//AOD mode\r
+ else if(gAnalysisLevel == "MC") {\r
AliMCEvent* gMCEvent = dynamic_cast<AliMCEvent*>(event);\r
//Calculating the multiplicity as the number of charged primaries\r
//within \pm 0.8 in eta and pT > 0.1 GeV/c\r
\r
//exclude non stable particles\r
if(!(gMCEvent->IsPhysicalPrimary(iParticle))) continue;\r
- if(track->Pt() < 0.1) continue;\r
-\r
+ \r
//++++++++++++++++\r
if (fMultiplicityEstimator == "V0M") {\r
if((track->Eta() > 5.1 || track->Eta() < 2.8)&&(track->Eta() < -3.7 || track->Eta() > -1.7)) \r
- continue;}\r
+ continue;}\r
else if (fMultiplicityEstimator == "V0A") {\r
if(track->Eta() > 5.1 || track->Eta() < 2.8) continue;}\r
else if (fMultiplicityEstimator == "V0C") {\r
if(track->Eta() > -1.7 || track->Eta() < -3.7) continue;}\r
else if (fMultiplicityEstimator == "TPC") {\r
- if(track->Eta() > 0.9 || track->Eta() < -0.9) continue;}\r
+ if(track->Eta() < fEtaMin || track->Eta() > fEtaMax) continue;\r
+ if(track->Pt() < fPtMin || track->Pt() > fPtMax) continue;\r
+ }\r
else{\r
+ if(track->Pt() < fPtMin || track->Pt() > fPtMax) continue;\r
if(track->Eta() < fEtaMin || track->Eta() > fEtaMax) continue;\r
}\r
//++++++++++++++++\r
-\r
+ \r
if(track->Charge() == 0) continue;\r
-\r
+ \r
gMultiplicity += 1;\r
}//loop over primaries\r
fHistMultiplicity->Fill(gMultiplicity);\r
- }//MC mode & multiplicity class\r
+ }//MC mode\r
+ else{\r
+ gMultiplicity = -1;\r
+ }\r
\r
+\r
+ // decide what should be returned only here\r
Double_t lReturnVal = -100;\r
if(fEventClass=="Multiplicity"){\r
lReturnVal = gMultiplicity;\r
Double_t AliAnalysisTaskBFPsi::GetReferenceMultiplicityFromAOD(AliVEvent *event){\r
//Function that returns the reference multiplicity from AODs (data or reco MC)\r
//Different ref. mult. implemented: V0M, V0A, V0C, TPC\r
- Double_t gRefMultiplicity = 0.;\r
+ Double_t gRefMultiplicity = 0., gRefMultiplicityTPC = 0.;\r
+ Double_t gRefMultiplicityVZERO = 0., gRefMultiplicityVZEROA = 0., gRefMultiplicityVZEROC = 0.;\r
\r
- if(fMultiplicityEstimator == "TPC") {\r
- // Loop over tracks in event\r
- for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {\r
- AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));\r
- if (!aodTrack) {\r
- AliError(Form("Could not receive track %d", iTracks));\r
- continue;\r
- }\r
- \r
- // AOD track cuts\r
- if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue;\r
+ AliAODHeader *header = dynamic_cast<AliAODHeader *>(event->GetHeader());\r
+ if(!header) {\r
+ Printf("ERROR: AOD header not available");\r
+ return -999;\r
+ }\r
+ Int_t gRunNumber = header->GetRunNumber();\r
+\r
+ // Loop over tracks in event\r
+ for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {\r
+ AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));\r
+ if (!aodTrack) {\r
+ AliError(Form("Could not receive track %d", iTracks));\r
+ continue;\r
+ }\r
+ \r
+ // AOD track cuts\r
+ if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue;\r
\r
- if(aodTrack->Charge() == 0) continue;\r
- // Kinematics cuts from ESD track cuts\r
- if( aodTrack->Pt() < fPtMin || aodTrack->Pt() > fPtMax) continue;\r
- if( aodTrack->Eta() < fEtaMin || aodTrack->Eta() > fEtaMax) continue;\r
+ if(aodTrack->Charge() == 0) continue;\r
+ // Kinematics cuts from ESD track cuts\r
+ if( aodTrack->Pt() < fPtMin || aodTrack->Pt() > fPtMax) continue;\r
+ if( aodTrack->Eta() < fEtaMin || aodTrack->Eta() > fEtaMax) continue;\r
+ \r
+ //=================PID (so far only for electron rejection)==========================//\r
+ if(fElectronRejection) {\r
+ // get the electron nsigma\r
+ Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kElectron));\r
\r
- //=================PID (so far only for electron rejection)==========================//\r
- if(fElectronRejection) {\r
- // get the electron nsigma\r
- Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kElectron));\r
- \r
- // check only for given momentum range\r
- if( aodTrack->Pt() > fElectronRejectionMinPt && aodTrack->Pt() < fElectronRejectionMaxPt ){\r
- //look only at electron nsigma\r
- if(!fElectronOnlyRejection) {\r
- //Make the decision based on the n-sigma of electrons\r
- if(nSigma < fElectronRejectionNSigma) continue;\r
- }\r
- else {\r
- Double_t nSigmaPions = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kPion));\r
- Double_t nSigmaKaons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kKaon));\r
- Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kProton));\r
- \r
- //Make the decision based on the n-sigma of electrons exclusively ( = track not in nsigma region for other species)\r
- if(nSigma < fElectronRejectionNSigma\r
- && nSigmaPions > fElectronRejectionNSigma\r
- && nSigmaKaons > fElectronRejectionNSigma\r
- && nSigmaProtons > fElectronRejectionNSigma ) continue;\r
- }\r
+ // check only for given momentum range\r
+ if( aodTrack->Pt() > fElectronRejectionMinPt && aodTrack->Pt() < fElectronRejectionMaxPt ){\r
+ //look only at electron nsigma\r
+ if(!fElectronOnlyRejection) {\r
+ //Make the decision based on the n-sigma of electrons\r
+ if(nSigma < fElectronRejectionNSigma) continue;\r
}\r
- }//electron rejection\r
- \r
- gRefMultiplicity += 1.0;\r
- }// track loop\r
- }// "TPC" multiplicity estimator\r
- else if((fMultiplicityEstimator == "V0M")||\r
- (fMultiplicityEstimator == "V0A")||\r
- (fMultiplicityEstimator == "V0C")) {\r
- //VZERO segmentation in two detectors (0-31: VZERO-C, 32-63: VZERO-A)\r
- for(Int_t i = 0; i < 64; i++) {\r
- fHistVZEROSignal->Fill(i,event->GetVZEROEqMultiplicity(i));\r
-\r
- if(fMultiplicityEstimator == "V0C") {\r
- if(i > 31) continue;}\r
- else if(fMultiplicityEstimator == "V0A") {\r
- if(i < 32) continue;}\r
- \r
- gRefMultiplicity += event->GetVZEROEqMultiplicity(i);\r
- }//loop over PMTs\r
- }\r
+ else {\r
+ Double_t nSigmaPions = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kPion));\r
+ Double_t nSigmaKaons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kKaon));\r
+ Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kProton));\r
+ \r
+ //Make the decision based on the n-sigma of electrons exclusively ( = track not in nsigma region for other species)\r
+ if(nSigma < fElectronRejectionNSigma\r
+ && nSigmaPions > fElectronRejectionNSigma\r
+ && nSigmaKaons > fElectronRejectionNSigma\r
+ && nSigmaProtons > fElectronRejectionNSigma ) continue;\r
+ }\r
+ }\r
+ }//electron rejection\r
+ \r
+ gRefMultiplicityTPC += 1.0;\r
+ }// track loop\r
+ \r
+ //VZERO segmentation in two detectors (0-31: VZERO-C, 32-63: VZERO-A)\r
+ for(Int_t iChannel = 0; iChannel < 64; iChannel++) {\r
+ fHistVZEROSignal->Fill(iChannel,event->GetVZEROEqMultiplicity(iChannel));\r
+ \r
+ if(iChannel < 32) \r
+ gRefMultiplicityVZEROC += event->GetVZEROEqMultiplicity(iChannel);\r
+ else if(iChannel >= 32) \r
+ gRefMultiplicityVZEROA += event->GetVZEROEqMultiplicity(iChannel);\r
+ }//loop over PMTs\r
+ \r
+ //Equalization of gain\r
+ Double_t gFactorA = GetEqualizationFactor(gRunNumber,"A");\r
+ if(gFactorA != 0)\r
+ gRefMultiplicityVZEROA /= gFactorA;\r
+ Double_t gFactorC = GetEqualizationFactor(gRunNumber,"C");\r
+ if(gFactorC != 0)\r
+ gRefMultiplicityVZEROC /= gFactorC;\r
+ if((gFactorA != 0)&&(gFactorC != 0)) \r
+ gRefMultiplicityVZERO = (gRefMultiplicityVZEROA/gFactorA)+(gRefMultiplicityVZEROC/gFactorC);\r
+ \r
+ if(fDebugLevel) \r
+ Printf("VZERO multiplicity: %.0f - TPC multiplicity: %.0f",gRefMultiplicityVZERO,gRefMultiplicityTPC);\r
+\r
+ fHistTPCvsVZEROMultiplicity->Fill(gRefMultiplicityVZERO,gRefMultiplicityTPC);\r
+\r
+ if(fMultiplicityEstimator == "TPC") \r
+ gRefMultiplicity = gRefMultiplicityTPC;\r
+ else if(fMultiplicityEstimator == "V0M")\r
+ gRefMultiplicity = gRefMultiplicityVZERO;\r
+ else if(fMultiplicityEstimator == "V0A")\r
+ gRefMultiplicity = gRefMultiplicityVZEROA;\r
+ else if(fMultiplicityEstimator == "V0C")\r
+ gRefMultiplicity = gRefMultiplicityVZEROC;\r
\r
return gRefMultiplicity;\r
}\r
\r
if(gAnalysisLevel == "AOD") { // handling of TPC only tracks different in AOD and ESD\r
// Loop over tracks in event\r
+ \r
for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {\r
AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));\r
if (!aodTrack) {\r
for(Int_t iTrackBit = 0; iTrackBit < 16; iTrackBit++){\r
fHistTrackStats->Fill(iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit));\r
}\r
+\r
if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue;\r
\r
vCharge = aodTrack->Charge();\r
Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kElectron));\r
\r
//Fill QA before the PID\r
- fHistdEdxVsPTPCbeforePID -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());\r
- fHistNSigmaTPCvsPtbeforePID -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma); \r
+ fHistdEdxVsPTPCbeforePIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());\r
+ fHistNSigmaTPCvsPtbeforePIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma); \r
//end of QA-before pid\r
\r
// check only for given momentum range\r
}\r
\r
//Fill QA after the PID\r
- fHistdEdxVsPTPCafterPID -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());\r
- fHistNSigmaTPCvsPtafterPID -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma); \r
+ fHistdEdxVsPTPCafterPIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());\r
+ fHistNSigmaTPCvsPtafterPIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma); \r
\r
}\r
//===========================end of PID (so far only for electron rejection)===============================//\r
\r
+ //+++++++++++++++++++++++++++++//\r
+ //===========================PID===============================// \r
+ if(fUsePID) {\r
+ Double_t prob[AliPID::kSPECIES]={0.};\r
+ Double_t probTPC[AliPID::kSPECIES]={0.};\r
+ Double_t probTOF[AliPID::kSPECIES]={0.};\r
+ Double_t probTPCTOF[AliPID::kSPECIES]={0.};\r
+ \r
+ Double_t nSigma = 0.;\r
+ Double_t nSigmaTPC = 0.; //++++\r
+ Double_t nSigmaTOF = 0.; //++++\r
+ Double_t nSigmaTPCTOF = 0.; //++++\r
+ UInt_t detUsedTPC = 0;\r
+ UInt_t detUsedTOF = 0;\r
+ UInt_t detUsedTPCTOF = 0;\r
+ \r
+ nSigmaTPC = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)fParticleOfInterest));\r
+ nSigmaTOF = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(aodTrack,(AliPID::EParticleType)fParticleOfInterest));\r
+ nSigmaTPCTOF = TMath::Sqrt(nSigmaTPC*nSigmaTPC + nSigmaTOF*nSigmaTOF);\r
+\r
+ //Decide what detector configuration we want to use\r
+ switch(fPidDetectorConfig) {\r
+ case kTPCpid:\r
+ fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);\r
+ nSigma = nSigmaTPC;\r
+ detUsedTPC = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTPC);\r
+ for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
+ prob[iSpecies] = probTPC[iSpecies];\r
+ break;\r
+ case kTOFpid:\r
+ fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);\r
+ nSigma = nSigmaTOF;\r
+ detUsedTOF = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTOF);\r
+ for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
+ prob[iSpecies] = probTOF[iSpecies];\r
+ break;\r
+ case kTPCTOF:\r
+ fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);\r
+ nSigma = nSigmaTPCTOF;\r
+ detUsedTPCTOF = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTPCTOF);\r
+ for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
+ prob[iSpecies] = probTPCTOF[iSpecies];\r
+ break;\r
+ default:\r
+ break;\r
+ }//end switch: define detector mask\r
+ \r
+ //Filling the PID QA\r
+ Double_t tofTime = -999., length = 999., tof = -999.;\r
+ Double_t c = TMath::C()*1.E-9;// m/ns\r
+ Double_t beta = -999.;\r
+ if ( (aodTrack->IsOn(AliAODTrack::kTOFin)) &&\r
+ (aodTrack->IsOn(AliAODTrack::kTIME)) ) { \r
+ tofTime = aodTrack->GetTOFsignal();//in ps\r
+ length = aodTrack->GetIntegratedLength();\r
+ tof = tofTime*1E-3; // ns \r
+ \r
+ if (tof <= 0) {\r
+ //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n");\r
+ continue;\r
+ }\r
+ if (length <= 0){\r
+ //printf("WARNING: track with negative length found!Skipping this track for PID checks\n");\r
+ continue;\r
+ }\r
+ \r
+ length = length*0.01; // in meters\r
+ tof = tof*c;\r
+ beta = length/tof;\r
+ \r
+ fHistBetavsPTOFbeforePID ->Fill(aodTrack->P()*aodTrack->Charge(),beta);\r
+ fHistProbTOFvsPtbeforePID ->Fill(aodTrack->Pt(),probTOF[fParticleOfInterest]);\r
+ fHistNSigmaTOFvsPtbeforePID ->Fill(aodTrack->Pt(),nSigmaTOF);\r
+ }//TOF signal \r
+ \r
+ fHistdEdxVsPTPCbeforePID -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());\r
+ fHistProbTPCvsPtbeforePID -> Fill(aodTrack->Pt(),probTPC[fParticleOfInterest]); \r
+ fHistNSigmaTPCvsPtbeforePID -> Fill(aodTrack->Pt(),nSigmaTPC);\r
+ \r
+ fHistProbTPCTOFvsPtbeforePID -> Fill(aodTrack->Pt(),probTPCTOF[fParticleOfInterest]);\r
+\r
+ //combined TPC&TOF\r
+ fHistBetaVsdEdXbeforePID->Fill(aodTrack->GetTPCsignal(),beta); //+++++++++ \r
+ fHistNSigmaTPCTOFvsPtbeforePID -> Fill(aodTrack->Pt(),nSigmaTPCTOF);\r
+ Printf("NSIGMA: %lf - nSigmaTOF: %lf - nSigmaTPC: %lf - nSigmaTPCTOF: %lf",nSigma,nSigmaTOF,nSigmaTPC,nSigmaTPCTOF);\r
+ \r
+ //end of QA-before pid\r
+ \r
+ if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) {\r
+ //Make the decision based on the n-sigma\r
+ if(fUsePIDnSigma) {\r
+ if(nSigma > fPIDNSigma) continue; \r
+ \r
+ fHistNSigmaTOFvsPtafterPID ->Fill(aodTrack->Pt(),nSigmaTOF);\r
+ fHistNSigmaTPCvsPtafterPID ->Fill(aodTrack->Pt(),nSigmaTPC);\r
+ fHistNSigmaTPCTOFvsPtafterPID ->Fill(aodTrack->Pt(),nSigmaTPCTOF);\r
+ }\r
+ //Make the decision based on the bayesian\r
+ else if(fUsePIDPropabilities) {\r
+ if(fParticleOfInterest != TMath::LocMax(AliPID::kSPECIES,prob)) continue;\r
+ if (prob[fParticleOfInterest] < fMinAcceptedPIDProbability) continue; \r
+ \r
+ fHistProbTOFvsPtafterPID ->Fill(aodTrack->Pt(),probTOF[fParticleOfInterest]);\r
+ fHistProbTPCvsPtafterPID ->Fill(aodTrack->Pt(),probTPC[fParticleOfInterest]); \r
+ fHistProbTPCTOFvsPtafterPID ->Fill(aodTrack->Pt(),probTPCTOF[fParticleOfInterest]);\r
+ \r
+ }\r
+ \r
+ //Fill QA after the PID\r
+ fHistBetavsPTOFafterPID ->Fill(aodTrack->P()*aodTrack->Charge(),beta);\r
+ fHistdEdxVsPTPCafterPID ->Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());\r
+ fHistBetaVsdEdXafterPID->Fill(aodTrack->GetTPCsignal(),beta); //+++++++++ \r
+ }\r
+ }\r
+ //===========================PID===============================//\r
+ //+++++++++++++++++++++++++++++//\r
+\r
+\r
Float_t dcaXY = aodTrack->DCA(); // this is the DCA from global track (not exactly what is cut on)\r
Float_t dcaZ = aodTrack->ZAtDCA(); // this is the DCA from global track (not exactly what is cut on)\r
\r
Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kElectron));\r
\r
//Fill QA before the PID\r
- fHistdEdxVsPTPCbeforePID -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());\r
- fHistNSigmaTPCvsPtbeforePID -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma); \r
+ fHistdEdxVsPTPCbeforePIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());\r
+ fHistNSigmaTPCvsPtbeforePIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma); \r
//end of QA-before pid\r
\r
// check only for given momentum range\r
}\r
\r
//Fill QA after the PID\r
- fHistdEdxVsPTPCafterPID -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());\r
- fHistNSigmaTPCvsPtafterPID -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma); \r
+ fHistdEdxVsPTPCafterPIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());\r
+ fHistNSigmaTPCvsPtafterPIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma); \r
\r
}\r
//===========================end of PID (so far only for electron rejection)===============================//\r
return tracksShuffled;\r
}\r
\r
+//________________________________________________________________________\r
+void AliAnalysisTaskBFPsi::SetVZEROCalibrationFile(const char* filename,\r
+ const char* lhcPeriod) {\r
+ //Function to setup the VZERO gain equalization\r
+ //============Get the equilization map============//\r
+ TFile *calibrationFile = TFile::Open(filename);\r
+ if((!calibrationFile)||(!calibrationFile->IsOpen())) {\r
+ Printf("No calibration file found!!!");\r
+ return;\r
+ }\r
+\r
+ TList *list = dynamic_cast<TList *>(calibrationFile->Get(lhcPeriod));\r
+ if(!list) {\r
+ Printf("Calibration TList not found!!!");\r
+ return;\r
+ }\r
+\r
+ fHistVZEROAGainEqualizationMap = dynamic_cast<TH1F *>(list->FindObject("gHistVZEROAGainEqualizationMap"));\r
+ if(!fHistVZEROAGainEqualizationMap) {\r
+ Printf("VZERO-A calibration object not found!!!");\r
+ return;\r
+ }\r
+ fHistVZEROCGainEqualizationMap = dynamic_cast<TH1F *>(list->FindObject("gHistVZEROCGainEqualizationMap"));\r
+ if(!fHistVZEROCGainEqualizationMap) {\r
+ Printf("VZERO-C calibration object not found!!!");\r
+ return;\r
+ }\r
+\r
+ fHistVZEROChannelGainEqualizationMap = dynamic_cast<TH2F *>(list->FindObject("gHistVZEROChannelGainEqualizationMap"));\r
+ if(!fHistVZEROChannelGainEqualizationMap) {\r
+ Printf("VZERO channel calibration object not found!!!");\r
+ return;\r
+ }\r
+}\r
+\r
+//________________________________________________________________________\r
+Double_t AliAnalysisTaskBFPsi::GetChannelEqualizationFactor(Int_t run, \r
+ Int_t channel) {\r
+ //\r
+ if(!fHistVZEROAGainEqualizationMap) return 1.0;\r
+\r
+ for(Int_t iBinX = 1; iBinX <= fHistVZEROChannelGainEqualizationMap->GetNbinsX(); iBinX++) {\r
+ Int_t gRunNumber = atoi(fHistVZEROChannelGainEqualizationMap->GetXaxis()->GetBinLabel(iBinX));\r
+ if(gRunNumber == run)\r
+ return fHistVZEROChannelGainEqualizationMap->GetBinContent(iBinX,channel+1);\r
+ }\r
+\r
+ return 1.0;\r
+}\r
+\r
+//________________________________________________________________________\r
+Double_t AliAnalysisTaskBFPsi::GetEqualizationFactor(Int_t run, \r
+ const char* side) {\r
+ //\r
+ if(!fHistVZEROAGainEqualizationMap) return 1.0;\r
+\r
+ TString gVZEROSide = side;\r
+ for(Int_t iBinX = 1; iBinX < fHistVZEROAGainEqualizationMap->GetNbinsX(); iBinX++) {\r
+ Int_t gRunNumber = atoi(fHistVZEROAGainEqualizationMap->GetXaxis()->GetBinLabel(iBinX));\r
+ //cout<<"Looking for run "<<run<<" - current run: "<<gRunNumber<<endl;\r
+ if(gRunNumber == run) {\r
+ if(gVZEROSide == "A") \r
+ return fHistVZEROAGainEqualizationMap->GetBinContent(iBinX);\r
+ else if(gVZEROSide == "C") \r
+ return fHistVZEROCGainEqualizationMap->GetBinContent(iBinX);\r
+ }\r
+ }\r
+\r
+ return 1.0;\r
+}\r
\r
//________________________________________________________________________\r
void AliAnalysisTaskBFPsi::FinishTaskOutput(){\r