[u/mrichter/AliRoot.git] / PWG2 / AliAnalysisTaskProtons.cxx

#include "TChain.h"
#include "TTree.h"
#include "TString.h"
#include "TList.h"
#include "TH2F.h"
#include "TH1I.h"
#include "TF1.h"
#include "TCanvas.h"
#include "TLorentzVector.h"

#include "AliAnalysisTask.h"
#include "AliAnalysisManager.h"

#include "AliESDEvent.h"
#include "AliESDInputHandler.h"
#include "AliAODEvent.h"
#include "AliAODInputHandler.h"
#include "AliMCEventHandler.h"
#include "AliMCEvent.h"
#include "AliStack.h"
#include "AliCFContainer.h"
#include "AliVertexerTracks.h"
#include "AliESDVertex.h"

#include "AliProtonAnalysis.h"
#include "AliAnalysisTaskProtons.h"

// Analysis task creating a the 2d y-p_t spectrum of p and antip
// Author: Panos Cristakoglou

ClassImp(AliAnalysisTaskProtons)

//________________________________________________________________________ 
AliAnalysisTaskProtons::AliAnalysisTaskProtons()
  : AliAnalysisTask(), fESD(0), fAOD(0), fMC(0), fAnalysisType("ESD"),
    fList(0), fProtonAnalysis(0),
    fElectronFunction(0), fMuonFunction(0),
    fPionFunction(0), fKaonFunction(0), fProtonFunction(0),
    fFunctionUsed(kFALSE) {
  //Dummy constructor
                                                                                                   
}

//________________________________________________________________________
AliAnalysisTaskProtons::AliAnalysisTaskProtons(const char *name) 
: AliAnalysisTask(name, ""), fESD(0), fAOD(0), fMC(0), fAnalysisType("ESD"),
  fList(0), fProtonAnalysis(0), 
  fElectronFunction(0), fMuonFunction(0),
  fPionFunction(0), fKaonFunction(0), fProtonFunction(0),
  fFunctionUsed(kFALSE),
  fTriggerMode(kMB2), fProtonAnalysisMode(kTPC),
  fVxMax(0), fVyMax(0), fVzMax(0) { 
  // Constructor

  // Define input and output slots here
  // Input slot #0 works with a TChain
  DefineInput(0, TChain::Class());
  // Output slot #0 writes into a TList container
  DefineOutput(0, TList::Class());
}

//________________________________________________________________________
void AliAnalysisTaskProtons::ConnectInputData(Option_t *) {
  // Connect ESD or AOD here
  // Called once

  TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
  if (!tree) {
    Printf("ERROR: Could not read chain from input slot 0");
  } else {
    if(fAnalysisType == "ESD") {
      AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
     
      if (!esdH) {
	Printf("ERROR: Could not get ESDInputHandler");
      } else
	fESD = esdH->GetEvent();
    }
    else if(fAnalysisType == "AOD") {
     AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
     
     if (!aodH) {
       Printf("ERROR: Could not get AODInputHandler");
     } else
       fAOD = aodH->GetEvent();
    }
    else if(fAnalysisType == "MC") {
     AliMCEventHandler* mcH = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
     if (!mcH) {
       Printf("ERROR: Could not retrieve MC event handler");
     }
     else
       fMC = mcH->MCEvent();
    }
    else
      Printf("Wrong analysis type: Only ESD, AOD and MC types are allowed!");
  }
}

//________________________________________________________________________
void AliAnalysisTaskProtons::CreateOutputObjects() {
  // Create histograms
  // Called once
  //Prior probabilities
  Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
  
  //proton analysis object
  fProtonAnalysis = new AliProtonAnalysis();

  if(fAnalysisType == "ESD") {
    //Use of TPConly tracks
    if(fProtonAnalysisMode == kTPC) {
      fProtonAnalysis->InitAnalysisHistograms(10, -0.5, 0.5, 16, 0.5, 0.9);
      fProtonAnalysis->UseTPCOnly();
      fProtonAnalysis->SetTPCpid();
      fProtonAnalysis->SetMinTPCClusters(100);
      fProtonAnalysis->SetMaxChi2PerTPCCluster(2.2);
      fProtonAnalysis->SetMaxCov11(0.5);
      fProtonAnalysis->SetMaxCov22(0.5);
      fProtonAnalysis->SetMaxCov33(0.5);
      fProtonAnalysis->SetMaxCov44(0.5);
      fProtonAnalysis->SetMaxCov55(0.5);
      fProtonAnalysis->SetMaxSigmaToVertexTPC(2.0);
      //fProtonAnalysis->SetMaxDCAXYTPC(1.5);
      //fProtonAnalysis->SetMaxDCAZTPC(1.5);
    }
    //Use of HybridTPC tracks
    else if(fProtonAnalysisMode == kHybrid) {
      fProtonAnalysis->InitAnalysisHistograms(10, -0.5, 0.5, 16, 0.5, 0.9);
      fProtonAnalysis->UseHybridTPC();
      fProtonAnalysis->SetTPCpid();
      fProtonAnalysis->SetMinTPCClusters(110);
      fProtonAnalysis->SetMaxChi2PerTPCCluster(2.2);
      fProtonAnalysis->SetMaxCov11(0.5);
      fProtonAnalysis->SetMaxCov22(0.5);
      fProtonAnalysis->SetMaxCov33(0.5);
      fProtonAnalysis->SetMaxCov44(0.5);
      fProtonAnalysis->SetMaxCov55(0.5);
      fProtonAnalysis->SetMaxSigmaToVertex(2.0);
      /*fProtonAnalysis->SetMaxDCAXY(1.5);
	fProtonAnalysis->SetMaxDCAZ(1.5);*/
      fProtonAnalysis->SetPointOnITSLayer6();
      fProtonAnalysis->SetPointOnITSLayer5();
      //fProtonAnalysis->SetPointOnITSLayer4();
      //fProtonAnalysis->SetPointOnITSLayer3();
      fProtonAnalysis->SetPointOnITSLayer2();
      fProtonAnalysis->SetPointOnITSLayer1();
      fProtonAnalysis->SetMinITSClusters(5);
    }
    //Combined tracking
    else if(fProtonAnalysisMode == kGlobal) {
      fProtonAnalysis->InitAnalysisHistograms(10, -0.5, 0.5, 16, 0.5, 0.9);
      fProtonAnalysis->SetMinTPCClusters(110);
      fProtonAnalysis->SetMaxChi2PerTPCCluster(2.2);
      fProtonAnalysis->SetMaxCov11(0.5);
      fProtonAnalysis->SetMaxCov22(0.5);
      fProtonAnalysis->SetMaxCov33(0.5);
      fProtonAnalysis->SetMaxCov44(0.5);
      fProtonAnalysis->SetMaxCov55(0.5);
      fProtonAnalysis->SetMaxSigmaToVertex(2.0);
      //fProtonAnalysis->SetMaxDCAXY(2.0);
      //fProtonAnalysis->SetMaxDCAZ(2.0);
      fProtonAnalysis->SetTPCRefit();
      fProtonAnalysis->SetPointOnITSLayer1();
      fProtonAnalysis->SetPointOnITSLayer2();
      //fProtonAnalysis->SetPointOnITSLayer3();
      //fProtonAnalysis->SetPointOnITSLayer4();
      fProtonAnalysis->SetPointOnITSLayer5();
      fProtonAnalysis->SetPointOnITSLayer6();
      fProtonAnalysis->SetMinITSClusters(5);
      fProtonAnalysis->SetITSRefit();
      fProtonAnalysis->SetESDpid();
    }

    if(fFunctionUsed)
      fProtonAnalysis->SetPriorProbabilityFunctions(fElectronFunction,
					      fMuonFunction,
					      fPionFunction,
					      fKaonFunction,
					      fProtonFunction);
    else
      fProtonAnalysis->SetPriorProbabilities(partFrac);
  }//ESD analysis

  fList = new TList();
  fList->Add(fProtonAnalysis->GetProtonYPtHistogram());
  fList->Add(fProtonAnalysis->GetAntiProtonYPtHistogram());
  fList->Add(fProtonAnalysis->GetEventHistogram());
  fList->Add(fProtonAnalysis->GetProtonContainer());
  fList->Add(fProtonAnalysis->GetAntiProtonContainer());
}

//________________________________________________________________________
void AliAnalysisTaskProtons::Exec(Option_t *) {
  // Main loop
  // Called for each event
  
  if(fAnalysisType == "ESD") {
    if (!fESD) {
      Printf("ERROR: fESD not available");
      return;
    }

    if(IsEventTriggered(fESD,fTriggerMode)) {
      const AliESDVertex *vertex = GetVertex(fESD,fProtonAnalysisMode,
					     fVxMax,fVyMax,fVzMax);
      if(vertex) {
	Printf("Proton ESD analysis task: There are %d tracks in this event", fESD->GetNumberOfTracks());
	fProtonAnalysis->Analyze(fESD,vertex);
      }//reconstructed vertex
    }//triggered event
  }//ESD analysis              
  
  else if(fAnalysisType == "AOD") {
    if (!fAOD) {
      Printf("ERROR: fAOD not available");
      return;
    }
    
    Printf("Proton AOD analysis task: There are %d tracks in this event", fAOD->GetNumberOfTracks());
    fProtonAnalysis->Analyze(fAOD);
  }//AOD analysis                                                                                                                                                                
  
  else if(fAnalysisType == "MC") {
    if (!fMC) {
      Printf("ERROR: Could not retrieve MC event");
      return;
    }

    AliStack* stack = fMC->Stack();
    if (!stack) {
      Printf("ERROR: Could not retrieve the stack");
      return;
    }
    Printf("Proton MC analysis task: There are %d primaries in this event", stack->GetNprimary());
    fProtonAnalysis->Analyze(stack);
  }//MC analysis                      

  // Post output data.
  PostData(0, fList);
}      

//________________________________________________________________________
void AliAnalysisTaskProtons::Terminate(Option_t *) {
  // Draw result to the screen
  // Called once at the end of the query
  
  fList = dynamic_cast<TList*> (GetOutputData(0));
  if (!fList) {
    Printf("ERROR: fList not available");
    return;
  }
   
  TH2F *fHistYPtProtons = (TH2F *)fList->At(0);
  TH2F *fHistYPtAntiProtons = (TH2F *)fList->At(1);
    
  TCanvas *c1 = new TCanvas("c1","p-\bar{p}",200,0,800,400);
  c1->SetFillColor(10); c1->SetHighLightColor(10); c1->Divide(2,1);

  c1->cd(1)->SetLeftMargin(0.15); c1->cd(1)->SetBottomMargin(0.15);  
  if (fHistYPtProtons) fHistYPtProtons->DrawCopy("colz");
  c1->cd(2)->SetLeftMargin(0.15); c1->cd(2)->SetBottomMargin(0.15);  
  if (fHistYPtAntiProtons) fHistYPtAntiProtons->DrawCopy("colz");
}

//________________________________________________________________________
Bool_t AliAnalysisTaskProtons::IsEventTriggered(const AliESDEvent *esd, 
						TriggerMode trigger) {
  // check if the event was triggered
  ULong64_t triggerMask = esd->GetTriggerMask();
  
  // definitions from p-p.cfg
  ULong64_t spdFO = (1 << 14);
  ULong64_t v0left = (1 << 11);
  ULong64_t v0right = (1 << 12);

  switch (trigger) {
  case kMB1: {
    if (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right)))
      return kTRUE;
    break;
  }
  case kMB2: {
    if (triggerMask & spdFO && ((triggerMask & v0left) || (triggerMask & v0right)))
      return kTRUE;
    break;
  }
  case kSPDFASTOR: {
    if (triggerMask & spdFO)
      return kTRUE;
    break;
  }
  }//switch

  return kFALSE;
}

//________________________________________________________________________
const AliESDVertex* AliAnalysisTaskProtons::GetVertex(AliESDEvent* esd, 
						      AnalysisMode mode,
						      Double_t gVxMax,
						      Double_t gVyMax,
						      Double_t gVzMax) {
  // Get the vertex from the ESD and returns it if the vertex is valid
  // Second argument decides which vertex is used (this selects
  // also the quality criteria that are applied)
  const AliESDVertex* vertex = 0;
  if(mode == kHybrid)
    vertex = esd->GetPrimaryVertexSPD();
  else if(mode == kTPC){
    Double_t kBz = esd->GetMagneticField();
    AliVertexerTracks vertexer(kBz);
    vertexer.SetTPCMode();
    AliESDVertex *vTPC = vertexer.FindPrimaryVertex(esd);
    esd->SetPrimaryVertexTPC(vTPC);
    for (Int_t i=0; i<esd->GetNumberOfTracks(); i++) {
      AliESDtrack *t = esd->GetTrack(i);
      t->RelateToVertexTPC(vTPC, kBz, kVeryBig);
    }
    delete vTPC;
    vertex = esd->GetPrimaryVertexTPC();
  }
  else if(mode == kGlobal)
    vertex = esd->GetPrimaryVertex();
  else
    Printf("GetVertex: ERROR: Invalid second argument %d", mode);
  
  if(!vertex) return 0;
  
  // check Ncontributors
  if(vertex->GetNContributors() <= 0) return 0;
  
  // check resolution
  Double_t zRes = vertex->GetZRes();
  if(zRes == 0) return 0;
  
  //check position
  if(TMath::Abs(vertex->GetXv()) > gVxMax) return 0;
  if(TMath::Abs(vertex->GetYv()) > gVyMax) return 0;
  if(TMath::Abs(vertex->GetZv()) > gVzMax) return 0;
  
  return vertex;
}
Commit	Line	Data
734d2c12	1	#include "TChain.h"
734d2c12	2	#include "TTree.h"
df01f00b	3	#include "TString.h"
734d2c12	4	#include "TList.h"
734d2c12	5	#include "TH2F.h"
df01f00b	6	#include "TH1I.h"
aafecd8b	7	#include "TF1.h"
734d2c12	8	#include "TCanvas.h"
	9	#include "TLorentzVector.h"
	10
	11	#include "AliAnalysisTask.h"
	12	#include "AliAnalysisManager.h"
	13
	14	#include "AliESDEvent.h"
	15	#include "AliESDInputHandler.h"
b620b667	16	#include "AliAODEvent.h"
b620b667	17	#include "AliAODInputHandler.h"
e4358d7f	18	#include "AliMCEventHandler.h"
	19	#include "AliMCEvent.h"
	20	#include "AliStack.h"
251e4034	21	#include "AliCFContainer.h"
6667f3a7	22	#include "AliVertexerTracks.h"
6667f3a7	23	#include "AliESDVertex.h"
734d2c12	24
5ae0977e	25	#include "AliProtonAnalysis.h"
734d2c12	26	#include "AliAnalysisTaskProtons.h"
	27
	28	// Analysis task creating a the 2d y-p_t spectrum of p and antip
	29	// Author: Panos Cristakoglou
	30
	31	ClassImp(AliAnalysisTaskProtons)
	32
df01f00b	33	//________________________________________________________________________
	34	AliAnalysisTaskProtons::AliAnalysisTaskProtons()
	35	: AliAnalysisTask(), fESD(0), fAOD(0), fMC(0), fAnalysisType("ESD"),
