]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGGA/GammaConv/AliV0ReaderV1.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Conversion Task able to run on AOD's, added different trigger selection, implemented...
[u/mrichter/AliRoot.git] / PWGGA / GammaConv / AliV0ReaderV1.cxx
CommitLineData
1528f6d1 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
11c1e680 3 * *
4 * Authors: Svein Lindal, Daniel Lohner *
5 * Version 1.0 *
6 * *
7 * *
8 * based on: on older version (see aliroot up to v5-04-42-AN) *
9 * AliV0Reader.cxx *
10 * Authors: Kathrin Koch, Kenneth Aamodt, Ana Marin *
11 * *
e5b6e8a6 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. *
1528f6d1 19 **************************************************************************/
20
21////////////////////////////////////////////////
e5b6e8a6 22//---------------------------------------------
1528f6d1 23// Class reconstructing conversion photons from V0s
24//---------------------------------------------
25////////////////////////////////////////////////
26
27#include "AliV0ReaderV1.h"
28#include "AliKFParticle.h"
29#include "AliAODv0.h"
30#include "AliESDv0.h"
31#include "AliAODEvent.h"
32#include "AliESDEvent.h"
0a2b2b4b 33#include "AliPID.h"
34#include "AliMCEvent.h"
35#include "AliStack.h"
36#include "AliMCEventHandler.h"
37#include "AliESDpid.h"
38#include "AliESDtrackCuts.h"
39#include "TRandom3.h"
40#include "AliGenCocktailEventHeader.h"
41#include "TList.h"
1528f6d1 42#include "AliKFConversionPhoton.h"
43#include "AliAODConversionPhoton.h"
44#include "AliConversionPhotonBase.h"
45#include "TVector.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"
53#include "TChain.h"
a280ac15 54#include "TFile.h"
a280ac15 55#include "TString.h"
56#include "TObjArray.h"
1528f6d1 57
58class iostream;
59
60
61using namespace std;
62
63ClassImp(AliV0ReaderV1)
64
65//________________________________________________________________________
66AliV0ReaderV1::AliV0ReaderV1(const char *name) : AliAnalysisTaskSE(name),
67 fConversionCuts(NULL),
68 fConversionGammas(NULL),
69 fUseImprovedVertex(kTRUE),
70 fUseOwnXYZCalculation(kTRUE),
71 fUseConstructGamma(kFALSE),
72 kUseAODConversionPhoton(kTRUE),
73 fCreateAOD(kFALSE),
74 fDeltaAODBranchName("GammaConv"),
75 fDeltaAODFilename("AliAODGammaConversion.root"),
ae947965 76 fCheckAODConsistenty(kFALSE),
a280ac15 77 fEventIsSelected(kFALSE),
78 fPeriodName("")
1528f6d1 79{
80 // Default constructor
81
82 DefineInput(0, TChain::Class());
83}
84
85//________________________________________________________________________
86AliV0ReaderV1::~AliV0ReaderV1()
87{
88 // default deconstructor
89
90 if(fConversionGammas){
91 fConversionGammas->Delete();// Clear Objects
92 delete fConversionGammas;
93 fConversionGammas=0x0;
94 }
95}
2bdd97ae 96/*
97//________________________________________________________________________
98AliV0ReaderV1::AliV0ReaderV1(AliV0ReaderV1 &original) : AliAnalysisTaskSE(original),
99 fConversionCuts(NULL),
100 fConversionGammas(NULL),
101 fUseImprovedVertex(original.fUseImprovedVertex),
102 fUseOwnXYZCalculation(original.fUseOwnXYZCalculation),
103 fUseConstructGamma(original.fUseConstructGamma),
104 kUseAODConversionPhoton(original.kUseAODConversionPhoton),
105 fCreateAOD(original.fCreateAOD),
106 fDeltaAODBranchName(original.fDeltaAODBranchName),
107 fDeltaAODFilename(original.fDeltaAODFilename),
108 fEventIsSelected(original.fEventIsSelected)
109{
110 // Default constructor
111
112 DefineInput(0, TChain::Class());
113}
114
115//____________________________________________________________
116AliV0ReaderV1 &AliV0ReaderV1::operator=(const AliV0ReaderV1 &ref){
117 //
118 // Assignment operator
119 // Only copies pointers, object is not the owner of the references
120 //
121 if(this != &ref){
122 AliAnalysisTaskSE::operator=(ref);
123 fUseImprovedVertex=ref.fUseImprovedVertex;
124 fUseOwnXYZCalculation=ref.fUseOwnXYZCalculation;
125 fUseConstructGamma=ref.fUseConstructGamma;
126 kUseAODConversionPhoton=ref.kUseAODConversionPhoton;
127 fCreateAOD=ref.fCreateAOD;
128 fDeltaAODBranchName=ref.fDeltaAODBranchName;
129 fDeltaAODFilename=ref.fDeltaAODFilename;
130 fEventIsSelected=ref.fEventIsSelected;
131 }
132 return *this;
133}
134*/
1528f6d1 135//________________________________________________________________________
136void AliV0ReaderV1::Init()
137{
138 // Initialize function to be called once before analysis
1528f6d1 139 if(fConversionCuts==NULL){
140 if(fConversionCuts==NULL)AliError("No Cut Selection initialized");
141 }
1528f6d1 142 if(fCreateAOD){kUseAODConversionPhoton=kTRUE;}
143
144 if(fConversionGammas != NULL){
145 delete fConversionGammas;
146 fConversionGammas=NULL;
147 }
148
149 if(fConversionGammas == NULL){
150 if(kUseAODConversionPhoton){
151 fConversionGammas = new TClonesArray("AliAODConversionPhoton",100);}
152 else{
153 fConversionGammas = new TClonesArray("AliKFConversionPhoton",100);}
154 }
155 fConversionGammas->Delete();//Reset the TClonesArray
2a1b1442 156}
157
158//________________________________________________________________________
159void AliV0ReaderV1::UserCreateOutputObjects()
160{
1528f6d1 161 // Create AODs
162
163 if(fCreateAOD){
164 if(fConversionCuts){
165 fDeltaAODBranchName.Append("_");
166 fDeltaAODBranchName.Append(fConversionCuts->GetCutNumber());
167 fDeltaAODBranchName.Append("_gamma");
168 }
169 fConversionGammas->SetName(fDeltaAODBranchName.Data());
170
171 AddAODBranch("TClonesArray", &fConversionGammas, fDeltaAODFilename.Data());
172 AliAnalysisManager::GetAnalysisManager()->RegisterExtraFile(fDeltaAODFilename.Data());
173 }
1528f6d1 174
1528f6d1 175}
1528f6d1 176//________________________________________________________________________
11c1e680 177Bool_t AliV0ReaderV1::Notify()
178{
a280ac15 179 if (fPeriodName.CompareTo("") == 0){
180 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
181 if(man) {
182 AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
183 if (inputHandler){
184 TTree* tree = (TTree*) inputHandler->GetTree();
185 TFile* file = (TFile*) tree->GetCurrentFile();
186 TString fileName(file->GetName());
187 TObjArray *arr = fileName.Tokenize("/");
188 for (Int_t i = 0; i < arr->GetEntriesFast();i++ ){
189 TObjString* testObjString = (TObjString*)arr->At(i);
190 if (testObjString->GetString().Contains("LHC")){
191 fPeriodName = testObjString->GetString();
192 i = arr->GetEntriesFast();
193 }
11c1e680 194 }
a280ac15 195 }
196 }
ccfa8c0d 197 if(!fConversionCuts->GetDoEtaShift()) return kTRUE; // No Eta Shift requested, continue
198 if(fConversionCuts->GetEtaShift() == 0.0){ // Eta Shift requested but not set, get shift automatically
199 fConversionCuts->GetCorrectEtaShiftFromPeriod(fPeriodName);
200 fConversionCuts->DoEtaShift(kFALSE); // Eta Shift Set, make sure that it is called only once
201 return kTRUE;
11c1e680 202 }
ccfa8c0d 203 else{
204 printf(" Gamma Conversion Reader %s :: Eta Shift Manually Set to %f \n\n",
205 (fConversionCuts->GetCutNumber()).Data(),fConversionCuts->GetEtaShift());
206 fConversionCuts->DoEtaShift(kFALSE); // Eta Shift Set, make sure that it is called only once
11c1e680 207 }
a280ac15 208 }
ccfa8c0d 209
11c1e680 210 return kTRUE;
211}
212//________________________________________________________________________
213void AliV0ReaderV1::UserExec(Option_t *option){
214
2a1b1442 215 // Check if correctly initialized
216 if(!fConversionGammas)Init();
217
1528f6d1 218 // User Exec
219 fEventIsSelected=ProcessEvent(fInputEvent,fMCEvent);
220}
221
222//________________________________________________________________________
223Bool_t AliV0ReaderV1::ProcessEvent(AliVEvent *inputEvent,AliMCEvent *mcEvent)
224{
ccfa8c0d 225
1528f6d1 226 //Reset the TClonesArray
227 fConversionGammas->Delete();
228
229 fInputEvent=inputEvent;
230 fMCEvent=mcEvent;
231
232 if(!fInputEvent){
233 AliError("No Input event");
234 return kFALSE;
235 }
236
237 if(!fConversionCuts){AliError("No ConversionCuts");return kFALSE;}
e5b6e8a6 238
1528f6d1 239 // Event Cuts
240 if(!fConversionCuts->EventIsSelected(fInputEvent,fMCEvent))return kFALSE;
241
242 // Set Magnetic Field
243 AliKFParticle::SetField(fInputEvent->GetMagneticField());
244
245 if(fInputEvent->IsA()==AliESDEvent::Class()){
246 ProcessESDV0s();
247 }
248 if(fInputEvent->IsA()==AliAODEvent::Class()){
249 GetAODConversionGammas();
250 }
251
252 // AOD Output
253 FillAODOutput();
254
255 return kTRUE;
256}
257///________________________________________________________________________
258void AliV0ReaderV1::FillAODOutput()
259{
260 // Fill AOD Output with reconstructed Photons
261
262 if(fInputEvent->IsA()==AliESDEvent::Class()){
263 ///Make sure delta aod is filled if standard aod is filled (for synchronization when reading aod with standard aod)
264 if(fCreateAOD) {
265 AliAODHandler * aodhandler = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
e5b6e8a6 266 if (aodhandler && aodhandler->GetFillAOD()) {
267 AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillExtension(kTRUE);
2a1b1442 268 //PostData(0, fConversionGammas);
e5b6e8a6 269
1528f6d1 270 }
271 }
272 }
273}
274
275///________________________________________________________________________
276const AliExternalTrackParam *AliV0ReaderV1::GetExternalTrackParam(AliESDv0 *fCurrentV0,Int_t &tracklabel,Int_t charge){
277
278 // Get External Track Parameter with given charge
279
280 if(!(charge==1||charge==-1)){AliError("Charge not defined");return 0x0;}
281
282 // Check for sign flip
283 if(fCurrentV0){
284 if(!fCurrentV0->GetParamN()||!fCurrentV0->GetParamP())return 0x0;
285 if(!fConversionCuts->GetTrack(fInputEvent,fCurrentV0->GetNindex())||!fConversionCuts->GetTrack(fInputEvent,fCurrentV0->GetPindex()))return 0x0;
286 if((fConversionCuts->GetTrack(fInputEvent,fCurrentV0->GetPindex()))->Charge()==charge){
287 tracklabel=fCurrentV0->GetPindex();
288 return fCurrentV0->GetParamP();}
289 if((fConversionCuts->GetTrack(fInputEvent,fCurrentV0->GetNindex()))->Charge()==charge){
290 tracklabel=fCurrentV0->GetNindex();
291 return fCurrentV0->GetParamN();}
292 }
293 return 0x0;
294}
295
296///________________________________________________________________________
297Bool_t AliV0ReaderV1::ProcessESDV0s()
298{
299 // Process ESD V0s for conversion photon reconstruction
1528f6d1 300 AliESDEvent *fESDEvent=dynamic_cast<AliESDEvent*>(fInputEvent);
301
302 AliKFConversionPhoton *fCurrentMotherKFCandidate=NULL;
e5b6e8a6 303
1528f6d1 304 if(fESDEvent){
305
306 for(Int_t currentV0Index=0;currentV0Index<fESDEvent->GetNumberOfV0s();currentV0Index++){
307 AliESDv0 *fCurrentV0=(AliESDv0*)(fESDEvent->GetV0(currentV0Index));
308 if(!fCurrentV0){
309 printf("Requested V0 does not exist");
ccfa8c0d 310 continue;
311 }
1528f6d1 312
313 fCurrentMotherKFCandidate=ReconstructV0(fCurrentV0,currentV0Index);
314
315 if(fCurrentMotherKFCandidate){
316
ae947965 317 //Bool_t aodV0Found = kFALSE;
318 if(AODEvent() && fCheckAODConsistenty){
319 // for(Int_t i = 0; i<AODEvent()->GetNumberOfV0s();i++){
320 // AliAODv0 *currebtAODV0 = AODEvent()->GetV0(i);
321 // cout<<currebtAODV0->GetID()<<endl;
322 // if(currebtAODV0->GetID() == currentV0Index)
323 // aodV0Found = kTRUE;
324 // }
325 // if(!aodV0Found)
326 // AliError(Form("AODV0 not Found belonging to ESDV0 %i",currentV0Index));
327
328
329 if(!(fConversionCuts->GetTrack(AODEvent(),fCurrentMotherKFCandidate->GetTrackLabelPositive())) || !(fConversionCuts->GetTrack(AODEvent(),fCurrentMotherKFCandidate->GetTrackLabelNegative()))){
330 fConversionGammas->Delete(); // Reset Gamma Array
331 AliError(Form("AOD Tracks not found for current ESD V0!!! V0 index %i, posESDtrack %i negESDtrack %i, Run Number: %i, Period Number: %i, NTracks: ESD %i AOD %i",
332 currentV0Index,
333 fCurrentMotherKFCandidate->GetTrackLabelPositive(),fCurrentMotherKFCandidate->GetTrackLabelNegative(),
334 fInputEvent->GetRunNumber(),fInputEvent->GetPeriodNumber(),
335 fInputEvent->GetNumberOfTracks(),AODEvent()->GetNumberOfTracks()));
336 return kTRUE;
337 }
338 }
1528f6d1 339 // Add Gamma to the TClonesArray
340
341 if(kUseAODConversionPhoton){
342 new((*fConversionGammas)[fConversionGammas->GetEntriesFast()]) AliAODConversionPhoton(fCurrentMotherKFCandidate);
ae947965 343 AliAODConversionPhoton * currentConversionPhoton = (AliAODConversionPhoton*)(fConversionGammas->At(fConversionGammas->GetEntriesFast()-1));
344 currentConversionPhoton->SetMass(fCurrentMotherKFCandidate->M());
345 currentConversionPhoton->SetMassToZero();
346
1528f6d1 347 }
348 else{
349 new((*fConversionGammas)[fConversionGammas->GetEntriesFast()]) AliKFConversionPhoton(*fCurrentMotherKFCandidate);
350 }
351
352 delete fCurrentMotherKFCandidate;
353 fCurrentMotherKFCandidate=NULL;
354 }
355 }
356 }
357 return kTRUE;
358}
359
360///________________________________________________________________________
361AliKFConversionPhoton *AliV0ReaderV1::ReconstructV0(AliESDv0 *fCurrentV0,Int_t currentV0Index)
362{
363 // Reconstruct conversion photon from ESD v0
364 fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kPhotonIn);
365
366 //checks if on the fly mode is set
367 if(!fConversionCuts->SelectV0Finder(fCurrentV0->GetOnFlyStatus())){
e5b6e8a6 368 fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kOnFly);
369 return 0x0;
1528f6d1 370 }
371
372 // TrackLabels
373 Int_t currentTrackLabels[2]={-1,-1};
374
375 // Get Daughter KF Particles
376
377 const AliExternalTrackParam *fCurrentExternalTrackParamPositive=GetExternalTrackParamP(fCurrentV0,currentTrackLabels[0]);
378 const AliExternalTrackParam *fCurrentExternalTrackParamNegative=GetExternalTrackParamN(fCurrentV0,currentTrackLabels[1]);
379
380 if(!fCurrentExternalTrackParamPositive||!fCurrentExternalTrackParamNegative)return 0x0;
381
382 // Apply some Cuts before Reconstruction
383
384 AliVTrack * posTrack = fConversionCuts->GetTrack(fInputEvent,currentTrackLabels[0]);
385 AliVTrack * negTrack = fConversionCuts->GetTrack(fInputEvent,currentTrackLabels[1]);
386
387 if(!negTrack || !posTrack) {
e5b6e8a6 388 fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kNoTracks);
389 return 0x0;
1528f6d1 390 }
391
392 // Track Cuts
393 if(!fConversionCuts->TracksAreSelected(negTrack, posTrack)){
e5b6e8a6 394 fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kTrackCuts);
395 return 0x0;
1528f6d1 396 }
397
398 // PID Cuts
e5b6e8a6 399 if(!fConversionCuts->dEdxCuts(negTrack) || !fConversionCuts->dEdxCuts(posTrack)) {
400 fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kdEdxCuts);
401 return 0x0;
402 }
1528f6d1 403
404 // Reconstruct Photon
405
406 AliKFConversionPhoton *fCurrentMotherKF=NULL;
407
408 AliKFParticle fCurrentNegativeKFParticle(*(fCurrentExternalTrackParamNegative),11);
409 AliKFParticle fCurrentPositiveKFParticle(*(fCurrentExternalTrackParamPositive),-11);
410
411 // Reconstruct Gamma
412
413 if(fUseConstructGamma){
ae947965 414 fCurrentMotherKF = new AliKFConversionPhoton();
415 fCurrentMotherKF->ConstructGamma(fCurrentNegativeKFParticle,fCurrentPositiveKFParticle);
1528f6d1 416 }else{
ae947965 417 fCurrentMotherKF = new AliKFConversionPhoton(fCurrentNegativeKFParticle,fCurrentPositiveKFParticle);
418 fCurrentMotherKF->SetMassConstraint(0,0.0001);
1528f6d1 419 }
420
421 // Set Track Labels
422
423 fCurrentMotherKF->SetTrackLabels(currentTrackLabels[0],currentTrackLabels[1]);
424
425 // Set V0 index
426
427 fCurrentMotherKF->SetV0Index(currentV0Index);
428
429 //Set MC Label
430
431 if(fMCEvent){
432
433 AliStack *fMCStack= fMCEvent->Stack();
434
e5b6e8a6 435 Int_t labelp=TMath::Abs(fConversionCuts->GetTrack(fInputEvent,fCurrentMotherKF->GetTrackLabelPositive())->GetLabel());
436 Int_t labeln=TMath::Abs(fConversionCuts->GetTrack(fInputEvent,fCurrentMotherKF->GetTrackLabelNegative())->GetLabel());
1528f6d1 437
438 TParticle *fNegativeMCParticle = fMCStack->Particle(labeln);
439 TParticle *fPositiveMCParticle = fMCStack->Particle(labelp);
440
441 if(fPositiveMCParticle&&fNegativeMCParticle){
442 fCurrentMotherKF->SetMCLabelPositive(labelp);
443 fCurrentMotherKF->SetMCLabelNegative(labeln);
444 }
445 }
446
447 // Update Vertex (moved for same eta compared to old)
448 if(fUseImprovedVertex == kTRUE){
449 AliKFVertex primaryVertexImproved(*fInputEvent->GetPrimaryVertex());
450 primaryVertexImproved+=*fCurrentMotherKF;
451 fCurrentMotherKF->SetProductionVertex(primaryVertexImproved);
452 }
453
454 // SetPsiPair
455
456 Double_t PsiPair=GetPsiPair(fCurrentV0,fCurrentExternalTrackParamPositive,fCurrentExternalTrackParamNegative);
457 fCurrentMotherKF->SetPsiPair(PsiPair);
458
e5b6e8a6 459
1528f6d1 460 // Recalculate ConversionPoint
461 Double_t dca[2]={0,0};
462 if(fUseOwnXYZCalculation){
463 Double_t convpos[3]={0,0,0};
464 if(!GetConversionPoint(fCurrentExternalTrackParamPositive,fCurrentExternalTrackParamNegative,convpos,dca)){
e5b6e8a6 465 fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kConvPointFail);
466 delete fCurrentMotherKF;
467 fCurrentMotherKF=NULL;
468 return 0x0;
1528f6d1 469 }
470
471 fCurrentMotherKF->SetConversionPoint(convpos);
472 }
473
474 if(fCurrentMotherKF->GetNDF() > 0.)
475 fCurrentMotherKF->SetChi2perNDF(fCurrentMotherKF->GetChi2()/fCurrentMotherKF->GetNDF()); //->Photon is created before all chi2 relevant changes are performed, set it "by hand"
e5b6e8a6 476
477
478 // Set Dilepton Mass (moved down for same eta compared to old)
479 fCurrentMotherKF->SetMass(fCurrentMotherKF->M());
480
1528f6d1 481 // Apply Photon Cuts
482
483 if(!fConversionCuts->PhotonCuts(fCurrentMotherKF,fInputEvent)){
484 fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kPhotonCuts);
485 delete fCurrentMotherKF;
486 fCurrentMotherKF=NULL;
487 return 0x0;
488 }
489
1528f6d1 490 fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kPhotonOut);
491 return fCurrentMotherKF;
492}
493
494///________________________________________________________________________
495Double_t AliV0ReaderV1::GetPsiPair(const AliESDv0* v0, const AliExternalTrackParam *positiveparam,const AliExternalTrackParam *negativeparam) const {
496 //
497 // Angle between daughter momentum plane and plane
498 //
499
a280ac15 500 AliExternalTrackParam nt(*negativeparam);
501 AliExternalTrackParam pt(*positiveparam);
1528f6d1 502
503 Float_t magField = fInputEvent->GetMagneticField();
504
505 Double_t xyz[3] = {0.,0.,0.};
506 v0->GetXYZ(xyz[0],xyz[1],xyz[2]);
e5b6e8a6 507
a280ac15 508 // Double_t pPlus[3] = {pt.Px(),pt.Py(),pt.Pz()};
509 // Double_t pMinus[3] = {nt.Px(),nt.Py(),nt.Pz()};
510
511 // Double_t u[3] = {pPlus[0]+pMinus[0],pPlus[1]+pMinus[1],pPlus[2]+pMinus[2]};
512 // Double_t normu = sqrt( (u[0]*u[0]) + (u[1]*u[1]) + (u[2]*u[2]) );
11c1e680 513
a280ac15 514 // u[0] = u[0] / normu;
515 // u[1] = u[1] / normu;
516 // u[2] = u[2] / normu;
517
518 // Double_t normpPlus = sqrt( (pPlus[0]*pPlus[0]) + (pPlus[1]*pPlus[1]) + (pPlus[2]*pPlus[2]) );
519 // Double_t normpMinus = sqrt( (pMinus[0]*pMinus[0]) + (pMinus[1]*pMinus[1]) + (pMinus[2]*pMinus[2]) );
520
521 // pPlus[0] = pPlus[0] / normpPlus;
522 // pPlus[1] = pPlus[1] / normpPlus;
523 // pPlus[2] = pPlus[2] / normpPlus;
524
525 // pMinus[0] = pMinus[0] / normpMinus;
526 // pMinus[1] = pMinus[1] / normpMinus;
527 // pMinus[2] = pMinus[2] / normpMinus;
528
529 // Double_t v[3] = {0,0,0}; // pPlus X pMinus
530 // v[0] = (pPlus[1]*pMinus[2]) - (pPlus[2]*pMinus[1]);
531 // v[1] = (pPlus[2]*pMinus[0]) - (pPlus[0]*pMinus[2]);
532 // v[2] = (pPlus[0]*pMinus[1]) - (pPlus[1]*pMinus[0]);
11c1e680 533
a280ac15 534 // Double_t w[3] = {0,0,0}; // u X v
535 // w[0] = (u[1]*v[2]) - (u[2]*v[1]);
536 // w[1] = (u[2]*v[0]) - (u[0]*v[2]);
537 // w[2] = (u[0]*v[1]) - (u[1]*v[0]);
538
539 // Double_t z[3] = {0,0,1};
540 // Double_t wc[3] = {0,0,0}; // u X v
541 // wc[0] = (u[1]*z[2]) - (u[2]*z[1]);
542 // wc[1] = (u[2]*z[0]) - (u[0]*z[2]);
543 // wc[2] = (u[0]*z[1]) - (u[1]*z[0]);
544
545 // Double_t PhiV = TMath::ACos((w[0]*wc[0]) + (w[1]*wc[1]) + (w[2]*wc[2]));
546 //return abs(PhiV);
547
548
549 // TVector3 pPlus(pt.Px(),pt.Py(),pt.Pz());
550 // TVector3 pMinus(nt.Px(),nt.Py(),nt.Pz());
551
552 // TVector3 u = pMinus + pPlus;
553 // u = u*(1/u.Mag());
554
555 // TVector3 pHPlus = pPlus*(1/pPlus.Mag());
556 // TVector3 pHMinus = pMinus*(1/pMinus.Mag());
557
558 // TVector3 v = pHPlus.Cross(pHMinus);
559 // TVector3 w = u.Cross(v);
560 // TVector3 z(0,0,1);
561 // TVector3 wc = u.Cross(z);
562
563 // Double_t PhiV = w * wc;
564
1528f6d1 565 Double_t mn[3] = {0,0,0};
566 Double_t mp[3] = {0,0,0};
e5b6e8a6 567
1528f6d1 568 v0->GetNPxPyPz(mn[0],mn[1],mn[2]);//reconstructed cartesian momentum components of negative daughter;
e5b6e8a6 569 v0->GetPPxPyPz(mp[0],mp[1],mp[2]);//reconstructed cartesian momentum components of positive daughter;
1528f6d1 570
571 Double_t deltat = 1.;
a280ac15 572 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
e5b6e8a6 573 Double_t radiussum = TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]) + 50;//radius to which tracks shall be propagated
574
1528f6d1 575 Double_t momPosProp[3] = {0,0,0};
576 Double_t momNegProp[3] = {0,0,0};
e5b6e8a6 577
1528f6d1 578 Double_t psiPair = 4.;
579 if(nt.PropagateTo(radiussum,magField) == 0) return psiPair; //propagate tracks to the outside -> Better Purity and Efficiency
1528f6d1 580 if(pt.PropagateTo(radiussum,magField) == 0) return psiPair; //propagate tracks to the outside -> Better Purity and Efficiency
e5b6e8a6 581
1528f6d1 582 pt.GetPxPyPz(momPosProp);//Get momentum vectors of tracks after propagation
583 nt.GetPxPyPz(momNegProp);
584
585 Double_t pEle =
586 TMath::Sqrt(momNegProp[0]*momNegProp[0]+momNegProp[1]*momNegProp[1]+momNegProp[2]*momNegProp[2]);//absolute momentum value of negative daughter
587
588 Double_t pPos =
589 TMath::Sqrt(momPosProp[0]*momPosProp[0]+momPosProp[1]*momPosProp[1]+momPosProp[2]*momPosProp[2]);//absolute momentum value of positive daughter
590
591 Double_t scalarproduct =
592 momPosProp[0]*momNegProp[0]+momPosProp[1]*momNegProp[1]+momPosProp[2]*momNegProp[2];//scalar product of propagated positive and negative daughters' momenta
593
594 Double_t chipair = TMath::ACos(scalarproduct/(pEle*pPos));//Angle between propagated daughter tracks
595
596 psiPair = TMath::Abs(TMath::ASin(deltat/chipair));
597
598 return psiPair;
599}
600
601///________________________________________________________________________
602Bool_t AliV0ReaderV1::GetHelixCenter(const AliExternalTrackParam *track,Double_t center[2]){
603
604 // Get Center of the helix track parametrization
605
606 Int_t charge=track->Charge();
607 Double_t b=fInputEvent->GetMagneticField();
608
609 Double_t helix[6];
610 track->GetHelixParameters(helix,b);
611
612 Double_t xpos = helix[5];
613 Double_t ypos = helix[0];
614 Double_t radius = TMath::Abs(1./helix[4]);
615 Double_t phi = helix[2];
616
617 if(phi < 0){
618 phi = phi + 2*TMath::Pi();
619 }
620
621 phi -= TMath::Pi()/2.;
622 Double_t xpoint = radius * TMath::Cos(phi);
623 Double_t ypoint = radius * TMath::Sin(phi);
624
625 if(b<0){
626 if(charge > 0){
627 xpoint = - xpoint;
628 ypoint = - ypoint;
629 }
630
631 if(charge < 0){
632 xpoint = xpoint;
633 ypoint = ypoint;
634 }
635 }
636 if(b>0){
637 if(charge > 0){
638 xpoint = xpoint;
639 ypoint = ypoint;
640 }
641
642 if(charge < 0){
643 xpoint = - xpoint;
644 ypoint = - ypoint;
645 }
646 }
647 center[0] = xpos + xpoint;
648 center[1] = ypos + ypoint;
649
650 return 1;
651}
652///________________________________________________________________________
653Bool_t AliV0ReaderV1::GetConversionPoint(const AliExternalTrackParam *pparam,const AliExternalTrackParam *nparam,Double_t convpos[3],Double_t dca[2]){
654
655 // Recalculate Conversion Point
656
657 if(!pparam||!nparam)return kFALSE;
e5b6e8a6 658
1528f6d1 659 Double_t helixcenterpos[2];
660 GetHelixCenter(pparam,helixcenterpos);
661
662 Double_t helixcenterneg[2];
663 GetHelixCenter(nparam,helixcenterneg);
664
665 Double_t helixpos[6];
666 pparam->GetHelixParameters(helixpos,fInputEvent->GetMagneticField());
667 Double_t posradius = TMath::Abs(1./helixpos[4]);
668
669 Double_t helixneg[6];
670 nparam->GetHelixParameters(helixneg,fInputEvent->GetMagneticField());
671 Double_t negradius = TMath::Abs(1./helixneg[4]);
672
673 // Calculate xy-position
674
675 Double_t xpos = helixcenterpos[0];
676 Double_t ypos = helixcenterpos[1];
677 Double_t xneg = helixcenterneg[0];
678 Double_t yneg = helixcenterneg[1];
679
680 convpos[0] = (xpos*negradius + xneg*posradius)/(negradius+posradius);
681 convpos[1] = (ypos*negradius+ yneg*posradius)/(negradius+posradius);
682
683
684 // Calculate Track XY vertex position
685
686 Double_t deltaXPos = convpos[0] - xpos;
687 Double_t deltaYPos = convpos[1] - ypos;
688
689 Double_t deltaXNeg = convpos[0] - xneg;
690 Double_t deltaYNeg = convpos[1] - yneg;
691
692 Double_t alphaPos = TMath::Pi() + TMath::ATan2(-deltaYPos,-deltaXPos);
693 Double_t alphaNeg = TMath::Pi() + TMath::ATan2(-deltaYNeg,-deltaXNeg);
694
695 Double_t vertexXNeg = xneg + TMath::Abs(negradius)*TMath::Cos(alphaNeg);
696 Double_t vertexYNeg = yneg + TMath::Abs(negradius)*TMath::Sin(alphaNeg);
697
698 Double_t vertexXPos = xpos + TMath::Abs(posradius)*TMath::Cos(alphaPos);
699 Double_t vertexYPos = ypos + TMath::Abs(posradius)*TMath::Sin(alphaPos);
700
701 AliExternalTrackParam p(*pparam);
702 AliExternalTrackParam n(*nparam);
703
704 TVector2 vertexPos(vertexXPos,vertexYPos);
705 TVector2 vertexNeg(vertexXNeg,vertexYNeg);
706
707 // Convert to local coordinate system
708 TVector2 vertexPosRot=vertexPos.Rotate(-p.GetAlpha());
709 TVector2 vertexNegRot=vertexNeg.Rotate(-n.GetAlpha());
710
711 // Propagate Track Params to Vertex
712
713 if(!p.PropagateTo(vertexPosRot.X(),fInputEvent->GetMagneticField()))return kFALSE;
714 if(!n.PropagateTo(vertexNegRot.X(),fInputEvent->GetMagneticField()))return kFALSE;
715
716 // Check whether propagation was sucessful
717
718 if(TMath::Abs(vertexPos.Mod()-TMath::Sqrt(p.GetX()*p.GetX()+p.GetY()*p.GetY()))>0.01)return kFALSE;
719 if(TMath::Abs(vertexNeg.Mod()-TMath::Sqrt(n.GetX()*n.GetX()+n.GetY()*n.GetY()))>0.01)return kFALSE;
720
721 // Calculate z position
722
723 convpos[2] = (p.GetZ()*negradius+n.GetZ()*posradius)/(negradius+posradius);
724
725 // Calculate DCA
726 TVector2 vdca=vertexPos-vertexNeg;
727 dca[0]=vdca.Mod();
728 dca[1]=TMath::Abs(n.GetZ()-p.GetZ());
729
730 return kTRUE;
731}
732//________________________________________________________________________
733Bool_t AliV0ReaderV1::GetAODConversionGammas(){
734
735 // Get reconstructed conversion photons from satellite AOD file
736
737 AliAODEvent *fAODEvent=dynamic_cast<AliAODEvent*>(fInputEvent);
738
739 if(fAODEvent){
740
e5b6e8a6 741 if(fConversionGammas == NULL){
742 fConversionGammas = new TClonesArray("AliAODConversionPhoton",100);
743 }
1528f6d1 744 fConversionGammas->Delete();//Reset the TClonesArray
745
746 //Get Gammas from satellite AOD gamma branch
747
748 AliAODConversionPhoton *gamma=0x0;
e5b6e8a6 749
1528f6d1 750 TClonesArray *fInputGammas=dynamic_cast<TClonesArray*>(fAODEvent->FindListObject(fDeltaAODBranchName.Data()));
751 if(!fInputGammas){
e5b6e8a6 752 FindDeltaAODBranchName();
753 fInputGammas=dynamic_cast<TClonesArray*>(fAODEvent->FindListObject(fDeltaAODBranchName.Data()));}
1528f6d1 754 if(!fInputGammas){AliError("No Gamma Satellites found");return kFALSE;}
755 // Apply Selection Cuts to Gammas and create local working copy
756 if(fInputGammas){
e5b6e8a6 757 for(Int_t i=0;i<fInputGammas->GetEntriesFast();i++){
758 gamma=dynamic_cast<AliAODConversionPhoton*>(fInputGammas->At(i));
759 if(gamma){
760 if(fConversionCuts->PhotonIsSelected(gamma,fInputEvent)){
761 new((*fConversionGammas)[fConversionGammas->GetEntriesFast()]) AliAODConversionPhoton(*gamma);}
762 }
763 }
1528f6d1 764 }
765 }
766
767 if(fConversionGammas->GetEntries()){return kTRUE;}
768
769 return kFALSE;
770}
771
772//________________________________________________________________________
773void AliV0ReaderV1::FindDeltaAODBranchName(){
774
775 // Find delta AOD branchname containing reconstructed photons
776
777 TList *list=fInputEvent->GetList();
778 for(Int_t ii=0;ii<list->GetEntries();ii++){
779 TString name((list->At(ii))->GetName());
780 if(name.BeginsWith(fDeltaAODBranchName)&&name.EndsWith("gamma")){
781 fDeltaAODBranchName=name;
782 AliInfo(Form("Set DeltaAOD BranchName to: %s",fDeltaAODBranchName.Data()));
783 return;
784 }
785 }
786}
1528f6d1 787//________________________________________________________________________
788void AliV0ReaderV1::Terminate(Option_t *)
789{
790
791}