X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=ANALYSIS%2FAliAnalysisTaskMCParticleFilter.cxx;h=1dc4f80fd4b266614a14eb028932e27b0b45f0d6;hb=d8a0be371a169cec5e84a3003c6973c481ac3e71;hp=2a13ecac2a098df6020c046e0aca8c44ecb41f84;hpb=7aad0c47882abec2d34b23d11dfb71c1bb2a3c7b;p=u%2Fmrichter%2FAliRoot.git diff --git a/ANALYSIS/AliAnalysisTaskMCParticleFilter.cxx b/ANALYSIS/AliAnalysisTaskMCParticleFilter.cxx index 2a13ecac2a0..1dc4f80fd4b 100644 --- a/ANALYSIS/AliAnalysisTaskMCParticleFilter.cxx +++ b/ANALYSIS/AliAnalysisTaskMCParticleFilter.cxx @@ -35,13 +35,18 @@ #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 "AliLog.h" @@ -51,8 +56,10 @@ ClassImp(AliAnalysisTaskMCParticleFilter) //____________________________________________________________________ AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter(): -AliAnalysisTaskSE(), + AliAnalysisTaskSE(), fTrackFilterMother(0x0), + fAODMcHeader(0x0), + fAODMcParticles(0x0), fHistList(0x0) { // Default constructor @@ -62,9 +69,10 @@ Bool_t AliAnalysisTaskMCParticleFilter::Notify() { // // Implemented Notify() to read the cross sections - // and number of trials from pyxsec.root + // from pyxsec.root // - TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree(); + AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); + TTree *tree = mgr->GetTree(); Double_t xsection = 0; UInt_t ntrials = 0; if(tree){ @@ -75,7 +83,11 @@ Bool_t AliAnalysisTaskMCParticleFilter::Notify() } TString fileName(curfile->GetName()); - if(fileName.Contains("AliESDs.root")){ + 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")){ @@ -95,15 +107,13 @@ Bool_t AliAnalysisTaskMCParticleFilter::Notify() } TTree *xtree = (TTree*)fxsec->Get("Xsection"); if(!xtree){ - Printf("%s:%d tree not found in the pyxsec.root",(char*)__FILE__,__LINE__); + 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); - ((TH1F*)(fHistList->FindObject("h1Trials")))->Fill("#sum{ntrials}",ntrials); - } return kTRUE; } @@ -114,10 +124,12 @@ Bool_t AliAnalysisTaskMCParticleFilter::Notify() AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter(const char* name): AliAnalysisTaskSE(name), fTrackFilterMother(0x0), + fAODMcHeader(0x0), + fAODMcParticles(0x0), fHistList(0x0) { // Default constructor - DefineOutput(1,TList::Class()); + DefineOutput(1, TList::Class()); } /* @@ -132,7 +144,14 @@ AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter(const AliAnalys //____________________________________________________________________ AliAnalysisTaskMCParticleFilter::~AliAnalysisTaskMCParticleFilter() { - // if( fTrackFilterMother ) delete fTrackFilterMother; + + if(fAODMcHeader){ + delete fAODMcHeader; + } + if(fAODMcParticles){ + fAODMcParticles->Delete(); + delete fAODMcParticles; + } } /* @@ -151,6 +170,8 @@ AliAnalysisTaskMCParticleFilter& AliAnalysisTaskMCParticleFilter::operator=(cons void AliAnalysisTaskMCParticleFilter::UserCreateOutputObjects() { // Create the output container + + if (OutputTree()&&fTrackFilterMother) OutputTree()->GetUserInfo()->Add(fTrackFilterMother); @@ -164,21 +185,19 @@ void AliAnalysisTaskMCParticleFilter::UserCreateOutputObjects() // mcparticles - TClonesArray *tca = new TClonesArray("AliAODMCParticle", 0); - tca->SetName(AliAODMCParticle::StdBranchName()); - AddAODBranch("TClonesArray",&tca); + fAODMcParticles = new TClonesArray("AliAODMCParticle", 0); + fAODMcParticles->SetName(AliAODMCParticle::StdBranchName()); + AddAODBranch("TClonesArray",&fAODMcParticles); // MC header... - AliAODMCHeader *mcHeader = new AliAODMCHeader(); - mcHeader->SetName(AliAODMCHeader::StdBranchName()); - AddAODBranch("AliAODMCHeader",&mcHeader); - - + fAODMcHeader = new AliAODMCHeader(); + fAODMcHeader->SetName(AliAODMCHeader::StdBranchName()); + AddAODBranch("AliAODMCHeader",&fAODMcHeader); AliMCEventHandler *mcH = (AliMCEventHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler()); AliAODHandler *aodH = dynamic_cast ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler()); if(!aodH){ - Printf("%s:&d Could not get AODHandler",(char*)__FILE__,__LINE__); + AliWarning("Could not get AODHandler"); return; } aodH->SetMCEventHandler(mcH); @@ -195,21 +214,23 @@ void AliAnalysisTaskMCParticleFilter::UserCreateOutputObjects() 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); + 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 - + // Execute analysis for current event + // + + // Fill AOD tracks from Kinematic stack + PostData(1,fHistList); + // get AliAOD Event AliAODEvent* aod = AODEvent(); if (!aod) { @@ -222,6 +243,13 @@ void AliAnalysisTaskMCParticleFilter::UserExec(Option_t */*option*/) AliWarning("No MC handler Found"); return; } + + AliAODHandler *aodH = dynamic_cast ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler()); + if(!aodH){ + AliWarning("Could not get AODHandler"); + return; + } + // fetch the output @@ -241,8 +269,58 @@ void AliAnalysisTaskMCParticleFilter::UserExec(Option_t */*option*/) 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(mcEH); + AliGenHijingEventHeader *hiH = 0; + AliCollisionGeometry *colG = 0; + AliGenDPMjetEventHeader *dpmH = 0; + // it can be only one save some casts + // assuming PYTHIA and HIJING are the most likely ones... + if(!pyH){ + hiH = dynamic_cast(mcEH); + if(!hiH){ + dpmH = dynamic_cast(mcEH); + } + } - // We take all real primaries + if (hiH || dpmH) colG = dynamic_cast(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(mcEH); + if (ccEH) { + TList *genHeaders = ccEH->GetHeaders(); + for (int imch=0; imchGetEntries(); imch++) { + if(!pyH)pyH = dynamic_cast(genHeaders->At(imch)); + if(!hiH)hiH = dynamic_cast(genHeaders->At(imch)); + if(!colG)colG = dynamic_cast(genHeaders->At(imch)); + if(!dpmH)dpmH = dynamic_cast(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())); + } + + Int_t j=0; @@ -259,18 +337,19 @@ void AliAnalysisTaskMCParticleFilter::UserExec(Option_t */*option*/) if (ip < nprim) { // Select the primary event write = kTRUE; - } else if (part->GetUniqueID() == 4) { + } 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() == 4)) { + 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 && Select(((AliMCParticle*) mcE->GetTrack(imo))->Particle(), rv, zv))write = kTRUE; - } else if (part->GetUniqueID() == 5) { + 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) { @@ -280,7 +359,7 @@ void AliAnalysisTaskMCParticleFilter::UserExec(Option_t */*option*/) // 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() == 4)) { + while((imo >= nprim) && (mother->GetUniqueID() == kPDecay)) { mother = (AliMCParticle*) mcE->GetTrack(imo); imo = mother->GetMother(); } @@ -289,55 +368,49 @@ void AliAnalysisTaskMCParticleFilter::UserExec(Option_t */*option*/) 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++; + j++; } } + /* - // check for charm daughters + // 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); + AliMCParticle* mcpart = (AliMCParticle*) mcE->GetTrack(ip); TParticle* part = mcpart->Particle(); - - // if((TMath::Abs(part->GetPdgCode())/400)==1){ + // debug info to check fro charm daugthers if((TMath::Abs(part->GetPdgCode()))==411){ - // cases iCharm++; - Int_t d0 = mcpart->GetFirstDaughter(); - Int_t d1 = mcpart->GetLastDaughter(); + Int_t d0 = part->GetFirstDaughter(); + Int_t d1 = part->GetLastDaughter(); if(d0>0&&d1>0){ for(int id = d0;id <= d1;id++){ - // AliMCParticle* daughter = mcE->GetTrack(id); iAll++; - if(mcH->IsParticleSelected(id))iTaken++; + 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)); } + */ - AliInfo(Form("Taken daughters %d/%d of %d charm",iTaken,iAll,iCharm)); - PostData(1,fHistList); + aodH->StoreMCParticles(); return; } @@ -348,18 +421,12 @@ void AliAnalysisTaskMCParticleFilter::Terminate(Option_t */*option*/){ // - 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; @@ -368,16 +435,6 @@ void AliAnalysisTaskMCParticleFilter::Terminate(Option_t */*option*/){ 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 *x = new TParameter("xsection",xsec); - tree->GetUserInfo()->Add(x); - TParameter *t = new TParameter("trials",trials); - tree->GetUserInfo()->Add(t); - */ - } @@ -413,3 +470,14 @@ Bool_t AliAnalysisTaskMCParticleFilter::Select(TParticle* part, Float_t rv, Floa } +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("---------------------------------------"); +}