[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(20, -1.0, 1.0, 48, 0.3, 1.5);
      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
  else if(fAnalysisType == "MC") 
    fProtonAnalysis->InitAnalysisHistograms(10, -0.5, 0.5, 16, 0.5, 0.9);
  
  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,kFALSE);//kTRUE in case of inclusive measurement
  }//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();
da8c4c1c	115	//fProtonAnalysis->SetTPCpid();
6667f3a7	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) {
da8c4c1c	152	fProtonAnalysis->InitAnalysisHistograms(20, -1.0, 1.0, 48, 0.3, 1.5);
6667f3a7	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
da8c4c1c	184	else if(fAnalysisType == "MC")
	185	fProtonAnalysis->InitAnalysisHistograms(10, -0.5, 0.5, 16, 0.5, 0.9);
	186
738619fd	187	fList = new TList();
6667f3a7	188	fList->Add(fProtonAnalysis->GetProtonYPtHistogram());
	189	fList->Add(fProtonAnalysis->GetAntiProtonYPtHistogram());
	190	fList->Add(fProtonAnalysis->GetEventHistogram());
	191	fList->Add(fProtonAnalysis->GetProtonContainer());
	192	fList->Add(fProtonAnalysis->GetAntiProtonContainer());
734d2c12	193	}
	194
	195	//________________________________________________________________________
	196	void AliAnalysisTaskProtons::Exec(Option_t *) {
	197	// Main loop
	198	// Called for each event
	199
b620b667	200	if(fAnalysisType == "ESD") {
	201	if (!fESD) {
	202	Printf("ERROR: fESD not available");
	203	return;
	204	}
df01f00b	205
6667f3a7	206	if(IsEventTriggered(fESD,fTriggerMode)) {
	207	const AliESDVertex *vertex = GetVertex(fESD,fProtonAnalysisMode,
	208	fVxMax,fVyMax,fVzMax);
	209	if(vertex) {
	210	Printf("Proton ESD analysis task: There are %d tracks in this event", fESD->GetNumberOfTracks());
	211	fProtonAnalysis->Analyze(fESD,vertex);
	212	}//reconstructed vertex
	213	}//triggered event
df01f00b	214	}//ESD analysis
df01f00b	215
b620b667	216	else if(fAnalysisType == "AOD") {
	217	if (!fAOD) {
	218	Printf("ERROR: fAOD not available");
	219	return;
	220	}
	221
e4358d7f	222	Printf("Proton AOD analysis task: There are %d tracks in this event", fAOD->GetNumberOfTracks());
6667f3a7	223	fProtonAnalysis->Analyze(fAOD);
df01f00b	224	}//AOD analysis
df01f00b	225
e4358d7f	226	else if(fAnalysisType == "MC") {
	227	if (!fMC) {
	228	Printf("ERROR: Could not retrieve MC event");
	229	return;
	230	}
df01f00b	231
e4358d7f	232	AliStack* stack = fMC->Stack();
	233	if (!stack) {
	234	Printf("ERROR: Could not retrieve the stack");
	235	return;
	236	}
	237	Printf("Proton MC analysis task: There are %d primaries in this event", stack->GetNprimary());
da8c4c1c	238	fProtonAnalysis->Analyze(stack,kFALSE);//kTRUE in case of inclusive measurement
df01f00b	239	}//MC analysis
734d2c12	240
	241	// Post output data.
	242	PostData(0, fList);
	243	}
	244
	245	//________________________________________________________________________
	246	void AliAnalysisTaskProtons::Terminate(Option_t *) {
	247	// Draw result to the screen
	248	// Called once at the end of the query
	249
	250	fList = dynamic_cast<TList*> (GetOutputData(0));
	251	if (!fList) {
	252	Printf("ERROR: fList not available");
	253	return;
	254	}
	255
	256	TH2F fHistYPtProtons = (TH2F )fList->At(0);
	257	TH2F fHistYPtAntiProtons = (TH2F )fList->At(1);
	258
df01f00b	259	TCanvas *c1 = new TCanvas("c1","p-\bar{p}",200,0,800,400);
df01f00b	260	c1->SetFillColor(10); c1->SetHighLightColor(10); c1->Divide(2,1);
734d2c12	261
df01f00b	262	c1->cd(1)->SetLeftMargin(0.15); c1->cd(1)->SetBottomMargin(0.15);
734d2c12	263	if (fHistYPtProtons) fHistYPtProtons->DrawCopy("colz");
df01f00b	264	c1->cd(2)->SetLeftMargin(0.15); c1->cd(2)->SetBottomMargin(0.15);
734d2c12	265	if (fHistYPtAntiProtons) fHistYPtAntiProtons->DrawCopy("colz");
	266	}
	267
6667f3a7	268	//________________________________________________________________________
	269	Bool_t AliAnalysisTaskProtons::IsEventTriggered(const AliESDEvent *esd,
	270	TriggerMode trigger) {
	271	// check if the event was triggered
	272	ULong64_t triggerMask = esd->GetTriggerMask();
	273
	274	// definitions from p-p.cfg
	275	ULong64_t spdFO = (1 << 14);
	276	ULong64_t v0left = (1 << 11);
	277	ULong64_t v0right = (1 << 12);
	278
	279	switch (trigger) {
	280	case kMB1: {
	281	if (triggerMask & spdFO \|\| ((triggerMask & v0left) \|\| (triggerMask & v0right)))
	282	return kTRUE;
	283	break;
	284	}
	285	case kMB2: {
	286	if (triggerMask & spdFO && ((triggerMask & v0left) \|\| (triggerMask & v0right)))
	287	return kTRUE;
	288	break;
	289	}
	290	case kSPDFASTOR: {
	291	if (triggerMask & spdFO)
	292	return kTRUE;
	293	break;
	294	}
	295	}//switch
	296
	297	return kFALSE;
	298	}
	299
	300	//________________________________________________________________________
	301	const AliESDVertex* AliAnalysisTaskProtons::GetVertex(AliESDEvent* esd,
	302	AnalysisMode mode,
	303	Double_t gVxMax,
	304	Double_t gVyMax,
	305	Double_t gVzMax) {
	306	// Get the vertex from the ESD and returns it if the vertex is valid
	307	// Second argument decides which vertex is used (this selects
	308	// also the quality criteria that are applied)
	309	const AliESDVertex* vertex = 0;
	310	if(mode == kHybrid)
	311	vertex = esd->GetPrimaryVertexSPD();
	312	else if(mode == kTPC){
	313	Double_t kBz = esd->GetMagneticField();
	314	AliVertexerTracks vertexer(kBz);
	315	vertexer.SetTPCMode();
	316	AliESDVertex *vTPC = vertexer.FindPrimaryVertex(esd);
	317	esd->SetPrimaryVertexTPC(vTPC);
	318	for (Int_t i=0; i<esd->GetNumberOfTracks(); i++) {
	319	AliESDtrack *t = esd->GetTrack(i);
	320	t->RelateToVertexTPC(vTPC, kBz, kVeryBig);
	321	}
	322	delete vTPC;
	323	vertex = esd->GetPrimaryVertexTPC();
	324	}
	325	else if(mode == kGlobal)
	326	vertex = esd->GetPrimaryVertex();
	327	else
	328	Printf("GetVertex: ERROR: Invalid second argument %d", mode);
	329
	330	if(!vertex) return 0;
	331
332	// check Ncontributors
333	if(vertex->GetNContributors() <= 0) return 0;
334
335	// check resolution
336	Double_t zRes = vertex->GetZRes();
337	if(zRes == 0) return 0;
338
339	//check position
340	if(TMath::Abs(vertex->GetXv()) > gVxMax) return 0;
341	if(TMath::Abs(vertex->GetYv()) > gVyMax) return 0;
342	if(TMath::Abs(vertex->GetZv()) > gVzMax) return 0;
343
344	return vertex;
345	}
346
734d2c12	347