+++ /dev/null
-/**************************************************************************
- * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- * *
- * Authors: Svein Lindal, Daniel Lohner *
- * Version 1.0 *
- * *
- * *
- * based on: on older version (see aliroot up to v5-04-42-AN) *
- * AliV0Reader.cxx *
- * Authors: Kathrin Koch, Kenneth Aamodt, Ana Marin *
- * *
- * 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. *
- **************************************************************************/
-
-////////////////////////////////////////////////
-//---------------------------------------------
-// Class reconstructing conversion photons from V0s
-//---------------------------------------------
-////////////////////////////////////////////////
-
-// The AliAODConversionPhotons will return the Position (== ESDTrack->GetID())of the positive (negative) ESDTrack
-// in the Array of ESDTracks stored in the ESDEvent via GetTrackLabelPositive() (GetTrackLabelNagative()).
-// Since it is the Position in the Array it is always positive.
-// After the conversion ESD --> AOD each AOD track will give you the former Position in the ESDArray via the
-// Function AODTrack->GetID(). AODTracks are stored for different TrackParameter (e.g. HybridTracks, TPC only ...). No Standard
-// AODTracks (not copies of AliExternalTrackParam) will be stored with negative values of AODTrack->GetID() (Formula: (ESDTrack->GetID()+1)*-1)
-// This leads to the possibility of more than one AOD track with the same negative ID. For that reason all AOD tracks are additionally flaged.
-// In this analysis we should therefore only use AODTracks with positive values of AODTrack->GetID().
-
-// If you want to find the AODTrack corresponding to the daugher track of a AliAODConversionPhoton you have to find the AODTrack with
-// AODTrack->GetID() == GetTrackLabelPositive() (GetTrackLabelNagative()).
-
-#include <TGeoGlobalMagField.h>
-
-#include "AliV0ReaderV1.h"
-#include "AliKFParticle.h"
-#include "AliAODv0.h"
-#include "AliESDv0.h"
-#include "AliAODEvent.h"
-#include "AliESDEvent.h"
-#include "AliPID.h"
-#include "AliMCEvent.h"
-#include "AliStack.h"
-#include "AliMCEventHandler.h"
-#include "AliESDpid.h"
-#include "AliESDtrackCuts.h"
-#include "TRandom3.h"
-#include "AliGenCocktailEventHeader.h"
-#include "TList.h"
-#include "AliKFConversionPhoton.h"
-#include "AliAODConversionPhoton.h"
-#include "AliConversionPhotonBase.h"
-#include "TVector.h"
-#include "AliKFVertex.h"
-#include "AliAODTrack.h"
-#include "AliESDtrack.h"
-#include "AliAnalysisManager.h"
-#include "AliInputEventHandler.h"
-#include "AliAODHandler.h"
-#include "AliPIDResponse.h"
-#include "TChain.h"
-#include "TFile.h"
-#include "TString.h"
-#include "TObjArray.h"
-
-class iostream;
-
-
-using namespace std;
-
-ClassImp(AliV0ReaderV1)
-
-//________________________________________________________________________
-AliV0ReaderV1::AliV0ReaderV1(const char *name) : AliAnalysisTaskSE(name),
- fConversionCuts(NULL),
- fEventCuts(NULL),
- fConversionGammas(NULL),
- fUseImprovedVertex(kTRUE),
- fUseOwnXYZCalculation(kTRUE),
- fUseConstructGamma(kFALSE),
- kUseAODConversionPhoton(kTRUE),
- fCreateAOD(kFALSE),
- fDeltaAODBranchName("GammaConv"),
- fDeltaAODFilename("AliAODGammaConversion.root"),
- fRelabelAODs(kFALSE),
- fEventIsSelected(kFALSE),
- fNumberOfPrimaryTracks(0),
- fPeriodName(""),
- fUseMassToZero(kTRUE)
-{
- // Default constructor
-
- DefineInput(0, TChain::Class());
-}
-
-//________________________________________________________________________
-AliV0ReaderV1::~AliV0ReaderV1()
-{
- // default deconstructor
-
- if(fConversionGammas){
- fConversionGammas->Delete();// Clear Objects
- delete fConversionGammas;
- fConversionGammas=0x0;
- }
-}
-
-/*
-//________________________________________________________________________
-AliV0ReaderV1::AliV0ReaderV1(AliV0ReaderV1 &original) : AliAnalysisTaskSE(original),
-fConversionCuts(NULL),
-fConversionGammas(NULL),
-fUseImprovedVertex(original.fUseImprovedVertex),
-fUseOwnXYZCalculation(original.fUseOwnXYZCalculation),
-fUseConstructGamma(original.fUseConstructGamma),
-kUseAODConversionPhoton(original.kUseAODConversionPhoton),
-fCreateAOD(original.fCreateAOD),
-fDeltaAODBranchName(original.fDeltaAODBranchName),
-fDeltaAODFilename(original.fDeltaAODFilename),
-fEventIsSelected(original.fEventIsSelected)
-{
-// Default constructor
-
-DefineInput(0, TChain::Class());
-}
-
-//____________________________________________________________
-AliV0ReaderV1 &AliV0ReaderV1::operator=(const AliV0ReaderV1 &ref){
-//
-// Assignment operator
-// Only copies pointers, object is not the owner of the references
-//
-if(this != &ref){
-AliAnalysisTaskSE::operator=(ref);
-fUseImprovedVertex=ref.fUseImprovedVertex;
-fUseOwnXYZCalculation=ref.fUseOwnXYZCalculation;
-fUseConstructGamma=ref.fUseConstructGamma;
-kUseAODConversionPhoton=ref.kUseAODConversionPhoton;
-fCreateAOD=ref.fCreateAOD;
-fDeltaAODBranchName=ref.fDeltaAODBranchName;
-fDeltaAODFilename=ref.fDeltaAODFilename;
-fEventIsSelected=ref.fEventIsSelected;
-}
-return *this;
-}
-*/
-
-//________________________________________________________________________
-void AliV0ReaderV1::Init()
-{
- // Initialize function to be called once before analysis
- if(fConversionCuts==NULL){
- if(fConversionCuts==NULL)AliError("No Conversion Cut Selection initialized");
- }
- if(fEventCuts==NULL){
- if(fEventCuts==NULL)AliError("No Event Cut Selection initialized");
- }
-
- if(fCreateAOD){kUseAODConversionPhoton=kTRUE;}
-
- if(fConversionGammas != NULL){
- delete fConversionGammas;
- fConversionGammas=NULL;
- }
-
- if(fConversionGammas == NULL){
- if(kUseAODConversionPhoton){
- fConversionGammas = new TClonesArray("AliAODConversionPhoton",100);}
- else{
- fConversionGammas = new TClonesArray("AliKFConversionPhoton",100);}
- }
- fConversionGammas->Delete();//Reset the TClonesArray
-}
-
-//________________________________________________________________________
-void AliV0ReaderV1::UserCreateOutputObjects()
-{
- // Create AODs
-
- if(fCreateAOD){
- if (fEventCuts){
- fDeltaAODBranchName.Append("_");
- fDeltaAODBranchName.Append(fEventCuts->GetCutNumber());
- }
- if(fConversionCuts){
- fDeltaAODBranchName.Append("_");
- fDeltaAODBranchName.Append(fConversionCuts->GetCutNumber());
- fDeltaAODBranchName.Append("_gamma");
- }
- fConversionGammas->SetName(fDeltaAODBranchName.Data());
-
- AddAODBranch("TClonesArray", &fConversionGammas, fDeltaAODFilename.Data());
- AliAnalysisManager::GetAnalysisManager()->RegisterExtraFile(fDeltaAODFilename.Data());
- }
-
-}
-//________________________________________________________________________
-Bool_t AliV0ReaderV1::Notify()
-{
- if (fPeriodName.CompareTo("") == 0){
- AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
- if(man) {
- AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
- if (inputHandler){
- TTree* tree = (TTree*) inputHandler->GetTree();
- TFile* file = (TFile*) tree->GetCurrentFile();
- TString fileName(file->GetName());
- TObjArray *arr = fileName.Tokenize("/");
- for (Int_t i = 0; i < arr->GetEntriesFast();i++ ){
- TObjString* testObjString = (TObjString*)arr->At(i);
- if (testObjString->GetString().Contains("LHC")){
- fPeriodName = testObjString->GetString();
- i = arr->GetEntriesFast();
- }
- }
- }
- }
- if(!fEventCuts->GetDoEtaShift()) return kTRUE; // No Eta Shift requested, continue
- if(fEventCuts->GetEtaShift() == 0.0){ // Eta Shift requested but not set, get shift automatically
- fEventCuts->GetCorrectEtaShiftFromPeriod(fPeriodName);
- fEventCuts->DoEtaShift(kFALSE); // Eta Shift Set, make sure that it is called only once
- return kTRUE;
- } else {
- printf(" Gamma Conversion Reader %s_%s :: Eta Shift Manually Set to %f \n\n",
- (fEventCuts->GetCutNumber()).Data(),(fConversionCuts->GetCutNumber()).Data(),fEventCuts->GetEtaShift());
- fEventCuts->DoEtaShift(kFALSE); // Eta Shift Set, make sure that it is called only once
- }
- }
- return kTRUE;
-}
-//________________________________________________________________________
-void AliV0ReaderV1::UserExec(Option_t *option){
-
- AliESDEvent * esdEvent = dynamic_cast<AliESDEvent*>(fInputEvent);
- if(esdEvent) {
- if (!TGeoGlobalMagField::Instance()->GetField()) esdEvent->InitMagneticField();
- }
-
- // Check if correctly initialized
- if(!fConversionGammas)Init();
-
- // User Exec
- fEventIsSelected=ProcessEvent(fInputEvent,fMCEvent);
-
- // AOD Output
- FillAODOutput();
-
-}
-
-//________________________________________________________________________
-Bool_t AliV0ReaderV1::ProcessEvent(AliVEvent *inputEvent,AliMCEvent *mcEvent)
-{
-
- //Reset the TClonesArray
- fConversionGammas->Delete();
-
- fInputEvent=inputEvent;
- fMCEvent=mcEvent;
-
- if(!fInputEvent){
- AliError("No Input event");
- return kFALSE;
- }
- if(!fEventCuts){AliError("No EventCuts");return kFALSE;}
- if(!fConversionCuts){AliError("No ConversionCuts");return kFALSE;}
-
- // Count Primary Tracks Event
- CountTracks();
-
- // Event Cuts
- if(!fEventCuts->EventIsSelected(fInputEvent,fMCEvent))return kFALSE;
-
- // Set Magnetic Field
- AliKFParticle::SetField(fInputEvent->GetMagneticField());
-
- if(fInputEvent->IsA()==AliESDEvent::Class()){
- ProcessESDV0s();
- }
- if(fInputEvent->IsA()==AliAODEvent::Class()){
- GetAODConversionGammas();
- }
-
- return kTRUE;
-}
-///________________________________________________________________________
-void AliV0ReaderV1::FillAODOutput()
-{
- // Fill AOD Output with reconstructed Photons
-
- if(fInputEvent->IsA()==AliESDEvent::Class()){
- ///Make sure delta aod is filled if standard aod is filled (for synchronization when reading aod with standard aod)
- if(fCreateAOD) {
- AliAODHandler * aodhandler = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
- if (aodhandler && aodhandler->GetFillAOD()) {
- AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillExtension(kTRUE);
- //PostData(0, fConversionGammas);
-
- }
- }
- }
-}
-
-///________________________________________________________________________
-const AliExternalTrackParam *AliV0ReaderV1::GetExternalTrackParam(AliESDv0 *fCurrentV0,Int_t &tracklabel,Int_t charge){
-
- // Get External Track Parameter with given charge
-
- if(!(charge==1||charge==-1)){AliError("Charge not defined");return 0x0;}
-
- // Check for sign flip
- if(fCurrentV0){
- if(!fCurrentV0->GetParamN()||!fCurrentV0->GetParamP())return 0x0;
- if(!fConversionCuts->GetTrack(fInputEvent,fCurrentV0->GetNindex())||!fConversionCuts->GetTrack(fInputEvent,fCurrentV0->GetPindex()))return 0x0;
- if((fConversionCuts->GetTrack(fInputEvent,fCurrentV0->GetPindex()))->Charge()==charge){
- tracklabel=fCurrentV0->GetPindex();
- return fCurrentV0->GetParamP();}
- if((fConversionCuts->GetTrack(fInputEvent,fCurrentV0->GetNindex()))->Charge()==charge){
- tracklabel=fCurrentV0->GetNindex();
- return fCurrentV0->GetParamN();}
- }
- return 0x0;
-}
-
-///________________________________________________________________________
-Bool_t AliV0ReaderV1::ProcessESDV0s()
-{
- // Process ESD V0s for conversion photon reconstruction
- AliESDEvent *fESDEvent=dynamic_cast<AliESDEvent*>(fInputEvent);
-
- AliKFConversionPhoton *fCurrentMotherKFCandidate=NULL;
-
- if(fESDEvent){
-
- for(Int_t currentV0Index=0;currentV0Index<fESDEvent->GetNumberOfV0s();currentV0Index++){
- AliESDv0 *fCurrentV0=(AliESDv0*)(fESDEvent->GetV0(currentV0Index));
- if(!fCurrentV0){
- printf("Requested V0 does not exist");
- continue;
- }
-
- fCurrentMotherKFCandidate=ReconstructV0(fCurrentV0,currentV0Index);
-
- if(fCurrentMotherKFCandidate){
- // Add Gamma to the TClonesArray
-
- if(kUseAODConversionPhoton){
- new((*fConversionGammas)[fConversionGammas->GetEntriesFast()]) AliAODConversionPhoton(fCurrentMotherKFCandidate);
- AliAODConversionPhoton * currentConversionPhoton = (AliAODConversionPhoton*)(fConversionGammas->At(fConversionGammas->GetEntriesFast()-1));
- currentConversionPhoton->SetMass(fCurrentMotherKFCandidate->M());
- if (fUseMassToZero) currentConversionPhoton->SetMassToZero();
- } else {
- new((*fConversionGammas)[fConversionGammas->GetEntriesFast()]) AliKFConversionPhoton(*fCurrentMotherKFCandidate);
- }
-
- delete fCurrentMotherKFCandidate;
- fCurrentMotherKFCandidate=NULL;
- }
- }
- }
- return kTRUE;
-}
-
-///________________________________________________________________________
-AliKFConversionPhoton *AliV0ReaderV1::ReconstructV0(AliESDv0 *fCurrentV0,Int_t currentV0Index)
-{
- // cout << currentV0Index << endl;
- // Reconstruct conversion photon from ESD v0
- fConversionCuts->FillPhotonCutIndex(AliConversionPhotonCuts::kPhotonIn);
-
- //checks if on the fly mode is set
- if(!fConversionCuts->SelectV0Finder(fCurrentV0->GetOnFlyStatus())){
- fConversionCuts->FillPhotonCutIndex(AliConversionPhotonCuts::kOnFly);
- return 0x0;
- }
-
- // TrackLabels
- Int_t currentTrackLabels[2]={-1,-1};
-
- // Get Daughter KF Particles
-
- const AliExternalTrackParam *fCurrentExternalTrackParamPositive=GetExternalTrackParamP(fCurrentV0,currentTrackLabels[0]);
- // cout << fCurrentExternalTrackParamPositive << "\t" << currentTrackLabels[0] << endl;
- const AliExternalTrackParam *fCurrentExternalTrackParamNegative=GetExternalTrackParamN(fCurrentV0,currentTrackLabels[1]);
- // cout << fCurrentExternalTrackParamNegative << "\t" << currentTrackLabels[1] << endl;
- if(!fCurrentExternalTrackParamPositive||!fCurrentExternalTrackParamNegative)return 0x0;
-
- // Apply some Cuts before Reconstruction
-
- AliVTrack * posTrack = fConversionCuts->GetTrack(fInputEvent,currentTrackLabels[0]);
- AliVTrack * negTrack = fConversionCuts->GetTrack(fInputEvent,currentTrackLabels[1]);
- if(!negTrack || !posTrack) {
- fConversionCuts->FillPhotonCutIndex(AliConversionPhotonCuts::kNoTracks);
- return 0x0;
- }
- // Track Cuts
- if(!fConversionCuts->TracksAreSelected(negTrack, posTrack)){
- fConversionCuts->FillPhotonCutIndex(AliConversionPhotonCuts::kTrackCuts);
- return 0x0;
- }
-
- fConversionCuts->FillV0EtaBeforedEdxCuts(fCurrentV0->Eta());
- if (!fConversionCuts->dEdxCuts(posTrack)) {
- fConversionCuts->FillPhotonCutIndex(AliConversionPhotonCuts::kdEdxCuts);
- return 0x0;
- }
- // PID Cuts
- if(!fConversionCuts->dEdxCuts(negTrack)) {
- fConversionCuts->FillPhotonCutIndex(AliConversionPhotonCuts::kdEdxCuts);
- return 0x0;
- }
- fConversionCuts->FillV0EtaAfterdEdxCuts(fCurrentV0->Eta());
- // Reconstruct Photon
- AliKFConversionPhoton *fCurrentMotherKF=NULL;
- // fUseConstructGamma = kFALSE;
- // cout << "construct gamma " << endl;
- AliKFParticle fCurrentNegativeKFParticle(*(fCurrentExternalTrackParamNegative),11);
- // cout << fCurrentExternalTrackParamNegative << "\t" << endl;
- AliKFParticle fCurrentPositiveKFParticle(*(fCurrentExternalTrackParamPositive),-11);
- // cout << fCurrentExternalTrackParamPositive << "\t" << endl;
- // cout << currentTrackLabels[0] << "\t" << currentTrackLabels[1] << endl;
- // cout << "construct gamma " <<fUseConstructGamma << endl;
-
- // Reconstruct Gamma
- if(fUseConstructGamma){
- fCurrentMotherKF = new AliKFConversionPhoton();
- fCurrentMotherKF->ConstructGamma(fCurrentNegativeKFParticle,fCurrentPositiveKFParticle);
- }else{
- fCurrentMotherKF = new AliKFConversionPhoton(fCurrentNegativeKFParticle,fCurrentPositiveKFParticle);
- fCurrentMotherKF->SetMassConstraint(0,0.0001);
- }
-
- // Set Track Labels
-
- fCurrentMotherKF->SetTrackLabels(currentTrackLabels[0],currentTrackLabels[1]);
-
- // Set V0 index
-
- fCurrentMotherKF->SetV0Index(currentV0Index);
-
- //Set MC Label
- if(fMCEvent){
-
- AliStack *fMCStack= fMCEvent->Stack();
-
- Int_t labelp=TMath::Abs(fConversionCuts->GetTrack(fInputEvent,fCurrentMotherKF->GetTrackLabelPositive())->GetLabel());
- Int_t labeln=TMath::Abs(fConversionCuts->GetTrack(fInputEvent,fCurrentMotherKF->GetTrackLabelNegative())->GetLabel());
-
- TParticle *fNegativeMCParticle = fMCStack->Particle(labeln);
- TParticle *fPositiveMCParticle = fMCStack->Particle(labelp);
-
- if(fPositiveMCParticle&&fNegativeMCParticle){
- fCurrentMotherKF->SetMCLabelPositive(labelp);
- fCurrentMotherKF->SetMCLabelNegative(labeln);
- }
- }
-
- // Update Vertex (moved for same eta compared to old)
- // cout << currentV0Index <<" \t before: \t" << fCurrentMotherKF->GetPx() << "\t" << fCurrentMotherKF->GetPy() << "\t" << fCurrentMotherKF->GetPz() << endl;
- if(fUseImprovedVertex == kTRUE){
- AliKFVertex primaryVertexImproved(*fInputEvent->GetPrimaryVertex());
- // cout << "Prim Vtx: " << primaryVertexImproved.GetX() << "\t" << primaryVertexImproved.GetY() << "\t" << primaryVertexImproved.GetZ() << endl;
- primaryVertexImproved+=*fCurrentMotherKF;
- fCurrentMotherKF->SetProductionVertex(primaryVertexImproved);
- }
- // SetPsiPair
-
- Double_t PsiPair=GetPsiPair(fCurrentV0,fCurrentExternalTrackParamPositive,fCurrentExternalTrackParamNegative);
- fCurrentMotherKF->SetPsiPair(PsiPair);
-
- // Recalculate ConversionPoint
- Double_t dca[2]={0,0};
- if(fUseOwnXYZCalculation){
- Double_t convpos[3]={0,0,0};
- if(!GetConversionPoint(fCurrentExternalTrackParamPositive,fCurrentExternalTrackParamNegative,convpos,dca)){
- fConversionCuts->FillPhotonCutIndex(AliConversionPhotonCuts::kConvPointFail);
- delete fCurrentMotherKF;
- fCurrentMotherKF=NULL;
- return 0x0;
- }
-
- fCurrentMotherKF->SetConversionPoint(convpos);
- }
-
- if(fCurrentMotherKF->GetNDF() > 0.)
- fCurrentMotherKF->SetChi2perNDF(fCurrentMotherKF->GetChi2()/fCurrentMotherKF->GetNDF()); //->Photon is created before all chi2 relevant changes are performed, set it "by hand"
-
-
- // Set Dilepton Mass (moved down for same eta compared to old)
- fCurrentMotherKF->SetMass(fCurrentMotherKF->M());
-
- // Apply Photon Cuts
-
- if(!fConversionCuts->PhotonCuts(fCurrentMotherKF,fInputEvent)){
- fConversionCuts->FillPhotonCutIndex(AliConversionPhotonCuts::kPhotonCuts);
- delete fCurrentMotherKF;
- fCurrentMotherKF=NULL;
- return 0x0;
- }
-
- // cout << currentV0Index <<" \t after: \t" <<fCurrentMotherKF->GetPx() << "\t" << fCurrentMotherKF->GetPy() << "\t" << fCurrentMotherKF->GetPz() << endl;
-
- fConversionCuts->FillPhotonCutIndex(AliConversionPhotonCuts::kPhotonOut);
- return fCurrentMotherKF;
-}
-
-///________________________________________________________________________
-Double_t AliV0ReaderV1::GetPsiPair(const AliESDv0* v0, const AliExternalTrackParam *positiveparam,const AliExternalTrackParam *negativeparam) const {
- //
- // Angle between daughter momentum plane and plane
- //
-
- AliExternalTrackParam nt(*negativeparam);
- AliExternalTrackParam pt(*positiveparam);
-
- Float_t magField = fInputEvent->GetMagneticField();
-
- Double_t xyz[3] = {0.,0.,0.};
- v0->GetXYZ(xyz[0],xyz[1],xyz[2]);
-
- // Double_t pPlus[3] = {pt.Px(),pt.Py(),pt.Pz()};
- // Double_t pMinus[3] = {nt.Px(),nt.Py(),nt.Pz()};
-
- // Double_t u[3] = {pPlus[0]+pMinus[0],pPlus[1]+pMinus[1],pPlus[2]+pMinus[2]};
- // Double_t normu = sqrt( (u[0]*u[0]) + (u[1]*u[1]) + (u[2]*u[2]) );
-
- // u[0] = u[0] / normu;
- // u[1] = u[1] / normu;
- // u[2] = u[2] / normu;
-
- // Double_t normpPlus = sqrt( (pPlus[0]*pPlus[0]) + (pPlus[1]*pPlus[1]) + (pPlus[2]*pPlus[2]) );
- // Double_t normpMinus = sqrt( (pMinus[0]*pMinus[0]) + (pMinus[1]*pMinus[1]) + (pMinus[2]*pMinus[2]) );
-
- // pPlus[0] = pPlus[0] / normpPlus;
- // pPlus[1] = pPlus[1] / normpPlus;
- // pPlus[2] = pPlus[2] / normpPlus;
-
- // pMinus[0] = pMinus[0] / normpMinus;
- // pMinus[1] = pMinus[1] / normpMinus;
- // pMinus[2] = pMinus[2] / normpMinus;
-
- // Double_t v[3] = {0,0,0}; // pPlus X pMinus
- // v[0] = (pPlus[1]*pMinus[2]) - (pPlus[2]*pMinus[1]);
- // v[1] = (pPlus[2]*pMinus[0]) - (pPlus[0]*pMinus[2]);
- // v[2] = (pPlus[0]*pMinus[1]) - (pPlus[1]*pMinus[0]);
-
- // Double_t w[3] = {0,0,0}; // u X v
- // w[0] = (u[1]*v[2]) - (u[2]*v[1]);
- // w[1] = (u[2]*v[0]) - (u[0]*v[2]);
- // w[2] = (u[0]*v[1]) - (u[1]*v[0]);
-
- // Double_t z[3] = {0,0,1};
- // Double_t wc[3] = {0,0,0}; // u X v
- // wc[0] = (u[1]*z[2]) - (u[2]*z[1]);
- // wc[1] = (u[2]*z[0]) - (u[0]*z[2]);
- // wc[2] = (u[0]*z[1]) - (u[1]*z[0]);
-
- // Double_t PhiV = TMath::ACos((w[0]*wc[0]) + (w[1]*wc[1]) + (w[2]*wc[2]));
- //return abs(PhiV);
-
-
- // TVector3 pPlus(pt.Px(),pt.Py(),pt.Pz());
- // TVector3 pMinus(nt.Px(),nt.Py(),nt.Pz());
-
- // TVector3 u = pMinus + pPlus;
- // u = u*(1/u.Mag());
-
- // TVector3 pHPlus = pPlus*(1/pPlus.Mag());
- // TVector3 pHMinus = pMinus*(1/pMinus.Mag());
-
- // TVector3 v = pHPlus.Cross(pHMinus);
- // TVector3 w = u.Cross(v);
- // TVector3 z(0,0,1);
- // TVector3 wc = u.Cross(z);
-
- // Double_t PhiV = w * wc;
-
- Double_t mn[3] = {0,0,0};
- Double_t mp[3] = {0,0,0};
-
- v0->GetNPxPyPz(mn[0],mn[1],mn[2]);//reconstructed cartesian momentum components of negative daughter;
- v0->GetPPxPyPz(mp[0],mp[1],mp[2]);//reconstructed cartesian momentum components of positive daughter;
-
- Double_t deltat = 1.;
- deltat = TMath::ATan(mp[2]/(TMath::Sqrt(mp[0]*mp[0] + mp[1]*mp[1])+1.e-13)) - TMath::ATan(mn[2]/(TMath::Sqrt(mn[0]*mn[0] + mn[1]*mn[1])+1.e-13));//difference of angles of the two daughter tracks with z-axis
- Double_t radiussum = TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]) + 50;//radius to which tracks shall be propagated
-
- Double_t momPosProp[3] = {0,0,0};
- Double_t momNegProp[3] = {0,0,0};
-
- Double_t psiPair = 4.;
- if(nt.PropagateTo(radiussum,magField) == 0) return psiPair; //propagate tracks to the outside -> Better Purity and Efficiency
- if(pt.PropagateTo(radiussum,magField) == 0) return psiPair; //propagate tracks to the outside -> Better Purity and Efficiency
-
- pt.GetPxPyPz(momPosProp);//Get momentum vectors of tracks after propagation
- nt.GetPxPyPz(momNegProp);
-
- Double_t pEle =
- TMath::Sqrt(momNegProp[0]*momNegProp[0]+momNegProp[1]*momNegProp[1]+momNegProp[2]*momNegProp[2]);//absolute momentum value of negative daughter
-
- Double_t pPos =
- TMath::Sqrt(momPosProp[0]*momPosProp[0]+momPosProp[1]*momPosProp[1]+momPosProp[2]*momPosProp[2]);//absolute momentum value of positive daughter
-
- Double_t scalarproduct =
- momPosProp[0]*momNegProp[0]+momPosProp[1]*momNegProp[1]+momPosProp[2]*momNegProp[2];//scalar product of propagated positive and negative daughters' momenta
-
- Double_t chipair = TMath::ACos(scalarproduct/(pEle*pPos));//Angle between propagated daughter tracks
-
- // psiPair = TMath::Abs(TMath::ASin(deltat/chipair));
- psiPair = TMath::ASin(deltat/chipair);
- return psiPair;
-}
-
-///________________________________________________________________________
-Bool_t AliV0ReaderV1::GetHelixCenter(const AliExternalTrackParam *track,Double_t center[2]){
-
- // Get Center of the helix track parametrization
-
- Int_t charge=track->Charge();
- Double_t b=fInputEvent->GetMagneticField();
-
- Double_t helix[6];
- track->GetHelixParameters(helix,b);
-
- Double_t xpos = helix[5];
- Double_t ypos = helix[0];
- Double_t radius = TMath::Abs(1./helix[4]);
- Double_t phi = helix[2];
-
- if(phi < 0){
- phi = phi + 2*TMath::Pi();
- }
-
- phi -= TMath::Pi()/2.;
- Double_t xpoint = radius * TMath::Cos(phi);
- Double_t ypoint = radius * TMath::Sin(phi);
-
- if(b<0){
- if(charge > 0){
- xpoint = - xpoint;
- ypoint = - ypoint;
- }
- }
- if(b>0){
- if(charge < 0){
- xpoint = - xpoint;
- ypoint = - ypoint;
- }
- }
- center[0] = xpos + xpoint;
- center[1] = ypos + ypoint;
-
- return 1;
-}
-///________________________________________________________________________
-Bool_t AliV0ReaderV1::GetConversionPoint(const AliExternalTrackParam *pparam,const AliExternalTrackParam *nparam,Double_t convpos[3],Double_t dca[2]){
-
- // Recalculate Conversion Point
-
- if(!pparam||!nparam)return kFALSE;
-
- Double_t helixcenterpos[2];
- GetHelixCenter(pparam,helixcenterpos);
-
- Double_t helixcenterneg[2];
- GetHelixCenter(nparam,helixcenterneg);
-
- Double_t helixpos[6];
- pparam->GetHelixParameters(helixpos,fInputEvent->GetMagneticField());
- Double_t posradius = TMath::Abs(1./helixpos[4]);
-
- Double_t helixneg[6];
- nparam->GetHelixParameters(helixneg,fInputEvent->GetMagneticField());
- Double_t negradius = TMath::Abs(1./helixneg[4]);
-
- // Calculate xy-position
-
- Double_t xpos = helixcenterpos[0];
- Double_t ypos = helixcenterpos[1];
- Double_t xneg = helixcenterneg[0];
- Double_t yneg = helixcenterneg[1];
-
- convpos[0] = (xpos*negradius + xneg*posradius)/(negradius+posradius);
- convpos[1] = (ypos*negradius+ yneg*posradius)/(negradius+posradius);
-
-
- // Calculate Track XY vertex position
-
- Double_t deltaXPos = convpos[0] - xpos;
- Double_t deltaYPos = convpos[1] - ypos;
-
- Double_t deltaXNeg = convpos[0] - xneg;
- Double_t deltaYNeg = convpos[1] - yneg;
-
- Double_t alphaPos = TMath::Pi() + TMath::ATan2(-deltaYPos,-deltaXPos);
- Double_t alphaNeg = TMath::Pi() + TMath::ATan2(-deltaYNeg,-deltaXNeg);
-
- Double_t vertexXNeg = xneg + TMath::Abs(negradius)*TMath::Cos(alphaNeg);
- Double_t vertexYNeg = yneg + TMath::Abs(negradius)*TMath::Sin(alphaNeg);
-
- Double_t vertexXPos = xpos + TMath::Abs(posradius)*TMath::Cos(alphaPos);
- Double_t vertexYPos = ypos + TMath::Abs(posradius)*TMath::Sin(alphaPos);
-
- AliExternalTrackParam p(*pparam);
- AliExternalTrackParam n(*nparam);
-
- TVector2 vertexPos(vertexXPos,vertexYPos);
- TVector2 vertexNeg(vertexXNeg,vertexYNeg);
-
- // Convert to local coordinate system
- TVector2 vertexPosRot=vertexPos.Rotate(-p.GetAlpha());
- TVector2 vertexNegRot=vertexNeg.Rotate(-n.GetAlpha());
-
- // Propagate Track Params to Vertex
-
- if(!p.PropagateTo(vertexPosRot.X(),fInputEvent->GetMagneticField()))return kFALSE;
- if(!n.PropagateTo(vertexNegRot.X(),fInputEvent->GetMagneticField()))return kFALSE;
-
- // Check whether propagation was sucessful
-
- if(TMath::Abs(vertexPos.Mod()-TMath::Sqrt(p.GetX()*p.GetX()+p.GetY()*p.GetY()))>0.01)return kFALSE;
- if(TMath::Abs(vertexNeg.Mod()-TMath::Sqrt(n.GetX()*n.GetX()+n.GetY()*n.GetY()))>0.01)return kFALSE;
-
- // Calculate z position
-
- convpos[2] = (p.GetZ()*negradius+n.GetZ()*posradius)/(negradius+posradius);
-
- // Calculate DCA
- TVector2 vdca=vertexPos-vertexNeg;
- dca[0]=vdca.Mod();
- dca[1]=TMath::Abs(n.GetZ()-p.GetZ());
-
- return kTRUE;
- }
- //________________________________________________________________________
- Bool_t AliV0ReaderV1::GetAODConversionGammas(){
-
- // Get reconstructed conversion photons from satellite AOD file
-
- AliAODEvent *fAODEvent=dynamic_cast<AliAODEvent*>(fInputEvent);
-
- if(fAODEvent){
-
- if(fConversionGammas == NULL){
- fConversionGammas = new TClonesArray("AliAODConversionPhoton",100);
- }
- fConversionGammas->Delete();//Reset the TClonesArray
-
- //Get Gammas from satellite AOD gamma branch
-
- AliAODConversionPhoton *gamma=0x0;
-
- TClonesArray *fInputGammas=dynamic_cast<TClonesArray*>(fAODEvent->FindListObject(fDeltaAODBranchName.Data()));
-
- if(!fInputGammas){
- FindDeltaAODBranchName();
- fInputGammas=dynamic_cast<TClonesArray*>(fAODEvent->FindListObject(fDeltaAODBranchName.Data()));}
- if(!fInputGammas){AliError("No Gamma Satellites found");return kFALSE;}
- // Apply Selection Cuts to Gammas and create local working copy
- if(fInputGammas){
- for(Int_t i=0;i<fInputGammas->GetEntriesFast();i++){
- gamma=dynamic_cast<AliAODConversionPhoton*>(fInputGammas->At(i));
- if(gamma){
- if(fRelabelAODs)RelabelAODPhotonCandidates(gamma);
- if(fConversionCuts->PhotonIsSelected(gamma,fInputEvent)){
- new((*fConversionGammas)[fConversionGammas->GetEntriesFast()]) AliAODConversionPhoton(*gamma);}
- }
- }
- }
- }
-
- if(fConversionGammas->GetEntries()){return kTRUE;}
-
- return kFALSE;
-}
-
-//________________________________________________________________________
-void AliV0ReaderV1::FindDeltaAODBranchName(){
-
- // Find delta AOD branchname containing reconstructed photons
-
- TList *list=fInputEvent->GetList();
- for(Int_t ii=0;ii<list->GetEntries();ii++){
- TString name((list->At(ii))->GetName());
- if(name.BeginsWith(fDeltaAODBranchName)&&name.EndsWith("gamma")){
- fDeltaAODBranchName=name;
- AliInfo(Form("Set DeltaAOD BranchName to: %s",fDeltaAODBranchName.Data()));
- return;
- }
- }
-}
-//________________________________________________________________________
-void AliV0ReaderV1::RelabelAODPhotonCandidates(AliAODConversionPhoton *PhotonCandidate){
-
- // Relabeling For AOD Event
- // ESDiD -> AODiD
- // MCLabel -> AODMCLabel
- Bool_t AODLabelPos = kFALSE;
- Bool_t AODLabelNeg = kFALSE;
-
- for(Int_t i = 0; i<fInputEvent->GetNumberOfTracks();i++){
- AliAODTrack *tempDaughter = static_cast<AliAODTrack*>(fInputEvent->GetTrack(i));
- if(!AODLabelPos){
- if( tempDaughter->GetID() == PhotonCandidate->GetTrackLabelPositive() ){
- PhotonCandidate->SetMCLabelPositive(abs(tempDaughter->GetLabel()));
- PhotonCandidate->SetLabelPositive(i);
- AODLabelPos = kTRUE;
- }
- }
- if(!AODLabelNeg){
- if( tempDaughter->GetID() == PhotonCandidate->GetTrackLabelNegative()){
- PhotonCandidate->SetMCLabelNegative(abs(tempDaughter->GetLabel()));
- PhotonCandidate->SetLabelNegative(i);
- AODLabelNeg = kTRUE;
- }
- }
- if(AODLabelNeg && AODLabelPos){
- return;
- }
- }
- if(!AODLabelPos || !AODLabelNeg){
- AliError(Form("NO AOD Daughters Found Pos: %i %i Neg: %i %i",AODLabelPos,PhotonCandidate->GetTrackLabelPositive(),AODLabelNeg,PhotonCandidate->GetTrackLabelNegative()));
- }
-
-}
-
-//************************************************************************
-// This function counts the number of primary tracks in the event, the
-// implementation for ESD and AOD had to be different due to the different
-// filters which are already applied on AOD level the so-called filter
-// bits, we tried to replicate the conditions in both but an exact match
-// could not be reached. The cuts are similar to the ones used by the
-// jet-group
-//************************************************************************
-//________________________________________________________________________
-void AliV0ReaderV1::CountTracks(){
-
- // In Case of MC count only MB tracks
- // if(fMCEvent && fInputEvent->IsA()==AliESDEvent::Class()){
- // fEventCuts->GetNotRejectedParticles(1,NULL,fMCEvent);
- // }
- // if(fMCEvent && fInputEvent->IsA()==AliAODEvent::Class()){
- // fEventCuts->GetNotRejectedParticles(1,NULL,fInputEvent);
- // }
-
- if(fInputEvent->IsA()==AliESDEvent::Class()){
- // Using standard function for setting Cuts
- Bool_t selectPrimaries=kTRUE;
- AliESDtrackCuts *EsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(selectPrimaries);
- EsdTrackCuts->SetMaxDCAToVertexZ(2);
- EsdTrackCuts->SetEtaRange(-0.8, 0.8);
- EsdTrackCuts->SetPtRange(0.15);
- fNumberOfPrimaryTracks = 0;
- for(Int_t iTracks = 0; iTracks < fInputEvent->GetNumberOfTracks(); iTracks++){
- AliESDtrack* curTrack = (AliESDtrack*) fInputEvent->GetTrack(iTracks);
- if(!curTrack) continue;
- if(!EsdTrackCuts->AcceptTrack(curTrack)) continue;
- //if(fMCEvent && !(fEventCuts->IsParticleFromBGEvent(abs(curTrack->GetLabel()),fMCEvent->Stack(),fInputEvent))) continue;
- fNumberOfPrimaryTracks++;
- }
- delete EsdTrackCuts;
- EsdTrackCuts=0x0;
- }
- else if(fInputEvent->IsA()==AliAODEvent::Class()){
- fNumberOfPrimaryTracks = 0;
- for(Int_t iTracks = 0; iTracks<fInputEvent->GetNumberOfTracks(); iTracks++){
- AliAODTrack* curTrack = (AliAODTrack*) fInputEvent->GetTrack(iTracks);
- if(curTrack->GetID()<0) continue; // Avoid double counting of tracks
- if(!curTrack->IsHybridGlobalConstrainedGlobal()) continue;
- if(abs(curTrack->Eta())>0.8) continue;
- if(curTrack->Pt()<0.15) continue;
- //if(fMCEvent && !(fEventCuts->IsParticleFromBGEvent(abs(curTrack->GetLabel()),NULL,fInputEvent))) continue;
- //if(abs(curTrack->ZAtDCA())>2) continue; // Only Set For TPCOnly tracks
- fNumberOfPrimaryTracks++;
- }
- }
-
- return;
-}
-
-//________________________________________________________________________
-void AliV0ReaderV1::Terminate(Option_t *)
-{
-
-}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff