1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Authors: Svein Lindal, Daniel Lohner *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
16 ////////////////////////////////////////////////
17 //---------------------------------------------
18 // Class reconstructing conversion photons from V0s
19 //---------------------------------------------
20 ////////////////////////////////////////////////
22 #include "AliV0ReaderV1.h"
23 #include "AliKFParticle.h"
26 #include "AliAODEvent.h"
27 #include "AliESDEvent.h"
28 #include "AliKFParticle.h"
29 #include "AliKFConversionPhoton.h"
30 #include "AliAODConversionPhoton.h"
31 #include "AliConversionPhotonBase.h"
33 #include "AliKFVertex.h"
34 #include "AliAODTrack.h"
35 #include "AliESDtrack.h"
36 #include "AliAnalysisManager.h"
37 #include "AliInputEventHandler.h"
38 #include "AliAODHandler.h"
39 #include "AliPIDResponse.h"
48 ClassImp(AliV0ReaderV1)
50 //________________________________________________________________________
51 AliV0ReaderV1::AliV0ReaderV1(const char *name) : AliAnalysisTaskSE(name),
52 fConversionCuts(NULL),
53 fConversionGammas(NULL),
54 fUseImprovedVertex(kTRUE),
55 fUseOwnXYZCalculation(kTRUE),
56 fUseConstructGamma(kFALSE),
57 kUseAODConversionPhoton(kTRUE),
59 fDeltaAODBranchName("GammaConv"),
60 fDeltaAODFilename("AliAODGammaConversion.root"),
61 fEventIsSelected(kFALSE)
63 // Default constructor
65 DefineInput(0, TChain::Class());
68 //________________________________________________________________________
69 AliV0ReaderV1::~AliV0ReaderV1()
71 // default deconstructor
73 if(fConversionGammas){
74 fConversionGammas->Delete();// Clear Objects
75 delete fConversionGammas;
76 fConversionGammas=0x0;
80 //________________________________________________________________________
81 AliV0ReaderV1::AliV0ReaderV1(AliV0ReaderV1 &original) : AliAnalysisTaskSE(original),
82 fConversionCuts(NULL),
83 fConversionGammas(NULL),
84 fUseImprovedVertex(original.fUseImprovedVertex),
85 fUseOwnXYZCalculation(original.fUseOwnXYZCalculation),
86 fUseConstructGamma(original.fUseConstructGamma),
87 kUseAODConversionPhoton(original.kUseAODConversionPhoton),
88 fCreateAOD(original.fCreateAOD),
89 fDeltaAODBranchName(original.fDeltaAODBranchName),
90 fDeltaAODFilename(original.fDeltaAODFilename),
91 fEventIsSelected(original.fEventIsSelected)
93 // Default constructor
95 DefineInput(0, TChain::Class());
98 //____________________________________________________________
99 AliV0ReaderV1 &AliV0ReaderV1::operator=(const AliV0ReaderV1 &ref){
101 // Assignment operator
102 // Only copies pointers, object is not the owner of the references
105 AliAnalysisTaskSE::operator=(ref);
106 fUseImprovedVertex=ref.fUseImprovedVertex;
107 fUseOwnXYZCalculation=ref.fUseOwnXYZCalculation;
108 fUseConstructGamma=ref.fUseConstructGamma;
109 kUseAODConversionPhoton=ref.kUseAODConversionPhoton;
110 fCreateAOD=ref.fCreateAOD;
111 fDeltaAODBranchName=ref.fDeltaAODBranchName;
112 fDeltaAODFilename=ref.fDeltaAODFilename;
113 fEventIsSelected=ref.fEventIsSelected;
118 //________________________________________________________________________
119 void AliV0ReaderV1::Init()
121 // Initialize function to be called once before analysis
123 if(fConversionCuts==NULL){
124 if(fConversionCuts==NULL)AliError("No Cut Selection initialized");
127 if(fCreateAOD){kUseAODConversionPhoton=kTRUE;}
129 if(fConversionGammas != NULL){
130 delete fConversionGammas;
131 fConversionGammas=NULL;
134 if(fConversionGammas == NULL){
135 if(kUseAODConversionPhoton){
136 fConversionGammas = new TClonesArray("AliAODConversionPhoton",100);}
138 fConversionGammas = new TClonesArray("AliKFConversionPhoton",100);}
140 fConversionGammas->Delete();//Reset the TClonesArray
144 //________________________________________________________________________
145 void AliV0ReaderV1::UserCreateOutputObjects()
151 fDeltaAODBranchName.Append("_");
152 fDeltaAODBranchName.Append(fConversionCuts->GetCutNumber());
153 fDeltaAODBranchName.Append("_gamma");
155 fConversionGammas->SetName(fDeltaAODBranchName.Data());
157 AddAODBranch("TClonesArray", &fConversionGammas, fDeltaAODFilename.Data());
158 AliAnalysisManager::GetAnalysisManager()->RegisterExtraFile(fDeltaAODFilename.Data());
163 //________________________________________________________________________
164 void AliV0ReaderV1::UserExec(Option_t *){
166 // Check if correctly initialized
167 if(!fConversionGammas)Init();
170 fEventIsSelected=ProcessEvent(fInputEvent,fMCEvent);
173 //________________________________________________________________________
174 Bool_t AliV0ReaderV1::ProcessEvent(AliVEvent *inputEvent,AliMCEvent *mcEvent)
176 //Reset the TClonesArray
177 fConversionGammas->Delete();
179 fInputEvent=inputEvent;
183 AliError("No Input event");
187 if(!fConversionCuts){AliError("No ConversionCuts");return kFALSE;}
190 if(!fConversionCuts->EventIsSelected(fInputEvent,fMCEvent))return kFALSE;
192 // Set Magnetic Field
193 AliKFParticle::SetField(fInputEvent->GetMagneticField());
195 if(fInputEvent->IsA()==AliESDEvent::Class()){
198 if(fInputEvent->IsA()==AliAODEvent::Class()){
199 GetAODConversionGammas();
207 ///________________________________________________________________________
208 void AliV0ReaderV1::FillAODOutput()
210 // Fill AOD Output with reconstructed Photons
212 if(fInputEvent->IsA()==AliESDEvent::Class()){
213 ///Make sure delta aod is filled if standard aod is filled (for synchronization when reading aod with standard aod)
215 AliAODHandler * aodhandler = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
216 if (aodhandler && aodhandler->GetFillAOD()) {
217 AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillExtension(kTRUE);
218 //PostData(0, fConversionGammas);
225 ///________________________________________________________________________
226 const AliExternalTrackParam *AliV0ReaderV1::GetExternalTrackParam(AliESDv0 *fCurrentV0,Int_t &tracklabel,Int_t charge){
228 // Get External Track Parameter with given charge
230 if(!(charge==1||charge==-1)){AliError("Charge not defined");return 0x0;}
232 // Check for sign flip
234 if(!fCurrentV0->GetParamN()||!fCurrentV0->GetParamP())return 0x0;
235 if(!fConversionCuts->GetTrack(fInputEvent,fCurrentV0->GetNindex())||!fConversionCuts->GetTrack(fInputEvent,fCurrentV0->GetPindex()))return 0x0;
236 if((fConversionCuts->GetTrack(fInputEvent,fCurrentV0->GetPindex()))->Charge()==charge){
237 tracklabel=fCurrentV0->GetPindex();
238 return fCurrentV0->GetParamP();}
239 if((fConversionCuts->GetTrack(fInputEvent,fCurrentV0->GetNindex()))->Charge()==charge){
240 tracklabel=fCurrentV0->GetNindex();
241 return fCurrentV0->GetParamN();}
246 ///________________________________________________________________________
247 Bool_t AliV0ReaderV1::ProcessESDV0s()
249 // Process ESD V0s for conversion photon reconstruction
251 AliESDEvent *fESDEvent=dynamic_cast<AliESDEvent*>(fInputEvent);
253 AliKFConversionPhoton *fCurrentMotherKFCandidate=NULL;
257 for(Int_t currentV0Index=0;currentV0Index<fESDEvent->GetNumberOfV0s();currentV0Index++){
258 AliESDv0 *fCurrentV0=(AliESDv0*)(fESDEvent->GetV0(currentV0Index));
260 printf("Requested V0 does not exist");
263 fCurrentMotherKFCandidate=ReconstructV0(fCurrentV0,currentV0Index);
265 if(fCurrentMotherKFCandidate){
267 // Add Gamma to the TClonesArray
269 if(kUseAODConversionPhoton){
270 new((*fConversionGammas)[fConversionGammas->GetEntriesFast()]) AliAODConversionPhoton(fCurrentMotherKFCandidate);
273 new((*fConversionGammas)[fConversionGammas->GetEntriesFast()]) AliKFConversionPhoton(*fCurrentMotherKFCandidate);
276 delete fCurrentMotherKFCandidate;
277 fCurrentMotherKFCandidate=NULL;
284 ///________________________________________________________________________
285 AliKFConversionPhoton *AliV0ReaderV1::ReconstructV0(AliESDv0 *fCurrentV0,Int_t currentV0Index)
287 // Reconstruct conversion photon from ESD v0
288 fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kPhotonIn);
290 //checks if on the fly mode is set
291 if(!fConversionCuts->SelectV0Finder(fCurrentV0->GetOnFlyStatus())){
292 fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kOnFly);
297 Int_t currentTrackLabels[2]={-1,-1};
299 // Get Daughter KF Particles
301 const AliExternalTrackParam *fCurrentExternalTrackParamPositive=GetExternalTrackParamP(fCurrentV0,currentTrackLabels[0]);
302 const AliExternalTrackParam *fCurrentExternalTrackParamNegative=GetExternalTrackParamN(fCurrentV0,currentTrackLabels[1]);
304 if(!fCurrentExternalTrackParamPositive||!fCurrentExternalTrackParamNegative)return 0x0;
306 // Apply some Cuts before Reconstruction
308 AliVTrack * posTrack = fConversionCuts->GetTrack(fInputEvent,currentTrackLabels[0]);
309 AliVTrack * negTrack = fConversionCuts->GetTrack(fInputEvent,currentTrackLabels[1]);
311 if(!negTrack || !posTrack) {
312 fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kNoTracks);
317 if(!fConversionCuts->TracksAreSelected(negTrack, posTrack)){
318 fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kTrackCuts);
323 if(!fConversionCuts->dEdxCuts(negTrack) || !fConversionCuts->dEdxCuts(posTrack)) {
324 fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kdEdxCuts);
328 // Reconstruct Photon
330 AliKFConversionPhoton *fCurrentMotherKF=NULL;
332 AliKFParticle fCurrentNegativeKFParticle(*(fCurrentExternalTrackParamNegative),11);
333 AliKFParticle fCurrentPositiveKFParticle(*(fCurrentExternalTrackParamPositive),-11);
337 if(fUseConstructGamma){
339 fCurrentMotherKF = new AliKFConversionPhoton();
340 fCurrentMotherKF->ConstructGamma(fCurrentNegativeKFParticle,fCurrentPositiveKFParticle);
342 fCurrentMotherKF = new AliKFConversionPhoton(fCurrentNegativeKFParticle,fCurrentPositiveKFParticle);
343 fCurrentMotherKF->SetMassConstraint(0,0);
348 fCurrentMotherKF->SetTrackLabels(currentTrackLabels[0],currentTrackLabels[1]);
352 fCurrentMotherKF->SetV0Index(currentV0Index);
358 AliStack *fMCStack= fMCEvent->Stack();
360 Int_t labeln=TMath::Abs(fConversionCuts->GetTrack(fInputEvent,fCurrentMotherKF->GetTrackLabelPositive())->GetLabel());
361 Int_t labelp=TMath::Abs(fConversionCuts->GetTrack(fInputEvent,fCurrentMotherKF->GetTrackLabelNegative())->GetLabel());
363 TParticle *fNegativeMCParticle = fMCStack->Particle(labeln);
364 TParticle *fPositiveMCParticle = fMCStack->Particle(labelp);
366 if(fPositiveMCParticle&&fNegativeMCParticle){
367 fCurrentMotherKF->SetMCLabelPositive(labelp);
368 fCurrentMotherKF->SetMCLabelNegative(labeln);
372 // Update Vertex (moved for same eta compared to old)
373 if(fUseImprovedVertex == kTRUE){
374 AliKFVertex primaryVertexImproved(*fInputEvent->GetPrimaryVertex());
375 primaryVertexImproved+=*fCurrentMotherKF;
376 fCurrentMotherKF->SetProductionVertex(primaryVertexImproved);
381 Double_t PsiPair=GetPsiPair(fCurrentV0,fCurrentExternalTrackParamPositive,fCurrentExternalTrackParamNegative);
382 fCurrentMotherKF->SetPsiPair(PsiPair);
385 // Recalculate ConversionPoint
386 Double_t dca[2]={0,0};
387 if(fUseOwnXYZCalculation){
388 Double_t convpos[3]={0,0,0};
389 if(!GetConversionPoint(fCurrentExternalTrackParamPositive,fCurrentExternalTrackParamNegative,convpos,dca)){
390 fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kConvPointFail);
391 delete fCurrentMotherKF;
392 fCurrentMotherKF=NULL;
396 fCurrentMotherKF->SetConversionPoint(convpos);
399 if(fCurrentMotherKF->GetNDF() > 0.)
400 fCurrentMotherKF->SetChi2perNDF(fCurrentMotherKF->GetChi2()/fCurrentMotherKF->GetNDF()); //->Photon is created before all chi2 relevant changes are performed, set it "by hand"
404 if(!fConversionCuts->PhotonCuts(fCurrentMotherKF,fInputEvent)){
405 fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kPhotonCuts);
406 delete fCurrentMotherKF;
407 fCurrentMotherKF=NULL;
411 // Set Dilepton Mass (moved down for same eta compared to old)
412 fCurrentMotherKF->SetMass(fCurrentMotherKF->M());
415 fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kPhotonOut);
416 return fCurrentMotherKF;
419 ///________________________________________________________________________
420 Double_t AliV0ReaderV1::GetPsiPair(const AliESDv0* v0, const AliExternalTrackParam *positiveparam,const AliExternalTrackParam *negativeparam) const {
422 // Angle between daughter momentum plane and plane
425 AliExternalTrackParam nt(*negativeparam);
426 AliExternalTrackParam pt(*positiveparam);
428 Float_t magField = fInputEvent->GetMagneticField();
430 Double_t xyz[3] = {0.,0.,0.};
431 v0->GetXYZ(xyz[0],xyz[1],xyz[2]);
433 Double_t mn[3] = {0,0,0};
434 Double_t mp[3] = {0,0,0};
436 v0->GetNPxPyPz(mn[0],mn[1],mn[2]);//reconstructed cartesian momentum components of negative daughter;
437 v0->GetPPxPyPz(mp[0],mp[1],mp[2]);//reconstructed cartesian momentum components of positive daughter;
439 Double_t deltat = 1.;
440 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
441 Double_t radiussum = TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]) + 50;//radius to which tracks shall be propagated
443 Double_t momPosProp[3] = {0,0,0};
444 Double_t momNegProp[3] = {0,0,0};
446 Double_t psiPair = 4.;
447 if(nt.PropagateTo(radiussum,magField) == 0) return psiPair; //propagate tracks to the outside -> Better Purity and Efficiency
449 if(pt.PropagateTo(radiussum,magField) == 0) return psiPair; //propagate tracks to the outside -> Better Purity and Efficiency
451 pt.GetPxPyPz(momPosProp);//Get momentum vectors of tracks after propagation
452 nt.GetPxPyPz(momNegProp);
455 TMath::Sqrt(momNegProp[0]*momNegProp[0]+momNegProp[1]*momNegProp[1]+momNegProp[2]*momNegProp[2]);//absolute momentum value of negative daughter
458 TMath::Sqrt(momPosProp[0]*momPosProp[0]+momPosProp[1]*momPosProp[1]+momPosProp[2]*momPosProp[2]);//absolute momentum value of positive daughter
460 Double_t scalarproduct =
461 momPosProp[0]*momNegProp[0]+momPosProp[1]*momNegProp[1]+momPosProp[2]*momNegProp[2];//scalar product of propagated positive and negative daughters' momenta
463 Double_t chipair = TMath::ACos(scalarproduct/(pEle*pPos));//Angle between propagated daughter tracks
465 psiPair = TMath::Abs(TMath::ASin(deltat/chipair));
470 ///________________________________________________________________________
471 Bool_t AliV0ReaderV1::GetHelixCenter(const AliExternalTrackParam *track,Double_t center[2]){
473 // Get Center of the helix track parametrization
475 Int_t charge=track->Charge();
476 Double_t b=fInputEvent->GetMagneticField();
479 track->GetHelixParameters(helix,b);
481 Double_t xpos = helix[5];
482 Double_t ypos = helix[0];
483 Double_t radius = TMath::Abs(1./helix[4]);
484 Double_t phi = helix[2];
487 phi = phi + 2*TMath::Pi();
490 phi -= TMath::Pi()/2.;
491 Double_t xpoint = radius * TMath::Cos(phi);
492 Double_t ypoint = radius * TMath::Sin(phi);
516 center[0] = xpos + xpoint;
517 center[1] = ypos + ypoint;
521 ///________________________________________________________________________
522 Bool_t AliV0ReaderV1::GetConversionPoint(const AliExternalTrackParam *pparam,const AliExternalTrackParam *nparam,Double_t convpos[3],Double_t dca[2]){
524 // Recalculate Conversion Point
526 if(!pparam||!nparam)return kFALSE;
528 Double_t helixcenterpos[2];
529 GetHelixCenter(pparam,helixcenterpos);
531 Double_t helixcenterneg[2];
532 GetHelixCenter(nparam,helixcenterneg);
534 Double_t helixpos[6];
535 pparam->GetHelixParameters(helixpos,fInputEvent->GetMagneticField());
536 Double_t posradius = TMath::Abs(1./helixpos[4]);
538 Double_t helixneg[6];
539 nparam->GetHelixParameters(helixneg,fInputEvent->GetMagneticField());
540 Double_t negradius = TMath::Abs(1./helixneg[4]);
542 // Calculate xy-position
544 Double_t xpos = helixcenterpos[0];
545 Double_t ypos = helixcenterpos[1];
546 Double_t xneg = helixcenterneg[0];
547 Double_t yneg = helixcenterneg[1];
549 convpos[0] = (xpos*negradius + xneg*posradius)/(negradius+posradius);
550 convpos[1] = (ypos*negradius+ yneg*posradius)/(negradius+posradius);
553 // Calculate Track XY vertex position
555 Double_t deltaXPos = convpos[0] - xpos;
556 Double_t deltaYPos = convpos[1] - ypos;
558 Double_t deltaXNeg = convpos[0] - xneg;
559 Double_t deltaYNeg = convpos[1] - yneg;
561 Double_t alphaPos = TMath::Pi() + TMath::ATan2(-deltaYPos,-deltaXPos);
562 Double_t alphaNeg = TMath::Pi() + TMath::ATan2(-deltaYNeg,-deltaXNeg);
564 Double_t vertexXNeg = xneg + TMath::Abs(negradius)*TMath::Cos(alphaNeg);
565 Double_t vertexYNeg = yneg + TMath::Abs(negradius)*TMath::Sin(alphaNeg);
567 Double_t vertexXPos = xpos + TMath::Abs(posradius)*TMath::Cos(alphaPos);
568 Double_t vertexYPos = ypos + TMath::Abs(posradius)*TMath::Sin(alphaPos);
570 AliExternalTrackParam p(*pparam);
571 AliExternalTrackParam n(*nparam);
573 TVector2 vertexPos(vertexXPos,vertexYPos);
574 TVector2 vertexNeg(vertexXNeg,vertexYNeg);
576 // Convert to local coordinate system
577 TVector2 vertexPosRot=vertexPos.Rotate(-p.GetAlpha());
578 TVector2 vertexNegRot=vertexNeg.Rotate(-n.GetAlpha());
580 // Propagate Track Params to Vertex
582 if(!p.PropagateTo(vertexPosRot.X(),fInputEvent->GetMagneticField()))return kFALSE;
583 if(!n.PropagateTo(vertexNegRot.X(),fInputEvent->GetMagneticField()))return kFALSE;
585 // Check whether propagation was sucessful
587 if(TMath::Abs(vertexPos.Mod()-TMath::Sqrt(p.GetX()*p.GetX()+p.GetY()*p.GetY()))>0.01)return kFALSE;
588 if(TMath::Abs(vertexNeg.Mod()-TMath::Sqrt(n.GetX()*n.GetX()+n.GetY()*n.GetY()))>0.01)return kFALSE;
590 // Calculate z position
592 convpos[2] = (p.GetZ()*negradius+n.GetZ()*posradius)/(negradius+posradius);
595 TVector2 vdca=vertexPos-vertexNeg;
597 dca[1]=TMath::Abs(n.GetZ()-p.GetZ());
601 //________________________________________________________________________
602 Bool_t AliV0ReaderV1::GetAODConversionGammas(){
604 // Get reconstructed conversion photons from satellite AOD file
606 AliAODEvent *fAODEvent=dynamic_cast<AliAODEvent*>(fInputEvent);
610 if(fConversionGammas == NULL){
611 fConversionGammas = new TClonesArray("AliAODConversionPhoton",100);
613 fConversionGammas->Delete();//Reset the TClonesArray
615 //Get Gammas from satellite AOD gamma branch
617 AliAODConversionPhoton *gamma=0x0;
619 TClonesArray *fInputGammas=dynamic_cast<TClonesArray*>(fAODEvent->FindListObject(fDeltaAODBranchName.Data()));
621 FindDeltaAODBranchName();
622 fInputGammas=dynamic_cast<TClonesArray*>(fAODEvent->FindListObject(fDeltaAODBranchName.Data()));}
623 if(!fInputGammas){AliError("No Gamma Satellites found");return kFALSE;}
624 // Apply Selection Cuts to Gammas and create local working copy
626 for(Int_t i=0;i<fInputGammas->GetEntriesFast();i++){
627 gamma=dynamic_cast<AliAODConversionPhoton*>(fInputGammas->At(i));
629 if(fConversionCuts->PhotonIsSelected(gamma,fInputEvent)){
630 new((*fConversionGammas)[fConversionGammas->GetEntriesFast()]) AliAODConversionPhoton(*gamma);}
636 if(fConversionGammas->GetEntries()){return kTRUE;}
641 //________________________________________________________________________
642 void AliV0ReaderV1::FindDeltaAODBranchName(){
644 // Find delta AOD branchname containing reconstructed photons
646 TList *list=fInputEvent->GetList();
647 for(Int_t ii=0;ii<list->GetEntries();ii++){
648 TString name((list->At(ii))->GetName());
649 if(name.BeginsWith(fDeltaAODBranchName)&&name.EndsWith("gamma")){
650 fDeltaAODBranchName=name;
651 AliInfo(Form("Set DeltaAOD BranchName to: %s",fDeltaAODBranchName.Data()));
658 //________________________________________________________________________
659 void AliV0ReaderV1::Terminate(Option_t *)