[u/mrichter/AliRoot.git] / CORRFW / test / AliCFV0Task.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

//-----------------------------------------------------------------------
// Example of task running on AliEn (CAF?)
// which provides standard way of calculating acceptance and efficiency
// between different steps of the procedure.
// The ouptut of the task is a AliCFContainer from which the efficiencies
// can be calculated
//-----------------------------------------------------------------------
// Author : R. Vernet, Consorzio Cometa - Catania (it)
//-----------------------------------------------------------------------


#ifndef ALICFV0TASK_CXX
#define ALICFV0TASK_CXX

#include "AliCFV0Task.h"
#include "TCanvas.h"
#include "AliStack.h"
#include "TParticle.h"
#include "TH1I.h"
#include "AliMCEvent.h"
#include "AliCFManager.h"
#include "AliCFContainer.h"
#include "AliESDtrack.h"
#include "AliLog.h"
#include "AliESDv0.h"
#include "AliV0vertexer.h"
#include "AliCFPair.h"
#include "AliESDEvent.h"
#include "AliAODEvent.h"
#include "TChain.h"
#include "AliCFParticleGenCuts.h"
#include "AliAODv0.h"
#include "TDatabasePDG.h"

//__________________________________________________________________________
AliCFV0Task::AliCFV0Task() :
  AliAnalysisTaskSE(),
  fRebuildV0s(0),
  fV0PDG(310),
  fCFManager(0x0),
  fHistEventsProcessed(0x0)
{
  //
  //Default ctor
  //
}
//___________________________________________________________________________
AliCFV0Task::AliCFV0Task(const Char_t* name) :
  AliAnalysisTaskSE(name),
  fRebuildV0s(0),
  fV0PDG(310),
  fCFManager(0x0),
  fHistEventsProcessed(0x0)
{
  //
  // Constructor. Initialization of Inputs and Outputs
  //
  Info("AliCFV0Task","Calling Constructor");
  /*
    DefineInput(0) and DefineOutput(0)
    are taken care of by AliAnalysisTaskSE constructor
  */
  DefineOutput(1,TH1I::Class());
  DefineOutput(2,AliCFContainer::Class());
//   DefineOutput(3,TList::Class());
}

//___________________________________________________________________________
AliCFV0Task& AliCFV0Task::operator=(const AliCFV0Task& c) 
{
  //
  // Assignment operator
  //
  if (this!=&c) {
    AliAnalysisTaskSE::operator=(c) ;
    fRebuildV0s = c.fRebuildV0s;
    fV0PDG      = c.fV0PDG;
    fCFManager  = c.fCFManager;
    fHistEventsProcessed = c.fHistEventsProcessed;
  }
  return *this;
}

//___________________________________________________________________________
AliCFV0Task::AliCFV0Task(const AliCFV0Task& c) :
  AliAnalysisTaskSE(c),
  fRebuildV0s(c.fRebuildV0s),
  fV0PDG(c.fV0PDG),
  fCFManager(c.fCFManager),
  fHistEventsProcessed(c.fHistEventsProcessed)
{
  //
  // Copy Constructor
  //
}

//___________________________________________________________________________
AliCFV0Task::~AliCFV0Task() {
  //
  //destructor
  //
  Info("~AliCFV0Task","Calling Destructor");
  if (fCFManager)           delete fCFManager ;
  if (fHistEventsProcessed) delete fHistEventsProcessed ;
}

