[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 "TString.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)
{
  // Default constructor
}

//____________________________________________________________________
AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter(const char* name):
    AliAnalysisTaskSE(name),
    fTrackFilterMother(0x0)
{
  // Default constructor
}


//____________________________________________________________________
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();
    Printf("AODMCHeader %p",mcHeader);
    Printf("AODMCHeader ** %p",&mcHeader);
    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);

    
}

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

  AliStack* stack = mcE->Stack();
  Int_t np    = mcE->GetNumberOfTracks();
  Int_t nprim = stack->GetNprimary();
  // 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
      TParticle* mother = part;
      Int_t imo = part->GetFirstMother();
      while((imo >= nprim) && (mother->GetUniqueID() == 4)) {
	mother =  mcE->Stack()->Particle(imo);
	imo =  mother->GetFirstMother();
      }
      // Select according to pseudorapidity and production point of primary ancestor
      if (imo < nprim && Select(mcE->Stack()->Particle(imo), rv, zv))write = kTRUE;         
    } else if (part->GetUniqueID() == 5) {
      // Now look for pair production
      Int_t imo = part->GetFirstMother();
      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
	TParticle* mother = mcE->Stack()->Particle(imo);
	imo = mother->GetFirstMother();
	while((imo >= nprim) && (mother->GetUniqueID() == 4)) {
	  mother =  mcE->Stack()->Particle(imo);
	  imo =  mother->GetFirstMother();
	}
	// Select according to pseudorapidity and production point 
	if (imo < nprim && Select(mother, 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++;
      Printf("Decay Mother %s",part->GetPDG()->GetName());
      Int_t d0 =  part->GetFirstDaughter();
      Int_t d1 =  part->GetLastDaughter();
      if(d0>0&&d1>0){
	for(int id = d0;id <= d1;id++){
	  TParticle* daughter =  mcE->Stack()->Particle(id);
	  Printf("Decay Daughter %s",daughter->GetPDG()->GetName());
	  iAll++;
	  if(mcH->IsParticleSelected(id))iTaken++;
	}
      }
    }
  }
  Printf("Taken daughters %d/%d of %d charm",iTaken,iAll,iCharm);


  return;
}

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"
	24	#include "TString.h"
	25
	26	#include "AliAnalysisTaskMCParticleFilter.h"
	27	#include "AliAnalysisManager.h"
	28	#include "AliAnalysisFilter.h"
	29	#include "AliHeader.h"
	30	#include "AliStack.h"
	31	#include "AliMCEvent.h"
	32	#include "AliMCEventHandler.h"
	33	#include "AliAODEvent.h"
	34	#include "AliAODHeader.h"
