STRANGENESS/Hypernuclei/AliAODRecoDecayLF2Prong.cxx
STRANGENESS/Hypernuclei/AliAODNuclExReplicator.cxx
STRANGENESS/Hypernuclei/AliAnalysisTaskESDNuclExFilter.cxx
+ STRANGENESS/Hypernuclei/AliAODMCNuclExReplicator.cxx
+ STRANGENESS/Hypernuclei/AliAnalysisTaskESDNuclExFilterMC.cxx
+ STRANGENESS/Hypernuclei/AliAnalysisTaskReadNuclexAOD.cxx
)
#pragma link C++ class AliAODRecoDecayLF2Prong+;
#pragma link C++ class AliAODNuclExReplicator+;
#pragma link C++ class AliAnalysisTaskESDNuclExFilter+;
-
+#pragma link C++ class AliAODMCNuclExReplicator+;
+#pragma link C++ class AliAnalysisTaskESDNuclExFilterMC+;
+#pragma link C++ class AliAnalysisTaskReadNuclexAOD+;
#endif
--- /dev/null
+AliAnalysisTask *AddTaskNuclexFilterMC( Bool_t onlyMuon = kTRUE,
+ Bool_t keepAllEvents = kTRUE,
+ Int_t mcMode = 0,
+ Int_t nsigmaTrk1 = 3,
+ Int_t partType1 = 2,
+ Int_t nsigmaTrk2 = 5,
+ Int_t partType2 = 7
+ )
+{
+
+ //get the current analysis manager
+
+ AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+ if (!mgr) {
+ Error("AddTaskNuclexFilter", "No analysis manager found.");
+ return 0;
+ }
+
+ //check for output aod handler
+ if (!mgr->GetOutputEventHandler()||mgr->GetOutputEventHandler()->IsA()!=AliAODHandler::Class()) {
+ Warning("AddTaskNuclExFilter","No AOD output handler available. Not adding the task!");
+ return 0;
+ }
+
+ //Do we have an MC handler?
+ // Bool_t hasMC=(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()!=0x0)||hasMC_aod;
+ /*
+ //Do we run on AOD?
+ Bool_t isAOD=mgr->GetInputEventHandler()->IsA()==AliAODInputHandler::Class();
+
+ //Allow merging of the filtered aods on grid trains
+ if(mgr->GetGridHandler()) {
+ printf(" SET MERGE FILTERED AODs \n");
+ mgr->GetGridHandler()->SetMergeAOD(kTRUE);
+ }
+
+ if(isAOD) {
+ //add options to AliAODHandler to duplicate input event
+ AliAODHandler *aodHandler = (AliAODHandler*)mgr->GetOutputEventHandler();
+ aodHandler->SetCreateNonStandardAOD();
+ }
+ */
+
+ Bool_t isAOD=mgr->GetInputEventHandler()->IsA()==AliAODInputHandler::Class();
+ if ( !isAOD) {
+ Printf("ERROR! This task can only run on AODs!");
+ }
+
+
+
+ //========= Add task to the ANALYSIS manager =====
+
+ cout<<"\n\n=============== addtask Nuclex =================== "<<endl<<endl;
+
+ //AliAnalysisTaskSE *esdnuclexfilter = new AliAnalysisTaskESDNuclExFilter("NuclEx Filter",onlyMuon,keepAllEvents,mcMode,nsigmaTrk1,partType1,nsigmaTrk2,partType2);
+
+ AliAnalysisTaskESDNuclExFilterMC *esdnuclexfilter = new AliAnalysisTaskESDNuclExFilterMC("NuclEx Filter MC",onlyMuon,keepAllEvents,mcMode,nsigmaTrk1,partType1,nsigmaTrk2,partType2);
+
+ esdnuclexfilter->SetWriteMuonAOD(kTRUE);
+
+ //================================================
+ // data containers
+ //================================================
+ // find input container
+
+ mgr->ConnectInput (esdnuclexfilter, 0, mgr->GetCommonInputContainer());
+ mgr->ConnectOutput (esdnuclexfilter, 0, mgr->GetCommonOutputContainer());
+
+
+ return esdnuclexfilter;
+}
--- /dev/null
+AliAnalysisTask *AddTaskReadNuclexAOD(TString name="name"){
+
+ //get the current analysis manager
+ AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+ if (!mgr) {
+ Error("AddTaskReadNuclexAOD", "No analysis manager found.");
+ return 0;
+ }
+
+ //========= Add task to the ANALYSIS manager =====
+
+ AliAnalysisTaskReadNuclexAOD *taskReadNuclexAOD = new AliAnalysisTaskReadNuclexAOD(name);
+
+ mgr->AddTask(taskReadNuclexAOD);
+
+ //================================================
+ // data containers
+ //================================================
+ // find input container
+
+ AliAnalysisDataContainer *cinput = mgr->GetCommonInputContainer();
+
+ TString outputFileName = AliAnalysisManager::GetCommonFileName();
+ //AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("Helium3Pi_tree", TTree::Class(), AliAnalysisManager::kOutputContainer, "AnalysisResults.root");
+ //AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("Helium3Pi_tree", TTree::Class(), AliAnalysisManager::kOutputContainer, "AnalysisResults.root");
+
+ AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("clisthistHyper" , TList::Class(), AliAnalysisManager::kOutputContainer, outputFileName);
+
+ AliAnalysisDataContainer *coutput2 = mgr->CreateContainer("treeAODrecoDecay", TTree::Class(),AliAnalysisManager::kOutputContainer,outputFileName);
+ AliAnalysisDataContainer *coutput3 = mgr->CreateContainer("treeMySecVert" , TTree::Class(),AliAnalysisManager::kOutputContainer,outputFileName);
+ // connect containers
+ mgr->ConnectInput (taskReadNuclexAOD, 0, cinput );
+ mgr->ConnectOutput (taskReadNuclexAOD, 1, coutput1);
+ mgr->ConnectOutput (taskReadNuclexAOD, 2, coutput2);
+ mgr->ConnectOutput (taskReadNuclexAOD, 3, coutput3);
+
+ return taskReadNuclexAOD;
+}
--- /dev/null
+/*************************************************************************
+* 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. *
+**************************************************************************/
+
+//
+// Implementation of a branch replicator
+// to produce aods with only few branches.
+//
+// This replicator is in charge of replicating the nuclei primary vertices
+// tracks identified as nuclei with Z>=2, secondary vertices in form of
+// AliAODRecoDecayLF2Prong and their daughter tracks.
+// These informations are stored into a reduced AODs (AliAOD.NuclEx.root)
+//
+// The vertices are filtered so that only the primary vertices make it
+// to the output aods.
+//
+// The secondary vertices are recreated here, as a AliAODRecoDecayLF2Prong
+// plus cuts that select secondary vericesoutside the primary vertex
+
+// Authors: S. Bufalino (stefania.bufalino@cern.ch)
+// R. Lea (ramona.lea@cern.ch)
+// Based on AliAODMuonReplicator.cxx
+
+//NOTE : This is a test on MC : no PID response only MC truth + select only 3LH
+// daughters : NOT INTENDED for any efficiency!!
+
+class AliESDv0;
+class AliESDVertex;
+class AliAODVertex;
+class AliAODRecoDecay;
+
+#include "AliStack.h"
+#include "AliAODEvent.h"
+#include "AliAODMCHeader.h"
+#include "AliAODMCParticle.h"
+#include "AliAODTZERO.h"
+#include "AliAODTrack.h"
+#include "AliAODVZERO.h"
+//#include "AliAnalysisCuts.h"
+#include "TF1.h"
+#include "AliExternalTrackParam.h"
+#include "AliESDv0.h"
+#include "AliAODv0.h"
+//#include "AliPIDResponse.h"
+#include <iostream>
+#include <cassert>
+#include "AliESDtrack.h"
+#include "TObjArray.h"
+#include "AliAnalysisFilter.h"
+#include "AliAODRecoDecay.h"
+#include "AliAODRecoDecayLF.h"
+#include "AliAODRecoDecayLF2Prong.h"
+
+#include <TFile.h>
+#include <TDatabasePDG.h>
+#include <TString.h>
+#include <TList.h>
+#include "AliLog.h"
+#include "AliVEvent.h"
+#include "AliVVertex.h"
+#include "AliVTrack.h"
+#include "AliVertexerTracks.h"
+#include "AliKFVertex.h"
+#include "AliESDEvent.h"
+#include "AliESDVertex.h"
+#include "AliESDtrackCuts.h"
+#include "AliAODEvent.h"
+#include "AliAnalysisFilter.h"
+//#include "AliAnalysisVertexingLF.h"
+#include "AliAnalysisManager.h"
+#include "AliAODMCNuclExReplicator.h"
+#include "TH1.h"
+#include "TCanvas.h"
+#include "AliInputEventHandler.h"
+
+using std::cout;
+using std::endl;
+
+ClassImp(AliAODMCNuclExReplicator)
+
+//_____________________________________________________________________________
+AliAODMCNuclExReplicator::AliAODMCNuclExReplicator(const char* name, const char* title,
+ Int_t mcMode,
+ Int_t nsigmaTrk1, Int_t partType1,
+ Int_t nsigmaTrk2, Int_t partType2
+ )
+ :AliAODBranchReplicator(name,title),
+ fBzkG(0.),
+ fCosMin(),
+ fDCAtracksMin(),
+ fRmax(),
+ fRmin(),
+ fDNmin(),
+ fDPmin(),
+ fHeader(0x0),
+ fVertices(0x0),
+ fNuclei(0x0),
+ fSecondaryVerices(0x0),
+ fDaughterTracks(0x0),
+ fList(0x0),
+ fMCParticles(0x0),
+ fMCHeader(0x0),
+ fMCMode(mcMode),
+ fLabelMap(),
+ fParticleSelected(),
+ fReplicateHeader(kTRUE), //replicateHeader //to be fixed
+ fnSigmaTrk1(nsigmaTrk1),
+ fnSigmaTrk2(nsigmaTrk2),
+ fpartType1(partType1),
+ fpartType2(partType2),
+ fSecVtxWithKF(kFALSE),
+ fVertexerTracks(0x0),
+ fV1(0x0),
+ fAODMapSize(0),
+ fAODMap(0)
+
+{
+ // default ctor
+}
+
+//_____________________________________________________________________________
+AliAODMCNuclExReplicator::~AliAODMCNuclExReplicator()
+{
+ // destructor
+ // delete fTrackCut;
+ // delete fVertexCut;
+ delete fList;
+}
+
+//_____________________________________________________________________________
+void AliAODMCNuclExReplicator::SelectParticle(Int_t i)
+{
+ // taking the absolute values here, need to take care
+ // of negative daughter and mother
+ // IDs when setting!
+
+ if (!IsParticleSelected(TMath::Abs(i)))
+ {
+ fParticleSelected.Add(TMath::Abs(i),1);
+ }
+}
+
+//_____________________________________________________________________________
+Bool_t AliAODMCNuclExReplicator::IsParticleSelected(Int_t i)
+{
+ // taking the absolute values here, need to take
+ // care with negative daughter and mother
+ // IDs when setting!
+ return (fParticleSelected.GetValue(TMath::Abs(i))==1);
+}
+
+
+//_____________________________________________________________________________
+void AliAODMCNuclExReplicator::CreateLabelMap(const AliAODEvent& source)
+{
+ //
+ // this should be called once all selections are done
+ //
+
+ fLabelMap.Delete();
+
+ TClonesArray* mcParticles = static_cast<TClonesArray*>(source.FindListObject(AliAODMCParticle::StdBranchName()));
+
+ Int_t i(0);
+ Int_t j(0);
+
+ TIter next(mcParticles);
+
+ while ( next() )
+ {
+ if (IsParticleSelected(i))
+ {
+ fLabelMap.Add(i,j++);
+ }
+ ++i;
+ }
+}
+
+//_____________________________________________________________________________
+Int_t AliAODMCNuclExReplicator::GetNewLabel(Int_t i)
+{
+ // Gets the label from the new created Map
+ // Call CreatLabelMap before
+ // otherwise only 0 returned
+ return fLabelMap.GetValue(TMath::Abs(i));
+}
+
+
+//_____________________________________________________________________________
+TList* AliAODMCNuclExReplicator::GetList() const
+{
+ // return (and build if not already done) our internal list of managed objects
+ if (!fList)
+ {
+
+ if ( fReplicateHeader )
+ {
+ fHeader = new AliAODHeader;
+ }
+
+
+ fSecondaryVerices = new TClonesArray("AliAODRecoDecayLF2Prong",30);
+ fSecondaryVerices->SetName("SecondaryVertices");
+
+ fVertices = new TClonesArray("AliAODVertex",2);
+ fVertices->SetName("vertices");
+
+ fNuclei = new TClonesArray("AliAODTrack",30);
+ fNuclei->SetName("Nuclei");
+
+ fDaughterTracks = new TClonesArray("AliAODTrack",30);
+ fDaughterTracks->SetName("DaughterTracks");
+
+
+ fList = new TList;
+ fList->SetOwner(kTRUE);
+
+ fList->Add(fHeader);
+ fList->Add(fVertices);
+ fList->Add(fNuclei);
+ fList->Add(fSecondaryVerices);
+ fList->Add(fDaughterTracks);
+
+
+ if ( fMCMode > 0 )
+ {
+ fMCHeader = new AliAODMCHeader;
+ fMCParticles = new TClonesArray("AliAODMCParticle",1000);
+ fMCParticles->SetName(AliAODMCParticle::StdBranchName());
+ fList->Add(fMCHeader);
+ fList->Add(fMCParticles);
+ }
+ }
+ return fList;
+}
+
+//_____________________________________________________________________________
+void AliAODMCNuclExReplicator::ReplicateAndFilter(const AliAODEvent& source)
+{
+ // Replicate (and filter if filters are there) the relevant parts we're interested in AODEvent
+
+ // cout<<"-------------------->QUESTO"<<endl;
+
+ //-----------------------------------------------
+ // AliPIDResponse
+
+ // AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
+ // AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
+
+ //--------------------------------------------------------
+
+ // printf("Execute NuclEx Replicator\n");
+
+ //---------------------------------
+
+ if (fReplicateHeader)
+ {
+ *fHeader = *(source.GetHeader());
+ }
+
+ fVertices->Clear("C");
+
+ fNuclei->Clear("C");
+
+ fSecondaryVerices->Clear("C");
+
+ fDaughterTracks->Clear("C");
+
+ //----------------------------------
+
+ //retrive MC infos
+
+ TClonesArray *arrayMC = 0;
+ AliAODMCHeader *mcHeader=0;
+ Int_t mumpdg=-100;
+
+ arrayMC = (TClonesArray*) source.GetList()->FindObject(AliAODMCParticle::StdBranchName());
+ if (!arrayMC) {
+ Printf("Error: MC particles branch not found!\n");
+ return;
+ }
+ // if(arrayMC)
+ // cout<<"Ho caricato array mc"<<endl;
+
+ mcHeader = (AliAODMCHeader*)source.GetList()->FindObject(AliAODMCHeader::StdBranchName());
+ if(!mcHeader) {
+ printf("AliAnalysisTaskSEDs::UserExec: MC header branch not found!\n");
+ return;
+ }
+ // if(mcHeader)
+ // cout<<"Ho caricato MC header"<<endl;
+
+ // cout<<"Centrality AOD source: "<<source.GetHeader()->GetCentrality()<<endl;
+
+ Int_t nsv(0);
+ // Int_t nnuclei(0);
+ Int_t ntracksD(0);
+
+ Int_t input(0);
+ Double_t xdummy,ydummy;
+
+ AliAODRecoDecayLF2Prong *io2Prong = 0;
+
+ TObjArray *twoTrackArray = new TObjArray(2);
+ Double_t dispersion;
+
+ // cout<<"Qui"<<endl;
+ //cout<<source.GetMagneticField()<<endl;
+
+ AliAODVertex *vtx = source.GetPrimaryVertex();
+
+ // cout<<"Source "<<source<<endl;
+ //cout<<"vtx: "<<vtx<<endl;
+
+ // A Set of very loose cut for a weak strange decay
+
+ fCosMin = 0.99;
+ fDCAtracksMin = 1;
+ fRmax = 200.;
+ fRmin = 0.1;
+ fDNmin = 0.05;
+ fDPmin = 0.05;
+
+ //----------------------------------------------------------
+
+ // Int_t nindices=0;
+ UShort_t *indices = 0;
+ const Int_t entries = source.GetNumberOfTracks();
+
+ Double_t pos[3],cov[6];
+ vtx->GetXYZ(pos);
+ vtx->GetCovarianceMatrix(cov);
+ fV1 = new AliESDVertex(pos,cov,100.,100,vtx->GetName());
+ // cout<<"fV1 pricipal loop: "<<fV1<<endl;
+
+ if(entries<=0) return;
+ indices = new UShort_t[entries];
+ memset(indices,0,sizeof(UShort_t)*entries);
+ fAODMapSize = 100000;
+ fAODMap = new Int_t[fAODMapSize];
+ memset(fAODMap,0,sizeof(Int_t)*fAODMapSize);
+ // cent=((AliAODEvent&)source)->GetCentrality();
+
+ //-------------------------------------------------------------
+
+ if(vtx->GetNContributors()<1) {
+
+ // SPD vertex cut
+ vtx =source.GetPrimaryVertexSPD();
+
+ if(vtx->GetNContributors()<1) {
+ Info("AliAnalysisTaskHelium3Pi","No good vertex, skip event");
+ return; // NO GOOD VERTEX, SKIP EVENT
+ }
+ }
+
+ Double_t xPrimaryVertex=0.,yPrimaryVertex=0.;
+ xPrimaryVertex=vtx->GetX();
+ yPrimaryVertex=vtx->GetY();
+
+ fBzkG=source.GetMagneticField();
+ fVertexerTracks=new AliVertexerTracks(fBzkG);
+
+ Double_t TrackNumber = source.GetNumberOfTracks();
+ Int_t label =-1;
+
+ //Tracks arrays
+
+ TArrayI Track0(TrackNumber); //Pions
+ Int_t nTrack0=0;
+
+ TArrayI Track1(TrackNumber); //Helium3
+ Int_t nTrack1=0;
+
+ for(Int_t j=0; j<TrackNumber; j++){
+
+ // cout<<"Inside loop tracks"<<endl;
+
+
+ AliVTrack *track = (AliVTrack*)source.GetTrack(j);
+
+ AliAODTrack *aodtrack =(AliAODTrack*)track;
+
+ //-----------------------------------------------------------
+ //Track cuts
+ if(aodtrack->GetTPCNcls() < 70 )continue;
+ if(aodtrack->Chi2perNDF() > 4 )continue;
+
+ if (!aodtrack->IsOn(AliAODTrack::kTPCrefit)) continue;
+ if (!aodtrack->IsOn(AliAODTrack::kTPCin)) continue;
+ if (aodtrack->IsOn(AliAODTrack::kITSpureSA)) continue;
+
+ //---------------------------------------------------------------
+
+ Double_t mom = aodtrack->P();
+
+ if(mom<0.150)continue;
+
+ label = TMath::Abs(aodtrack->GetLabel());
+
+ AliAODMCParticle *part = (AliAODMCParticle*) arrayMC->At(label);
+
+ Int_t PDGCode=part->GetPdgCode();
+
+ Int_t mumid = part->GetMother();
+
+ if(mumid>-1){
+ AliAODMCParticle *mother = (AliAODMCParticle*) arrayMC->At(mumid);
+ mumpdg = mother->GetPdgCode();
+ }
+
+ // if(mumpdg == 1010010030 ||mumpdg == -1010010030 ){
+
+ if(mumpdg == 1010010030){
+
+ // if(PDGCode==-211 || PDGCode==+211){
+ if(PDGCode==-211){
+ Track0[nTrack0++]=j;
+ }
+
+ // if(PDGCode==1000020030 ||PDGCode==-1000020030 ){
+ if(PDGCode==1000020030){
+ Track1[nTrack1++]=j;
+ // new((*fNuclei)[nnuclei++]) AliAODTrack(*aodtrack);
+ }
+ }
+ }
+
+ //Set Track Daughters
+ Track0.Set(nTrack0);
+ Track1.Set(nTrack1);
+
+
+ // cout<<"Track loops..."<<endl;
+ // cout<<"npos "<<nTrack1<<endl;
+ // cout<<"nneg "<<nTrack0<<endl;
+
+
+ AliAODTrack *postrackAOD = 0;
+ AliAODTrack *negtrackAOD = 0;
+
+ AliESDtrack *postrack = 0;
+ AliESDtrack *negtrack = 0;
+
+ Bool_t isOk=kFALSE;
+
+ for (Int_t i=0; i<nTrack1; i++){ //! He Tracks Loop
+
+ Int_t Track1idx=Track1[i];
+
+ AliVTrack *trackpos = (AliVTrack*)source.GetTrack(Track1idx);
+
+ postrackAOD = (AliAODTrack*)trackpos;
+ postrack = new AliESDtrack(trackpos);
+
+ //--- MC infos
+
+ Int_t labelpos = TMath::Abs(postrack->GetLabel());
+ AliAODMCParticle *partPos = (AliAODMCParticle*) arrayMC->At(labelpos);
+ Int_t mumidPos = partPos->GetMother();
+
+ //------------------------------
+
+ for (Int_t k=0; k <nTrack0 ; k++) { //! Pion Tracks Loop
+
+ Int_t Track0idx=Track0[k];
+
+ AliVTrack *trackneg = (AliVTrack*)source.GetTrack(Track0idx);
+ negtrackAOD =(AliAODTrack*)trackneg;
+ negtrack = new AliESDtrack(trackneg);
+
+ //--- MC infos
+
+ Int_t labelneg = TMath::Abs(negtrack->GetLabel());
+ AliAODMCParticle *partNeg = (AliAODMCParticle*) arrayMC->At(labelneg);
+ Int_t mumidNeg = partNeg->GetMother();
+
+
+ //------------------------------
+ // if(mumidPos == mumidNeg && mumidNeg > 0){
+ isOk=kFALSE;
+
+ if(mumidPos == mumidNeg && mumidNeg > 0)
+ isOk = kTRUE;
+
+ twoTrackArray->AddAt(negtrack,0);
+ twoTrackArray->AddAt(postrack,1);
+
+ Double_t dcap1n1 = postrack->GetDCA(negtrack,fBzkG,xdummy,ydummy);
+
+ Double_t dcap1toPV = TMath::Abs(postrack->GetD(xPrimaryVertex, yPrimaryVertex,fBzkG));
+ Double_t dcan1toPV = TMath::Abs(negtrack->GetD(xPrimaryVertex, yPrimaryVertex,fBzkG));
+
+ if(dcap1n1>fDCAtracksMin)continue;
+ if((xdummy+ydummy)>fRmax )continue;
+ if((xdummy+ydummy)< fRmin)continue;
+
+ if ( dcan1toPV < fDNmin)
+ if ( dcap1toPV < fDPmin) continue;
+
+ // cout<<"dcap1n1: "<<dcap1n1<<endl;
+
+ AliAODVertex *vertexp1n1 = ReconstructSecondaryVertex(twoTrackArray,dispersion,kTRUE);
+
+ if(!vertexp1n1) {
+
+ twoTrackArray->Clear();
+ delete negtrack;
+ negtrack=NULL;
+ continue;
+ }
+
+ io2Prong = Make2Prong(twoTrackArray,source,vertexp1n1,dcap1n1);
+
+ if(io2Prong->CosPointingAngle()<fCosMin)continue;
+
+
+ // AliAODTrack *trk0 = (AliAODTrack*)io2Prong->GetDaughter(0);
+ // AliAODTrack *trk1 = (AliAODTrack*)io2Prong->GetDaughter(1);
+
+ // cout<<"**********************************************"<<endl;
+ // cout<<trk0/*->GetID()*/<<" "<<negtrackAOD->GetID()<<endl;
+ // cout<<trk1/*->GetID()*/<<" "<<postrackAOD->GetID()<<endl;
+ // cout<<"d0 io2Prong: "<<io2Prong->GetProngID(1)<<endl;
+ // cout<<"d1 io2Prong: "<<io2Prong->GetProngID(0)<<endl;
+ // cout<<"**********************************************"<<endl;
+
+ // rd = new((*fSecondaryVerices)[nsv++]) AliAODRecoDecayLF2Prong(*io2Prong);
+ if(isOk){
+ new((*fSecondaryVerices)[nsv++]) AliAODRecoDecayLF2Prong(*io2Prong);
+ new((*fDaughterTracks)[ntracksD++]) AliAODTrack(*negtrackAOD);
+ new((*fDaughterTracks)[ntracksD++]) AliAODTrack(*postrackAOD);
+
+ }
+ // rd->SetSecondaryVtx(vertexp1n1);
+ // vertexp1n1->SetParent(rd);
+ // AddRefs(vertexp1n1,rd,source,twoTrackArray);
+
+ delete negtrack;
+ negtrack = NULL;
+
+ delete vertexp1n1;
+ vertexp1n1= NULL;
+ continue;
+ // }
+ }
+
+ delete postrack;
+ postrack = NULL;
+
+ }
+
+ //----------------------------------------------------------
+
+ assert(fVertices!=0x0);
+ fVertices->Clear("C");
+ TIter nextV(source.GetVertices());
+ AliAODVertex* v;
+ Int_t nvertices(0);
+
+ while ( ( v = static_cast<AliAODVertex*>(nextV()) ) )
+ {
+ if ( v->GetType() == AliAODVertex::kPrimary ||
+ v->GetType() == AliAODVertex::kMainSPD ||
+ v->GetType() == AliAODVertex::kPileupSPD ||
+ v->GetType() == AliAODVertex::kPileupTracks||
+ v->GetType() == AliAODVertex::kMainTPC )
+ {
+
+ AliAODVertex* tmp = v->CloneWithoutRefs();
+ //AliAODVertex* copiedVertex = new((*fVertices)[nvertices++]) AliAODVertex(*tmp);
+ new((*fVertices)[nvertices++]) AliAODVertex(*tmp);
+
+ // to insure the main vertex retains the ncontributors information
+ // (which is otherwise computed dynamically from
+ // references to tracks, which we do not keep in muon aods...)
+ // we set it here
+
+ //copiedVertex->SetNContributors(v->GetNContributors());
+
+ // fVertices->Delete();
+ // delete copiedVertex;
+ delete tmp;
+ // printf("....Prendo il vertice primario...\n");
+ }
+ // printf("....Prendo il vertice primario...\n");
+ }
+
+ //printf("....Done NuclEx Replicator...\n");
+
+ AliDebug(1,Form("input mu tracks=%d tracks=%d vertices=%d nnuclei=%d",
+ input,fSecondaryVerices->GetEntries(),fVertices->GetEntries(),fNuclei->GetEntries()));
+
+ // cout<<"Delete..."<<endl;
+ // delete foPion;
+ // foPion = NULL;
+ //cout<<"Delete 1"<< endl;
+
+ if(io2Prong) {delete io2Prong; io2Prong=NULL;}
+ //cout<<"Delete 2"<< endl;
+ twoTrackArray->Delete(); delete twoTrackArray;
+ // cout<<"Delete 3"<< endl;
+ // vtx->Delete(); delete vtx;
+ // cout<<"Delete 4"<< endl;
+ if(fV1) { delete fV1; fV1=NULL; }
+ // cout<<"Delete 5"<< endl;
+ if(fAODMap) { delete [] fAODMap; fAODMap=NULL; }
+ delete []indices;
+ // cout<<"Delete 6"<< endl;
+ delete fVertexerTracks;
+
+ // cout<<"Mi rompo alla fine. Perche???"<<endl;
+
+}
+
+//-----------------------------------------------------------------------------
+
+AliAODVertex *AliAODMCNuclExReplicator::ReconstructSecondaryVertex(TObjArray *trkArray,
+ Double_t &dispersion,Bool_t useTRefArray) const
+
+{
+ // Secondary vertex reconstruction with AliVertexerTracks or AliKFParticle
+ //AliCodeTimerAuto("",0);
+
+ AliESDVertex *vertexESD = 0;
+ AliAODVertex *vertexAOD = 0;
+
+ if(!fSecVtxWithKF) { // AliVertexerTracks
+
+ fVertexerTracks->SetVtxStart(fV1);
+ vertexESD = (AliESDVertex*)fVertexerTracks->VertexForSelectedESDTracks(trkArray);
+
+ if(!vertexESD) return vertexAOD;
+
+
+ Double_t vertRadius2=vertexESD->GetXv()*vertexESD->GetXv()+vertexESD->GetYv()*vertexESD->GetYv();
+ if(vertRadius2>200.){
+ // vertex outside beam pipe, reject candidate to avoid propagation through material
+ delete vertexESD; vertexESD=NULL;
+ return vertexAOD;
+ }
+
+ } else { // Kalman Filter vertexer (AliKFParticle)
+
+ AliKFParticle::SetField(fBzkG);
+
+ AliKFVertex vertexKF;
+
+ Int_t nTrks = trkArray->GetEntriesFast();
+ for(Int_t i=0; i<nTrks; i++) {
+ AliESDtrack *esdTrack = (AliESDtrack*)trkArray->At(i);
+ AliKFParticle daughterKF(*esdTrack,211);
+ vertexKF.AddDaughter(daughterKF);
+ }
+ vertexESD = new AliESDVertex(vertexKF.Parameters(),
+ vertexKF.CovarianceMatrix(),
+ vertexKF.GetChi2(),
+ vertexKF.GetNContributors());
+
+ }
+ // convert to AliAODVertex
+ Double_t pos[3],cov[6],chi2perNDF;
+ vertexESD->GetXYZ(pos); // position
+ vertexESD->GetCovMatrix(cov); //covariance matrix
+ chi2perNDF = vertexESD->GetChi2toNDF();
+ dispersion = vertexESD->GetDispersion();
+ delete vertexESD;
+ vertexESD=NULL;
+
+ Int_t nprongs= (useTRefArray ? 0 : trkArray->GetEntriesFast());
+ vertexAOD = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,nprongs);
+
+ // cout<<"------------------> Reconstruct vertexAOD: "<<vertexAOD<<endl;
+
+ return vertexAOD;
+
+
+}
+
+//-----------------------------------------------------------------------------
+
+AliAODRecoDecayLF2Prong* AliAODMCNuclExReplicator::Make2Prong(TObjArray *twoTrackArray,const AliAODEvent &evento,
+ AliAODVertex *secVert,Double_t dca)
+
+
+{
+
+ //cout<<"Inside make 2 prong"<<endl;
+
+ Int_t charge[2];
+ charge[0]=1; //it was -1 //Ramona
+ charge[1]=2;
+
+ // From AliAnalysisVertexingLF.cxx Method:Make2Prongs
+
+ fBzkG = evento.GetMagneticField();
+
+ Double_t px[2],py[2],pz[2],d0[2],d0err[2];
+ // Also this has been changed //Ramona
+ AliESDtrack *negtrack = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
+ AliESDtrack *postrack = (AliESDtrack*)twoTrackArray->UncheckedAt(1);
+
+ //cout<<"negtrack: "<<negtrack<<" postrack: "<<postrack<<endl;
+ //cout<<"kVeryBig: "<<kVeryBig<<endl;
+ //cout<<"secVert: "<<secVert<<endl;
+
+ // // propagate tracks to secondary vertex, to compute inv. mass
+
+ postrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
+ negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
+
+ Double_t momentum[3];
+
+ negtrack->GetPxPyPz(momentum);
+ px[0] = charge[0]*momentum[0];
+ py[0] = charge[0]*momentum[1];
+ pz[0] = charge[0]*momentum[2];
+
+ postrack->GetPxPyPz(momentum);
+ px[1] = charge[1]*momentum[0];
+ py[1] = charge[1]*momentum[1];
+ pz[1] = charge[1]*momentum[2];
+
+ //cout<< px[0] <<" "<< " "<< py[0] << " "<< pz[0]<<endl;
+ // px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
+ //cout<< px[1] <<" "<< " "<< py[1] << " "<< pz[1]<<endl;
+ //px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
+
+
+ // primary vertex to be used by this candidate
+ AliAODVertex *primVertexAOD = evento.GetPrimaryVertex();
+ //cout<<"primVertexAOD "<<primVertexAOD<<endl;
+ if(!primVertexAOD) return 0x0;
+
+ Double_t d0z0[2],covd0z0[3];
+
+ d0z0[0] = -999.;
+ d0z0[1] = -999.;
+ covd0z0[0] = -999.;
+ covd0z0[1] = -999.;
+ covd0z0[2] = -999.;
+
+ d0[0] = d0z0[0];
+ d0err[0] = TMath::Sqrt(TMath::Abs(covd0z0[0]));
+
+ d0[1] = d0z0[0];
+ d0err[1] = TMath::Sqrt(TMath::Abs(covd0z0[0]));
+
+ // create the object AliAODRecoDecayLF2Prong
+ // AliAODRecoDecayLF2Prong *the2Prong = new AliAODRecoDecayLF2Prong(secVert,px,py,pz,d0,d0err,dcap1n1);
+ AliAODRecoDecayLF2Prong *the2Prong = new AliAODRecoDecayLF2Prong(secVert,px,py,pz,d0,d0err,dca);
+ the2Prong->SetOwnPrimaryVtx(primVertexAOD);
+
+ // the2Prong->SetProngIDs(2,{-999,999});
+ // UShort_t id0[2]={99999,999999};
+ // the2Prong->SetProngIDs(2,id0);
+
+ UShort_t id[2]={(UShort_t)postrack->GetID(),(UShort_t)negtrack->GetID()};
+ the2Prong->SetProngIDs(2,id);
+
+ // cout<<"\n\n\nMake 2 Prong: id[0]"<<id[0]<<" id[1]: "<<id[1]<<endl;
+ // cout<<"Get: 1 "<<the2Prong->GetProngID(0)<<" 2 "<<the2Prong->GetProngID(1)<<endl;
+ // cout<<"Get: 3 "<<the2Prong->GetProngID(2)<<" 4 "<<the2Prong->GetProngID(3)<<endl<<endl<<endl;
+ //delete primVertexAOD; primVertexAOD=NULL;
+
+ if(postrack->Charge()!=0 && negtrack->Charge()!=0) {
+
+ AddDaughterRefs(secVert,(AliAODEvent&)evento,twoTrackArray);
+ // AddDaughterRefs(secVert,(AliAODEvent*)evento,twoTrackArray);
+
+ }
+
+ return the2Prong;
+
+ delete the2Prong;
+}
+
+//----------------------------------------------------------------------------
+void AliAODMCNuclExReplicator::AddDaughterRefs(AliAODVertex *v,
+ const AliAODEvent &event,
+ const TObjArray *trkArray) const
+
+{
+ // Add the AOD tracks as daughters of the vertex (TRef)
+ // AliCodeTimerAuto("",0);
+ // cout<<"Inside AddDaughterRefs "<<endl;
+
+ Int_t nDg = v->GetNDaughters();
+
+ TObject *dg = 0;
+ if(nDg) dg = v->GetDaughter(0);
+ if(dg) return; // daughters already added
+
+ Int_t nTrks = trkArray->GetEntriesFast();
+
+ AliExternalTrackParam *track = 0;
+ AliAODTrack *aodTrack = 0;
+ Int_t id;
+
+ for(Int_t i=0; i<nTrks; i++) {
+ track = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
+
+ id = (Int_t)track->GetID();
+ // printf("---> %d\n",id);
+ if(id<0) continue; // this track is a AliAODRecoDecay
+
+ aodTrack = (AliAODTrack*)event.GetTrack(fAODMap[id]);
+ v->AddDaughter(aodTrack);
+ }
+ return;
+}
+//----------------------------------------------------------------------------
+
+void AliAODMCNuclExReplicator::AddRefs(AliAODVertex *v,AliAODRecoDecayLF *rd,
+ const AliAODEvent &event,
+ const TObjArray *trkArray) const
+
+{
+ // Add the AOD tracks as daughters of the vertex (TRef)
+ // Also add the references to the primary vertex and to the cuts
+ //AliCodeTimerAuto("",0);
+
+ AddDaughterRefs(v,event,trkArray);
+ rd->SetPrimaryVtxRef((AliAODVertex*)event.GetPrimaryVertex());
+ return;
+}
+
+//---------------------------------------------------------------------------
+
+void AliAODMCNuclExReplicator::Terminate(){
+
+}
--- /dev/null
+#ifndef ALIAODMCNUCLEXREPLICATOR_H
+#define ALIAODMCNUCLEXREPLICATOR_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+* See cxx source for full Copyright notice */
+
+#ifndef ALIDAODBRANCHREPLICATOR_H
+# include "AliAODBranchReplicator.h"
+#endif
+#ifndef ROOT_TExMap
+# include "TExMap.h"
+#endif
+
+
+//
+// Implementation of a branch replicator
+// to produce aods with only few branches.
+//
+// Authors: S. Bufalino (stefania.bufalino@cern.ch)
+// R. Lea (ramona.lea@cern.ch)
+//
+// Based on AliAODMuonReplicator.h (L. Aphecetche (Subatech))
+
+//class AliAnalysisCuts;
+class TClonesArray;
+class AliAODMCHeader;
+class AliAODVZERO;
+class AliAODTZERO;
+//class AliPIDResponse;
+class AliESDv0;
+class TArrayI;
+class AliAODv0;
+class TRefArray;
+class AliAODRecoDecay;
+class AliAODRecoDecayLF;
+class AliAODRecoDecayLF2Prong;
+class AliVertexerTracks;
+class AliAODHeader;
+
+class AliESDVertex;
+class AliESDtrack;
+class AliVEvent;
+class AliAODVertex;
+class AliVertexerTracks;
+class AliESDv0;
+class AliAODv0;
+
+class TH1F;
+
+// TODO add new constructor for the 3 prong case (it will write an AliAODRecoDecayLF3Prong
+
+class AliAODMCNuclExReplicator : public AliAODBranchReplicator
+{
+ public:
+
+ AliAODMCNuclExReplicator(const char* name="AliAODMCNuclExReplicator",
+ const char* title="Branch Replicator for muon related branches",
+ /* AliAnalysisCuts* trackCut=0x0, */
+ /* AliAnalysisCuts* vertexCut=0x0, */
+ Int_t mcMode=0,
+ Int_t nsigmaTrk1=3, Int_t partType1 = 2,
+ Int_t nsigmaTrk2=3, Int_t partType2 = 7
+ );
+ // Int_t partType; // 0 = e; 1 = mu; 2 = pi; 3 = K; 4= p; 5 = d; 6 =t ; 7 = 3He; 8=4He;
+
+ virtual ~AliAODMCNuclExReplicator();
+
+ virtual TList* GetList() const;
+
+ virtual void ReplicateAndFilter(const AliAODEvent& source);
+
+ virtual void Terminate();
+
+ private:
+
+ // TO DO : Implemet MC
+ // void FilterMC(const AliAODEvent& source);
+ void SelectParticle(Int_t i);
+ Bool_t IsParticleSelected(Int_t i);
+ void CreateLabelMap(const AliAODEvent& source);
+ Int_t GetNewLabel(Int_t i);
+
+
+ Double_t fBzkG; // z componenent of field in kG
+
+ Double_t fCosMin ;
+ Double_t fDCAtracksMin ;
+ Double_t fRmax ;
+ Double_t fRmin ;
+ Double_t fDNmin ;
+ Double_t fDPmin ;
+
+ private:
+
+ mutable AliAODHeader* fHeader; //! internal header object
+
+ // AliAnalysisCuts* fVertexCut; // decides which vertices to keep
+ mutable TClonesArray* fVertices; //! internal array of vertices
+
+ mutable TClonesArray* fNuclei; //! internal array of nuclei tracks
+
+ // AliAnalysisCuts* fTrackCut; // decides which tracks to keep
+ mutable TClonesArray* fSecondaryVerices; //! internal array of secondary vertices canditates
+
+ mutable TClonesArray* fDaughterTracks; //! internal SV daughter tracks
+
+
+
+ mutable TList* fList; //! internal list of managed objects (fVertices and fTracks)
+
+ mutable TClonesArray* fMCParticles; //! internal array of MC particles
+ mutable AliAODMCHeader* fMCHeader; //! internal array of MC header
+ Int_t fMCMode; // MC filtering switch (0=none=no mc information,1=normal=simple copy,>=2=aggressive=filter out)
+
+ TExMap fLabelMap; //! for MC label remapping (in case of aggressive filtering)
+ TExMap fParticleSelected; //! List of selected MC particles
+
+ Bool_t fReplicateHeader; //! whether or not the replicate the AOD Header
+
+ Int_t fnSigmaTrk1;
+ Int_t fnSigmaTrk2;
+ Int_t fpartType1;
+ Int_t fpartType2;
+
+
+
+ Bool_t fSecVtxWithKF; // if kTRUE use KF vertexer, else AliVertexerTracks
+ AliVertexerTracks* fVertexerTracks; // vertexer, to compute secondary vertices
+ AliESDVertex *fV1; // primary vertex
+
+ Int_t fAODMapSize; // size of fAODMap
+ Int_t *fAODMap; //[fAODMapSize] map between index and ID for AOD tracks
+
+ AliAODVertex* ReconstructSecondaryVertex(TObjArray *trkArray,Double_t &dispersion,Bool_t useTRefArray=kTRUE) const;
+
+ AliAODRecoDecayLF2Prong* Make2Prong(TObjArray *twoTrackArray,const AliAODEvent &evento,
+ AliAODVertex *secVert,Double_t dca);
+
+ /* AliAODRecoDecayLF2Prong* Make2Prong(TObjArray *twoTrackArray,AliAODEvent *evento, */
+ /* AliAODVertex *secVert,Double_t dca); */
+
+
+ void AddDaughterRefs(AliAODVertex *v, const AliAODEvent &event, const TObjArray *trkArray) const;
+ /* void AddDaughterRefs(AliAODVertex *v, AliAODEvent *event, const TObjArray *trkArray) const; */
+ void AddRefs(AliAODVertex *v,AliAODRecoDecayLF *rd, const AliAODEvent &event, const TObjArray *trkArray) const;
+ /* void AddRefs(AliAODVertex *v,AliAODRecoDecayLF *rd, AliAODEvent *event, const TObjArray *trkArray) const; */
+
+ private:
+
+ // AliPIDResponse *fPIDResponse; //! PID response object
+
+ AliAODMCNuclExReplicator(const AliAODMCNuclExReplicator&);
+ AliAODMCNuclExReplicator& operator=(const AliAODMCNuclExReplicator&);
+
+ ClassDef(AliAODMCNuclExReplicator,5) // Branch replicator for ESD to muon AOD.
+ };
+
+#endif
--- /dev/null
+/**************************************************************************
+ * 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. *
+ **************************************************************************/
+
+//
+// Create a new AOD starting from the general AOD. This Task can be used also strating
+//from ESD changing the input handler. (Method to be testeted on the grid)
+// filtering of the ESD.
+//
+// Authors: S. Bufalino (stefania.bufalino@cern.ch)
+// R. Lea (ramona.lea@cern.ch)
+// Based on AliAnalysisTaskESDMuonFilter.cxx
+//
+// (see AddFilteredAOD method)
+//
+
+#include "AliAnalysisTaskESDNuclExFilterMC.h"
+
+#include "AliAODEvent.h"
+#include "AliAODHandler.h"
+#include "AliAODExtension.h"
+#include "AliAODMCParticle.h"
+#include "AliAODMCNuclExReplicator.h"
+#include "AliAODVertex.h"
+#include "AliAnalysisFilter.h"
+#include "AliAnalysisManager.h"
+#include "AliCodeTimer.h"
+#include "AliESDEvent.h"
+#include "AliESDInputHandler.h"
+#include "AliESDMuonTrack.h"
+#include "AliESDVertex.h"
+#include "AliESDtrack.h"
+#include "AliLog.h"
+#include "AliMCEvent.h"
+#include "AliMCEventHandler.h"
+#include "AliMultiplicity.h"
+#include <TChain.h>
+#include <TFile.h>
+#include <TParticle.h>
+#include "AliESDtrackCuts.h"
+#include "AliAODTrack.h"
+#include "AliAODPid.h"
+#include "TF1.h"
+//#include "AliPIDResponse.h"
+
+using std::cout;
+using std::endl;
+ClassImp(AliAnalysisTaskESDNuclExFilterMC)
+
+////////////////////////////////////////////////////////////////////////
+
+AliAnalysisTaskESDNuclExFilterMC::AliAnalysisTaskESDNuclExFilterMC(Bool_t onlyMuon, Bool_t keepAllEvents, Int_t mcMode, Int_t nsigmaTrk1,Int_t nsigmaTrk2, Int_t partType1,Int_t partType2):
+ AliAnalysisTaskSE(),
+ fTrackFilter(0x0),
+ fEnableMuonAOD(kTRUE),
+ fEnableDimuonAOD(kTRUE),
+ fOnlyMuon(onlyMuon),
+ fKeepAllEvents(keepAllEvents),
+ fMCMode(mcMode),
+ fnSigmaTrk1(nsigmaTrk1),
+ fnSigmaTrk2(nsigmaTrk2),
+ fpartType1(partType1),
+ fpartType2(partType2),
+ murep(0x0)
+{
+ // Default constructor
+}
+
+AliAnalysisTaskESDNuclExFilterMC::AliAnalysisTaskESDNuclExFilterMC(const char* name, Bool_t onlyMuon, Bool_t keepAllEvents, Int_t mcMode, Int_t nsigmaTrk1,Int_t nsigmaTrk2, Int_t partType1,Int_t partType2):
+ AliAnalysisTaskSE(name),
+ fTrackFilter(0x0),
+ fEnableMuonAOD(kTRUE),
+ fEnableDimuonAOD(kTRUE),
+ fOnlyMuon(onlyMuon),
+ fKeepAllEvents(keepAllEvents),
+ fMCMode(mcMode),
+ fnSigmaTrk1(nsigmaTrk1),
+ fnSigmaTrk2(nsigmaTrk2),
+ fpartType1(partType1),
+ fpartType2(partType2),
+ murep(0x0)
+
+{
+ // Constructor
+}
+
+//______________________________________________________________________________
+void AliAnalysisTaskESDNuclExFilterMC::UserCreateOutputObjects()
+{
+
+ // Create the output container
+ if (fTrackFilter) OutputTree()->GetUserInfo()->Add(fTrackFilter);
+
+}
+
+//______________________________________________________________________________
+void AliAnalysisTaskESDNuclExFilterMC::PrintTask(Option_t *option, Int_t indent) const
+{
+ // Specify how we are configured
+
+ AliAnalysisTaskSE::PrintTask(option,indent);
+
+ TString spaces(' ',indent+3);
+
+ if ( fOnlyMuon )
+ {
+ cout << spaces.Data() << "Keep only muon information " << endl;
+ }
+ else
+ {
+ cout << spaces.Data() << "Keep all information from standard AOD" << endl;
+ }
+
+ if ( fKeepAllEvents )
+ {
+ cout << spaces.Data() << "Keep all events, regardless of number of muons" << endl;
+ }
+ else
+ {
+ cout << spaces.Data() << "Keep only events with at least one muon" << endl;
+ }
+
+ if ( fMCMode > 0 )
+ {
+ cout << spaces.Data() << "Assuming work on MC data (i.e. will transmit MC branches)" << endl;
+ }
+}
+
+//______________________________________________________________________________
+void AliAnalysisTaskESDNuclExFilterMC::AddFilteredAOD(const char* aodfilename, const char* title, Bool_t toMerge)
+{
+
+ //cout<<"Entro ne ADDFILTETEDAOD"<<endl;
+
+ AliAODHandler *aodH = (AliAODHandler*)((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
+ if (!aodH) Fatal("UserCreateOutputObjects", "No AOD handler");
+
+ if(aodH){
+
+ AliAODExtension* ext = aodH->AddFilteredAOD(aodfilename,title,toMerge);
+
+ if (!ext) return;
+
+ if ( fOnlyMuon )
+ {
+
+ if(!murep)delete murep;
+
+ murep = new AliAODMCNuclExReplicator("NuclExReplicator",
+ "remove non interesting tracks",
+ fMCMode,fnSigmaTrk1,fnSigmaTrk2,fpartType1,fpartType2);
+
+
+ ext->DropUnspecifiedBranches(); // all branches not part of a FilterBranch call (below) will be dropped
+
+ ext->FilterBranch("header",murep);
+ ext->FilterBranch("vertices",murep);
+ ext->FilterBranch("nuclei",murep);
+ ext->FilterBranch("secvertices",murep); //per test
+ ext->FilterBranch("daughtertracks",murep);
+
+ //cout<<"add filterd aod"<<endl;
+
+ if ( fMCMode > 0 )
+ {
+ // MC branches will be copied (if present), as they are, but only
+ // for events with at least one muon.
+ // For events w/o muon, mcparticles array will be empty and mcheader will be dummy
+ // (e.g. strlen(GetGeneratorName())==0)
+
+ ext->FilterBranch("mcparticles",murep);
+ ext->FilterBranch("mcHeader",murep);
+ }
+ }
+
+ }
+ //cout<<"fine add filterd"<<endl;
+}
+
+//______________________________________________________________________________
+void AliAnalysisTaskESDNuclExFilterMC::Init()
+{
+
+ // Initialization
+ if(fEnableMuonAOD)
+ AddFilteredAOD("AliAOD.NuclEx.root", "NuclexFilteredEvents",kTRUE);
+}
+
+
+//______________________________________________________________________________
+void AliAnalysisTaskESDNuclExFilterMC::UserExec(Option_t */*option*/)
+{
+ // Execute analysis for current event
+
+ Long64_t ientry = Entry();
+ if(fDebug)printf("Muon Filter: Analysing event # %5d\n", (Int_t) ientry);
+
+ //***************************************************
+
+ ConvertESDtoAOD();
+
+ //*************************************************
+
+}
+
+
+void AliAnalysisTaskESDNuclExFilterMC::ConvertESDtoAOD()
+
+{
+
+ //cout<<"========================> CONVERT ESD TO AOD <============================="<<endl;
+
+
+ AliAODEvent *lAODevent=(AliAODEvent*)InputEvent();
+
+ // Read primary vertex from AOD event
+ // AliAODVertex *primary = *(AODEvent()->GetPrimaryVertex());
+
+ AliAODVertex *primary = lAODevent->GetPrimaryVertex();
+ if (fDebug && primary) primary->Print();
+ //cout<<"Primary vtx x: "<<primary->GetX()<<" "<<primary->GetY()<<" "<<primary->GetZ()<<endl;
+
+ AliAODHandler* handler = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
+
+ if ( handler ){
+
+ AliAODExtension *extNuclEx = handler->GetFilteredAOD("AliAOD.NuclEx.root");
+
+ if ( extNuclEx ) {
+
+ extNuclEx->SetEvent(lAODevent);
+ extNuclEx->SelectEvent();
+ extNuclEx->FinishEvent();
+
+ }
+ }
+
+}
+//------------------------------------------------
+void AliAnalysisTaskESDNuclExFilterMC::Terminate(Option_t */*option*/)
+{
+ // Terminate analysis
+ //
+ // delete murep;
+
+ if (fDebug > 1) printf("AnalysisESDfilter: Terminate() \n");
+}
--- /dev/null
+#ifndef ALIANALYSISTASKESDNUCLEXFILTERMC_H
+#define ALIANALYSISTASKESDNUCLEXFILTERMC_H
+
+
+// Create a new AOD starting from the general AOD. This Task can be used also strating
+//from ESD changing the input handler. (Method to be testeted on the grid)
+// filtering of the ESD.
+//
+// Authors: S. Bufalino (stefania.bufalino@cern.ch)
+// R. Lea (ramona.lea@cern.ch)
+// Based on AliAnalysisTaskESDMuonFilter.h
+
+#ifndef ALIANALYSISTASKSE_H
+# include "AliAnalysisTaskSE.h"
+# include "AliESDtrack.h"
+# include "AliAODTrack.h"
+# include "AliAODPid.h"
+# include "AliAODMCNuclExReplicator.h"
+//#include "AliAOD3LH.h"
+
+#endif
+#ifndef ALIANALYSISCUTS_H
+# include "AliAnalysisCuts.h"
+#endif
+
+class AliAnalysisFilter;
+class AliESDtrack;
+//class AliPIDResponse;
+
+class AliAnalysisTaskESDNuclExFilterMC : public AliAnalysisTaskSE
+{
+ public:
+ AliAnalysisTaskESDNuclExFilterMC(Bool_t onlyMuon=kTRUE, Bool_t keepAllEvents=kTRUE, Int_t mcMode=0, Int_t nsigmaTrk1=3 ,Int_t nsigmaTrk2 = 3, Int_t partType1 = 2,Int_t partType2 = 7);
+ AliAnalysisTaskESDNuclExFilterMC(const char* name, Bool_t onlyMuon=kTRUE, Bool_t keepAllEvents=kTRUE, Int_t mcMode=0, Int_t nsigmaTrk1=3 ,Int_t nsigmaTrk2 = 3, Int_t partType1 = 2,Int_t partType2 = 7);
+ virtual ~AliAnalysisTaskESDNuclExFilterMC() {;}
+ // Implementation of interface methods
+ virtual void UserCreateOutputObjects();
+ virtual void Init();
+ virtual void LocalInit() {Init();}
+ virtual void UserExec(Option_t *option);
+ virtual void Terminate(Option_t *option);
+
+ virtual void ConvertESDtoAOD();
+
+
+ // Setters
+ virtual void SetTrackFilter(AliAnalysisFilter* trackF) {fTrackFilter = trackF;}
+ void SetWriteMuonAOD(Bool_t enableMuonAOD){fEnableMuonAOD = enableMuonAOD;}
+ void SetWriteDimuonAOD(Bool_t enableDimuonAOD){fEnableDimuonAOD = enableDimuonAOD;}
+
+ void PrintTask(Option_t *option="", Int_t indent=0) const;
+
+
+ private:
+ AliAnalysisTaskESDNuclExFilterMC(const AliAnalysisTaskESDNuclExFilterMC&);
+ AliAnalysisTaskESDNuclExFilterMC& operator=(const AliAnalysisTaskESDNuclExFilterMC&);
+ void AddFilteredAOD(const char* aodfilename, const char* title, Bool_t toMerge);
+
+ AliAnalysisFilter* fTrackFilter; // Track Filter
+ Bool_t fEnableMuonAOD; // flag for enabling Muon AOD production
+ Bool_t fEnableDimuonAOD; // flag for enabling Dimuon AOD production
+ Bool_t fOnlyMuon; // flag for disabling branches irrelevant for (most) muon analyses
+ Bool_t fKeepAllEvents; // keep even events where there's no muons (to get e.g. unbiased vertex distribution)
+ Int_t fMCMode; // whether and how we're filtering MC data
+
+ Int_t fnSigmaTrk1;
+ Int_t fnSigmaTrk2;
+ Int_t fpartType1;
+ Int_t fpartType2;
+
+ AliAODMCNuclExReplicator* murep;
+
+ /* AliPIDResponse *fPIDResponse; //! PID response object */
+
+ ClassDef(AliAnalysisTaskESDNuclExFilterMC, 5); // Analysis task for standard ESD filtering
+};
+
+/* class AliAnalysisNonMuonTrackCuts : public AliAnalysisCuts */
+/* { */
+/* public: */
+/* AliAnalysisNonMuonTrackCuts(); */
+/* virtual ~AliAnalysisNonMuonTrackCuts() {} */
+/* /\* virtual Bool_t IsSelected(TObject* obj); *\/ */
+/* /\* virtual Bool_t IsSelected(TList* /\\* list *\\/ ) { return kTRUE; } *\/ */
+
+/* ClassDef(AliAnalysisNonMuonTrackCuts,1); // Select muon spectrometer tracks */
+/* }; */
+
+/* class AliAnalysisNonPrimaryVertices : public AliAnalysisCuts */
+/* { */
+/* public: */
+/* AliAnalysisNonPrimaryVertices(); */
+/* virtual ~AliAnalysisNonPrimaryVertices() {} */
+/* /\* virtual Bool_t IsSelected(TObject* obj); *\/ */
+/* /\* virtual Bool_t IsSelected(TList* /\\* list *\\/ ) { return kTRUE; } *\/ */
+
+/* ClassDef(AliAnalysisNonPrimaryVertices,1); // Select primary vertices */
+/* }; */
+
+#endif
--- /dev/null
+/**************************************************************************
+ * Contributors are not mentioned at all. *
+ * *
+ * 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. *
+ **************************************************************************/
+
+//
+//-----------------------------------------------------------------
+// AliAnalysisTaskNUclexAOD class AOD
+//-----------------------------------------------------------------
+
+
+class TTree;
+class TParticle;
+class TVector3;
+
+#include "AliAnalysisManager.h"
+#include <AliMCEventHandler.h>
+#include <AliMCEvent.h>
+#include <AliStack.h>
+
+class AliESDVertex;
+class AliAODVertex;
+class AliESDv0;
+class AliAODv0;
+class AliCascadeVertexer;
+
+#include <iostream>
+#include "AliAnalysisTaskSE.h"
+#include "TList.h"
+#include "TH1.h"
+#include "TH2.h"
+#include "TH3.h"
+#include "TNtuple.h"
+#include "TGraph.h"
+#include "TCutG.h"
+#include "TF1.h"
+#include "TCanvas.h"
+#include "TMath.h"
+#include "TChain.h"
+#include "Riostream.h"
+#include "AliLog.h"
+#include "AliCascadeVertexer.h"
+#include "AliESDEvent.h"
+#include "AliESDtrack.h"
+#include "AliExternalTrackParam.h"
+#include "AliAODEvent.h"
+#include "AliInputEventHandler.h"
+#include "AliESDcascade.h"
+#include "AliAODcascade.h"
+#include "AliAnalysisTaskReadNuclexAOD.h"
+#include "AliESDtrackCuts.h"
+#include "AliCentrality.h"
+#include "TString.h"
+#include <TDatime.h>
+#include <TRandom3.h>
+#include <AliAODPid.h>
+#include <AliAODHandler.h>
+#include <AliAODRecoDecay.h>
+#include <AliAODRecoDecayLF.h>
+#include <AliAODRecoDecayLF2Prong.h>
+
+ClassImp(AliAnalysisTaskReadNuclexAOD)
+
+//________________________________________________________________________
+AliAnalysisTaskReadNuclexAOD::AliAnalysisTaskReadNuclexAOD()
+ : AliAnalysisTaskSE(),
+ fAnalysisType("AOD"),
+ fCollidingSystems(0),
+ fDataType("REAL"),
+ fListHist(0),
+ fHistEventMultiplicity(0),
+ fHistTrackMultiplicity(0),
+ fHistTrackMultiplicityCent(0),
+ fHistTrackMultiplicitySemiCent(0),
+ fHistTrackMultiplicityMB(0),
+ fHistTrackMultiplicityPVCent(0),
+ fHistTrackMultiplicityPVSemiCent(0),
+ fHistTrackMultiplicityPVMB(0),
+ fhBB(0),
+ fhBBPions(0),
+ fhBBHe(0),
+ aodTree(0x0),
+ Nuclei(0x0),
+ SecondaryVertices(0x0),
+ DaughterTracks(0x0),
+ nevent(0x0),
+ fNtuple1(0),
+ trunNumber(0),
+ tbunchcross(0),
+ torbit(0),
+ tperiod(0),
+ teventtype(0),
+ tTrackNumber(0),
+ tpercentile(0),
+ txPrimaryVertex(0),
+ tyPrimaryVertex(0),
+ tzPrimaryVertex(0),
+ txSecondaryVertex(0),
+ tySecondaryVertex(0),
+ tzSecondaryVertex(0),
+ tdcaTracks(0),
+ tCosPointingAngle(0),
+ tDCAV0toPrimaryVertex(0),
+ tHeSign(0),
+ tHepInTPC(0),
+ tHeTPCsignal(0),
+ tDcaHeToPrimVertex(0),
+ tHeEta(0),
+ tmomHex(0),
+ tmomHey(0),
+ tmomHez(0),
+ tmomHeAtSVx(0),
+ tmomHeAtSVy(0),
+ tmomHeAtSVz(0),
+ tHeTPCNcls(0),
+ tHeimpactXY(0),
+ tHeimpactZ(0),
+ tHeITSClusterMap(0),
+ tIsHeITSRefit(0),
+ tPionSign(0),
+ tPionpInTPC(0),
+ tPionTPCsignal(0),
+ tDcaPionToPrimVertex(0),
+ tPionEta(0),
+ tmomPionx(0),
+ tmomPiony(0),
+ tmomPionz(0),
+ tmomNegPionAtSVx(0),
+ tmomNegPionAtSVy(0),
+ tmomNegPionAtSVz(0),
+ tPionTPCNcls(0),
+ tPionimpactXY(0),
+ tPionimpactZ(0),
+ tPionITSClusterMap(0),
+ tIsPiITSRefit(0),
+ txn(0),
+ txp(0),
+ tchi2He(0),
+ tchi2Pi(0),
+ fNtuple4(0),
+ t3LHlrunNumber(0),
+ t3LHlBCNumber(0),
+ t3LHlOrbitNumber(0),
+ t3LHlPeriodNumber(0),
+ t3LHleventtype(0),
+ t3LHlpercentile(0),
+ t3LHlPx(0),
+ t3LHlPy(0),
+ t3LHlPz(0),
+ t3LHlEta(0),
+ t3LHlY(0),
+ t3LHlM(0),
+ t3LHlxPrimaryVertex(0),
+ t3LHlyPrimaryVertex(0),
+ t3LHlzPrimaryVertex(0),
+ t3LHlp0px(0),
+ t3LHlp0py(0),
+ t3LHlp0pz(0),
+ t3LHlp0ch(0),
+ t3LHlp1px(0),
+ t3LHlp1py(0),
+ t3LHlp1pz(0),
+ t3LHlp1ch(0)
+
+{
+ // Dummy Constructor
+}
+
+//________________________________________________________________________
+AliAnalysisTaskReadNuclexAOD::AliAnalysisTaskReadNuclexAOD(const char *name)
+ : AliAnalysisTaskSE(name),
+ fAnalysisType("AOD"),
+ fCollidingSystems(0),
+ fDataType("REAL"),
+ fListHist(0),
+ fHistEventMultiplicity(0),
+ fHistTrackMultiplicity(0),
+ fHistTrackMultiplicityCent(0),
+ fHistTrackMultiplicitySemiCent(0),
+ fHistTrackMultiplicityMB(0),
+ fHistTrackMultiplicityPVCent(0),
+ fHistTrackMultiplicityPVSemiCent(0),
+ fHistTrackMultiplicityPVMB(0),
+ fhBB(0),
+ fhBBPions(0),
+ fhBBHe(0),
+ aodTree(0x0),
+ Nuclei(0x0),
+ SecondaryVertices(0x0),
+ DaughterTracks(0x0),
+ nevent(0x0),
+ fNtuple1(0),
+ trunNumber(0),
+ tbunchcross(0),
+ torbit(0),
+ tperiod(0),
+ teventtype(0),
+ tTrackNumber(0),
+ tpercentile(0),
+ txPrimaryVertex(0),
+ tyPrimaryVertex(0),
+ tzPrimaryVertex(0),
+ txSecondaryVertex(0),
+ tySecondaryVertex(0),
+ tzSecondaryVertex(0),
+ tdcaTracks(0),
+ tCosPointingAngle(0),
+ tDCAV0toPrimaryVertex(0),
+ tHeSign(0),
+ tHepInTPC(0),
+ tHeTPCsignal(0),
+ tDcaHeToPrimVertex(0),
+ tHeEta(0),
+ tmomHex(0),
+ tmomHey(0),
+ tmomHez(0),
+ tmomHeAtSVx(0),
+ tmomHeAtSVy(0),
+ tmomHeAtSVz(0),
+ tHeTPCNcls(0),
+ tHeimpactXY(0),
+ tHeimpactZ(0),
+ tHeITSClusterMap(0),
+ tIsHeITSRefit(0),
+ tPionSign(0),
+ tPionpInTPC(0),
+ tPionTPCsignal(0),
+ tDcaPionToPrimVertex(0),
+ tPionEta(0),
+ tmomPionx(0),
+ tmomPiony(0),
+ tmomPionz(0),
+ tmomNegPionAtSVx(0),
+ tmomNegPionAtSVy(0),
+ tmomNegPionAtSVz(0),
+ tPionTPCNcls(0),
+ tPionimpactXY(0),
+ tPionimpactZ(0),
+ tPionITSClusterMap(0),
+ tIsPiITSRefit(0),
+ txn(0),
+ txp(0),
+ tchi2He(0),
+ tchi2Pi(0),
+ fNtuple4(0),
+ t3LHlrunNumber(0),
+ t3LHlBCNumber(0),
+ t3LHlOrbitNumber(0),
+ t3LHlPeriodNumber(0),
+ t3LHleventtype(0),
+ t3LHlpercentile(0),
+ t3LHlPx(0),
+ t3LHlPy(0),
+ t3LHlPz(0),
+ t3LHlEta(0),
+ t3LHlY(0),
+ t3LHlM(0),
+ t3LHlxPrimaryVertex(0),
+ t3LHlyPrimaryVertex(0),
+ t3LHlzPrimaryVertex(0),
+ t3LHlp0px(0),
+ t3LHlp0py(0),
+ t3LHlp0pz(0),
+ t3LHlp0ch(0),
+ t3LHlp1px(0),
+ t3LHlp1py(0),
+ t3LHlp1pz(0),
+ t3LHlp1ch(0)
+
+{
+ // Define input and output slots here
+ // Input slot #0 works with a TChain
+ //DefineInput(0, TChain::Class());
+ // Output slot #0 writes into a TList container (Cascade)
+ DefineOutput(1, TList::Class());
+
+ DefineInput(0, TChain::Class());
+
+ DefineOutput(1, TList::Class());
+ DefineOutput(2, TTree::Class());
+ DefineOutput(3, TTree::Class());
+
+}
+//_______________________________________________________
+AliAnalysisTaskReadNuclexAOD::~AliAnalysisTaskReadNuclexAOD()
+{
+ // Destructor
+ if (fListHist) {
+ delete fListHist;
+ fListHist = 0;
+ }
+
+ if(fNtuple1) delete fNtuple1;
+ if(fNtuple4) delete fNtuple4;
+}
+//=================DEFINITION BETHE BLOCH==============================
+
+Double_t AliAnalysisTaskReadNuclexAOD::BetheBloch(Double_t betaGamma,Double_t charge,Bool_t isPbPb) {
+
+ Double_t kp1, kp2, kp3, kp4, kp5;
+
+ if(isPbPb){
+
+
+
+ //pass2 2011
+ kp1 = 4.7*charge*charge;
+ kp2 = 8.98482806165147636e+00;
+ kp3 = 1.54000000000000005e-05;
+ kp4 = 2.30445734159456084e+00;
+ kp5 = 2.25624744086878559e+00;
+
+ }
+
+ else{
+
+
+ //pass2 2011
+ kp1 = 4.7*charge*charge;
+ kp2 = 8.98482806165147636e+00;
+ kp3 = 1.54000000000000005e-05;
+ kp4 = 2.30445734159456084e+00;
+ kp5 = 2.25624744086878559e+00;
+
+
+ }
+
+ Double_t beta = betaGamma / TMath::Sqrt(1.0 + betaGamma * betaGamma);
+
+ Double_t aa = TMath::Power(beta, kp4);
+ Double_t bb = TMath::Power(1.0 / betaGamma, kp5);
+
+ bb = TMath::Log(kp3 + bb);
+
+ Double_t out = (kp2 - aa - bb) * kp1 / aa;
+
+ return out;
+
+}
+
+void AliAnalysisTaskReadNuclexAOD::Init()
+{
+ // cout<<"----------------------------------> Vediamo"<<endl;
+
+ // AliAODHandler* aodHandler = (AliAODHandler*)((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
+ // cout<<aodHandler<<endl;
+ // TTree *aodTree = aodHandler->GetTree();
+ // // aodTree = aodHandler->GetTree();
+ // cout<<"aodTree: "<<aodTree<<endl;
+ // aodTree->SetBranchAddress("Nuclei", &Nuclei);
+ // aodTree->SetBranchAddress("SecondaryVertices",&SecondaryVertices);
+ // aodTree->SetBranchAddress("DaughterTracks", &DaughterTracks);
+
+
+}
+
+//==================DEFINITION OF OUTPUT OBJECTS==============================
+
+void AliAnalysisTaskReadNuclexAOD::UserCreateOutputObjects()
+{
+ fListHist = new TList();
+ fListHist->SetOwner(); // IMPORTANT!
+
+
+ if(! fHistEventMultiplicity ){
+
+ fHistEventMultiplicity = new TH1F( "fHistEventMultiplicity" , "Nb of Events" , 13 , 0, 13);
+
+ fHistEventMultiplicity->GetXaxis()->SetBinLabel(1,"All Events");
+ fHistEventMultiplicity->GetXaxis()->SetBinLabel(2,"Events w/PV");
+ fHistEventMultiplicity->GetXaxis()->SetBinLabel(3,"Events w/|Vz|<10cm");
+ fHistEventMultiplicity->GetXaxis()->SetBinLabel(4,"Central Events");
+ fHistEventMultiplicity->GetXaxis()->SetBinLabel(5,"SemiCentral Events");
+ fHistEventMultiplicity->GetXaxis()->SetBinLabel(6,"MB Events");
+ fHistEventMultiplicity->GetXaxis()->SetBinLabel(7,"Central Events w/|Vz|<10cm");
+ fHistEventMultiplicity->GetXaxis()->SetBinLabel(8,"SemiCentral Events w/|Vz|<10cm");
+ fHistEventMultiplicity->GetXaxis()->SetBinLabel(9,"MB Events w/|Vz|<10cm");
+ fHistEventMultiplicity->GetXaxis()->SetBinLabel(10,"Sec. Vert. Tree");
+ fHistEventMultiplicity->GetXaxis()->SetBinLabel(11,"Sec. Vert. Rec");
+ fHistEventMultiplicity->GetXaxis()->SetBinLabel(12,"Event Not taken");
+ fListHist->Add(fHistEventMultiplicity);
+ }
+
+
+ if(! fHistTrackMultiplicity ){
+ fHistTrackMultiplicity = new TH2F( "fHistTrackMultiplicity" , "Nb of Tracks", 2500,0, 25000,210,-1,104);
+ fHistTrackMultiplicity->GetXaxis()->SetTitle("Number of tracks");
+ fHistTrackMultiplicity->GetYaxis()->SetTitle("Percentile");
+ fListHist->Add(fHistTrackMultiplicity);
+ }
+
+ if(! fHistTrackMultiplicityCent ){
+ fHistTrackMultiplicityCent = new TH2F( "fHistTrackMultiplicityCent", "Nb of Tracks Central Events", 2500,0, 25000,210,-1,104 );
+ fHistTrackMultiplicityCent->GetXaxis()->SetTitle("Number of tracks");
+ fHistTrackMultiplicityCent->GetYaxis()->SetTitle("Percentile");
+ fListHist->Add(fHistTrackMultiplicityCent);
+ }
+
+ if(! fHistTrackMultiplicitySemiCent ){
+ fHistTrackMultiplicitySemiCent = new TH2F( "fHistTrackMultiplicitySemiCent" , "Nb of Tracks SemiCentral Events", 2500,0, 25000 ,210,-1,104);
+ fHistTrackMultiplicitySemiCent->GetXaxis()->SetTitle("Number of tracks");
+ fHistTrackMultiplicitySemiCent->GetYaxis()->SetTitle("Percentile");
+ fListHist->Add(fHistTrackMultiplicitySemiCent);
+ }
+
+ if(! fHistTrackMultiplicityMB ){
+ fHistTrackMultiplicityMB = new TH2F( "fHistTrackMultiplicityMB" , "Nb of Tracks MBral Events", 2500,0, 25000,210,-1,104 );
+ fHistTrackMultiplicityMB->GetXaxis()->SetTitle("Number of tracks");
+ fHistTrackMultiplicityMB->GetYaxis()->SetTitle("Percentile");
+ fListHist->Add(fHistTrackMultiplicityMB);
+ }
+
+ if(! fHistTrackMultiplicityPVCent ){
+ fHistTrackMultiplicityPVCent = new TH2F( "fHistTrackMultiplicityPVCent" , "Nb of Tracks Central Events", 2500,0, 25000,210,-1,104 );
+ fHistTrackMultiplicityPVCent->GetXaxis()->SetTitle("Number of tracks");
+ fHistTrackMultiplicityPVCent->GetYaxis()->SetTitle("Percentile");
+ fListHist->Add(fHistTrackMultiplicityPVCent);
+ }
+
+ if(! fHistTrackMultiplicityPVSemiCent ){
+ fHistTrackMultiplicityPVSemiCent = new TH2F( "fHistTrackMultiplicityPVSemiCent" , "Nb of Tracks SemiCentral Events", 2500,0, 25000 ,210,-1,104);
+ fHistTrackMultiplicityPVSemiCent->GetXaxis()->SetTitle("Number of tracks");
+ fHistTrackMultiplicityPVSemiCent->GetYaxis()->SetTitle("Percentile");
+ fListHist->Add(fHistTrackMultiplicityPVSemiCent);
+ }
+
+ if(! fHistTrackMultiplicityPVMB ){
+ fHistTrackMultiplicityPVMB = new TH2F( "fHistTrackMultiplicityPVMB" , "Nb of Tracks MBral Events", 2500,0, 25000,210,-1,104 );
+ fHistTrackMultiplicityPVMB->GetXaxis()->SetTitle("Number of tracks");
+ fHistTrackMultiplicityPVMB->GetYaxis()->SetTitle("Percentile");
+ fListHist->Add(fHistTrackMultiplicityPVMB);
+ }
+
+ if(! fhBB ){
+ fhBB = new TH2F( "fhBB" , "BetheBlochTPC" , 120,-6,6,150,0,1500);
+ fhBB->GetXaxis()->SetTitle("p/z (GeV/#it{c})");
+ fhBB->GetYaxis()->SetTitle("TPC Signal");
+ fListHist->Add(fhBB);
+ }
+
+ if(! fhBBPions ){
+ fhBBPions = new TH2F( "fhBBPions" , "Bethe-Bloch TPC Pions" , 120,-6,6,150,0,1500);
+ fhBBPions->GetXaxis()->SetTitle("p/z (GeV/#it{c})");
+ fhBBPions->GetYaxis()->SetTitle("TPC Signal");
+ fListHist->Add(fhBBPions);
+ }
+
+ if(! fhBBHe ){
+ fhBBHe = new TH2F( "fhBBHe" , "Bethe-Bloch TPC He" , 120,-6,6,150,0,1500);
+ fhBBHe->GetXaxis()->SetTitle("p/z (GeV/#it{c})");
+ fhBBHe->GetYaxis()->SetTitle("TPC Signal");
+ fListHist->Add(fhBBHe);
+ }
+
+
+
+ if(! fNtuple1 ) {
+
+ fNtuple1 = new TTree("fNtuple1","fNtuple1");
+
+ fNtuple1->Branch("trunNumber" ,&trunNumber ,"trunNumber/F");
+ fNtuple1->Branch("tbunchcross" ,&tbunchcross ,"tbunchcross/F");
+ fNtuple1->Branch("torbit" ,&torbit ,"torbit/F");
+ fNtuple1->Branch("tperiod" ,&tperiod ,"tperiod/F");
+ fNtuple1->Branch("teventtype" ,&teventtype ,"teventtype/F");
+ fNtuple1->Branch("tTrackNumber" ,&tTrackNumber ,"tTrackNumber/F");
+ fNtuple1->Branch("tpercentile" ,&tpercentile ,"tpercentile/F") ;
+ fNtuple1->Branch("txPrimaryVertex" ,&txPrimaryVertex ,"txPrimaryVertex/F");
+ fNtuple1->Branch("tyPrimaryVertex" ,&tyPrimaryVertex ,"tyPrimaryVertex/F");
+ fNtuple1->Branch("tzPrimaryVertex" ,&tzPrimaryVertex ,"tzPrimaryVertex/F");
+ fNtuple1->Branch("txSecondaryVertex" ,&txSecondaryVertex ,"txSecondaryVertex/F");
+ fNtuple1->Branch("tySecondaryVertex" ,&tySecondaryVertex ,"tySecondaryVertex/F");
+ fNtuple1->Branch("tzSecondaryVertex" ,&tzSecondaryVertex ,"tzSecondaryVertex/F");
+ fNtuple1->Branch("tdcaTracks" ,&tdcaTracks ,"tdcaTracks/F");
+ fNtuple1->Branch("tCosPointingAngle" ,&tCosPointingAngle ,"tCosPointingAngle/F");
+ fNtuple1->Branch("tDCAV0toPrimaryVertex",&tDCAV0toPrimaryVertex,"tDCAV0toPrimaryVertex/F");
+ fNtuple1->Branch("tHeSign" ,&tHeSign ,"tHeSign/F");
+ fNtuple1->Branch("tHepInTPC" ,&tHepInTPC ,"tHepInTPC/F");
+ fNtuple1->Branch("tHeTPCsignal" ,&tHeTPCsignal ,"tHeTPCsignal/F");
+ fNtuple1->Branch("tDcaHeToPrimVertex" ,&tDcaHeToPrimVertex ,"tDcaHeToPrimVertex/F");
+ fNtuple1->Branch("tHeEta" ,&tHeEta ,"tHeEta/F");
+ fNtuple1->Branch("tmomHex" ,&tmomHex ,"tmomHex/F");
+ fNtuple1->Branch("tmomHey" ,&tmomHey ,"tmomHey/F");
+ fNtuple1->Branch("tmomHez" ,&tmomHez ,"tmomHez/F");
+ fNtuple1->Branch("tmomHeAtSVx" ,&tmomHeAtSVx ,"tmomHeAtSVx/F");
+ fNtuple1->Branch("tmomHeAtSVy" ,&tmomHeAtSVy ,"tmomHeAtSVy/F");
+ fNtuple1->Branch("tmomHeAtSVz" ,&tmomHeAtSVz ,"tmomHeAtSVz/F");
+ fNtuple1->Branch("tHeTPCNcls" ,&tHeTPCNcls ,"tHeTPCNcls/F");
+ fNtuple1->Branch("tHeimpactXY" ,&tHeimpactXY ,"tHeimpactXY/F");
+ fNtuple1->Branch("tHeimpactZ" ,&tHeimpactZ ,"tHeimpactZ/F");
+ fNtuple1->Branch("tHeITSClusterMap" ,&tHeITSClusterMap ,"tHeITSClusterMap/F");
+ fNtuple1->Branch("tIsHeITSRefit" ,&tIsHeITSRefit ,"tIsHeITSRefit/F");
+ fNtuple1->Branch("tPionSign" ,&tPionSign ,"tPionSign/F");
+ fNtuple1->Branch("tPionpInTPC" ,&tPionpInTPC ,"tPionpInTPC/F");
+ fNtuple1->Branch("tPionTPCsignal" ,&tPionTPCsignal ,"tPionTPCsignal/F");
+ fNtuple1->Branch("tDcaPionToPrimVertex" ,&tDcaPionToPrimVertex ,"tDcaPionToPrimVertex/F");
+ fNtuple1->Branch("tPionEta" ,&tPionEta ,"tPionEta/F");
+ fNtuple1->Branch("tmomPionx" ,&tmomPionx ,"tmomPionx/F");
+ fNtuple1->Branch("tmomPiony" ,&tmomPiony ,"tmomPiony/F");
+ fNtuple1->Branch("tmomPionz" ,&tmomPionz ,"tmomPionz/F");
+ fNtuple1->Branch("tmomNegPionAtSVx" ,&tmomNegPionAtSVx ,"tmomNegPionAtSVx/F");
+ fNtuple1->Branch("tmomNegPionAtSVy" ,&tmomNegPionAtSVy ,"tmomNegPionAtSVy/F");
+ fNtuple1->Branch("tmomNegPionAtSVz" ,&tmomNegPionAtSVz ,"tmomNegPionAtSVz/F");
+ fNtuple1->Branch("tPionTPCNcls" ,&tPionTPCNcls ,"tPionTPCNcls/F");
+ fNtuple1->Branch("tPionimpactXY" ,&tPionimpactXY ,"tPionimpactXY/F");
+ fNtuple1->Branch("tPionimpactZ" ,&tPionimpactZ ,"tPionimpactZ/F");
+ fNtuple1->Branch("tPionITSClusterMap" ,&tPionITSClusterMap ,"tPionITSClusterMap/F");
+ fNtuple1->Branch("tIsPiITSRefit" ,&tIsPiITSRefit ,"tIsPiITSRefit/F");
+ fNtuple1->Branch("txn" ,&txn ,"txn/F");
+ fNtuple1->Branch("txp" ,&txp ,"txp/F");
+ fNtuple1->Branch("tchi2He" ,&tchi2He ,"tchi2He/F");
+ fNtuple1->Branch("tchi2Pi" ,&tchi2Pi ,"tchi2Pi/F");
+
+ }
+
+ if(! fNtuple4 ) {
+
+ fNtuple4 = new TTree("fNtuple4","fNtuple4");
+
+ fNtuple4->Branch("t3LHlrunNumber" ,&t3LHlrunNumber ,"t3LHlrunNumber/F");
+ fNtuple4->Branch("t3LHlBCNumber" ,&t3LHlBCNumber ,"t3LHlBCNumber/F");
+ fNtuple4->Branch("t3LHlOrbitNumber" ,&t3LHlOrbitNumber ,"t3LHlOrbitNumber/F");
+ fNtuple4->Branch("t3LHlPeriodNumber" ,&t3LHlPeriodNumber ,"t3LHlPeriodNumber/F");
+ fNtuple4->Branch("t3LHleventtype" ,&t3LHleventtype ,"t3LHleventtype/F");
+ fNtuple4->Branch("t3LHlpercentile" ,&t3LHlpercentile ,"t3LHlpercentile/F");
+ fNtuple4->Branch("t3LHlPx" ,&t3LHlPx ,"t3LHlPx/F");
+ fNtuple4->Branch("t3LHlPy" ,&t3LHlPy ,"t3LHlPy/F");
+ fNtuple4->Branch("t3LHlPz" ,&t3LHlPz ,"t3LHlPz/F");
+ fNtuple4->Branch("t3LHlEta" ,&t3LHlEta ,"t3LHlEta/F");
+ fNtuple4->Branch("t3LHlY" ,&t3LHlY ,"t3LHlY/F");
+ fNtuple4->Branch("t3LHlM" ,&t3LHlM ,"t3LHlM/F");
+ fNtuple4->Branch("t3LHlxPrimaryVertex" ,&t3LHlxPrimaryVertex ,"t3LHlxPrimaryVertex/F");
+ fNtuple4->Branch("t3LHlyPrimaryVertex" ,&t3LHlyPrimaryVertex ,"t3LHlyPrimaryVertex/F");
+ fNtuple4->Branch("t3LHlzPrimaryVertex" ,&t3LHlzPrimaryVertex ,"t3LHlzPrimaryVertex/F");
+ fNtuple4->Branch("t3LHlp0px" ,&t3LHlp0px ,"t3LHlp0px/F");
+ fNtuple4->Branch("t3LHlp0py" ,&t3LHlp0py ,"t3LHlp0py/F");
+ fNtuple4->Branch("t3LHlp0pz" ,&t3LHlp0pz ,"t3LHlp0pz/F");
+ fNtuple4->Branch("t3LHlp0ch" ,&t3LHlp0ch ,"t3LHlp0ch/F");
+ fNtuple4->Branch("t3LHlp1px" ,&t3LHlp1px ,"t3LHlp1px/F");
+ fNtuple4->Branch("t3LHlp1py" ,&t3LHlp1py ,"t3LHlp1py/F");
+ fNtuple4->Branch("t3LHlp1pz" ,&t3LHlp1pz ,"t3LHlp1pz/F");
+ fNtuple4->Branch("t3LHlp1ch" ,&t3LHlp1ch ,"t3LHlp1ch/F");
+
+
+ }
+
+
+ PostData(1, fListHist);
+ PostData(2, fNtuple1);
+ PostData(3, fNtuple4);
+}// end UserCreateOutputObjects
+
+
+
+//====================== USER EXEC ========================
+
+void AliAnalysisTaskReadNuclexAOD::UserExec(Option_t *)
+{
+ //_______________________________________________________________________
+
+ //!*********************!//
+ //! Define variables !//
+ //!*********************!//
+
+ Double_t xPrimaryVertex=0.,yPrimaryVertex=0.,zPrimaryVertex=0.;
+
+ ULong_t statusT;
+ ULong_t statusPi;
+
+ Bool_t IsHeITSRefit=kFALSE,IsPiITSRefit=kFALSE ;
+
+ Double_t fPos[3]={0.,0.,0.};
+
+ Double_t runNumber=0.;
+ Double_t BCNumber=0.;
+ Double_t OrbitNumber=0.;
+ Double_t PeriodNumber=0.;
+
+ Double_t Helium3Mass = 2.80839;
+ Double_t PionMass = 0.13957;
+
+ // TLORENTZ vectors
+
+ TLorentzVector vPion,vHelium,vSum;
+ TLorentzVector vHyp;
+ TLorentzVector vp0,vp1,vS;
+ //!----------------------------------------------------------------
+ //! A set of very loose parameters for cuts
+
+ // Double_t fgChi2max=33.; //! max chi2
+ Double_t fgDNmin=0.05; //! min imp parameter for the 1st daughter = 500um
+ // Double_t fgDPmin=0.05; //! min imp parameter for the 2nd daughter = 500um
+ Double_t fgDCAmax=1.; //! max DCA between the daughter tracks in cm
+ // Double_t fgCPAmin=0.99; //! min cosine of V0's pointing angle
+ Double_t fgRmin=0.1; //! min radius of the fiducial volume = 1 mm
+ Double_t fgRmax=200.; //! max radius of the fiducial volume = 2 m
+
+ //!-----------------------------------------------------------------
+
+ // Main loop
+ // Called for EACH event
+
+ nevent++; //--> This is important
+
+ AliVEvent *event = InputEvent();
+ if (!event) { Printf("ERROR: Could not retrieve event"); return; }
+ Info("AliAnalysisTaskReadNuclexAOD","Starting UserExec");
+
+ SetDataType("REAL");
+
+ AliAODEvent *lAODevent=(AliAODEvent *)InputEvent();
+ if (!lAODevent) {
+ Printf("ERROR: aod not available");
+ // fHistLog->Fill(98);
+ return;
+ }
+
+ fHistEventMultiplicity->Fill(0.5);
+
+ //*************************************************************
+
+ runNumber = lAODevent->GetRunNumber();
+ BCNumber = lAODevent->GetBunchCrossNumber();
+ OrbitNumber = lAODevent->GetOrbitNumber();
+ PeriodNumber = lAODevent->GetPeriodNumber();
+
+ AliCentrality *centrality = lAODevent->GetCentrality();
+ centrality = lAODevent->GetCentrality();
+
+ Float_t percentile=centrality->GetCentralityPercentile("V0M");
+ Float_t refMult = lAODevent->GetHeader()->GetRefMultiplicity();
+
+ fHistTrackMultiplicity->Fill(refMult,percentile); //tracce per evento
+
+ // cout<<"Centrality: "<<percentile<<" number of track: "<<lAODevent->GetHeader()->GetRefMultiplicity()<<endl;
+ //*************************************************************
+
+ AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
+ AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
+
+ Int_t eventtype=-99;
+
+ Bool_t isSelectedCentral = (inputHandler->IsEventSelected() & AliVEvent::kCentral);
+ Bool_t isSelectedSemiCentral = (inputHandler->IsEventSelected() & AliVEvent::kSemiCentral);
+ Bool_t isSelectedMB = (inputHandler->IsEventSelected() & AliVEvent::kMB);
+
+ if(isSelectedCentral){
+ fHistEventMultiplicity->Fill(3.5);
+ fHistTrackMultiplicityCent->Fill(refMult,percentile);
+ eventtype=1;
+ }
+
+ if(isSelectedSemiCentral){
+ fHistEventMultiplicity->Fill(4.5);
+ fHistTrackMultiplicitySemiCent->Fill(refMult,percentile);
+ eventtype=2;
+ }
+
+ if(isSelectedMB){
+ fHistEventMultiplicity->Fill(5.5);
+ fHistTrackMultiplicityMB->Fill(refMult,percentile);
+ eventtype=3;
+ }
+
+ if(eventtype==-99)return;
+ // cout<<"eventtype: "<<eventtype<<endl;
+
+ // if(!isSelectedCentral)return;
+ // if(!isSelectedSemiCentral)return;
+ // if(!isSelectedMB)return;
+
+ cout<<"eventtype: "<<eventtype<<endl;
+
+ //-------------------------------
+
+ Double_t lBestPrimaryVtxPos[3] = {-100.0, -100.0, -100.0};
+ Double_t lMagneticField = -10.;
+
+ AliAODVertex* primvtx = lAODevent->GetPrimaryVertex();
+ if(!primvtx){
+ cout<<"No PV"<<endl;
+ return;
+ }
+ fHistEventMultiplicity->Fill(1.5);
+ primvtx->GetXYZ( lBestPrimaryVtxPos );
+
+ if(TMath::Abs(lBestPrimaryVtxPos[2]) > 10.0 ) {
+ AliWarning("Pb / | Z position of Best Prim Vtx | > 10.0 cm ... return !");
+ return;
+ }
+ fHistEventMultiplicity->Fill(2.5);
+
+
+ if(eventtype==1){
+ fHistTrackMultiplicityPVCent->Fill(refMult,percentile);
+ fHistEventMultiplicity->Fill(6.5);
+ }
+
+ if(eventtype==2){
+ fHistTrackMultiplicityPVSemiCent->Fill(refMult,percentile);
+ fHistEventMultiplicity->Fill(7.5);
+ }
+
+ if(eventtype==3){
+ fHistTrackMultiplicityPVMB->Fill(refMult,percentile);
+ fHistEventMultiplicity->Fill(8.5);
+ }
+
+ lMagneticField = lAODevent->GetMagneticField();
+
+ xPrimaryVertex = lBestPrimaryVtxPos[0];
+ yPrimaryVertex = lBestPrimaryVtxPos[1];
+ zPrimaryVertex = lBestPrimaryVtxPos[2];
+
+ //--------------------------------------------------------
+
+ AliAODHandler* aodHandler = (AliAODHandler*)((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
+
+ aodTree = aodHandler->GetTree();
+ aodTree->SetBranchAddress("Nuclei",&Nuclei);
+ aodTree->SetBranchAddress("SecondaryVertices",&SecondaryVertices);
+ aodTree->SetBranchAddress("DaughterTracks",&DaughterTracks);
+
+ //---------------------------------------------------------
+
+ aodTree->GetEntry(nevent);
+
+ if(Nuclei-> GetEntriesFast()>0)
+ fHistEventMultiplicity->Fill(2);
+
+ if(SecondaryVertices-> GetEntriesFast()>0)
+ fHistEventMultiplicity->Fill(3);
+
+ for(Int_t i=0; i < SecondaryVertices-> GetEntriesFast() ; i++){
+
+ AliAODRecoDecayLF2Prong *d = (AliAODRecoDecayLF2Prong*)SecondaryVertices->UncheckedAt(i);
+ //cout<<"\n\n\nPT: "<<d->Pt()<<endl<<endl<<endl<<endl;
+ if(!d) continue;
+
+ fHistEventMultiplicity->Fill(10.5);
+
+ // Double_t tDecayVertexV0[3];
+ // if(d->GetXYZ(tDecayVertexV0))
+
+ Double_t tV0mom[3];
+ d->PxPyPz(tV0mom);
+
+ // TLORENTZ vectors
+
+ // cout<<"px prong 1: "<<d->PxProng(0)<<endl;
+
+ vHyp.SetXYZM(tV0mom[0],tV0mom[1],tV0mom[2],2.991);
+
+
+ vp0.SetXYZM(d->PxProng(0),d->PyProng(0),d->PzProng(0),PionMass);
+ vp1.SetXYZM(d->PxProng(1),d->PyProng(1),d->PzProng(1),Helium3Mass);
+ vS = vp0 + vp1;
+
+ // cout<<"Momemnto pt: "<<lPt<<endl;
+
+ t3LHlrunNumber =(Float_t)runNumber;
+ t3LHlBCNumber =(Float_t)BCNumber;
+ t3LHlOrbitNumber =(Float_t)OrbitNumber;
+ t3LHlPeriodNumber =(Float_t)PeriodNumber;
+ t3LHleventtype =(Float_t)eventtype;
+ t3LHlpercentile =(Float_t)percentile;
+ t3LHlPx =(Float_t)tV0mom[0];
+ t3LHlPy =(Float_t)tV0mom[1];
+ t3LHlPz =(Float_t)tV0mom[2];
+ t3LHlEta =(Float_t)d->Eta();
+ t3LHlY =(Float_t)vHyp.Rapidity();
+ t3LHlM =(Float_t)vS.M();
+ t3LHlxPrimaryVertex =(Float_t)xPrimaryVertex;
+ t3LHlyPrimaryVertex =(Float_t)yPrimaryVertex;
+ t3LHlzPrimaryVertex =(Float_t)zPrimaryVertex;
+
+ t3LHlp0px =(Float_t)d->PxProng(0);
+ t3LHlp0py =(Float_t)d->PyProng(0);
+ t3LHlp0pz =(Float_t)d->PzProng(0);
+ t3LHlp0ch =(Float_t)d->ChargeProng(0);
+ t3LHlp1px =(Float_t)d->PxProng(1);
+ t3LHlp1py =(Float_t)d->PyProng(1);
+ t3LHlp1pz =(Float_t)d->PzProng(1);
+ t3LHlp1ch =(Float_t)d->ChargeProng(1);
+
+ fNtuple4->Fill();
+
+ }
+
+
+ //----------------------------------------
+ // Loop on tracks like-AnalysisTaskHelium3Pi.cxx
+
+ Int_t TrackNumber= DaughterTracks-> GetEntriesFast();
+
+ TArrayI Track0(TrackNumber); //Pions
+ Int_t nTrack0=0;
+
+ TArrayI Track1(TrackNumber); //Helium3
+ Int_t nTrack1=0;
+
+
+ Float_t impactXY=-999, impactZ=-999;
+ Float_t impactXYpi=-999, impactZpi=-999;
+
+ for(Int_t j=0; j<TrackNumber; j++){
+ // cout<<"eventtype inside: "<<eventtype<<endl;
+ // cout<<"Inside loop tracks"<<endl;
+ AliVTrack *track = (AliVTrack*)DaughterTracks->UncheckedAt(j);
+ AliAODTrack *aodtrack =(AliAODTrack*)track;
+
+ fhBB->Fill(aodtrack->GetDetPid()->GetTPCmomentum()*aodtrack->Charge(),aodtrack->GetTPCsignal());
+
+ //-----------------------------------------------------------
+ //Track cuts
+ // if(aodtrack->GetTPCNcls() < 80 )continue;
+ // if(aodtrack->Chi2perNDF() > 5 )continue;
+
+ // if (!aodtrack->IsOn(AliAODTrack::kTPCrefit)) continue;
+ // if (!aodtrack->IsOn(AliAODTrack::kTPCin)) continue;
+ // if (aodtrack->IsOn(AliAODTrack::kITSpureSA)) continue;
+
+ //---------------------------------------------------------------
+
+ if (aodtrack->GetTPCsignal()<100){
+ Track0[nTrack0++]=j;
+ fhBBPions->Fill(aodtrack->GetDetPid()->GetTPCmomentum()*aodtrack->Charge(),aodtrack->GetTPCsignal());
+ }
+
+ else{
+ Track1[nTrack1++]=j;
+ fhBBHe->Fill(aodtrack->GetDetPid()->GetTPCmomentum()*aodtrack->Charge(),aodtrack->GetTPCsignal());
+ }
+
+ }
+
+ Track0.Set(nTrack0);
+ Track1.Set(nTrack1);
+
+
+ //cout<<"nTrack0: "<<nTrack0<<endl;
+ //cout<<"nTrack1: "<<nTrack1<<endl;
+
+
+ //cout<<"Track loops..."<<endl;
+
+ Double_t DcaHeToPrimVertex=0;
+ Double_t DcaPionToPrimVertex=0;
+
+ impactXY=-999, impactZ=-999;
+ impactXYpi=-999, impactZpi=-999;
+
+
+ AliESDtrack *PionTrack = 0x0;
+ AliESDtrack *HeTrack = 0x0;
+
+
+ // Vettors for il PxPyPz
+
+ Double_t momPionVett[3];
+ for(Int_t i=0;i<3;i++)momPionVett[i]=0;
+
+ Double_t momHeVett[3];
+ for(Int_t i=0;i<3;i++)momHeVett[i]=0;
+
+ //At SV
+
+ Double_t momPionVettAt[3];
+ for(Int_t i=0;i<3;i++)momPionVettAt[i]=0;
+
+ Double_t momHeVettAt[3];
+ for(Int_t i=0;i<3;i++)momHeVettAt[i]=0;
+
+
+ //Loop on the tracks
+
+ for (Int_t k=0; k <nTrack0 ; k++) { //! Pion Tracks Loop
+
+ DcaPionToPrimVertex=0.;
+ DcaHeToPrimVertex=0;
+
+ Int_t Track0idx=Track0[k];
+
+ AliVTrack *trackpi = (AliVTrack*)DaughterTracks->UncheckedAt(Track0idx);
+ PionTrack = new AliESDtrack(trackpi);
+
+ statusPi = (ULong_t)trackpi->GetStatus();
+ IsPiITSRefit = ((statusPi) & (AliESDtrack::kITSrefit));
+
+ if (PionTrack)
+ DcaPionToPrimVertex = TMath::Abs(PionTrack->GetD(xPrimaryVertex, yPrimaryVertex,lMagneticField));
+
+ //cout<<"DcaPionToPrimVertex: "<<DcaPionToPrimVertex<<endl;
+
+ if(DcaPionToPrimVertex<0.2)continue;
+ AliExternalTrackParam trackInPion(*PionTrack);
+
+ for (Int_t i=0; i<nTrack1; i++){ //! He Tracks Loop
+
+ Int_t Track1idx=Track1[i];
+
+ AliVTrack *trackhe = (AliVTrack*)DaughterTracks->UncheckedAt(Track1idx);
+ HeTrack = new AliESDtrack(trackhe);
+
+ statusT= (ULong_t)HeTrack->GetStatus();
+ IsHeITSRefit = (statusT & AliESDtrack::kITSrefit);
+
+ if (HeTrack)
+ DcaHeToPrimVertex = TMath::Abs(HeTrack->GetD(xPrimaryVertex, yPrimaryVertex,lMagneticField)); //OK
+
+ AliExternalTrackParam trackInHe(*HeTrack);
+
+ if ( DcaPionToPrimVertex < fgDNmin) //OK
+ if ( DcaHeToPrimVertex < fgDNmin) continue; //OK
+
+ Double_t xn, xp;
+ Double_t dca=0.;
+
+ dca= PionTrack->GetDCA(HeTrack,lMagneticField,xn,xp); //!dca (Neg to Pos)
+
+ //cout<<"dca tracks: "<<dca<<endl;
+
+ if (dca > fgDCAmax) continue;
+ if ((xn+xp) > 2*fgRmax) continue;
+ if ((xn+xp) < 2*fgRmin) continue;
+
+ Bool_t corrected=kFALSE;
+ if ((trackInPion.GetX() > 3.) && (xn < 3.)) {
+ //correct for the beam pipe material
+ corrected=kTRUE;
+ }
+ if ((trackInHe.GetX() > 3.) && (xp < 3.)) {
+ //correct for the beam pipe material
+ corrected=kTRUE;
+ }
+ if (corrected) {
+ dca=trackInPion.GetDCA(&trackInHe,lMagneticField,xn,xp);
+ if (dca > fgDCAmax) continue;
+ if ((xn+xp) > 2*fgRmax) continue;
+ if ((xn+xp) < 2*fgRmin) continue;
+ }
+
+ //=============================================//
+ // Make "V0" with found tracks //
+ //=============================================//
+
+ trackInPion.PropagateTo(xn,lMagneticField);
+ trackInHe.PropagateTo(xp,lMagneticField);
+
+ AliESDv0 vertex(trackInPion,Track0idx,trackInHe,Track1idx);
+ // if (vertex.GetChi2V0() > fgChi2max) continue;
+
+ //cout<<"Made v0"<<endl;
+
+ Float_t CosPointingAngle=vertex.GetV0CosineOfPointingAngle(xPrimaryVertex,yPrimaryVertex,zPrimaryVertex); //PointingAngle
+ //cout<<"cpa: "<<CosPointingAngle<<endl;
+ // if (CosPointingAngle < fgCPAmin) continue;
+
+ vertex.SetDcaV0Daughters(dca);
+ vertex.SetV0CosineOfPointingAngle(CosPointingAngle);
+
+ fPos[0]=vertex.Xv();
+ fPos[1]=vertex.Yv();
+ fPos[2]=vertex.Zv();
+
+ HeTrack->PxPyPz(momHeVett);
+ PionTrack->PxPyPz(momPionVett);
+
+ Double_t raggio=TMath::Sqrt(fPos[0]*fPos[0]+fPos[1]*fPos[1]+fPos[2]*fPos[2]);
+ HeTrack->GetPxPyPzAt(raggio,lMagneticField,momHeVettAt);
+ PionTrack->GetPxPyPzAt(raggio,lMagneticField,momPionVettAt);
+
+ //------------------------------------------------------------------------//
+
+ HeTrack->GetImpactParameters(impactXY, impactZ);
+ PionTrack->GetImpactParameters(impactXYpi, impactZpi);
+
+ // if(vertex.GetD(xPrimaryVertex,yPrimaryVertex,zPrimaryVertex)>3) continue;
+
+ //salvo solo fino a 3.1 GeV/c2
+
+ vHelium.SetXYZM(2*momHeVettAt[0],2*momHeVettAt[1],2*momHeVettAt[2],Helium3Mass);
+ vPion.SetXYZM(momPionVettAt[0],momPionVettAt[1],momPionVettAt[2],PionMass);
+ vSum=vHelium+vPion;
+
+ // if(vSum.M()>3.2)
+ // continue;
+
+ Int_t fIdxInt[200]; //dummy array
+ Int_t nClustersTPCHe = HeTrack->GetTPCclusters(fIdxInt);
+ Int_t nClustersTPCPi = PionTrack->GetTPCclusters(fIdxInt);
+
+ //cout<<"Cos PA: "<<CosPointingAngle<<endl;
+ fHistEventMultiplicity->Fill(11.5);
+ //-------------------------------------------------------------
+
+ trunNumber =(Float_t)runNumber;
+ tbunchcross =(Float_t)BCNumber;
+ torbit =(Float_t)OrbitNumber;
+ tperiod =(Float_t)PeriodNumber;
+ teventtype =(Float_t)eventtype;
+ tTrackNumber =(Float_t)TrackNumber;
+ tpercentile =(Float_t)percentile;
+ txPrimaryVertex =(Float_t)xPrimaryVertex; //PRIMARY
+ tyPrimaryVertex =(Float_t)yPrimaryVertex;
+ tzPrimaryVertex =(Float_t)zPrimaryVertex;
+ txSecondaryVertex =(Float_t)fPos[0]; //SECONDARY
+ tySecondaryVertex =(Float_t)fPos[1];
+ tzSecondaryVertex =(Float_t)fPos[2];
+ tdcaTracks =(Float_t)dca; //between 2 tracks
+ tCosPointingAngle =(Float_t)CosPointingAngle; //cosPointingAngle da V0
+ tDCAV0toPrimaryVertex =(Float_t)vertex.GetD(xPrimaryVertex,yPrimaryVertex,zPrimaryVertex);
+ tHeSign =(Float_t)HeTrack->GetSign(); //He
+ tHepInTPC =(Float_t)trackInHe.GetP();
+ tHeTPCsignal =(Float_t)HeTrack->GetTPCsignal();
+ tDcaHeToPrimVertex =(Float_t)DcaHeToPrimVertex;
+ tHeEta =(Float_t)HeTrack->Eta();
+ tmomHex =(Float_t)momHeVett[0];
+ tmomHey =(Float_t)momHeVett[1];
+ tmomHez =(Float_t)momHeVett[2];
+ tmomHeAtSVx =(Float_t)momHeVettAt[0];
+ tmomHeAtSVy =(Float_t)momHeVettAt[1];
+ tmomHeAtSVz =(Float_t)momHeVettAt[2];
+ tHeTPCNcls =(Float_t)HeTrack->GetTPCNcls();
+ tHeimpactXY =(Float_t)impactXY;
+ tHeimpactZ =(Float_t)impactZ;
+ tHeITSClusterMap =(Float_t)HeTrack->GetITSClusterMap();
+ tIsHeITSRefit =(Float_t)IsHeITSRefit;
+ tPionSign =(Float_t)PionTrack->GetSign(); //Pion
+ tPionpInTPC =(Float_t)trackInPion.GetP();
+ tPionTPCsignal =(Float_t)PionTrack->GetTPCsignal();
+ tDcaPionToPrimVertex =(Float_t)DcaPionToPrimVertex;
+ tPionEta =(Float_t)PionTrack->Eta();
+ tmomPionx =(Float_t)momPionVett[0];
+ tmomPiony =(Float_t)momPionVett[1];
+ tmomPionz =(Float_t)momPionVett[2];
+ tmomNegPionAtSVx =(Float_t)momPionVettAt[0];
+ tmomNegPionAtSVy =(Float_t)momPionVettAt[1];
+ tmomNegPionAtSVz =(Float_t)momPionVettAt[2];
+ tPionTPCNcls =(Float_t)PionTrack->GetTPCNcls();
+ tPionimpactXY =(Float_t)impactXYpi;
+ tPionimpactZ =(Float_t)impactZpi;
+ tPionITSClusterMap =(Float_t)PionTrack->GetITSClusterMap();
+ tIsPiITSRefit =(Float_t)IsPiITSRefit;
+ txn =(Float_t)xn;
+ txp =(Float_t)xp;
+ tchi2He =(Float_t)HeTrack->GetTPCchi2()/(Float_t)(nClustersTPCHe);
+ tchi2Pi =(Float_t)PionTrack->GetTPCchi2()/(Float_t)(nClustersTPCPi);
+
+ fNtuple1->Fill();
+ vertex.Delete();
+
+ }
+
+ }
+
+ PostData(1,fListHist);
+ PostData(2,fNtuple1);
+ PostData(3,fNtuple4);
+
+} //end userexec
+
+
+//________________________________________________________________________
+
+void AliAnalysisTaskReadNuclexAOD::Terminate(Option_t *)
+{
+ // Draw result to the screen
+ // Called once at the end of the query
+}
+
+
--- /dev/null
+#ifndef ALIANALYSISTASKREASNUCLEXAOD_H
+#define ALIANALYSISTASKREASNUCLEXAOD_H
+
+/* See cxx source for full Copyright notice */
+
+//-----------------------------------------------------------------
+// AliAnalysisTaskReadNuclexAOD class
+//-----------------------------------------------------------------
+
+class TList;
+class TH1F;
+class TH2F;
+class TH3F;
+class TNtuple;
+class AliAODPid;
+class AliESDcascade;
+//class AliCascadeVertexer;
+
+#include "TString.h"
+
+#include "AliAnalysisTaskSE.h"
+
+class AliAnalysisTaskReadNuclexAOD : public AliAnalysisTaskSE {
+ public:
+ AliAnalysisTaskReadNuclexAOD();
+ AliAnalysisTaskReadNuclexAOD(const char *name);
+ virtual ~AliAnalysisTaskReadNuclexAOD();
+
+ virtual void UserCreateOutputObjects();
+ virtual void Init();
+ virtual void LocalInit() {Init();}
+ virtual void UserExec(Option_t *option);
+ virtual void Terminate(Option_t *);
+
+ void SetCollidingSystems(Short_t collidingSystems = 0) {fCollidingSystems = collidingSystems;}
+ void SetAnalysisType (const char* analysisType = "ESD") {fAnalysisType = analysisType;}
+ void SetDataType (const char* dataType = "REAL") {fDataType = dataType;}
+
+ //Double_t BetheBloch(Double_t bg,Double_t Charge,Bool_t optMC);
+ Double_t BetheBloch(Double_t bg,Double_t Charge,Bool_t isPbPb);
+
+
+ private:
+
+ TString fAnalysisType; //! "ESD" or "AOD" analysis type
+
+ Short_t fCollidingSystems; //! 0 = pp collisions or 1 = AA collisions
+ TString fDataType; //! "REAL" or "SIM" data type
+
+
+ TList *fListHist; //! List of histograms
+
+
+ TH1F *fHistEventMultiplicity;
+ TH2F *fHistTrackMultiplicity;
+ TH2F *fHistTrackMultiplicityCent;
+ TH2F *fHistTrackMultiplicitySemiCent;
+ TH2F *fHistTrackMultiplicityMB;
+ TH2F *fHistTrackMultiplicityPVCent;
+ TH2F *fHistTrackMultiplicityPVSemiCent;
+ TH2F *fHistTrackMultiplicityPVMB;
+ TH2F *fhBB;
+ TH2F *fhBBPions;
+ TH2F *fhBBHe;
+
+ TTree *aodTree;
+ TClonesArray *Nuclei;
+ TClonesArray *SecondaryVertices;
+ TClonesArray *DaughterTracks;
+
+ Int_t nevent; //event of the reduced tree
+
+ TTree *fNtuple1; //! Tree Pairs Pi/Proton "standard"
+
+ Float_t trunNumber;
+ Float_t tbunchcross;
+ Float_t torbit;
+ Float_t tperiod;
+ Float_t teventtype;
+ Float_t tTrackNumber;
+ Float_t tpercentile;
+ Float_t txPrimaryVertex;
+ Float_t tyPrimaryVertex;
+ Float_t tzPrimaryVertex;
+ Float_t txSecondaryVertex;
+ Float_t tySecondaryVertex;
+ Float_t tzSecondaryVertex;
+ Float_t tdcaTracks;
+ Float_t tCosPointingAngle;
+ Float_t tDCAV0toPrimaryVertex;
+ Float_t tHeSign;
+ Float_t tHepInTPC;
+ Float_t tHeTPCsignal;
+ Float_t tDcaHeToPrimVertex;
+ Float_t tHeEta;
+ Float_t tmomHex;
+ Float_t tmomHey;
+ Float_t tmomHez;
+ Float_t tmomHeAtSVx;
+ Float_t tmomHeAtSVy;
+ Float_t tmomHeAtSVz;
+ Float_t tHeTPCNcls;
+ Float_t tHeimpactXY;
+ Float_t tHeimpactZ;
+ Float_t tHeITSClusterMap;
+ Float_t tIsHeITSRefit;
+ Float_t tPionSign;
+ Float_t tPionpInTPC;
+ Float_t tPionTPCsignal;
+ Float_t tDcaPionToPrimVertex;
+ Float_t tPionEta;
+ Float_t tmomPionx;
+ Float_t tmomPiony;
+ Float_t tmomPionz;
+ Float_t tmomNegPionAtSVx;
+ Float_t tmomNegPionAtSVy;
+ Float_t tmomNegPionAtSVz;
+ Float_t tPionTPCNcls;
+ Float_t tPionimpactXY;
+ Float_t tPionimpactZ;
+ Float_t tPionITSClusterMap;
+ Float_t tIsPiITSRefit;
+ Float_t txn;
+ Float_t txp;
+ Float_t tchi2He;
+ Float_t tchi2Pi;
+
+
+ TTree *fNtuple4;
+
+ Float_t t3LHlrunNumber;
+ Float_t t3LHlBCNumber;
+ Float_t t3LHlOrbitNumber;
+ Float_t t3LHlPeriodNumber;
+ Float_t t3LHleventtype;
+ Float_t t3LHlpercentile;
+ Float_t t3LHlPx;
+ Float_t t3LHlPy;
+ Float_t t3LHlPz;
+ Float_t t3LHlEta;
+ Float_t t3LHlY;
+ Float_t t3LHlM;
+ Float_t t3LHlxPrimaryVertex;
+ Float_t t3LHlyPrimaryVertex;
+ Float_t t3LHlzPrimaryVertex;
+ Float_t t3LHlp0px;
+ Float_t t3LHlp0py;
+ Float_t t3LHlp0pz;
+ Float_t t3LHlp0ch;
+ Float_t t3LHlp1px;
+ Float_t t3LHlp1py;
+ Float_t t3LHlp1pz;
+ Float_t t3LHlp1ch;
+
+
+
+
+ static const Int_t fgNrot;
+
+
+ AliAnalysisTaskReadNuclexAOD(const AliAnalysisTaskReadNuclexAOD&); // not implemented
+ AliAnalysisTaskReadNuclexAOD& operator=(const AliAnalysisTaskReadNuclexAOD&); // not implemented
+
+ ClassDef(AliAnalysisTaskReadNuclexAOD, 0);
+};
+
+#endif
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff