[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 "AliESDInputHandler.h"
#include "AliAODEvent.h"
#include "AliAODHeader.h"
#include "AliAODMCHeader.h"
#include "AliAODHandler.h"
#include "AliAODVertex.h"
#include "AliAODMCParticle.h"
#include "AliCollisionGeometry.h"
#include "AliGenDPMjetEventHeader.h"
#include "AliGenHijingEventHeader.h"
#include "AliGenPythiaEventHeader.h"
#include "AliGenCocktailEventHeader.h"
#include "AliGenEventHeaderTunedPbPb.h"
#include "AliESDtrack.h"
#include "AliAODTrack.h"
#include "AliAODPid.h"
#include "AliESDpid.h"

#include "AliLog.h"

ClassImp(AliAnalysisTaskMCParticleFilter)

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

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

Bool_t AliAnalysisTaskMCParticleFilter::Notify()
{
  //
  // Implemented Notify() to read the cross sections
  // from pyxsec.root
  // 
  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
  TTree *tree = mgr->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());
    TString datafile = mgr->GetInputEventHandler()->GetInputFileName();
    if (fileName.Contains(datafile)) {
        fileName.ReplaceAll(datafile, "pyxsec.root");
    }
    else 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){
      AliWarning(Form("%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);
  }
  return kTRUE;
}


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

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

  if(fAODMcHeader){
    delete fAODMcHeader;
  }
  if(fAODMcParticles){
    fAODMcParticles->Delete();
    delete fAODMcParticles;
  }
}