dce1b636	35	#include "AliAODMCHeader.h"
da97a08a	36	#include "AliAODHandler.h"
	37	#include "AliAODVertex.h"
	38	#include "AliAODMCParticle.h"
	39
	40
	41	#include "AliLog.h"
	42
	43	ClassImp(AliAnalysisTaskMCParticleFilter)
	44
	45	////////////////////////////////////////////////////////////////////////
	46
	47	//____________________________________________________________________
	48	AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter():
	49	AliAnalysisTaskSE(),
	50	fTrackFilterMother(0x0)
	51	{
	52	// Default constructor
	53	}
	54
	55	//____________________________________________________________________
	56	AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter(const char* name):
	57	AliAnalysisTaskSE(name),
	58	fTrackFilterMother(0x0)
	59	{
	60	// Default constructor
	61	}
	62
	63
	64	//____________________________________________________________________
	65	AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter(const AliAnalysisTaskMCParticleFilter& obj):
	66	AliAnalysisTaskSE(obj),
	67	fTrackFilterMother(obj.fTrackFilterMother)
	68	{
	69	// Copy constructor
	70	}
	71
	72	//____________________________________________________________________
	73	AliAnalysisTaskMCParticleFilter::~AliAnalysisTaskMCParticleFilter()
	74	{
	75	// if( fTrackFilterMother ) delete fTrackFilterMother;
	76	}
	77
	78
	79	//____________________________________________________________________
	80	AliAnalysisTaskMCParticleFilter& AliAnalysisTaskMCParticleFilter::operator=(const AliAnalysisTaskMCParticleFilter& other)
	81	{
	82	// Assignment
	83	if(this!=&other) {
	84	AliAnalysisTaskSE::operator=(other);
	85	fTrackFilterMother = other.fTrackFilterMother;
	86	}
	87	return *this;
	88	}
	89
	90	//____________________________________________________________________
	91	void AliAnalysisTaskMCParticleFilter::UserCreateOutputObjects()
	92	{
	93	// Create the output container
	94	if (OutputTree()&&fTrackFilterMother)
	95	OutputTree()->GetUserInfo()->Add(fTrackFilterMother);
d2740fba	96
	97	// this part is mainly needed to set the MCEventHandler
	98	// to the AODHandler, this will not be needed when
	99	// AODHandler goes to ANALYSISalice
	100	// setting in the steering macro will not work on proof :(
	101	// so we have to do it in a task
	102
	103	// the branch booking can also go into the AODHandler then
	104
dce1b636	105
dce1b636	106	// mcparticles
da97a08a	107	TClonesArray *tca = new TClonesArray("AliAODMCParticle", 0);
	108	tca->SetName(AliAODMCParticle::StdBranchName());
	109	AddAODBranch("TClonesArray",&tca);
da97a08a	110
dce1b636	111	// MC header...
	112	AliAODMCHeader *mcHeader = new AliAODMCHeader();
	113	Printf("AODMCHeader %p",mcHeader);
	114	Printf("AODMCHeader ** %p",&mcHeader);
	115	mcHeader->SetName(AliAODMCHeader::StdBranchName());
	116	AddAODBranch("AliAODMCHeader",&mcHeader);
	117
	118
	119
	120	AliMCEventHandler mcH = (AliMCEventHandler) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
da97a08a	121	AliAODHandler aodH = dynamic_cast<AliAODHandler> ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
	122	if(!aodH){
	123	Printf("%s:&d Could not get AODHandler",(char*)__FILE__,__LINE__);
	124	return;
	125	}
	126	aodH->SetMCEventHandler(mcH);
dce1b636	127
da97a08a	128
	129	}
	130
	131	//____________________________________________________________________
	132	void AliAnalysisTaskMCParticleFilter::UserExec(Option_t /option*/)
	133	{
	134	// Execute analysis for current event
	135	//
	136
	137	// Fill AOD tracks from Kinematic stack
	138
	139	// get AliAOD Event
	140	AliAODEvent* aod = AODEvent();
	141	if (!aod) {
	142	AliWarning("No Output Handler connected, doing nothing !") ;
	143	return;
	144	}
	145
	146	AliMCEventHandler mcH = (AliMCEventHandler) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
	147	if(!mcH){
	148	AliWarning("No MC handler Found");
	149	return;
	150	}
	151
	152
	153	// fetch the output
	154	// Fill control histos
	155
	156	// tmp array for holding the mctracks
	157	// need to set mother and duagther __before__
	158	// filling in the tree
	159
	160	AliMCEvent *mcE = MCEvent();
	161
	162	if(!mcE){
	163	AliWarning("No MC event Found");
	164	return;
	165	}
	166
	167	AliStack* stack = mcE->Stack();
	168	Int_t np = mcE->GetNumberOfTracks();
	169	Int_t nprim = stack->GetNprimary();
	170	// TODO ADD MC VERTEX
	171
	172	// We take all real primaries
	173
	174
da97a08a	175	Int_t j=0;
	176	for (Int_t ip = 0; ip < np; ip++){
	177	AliMCParticle* mcpart = mcE->GetTrack(ip);
	178	TParticle* part = mcpart->Particle();
da97a08a	179	Float_t xv = part->Vx();
	180	Float_t yv = part->Vy();
	181	Float_t zv = part->Vz();
	182	Float_t rv = TMath::Sqrt(xv * xv + yv * yv);
	183
	184	Bool_t write = kFALSE;
63892dd0	185
da97a08a	186	if (ip < nprim) {
	187	// Select the primary event
	188	write = kTRUE;
	189	} else if (part->GetUniqueID() == 4) {
	190	// Particles from decay
	191	// Check that the decay chain ends at a primary particle
	192	TParticle* mother = part;
	193	Int_t imo = part->GetFirstMother();
	194	while((imo >= nprim) && (mother->GetUniqueID() == 4)) {
	195	mother = mcE->Stack()->Particle(imo);
	196	imo = mother->GetFirstMother();
	197	}
63892dd0	198	// Select according to pseudorapidity and production point of primary ancestor
63892dd0	199	if (imo < nprim && Select(mcE->Stack()->Particle(imo), rv, zv))write = kTRUE;
da97a08a	200	} else if (part->GetUniqueID() == 5) {
da97a08a	201	// Now look for pair production
da97a08a	202	Int_t imo = part->GetFirstMother();
	203	if (imo < nprim) {
	204	// Select, if the gamma is a primary
	205	write = kTRUE;
	206	} else {
	207	// Check if the gamma comes from the decay chain of a primary particle
	208	TParticle* mother = mcE->Stack()->Particle(imo);
	209	imo = mother->GetFirstMother();
	210	while((imo >= nprim) && (mother->GetUniqueID() == 4)) {
	211	mother = mcE->Stack()->Particle(imo);
	212	imo = mother->GetFirstMother();
	213	}
	214	// Select according to pseudorapidity and production point
	215	if (imo < nprim && Select(mother, rv, zv))
	216	write = kTRUE;
	217	}
	218	}
	219	/*
	220	else if (part->GetUniqueID() == 13){
	221	// Evaporation
	222	// Check that we end at a primary particle
	223	TParticle* mother = part;
	224	Int_t imo = part->GetFirstMother();
	225	while((imo >= nprim) && ((mother->GetUniqueID() == 4 ) \|\| ( mother->GetUniqueID() == 13))) {
	226	mother = mcE->Stack()->Particle(imo);
	227	imo = mother->GetFirstMother();
	228	}
	229	// Select according to pseudorapidity and production point
	230	if (imo < nprim && Select(mother, rv, zv))
	231	write = kTRUE;
	232	}
	233	*/
	234	if (write) {
	235	if(mcH)mcH->SelectParticle(ip);
	236	j++;
	237	}
	238	}
	239
63892dd0	240
	241	// check for charm daughters
	242	static int iTaken = 0;
	243	static int iAll = 0;
	244	static int iCharm = 0;
	245	for (Int_t ip = 0; ip < np; ip++){
	246	AliMCParticle* mcpart = mcE->GetTrack(ip);
	247	TParticle* part = mcpart->Particle();
	248
	249	// if((TMath::Abs(part->GetPdgCode())/400)==1){
	250	if((TMath::Abs(part->GetPdgCode()))==411){
	251	// cases
	252	iCharm++;
	253	Printf("Decay Mother %s",part->GetPDG()->GetName());
	254	Int_t d0 = part->GetFirstDaughter();
	255	Int_t d1 = part->GetLastDaughter();
	256	if(d0>0&&d1>0){
	257	for(int id = d0;id <= d1;id++){
	258	TParticle* daughter = mcE->Stack()->Particle(id);
	259	Printf("Decay Daughter %s",daughter->GetPDG()->GetName());
	260	iAll++;
	261	if(mcH->IsParticleSelected(id))iTaken++;
	262	}
	263	}
	264	}
	265	}
	266	Printf("Taken daughters %d/%d of %d charm",iTaken,iAll,iCharm);
	267
	268
da97a08a	269	return;
	270	}
	271
	272	Bool_t AliAnalysisTaskMCParticleFilter::Select(TParticle* part, Float_t rv, Float_t zv)
	273	{
	274	// Selection accoring to eta of the mother and production point
	275	// This has to be refined !!!!!!
	276
	277	// Esp if we don't have collisison in the central barrel but e.g. beam gas
	278	// return kTRUE;
	279
	280	Float_t eta = part->Eta();
	281
	282	// central barrel consider everything in the ITS...
	283	// large eta window for smeared vertex and SPD acceptance (2 at z = 0)
	284	// larger for V0s in the TPC
	285	// if (TMath::Abs(eta) < 2.5 && rv < 250. && TMath::Abs(zv)<255)return kTRUE;
	286
63892dd0	287	if (TMath::Abs(eta) < 2.5 && rv < 170)return kTRUE;
63892dd0	288
da97a08a	289	// Andreas' Cuts
63892dd0	290	// if (TMath::Abs(eta) < 1. && rv < 170)return kTRUE;
	291
	292
da97a08a	293
	294	// Muon arm
	295	if(eta > -4.2 && eta < -2.3 && zv > -500.)return kTRUE; // Muon arms
	296
	297	// PMD acceptance 2.3 <= eta < = 3.5
	298	// if(eta>2.0&&eta<3.7)return kTRUE;
	299
	300	return kFALSE;
	301
	302	}
	303