6667f3a7	36	fList(0), fProtonAnalysis(0),
df01f00b	37	fElectronFunction(0), fMuonFunction(0),
	38	fPionFunction(0), fKaonFunction(0), fProtonFunction(0),
	39	fFunctionUsed(kFALSE) {
db10bcb0	40	//Dummy constructor
df01f00b	41
db10bcb0	42	}
db10bcb0	43
734d2c12	44	//________________________________________________________________________
734d2c12	45	AliAnalysisTaskProtons::AliAnalysisTaskProtons(const char *name)
df01f00b	46	: AliAnalysisTask(name, ""), fESD(0), fAOD(0), fMC(0), fAnalysisType("ESD"),
6667f3a7	47	fList(0), fProtonAnalysis(0),
aafecd8b	48	fElectronFunction(0), fMuonFunction(0),
aafecd8b	49	fPionFunction(0), fKaonFunction(0), fProtonFunction(0),
6667f3a7	50	fFunctionUsed(kFALSE),
	51	fTriggerMode(kMB2), fProtonAnalysisMode(kTPC),
	52	fVxMax(0), fVyMax(0), fVzMax(0) {
734d2c12	53	// Constructor
	54
	55	// Define input and output slots here
	56	// Input slot #0 works with a TChain
	57	DefineInput(0, TChain::Class());
	58	// Output slot #0 writes into a TList container
	59	DefineOutput(0, TList::Class());
	60	}
	61
	62	//________________________________________________________________________
	63	void AliAnalysisTaskProtons::ConnectInputData(Option_t *) {
	64	// Connect ESD or AOD here
	65	// Called once
	66
	67	TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
	68	if (!tree) {
	69	Printf("ERROR: Could not read chain from input slot 0");
	70	} else {
b620b667	71	if(fAnalysisType == "ESD") {
df01f00b	72	AliESDInputHandler esdH = dynamic_cast<AliESDInputHandler> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
df01f00b	73
b620b667	74	if (!esdH) {
	75	Printf("ERROR: Could not get ESDInputHandler");
	76	} else
	77	fESD = esdH->GetEvent();
	78	}
	79	else if(fAnalysisType == "AOD") {
df01f00b	80	AliAODInputHandler aodH = dynamic_cast<AliAODInputHandler> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
	81
	82	if (!aodH) {
	83	Printf("ERROR: Could not get AODInputHandler");
	84	} else
	85	fAOD = aodH->GetEvent();
b620b667	86	}
e4358d7f	87	else if(fAnalysisType == "MC") {
df01f00b	88	AliMCEventHandler* mcH = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
	89	if (!mcH) {
	90	Printf("ERROR: Could not retrieve MC event handler");
	91	}
	92	else
	93	fMC = mcH->MCEvent();
e4358d7f	94	}
df01f00b	95	else
e4358d7f	96	Printf("Wrong analysis type: Only ESD, AOD and MC types are allowed!");
734d2c12	97	}
	98	}
	99
	100	//________________________________________________________________________
	101	void AliAnalysisTaskProtons::CreateOutputObjects() {
	102	// Create histograms
	103	// Called once
df01f00b	104	//Prior probabilities
734d2c12	105	Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
df01f00b	106
df01f00b	107	//proton analysis object
6667f3a7	108	fProtonAnalysis = new AliProtonAnalysis();
df01f00b	109
b620b667	110	if(fAnalysisType == "ESD") {
251e4034	111	//Use of TPConly tracks
6667f3a7	112	if(fProtonAnalysisMode == kTPC) {
	113	fProtonAnalysis->InitAnalysisHistograms(10, -0.5, 0.5, 16, 0.5, 0.9);
	114	fProtonAnalysis->UseTPCOnly();
	115	fProtonAnalysis->SetTPCpid();
	116	fProtonAnalysis->SetMinTPCClusters(100);
	117	fProtonAnalysis->SetMaxChi2PerTPCCluster(2.2);
	118	fProtonAnalysis->SetMaxCov11(0.5);
	119	fProtonAnalysis->SetMaxCov22(0.5);
	120	fProtonAnalysis->SetMaxCov33(0.5);
	121	fProtonAnalysis->SetMaxCov44(0.5);
	122	fProtonAnalysis->SetMaxCov55(0.5);
	123	fProtonAnalysis->SetMaxSigmaToVertexTPC(2.0);
	124	//fProtonAnalysis->SetMaxDCAXYTPC(1.5);
	125	//fProtonAnalysis->SetMaxDCAZTPC(1.5);
	126	}
	127	//Use of HybridTPC tracks
	128	else if(fProtonAnalysisMode == kHybrid) {
	129	fProtonAnalysis->InitAnalysisHistograms(10, -0.5, 0.5, 16, 0.5, 0.9);
	130	fProtonAnalysis->UseHybridTPC();
	131	fProtonAnalysis->SetTPCpid();
	132	fProtonAnalysis->SetMinTPCClusters(110);
	133	fProtonAnalysis->SetMaxChi2PerTPCCluster(2.2);
	134	fProtonAnalysis->SetMaxCov11(0.5);
	135	fProtonAnalysis->SetMaxCov22(0.5);
	136	fProtonAnalysis->SetMaxCov33(0.5);
	137	fProtonAnalysis->SetMaxCov44(0.5);
	138	fProtonAnalysis->SetMaxCov55(0.5);
	139	fProtonAnalysis->SetMaxSigmaToVertex(2.0);
	140	/*fProtonAnalysis->SetMaxDCAXY(1.5);
	141	fProtonAnalysis->SetMaxDCAZ(1.5);*/
	142	fProtonAnalysis->SetPointOnITSLayer6();
	143	fProtonAnalysis->SetPointOnITSLayer5();
	144	//fProtonAnalysis->SetPointOnITSLayer4();
	145	//fProtonAnalysis->SetPointOnITSLayer3();
	146	fProtonAnalysis->SetPointOnITSLayer2();
	147	fProtonAnalysis->SetPointOnITSLayer1();
	148	fProtonAnalysis->SetMinITSClusters(5);
	149	}
	150	//Combined tracking
	151	else if(fProtonAnalysisMode == kGlobal) {
	152	fProtonAnalysis->InitAnalysisHistograms(10, -0.5, 0.5, 16, 0.5, 0.9);
	153	fProtonAnalysis->SetMinTPCClusters(110);
	154	fProtonAnalysis->SetMaxChi2PerTPCCluster(2.2);
	155	fProtonAnalysis->SetMaxCov11(0.5);
	156	fProtonAnalysis->SetMaxCov22(0.5);
	157	fProtonAnalysis->SetMaxCov33(0.5);
	158	fProtonAnalysis->SetMaxCov44(0.5);
	159	fProtonAnalysis->SetMaxCov55(0.5);
	160	fProtonAnalysis->SetMaxSigmaToVertex(2.0);
	161	//fProtonAnalysis->SetMaxDCAXY(2.0);
	162	//fProtonAnalysis->SetMaxDCAZ(2.0);
	163	fProtonAnalysis->SetTPCRefit();
	164	fProtonAnalysis->SetPointOnITSLayer1();
	165	fProtonAnalysis->SetPointOnITSLayer2();
	166	//fProtonAnalysis->SetPointOnITSLayer3();
	167	//fProtonAnalysis->SetPointOnITSLayer4();
	168	fProtonAnalysis->SetPointOnITSLayer5();
	169	fProtonAnalysis->SetPointOnITSLayer6();
	170	fProtonAnalysis->SetMinITSClusters(5);
	171	fProtonAnalysis->SetITSRefit();
	172	fProtonAnalysis->SetESDpid();
	173	}
df01f00b	174
6667f3a7	175	if(fFunctionUsed)
	176	fProtonAnalysis->SetPriorProbabilityFunctions(fElectronFunction,
	177	fMuonFunction,
	178	fPionFunction,
	179	fKaonFunction,
	180	fProtonFunction);
	181	else
	182	fProtonAnalysis->SetPriorProbabilities(partFrac);
	183	}//ESD analysis
738619fd	184
738619fd	185	fList = new TList();
6667f3a7	186	fList->Add(fProtonAnalysis->GetProtonYPtHistogram());
	187	fList->Add(fProtonAnalysis->GetAntiProtonYPtHistogram());
	188	fList->Add(fProtonAnalysis->GetEventHistogram());
	189	fList->Add(fProtonAnalysis->GetProtonContainer());
	190	fList->Add(fProtonAnalysis->GetAntiProtonContainer());
734d2c12	191	}
	192
	193	//________________________________________________________________________
	194	void AliAnalysisTaskProtons::Exec(Option_t *) {
	195	// Main loop
	196	// Called for each event
	197
b620b667	198	if(fAnalysisType == "ESD") {
	199	if (!fESD) {
	200	Printf("ERROR: fESD not available");
	201	return;
	202	}
df01f00b	203
6667f3a7	204	if(IsEventTriggered(fESD,fTriggerMode)) {
	205	const AliESDVertex *vertex = GetVertex(fESD,fProtonAnalysisMode,
	206	fVxMax,fVyMax,fVzMax);
	207	if(vertex) {
	208	Printf("Proton ESD analysis task: There are %d tracks in this event", fESD->GetNumberOfTracks());
	209	fProtonAnalysis->Analyze(fESD,vertex);
	210	}//reconstructed vertex
	211	}//triggered event
df01f00b	212	}//ESD analysis
df01f00b	213
b620b667	214	else if(fAnalysisType == "AOD") {
	215	if (!fAOD) {
	216	Printf("ERROR: fAOD not available");
	217	return;
	218	}
	219
e4358d7f	220	Printf("Proton AOD analysis task: There are %d tracks in this event", fAOD->GetNumberOfTracks());
6667f3a7	221	fProtonAnalysis->Analyze(fAOD);
df01f00b	222	}//AOD analysis
df01f00b	223
e4358d7f	224	else if(fAnalysisType == "MC") {
	225	if (!fMC) {
	226	Printf("ERROR: Could not retrieve MC event");
	227	return;
	228	}
df01f00b	229
e4358d7f	230	AliStack* stack = fMC->Stack();
	231	if (!stack) {
	232	Printf("ERROR: Could not retrieve the stack");
	233	return;
	234	}
	235	Printf("Proton MC analysis task: There are %d primaries in this event", stack->GetNprimary());
6667f3a7	236	fProtonAnalysis->Analyze(stack);
df01f00b	237	}//MC analysis
734d2c12	238
	239	// Post output data.
	240	PostData(0, fList);
	241	}
	242
	243	//________________________________________________________________________
	244	void AliAnalysisTaskProtons::Terminate(Option_t *) {
	245	// Draw result to the screen
	246	// Called once at the end of the query
	247
	248	fList = dynamic_cast<TList*> (GetOutputData(0));
	249	if (!fList) {
	250	Printf("ERROR: fList not available");
	251	return;
	252	}
	253
	254	TH2F fHistYPtProtons = (TH2F )fList->At(0);
	255	TH2F fHistYPtAntiProtons = (TH2F )fList->At(1);
	256
df01f00b	257	TCanvas *c1 = new TCanvas("c1","p-\bar{p}",200,0,800,400);
df01f00b	258	c1->SetFillColor(10); c1->SetHighLightColor(10); c1->Divide(2,1);
734d2c12	259
df01f00b	260	c1->cd(1)->SetLeftMargin(0.15); c1->cd(1)->SetBottomMargin(0.15);
734d2c12	261	if (fHistYPtProtons) fHistYPtProtons->DrawCopy("colz");
df01f00b	262	c1->cd(2)->SetLeftMargin(0.15); c1->cd(2)->SetBottomMargin(0.15);
734d2c12	263	if (fHistYPtAntiProtons) fHistYPtAntiProtons->DrawCopy("colz");
	264	}
	265
6667f3a7	266	//________________________________________________________________________
	267	Bool_t AliAnalysisTaskProtons::IsEventTriggered(const AliESDEvent *esd,
	268	TriggerMode trigger) {
	269	// check if the event was triggered
	270	ULong64_t triggerMask = esd->GetTriggerMask();
	271
	272	// definitions from p-p.cfg
	273	ULong64_t spdFO = (1 << 14);
	274	ULong64_t v0left = (1 << 11);
	275	ULong64_t v0right = (1 << 12);
	276
	277	switch (trigger) {
	278	case kMB1: {
	279	if (triggerMask & spdFO \|\| ((triggerMask & v0left) \|\| (triggerMask & v0right)))
	280	return kTRUE;
	281	break;
	282	}
	283	case kMB2: {
	284	if (triggerMask & spdFO && ((triggerMask & v0left) \|\| (triggerMask & v0right)))
	285	return kTRUE;
	286	break;
	287	}
	288	case kSPDFASTOR: {
	289	if (triggerMask & spdFO)
	290	return kTRUE;
	291	break;
	292	}
	293	}//switch
	294
	295	return kFALSE;
	296	}
	297
	298	//________________________________________________________________________
	299	const AliESDVertex* AliAnalysisTaskProtons::GetVertex(AliESDEvent* esd,
	300	AnalysisMode mode,
	301	Double_t gVxMax,
	302	Double_t gVyMax,
	303	Double_t gVzMax) {
	304	// Get the vertex from the ESD and returns it if the vertex is valid
	305	// Second argument decides which vertex is used (this selects
	306	// also the quality criteria that are applied)
	307	const AliESDVertex* vertex = 0;
	308	if(mode == kHybrid)
	309	vertex = esd->GetPrimaryVertexSPD();
	310	else if(mode == kTPC){
	311	Double_t kBz = esd->GetMagneticField();
	312	AliVertexerTracks vertexer(kBz);
	313	vertexer.SetTPCMode();
	314	AliESDVertex *vTPC = vertexer.FindPrimaryVertex(esd);
	315	esd->SetPrimaryVertexTPC(vTPC);
	316	for (Int_t i=0; i<esd->GetNumberOfTracks(); i++) {
	317	AliESDtrack *t = esd->GetTrack(i);
	318	t->RelateToVertexTPC(vTPC, kBz, kVeryBig);
	319	}
	320	delete vTPC;
	321	vertex = esd->GetPrimaryVertexTPC();
	322	}
	323	else if(mode == kGlobal)
	324	vertex = esd->GetPrimaryVertex();
	325	else
	326	Printf("GetVertex: ERROR: Invalid second argument %d", mode);
	327
	328	if(!vertex) return 0;
	329
330	// check Ncontributors
331	if(vertex->GetNContributors() <= 0) return 0;
332
333	// check resolution
334	Double_t zRes = vertex->GetZRes();
335	if(zRes == 0) return 0;
336
337	//check position
338	if(TMath::Abs(vertex->GetXv()) > gVxMax) return 0;
339	if(TMath::Abs(vertex->GetYv()) > gVyMax) return 0;
340	if(TMath::Abs(vertex->GetZv()) > gVzMax) return 0;
341
342	return vertex;
343	}
344
734d2c12	345