/*
//____________________________________________________________________
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
    fAODMcParticles = new TClonesArray("AliAODMCParticle", 0);
    fAODMcParticles->SetName(AliAODMCParticle::StdBranchName());
    AddAODBranch("TClonesArray",&fAODMcParticles);

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

    AliMCEventHandler *mcH = (AliMCEventHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler()); 
    AliAODHandler *aodH = dynamic_cast<AliAODHandler*> ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
    if(!aodH){
      AliWarning("Could not get AODHandler");
      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 MC header",1,0,1);
    h1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
    fHistList->Add(h1Trials);

    PostData(1,fHistList);
    
}

//____________________________________________________________________
void AliAnalysisTaskMCParticleFilter::UserExec(Option_t */*option*/)
{
  // Execute analysis for current event
  //
  
  // Fill AOD tracks from Kinematic stack
  PostData(1,fHistList);
  
  // 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;
  }
  
  AliAODHandler *aodH = dynamic_cast<AliAODHandler*> ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
  if(!aodH){
    AliWarning("Could not get AODHandler");
    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

  // Get the proper MC Collision Geometry
  AliGenEventHeader* mcEH = mcE->GenEventHeader();

  AliGenPythiaEventHeader *pyH  = dynamic_cast<AliGenPythiaEventHeader*>(mcEH);
  AliGenHijingEventHeader *hiH  = 0;
  AliCollisionGeometry    *colG = 0;
  AliGenDPMjetEventHeader *dpmH = 0;
  AliGenEventHeaderTunedPbPb *tunedH = 0;

  // it can be only one save some casts
  // assuming PYTHIA and HIJING are the most likely ones...
  if(!pyH){
      hiH = dynamic_cast<AliGenHijingEventHeader*>(mcEH);
      if(!hiH){
	  dpmH = dynamic_cast<AliGenDPMjetEventHeader*>(mcEH);
	  if(!dpmH){
	    tunedH = dynamic_cast<AliGenEventHeaderTunedPbPb*>(mcEH);
	  }
      }
  }
  
  if (hiH || dpmH) colG = dynamic_cast<AliCollisionGeometry*>(mcEH);

  // fetch the trials on a event by event basis, not from pyxsec.root otherwise 
  // we will get a problem when running on proof since Notify may be called 
  // more than once per file
  // consider storing this information in the AOD output via AliAODHandler
  Float_t ntrials = 0;
  if (!colG) {
    AliGenCocktailEventHeader *ccEH = dynamic_cast<AliGenCocktailEventHeader *>(mcEH);
    if (ccEH) {
      TList *genHeaders = ccEH->GetHeaders();
      for (int imch=0; imch<genHeaders->GetEntries(); imch++) {
	if(!pyH)pyH = dynamic_cast<AliGenPythiaEventHeader*>(genHeaders->At(imch));
	if(!hiH)hiH = dynamic_cast<AliGenHijingEventHeader*>(genHeaders->At(imch));
	if(!colG)colG = dynamic_cast<AliCollisionGeometry *>(genHeaders->At(imch));
	if(!dpmH)dpmH = dynamic_cast<AliGenDPMjetEventHeader*>(genHeaders->At(imch));
      }
    }
  }

  // take the trials from the p+p event
  if(hiH)ntrials = hiH->Trials();
  if(dpmH)ntrials = dpmH->Trials();
  if(pyH)ntrials = pyH->Trials();
  if(ntrials)((TH1F*)(fHistList->FindObject("h1Trials")))->Fill("#sum{ntrials}",ntrials); 
  

  if (colG) {
    fAODMcHeader->SetReactionPlaneAngle(colG->ReactionPlaneAngle());
    AliInfo(Form("Found Collision Geometry. Got Reaction Plane %lf\n", colG->ReactionPlaneAngle()));
  }
  else if(tunedH) {
    fAODMcHeader->SetReactionPlaneAngle(tunedH->GetPsi2());
    fAODMcHeader->SetCrossSection(tunedH->GetCentrality());
  }


  Int_t j=0;
  for (Int_t ip = 0; ip < np; ip++){
    AliMCParticle* mcpart = (AliMCParticle*) 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() == kPDecay) {
      // 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() == kPDecay)) {
	mother =  (AliMCParticle*) mcE->GetTrack(imo);
	imo =  mother->GetMother();
      }
      // Select according to pseudorapidity and production point of primary ancestor
      if (imo < nprim)write = kTRUE;         
      // if(!Select(((AliMCParticle*) mcE->GetTrack(imo))->Particle(), rv, zv))write = kFALSE; // selection on eta and phi of the mother
    } else if (part->GetUniqueID() == kPPair) {
      // 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 =  (AliMCParticle*) mcE->GetTrack(imo);
	imo = mother->GetMother();
	while((imo >= nprim) && (mother->GetUniqueID() == kPDecay)) {
	  mother =   (AliMCParticle*) mcE->GetTrack(imo);
	  imo =  mother->GetMother();
	}
	// Select according to pseudorapidity and production point 
	if (imo < nprim && Select(mother->Particle(), rv, zv)) 
	  write = kTRUE;
      }
    }

    if (write) {
      if(mcH)mcH->SelectParticle(ip);
      j++;      
    }
  }

  /*

  // check varibales for charm need all daughters
  static int  iTaken = 0;
  static int  iAll = 0;
  static int  iCharm = 0;

  for (Int_t ip = 0; ip < np; ip++){
    AliMCParticle* mcpart = (AliMCParticle*) mcE->GetTrack(ip);
    TParticle* part = mcpart->Particle();
    // debug info to check fro charm daugthers
    if((TMath::Abs(part->GetPdgCode()))==411){
      iCharm++;
      Int_t d0 =  part->GetFirstDaughter();
      Int_t d1 =  part->GetLastDaughter();
      if(d0>0&&d1>0){
	for(int id = d0;id <= d1;id++){
	  iAll++;
	  if(mcH->IsParticleSelected(id)){
	    iTaken++;
	    Printf("> %d/%d Taken",id,nprim);
	    PrintMCParticle( (AliMCParticle*)mcE->GetTrack(id),id);
	  }
	  else{
	    Printf("> %d/%d NOT Taken",id,nprim);
	    PrintMCParticle( (AliMCParticle*)mcE->GetTrack(id),id);
	  }
	}
      }
    }// if charm
    AliInfo(Form("Taken daughters %d/%d of %d charm",iTaken,iAll,iCharm));
  }
  */


  aodH->StoreMCParticles();

  return;
}

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


  fHistList = (TList*)GetOutputData(1);
  if(!fHistList){
    Printf("%s:%d Output list 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));

}

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;
 
}