//_________________________________________________
void AliCFV0Task::UserExec(Option_t *)
{
  //
  // Main loop function
  //
  Info("UserExec","") ;

  if (!fInputEvent) {
    Error("UserExec","NO EVENT FOUND!");
    return;
  }

  if (!fMCEvent) {
    Error("UserExec","NO MC INFO FOUND!");
    return;
  }

  fCFManager->SetMCEventInfo (fMCEvent);
  fCFManager->SetRecEventInfo(fInputEvent);

  Bool_t isESDEvent = strcmp(fInputEvent->ClassName(),"AliESDEvent") == 0 ? kTRUE : kFALSE ;
  Bool_t isAODEvent = strcmp(fInputEvent->ClassName(),"AliAODEvent") == 0 ? kTRUE : kFALSE ;

  AliESDEvent* fESD = dynamic_cast<AliESDEvent*>(fInputEvent);
  AliAODEvent* fAOD = dynamic_cast<AliAODEvent*>(fInputEvent);
  

  // MC-event selection
  Double_t containerInput[2] ;
        
  //loop on the MC event
  Info("UserExec","Looping on MC event");
  for (Int_t ipart=0; ipart<fMCEvent->GetNumberOfTracks(); ipart++) { 
    AliMCParticle *mcPart  = (AliMCParticle*)fMCEvent->GetTrack(ipart);

    //check the MC-level cuts
    if (!fCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,mcPart)) continue;
    containerInput[0] = mcPart->Pt();
    containerInput[1] = mcPart->Y() ;
    //fill the container for Gen-level selection
    fCFManager->GetParticleContainer()->Fill(containerInput,kStepGenerated);
    
    //check the Acceptance-level cuts
    if (!fCFManager->CheckParticleCuts(AliCFManager::kPartAccCuts,mcPart)) continue;
    //fill the container for Acceptance-level selection
    fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructible);
  }    


  //Now go to rec level
  Info("UserExec","Looping on %s",fInputEvent->ClassName());
 
  if (isESDEvent && fRebuildV0s) RebuildV0s(fESD) ;

  Info("UserExec","There are %d V0s in event",fInputEvent->GetNumberOfV0s());
  
  AliESDv0*  esdV0 = 0x0;
  AliAODv0*  aodV0 = 0x0;
  AliCFPair* pair  = 0x0;
  Int_t      labMCV0 = 0;

  for (Int_t iV0 = 0; iV0<fInputEvent->GetNumberOfV0s(); iV0++) {

    if (isESDEvent) {
      esdV0 = fESD->GetV0(iV0);
      //check if mother reconstructed V0 can be associated to a MC V0
      labMCV0 = IsMcV0(esdV0);
      if (labMCV0 <0) continue;
      pair = new AliCFPair(esdV0,fESD);
    }
    else if (isAODEvent) {
      aodV0 = fAOD->GetV0(iV0);
      labMCV0 = IsMcV0(aodV0);
      if (labMCV0 <0) continue;
      pair = new AliCFPair(aodV0);
    }
    else {
      Error("UserExec","Error: input data file is not an ESD nor an AOD");
      return;
    }
    
    pair->SetV0PDG(fV0PDG);
    pair->SetLabel(labMCV0);
    if (!fCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pair)) continue;
    
    //check if associated MC v0 passes the cuts
    AliMCParticle* mcV0 = (AliMCParticle*)fMCEvent->GetTrack(labMCV0);
    if (!mcV0) continue;

    if (!fCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,mcV0)) continue; 
    
    //fill the container
    Double_t pt, rapidity;
    
    if (isESDEvent) {
      Double_t mom[3];
      esdV0->GetPxPyPz(mom[0],mom[1],mom[2]);
      Double_t pt2 = mom[0]*mom[0]+mom[1]*mom[1] ;
      pt  = TMath::Sqrt(pt2);
      Double_t energy  = TMath::Sqrt(pt2 + mom[2]*mom[2] + TMath::Power(TDatabasePDG::Instance()->GetParticle(fV0PDG)->Mass(),2));
      rapidity = GetRapidity(energy,mom[2]) ;
    }
    else {
      pt = aodV0->Pt();
      rapidity = aodV0->Y(fV0PDG);
    }
    
    containerInput[0] = pt ;
    containerInput[1] = rapidity ;
    fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed) ;   
    if (!fCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,pair)) continue ;
    fCFManager->GetParticleContainer()->Fill(containerInput,kStepSelected);
    
    delete pair;
  }
  
  fHistEventsProcessed->Fill(0);
  /* PostData(0) is taken care of by AliAnalysisTaskSE */
  PostData(1,fHistEventsProcessed) ;
  PostData(2,fCFManager->GetParticleContainer()) ;
}


//___________________________________________________________________________
void AliCFV0Task::Terminate(Option_t*)
{
  // The Terminate() function is the last function to be called during
  // a query. It always runs on the client, it can be used to present
  // the results graphically or save the results to file.

  Info("Terminate","");
  AliAnalysisTaskSE::Terminate();


  //draw some example plots....

  AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));

  TH1D* h00 =   cont->ShowProjection(0,0) ;
  TH1D* h01 =   cont->ShowProjection(0,1) ;
  TH1D* h02 =   cont->ShowProjection(0,2) ;
  TH1D* h03 =   cont->ShowProjection(0,3) ;

  TH1D* h10 =   cont->ShowProjection(1,0) ;
  TH1D* h11 =   cont->ShowProjection(1,1) ;
  TH1D* h12 =   cont->ShowProjection(1,2) ;
  TH1D* h13 =   cont->ShowProjection(1,3) ;

  Double_t max1 = h00->GetMaximum();
  Double_t max2 = h10->GetMaximum();

  h00->GetYaxis()->SetRangeUser(0,max1*1.2);
  h01->GetYaxis()->SetRangeUser(0,max1*1.2);
  h02->GetYaxis()->SetRangeUser(0,max1*1.2);
  h03->GetYaxis()->SetRangeUser(0,max1*1.2);

  h10->GetYaxis()->SetRangeUser(0,max2*1.2);
  h11->GetYaxis()->SetRangeUser(0,max2*1.2);
  h12->GetYaxis()->SetRangeUser(0,max2*1.2);
  h13->GetYaxis()->SetRangeUser(0,max2*1.2);

  h00->SetMarkerStyle(23) ;
  h01->SetMarkerStyle(24) ;
  h02->SetMarkerStyle(25) ;
  h03->SetMarkerStyle(26) ;

  h10->SetMarkerStyle(23) ;
  h11->SetMarkerStyle(24) ;
  h12->SetMarkerStyle(25) ;
  h13->SetMarkerStyle(26) ;

  TCanvas * c =new TCanvas("c","",1400,800);
  c->Divide(4,2);

  c->cd(1);
  h00->Draw("p");
  c->cd(2);
  h01->Draw("p");
  c->cd(3);
  h02->Draw("p");
  c->cd(4);
  h03->Draw("p");
  c->cd(5);
  h10->Draw("p");
  c->cd(6);
  h11->Draw("p");
  c->cd(7);
  h12->Draw("p");
  c->cd(8);
  h13->Draw("p");

  c->SaveAs("plots.eps");
}

//___________________________________________________________________________
void AliCFV0Task::UserCreateOutputObjects() {
  //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
  //TO BE SET BEFORE THE EXECUTION OF THE TASK
  //
  Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());

  //slot #1
  OpenFile(1);
  fHistEventsProcessed = new TH1I("fHistEventsProcessed","",1,0,1) ;
}

//___________________________________________________________________________
Int_t AliCFV0Task::GetV0Label(UInt_t lab1, UInt_t lab2) const {
  //
  // returns the label of the V0, given the labels of the 2 daughter tracks
  // returns -1 if the V0 is fake
  //
  
  AliStack* stack = fMCEvent->Stack();
  TParticle* part1 = stack->Particle(lab1) ;
  TParticle* part2 = stack->Particle(lab2) ;
  
  Int_t part1MotherLab=part1->GetFirstMother();
  Int_t part2MotherLab=part2->GetFirstMother();

  if (part1MotherLab==-1 || part2MotherLab==-1) return -1 ;
  if (part1MotherLab != part2MotherLab )        return -1 ;
  
  if (stack->Particle(part1MotherLab)->GetPdgCode() != fV0PDG ) return -1 ;

  switch (fV0PDG) {
  case kK0Short : 
    if ( (part1->GetPdgCode()==  -211 && part2->GetPdgCode()==  211)  || 
	 (part1->GetPdgCode()==   211 && part2->GetPdgCode()== -211)   )   return part1MotherLab ;
    break ;
  case kLambda0 :
    if ( (part1->GetPdgCode()==  -211 && part2->GetPdgCode()== 2212)  || 
	 (part1->GetPdgCode()==  2212 && part2->GetPdgCode()== -211)   )   return part1MotherLab ;
    break ;
  case kLambda0Bar :
    if ( (part1->GetPdgCode()==   211 && part2->GetPdgCode()== -2212) || 
	 (part1->GetPdgCode()== -2212 && part2->GetPdgCode()==   211)  )   return part1MotherLab ;
    break ;
  default :
    return -1;
    break ;
  }
  
  return -1 ;
}

//___________________________________________________________________________
Double_t AliCFV0Task::GetRapidity(Double_t energy, Double_t pz) {
  //
  // calculates rapidity, checking energy is larger than pz, otherwise returns 999.
  //
  
  if (energy == pz || energy == -pz) {
    printf("GetRapidity : ERROR : rapidity for 4-vector with E = Pz -- infinite result");
    return 999;
  }
  if (energy < pz) {
    printf("GetRapidity : ERROR : rapidity for 4-vector with E = Pz -- infinite result");
    return 999;
  }
  Double_t y = 0.5 * log((energy + pz) / (energy - pz));
  return y;
} 

//___________________________________________________________________________
Int_t AliCFV0Task::IsMcV0(AliESDv0* v0) const {
  //
  // check if the passed V0 is associated to a MC one, 
  //   and returns the corresponding geant label.
  // returns -1 if the V0 is fake (i.e. label<0).
  //

  AliESDEvent* fESD = dynamic_cast<AliESDEvent*>(fInputEvent);

  Int_t nindex = v0->GetNindex();
  Int_t pindex = v0->GetPindex();
  AliESDtrack *nTrack = fESD->GetTrack(nindex) ;
  AliESDtrack *pTrack = fESD->GetTrack(pindex) ;
  
  if (!nTrack || !pTrack) return -1 ;

  Int_t nlab  = nTrack->GetLabel() ;
  Int_t plab  = pTrack->GetLabel() ;
  
  if (nlab <0 || plab <0) return -1 ;

  return GetV0Label((UInt_t)nlab,(UInt_t)plab) ;
}

//___________________________________________________________________________
Int_t AliCFV0Task::IsMcV0(AliAODv0* v0) const {
  //
  // check if the passed V0 is associated to a MC one, 
  //   and returns the corresponding geant label.
  // returns -1 if the V0 is fake (i.e. label<0).
  //

  AliAODVertex * vtx = v0->GetSecondaryVtx();

  AliAODTrack *nTrack = (AliAODTrack*)vtx->GetDaughter(1); //neg is filled after pos in AliAnalysisTaskStrange
  AliAODTrack *pTrack = (AliAODTrack*)vtx->GetDaughter(0);
  
  if (!nTrack || !pTrack) return -1 ;

  Int_t nlab  = nTrack->GetLabel() ;
  Int_t plab  = pTrack->GetLabel() ;

  if (nlab <0 || plab <0) return -1 ;

  return GetV0Label((UInt_t)nlab,(UInt_t)plab) ;
}


//___________________________________________________________________________
void AliCFV0Task::RebuildV0s(AliESDEvent* fESD) {
  
  fESD->ResetV0s();

  //These are pp cuts : to change if Pb+Pb !
  Double_t cuts[]={33,  // max. allowed chi2
		   0.01,// min. allowed negative daughter's impact parameter 
		   0.01,// min. allowed positive daughter's impact parameter 
		   0.5,// max. allowed DCA between the daughter tracks
		   0.98,// max. allowed cosine of V0's pointing angle
		   0.2,  // min. radius of the fiducial volume
		   100.   // max. radius of the fiducial volume
  };
  //   AliV0vertexer* v0Vertexer = new AliV0vertexer(cuts);
  //   v0Vertexer->SetVertex(primVertexPos);
  AliV0vertexer* v0Vertexer = new AliV0vertexer();
  v0Vertexer->SetCuts(cuts) ;
  v0Vertexer->Tracks2V0vertices(fESD);
  delete v0Vertexer ;
}

#endif
Commit	Line	Data
2fbc0b17	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	//-----------------------------------------------------------------------
	17	// Example of task running on AliEn (CAF?)
	18	// which provides standard way of calculating acceptance and efficiency
	19	// between different steps of the procedure.
	20	// The ouptut of the task is a AliCFContainer from which the efficiencies
	21	// can be calculated
	22	//-----------------------------------------------------------------------
	23	// Author : R. Vernet, Consorzio Cometa - Catania (it)
	24	//-----------------------------------------------------------------------
	25
	26
	27	#ifndef ALICFV0TASK_CXX
	28	#define ALICFV0TASK_CXX
2fbc0b17	29
	30	#include "AliCFV0Task.h"
	31	#include "TCanvas.h"
	32	#include "AliStack.h"
	33	#include "TParticle.h"
	34	#include "TH1I.h"
2fbc0b17	35	#include "AliMCEvent.h"
2fbc0b17	36	#include "AliCFManager.h"
2fbc0b17	37	#include "AliCFContainer.h"
2fbc0b17	38	#include "AliESDtrack.h"
	39	#include "AliLog.h"
	40	#include "AliESDv0.h"
	41	#include "AliV0vertexer.h"
	42	#include "AliCFPair.h"
318a0e1f	43	#include "AliESDEvent.h"
86c32a36	44	#include "AliAODEvent.h"
318a0e1f	45	#include "TChain.h"
318a0e1f	46	#include "AliCFParticleGenCuts.h"
86c32a36	47	#include "AliAODv0.h"
0c01ae65	48	#include "TDatabasePDG.h"
2fbc0b17	49
	50	//__________________________________________________________________________
	51	AliCFV0Task::AliCFV0Task() :
318a0e1f	52	AliAnalysisTaskSE(),
2fbc0b17	53	fRebuildV0s(0),
2fbc0b17	54	fV0PDG(310),
2fbc0b17	55	fCFManager(0x0),
	56	fHistEventsProcessed(0x0)
	57	{
318a0e1f	58	//
	59	//Default ctor
	60	//
2fbc0b17	61	}
	62	//___________________________________________________________________________
	63	AliCFV0Task::AliCFV0Task(const Char_t* name) :
318a0e1f	64	AliAnalysisTaskSE(name),
2fbc0b17	65	fRebuildV0s(0),
2fbc0b17	66	fV0PDG(310),
2fbc0b17	67	fCFManager(0x0),
	68	fHistEventsProcessed(0x0)
	69	{
	70	//
	71	// Constructor. Initialization of Inputs and Outputs
	72	//
	73	Info("AliCFV0Task","Calling Constructor");
318a0e1f	74	/*
	75	DefineInput(0) and DefineOutput(0)
	76	are taken care of by AliAnalysisTaskSE constructor
	77	*/
	78	DefineOutput(1,TH1I::Class());
	79	DefineOutput(2,AliCFContainer::Class());
	80	// DefineOutput(3,TList::Class());
2fbc0b17	81	}
	82
	83	//___________________________________________________________________________
	84	AliCFV0Task& AliCFV0Task::operator=(const AliCFV0Task& c)
	85	{
	86	//
	87	// Assignment operator
	88	//
	89	if (this!=&c) {
318a0e1f	90	AliAnalysisTaskSE::operator=(c) ;
2fbc0b17	91	fRebuildV0s = c.fRebuildV0s;
2fbc0b17	92	fV0PDG = c.fV0PDG;
2fbc0b17	93	fCFManager = c.fCFManager;
	94	fHistEventsProcessed = c.fHistEventsProcessed;
	95	}
	96	return *this;
	97	}
	98
	99	//___________________________________________________________________________
	100	AliCFV0Task::AliCFV0Task(const AliCFV0Task& c) :
318a0e1f	101	AliAnalysisTaskSE(c),
2fbc0b17	102	fRebuildV0s(c.fRebuildV0s),
2fbc0b17	103	fV0PDG(c.fV0PDG),
2fbc0b17	104	fCFManager(c.fCFManager),
	105	fHistEventsProcessed(c.fHistEventsProcessed)
	106	{
	107	//
	108	// Copy Constructor
	109	//
	110	}
	111
	112	//___________________________________________________________________________
	113	AliCFV0Task::~AliCFV0Task() {
	114	//
	115	//destructor
	116	//
	117	Info("~AliCFV0Task","Calling Destructor");
2fbc0b17	118	if (fCFManager) delete fCFManager ;
	119	if (fHistEventsProcessed) delete fHistEventsProcessed ;
	120	}
	121
2fbc0b17	122	//_________________________________________________
318a0e1f	123	void AliCFV0Task::UserExec(Option_t *)
2fbc0b17	124	{
	125	//
	126	// Main loop function
	127	//
318a0e1f	128	Info("UserExec","") ;
2fbc0b17	129
86c32a36	130	if (!fInputEvent) {
86c32a36	131	Error("UserExec","NO EVENT FOUND!");
318a0e1f	132	return;
318a0e1f	133	}
2fbc0b17	134
0c01ae65	135	if (!fMCEvent) {
	136	Error("UserExec","NO MC INFO FOUND!");
	137	return;
	138	}
	139
	140	fCFManager->SetMCEventInfo (fMCEvent);
	141	fCFManager->SetRecEventInfo(fInputEvent);
2fbc0b17	142
86c32a36	143	Bool_t isESDEvent = strcmp(fInputEvent->ClassName(),"AliESDEvent") == 0 ? kTRUE : kFALSE ;
	144	Bool_t isAODEvent = strcmp(fInputEvent->ClassName(),"AliAODEvent") == 0 ? kTRUE : kFALSE ;
	145
	146	AliESDEvent* fESD = dynamic_cast<AliESDEvent*>(fInputEvent);
	147	AliAODEvent* fAOD = dynamic_cast<AliAODEvent*>(fInputEvent);
	148
	149
2fbc0b17	150	// MC-event selection
	151	Double_t containerInput[2] ;
	152
	153	//loop on the MC event
318a0e1f	154	Info("UserExec","Looping on MC event");
318a0e1f	155	for (Int_t ipart=0; ipart<fMCEvent->GetNumberOfTracks(); ipart++) {
0c01ae65	156	AliMCParticle mcPart = (AliMCParticle)fMCEvent->GetTrack(ipart);
2fbc0b17	157
	158	//check the MC-level cuts
	159	if (!fCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,mcPart)) continue;
	160	containerInput[0] = mcPart->Pt();
1d519a00	161	containerInput[1] = mcPart->Y() ;
2fbc0b17	162	//fill the container for Gen-level selection
	163	fCFManager->GetParticleContainer()->Fill(containerInput,kStepGenerated);
	164
	165	//check the Acceptance-level cuts
	166	if (!fCFManager->CheckParticleCuts(AliCFManager::kPartAccCuts,mcPart)) continue;
	167	//fill the container for Acceptance-level selection
	168	fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructible);
	169	}
	170
	171
	172	//Now go to rec level
