1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Authors: Svein Lindal, Daniel Lohner *
8 * based on: on older version (see aliroot up to v5-04-42-AN) *
10 * Authors: Kathrin Koch, Kenneth Aamodt, Ana Marin *
12 * Permission to use, copy, modify and distribute this software and its *
13 * documentation strictly for non-commercial purposes is hereby granted *
14 * without fee, provided that the above copyright notice appears in all *
15 * copies and that both the copyright notice and this permission notice *
16 * appear in the supporting documentation. The authors make no claims *
17 * about the suitability of this software for any purpose. It is *
18 * provided "as is" without express or implied warranty. *
19 **************************************************************************/
21 ////////////////////////////////////////////////
22 //---------------------------------------------
23 // Class reconstructing conversion photons from V0s
24 //---------------------------------------------
25 ////////////////////////////////////////////////
27 #include "AliV0ReaderV1.h"
28 #include "AliKFParticle.h"
31 #include "AliAODEvent.h"
32 #include "AliESDEvent.h"
34 #include "AliMCEvent.h"
36 #include "AliMCEventHandler.h"
37 #include "AliESDpid.h"
38 #include "AliESDtrackCuts.h"
40 #include "AliGenCocktailEventHeader.h"
42 #include "AliKFConversionPhoton.h"
43 #include "AliAODConversionPhoton.h"
44 #include "AliConversionPhotonBase.h"
46 #include "AliKFVertex.h"
47 #include "AliAODTrack.h"
48 #include "AliESDtrack.h"
49 #include "AliAnalysisManager.h"
50 #include "AliInputEventHandler.h"
51 #include "AliAODHandler.h"
52 #include "AliPIDResponse.h"
56 #include "TObjArray.h"
63 ClassImp(AliV0ReaderV1)
65 //________________________________________________________________________
66 AliV0ReaderV1::AliV0ReaderV1(const char *name) : AliAnalysisTaskSE(name),
67 fConversionCuts(NULL),
68 fConversionGammas(NULL),
69 fUseImprovedVertex(kTRUE),
70 fUseOwnXYZCalculation(kTRUE),
71 fUseConstructGamma(kFALSE),
72 kUseAODConversionPhoton(kTRUE),
74 fDeltaAODBranchName("GammaConv"),
75 fDeltaAODFilename("AliAODGammaConversion.root"),
76 fEventIsSelected(kFALSE),
79 // Default constructor
81 DefineInput(0, TChain::Class());
84 //________________________________________________________________________
85 AliV0ReaderV1::~AliV0ReaderV1()
87 // default deconstructor
89 if(fConversionGammas){
90 fConversionGammas->Delete();// Clear Objects
91 delete fConversionGammas;
92 fConversionGammas=0x0;
96 //________________________________________________________________________
97 AliV0ReaderV1::AliV0ReaderV1(AliV0ReaderV1 &original) : AliAnalysisTaskSE(original),
98 fConversionCuts(NULL),
99 fConversionGammas(NULL),
100 fUseImprovedVertex(original.fUseImprovedVertex),
101 fUseOwnXYZCalculation(original.fUseOwnXYZCalculation),
102 fUseConstructGamma(original.fUseConstructGamma),
103 kUseAODConversionPhoton(original.kUseAODConversionPhoton),
104 fCreateAOD(original.fCreateAOD),
105 fDeltaAODBranchName(original.fDeltaAODBranchName),
106 fDeltaAODFilename(original.fDeltaAODFilename),
107 fEventIsSelected(original.fEventIsSelected)
109 // Default constructor
111 DefineInput(0, TChain::Class());
114 //____________________________________________________________
115 AliV0ReaderV1 &AliV0ReaderV1::operator=(const AliV0ReaderV1 &ref){
117 // Assignment operator
118 // Only copies pointers, object is not the owner of the references
121 AliAnalysisTaskSE::operator=(ref);
122 fUseImprovedVertex=ref.fUseImprovedVertex;
123 fUseOwnXYZCalculation=ref.fUseOwnXYZCalculation;
124 fUseConstructGamma=ref.fUseConstructGamma;
125 kUseAODConversionPhoton=ref.kUseAODConversionPhoton;
126 fCreateAOD=ref.fCreateAOD;
127 fDeltaAODBranchName=ref.fDeltaAODBranchName;
128 fDeltaAODFilename=ref.fDeltaAODFilename;
129 fEventIsSelected=ref.fEventIsSelected;
134 //________________________________________________________________________
135 void AliV0ReaderV1::Init()
137 // Initialize function to be called once before analysis
139 if(fConversionCuts==NULL){
140 if(fConversionCuts==NULL)AliError("No Cut Selection initialized");
143 if(fCreateAOD){kUseAODConversionPhoton=kTRUE;}
145 if(fConversionGammas != NULL){
146 delete fConversionGammas;
147 fConversionGammas=NULL;
150 if(fConversionGammas == NULL){
151 if(kUseAODConversionPhoton){
152 fConversionGammas = new TClonesArray("AliAODConversionPhoton",100);}
154 fConversionGammas = new TClonesArray("AliKFConversionPhoton",100);}
156 fConversionGammas->Delete();//Reset the TClonesArray
159 //________________________________________________________________________
160 void AliV0ReaderV1::UserCreateOutputObjects()
166 fDeltaAODBranchName.Append("_");
167 fDeltaAODBranchName.Append(fConversionCuts->GetCutNumber());
168 fDeltaAODBranchName.Append("_gamma");
170 fConversionGammas->SetName(fDeltaAODBranchName.Data());
172 AddAODBranch("TClonesArray", &fConversionGammas, fDeltaAODFilename.Data());
173 AliAnalysisManager::GetAnalysisManager()->RegisterExtraFile(fDeltaAODFilename.Data());
177 //________________________________________________________________________
178 Bool_t AliV0ReaderV1::Notify()
180 if (fPeriodName.CompareTo("") == 0){
181 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
183 AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
185 TTree* tree = (TTree*) inputHandler->GetTree();
186 TFile* file = (TFile*) tree->GetCurrentFile();
187 TString fileName(file->GetName());
188 TObjArray *arr = fileName.Tokenize("/");
189 for (Int_t i = 0; i < arr->GetEntriesFast();i++ ){
190 TObjString* testObjString = (TObjString*)arr->At(i);
191 if (testObjString->GetString().Contains("LHC")){
192 fPeriodName = testObjString->GetString();
193 i = arr->GetEntriesFast();
199 if(fConversionCuts->GetDoEtaShift()){
200 if(fPeriodName.CompareTo("LHC12g") == 0 || //pilot run 2012
201 fPeriodName.CompareTo("LHC13b") == 0 || //mainly minimum bias
202 fPeriodName.CompareTo("LHC13c") == 0 || //mainly minimum bias
203 fPeriodName.CompareTo("LHC13d") == 0 || //mainly triggered
204 fPeriodName.CompareTo("LHC13e") == 0 || //mainly triggered
205 fPeriodName.CompareTo("LHC13c3") == 0 || //MC Starlight, anchor LHC13d+e
206 fPeriodName.CompareTo("LHC13c2") == 0 || //MC Starlight, coherent J/Psi, UPC muon anchor LHC13d+e
207 fPeriodName.CompareTo("LHC13b4") == 0 || //MC Pythia 6 (Jet-Jet), anchor LHC13b
208 fPeriodName.CompareTo("LHC13b2_fix_1") == 0 || //MC DPMJET, anchr LHC13b+c
209 fPeriodName.CompareTo("LHC13b3") == 0 || //MC HIJING, weighted to number of events per run, anchor LHC13b
210 fPeriodName.CompareTo("LHC13b2") == 0 || // MC DPMJET, wrong energy, anchor LHC13b
211 fPeriodName.CompareTo("LHC13b2_plus") == 0 || // MC DPMJET, weighted to number event per run, anchor LHC13b
212 fPeriodName.CompareTo("LHC13c1_bis") == 0 || // MC AMPT fast generation, pT hardbin, anchor ?
213 fPeriodName.CompareTo("LHC13c1") == 0 || // MC AMPT fast generation, anchor ?
214 fPeriodName.CompareTo("LHC13b1") == 0 || // MC DPMJET, fragments, with fixed label 0, anchor LHC12g
215 fPeriodName.CompareTo("LHC12g4b_fix") == 0 || // MC DPMJET, with fixed label 0, anchor LHC12g
216 fPeriodName.CompareTo("LHC12g1_fix") == 0 || // MC ?, with fixed label 0, anchor LHC12g
217 fPeriodName.CompareTo("LHC12g4c") == 0 || // MC DPMJET, shifted vertex runs, anchor LHC12g
218 fPeriodName.CompareTo("LHC12h6") == 0 || // MC muon cocktail, anchor LHC12g
219 fPeriodName.CompareTo("LHC12g4b") == 0 || // MC DPMJET 3rd iteration, anchor LHC12g
220 fPeriodName.CompareTo("LHC12g4a") == 0 || // MC DPMJET improved, anchor LHC12g
221 fPeriodName.CompareTo("LHC12g4") == 0 || // MC DPMJET, anchor LHC12g
222 fPeriodName.CompareTo("LHC12g5") == 0 || // MC PHOJET, anchor LHC12g
223 fPeriodName.CompareTo("LHC12g2") == 0 || // MC Starlight background, anchor LHC12g
224 fPeriodName.CompareTo("LHC12g1") == 0 // MC ?, anchor LHC12g
226 cout<<"AliV0ReaderV1 --> pPb Run!!! Eta Shift of "<<-0.465<<endl;
227 fConversionCuts->SetEtaShift(-0.465);
229 else if(fPeriodName.CompareTo("LHC13f") == 0 ||
230 fPeriodName.CompareTo("LHC13c6b") == 0 ||// MC Jpsi -> mumu, anchor LHC13f
231 fPeriodName.CompareTo("LHC13c5") == 0 || //MC Starlight, gamma gamma UPC muon, anchor LHC13f
232 fPeriodName.CompareTo("LHC13c4") == 0 //MC Starlight, coherent JPsi, UPC muon, anchor LHC13f
234 cout<<"AliV0ReaderV1 --> Pbp Run!!! Eta Shift of +"<<0.465<<endl;
235 fConversionCuts->SetEtaShift(+0.465);
237 else cout<<"AliV0ReaderV1 --> Eta Shift Requested but Period not Known"<<endl;
239 if(fConversionCuts->IsEtaShiftForced() == 1){
240 cout<<"AliV0ReaderV1 --> Force Eta Shift for non Pbp or pPb Run!!! Eta Shift of "<<-0.465<<endl;
241 fConversionCuts->SetEtaShift(-0.465);
243 else if(fConversionCuts->IsEtaShiftForced() == 2){
244 cout<<"AliV0ReaderV1 --> Force Eta Shift for non Pbp or pPb Run!!! Eta Shift of +"<<0.465<<endl;
245 fConversionCuts->SetEtaShift(0.465);
251 //________________________________________________________________________
252 void AliV0ReaderV1::UserExec(Option_t *option){
255 // Check if correctly initialized
256 if(!fConversionGammas)Init();
259 fEventIsSelected=ProcessEvent(fInputEvent,fMCEvent);
262 //________________________________________________________________________
263 Bool_t AliV0ReaderV1::ProcessEvent(AliVEvent *inputEvent,AliMCEvent *mcEvent)
265 //Reset the TClonesArray
266 fConversionGammas->Delete();
268 fInputEvent=inputEvent;
272 AliError("No Input event");
276 if(!fConversionCuts){AliError("No ConversionCuts");return kFALSE;}
279 if(!fConversionCuts->EventIsSelected(fInputEvent,fMCEvent))return kFALSE;
281 // Set Magnetic Field
282 AliKFParticle::SetField(fInputEvent->GetMagneticField());
284 if(fInputEvent->IsA()==AliESDEvent::Class()){
287 if(fInputEvent->IsA()==AliAODEvent::Class()){
288 GetAODConversionGammas();
296 ///________________________________________________________________________
297 void AliV0ReaderV1::FillAODOutput()
299 // Fill AOD Output with reconstructed Photons
301 if(fInputEvent->IsA()==AliESDEvent::Class()){
302 ///Make sure delta aod is filled if standard aod is filled (for synchronization when reading aod with standard aod)
304 AliAODHandler * aodhandler = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
305 if (aodhandler && aodhandler->GetFillAOD()) {
306 AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillExtension(kTRUE);
307 //PostData(0, fConversionGammas);
314 ///________________________________________________________________________
315 const AliExternalTrackParam *AliV0ReaderV1::GetExternalTrackParam(AliESDv0 *fCurrentV0,Int_t &tracklabel,Int_t charge){
317 // Get External Track Parameter with given charge
319 if(!(charge==1||charge==-1)){AliError("Charge not defined");return 0x0;}
321 // Check for sign flip
323 if(!fCurrentV0->GetParamN()||!fCurrentV0->GetParamP())return 0x0;
324 if(!fConversionCuts->GetTrack(fInputEvent,fCurrentV0->GetNindex())||!fConversionCuts->GetTrack(fInputEvent,fCurrentV0->GetPindex()))return 0x0;
325 if((fConversionCuts->GetTrack(fInputEvent,fCurrentV0->GetPindex()))->Charge()==charge){
326 tracklabel=fCurrentV0->GetPindex();
327 return fCurrentV0->GetParamP();}
328 if((fConversionCuts->GetTrack(fInputEvent,fCurrentV0->GetNindex()))->Charge()==charge){
329 tracklabel=fCurrentV0->GetNindex();
330 return fCurrentV0->GetParamN();}
335 ///________________________________________________________________________
336 Bool_t AliV0ReaderV1::ProcessESDV0s()
338 // Process ESD V0s for conversion photon reconstruction
340 AliESDEvent *fESDEvent=dynamic_cast<AliESDEvent*>(fInputEvent);
342 AliKFConversionPhoton *fCurrentMotherKFCandidate=NULL;
346 for(Int_t currentV0Index=0;currentV0Index<fESDEvent->GetNumberOfV0s();currentV0Index++){
347 AliESDv0 *fCurrentV0=(AliESDv0*)(fESDEvent->GetV0(currentV0Index));
349 printf("Requested V0 does not exist");
352 fCurrentMotherKFCandidate=ReconstructV0(fCurrentV0,currentV0Index);
354 if(fCurrentMotherKFCandidate){
356 // Add Gamma to the TClonesArray
358 if(kUseAODConversionPhoton){
359 new((*fConversionGammas)[fConversionGammas->GetEntriesFast()]) AliAODConversionPhoton(fCurrentMotherKFCandidate);
362 new((*fConversionGammas)[fConversionGammas->GetEntriesFast()]) AliKFConversionPhoton(*fCurrentMotherKFCandidate);
365 delete fCurrentMotherKFCandidate;
366 fCurrentMotherKFCandidate=NULL;
373 ///________________________________________________________________________
374 AliKFConversionPhoton *AliV0ReaderV1::ReconstructV0(AliESDv0 *fCurrentV0,Int_t currentV0Index)
376 // Reconstruct conversion photon from ESD v0
377 fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kPhotonIn);
379 //checks if on the fly mode is set
380 if(!fConversionCuts->SelectV0Finder(fCurrentV0->GetOnFlyStatus())){
381 fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kOnFly);
386 Int_t currentTrackLabels[2]={-1,-1};
388 // Get Daughter KF Particles
390 const AliExternalTrackParam *fCurrentExternalTrackParamPositive=GetExternalTrackParamP(fCurrentV0,currentTrackLabels[0]);
391 const AliExternalTrackParam *fCurrentExternalTrackParamNegative=GetExternalTrackParamN(fCurrentV0,currentTrackLabels[1]);
393 if(!fCurrentExternalTrackParamPositive||!fCurrentExternalTrackParamNegative)return 0x0;
395 // Apply some Cuts before Reconstruction
397 AliVTrack * posTrack = fConversionCuts->GetTrack(fInputEvent,currentTrackLabels[0]);
398 AliVTrack * negTrack = fConversionCuts->GetTrack(fInputEvent,currentTrackLabels[1]);
400 if(!negTrack || !posTrack) {
401 fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kNoTracks);
406 if(!fConversionCuts->TracksAreSelected(negTrack, posTrack)){
407 fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kTrackCuts);
412 if(!fConversionCuts->dEdxCuts(negTrack) || !fConversionCuts->dEdxCuts(posTrack)) {
413 fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kdEdxCuts);
417 // Reconstruct Photon
419 AliKFConversionPhoton *fCurrentMotherKF=NULL;
421 AliKFParticle fCurrentNegativeKFParticle(*(fCurrentExternalTrackParamNegative),11);
422 AliKFParticle fCurrentPositiveKFParticle(*(fCurrentExternalTrackParamPositive),-11);
426 if(fUseConstructGamma){
427 fCurrentMotherKF = new AliKFConversionPhoton();
428 fCurrentMotherKF->ConstructGamma(fCurrentNegativeKFParticle,fCurrentPositiveKFParticle);
430 fCurrentMotherKF = new AliKFConversionPhoton(fCurrentNegativeKFParticle,fCurrentPositiveKFParticle);
431 fCurrentMotherKF->SetMassConstraint(0,0);
436 fCurrentMotherKF->SetTrackLabels(currentTrackLabels[0],currentTrackLabels[1]);
440 fCurrentMotherKF->SetV0Index(currentV0Index);
446 AliStack *fMCStack= fMCEvent->Stack();
448 Int_t labelp=TMath::Abs(fConversionCuts->GetTrack(fInputEvent,fCurrentMotherKF->GetTrackLabelPositive())->GetLabel());
449 Int_t labeln=TMath::Abs(fConversionCuts->GetTrack(fInputEvent,fCurrentMotherKF->GetTrackLabelNegative())->GetLabel());
451 TParticle *fNegativeMCParticle = fMCStack->Particle(labeln);
452 TParticle *fPositiveMCParticle = fMCStack->Particle(labelp);
454 if(fPositiveMCParticle&&fNegativeMCParticle){
455 fCurrentMotherKF->SetMCLabelPositive(labelp);
456 fCurrentMotherKF->SetMCLabelNegative(labeln);
460 // Update Vertex (moved for same eta compared to old)
461 if(fUseImprovedVertex == kTRUE){
462 AliKFVertex primaryVertexImproved(*fInputEvent->GetPrimaryVertex());
463 primaryVertexImproved+=*fCurrentMotherKF;
464 fCurrentMotherKF->SetProductionVertex(primaryVertexImproved);
469 Double_t PsiPair=GetPsiPair(fCurrentV0,fCurrentExternalTrackParamPositive,fCurrentExternalTrackParamNegative);
470 fCurrentMotherKF->SetPsiPair(PsiPair);
473 // Recalculate ConversionPoint
474 Double_t dca[2]={0,0};
475 if(fUseOwnXYZCalculation){
476 Double_t convpos[3]={0,0,0};
477 if(!GetConversionPoint(fCurrentExternalTrackParamPositive,fCurrentExternalTrackParamNegative,convpos,dca)){
478 fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kConvPointFail);
479 delete fCurrentMotherKF;
480 fCurrentMotherKF=NULL;
484 fCurrentMotherKF->SetConversionPoint(convpos);
487 if(fCurrentMotherKF->GetNDF() > 0.)
488 fCurrentMotherKF->SetChi2perNDF(fCurrentMotherKF->GetChi2()/fCurrentMotherKF->GetNDF()); //->Photon is created before all chi2 relevant changes are performed, set it "by hand"
491 // Set Dilepton Mass (moved down for same eta compared to old)
492 fCurrentMotherKF->SetMass(fCurrentMotherKF->M());
496 if(!fConversionCuts->PhotonCuts(fCurrentMotherKF,fInputEvent)){
497 fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kPhotonCuts);
498 delete fCurrentMotherKF;
499 fCurrentMotherKF=NULL;
503 fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kPhotonOut);
504 return fCurrentMotherKF;
507 ///________________________________________________________________________
508 Double_t AliV0ReaderV1::GetPsiPair(const AliESDv0* v0, const AliExternalTrackParam *positiveparam,const AliExternalTrackParam *negativeparam) const {
510 // Angle between daughter momentum plane and plane
513 AliExternalTrackParam nt(*negativeparam);
514 AliExternalTrackParam pt(*positiveparam);
516 Float_t magField = fInputEvent->GetMagneticField();
518 Double_t xyz[3] = {0.,0.,0.};
519 v0->GetXYZ(xyz[0],xyz[1],xyz[2]);
521 // Double_t pPlus[3] = {pt.Px(),pt.Py(),pt.Pz()};
522 // Double_t pMinus[3] = {nt.Px(),nt.Py(),nt.Pz()};
524 // Double_t u[3] = {pPlus[0]+pMinus[0],pPlus[1]+pMinus[1],pPlus[2]+pMinus[2]};
525 // Double_t normu = sqrt( (u[0]*u[0]) + (u[1]*u[1]) + (u[2]*u[2]) );
527 // u[0] = u[0] / normu;
528 // u[1] = u[1] / normu;
529 // u[2] = u[2] / normu;
531 // Double_t normpPlus = sqrt( (pPlus[0]*pPlus[0]) + (pPlus[1]*pPlus[1]) + (pPlus[2]*pPlus[2]) );
532 // Double_t normpMinus = sqrt( (pMinus[0]*pMinus[0]) + (pMinus[1]*pMinus[1]) + (pMinus[2]*pMinus[2]) );
534 // pPlus[0] = pPlus[0] / normpPlus;
535 // pPlus[1] = pPlus[1] / normpPlus;
536 // pPlus[2] = pPlus[2] / normpPlus;
538 // pMinus[0] = pMinus[0] / normpMinus;
539 // pMinus[1] = pMinus[1] / normpMinus;
540 // pMinus[2] = pMinus[2] / normpMinus;
542 // Double_t v[3] = {0,0,0}; // pPlus X pMinus
543 // v[0] = (pPlus[1]*pMinus[2]) - (pPlus[2]*pMinus[1]);
544 // v[1] = (pPlus[2]*pMinus[0]) - (pPlus[0]*pMinus[2]);
545 // v[2] = (pPlus[0]*pMinus[1]) - (pPlus[1]*pMinus[0]);
547 // Double_t w[3] = {0,0,0}; // u X v
548 // w[0] = (u[1]*v[2]) - (u[2]*v[1]);
549 // w[1] = (u[2]*v[0]) - (u[0]*v[2]);
550 // w[2] = (u[0]*v[1]) - (u[1]*v[0]);
552 // Double_t z[3] = {0,0,1};
553 // Double_t wc[3] = {0,0,0}; // u X v
554 // wc[0] = (u[1]*z[2]) - (u[2]*z[1]);
555 // wc[1] = (u[2]*z[0]) - (u[0]*z[2]);
556 // wc[2] = (u[0]*z[1]) - (u[1]*z[0]);
558 // Double_t PhiV = TMath::ACos((w[0]*wc[0]) + (w[1]*wc[1]) + (w[2]*wc[2]));
562 // TVector3 pPlus(pt.Px(),pt.Py(),pt.Pz());
563 // TVector3 pMinus(nt.Px(),nt.Py(),nt.Pz());
565 // TVector3 u = pMinus + pPlus;
566 // u = u*(1/u.Mag());
568 // TVector3 pHPlus = pPlus*(1/pPlus.Mag());
569 // TVector3 pHMinus = pMinus*(1/pMinus.Mag());
571 // TVector3 v = pHPlus.Cross(pHMinus);
572 // TVector3 w = u.Cross(v);
573 // TVector3 z(0,0,1);
574 // TVector3 wc = u.Cross(z);
576 // Double_t PhiV = w * wc;
578 Double_t mn[3] = {0,0,0};
579 Double_t mp[3] = {0,0,0};
581 v0->GetNPxPyPz(mn[0],mn[1],mn[2]);//reconstructed cartesian momentum components of negative daughter;
582 v0->GetPPxPyPz(mp[0],mp[1],mp[2]);//reconstructed cartesian momentum components of positive daughter;
584 Double_t deltat = 1.;
585 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
586 Double_t radiussum = TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]) + 50;//radius to which tracks shall be propagated
588 Double_t momPosProp[3] = {0,0,0};
589 Double_t momNegProp[3] = {0,0,0};
591 Double_t psiPair = 4.;
592 if(nt.PropagateTo(radiussum,magField) == 0) return psiPair; //propagate tracks to the outside -> Better Purity and Efficiency
593 if(pt.PropagateTo(radiussum,magField) == 0) return psiPair; //propagate tracks to the outside -> Better Purity and Efficiency
595 pt.GetPxPyPz(momPosProp);//Get momentum vectors of tracks after propagation
596 nt.GetPxPyPz(momNegProp);
599 TMath::Sqrt(momNegProp[0]*momNegProp[0]+momNegProp[1]*momNegProp[1]+momNegProp[2]*momNegProp[2]);//absolute momentum value of negative daughter
602 TMath::Sqrt(momPosProp[0]*momPosProp[0]+momPosProp[1]*momPosProp[1]+momPosProp[2]*momPosProp[2]);//absolute momentum value of positive daughter
604 Double_t scalarproduct =
605 momPosProp[0]*momNegProp[0]+momPosProp[1]*momNegProp[1]+momPosProp[2]*momNegProp[2];//scalar product of propagated positive and negative daughters' momenta
607 Double_t chipair = TMath::ACos(scalarproduct/(pEle*pPos));//Angle between propagated daughter tracks
609 psiPair = TMath::Abs(TMath::ASin(deltat/chipair));
614 ///________________________________________________________________________
615 Bool_t AliV0ReaderV1::GetHelixCenter(const AliExternalTrackParam *track,Double_t center[2]){
617 // Get Center of the helix track parametrization
619 Int_t charge=track->Charge();
620 Double_t b=fInputEvent->GetMagneticField();
623 track->GetHelixParameters(helix,b);
625 Double_t xpos = helix[5];
626 Double_t ypos = helix[0];
627 Double_t radius = TMath::Abs(1./helix[4]);
628 Double_t phi = helix[2];
631 phi = phi + 2*TMath::Pi();
634 phi -= TMath::Pi()/2.;
635 Double_t xpoint = radius * TMath::Cos(phi);
636 Double_t ypoint = radius * TMath::Sin(phi);
660 center[0] = xpos + xpoint;
661 center[1] = ypos + ypoint;
665 ///________________________________________________________________________
666 Bool_t AliV0ReaderV1::GetConversionPoint(const AliExternalTrackParam *pparam,const AliExternalTrackParam *nparam,Double_t convpos[3],Double_t dca[2]){
668 // Recalculate Conversion Point
670 if(!pparam||!nparam)return kFALSE;
672 Double_t helixcenterpos[2];
673 GetHelixCenter(pparam,helixcenterpos);
675 Double_t helixcenterneg[2];
676 GetHelixCenter(nparam,helixcenterneg);
678 Double_t helixpos[6];
679 pparam->GetHelixParameters(helixpos,fInputEvent->GetMagneticField());
680 Double_t posradius = TMath::Abs(1./helixpos[4]);
682 Double_t helixneg[6];
683 nparam->GetHelixParameters(helixneg,fInputEvent->GetMagneticField());
684 Double_t negradius = TMath::Abs(1./helixneg[4]);
686 // Calculate xy-position
688 Double_t xpos = helixcenterpos[0];
689 Double_t ypos = helixcenterpos[1];
690 Double_t xneg = helixcenterneg[0];
691 Double_t yneg = helixcenterneg[1];
693 convpos[0] = (xpos*negradius + xneg*posradius)/(negradius+posradius);
694 convpos[1] = (ypos*negradius+ yneg*posradius)/(negradius+posradius);
697 // Calculate Track XY vertex position
699 Double_t deltaXPos = convpos[0] - xpos;
700 Double_t deltaYPos = convpos[1] - ypos;
702 Double_t deltaXNeg = convpos[0] - xneg;
703 Double_t deltaYNeg = convpos[1] - yneg;
705 Double_t alphaPos = TMath::Pi() + TMath::ATan2(-deltaYPos,-deltaXPos);
706 Double_t alphaNeg = TMath::Pi() + TMath::ATan2(-deltaYNeg,-deltaXNeg);
708 Double_t vertexXNeg = xneg + TMath::Abs(negradius)*TMath::Cos(alphaNeg);
709 Double_t vertexYNeg = yneg + TMath::Abs(negradius)*TMath::Sin(alphaNeg);
711 Double_t vertexXPos = xpos + TMath::Abs(posradius)*TMath::Cos(alphaPos);
712 Double_t vertexYPos = ypos + TMath::Abs(posradius)*TMath::Sin(alphaPos);
714 AliExternalTrackParam p(*pparam);
715 AliExternalTrackParam n(*nparam);
717 TVector2 vertexPos(vertexXPos,vertexYPos);
718 TVector2 vertexNeg(vertexXNeg,vertexYNeg);
720 // Convert to local coordinate system
721 TVector2 vertexPosRot=vertexPos.Rotate(-p.GetAlpha());
722 TVector2 vertexNegRot=vertexNeg.Rotate(-n.GetAlpha());
724 // Propagate Track Params to Vertex
726 if(!p.PropagateTo(vertexPosRot.X(),fInputEvent->GetMagneticField()))return kFALSE;
727 if(!n.PropagateTo(vertexNegRot.X(),fInputEvent->GetMagneticField()))return kFALSE;
729 // Check whether propagation was sucessful
731 if(TMath::Abs(vertexPos.Mod()-TMath::Sqrt(p.GetX()*p.GetX()+p.GetY()*p.GetY()))>0.01)return kFALSE;
732 if(TMath::Abs(vertexNeg.Mod()-TMath::Sqrt(n.GetX()*n.GetX()+n.GetY()*n.GetY()))>0.01)return kFALSE;
734 // Calculate z position
736 convpos[2] = (p.GetZ()*negradius+n.GetZ()*posradius)/(negradius+posradius);
739 TVector2 vdca=vertexPos-vertexNeg;
741 dca[1]=TMath::Abs(n.GetZ()-p.GetZ());
745 //________________________________________________________________________
746 Bool_t AliV0ReaderV1::GetAODConversionGammas(){
748 // Get reconstructed conversion photons from satellite AOD file
750 AliAODEvent *fAODEvent=dynamic_cast<AliAODEvent*>(fInputEvent);
754 if(fConversionGammas == NULL){
755 fConversionGammas = new TClonesArray("AliAODConversionPhoton",100);
757 fConversionGammas->Delete();//Reset the TClonesArray
759 //Get Gammas from satellite AOD gamma branch
761 AliAODConversionPhoton *gamma=0x0;
763 TClonesArray *fInputGammas=dynamic_cast<TClonesArray*>(fAODEvent->FindListObject(fDeltaAODBranchName.Data()));
765 FindDeltaAODBranchName();
766 fInputGammas=dynamic_cast<TClonesArray*>(fAODEvent->FindListObject(fDeltaAODBranchName.Data()));}
767 if(!fInputGammas){AliError("No Gamma Satellites found");return kFALSE;}
768 // Apply Selection Cuts to Gammas and create local working copy
770 for(Int_t i=0;i<fInputGammas->GetEntriesFast();i++){
771 gamma=dynamic_cast<AliAODConversionPhoton*>(fInputGammas->At(i));
773 if(fConversionCuts->PhotonIsSelected(gamma,fInputEvent)){
774 new((*fConversionGammas)[fConversionGammas->GetEntriesFast()]) AliAODConversionPhoton(*gamma);}
780 if(fConversionGammas->GetEntries()){return kTRUE;}
785 //________________________________________________________________________
786 void AliV0ReaderV1::FindDeltaAODBranchName(){
788 // Find delta AOD branchname containing reconstructed photons
790 TList *list=fInputEvent->GetList();
791 for(Int_t ii=0;ii<list->GetEntries();ii++){
792 TString name((list->At(ii))->GetName());
793 if(name.BeginsWith(fDeltaAODBranchName)&&name.EndsWith("gamma")){
794 fDeltaAODBranchName=name;
795 AliInfo(Form("Set DeltaAOD BranchName to: %s",fDeltaAODBranchName.Data()));
802 //________________________________________________________________________
803 void AliV0ReaderV1::Terminate(Option_t *)