void AliAnalysisTaskMCParticleFilter::PrintMCParticle(const AliMCParticle *mcp,Int_t np){
  
  const TParticle *p = mcp->Particle();
  Printf("Nr.%d --- Status %d ---- Mother1 %d Mother2 %d Daughter1 %d Daughte\
r2 %d ",
	 np,p->GetStatusCode(),p->GetMother(0),p->GetMother(1),p->GetDaughter \
	 (0),p->GetDaughter(1));
  Printf("Eta %3.3f Phi %3.3f ",p->Eta(),p->Phi());
  p->Print();
  Printf("---------------------------------------");
}
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"
fdbaa4ce	38	#include "AliESDInputHandler.h"
da97a08a	39	#include "AliAODEvent.h"
da97a08a	40	#include "AliAODHeader.h"
dce1b636	41	#include "AliAODMCHeader.h"
da97a08a	42	#include "AliAODHandler.h"
	43	#include "AliAODVertex.h"
	44	#include "AliAODMCParticle.h"
e988332e	45	#include "AliCollisionGeometry.h"
c4b2bb17	46	#include "AliGenDPMjetEventHeader.h"
	47	#include "AliGenHijingEventHeader.h"
	48	#include "AliGenPythiaEventHeader.h"
e988332e	49	#include "AliGenCocktailEventHeader.h"
995e33f0	50	#include "AliGenEventHeaderTunedPbPb.h"
63825dc4	51	#include "AliESDtrack.h"
	52	#include "AliAODTrack.h"
	53	#include "AliAODPid.h"
	54	#include "AliESDpid.h"
da97a08a	55
	56	#include "AliLog.h"
	57
	58	ClassImp(AliAnalysisTaskMCParticleFilter)
	59
	60	////////////////////////////////////////////////////////////////////////
	61
	62	//____________________________________________________________________
	63	AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter():
3189f57d	64	AliAnalysisTaskSE(),
a3b51fd2	65	fTrackFilterMother(0x0),
3189f57d	66	fAODMcHeader(0x0),
3189f57d	67	fAODMcParticles(0x0),
a3b51fd2	68	fHistList(0x0)
da97a08a	69	{
	70	// Default constructor
	71	}
	72
a3b51fd2	73	Bool_t AliAnalysisTaskMCParticleFilter::Notify()
	74	{
	75	//
	76	// Implemented Notify() to read the cross sections
c4b2bb17	77	// from pyxsec.root
a3b51fd2	78	//
fdbaa4ce	79	AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
fdbaa4ce	80	TTree *tree = mgr->GetTree();
a3b51fd2	81	Double_t xsection = 0;
	82	UInt_t ntrials = 0;
	83	if(tree){
	84	TFile *curfile = tree->GetCurrentFile();
	85	if (!curfile) {
	86	Error("Notify","No current file");
	87	return kFALSE;
	88	}
	89
	90	TString fileName(curfile->GetName());
fdbaa4ce	91	TString datafile = mgr->GetInputEventHandler()->GetInputFileName();
	92	if (fileName.Contains(datafile)) {
	93	fileName.ReplaceAll(datafile, "pyxsec.root");
	94	}
	95	else if(fileName.Contains("AliESDs.root")){
a3b51fd2	96	fileName.ReplaceAll("AliESDs.root", "pyxsec.root");
	97	}
	98	else if(fileName.Contains("AliAOD.root")){
	99	fileName.ReplaceAll("AliAOD.root", "pyxsec.root");
	100	}
	101	else if(fileName.Contains("galice.root")){
	102	// for running with galice and kinematics alone...
	103	fileName.ReplaceAll("galice.root", "pyxsec.root");
	104	}
	105
	106
	107	TFile *fxsec = TFile::Open(fileName.Data());
	108	if(!fxsec){
	109	AliInfo(Form("%s:%d %s not found in the Input",(char*)__FILE__,__LINE__,fileName.Data()));
	110	// not a severe condition we just do not have the information...
	111	return kTRUE;
	112	}
	113	TTree xtree = (TTree)fxsec->Get("Xsection");
	114	if(!xtree){
cd0ecb97	115	AliWarning(Form("%s:%d tree not found in the pyxsec.root",(char*)__FILE__,__LINE__));
a3b51fd2	116	return kTRUE;
	117	}
	118	xtree->SetBranchAddress("xsection",&xsection);
	119	xtree->SetBranchAddress("ntrials",&ntrials);
	120	xtree->GetEntry(0);
	121	((TProfile*)(fHistList->FindObject("h1Xsec")))->Fill("<#sigma>",xsection);
a3b51fd2	122	}
	123	return kTRUE;
	124	}
	125
	126
	127
