[u/mrichter/AliRoot.git] / ANALYSIS / AliAnalysisTaskMCParticleFilter.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.                  *
 **************************************************************************/

//
//  Analysis task for Kinematic filtering
//  Fill AODtrackMC tracks from Kinematic stack
//
 
#include <TChain.h>
#include <TFile.h>
#include "TParticle.h"
#include "TParameter.h"
#include "TString.h"
#include "TList.h"
#include "TProfile.h"
#include "TH1F.h"


#include "AliAnalysisTaskMCParticleFilter.h"
#include "AliAnalysisManager.h"
#include "AliAnalysisFilter.h"
#include "AliHeader.h"
#include "AliStack.h"
#include "AliMCEvent.h"
#include "AliMCEventHandler.h"
#include "AliAODEvent.h"
#include "AliAODHeader.h"
#include "AliAODMCHeader.h"
#include "AliAODHandler.h"
#include "AliAODVertex.h"
#include "AliAODMCParticle.h"


#include "AliLog.h"

ClassImp(AliAnalysisTaskMCParticleFilter)

////////////////////////////////////////////////////////////////////////

//____________________________________________________________________
AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter():
AliAnalysisTaskSE(),
  fTrackFilterMother(0x0),
  fHistList(0x0)
{
  // Default constructor
}

Bool_t AliAnalysisTaskMCParticleFilter::Notify()
{
  //
  // Implemented Notify() to read the cross sections
  // and number of trials from pyxsec.root
  // 
  TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
  Double_t xsection = 0;
  UInt_t   ntrials  = 0;
  if(tree){
    TFile *curfile = tree->GetCurrentFile();
    if (!curfile) {
      Error("Notify","No current file");
      return kFALSE;
    }

    TString fileName(curfile->GetName());
    if(fileName.Contains("AliESDs.root")){
        fileName.ReplaceAll("AliESDs.root", "pyxsec.root");
    }
    else if(fileName.Contains("AliAOD.root")){
        fileName.ReplaceAll("AliAOD.root", "pyxsec.root");
    }
    else if(fileName.Contains("galice.root")){
        // for running with galice and kinematics alone...                      
        fileName.ReplaceAll("galice.root", "pyxsec.root");
    }


    TFile *fxsec = TFile::Open(fileName.Data());
    if(!fxsec){
      AliInfo(Form("%s:%d %s not found in the Input",(char*)__FILE__,__LINE__,fileName.Data()));
      // not a severe condition we just do not have the information...
      return kTRUE;
    }
    TTree *xtree = (TTree*)fxsec->Get("Xsection");
    if(!xtree){
      Printf("%s:%d tree not found in the pyxsec.root",(char*)__FILE__,__LINE__);
      return kTRUE;
    }
    xtree->SetBranchAddress("xsection",&xsection);
    xtree->SetBranchAddress("ntrials",&ntrials);
    xtree->GetEntry(0);
    ((TProfile*)(fHistList->FindObject("h1Xsec")))->Fill("<#sigma>",xsection);
    ((TH1F*)(fHistList->FindObject("h1Trials")))->Fill("#sum{ntrials}",ntrials); 

  }
  return kTRUE;
}


//____________________________________________________________________
AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter(const char* name):
    AliAnalysisTaskSE(name),
    fTrackFilterMother(0x0),
    fHistList(0x0)
{
  // Default constructor
  DefineOutput(1,TList::Class());  
}

/*
//____________________________________________________________________
AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter(const AliAnalysisTaskMCParticleFilter& obj):
    AliAnalysisTaskSE(obj),
    fTrackFilterMother(obj.fTrackFilterMother)
{
  // Copy constructor
}
*/
//____________________________________________________________________
AliAnalysisTaskMCParticleFilter::~AliAnalysisTaskMCParticleFilter()
{
  //  if( fTrackFilterMother ) delete fTrackFilterMother;
}

/*
//____________________________________________________________________
AliAnalysisTaskMCParticleFilter& AliAnalysisTaskMCParticleFilter::operator=(const AliAnalysisTaskMCParticleFilter& other)
{
// Assignment
  if(this!=&other) {
    AliAnalysisTaskSE::operator=(other);
    fTrackFilterMother = other.fTrackFilterMother;
  }
  return *this;
}
*/
//____________________________________________________________________
void AliAnalysisTaskMCParticleFilter::UserCreateOutputObjects()
{
  // Create the output container
    if (OutputTree()&&fTrackFilterMother) 
	OutputTree()->GetUserInfo()->Add(fTrackFilterMother);

    // this part is mainly needed to set the MCEventHandler
    // to the AODHandler, this will not be needed when
    // AODHandler goes to ANALYSISalice
    // setting in the steering macro will not work on proof :(
    // so we have to do it in a task

    // the branch booking can also go into the AODHandler then


    // mcparticles
    TClonesArray *tca = new TClonesArray("AliAODMCParticle", 0);
    tca->SetName(AliAODMCParticle::StdBranchName());
    AddAODBranch("TClonesArray",&tca);

    // MC header...
    AliAODMCHeader *mcHeader = new AliAODMCHeader();
    mcHeader->SetName(AliAODMCHeader::StdBranchName());
    AddAODBranch("AliAODMCHeader",&mcHeader);    

    
    AliMCEventHandler *mcH = (AliMCEventHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler()); 
    AliAODHandler *aodH = dynamic_cast<AliAODHandler*> ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
    if(!aodH){
      Printf("%s:&d Could not get AODHandler",(char*)__FILE__,__LINE__);
      return;
    }
    aodH->SetMCEventHandler(mcH);


    // these are histograms, for averaging and summing
    // do it via histograms to be PROOF-proof
    // which merges the results from different workers
    // histos are not saved reside only in memory
    
    
    fHistList = new TList();
    TProfile *h1Xsec = new TProfile("h1Xsec","xsec from pyxsec.root",1,0,1);
    h1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
    fHistList->Add(h1Xsec);
    TH1F* h1Trials = new TH1F("h1Trials","trials from pyxsec.root",1,0,1);
    h1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
    fHistList->Add(h1Trials);

    
}

//____________________________________________________________________
void AliAnalysisTaskMCParticleFilter::UserExec(Option_t */*option*/)
{
// Execute analysis for current event
//

// Fill AOD tracks from Kinematic stack
    
  // get AliAOD Event 
  AliAODEvent* aod = AODEvent();
  if (!aod) {
      AliWarning("No Output Handler connected, doing nothing !") ;
      return;
  }

  AliMCEventHandler *mcH = (AliMCEventHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler()); 
  if(!mcH){ 
    AliWarning("No MC handler Found");
    return;
  }


  // fetch the output 
  // Fill control histos

  // tmp array for holding the mctracks
  // need to set mother and duagther __before__
  // filling in the tree

  AliMCEvent *mcE = MCEvent();

  if(!mcE){
    AliWarning("No MC event Found");
    return;
  }

  Int_t np    = mcE->GetNumberOfTracks();
  Int_t nprim = mcE->GetNumberOfPrimaries();
  // TODO ADD MC VERTEX
  
  // We take all real primaries


  Int_t j=0;
  for (Int_t ip = 0; ip < np; ip++){
    AliMCParticle* mcpart =  mcE->GetTrack(ip);
    TParticle* part = mcpart->Particle();
    Float_t xv = part->Vx();
    Float_t yv = part->Vy();
    Float_t zv = part->Vz();
    Float_t rv = TMath::Sqrt(xv * xv + yv * yv);
      
    Bool_t write = kFALSE;
    
    if (ip < nprim) {
      // Select the primary event
      write = kTRUE;
    } else if (part->GetUniqueID() == 4) {
      // Particles from decay
      // Check that the decay chain ends at a primary particle
      AliMCParticle* mother = mcpart;
      Int_t imo = mcpart->GetMother();
      while((imo >= nprim) && (mother->GetUniqueID() == 4)) {
	mother = mcE->GetTrack(imo);
	imo =  mother->GetMother();
      }
      // Select according to pseudorapidity and production point of primary ancestor
      if (imo < nprim && Select(mcE->GetTrack(imo)->Particle(), rv, zv))write = kTRUE;         
    } else if (part->GetUniqueID() == 5) {
      // Now look for pair production
      Int_t imo = mcpart->GetMother();
      if (imo < nprim) {
	// Select, if the gamma is a primary
	write = kTRUE;
      } else {
	// Check if the gamma comes from the decay chain of a primary particle
	AliMCParticle* mother = mcE->GetTrack(imo);
	imo = mother->GetMother();
	while((imo >= nprim) && (mother->GetUniqueID() == 4)) {
	  mother =  mcE->GetTrack(imo);
	  imo =  mother->GetMother();
	}
	// Select according to pseudorapidity and production point 
	if (imo < nprim && Select(mother->Particle(), rv, zv)) 
	  write = kTRUE;
      }
    }
    /*
    else if (part->GetUniqueID() == 13){
      // Evaporation
      // Check that we end at a primary particle
      TParticle* mother = part;
      Int_t imo = part->GetFirstMother();
      while((imo >= nprim) && ((mother->GetUniqueID() == 4 ) || ( mother->GetUniqueID() == 13))) {
	mother =  mcE->Stack()->Particle(imo);
	imo =  mother->GetFirstMother();
      }
      // Select according to pseudorapidity and production point 
	if (imo < nprim && Select(mother, rv, zv)) 
	  write = kTRUE;
    }
    */    
    if (write) {
      if(mcH)mcH->SelectParticle(ip);
      j++;
    }
  }


    // check for charm daughters
  static int  iTaken = 0;
  static int  iAll = 0;
  static int  iCharm = 0;
  for (Int_t ip = 0; ip < np; ip++){
    AliMCParticle* mcpart =  mcE->GetTrack(ip);
    TParticle* part = mcpart->Particle();

    //    if((TMath::Abs(part->GetPdgCode())/400)==1){
    if((TMath::Abs(part->GetPdgCode()))==411){
      // cases 
      iCharm++;
      Int_t d0 =  mcpart->GetFirstDaughter();
      Int_t d1 =  mcpart->GetLastDaughter();
      if(d0>0&&d1>0){
	for(int id = d0;id <= d1;id++){
	  // AliMCParticle* daughter =  mcE->GetTrack(id);
	  iAll++;
	  if(mcH->IsParticleSelected(id))iTaken++;
	}
      }
    }
  }

  AliInfo(Form("Taken daughters %d/%d of %d charm",iTaken,iAll,iCharm));

  PostData(1,fHistList);

  return;
}

void AliAnalysisTaskMCParticleFilter::Terminate(Option_t */*option*/){
  //
  // Terminate the execution save the PYTHIA x_section to the UserInfo()
  //


  TTree *tree = (TTree*)GetOutputData(0);
  fHistList = (TList*)GetOutputData(1);
  if(!fHistList){
    Printf("%s:%d Output list not found",(char*)__FILE__,__LINE__);
    return;
  }

  if(!tree){
    Printf("%s:%d Output tree not found",(char*)__FILE__,__LINE__);
    return;
  }

  TProfile *hXsec = ((TProfile*)(fHistList->FindObject("h1Xsec")));
  TH1F* hTrials  = ((TH1F*)(fHistList->FindObject("h1Trials")));
  if(!hXsec)return;
  if(!hTrials)return;

  Float_t xsec = hXsec->GetBinContent(1);
  Float_t trials = hTrials->Integral();
  AliInfo(Form("Average -section %.4E and total number of trials %E",xsec,trials));
  
  // This only works if the tree is not a special output
  // in this case it is already written
  /*
  TParameter<float> *x = new TParameter<float>("xsection",xsec);
  tree->GetUserInfo()->Add(x);
  TParameter<float> *t = new TParameter<float>("trials",trials);
  tree->GetUserInfo()->Add(t);
  */
  

}

Bool_t AliAnalysisTaskMCParticleFilter::Select(TParticle* part, Float_t rv, Float_t zv)
{
  // Selection accoring to eta of the mother and production point
  // This has to be refined !!!!!!

  // Esp if we don't have collisison in the central barrel but e.g. beam gas
  //  return kTRUE;

  Float_t eta = part->Eta();

  // central barrel consider everything in the ITS...
  // large eta window for smeared vertex and SPD acceptance (2 at z = 0)
  // larger for V0s in the TPC
  //  if (TMath::Abs(eta) < 2.5 && rv < 250. && TMath::Abs(zv)<255)return kTRUE;

  if (TMath::Abs(eta) < 2.5 && rv < 170)return kTRUE;   

  // Andreas' Cuts
  //  if (TMath::Abs(eta) < 1. && rv < 170)return kTRUE;   


  // Muon arm
  if(eta > -4.2 && eta < -2.3 && zv > -500.)return kTRUE; // Muon arms

  // PMD acceptance 2.3 <= eta < = 3.5
  //  if(eta>2.0&&eta<3.7)return kTRUE; 

  return kFALSE;
 
}
Commit	Line	Data
da97a08a	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	// Analysis task for Kinematic filtering
	18	// Fill AODtrackMC tracks from Kinematic stack
	19	//
	20
	21	#include <TChain.h>
	22	#include <TFile.h>
	23	#include "TParticle.h"
a3b51fd2	24	#include "TParameter.h"
da97a08a	25	#include "TString.h"
a3b51fd2	26	#include "TList.h"
	27	#include "TProfile.h"
	28	#include "TH1F.h"
	29
da97a08a	30
	31	#include "AliAnalysisTaskMCParticleFilter.h"
	32	#include "AliAnalysisManager.h"
	33	#include "AliAnalysisFilter.h"
	34	#include "AliHeader.h"
	35	#include "AliStack.h"
	36	#include "AliMCEvent.h"
	37	#include "AliMCEventHandler.h"
	38	#include "AliAODEvent.h"
	39	#include "AliAODHeader.h"
dce1b636	40	#include "AliAODMCHeader.h"
da97a08a	41	#include "AliAODHandler.h"
	42	#include "AliAODVertex.h"
	43	#include "AliAODMCParticle.h"
	44
	45
	46	#include "AliLog.h"
	47
	48	ClassImp(AliAnalysisTaskMCParticleFilter)
	49
	50	////////////////////////////////////////////////////////////////////////
	51
	52	//____________________________________________________________________
	53	AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter():
	54	AliAnalysisTaskSE(),
a3b51fd2	55	fTrackFilterMother(0x0),
a3b51fd2	56	fHistList(0x0)
da97a08a	57	{
	58	// Default constructor
	59	}
	60
a3b51fd2	61	Bool_t AliAnalysisTaskMCParticleFilter::Notify()
	62	{
	63	//
	64	// Implemented Notify() to read the cross sections
	65	// and number of trials from pyxsec.root
	66	//
	67	TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
	68	Double_t xsection = 0;
	69	UInt_t ntrials = 0;
	70	if(tree){
	71	TFile *curfile = tree->GetCurrentFile();
	72	if (!curfile) {
	73	Error("Notify","No current file");
	74	return kFALSE;
	75	}
	76
	77	TString fileName(curfile->GetName());
	78	if(fileName.Contains("AliESDs.root")){
	79	fileName.ReplaceAll("AliESDs.root", "pyxsec.root");
	80	}
	81	else if(fileName.Contains("AliAOD.root")){
	82	fileName.ReplaceAll("AliAOD.root", "pyxsec.root");
	83	}
	84	else if(fileName.Contains("galice.root")){
	85	// for running with galice and kinematics alone...
	86	fileName.ReplaceAll("galice.root", "pyxsec.root");
	87	}
	88
	89
	90	TFile *fxsec = TFile::Open(fileName.Data());
	91	if(!fxsec){
	92	AliInfo(Form("%s:%d %s not found in the Input",(char*)__FILE__,__LINE__,fileName.Data()));
	93	// not a severe condition we just do not have the information...
	94	return kTRUE;
	95	}
	96	TTree xtree = (TTree)fxsec->Get("Xsection");
	97	if(!xtree){
	98	Printf("%s:%d tree not found in the pyxsec.root",(char*)__FILE__,__LINE__);
	99	return kTRUE;
	100	}
	101	xtree->SetBranchAddress("xsection",&xsection);
	102	xtree->SetBranchAddress("ntrials",&ntrials);
	103	xtree->GetEntry(0);
	104	((TProfile*)(fHistList->FindObject("h1Xsec")))->Fill("<#sigma>",xsection);
	105	((TH1F*)(fHistList->FindObject("h1Trials")))->Fill("#sum{ntrials}",ntrials);
	106
	107	}
	108	return kTRUE;
	109	}
	110
	111
	112
da97a08a	113	//____________________________________________________________________
	114	AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter(const char* name):
	115	AliAnalysisTaskSE(name),
a3b51fd2	116	fTrackFilterMother(0x0),
a3b51fd2	117	fHistList(0x0)
da97a08a	118	{
da97a08a	119	// Default constructor
a3b51fd2	120	DefineOutput(1,TList::Class());
da97a08a	121	}
da97a08a	122
a3b51fd2	123	/*
da97a08a	124	//____________________________________________________________________
	125	AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter(const AliAnalysisTaskMCParticleFilter& obj):
	126	AliAnalysisTaskSE(obj),
	127	fTrackFilterMother(obj.fTrackFilterMother)
	128	{
	129	// Copy constructor
	130	}
a3b51fd2	131	*/
da97a08a	132	//____________________________________________________________________
	133	AliAnalysisTaskMCParticleFilter::~AliAnalysisTaskMCParticleFilter()
	134	{
	135	// if( fTrackFilterMother ) delete fTrackFilterMother;
	136	}
	137
a3b51fd2	138	/*
da97a08a	139	//____________________________________________________________________
	140	AliAnalysisTaskMCParticleFilter& AliAnalysisTaskMCParticleFilter::operator=(const AliAnalysisTaskMCParticleFilter& other)
	141	{
	142	// Assignment
	143	if(this!=&other) {
	144	AliAnalysisTaskSE::operator=(other);
	145	fTrackFilterMother = other.fTrackFilterMother;
	146	}
	147	return *this;
	148	}
a3b51fd2	149	*/
da97a08a	150	//____________________________________________________________________
	151	void AliAnalysisTaskMCParticleFilter::UserCreateOutputObjects()
	152	{
	153	// Create the output container
	154	if (OutputTree()&&fTrackFilterMother)
	155	OutputTree()->GetUserInfo()->Add(fTrackFilterMother);
d2740fba	156
	157	// this part is mainly needed to set the MCEventHandler
	158	// to the AODHandler, this will not be needed when
	159	// AODHandler goes to ANALYSISalice
	160	// setting in the steering macro will not work on proof :(
	161	// so we have to do it in a task
	162
	163	// the branch booking can also go into the AODHandler then
	164
dce1b636	165
dce1b636	166	// mcparticles
da97a08a	167	TClonesArray *tca = new TClonesArray("AliAODMCParticle", 0);
	168	tca->SetName(AliAODMCParticle::StdBranchName());
	169	AddAODBranch("TClonesArray",&tca);
da97a08a	170
dce1b636	171	// MC header...
dce1b636	172	AliAODMCHeader *mcHeader = new AliAODMCHeader();
dce1b636	173	mcHeader->SetName(AliAODMCHeader::StdBranchName());
	174	AddAODBranch("AliAODMCHeader",&mcHeader);
	175
	176
	177
	178	AliMCEventHandler mcH = (AliMCEventHandler) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
da97a08a	179	AliAODHandler aodH = dynamic_cast<AliAODHandler> ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
	180	if(!aodH){
	181	Printf("%s:&d Could not get AODHandler",(char*)__FILE__,__LINE__);
	182	return;
	183	}
	184	aodH->SetMCEventHandler(mcH);
dce1b636	185
a3b51fd2	186
	187	// these are histograms, for averaging and summing
	188	// do it via histograms to be PROOF-proof
	189	// which merges the results from different workers
	190	// histos are not saved reside only in memory
	191
	192
	193
	194	fHistList = new TList();
	195	TProfile *h1Xsec = new TProfile("h1Xsec","xsec from pyxsec.root",1,0,1);
	196	h1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
	197	fHistList->Add(h1Xsec);
	198	TH1F* h1Trials = new TH1F("h1Trials","trials from pyxsec.root",1,0,1);
	199	h1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
	200	fHistList->Add(h1Trials);
	201
da97a08a	202
	203	}
	204
	205	//____________________________________________________________________
	206	void AliAnalysisTaskMCParticleFilter::UserExec(Option_t /option*/)
	207	{
	208	// Execute analysis for current event
	209	//
	210
	211	// Fill AOD tracks from Kinematic stack
	212
	213	// get AliAOD Event
	214	AliAODEvent* aod = AODEvent();
	215	if (!aod) {
	216	AliWarning("No Output Handler connected, doing nothing !") ;
	217	return;
	218	}
	219
	220	AliMCEventHandler mcH = (AliMCEventHandler) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
	221	if(!mcH){
	222	AliWarning("No MC handler Found");
	223	return;
	224	}
	225
	226
	227	// fetch the output
	228	// Fill control histos
	229
	230	// tmp array for holding the mctracks
	231	// need to set mother and duagther __before__
	232	// filling in the tree
	233
	234	AliMCEvent *mcE = MCEvent();
	235
	236	if(!mcE){
	237	AliWarning("No MC event Found");
	238	return;
	239	}
	240
da97a08a	241	Int_t np = mcE->GetNumberOfTracks();
93836e1b	242	Int_t nprim = mcE->GetNumberOfPrimaries();
da97a08a	243	// TODO ADD MC VERTEX
	244
	245	// We take all real primaries
	246
	247
da97a08a	248	Int_t j=0;
	249	for (Int_t ip = 0; ip < np; ip++){
	250	AliMCParticle* mcpart = mcE->GetTrack(ip);
	251	TParticle* part = mcpart->Particle();
da97a08a	252	Float_t xv = part->Vx();
	253	Float_t yv = part->Vy();
	254	Float_t zv = part->Vz();
	255	Float_t rv = TMath::Sqrt(xv * xv + yv * yv);
	256
	257	Bool_t write = kFALSE;
63892dd0	258
da97a08a	259	if (ip < nprim) {
	260	// Select the primary event
	261	write = kTRUE;
	262	} else if (part->GetUniqueID() == 4) {
	263	// Particles from decay
	264	// Check that the decay chain ends at a primary particle
93836e1b	265	AliMCParticle* mother = mcpart;
93836e1b	266	Int_t imo = mcpart->GetMother();
da97a08a	267	while((imo >= nprim) && (mother->GetUniqueID() == 4)) {
93836e1b	268	mother = mcE->GetTrack(imo);
93836e1b	269	imo = mother->GetMother();
da97a08a	270	}
63892dd0	271	// Select according to pseudorapidity and production point of primary ancestor
93836e1b	272	if (imo < nprim && Select(mcE->GetTrack(imo)->Particle(), rv, zv))write = kTRUE;
da97a08a	273	} else if (part->GetUniqueID() == 5) {
da97a08a	274	// Now look for pair production
93836e1b	275	Int_t imo = mcpart->GetMother();
da97a08a	276	if (imo < nprim) {
	277	// Select, if the gamma is a primary
	278	write = kTRUE;
	279	} else {
	280	// Check if the gamma comes from the decay chain of a primary particle
93836e1b	281	AliMCParticle* mother = mcE->GetTrack(imo);
93836e1b	282	imo = mother->GetMother();
da97a08a	283	while((imo >= nprim) && (mother->GetUniqueID() == 4)) {
93836e1b	284	mother = mcE->GetTrack(imo);
93836e1b	285	imo = mother->GetMother();
da97a08a	286	}
da97a08a	287	// Select according to pseudorapidity and production point
93836e1b	288	if (imo < nprim && Select(mother->Particle(), rv, zv))
da97a08a	289	write = kTRUE;
	290	}
	291	}
	292	/*
	293	else if (part->GetUniqueID() == 13){
	294	// Evaporation
	295	// Check that we end at a primary particle
	296	TParticle* mother = part;
	297	Int_t imo = part->GetFirstMother();
	298	while((imo >= nprim) && ((mother->GetUniqueID() == 4 ) \|\| ( mother->GetUniqueID() == 13))) {
	299	mother = mcE->Stack()->Particle(imo);
	300	imo = mother->GetFirstMother();
	301	}
	302	// Select according to pseudorapidity and production point
	303	if (imo < nprim && Select(mother, rv, zv))
	304	write = kTRUE;
	305	}
	306	*/
	307	if (write) {
	308	if(mcH)mcH->SelectParticle(ip);
	309	j++;
	310	}
	311	}
	312
63892dd0	313
	314	// check for charm daughters
	315	static int iTaken = 0;
	316	static int iAll = 0;
	317	static int iCharm = 0;
	318	for (Int_t ip = 0; ip < np; ip++){
	319	AliMCParticle* mcpart = mcE->GetTrack(ip);
	320	TParticle* part = mcpart->Particle();
	321
	322	// if((TMath::Abs(part->GetPdgCode())/400)==1){
	323	if((TMath::Abs(part->GetPdgCode()))==411){
	324	// cases
	325	iCharm++;
93836e1b	326	Int_t d0 = mcpart->GetFirstDaughter();
93836e1b	327	Int_t d1 = mcpart->GetLastDaughter();
63892dd0	328	if(d0>0&&d1>0){
63892dd0	329	for(int id = d0;id <= d1;id++){
a3b51fd2	330	// AliMCParticle* daughter = mcE->GetTrack(id);
63892dd0	331	iAll++;
	332	if(mcH->IsParticleSelected(id))iTaken++;
	333	}
	334	}
	335	}
	336	}
a3c57861	337
a3c57861	338	AliInfo(Form("Taken daughters %d/%d of %d charm",iTaken,iAll,iCharm));
63892dd0	339
a3b51fd2	340	PostData(1,fHistList);
63892dd0	341
da97a08a	342	return;
	343	}
	344
a3b51fd2	345	void AliAnalysisTaskMCParticleFilter::Terminate(Option_t /option*/){
	346	//
	347	// Terminate the execution save the PYTHIA x_section to the UserInfo()
	348	//
	349
	350
	351	TTree tree = (TTree)GetOutputData(0);
	352	fHistList = (TList*)GetOutputData(1);
	353	if(!fHistList){
	354	Printf("%s:%d Output list not found",(char*)__FILE__,__LINE__);
	355	return;
	356	}
	357
	358	if(!tree){
	359	Printf("%s:%d Output tree not found",(char*)__FILE__,__LINE__);
	360	return;
	361	}
	362
	363	TProfile hXsec = ((TProfile)(fHistList->FindObject("h1Xsec")));
	364	TH1F* hTrials = ((TH1F*)(fHistList->FindObject("h1Trials")));
	365	if(!hXsec)return;
	366	if(!hTrials)return;
	367
	368	Float_t xsec = hXsec->GetBinContent(1);
	369	Float_t trials = hTrials->Integral();
	370	AliInfo(Form("Average -section %.4E and total number of trials %E",xsec,trials));
	371
	372	// This only works if the tree is not a special output
	373	// in this case it is already written
	374	/*
	375	TParameter<float> *x = new TParameter<float>("xsection",xsec);
	376	tree->GetUserInfo()->Add(x);
	377	TParameter<float> *t = new TParameter<float>("trials",trials);
	378	tree->GetUserInfo()->Add(t);
	379	*/
	380
	381
	382	}
	383
da97a08a	384	Bool_t AliAnalysisTaskMCParticleFilter::Select(TParticle* part, Float_t rv, Float_t zv)
	385	{
	386	// Selection accoring to eta of the mother and production point
	387	// This has to be refined !!!!!!
	388
	389	// Esp if we don't have collisison in the central barrel but e.g. beam gas
	390	// return kTRUE;
	391
	392	Float_t eta = part->Eta();
	393
	394	// central barrel consider everything in the ITS...
	395	// large eta window for smeared vertex and SPD acceptance (2 at z = 0)
	396	// larger for V0s in the TPC
	397	// if (TMath::Abs(eta) < 2.5 && rv < 250. && TMath::Abs(zv)<255)return kTRUE;
	398
63892dd0	399	if (TMath::Abs(eta) < 2.5 && rv < 170)return kTRUE;
63892dd0	400
da97a08a	401	// Andreas' Cuts
63892dd0	402	// if (TMath::Abs(eta) < 1. && rv < 170)return kTRUE;
	403
	404
da97a08a	405
	406	// Muon arm
	407	if(eta > -4.2 && eta < -2.3 && zv > -500.)return kTRUE; // Muon arms
	408
	409	// PMD acceptance 2.3 <= eta < = 3.5
	410	// if(eta>2.0&&eta<3.7)return kTRUE;
	411
	412	return kFALSE;
	413
	414	}
	415