86c32a36	173	Info("UserExec","Looping on %s",fInputEvent->ClassName());
86c32a36	174
86c32a36	175	if (isESDEvent && fRebuildV0s) RebuildV0s(fESD) ;
2fbc0b17	176
86c32a36	177	Info("UserExec","There are %d V0s in event",fInputEvent->GetNumberOfV0s());
	178
	179	AliESDv0* esdV0 = 0x0;
	180	AliAODv0* aodV0 = 0x0;
	181	AliCFPair* pair = 0x0;
	182	Int_t labMCV0 = 0;
	183
	184	for (Int_t iV0 = 0; iV0<fInputEvent->GetNumberOfV0s(); iV0++) {
	185
	186	if (isESDEvent) {
	187	esdV0 = fESD->GetV0(iV0);
	188	//check if mother reconstructed V0 can be associated to a MC V0
	189	labMCV0 = IsMcV0(esdV0);
	190	if (labMCV0 <0) continue;
	191	pair = new AliCFPair(esdV0,fESD);
	192	}
	193	else if (isAODEvent) {
	194	aodV0 = fAOD->GetV0(iV0);
	195	labMCV0 = IsMcV0(aodV0);
	196	if (labMCV0 <0) continue;
	197	pair = new AliCFPair(aodV0);
	198	}
	199	else {
	200	Error("UserExec","Error: input data file is not an ESD nor an AOD");
	201	return;
	202	}
	203
	204	pair->SetV0PDG(fV0PDG);
	205	pair->SetLabel(labMCV0);
	206	if (!fCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pair)) continue;
	207
2fbc0b17	208	//check if associated MC v0 passes the cuts
0c01ae65	209	AliMCParticle* mcV0 = (AliMCParticle*)fMCEvent->GetTrack(labMCV0);
2fbc0b17	210	if (!mcV0) continue;
	211
	212	if (!fCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,mcV0)) continue;
	213
	214	//fill the container
86c32a36	215	Double_t pt, rapidity;
	216
	217	if (isESDEvent) {
	218	Double_t mom[3];
	219	esdV0->GetPxPyPz(mom[0],mom[1],mom[2]);
	220	Double_t pt2 = mom[0]mom[0]+mom[1]mom[1] ;
	221	pt = TMath::Sqrt(pt2);
	222	Double_t energy = TMath::Sqrt(pt2 + mom[2]*mom[2] + TMath::Power(TDatabasePDG::Instance()->GetParticle(fV0PDG)->Mass(),2));
	223	rapidity = GetRapidity(energy,mom[2]) ;
	224	}
	225	else {
	226	pt = aodV0->Pt();
	227	rapidity = aodV0->Y(fV0PDG);
	228	}
	229
2fbc0b17	230	containerInput[0] = pt ;
86c32a36	231	containerInput[1] = rapidity ;
2fbc0b17	232	fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed) ;
86c32a36	233	if (!fCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,pair)) continue ;
2fbc0b17	234	fCFManager->GetParticleContainer()->Fill(containerInput,kStepSelected);
86c32a36	235
86c32a36	236	delete pair;
2fbc0b17	237	}
	238
	239	fHistEventsProcessed->Fill(0);
318a0e1f	240	/* PostData(0) is taken care of by AliAnalysisTaskSE */
	241	PostData(1,fHistEventsProcessed) ;
	242	PostData(2,fCFManager->GetParticleContainer()) ;
2fbc0b17	243	}
	244
	245
	246	//___________________________________________________________________________
	247	void AliCFV0Task::Terminate(Option_t*)
	248	{
	249	// The Terminate() function is the last function to be called during
	250	// a query. It always runs on the client, it can be used to present
	251	// the results graphically or save the results to file.
	252
	253	Info("Terminate","");
318a0e1f	254	AliAnalysisTaskSE::Terminate();
	255
	256
	257	//draw some example plots....
	258
	259	AliCFContainer cont= dynamic_cast<AliCFContainer> (GetOutputData(2));
	260
	261	TH1D* h00 = cont->ShowProjection(0,0) ;
	262	TH1D* h01 = cont->ShowProjection(0,1) ;
	263	TH1D* h02 = cont->ShowProjection(0,2) ;
	264	TH1D* h03 = cont->ShowProjection(0,3) ;
2fbc0b17	265
318a0e1f	266	TH1D* h10 = cont->ShowProjection(1,0) ;
	267	TH1D* h11 = cont->ShowProjection(1,1) ;
	268	TH1D* h12 = cont->ShowProjection(1,2) ;
	269	TH1D* h13 = cont->ShowProjection(1,3) ;
2fbc0b17	270
318a0e1f	271	Double_t max1 = h00->GetMaximum();
318a0e1f	272	Double_t max2 = h10->GetMaximum();
2fbc0b17	273
318a0e1f	274	h00->GetYaxis()->SetRangeUser(0,max1*1.2);
	275	h01->GetYaxis()->SetRangeUser(0,max1*1.2);
	276	h02->GetYaxis()->SetRangeUser(0,max1*1.2);
	277	h03->GetYaxis()->SetRangeUser(0,max1*1.2);
2fbc0b17	278
318a0e1f	279	h10->GetYaxis()->SetRangeUser(0,max2*1.2);
	280	h11->GetYaxis()->SetRangeUser(0,max2*1.2);
	281	h12->GetYaxis()->SetRangeUser(0,max2*1.2);
	282	h13->GetYaxis()->SetRangeUser(0,max2*1.2);
2fbc0b17	283
318a0e1f	284	h00->SetMarkerStyle(23) ;
	285	h01->SetMarkerStyle(24) ;
	286	h02->SetMarkerStyle(25) ;
	287	h03->SetMarkerStyle(26) ;
2fbc0b17	288
318a0e1f	289	h10->SetMarkerStyle(23) ;
	290	h11->SetMarkerStyle(24) ;
	291	h12->SetMarkerStyle(25) ;
	292	h13->SetMarkerStyle(26) ;
2fbc0b17	293
	294	TCanvas * c =new TCanvas("c","",1400,800);
	295	c->Divide(4,2);
	296
	297	c->cd(1);
318a0e1f	298	h00->Draw("p");
2fbc0b17	299	c->cd(2);
318a0e1f	300	h01->Draw("p");
2fbc0b17	301	c->cd(3);
318a0e1f	302	h02->Draw("p");
2fbc0b17	303	c->cd(4);
318a0e1f	304	h03->Draw("p");
2fbc0b17	305	c->cd(5);
318a0e1f	306	h10->Draw("p");
2fbc0b17	307	c->cd(6);
318a0e1f	308	h11->Draw("p");
2fbc0b17	309	c->cd(7);
318a0e1f	310	h12->Draw("p");
2fbc0b17	311	c->cd(8);
318a0e1f	312	h13->Draw("p");
2fbc0b17	313
2fbc0b17	314	c->SaveAs("plots.eps");
2fbc0b17	315	}
	316
	317	//___________________________________________________________________________
318a0e1f	318	void AliCFV0Task::UserCreateOutputObjects() {
2fbc0b17	319	//HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
	320	//TO BE SET BEFORE THE EXECUTION OF THE TASK
	321	//
318a0e1f	322	Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
2fbc0b17	323
318a0e1f	324	//slot #1
318a0e1f	325	OpenFile(1);
2fbc0b17	326	fHistEventsProcessed = new TH1I("fHistEventsProcessed","",1,0,1) ;
	327	}
	328
	329	//___________________________________________________________________________
318a0e1f	330	Int_t AliCFV0Task::GetV0Label(UInt_t lab1, UInt_t lab2) const {
2fbc0b17	331	//
	332	// returns the label of the V0, given the labels of the 2 daughter tracks
	333	// returns -1 if the V0 is fake
	334	//
	335
318a0e1f	336	AliStack* stack = fMCEvent->Stack();
2fbc0b17	337	TParticle* part1 = stack->Particle(lab1) ;
	338	TParticle* part2 = stack->Particle(lab2) ;
	339
	340	Int_t part1MotherLab=part1->GetFirstMother();
	341	Int_t part2MotherLab=part2->GetFirstMother();
318a0e1f	342
2fbc0b17	343	if (part1MotherLab==-1 \|\| part2MotherLab==-1) return -1 ;
2fbc0b17	344	if (part1MotherLab != part2MotherLab ) return -1 ;
318a0e1f	345
2fbc0b17	346	if (stack->Particle(part1MotherLab)->GetPdgCode() != fV0PDG ) return -1 ;
	347
	348	switch (fV0PDG) {
	349	case kK0Short :
	350	if ( (part1->GetPdgCode()== -211 && part2->GetPdgCode()== 211) \|\|
	351	(part1->GetPdgCode()== 211 && part2->GetPdgCode()== -211) ) return part1MotherLab ;
	352	break ;
	353	case kLambda0 :
	354	if ( (part1->GetPdgCode()== -211 && part2->GetPdgCode()== 2212) \|\|
	355	(part1->GetPdgCode()== 2212 && part2->GetPdgCode()== -211) ) return part1MotherLab ;
	356	break ;
	357	case kLambda0Bar :
	358	if ( (part1->GetPdgCode()== 211 && part2->GetPdgCode()== -2212) \|\|
	359	(part1->GetPdgCode()== -2212 && part2->GetPdgCode()== 211) ) return part1MotherLab ;
	360	break ;
	361	default :
	362	return -1;
	363	break ;
	364	}
	365
	366	return -1 ;
	367	}
	368
	369	//___________________________________________________________________________
	370	Double_t AliCFV0Task::GetRapidity(Double_t energy, Double_t pz) {
86c32a36	371	//
	372	// calculates rapidity, checking energy is larger than pz, otherwise returns 999.
	373	//
	374
2fbc0b17	375	if (energy == pz \|\| energy == -pz) {
	376	printf("GetRapidity : ERROR : rapidity for 4-vector with E = Pz -- infinite result");
	377	return 999;
	378	}
	379	if (energy < pz) {
	380	printf("GetRapidity : ERROR : rapidity for 4-vector with E = Pz -- infinite result");
	381	return 999;
	382	}
	383	Double_t y = 0.5 * log((energy + pz) / (energy - pz));
	384	return y;
	385	}
	386
	387	//___________________________________________________________________________