da97a08a	128	//____________________________________________________________________
	129	AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter(const char* name):
	130	AliAnalysisTaskSE(name),
a3b51fd2	131	fTrackFilterMother(0x0),
3189f57d	132	fAODMcHeader(0x0),
3189f57d	133	fAODMcParticles(0x0),
a3b51fd2	134	fHistList(0x0)
da97a08a	135	{
da97a08a	136	// Default constructor
cc008a4b	137	DefineOutput(1, TList::Class());
da97a08a	138	}
da97a08a	139
a3b51fd2	140	/*
da97a08a	141	//____________________________________________________________________
	142	AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter(const AliAnalysisTaskMCParticleFilter& obj):
	143	AliAnalysisTaskSE(obj),
	144	fTrackFilterMother(obj.fTrackFilterMother)
	145	{
	146	// Copy constructor
	147	}
a3b51fd2	148	*/
da97a08a	149	//____________________________________________________________________
	150	AliAnalysisTaskMCParticleFilter::~AliAnalysisTaskMCParticleFilter()
	151	{
3189f57d	152
	153	if(fAODMcHeader){
	154	delete fAODMcHeader;
	155	}
	156	if(fAODMcParticles){
	157	fAODMcParticles->Delete();
	158	delete fAODMcParticles;
	159	}
da97a08a	160	}
da97a08a	161
a3b51fd2	162	/*
da97a08a	163	//____________________________________________________________________
	164	AliAnalysisTaskMCParticleFilter& AliAnalysisTaskMCParticleFilter::operator=(const AliAnalysisTaskMCParticleFilter& other)
	165	{
	166	// Assignment
	167	if(this!=&other) {
	168	AliAnalysisTaskSE::operator=(other);
	169	fTrackFilterMother = other.fTrackFilterMother;
	170	}
	171	return *this;
	172	}
a3b51fd2	173	*/
da97a08a	174	//____________________________________________________________________
	175	void AliAnalysisTaskMCParticleFilter::UserCreateOutputObjects()
	176	{
	177	// Create the output container
5c66ad80	178
a7cd6336	179
da97a08a	180	if (OutputTree()&&fTrackFilterMother)
da97a08a	181	OutputTree()->GetUserInfo()->Add(fTrackFilterMother);
d2740fba	182
	183	// this part is mainly needed to set the MCEventHandler
	184	// to the AODHandler, this will not be needed when
	185	// AODHandler goes to ANALYSISalice
	186	// setting in the steering macro will not work on proof :(
	187	// so we have to do it in a task
	188
	189	// the branch booking can also go into the AODHandler then
	190
dce1b636	191
dce1b636	192	// mcparticles
3189f57d	193	fAODMcParticles = new TClonesArray("AliAODMCParticle", 0);
	194	fAODMcParticles->SetName(AliAODMCParticle::StdBranchName());
	195	AddAODBranch("TClonesArray",&fAODMcParticles);
da97a08a	196
dce1b636	197	// MC header...
3189f57d	198	fAODMcHeader = new AliAODMCHeader();
	199	fAODMcHeader->SetName(AliAODMCHeader::StdBranchName());
	200	AddAODBranch("AliAODMCHeader",&fAODMcHeader);
dce1b636	201
dce1b636	202	AliMCEventHandler mcH = (AliMCEventHandler) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
da97a08a	203	AliAODHandler aodH = dynamic_cast<AliAODHandler> ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
da97a08a	204	if(!aodH){
3827f93a	205	AliWarning("Could not get AODHandler");
da97a08a	206	return;
	207	}
	208	aodH->SetMCEventHandler(mcH);
dce1b636	209
a3b51fd2	210
	211	// these are histograms, for averaging and summing
	212	// do it via histograms to be PROOF-proof
	213	// which merges the results from different workers
	214	// histos are not saved reside only in memory
	215
	216
	217
	218	fHistList = new TList();
	219	TProfile *h1Xsec = new TProfile("h1Xsec","xsec from pyxsec.root",1,0,1);
	220	h1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
	221	fHistList->Add(h1Xsec);
c4b2bb17	222	TH1F* h1Trials = new TH1F("h1Trials","trials from MC header",1,0,1);
a3b51fd2	223	h1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
	224	fHistList->Add(h1Trials);
	225
5c66ad80	226	PostData(1,fHistList);
da97a08a	227
	228	}
	229
	230	//____________________________________________________________________
	231	void AliAnalysisTaskMCParticleFilter::UserExec(Option_t /option*/)
	232	{
5c66ad80	233	// Execute analysis for current event
	234	//
	235
	236	// Fill AOD tracks from Kinematic stack
	237	PostData(1,fHistList);
	238
da97a08a	239	// get AliAOD Event
	240	AliAODEvent* aod = AODEvent();
	241	if (!aod) {
	242	AliWarning("No Output Handler connected, doing nothing !") ;
	243	return;
	244	}
	245
	246	AliMCEventHandler mcH = (AliMCEventHandler) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
	247	if(!mcH){
	248	AliWarning("No MC handler Found");
	249	return;
	250	}
3827f93a	251
	252	AliAODHandler aodH = dynamic_cast<AliAODHandler> ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
	253	if(!aodH){
	254	AliWarning("Could not get AODHandler");
	255	return;
	256	}
	257
da97a08a	258
	259
	260	// fetch the output
	261	// Fill control histos
	262
	263	// tmp array for holding the mctracks
	264	// need to set mother and duagther __before__
	265	// filling in the tree
	266
	267	AliMCEvent *mcE = MCEvent();
	268
	269	if(!mcE){
	270	AliWarning("No MC event Found");
	271	return;
	272	}
	273
da97a08a	274	Int_t np = mcE->GetNumberOfTracks();
93836e1b	275	Int_t nprim = mcE->GetNumberOfPrimaries();
da97a08a	276	// TODO ADD MC VERTEX
e988332e	277
e988332e	278	// Get the proper MC Collision Geometry
e988332e	279	AliGenEventHeader* mcEH = mcE->GenEventHeader();
e988332e	280
cc008a4b	281	AliGenPythiaEventHeader pyH = dynamic_cast<AliGenPythiaEventHeader>(mcEH);
	282	AliGenHijingEventHeader *hiH = 0;
	283	AliCollisionGeometry *colG = 0;
c4b2bb17	284	AliGenDPMjetEventHeader *dpmH = 0;
995e33f0	285	AliGenEventHeaderTunedPbPb *tunedH = 0;
995e33f0	286
c4b2bb17	287	// it can be only one save some casts
	288	// assuming PYTHIA and HIJING are the most likely ones...
	289	if(!pyH){
cc008a4b	290	hiH = dynamic_cast<AliGenHijingEventHeader*>(mcEH);
	291	if(!hiH){
	292	dpmH = dynamic_cast<AliGenDPMjetEventHeader*>(mcEH);
995e33f0	293	if(!dpmH){
	294	tunedH = dynamic_cast<AliGenEventHeaderTunedPbPb*>(mcEH);
	295	}
cc008a4b	296	}
c4b2bb17	297	}
cc008a4b	298
cc008a4b	299	if (hiH \|\| dpmH) colG = dynamic_cast<AliCollisionGeometry*>(mcEH);
c4b2bb17	300
	301	// fetch the trials on a event by event basis, not from pyxsec.root otherwise
	302	// we will get a problem when running on proof since Notify may be called
	303	// more than once per file
	304	// consider storing this information in the AOD output via AliAODHandler
	305	Float_t ntrials = 0;
e988332e	306	if (!colG) {
	307	AliGenCocktailEventHeader ccEH = dynamic_cast<AliGenCocktailEventHeader >(mcEH);
	308	if (ccEH) {
	309	TList *genHeaders = ccEH->GetHeaders();
	310	for (int imch=0; imch<genHeaders->GetEntries(); imch++) {
04152663	311	if(!pyH)pyH = dynamic_cast<AliGenPythiaEventHeader*>(genHeaders->At(imch));
04152663	312	if(!hiH)hiH = dynamic_cast<AliGenHijingEventHeader*>(genHeaders->At(imch));
c4b2bb17	313	if(!colG)colG = dynamic_cast<AliCollisionGeometry *>(genHeaders->At(imch));
04152663	314	if(!dpmH)dpmH = dynamic_cast<AliGenDPMjetEventHeader*>(genHeaders->At(imch));
e988332e	315	}
	316	}
	317	}
	318
c4b2bb17	319	// take the trials from the p+p event
c4b2bb17	320	if(hiH)ntrials = hiH->Trials();
5ae590c4	321	if(dpmH)ntrials = dpmH->Trials();
c4b2bb17	322	if(pyH)ntrials = pyH->Trials();
	323	if(ntrials)((TH1F*)(fHistList->FindObject("h1Trials")))->Fill("#sum{ntrials}",ntrials);
	324
	325
	326
	327
e988332e	328	if (colG) {
3189f57d	329	fAODMcHeader->SetReactionPlaneAngle(colG->ReactionPlaneAngle());
cd0ecb97	330	AliInfo(Form("Found Collision Geometry. Got Reaction Plane %lf\n", colG->ReactionPlaneAngle()));
e988332e	331	}
995e33f0	332	else if(tunedH) {
	333	fAODMcHeader->SetReactionPlaneAngle(tunedH->GetPsi2());
	334	fAODMcHeader->SetCrossSection(tunedH->GetCentrality());
	335	}
de514cf0	336
da97a08a	337
da97a08a	338	Int_t j=0;
da97a08a	339	for (Int_t ip = 0; ip < np; ip++){
7aad0c47	340	AliMCParticle* mcpart = (AliMCParticle*) mcE->GetTrack(ip);
da97a08a	341	TParticle* part = mcpart->Particle();
da97a08a	342	Float_t xv = part->Vx();
	343	Float_t yv = part->Vy();
	344	Float_t zv = part->Vz();
	345	Float_t rv = TMath::Sqrt(xv * xv + yv * yv);
	346
	347	Bool_t write = kFALSE;
63892dd0	348
da97a08a	349	if (ip < nprim) {
	350	// Select the primary event
	351	write = kTRUE;
43d79481	352	} else if (part->GetUniqueID() == kPDecay) {
da97a08a	353	// Particles from decay
da97a08a	354	// Check that the decay chain ends at a primary particle
93836e1b	355	AliMCParticle* mother = mcpart;
93836e1b	356	Int_t imo = mcpart->GetMother();
43d79481	357	while((imo >= nprim) && (mother->GetUniqueID() == kPDecay)) {
7aad0c47	358	mother = (AliMCParticle*) mcE->GetTrack(imo);
93836e1b	359	imo = mother->GetMother();
da97a08a	360	}
63892dd0	361	// Select according to pseudorapidity and production point of primary ancestor
c7ef8308	362	if (imo < nprim)write = kTRUE;
c7ef8308	363	// if(!Select(((AliMCParticle*) mcE->GetTrack(imo))->Particle(), rv, zv))write = kFALSE; // selection on eta and phi of the mother
43d79481	364	} else if (part->GetUniqueID() == kPPair) {
da97a08a	365	// Now look for pair production
93836e1b	366	Int_t imo = mcpart->GetMother();
da97a08a	367	if (imo < nprim) {
	368	// Select, if the gamma is a primary
	369	write = kTRUE;
	370	} else {
	371	// Check if the gamma comes from the decay chain of a primary particle
7aad0c47	372	AliMCParticle* mother = (AliMCParticle*) mcE->GetTrack(imo);
93836e1b	373	imo = mother->GetMother();
43d79481	374	while((imo >= nprim) && (mother->GetUniqueID() == kPDecay)) {
7aad0c47	375	mother = (AliMCParticle*) mcE->GetTrack(imo);
93836e1b	376	imo = mother->GetMother();
da97a08a	377	}
da97a08a	378	// Select according to pseudorapidity and production point
93836e1b	379	if (imo < nprim && Select(mother->Particle(), rv, zv))
da97a08a	380	write = kTRUE;
	381	}
	382	}
c7ef8308	383
da97a08a	384	if (write) {
da97a08a	385	if(mcH)mcH->SelectParticle(ip);
c7ef8308	386	j++;
	387	}
	388	}
	389
	390	/*
b7aedb79	391
	392	// check varibales for charm need all daughters
	393	static int iTaken = 0;
	394	static int iAll = 0;
	395	static int iCharm = 0;
	396
c7ef8308	397	for (Int_t ip = 0; ip < np; ip++){
	398	AliMCParticle* mcpart = (AliMCParticle*) mcE->GetTrack(ip);
	399	TParticle* part = mcpart->Particle();
	400	// debug info to check fro charm daugthers
	401	if((TMath::Abs(part->GetPdgCode()))==411){
	402	iCharm++;
	403	Int_t d0 = part->GetFirstDaughter();
	404	Int_t d1 = part->GetLastDaughter();
	405	if(d0>0&&d1>0){
	406	for(int id = d0;id <= d1;id++){
	407	iAll++;
	408	if(mcH->IsParticleSelected(id)){
	409	iTaken++;
	410	Printf("> %d/%d Taken",id,nprim);
	411	PrintMCParticle( (AliMCParticle*)mcE->GetTrack(id),id);
	412	}
	413	else{
	414	Printf("> %d/%d NOT Taken",id,nprim);
	415	PrintMCParticle( (AliMCParticle*)mcE->GetTrack(id),id);
de514cf0	416	}
63892dd0	417	}
c7ef8308	418	}
	419	}// if charm
	420	AliInfo(Form("Taken daughters %d/%d of %d charm",iTaken,iAll,iCharm));
63892dd0	421	}
c7ef8308	422	*/
a3c57861	423
63892dd0	424
3827f93a	425	aodH->StoreMCParticles();
63892dd0	426
da97a08a	427	return;
	428	}
	429
a3b51fd2	430	void AliAnalysisTaskMCParticleFilter::Terminate(Option_t /option*/){
	431	//
	432	// Terminate the execution save the PYTHIA x_section to the UserInfo()
	433	//
	434
	435
a3b51fd2	436	fHistList = (TList*)GetOutputData(1);
	437	if(!fHistList){
	438	Printf("%s:%d Output list not found",(char*)__FILE__,__LINE__);
	439	return;
	440	}
	441
a3b51fd2	442	TProfile hXsec = ((TProfile)(fHistList->FindObject("h1Xsec")));
	443	TH1F* hTrials = ((TH1F*)(fHistList->FindObject("h1Trials")));
	444	if(!hXsec)return;
	445	if(!hTrials)return;
	446
	447	Float_t xsec = hXsec->GetBinContent(1);
	448	Float_t trials = hTrials->Integral();
	449	AliInfo(Form("Average -section %.4E and total number of trials %E",xsec,trials));
a3b51fd2	450
	451	}
	452
da97a08a	453	Bool_t AliAnalysisTaskMCParticleFilter::Select(TParticle* part, Float_t rv, Float_t zv)
	454	{
	455	// Selection accoring to eta of the mother and production point
	456	// This has to be refined !!!!!!
	457
	458	// Esp if we don't have collisison in the central barrel but e.g. beam gas
	459	// return kTRUE;
	460
	461	Float_t eta = part->Eta();
	462
	463	// central barrel consider everything in the ITS...
	464	// large eta window for smeared vertex and SPD acceptance (2 at z = 0)
	465	// larger for V0s in the TPC
	466	// if (TMath::Abs(eta) < 2.5 && rv < 250. && TMath::Abs(zv)<255)return kTRUE;
	467
63892dd0	468	if (TMath::Abs(eta) < 2.5 && rv < 170)return kTRUE;
63892dd0	469
da97a08a	470	// Andreas' Cuts
63892dd0	471	// if (TMath::Abs(eta) < 1. && rv < 170)return kTRUE;
	472
	473
da97a08a	474
	475	// Muon arm
	476	if(eta > -4.2 && eta < -2.3 && zv > -500.)return kTRUE; // Muon arms
	477
	478	// PMD acceptance 2.3 <= eta < = 3.5
	479	// if(eta>2.0&&eta<3.7)return kTRUE;
	480
	481	return kFALSE;
	482
	483	}
	484
c7ef8308	485	void AliAnalysisTaskMCParticleFilter::PrintMCParticle(const AliMCParticle *mcp,Int_t np){
	486
	487	const TParticle *p = mcp->Particle();
	488	Printf("Nr.%d --- Status %d ---- Mother1 %d Mother2 %d Daughter1 %d Daughte\
	489	r2 %d ",
	490	np,p->GetStatusCode(),p->GetMother(0),p->GetMother(1),p->GetDaughter \
	491	(0),p->GetDaughter(1));
	492	Printf("Eta %3.3f Phi %3.3f ",p->Eta(),p->Phi());
	493	p->Print();
	494	Printf("---------------------------------------");
	495	}