318a0e1f	388	Int_t AliCFV0Task::IsMcV0(AliESDv0* v0) const {
86c32a36	389	//
	390	// check if the passed V0 is associated to a MC one,
	391	// and returns the corresponding geant label.
	392	// returns -1 if the V0 is fake (i.e. label<0).
	393	//
318a0e1f	394
	395	AliESDEvent* fESD = dynamic_cast<AliESDEvent*>(fInputEvent);
	396
2fbc0b17	397	Int_t nindex = v0->GetNindex();
2fbc0b17	398	Int_t pindex = v0->GetPindex();
318a0e1f	399	AliESDtrack *nTrack = fESD->GetTrack(nindex) ;
318a0e1f	400	AliESDtrack *pTrack = fESD->GetTrack(pindex) ;
2fbc0b17	401
	402	if (!nTrack \|\| !pTrack) return -1 ;
	403
	404	Int_t nlab = nTrack->GetLabel() ;
	405	Int_t plab = pTrack->GetLabel() ;
	406
	407	if (nlab <0 \|\| plab <0) return -1 ;
	408
318a0e1f	409	return GetV0Label((UInt_t)nlab,(UInt_t)plab) ;
2fbc0b17	410	}
2fbc0b17	411
86c32a36	412	//___________________________________________________________________________
	413	Int_t AliCFV0Task::IsMcV0(AliAODv0* v0) const {
	414	//
	415	// check if the passed V0 is associated to a MC one,
	416	// and returns the corresponding geant label.
	417	// returns -1 if the V0 is fake (i.e. label<0).
	418	//
	419
	420	AliAODVertex * vtx = v0->GetSecondaryVtx();
	421
	422	AliAODTrack nTrack = (AliAODTrack)vtx->GetDaughter(1); //neg is filled after pos in AliAnalysisTaskStrange
	423	AliAODTrack pTrack = (AliAODTrack)vtx->GetDaughter(0);
	424
	425	if (!nTrack \|\| !pTrack) return -1 ;
	426
	427	Int_t nlab = nTrack->GetLabel() ;
	428	Int_t plab = pTrack->GetLabel() ;
	429
	430	if (nlab <0 \|\| plab <0) return -1 ;
	431
	432	return GetV0Label((UInt_t)nlab,(UInt_t)plab) ;
	433	}
	434
2fbc0b17	435
2fbc0b17	436	//___________________________________________________________________________
318a0e1f	437	void AliCFV0Task::RebuildV0s(AliESDEvent* fESD) {
2fbc0b17	438
	439	fESD->ResetV0s();
	440
	441	//These are pp cuts : to change if Pb+Pb !
	442	Double_t cuts[]={33, // max. allowed chi2
318a0e1f	443	0.01,// min. allowed negative daughter's impact parameter
	444	0.01,// min. allowed positive daughter's impact parameter
	445	0.5,// max. allowed DCA between the daughter tracks
	446	0.98,// max. allowed cosine of V0's pointing angle
	447	0.2, // min. radius of the fiducial volume
2fbc0b17	448	100. // max. radius of the fiducial volume
	449	};
	450	// AliV0vertexer* v0Vertexer = new AliV0vertexer(cuts);
	451	// v0Vertexer->SetVertex(primVertexPos);
	452	AliV0vertexer* v0Vertexer = new AliV0vertexer();
	453	v0Vertexer->SetCuts(cuts) ;
	454	v0Vertexer->Tracks2V0vertices(fESD);
	455	delete v0Vertexer ;
	456	}
	457
	458	#endif