]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG3/muon/AliAnalysisTaskMuonTreeBuilder.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
some fixes wrt to the AliAODExtension behavior (Laurent)
[u/mrichter/AliRoot.git] / PWG3 / muon / AliAnalysisTaskMuonTreeBuilder.cxx
CommitLineData
9fe33e0b 1#ifndef ALIANALYSISTASKMUONTREEBUILDER_CXX
2#define ALIANALYSISTASKMUONTREEBUILDER_CXX
3
4#include "AliAnalysisTaskMuonTreeBuilder.h"
5#include "AliMCEvent.h"
6#include "AliESDMuonTrack.h"
7#include "AliESDVertex.h"
8#include "AliStack.h"
9#include "AliHeader.h"
10#include "AliESDHeader.h"
11#include "TParticle.h"
12#include "TLorentzVector.h"
13#include "TFile.h"
14#include "TH1I.h"
15#include "AliAnalysisManager.h"
16#include "AliESDEvent.h"
17#include "AliAODEvent.h"
18#include "TChain.h"
19#include "AliESDtrack.h"
20#include "AliLog.h"
21#include "AliESDtrack.h"
22#include "AliESDInputHandler.h"
23#include "TCanvas.h"
24#include "AliPhysicsSelection.h"
25#include "AliMultiplicity.h"
26
27// Analysis task for muon-dimuon analysis
28// Works for real and MC events
29// Works with the corresponding AddTask macro
c9726e00 30// Includes a flag for physics selection
9fe33e0b 31//
32// author: L. Bianchi - Universita' & INFN Torino
33
34ClassImp(AliAnalysisTaskMuonTreeBuilder)
35
36//__________________________________________________________________________
37AliAnalysisTaskMuonTreeBuilder::AliAnalysisTaskMuonTreeBuilder() :
38 fNevt(0),
39 fBeamEnergy(5000.),
40 fOutput(0x0),
41 fOutputTree(0x0),
42 fIsMC(kFALSE),
43 fIsSelected(kFALSE),
44 fNumMuonTracks(0),
45 fNumSPDTracklets(0),
46 fNumContributors(0),
47 fNumDimuons(0)
48{
49 //
50 //Default ctor
51 //
52}
53//___________________________________________________________________________
54AliAnalysisTaskMuonTreeBuilder::AliAnalysisTaskMuonTreeBuilder(const Char_t* name) :
55 AliAnalysisTaskSE(name),
56 fNevt(0),
57 fBeamEnergy(5000.),
58 fOutput(0x0),
59 fOutputTree(0x0),
60 fIsMC(kFALSE),
61 fIsSelected(kFALSE),
62 fNumMuonTracks(0),
63 fNumSPDTracklets(0),
64 fNumContributors(0),
65 fNumDimuons(0)
66{
67 //
68 // Constructor. Initialization of Inputs and Outputs
69 //
70 Info("AliAnalysisTaskMuonTreeBuilder","Calling Constructor");
71
72 DefineOutput(1,TTree::Class());
73}
74
75//___________________________________________________________________________
76AliAnalysisTaskMuonTreeBuilder& AliAnalysisTaskMuonTreeBuilder::operator=(const AliAnalysisTaskMuonTreeBuilder& c)
77{
78 //
79 // Assignment operator
80 //
81 if (this!=&c) {
82 AliAnalysisTaskSE::operator=(c) ;
83 fNevt = c.fNevt ;
84 }
85 return *this;
86}
87
88//___________________________________________________________________________
89AliAnalysisTaskMuonTreeBuilder::AliAnalysisTaskMuonTreeBuilder(const AliAnalysisTaskMuonTreeBuilder& c) :
90 AliAnalysisTaskSE(c),
91 fNevt(c.fNevt),
92 fBeamEnergy(c.fBeamEnergy),
93 fOutput(c.fOutput),
94 fOutputTree(c.fOutputTree),
95 fIsMC(kFALSE),
96 fIsSelected(kFALSE),
97 fNumMuonTracks(c.fNumMuonTracks),
98 fNumSPDTracklets(c.fNumSPDTracklets),
99 fNumContributors(c.fNumContributors),
100 fNumDimuons(c.fNumDimuons)
101{
102 //
103 // Copy Constructor FIDUCIAL REGIONS?
104 //
105}
106
107//___________________________________________________________________________
108AliAnalysisTaskMuonTreeBuilder::~AliAnalysisTaskMuonTreeBuilder() {
109 //
110 //destructor
111 //
112 Info("~AliAnalysisTaskMuonTreeBuilder","Calling Destructor");
113}
114
115//___________________________________________________________________________
116void AliAnalysisTaskMuonTreeBuilder::UserCreateOutputObjects(){
117
118//
119// Creating User-Defined Output Objects
120//
121
122 // TREE OUTPUT----------------------------------------------------------
123 OpenFile(1);
124 fOutputTree = new TTree("krec","Tree of reconstructed muons");
125
126 fOutputTree->Branch("IsSelected",&fIsSelected,"IsSelected/O");
127 fOutputTree->Branch("FiredTriggerClasses",fTrigClass,"FiredTriggerClasses/C");
128
129 fOutputTree->Branch("NumMuonTracks",&fNumMuonTracks,"NumMuonTracks/I");
130 fOutputTree->Branch("NumContributors",&fNumContributors,"NumContributors/I");
131 fOutputTree->Branch("NumSPDTracklets",&fNumSPDTracklets,"NumSPDTraclets/I");
132 fOutputTree->Branch("Vertex",fVertex,"Vertex[3]/D");
133 fOutputTree->Branch("pT",fpT,"pT[10]/D");
134 fOutputTree->Branch("E",fE,"E[10]/D");
135 fOutputTree->Branch("px",fpx,"px[10]/D");
136 fOutputTree->Branch("py",fpy,"py[10]/D");
137 fOutputTree->Branch("pz",fpz,"pz[10]/D");
138 fOutputTree->Branch("pxUncorr",fpxUncorr,"pxUncorr[10]/D");
139 fOutputTree->Branch("pyUncorr",fpyUncorr,"pyUncorr[10]/D");
140 fOutputTree->Branch("pzUncorr",fpzUncorr,"pzUncorr[10]/D");
141 fOutputTree->Branch("y",fy,"y[10]/D");
142 fOutputTree->Branch("eta",feta,"eta[10]/D");
143 fOutputTree->Branch("phi",fphi,"phi[10]/D");
144 fOutputTree->Branch("MatchTrig",fMatchTrig,"MatchTrig[10]/I");
145 fOutputTree->Branch("TrackChi2",fTrackChi2,"TrackChi2[10]/D");
146 fOutputTree->Branch("MatchTrigChi2",fMatchTrigChi2,"MatchTrigChi2[10]/D");
147 fOutputTree->Branch("DCA",fDCA,"DCA[10]/D");
148 fOutputTree->Branch("Charge",fCharge,"Charge[10]/S");
149 fOutputTree->Branch("MuFamily",fMuFamily,"MuFamily[10]/I");
150 fOutputTree->Branch("RAtAbsEnd",fRAtAbsEnd,"RAtAbsEnd[10]/D");
151
152 fOutputTree->Branch("NumDimuons",&fNumDimuons,"NumDimuons/I");
153 fOutputTree->Branch("DimuConstituent",fDimuonConstituent,"DimuonConstituent[45][2]/I");
154 fOutputTree->Branch("pTdimuon",fpTdimuon,"pTdimuon[45]/D");
155 fOutputTree->Branch("pxdimuon",fpxdimuon,"pxdimuon[45]/D");
156 fOutputTree->Branch("pydimuon",fpydimuon,"pydimuon[45]/D");
157 fOutputTree->Branch("pzdimuon",fpzdimuon,"pzdimuon[45]/D");
158 fOutputTree->Branch("ydimuon",fydimuon,"ydimuon[45]/D");
159 fOutputTree->Branch("Imassdimuon",fiMassdimuon,"iMassdimuon[45]/D");
160 fOutputTree->Branch("costCS",fcostCS,"costCS[45]/D");
161 fOutputTree->Branch("costHE",fcostHE,"costHE[45]/D");
162 fOutputTree->Branch("phiCS",fphiCS,"phiCS[45]/D");
163 fOutputTree->Branch("phiHE",fphiHE,"phiHE[45]/D");
164
165 fOutputTree->Branch("PDG",fPDG,"PDG[10]/I");
166 fOutputTree->Branch("PDGmother",fPDGmother,"PDGmother[10]/I");
167 fOutputTree->Branch("PDGdimu",fPDGdimu,"PDGdimu[45]/I");
168
169}
170
171
172
173//_________________________________________________
174void AliAnalysisTaskMuonTreeBuilder::UserExec(Option_t *)
175{
176//
177// User Exec
178//
179
180 fNevt++;
181
182
183 fNumMuonTracks=0;
184 fNumSPDTracklets=666.;
185 fNumContributors=666.;
186 fNumDimuons=0;
187 fIsSelected=kFALSE;
188 fVertex[0]=666.; fVertex[1]=666.; fVertex[2]=666.;
189 for(Int_t i=0; i<10;i++){
190 fpT[i]=666.;
191 fE[i]=666.;
192 fpx[i]=666;
193 fpy[i]=666;
194 fpz[i]=666;
195 fpxUncorr[i]=666;
196 fpyUncorr[i]=666;
197 fpzUncorr[i]=666;
198 fy[i]=666.;
199 feta[i]=666.;
200 fphi[i]=666.;
201 fMatchTrig[i]=666.;
202 fTrackChi2[i]=666.;
203 fMatchTrigChi2[i]=666.;
204 fDCA[i]=666.;
205 fPDG[i]=666;
206 fPDGmother[i]=666;
207 fCharge[i]=666;
208 fMuFamily[i]=666;
209 fRAtAbsEnd[i]=666;
210 }
211 for(Int_t i=0; i<45;i++){
212 fDimuonConstituent[i][0]=666; fDimuonConstituent[i][1]=666;
213 fpTdimuon[i]=666.;
214 fpxdimuon[i]=666.;
215 fpydimuon[i]=666.;
216 fpzdimuon[i]=666.;
217 fydimuon[i]=666.;
218 fiMassdimuon[i]=666.;
219 fcostCS[i]=666.;
220 fcostHE[i]=666.;
221 fphiCS[i]=666.;
222 fphiHE[i]=666.;
223 fPDGdimu[i]=666;
224 }
225
226
227////////
228//// ESD
229////////
230
231 AliESDEvent *fESD = 0x0;
232 AliMCEvent* mcEvent = 0x0;
233
234 if(fIsMC){
235 if (!fMCEvent) {
236 Error("UserExec","NO MC EVENT FOUND!");
237 return;
238 }
239 }
240
241 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
242
243 fIsSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
244
245 fNumMuonTracks = fESD->GetNumberOfMuonTracks() ;
246 Int_t loopEnd = fNumMuonTracks;
247 if(!fIsMC) {
248 TString cla = fESD->GetFiredTriggerClasses();
249 sprintf(fTrigClass,"%s",cla.Data());
250 }
251
252 if(fNumMuonTracks>0 && fIsMC){
253 mcEvent = MCEvent();
254 }
255 fNumContributors = fESD->GetPrimaryVertexSPD()->GetNContributors();
256 fNumSPDTracklets = fESD->GetMultiplicity()->GetNumberOfTracklets();
257 fVertex[0]=fESD->GetPrimaryVertexSPD()->GetX();
258 fVertex[1]=fESD->GetPrimaryVertexSPD()->GetY();
259 fVertex[2]=fESD->GetPrimaryVertexSPD()->GetZ();
260 printf("fVertex : %f - %f - %f\n",fVertex[0],fVertex[1],fVertex[2]);
261
262 Int_t numdimu = 0;
263 for (Int_t j = 0; j<loopEnd; j++) {
264 AliESDMuonTrack* mu1 = new AliESDMuonTrack(*(fESD->GetMuonTrack(j)));
265 if (!mu1->ContainTrackerData()) {fNumMuonTracks=fNumMuonTracks-1; continue;}
266 Double_t charge1 = mu1->Charge();
267 fCharge[j] = mu1->Charge();
268 fpT[j] = mu1->Pt();
269 fpx[j] = mu1->Px();
270 fpy[j] = mu1->Py();
271 fpz[j] = mu1->Pz();
272 fpxUncorr[j] = mu1->PxUncorrected();
273 fpyUncorr[j] = mu1->PyUncorrected();
274 fpzUncorr[j] = mu1->PzUncorrected();
275 fy[j] = mu1->Y();
276 feta[j] = mu1->Eta();
277 fphi[j] = Phideg(mu1->Phi());
278 Double_t emu1 = mu1->E();
279 fE[j] = emu1;
280 fDCA[j] = mu1->GetDCA();
281 fTrackChi2[j] = mu1->GetChi2();
282 fMatchTrig[j] = mu1->GetMatchTrigger();
283 fMatchTrigChi2[j]= mu1->GetChi2MatchTrigger();
284 fRAtAbsEnd[j]=mu1->GetRAtAbsorberEnd();
285
286 AliMCParticle* mcTrack = 0x0;
287 if(fIsMC){
288 if(mu1->GetLabel()==-1) continue;
289 mcTrack = (AliMCParticle*)mcEvent->GetTrack(mu1->GetLabel());
290 fPDG[j] = mcTrack->PdgCode();
291 if (mcTrack->GetMother()==-1) continue;
292 fPDGmother[j] = ((AliMCParticle*)mcEvent->GetTrack(mcTrack->GetMother()))->PdgCode();
293 if (TMath::Abs(fPDG[j])==13) fMuFamily[j] = FindMuFamily(mcTrack,mcEvent);
294 }
295 for (Int_t jj = j+1; jj<loopEnd; jj++) {
296 AliESDMuonTrack* mu2 = new AliESDMuonTrack(*(fESD->GetMuonTrack(jj)));
297 if (!mu2->ContainTrackerData()) continue;
298
299 Double_t pxmu2 = mu2->Px();
300 Double_t pymu2 = mu2->Py();
301 Double_t pzmu2 = mu2->Pz();
302 Double_t emu2 = mu2->E();
303 //Double_t charge2= mu2->Charge();
304
305 fpTdimuon[numdimu] = TMath::Sqrt((fpx[j]+pxmu2)*(fpx[j]+pxmu2)+(fpy[j]+pymu2)*(fpy[j]+pymu2));
306 fpxdimuon[numdimu] = fpx[j]+pxmu2;
307 fpydimuon[numdimu] = fpy[j]+pymu2;
308 fpzdimuon[numdimu] = fpz[j]+pzmu2;
309 fydimuon[numdimu] = Rap((emu1+emu2),(fpz[j]+pzmu2));
310 fiMassdimuon[numdimu] = Imass(emu1,fpx[j],fpy[j],fpz[j],emu2,pxmu2,pymu2,pzmu2);
311 fcostCS[numdimu]=CostCS(fpx[j],fpy[j],fpz[j],emu1,charge1,pxmu2,pymu2,pzmu2,emu2);
312 fcostHE[numdimu]=CostHE(fpx[j],fpy[j],fpz[j],emu1,charge1,pxmu2,pymu2,pzmu2,emu2);
313 fphiCS[numdimu] = PhiCS(fpx[j],fpy[j],fpz[j],emu1,charge1,pxmu2,pymu2,pzmu2,emu2);
314 fphiHE[numdimu] = PhiHE(fpx[j],fpy[j],fpz[j],emu1,charge1,pxmu2,pymu2,pzmu2,emu2);
315
316 fDimuonConstituent[numdimu][0]=j; fDimuonConstituent[numdimu][1]=jj;
317
318 numdimu++;
319
320 if(fIsMC){
321 if(mu2->GetLabel()==-1) continue;
322 if(TMath::Abs(fPDG[j])==13 && TMath::Abs(((AliMCParticle*)mcEvent->GetTrack(mu2->GetLabel()))->PdgCode())==13) fPDGdimu[numdimu-1]=FindDimuFamily(mcTrack,(AliMCParticle*)mcEvent->GetTrack(mu2->GetLabel()),mcEvent);
323 else fPDGdimu[numdimu-1]=-3;
324 }
325
326 delete mu2;
327 }
328 fNumDimuons=numdimu;
329
330
331 delete mu1;
332 }
333 fOutputTree->Fill();
334 PostData(1,fOutputTree);
335}
336
337//________________________________________________________________________
338Int_t AliAnalysisTaskMuonTreeBuilder::FindDimuFamily(AliMCParticle* mcTrack1,AliMCParticle* mcTrack2, AliMCEvent* mcEvent) const
339{
340 // finds the family of the dimuon (works only if the 2 muons are real muons (not hadrons))
341 Int_t familynumber;
342
343 if(mcTrack1->GetMother()==mcTrack2->GetMother()) familynumber = TMath::Abs(((AliMCParticle*)mcEvent->GetTrack(mcTrack1->GetMother()))->PdgCode());
344 else{
345 Int_t familymu1 = FindMuFamily(mcTrack1,mcEvent);
346 Int_t familymu2 = FindMuFamily(mcTrack2,mcEvent);
347 if(familymu1==5 && familymu2==5) familynumber=5; //bb dimuon
348 else if(familymu1==4 && familymu2==4) familynumber=4; //cc dimuon
349 else if((familymu1==4 && familymu2==5)||(familymu2==4 && familymu1==5)) familynumber=45; //bc dimuon
350 else if (familymu1==-2 || familymu2==-2 || familymu1==-3 || familymu2==-3) familynumber=-2; //hadron dimuon (at least 1 hadron involved)
351 else familynumber=-1;
352 }
353
354 return familynumber;
355}
356
357//________________________________________________________________________
358Int_t AliAnalysisTaskMuonTreeBuilder::FindMuFamily(AliMCParticle* mcTrack, AliMCEvent* mcEvent) const
359{
360 // finds the family of the muon
361 Int_t imother = mcTrack->GetMother();
362 if ( imother<0 ) return -1; // Drell-Yan Muon
363
364 Int_t igrandma = imother;
365
366 AliMCParticle* motherPart = (AliMCParticle*)mcEvent->GetTrack(imother);
367 Int_t motherPdg = motherPart->PdgCode();
368
369 // Track is an heavy flavor muon
370 Int_t absPdg = TMath::Abs(motherPdg);
371 if(absPdg/100==5 || absPdg/1000==5) return 5;
372 if(absPdg/100==4 || absPdg/1000==4){
373 Int_t newMother = -1;
374 igrandma = imother;
375 Int_t absGrandMotherPdg = TMath::Abs(motherPart->PdgCode());
376 while ( absGrandMotherPdg > 10 ) {
377 igrandma = ((AliMCParticle*)mcEvent->GetTrack(igrandma))->GetMother();
378 if( igrandma < 0 ) break;
379 absGrandMotherPdg = TMath::Abs(((AliMCParticle*)mcEvent->GetTrack(igrandma))->PdgCode());
380 }
381
382 if (absGrandMotherPdg==5) newMother = 5; // Charm from beauty
383 else if (absGrandMotherPdg==4) newMother = 4;
384
385 if(newMother<0) {
386 //AliWarning("Mother not correctly found! Set to charm!\n");
387 newMother = 4;
388 }
389
390 return newMother;
391 }
392
393 Int_t nPrimaries = mcEvent->Stack()->GetNprimary();
394 // Track is a bkg. muon
395 if (imother<nPrimaries) {
396 return -2; //is a primary
397 }
398 else {
399 return -3; //is a secondary
400 }
401
402}
403
404//________________________________________________________________________
405Double_t AliAnalysisTaskMuonTreeBuilder::Imass(Double_t e1, Double_t px1, Double_t py1, Double_t pz1,
406 Double_t e2, Double_t px2, Double_t py2, Double_t pz2) const
407{
408// invariant mass calculation
409 Double_t imassrec = TMath::Sqrt((e1+e2)*(e1+e2)-((px1+px2)*(px1+px2)+
410 (py1+py2)*(py1+py2)+(pz1+pz2)*(pz1+pz2)));
411 return imassrec;
412}
413
414//________________________________________________________________________
415Double_t AliAnalysisTaskMuonTreeBuilder::Rap(Double_t e, Double_t pz) const
416{
417// calculate rapidity
418 Double_t rap;
c9726e00 419 if(e>TMath::Abs(pz)){
9fe33e0b 420 rap = 0.5*TMath::Log((e+pz)/(e-pz));
421 return rap;
422 }
423 else{
c9726e00 424 rap = 666.;
9fe33e0b 425 return rap;
426 }
427}
428
429//________________________________________________________________________
430Double_t AliAnalysisTaskMuonTreeBuilder::Phideg(Double_t phi) const
431{
432// calculate Phi in range [-180,180]
433 Double_t phideg;
434
435 phideg = phi-TMath::Pi();
436 phideg = phideg*57.296;
437 return phideg;
438}
439
440//________________________________________________________________________
441Double_t AliAnalysisTaskMuonTreeBuilder::CostCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
442Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
443{
c9726e00 444// Cosine of the theta decay angle (mu+) in the Collins-Soper frame
445
446 TLorentzVector pMu1CM, pMu2CM, pProjCM, pTargCM, pDimuCM; // In the CM. frame
447 TLorentzVector pMu1Dimu, pMu2Dimu, pProjDimu, pTargDimu; // In the dimuon rest frame
9fe33e0b 448 TVector3 beta,zaxisCS;
449 Double_t mp=0.93827231;
450 //
451 // --- Fill the Lorentz vector for projectile and target in the CM frame
452 //
c9726e00 453 pProjCM.SetPxPyPzE(0.,0.,-fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp));
454 pTargCM.SetPxPyPzE(0.,0.,fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp));
9fe33e0b 455 //
456 // --- Get the muons parameters in the CM frame
457 //
c9726e00 458 pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
459 pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
9fe33e0b 460 //
461 // --- Obtain the dimuon parameters in the CM frame
462 //
c9726e00 463 pDimuCM=pMu1CM+pMu2CM;
9fe33e0b 464 //
465 // --- Translate the dimuon parameters in the dimuon rest frame
466 //
c9726e00 467 beta=(-1./pDimuCM.E())*pDimuCM.Vect();
468 if(beta.Mag()>=1) return 666.;
469 pMu1Dimu=pMu1CM;
470 pMu2Dimu=pMu2CM;
471 pProjDimu=pProjCM;
472 pTargDimu=pTargCM;
473 pMu1Dimu.Boost(beta);
474 pMu2Dimu.Boost(beta);
475 pProjDimu.Boost(beta);
476 pTargDimu.Boost(beta);
477
478 //Debugging part -------------------------------------
479 Double_t debugProj[4]={0.,0.,0.,0.};
480 Double_t debugTarg[4]={0.,0.,0.,0.};
481 Double_t debugMu1[4]={0.,0.,0.,0.};
482 Double_t debugMu2[4]={0.,0.,0.,0.};
483 pMu1Dimu.GetXYZT(debugMu1);
484 pMu2Dimu.GetXYZT(debugMu2);
485 pProjDimu.GetXYZT(debugProj);
486 pTargDimu.GetXYZT(debugTarg);
487 if (debugProj[0]!=debugProj[0] ||debugProj[1]!=debugProj[1] || debugProj[2]!=debugProj[2] ||debugProj[3]!=debugProj[3]) return 666;
488 if (debugTarg[0]!=debugTarg[0] ||debugTarg[1]!=debugTarg[1] || debugTarg[2]!=debugTarg[2] ||debugTarg[3]!=debugTarg[3]) return 666;
489 if (debugMu1[0]!=debugMu1[0] ||debugMu1[1]!=debugMu1[1] || debugMu1[2]!=debugMu1[2] ||debugMu1[3]!=debugMu1[3]) return 666;
490 if (debugMu2[0]!=debugMu2[0] ||debugMu2[1]!=debugMu2[1] || debugMu2[2]!=debugMu2[2] ||debugMu2[3]!=debugMu2[3]) return 666;
491 //----------------------------------------------------
492
9fe33e0b 493 // --- Determine the z axis for the CS angle
c9726e00 494 zaxisCS=(((pProjDimu.Vect()).Unit())-((pTargDimu.Vect()).Unit())).Unit();
495
9fe33e0b 496 // --- Determine the CS angle (angle between mu+ and the z axis defined above)
497 Double_t cost;
498
c9726e00 499 if(charge1>0) {cost = zaxisCS.Dot((pMu1Dimu.Vect()).Unit());}
500 else {cost = zaxisCS.Dot((pMu2Dimu.Vect()).Unit());}
501
9fe33e0b 502 return cost;
503}
504
505//________________________________________________________________________
506Double_t AliAnalysisTaskMuonTreeBuilder::CostHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
507Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
508{
c9726e00 509// Cosine of the theta decay angle (mu+) in the Helicity frame
510
511 TLorentzVector pMu1CM, pMu2CM, pDimuCM; // In the CM frame
512 TLorentzVector pMu1Dimu, pMu2Dimu; // In the dimuon rest frame
9fe33e0b 513 TVector3 beta,zaxisCS;
514 //
515 // --- Get the muons parameters in the CM frame
516 //
c9726e00 517 pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
518 pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
9fe33e0b 519 //
520 // --- Obtain the dimuon parameters in the CM frame
521 //
c9726e00 522 pDimuCM=pMu1CM+pMu2CM;
9fe33e0b 523 //
524 // --- Translate the muon parameters in the dimuon rest frame
525 //
c9726e00 526 beta=(-1./pDimuCM.E())*pDimuCM.Vect();
527 if(beta.Mag()>=1) return 666.;
528 pMu1Dimu=pMu1CM;
529 pMu2Dimu=pMu2CM;
530 pMu1Dimu.Boost(beta);
531 pMu2Dimu.Boost(beta);
532
533 //Debugging part -------------------------------------
534 Double_t debugMu1[4]={0.,0.,0.,0.};
535 Double_t debugMu2[4]={0.,0.,0.,0.};
536 pMu1Dimu.GetXYZT(debugMu1);
537 pMu2Dimu.GetXYZT(debugMu2);
538 if (debugMu1[0]!=debugMu1[0] ||debugMu1[1]!=debugMu1[1] || debugMu1[2]!=debugMu1[2] ||debugMu1[3]!=debugMu1[3]) return 666;
539 if (debugMu2[0]!=debugMu2[0] ||debugMu2[1]!=debugMu2[1] || debugMu2[2]!=debugMu2[2] ||debugMu2[3]!=debugMu2[3]) return 666;
540 //----------------------------------------------------
541
9fe33e0b 542 // --- Determine the z axis for the calculation of the polarization angle (i.e. the direction of the dimuon in the CM system)
9fe33e0b 543 TVector3 zaxis;
c9726e00 544 zaxis=(pDimuCM.Vect()).Unit();
9fe33e0b 545
546 // --- Calculation of the polarization angle (angle between mu+ and the z axis defined above)
547 Double_t cost;
c9726e00 548 if(charge1>0) {cost = zaxis.Dot((pMu1Dimu.Vect()).Unit());}
549 else {cost = zaxis.Dot((pMu2Dimu.Vect()).Unit());}
9fe33e0b 550 return cost;
551}
552
553//________________________________________________________________________
554Double_t AliAnalysisTaskMuonTreeBuilder::PhiCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
555Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
556{
c9726e00 557// Phi decay angle (mu+) in the Collins-Soper frame
558
559 TLorentzVector pMu1CM, pMu2CM, pProjCM, pTargCM, pDimuCM; // In the CM frame
560 TLorentzVector pMu1Dimu, pMu2Dimu, pProjDimu, pTargDimu; // In the dimuon rest frame
9fe33e0b 561 TVector3 beta,yaxisCS, xaxisCS, zaxisCS;
562 Double_t mp=0.93827231;
c9726e00 563
9fe33e0b 564 // --- Fill the Lorentz vector for projectile and target in the CM frame
c9726e00 565 pProjCM.SetPxPyPzE(0.,0.,-fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp));
566 pTargCM.SetPxPyPzE(0.,0.,fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp));
567
9fe33e0b 568 // --- Get the muons parameters in the CM frame
c9726e00 569 pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
570 pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
571
9fe33e0b 572 // --- Obtain the dimuon parameters in the CM frame
c9726e00 573 pDimuCM=pMu1CM+pMu2CM;
574
9fe33e0b 575 // --- Translate the dimuon parameters in the dimuon rest frame
c9726e00 576 beta=(-1./pDimuCM.E())*pDimuCM.Vect();
577 if(beta.Mag()>=1) return 666.;
578 pMu1Dimu=pMu1CM;
579 pMu2Dimu=pMu2CM;
580 pProjDimu=pProjCM;
581 pTargDimu=pTargCM;
582 pMu1Dimu.Boost(beta);
583 pMu2Dimu.Boost(beta);
584 pProjDimu.Boost(beta);
585 pTargDimu.Boost(beta);
586
587 //Debugging part -------------------------------------
588 Double_t debugProj[4]={0.,0.,0.,0.};
589 Double_t debugTarg[4]={0.,0.,0.,0.};
590 Double_t debugMu1[4]={0.,0.,0.,0.};
591 Double_t debugMu2[4]={0.,0.,0.,0.};
592 pMu1Dimu.GetXYZT(debugMu1);
593 pMu2Dimu.GetXYZT(debugMu2);
594 pProjDimu.GetXYZT(debugProj);
595 pTargDimu.GetXYZT(debugTarg);
596 if (debugProj[0]!=debugProj[0] ||debugProj[1]!=debugProj[1] || debugProj[2]!=debugProj[2] ||debugProj[3]!=debugProj[3]) return 666;
597 if (debugTarg[0]!=debugTarg[0] ||debugTarg[1]!=debugTarg[1] || debugTarg[2]!=debugTarg[2] ||debugTarg[3]!=debugTarg[3]) return 666;
598 if (debugMu1[0]!=debugMu1[0] ||debugMu1[1]!=debugMu1[1] || debugMu1[2]!=debugMu1[2] ||debugMu1[3]!=debugMu1[3]) return 666;
599 if (debugMu2[0]!=debugMu2[0] ||debugMu2[1]!=debugMu2[1] || debugMu2[2]!=debugMu2[2] ||debugMu2[3]!=debugMu2[3]) return 666;
600 //----------------------------------------------------
601
9fe33e0b 602 // --- Determine the z axis for the CS angle
c9726e00 603 zaxisCS=(((pProjDimu.Vect()).Unit())-((pTargDimu.Vect()).Unit())).Unit();
604 yaxisCS=(((pProjDimu.Vect()).Unit()).Cross((pTargDimu.Vect()).Unit())).Unit();
9fe33e0b 605 xaxisCS=(yaxisCS.Cross(zaxisCS)).Unit();
606
607 Double_t phi=0.;
608 if(charge1>0) {
c9726e00 609 phi = TMath::ATan2((pMu1Dimu.Vect()).Dot(yaxisCS),((pMu1Dimu.Vect()).Dot(xaxisCS)));
9fe33e0b 610 } else {
c9726e00 611 phi = TMath::ATan2((pMu2Dimu.Vect()).Dot(yaxisCS),((pMu2Dimu.Vect()).Dot(xaxisCS)));
9fe33e0b 612 }
613 if (phi>TMath::Pi()) phi=phi-TMath::Pi();
c9726e00 614
9fe33e0b 615 return phi;
616}
617
618//________________________________________________________________________
619Double_t AliAnalysisTaskMuonTreeBuilder::PhiHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
620Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
621{
c9726e00 622// Phi decay angle (mu+) in the Helicity frame
623 TLorentzVector pMu1Lab, pMu2Lab, pProjLab, pTargLab, pDimuLab; // In the lab. frame
624 TLorentzVector pMu1Dimu, pMu2Dimu, pProjDimu, pTargDimu; // In the dimuon rest frame
9fe33e0b 625 TVector3 beta,xaxis, yaxis,zaxis;
626 Double_t mp=0.93827231;
627
9fe33e0b 628 // --- Get the muons parameters in the LAB frame
c9726e00 629 pMu1Lab.SetPxPyPzE(px1,py1,pz1,e1);
630 pMu2Lab.SetPxPyPzE(px2,py2,pz2,e2);
631
9fe33e0b 632 // --- Obtain the dimuon parameters in the LAB frame
c9726e00 633 pDimuLab=pMu1Lab+pMu2Lab;
634 zaxis=(pDimuLab.Vect()).Unit();
9fe33e0b 635
9fe33e0b 636 // --- Translate the muon parameters in the dimuon rest frame
c9726e00 637 beta=(-1./pDimuLab.E())*pDimuLab.Vect();
638 if(beta.Mag()>=1.) return 666.;
9fe33e0b 639
c9726e00 640 pProjLab.SetPxPyPzE(0.,0.,-fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp)); // proiettile
641 pTargLab.SetPxPyPzE(0.,0.,fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp)); // bersaglio
9fe33e0b 642
c9726e00 643 pProjDimu=pProjLab;
644 pTargDimu=pTargLab;
9fe33e0b 645
c9726e00 646 pProjDimu.Boost(beta);
647 pTargDimu.Boost(beta);
9fe33e0b 648
c9726e00 649 yaxis=((pProjDimu.Vect()).Cross(pTargDimu.Vect())).Unit();
9fe33e0b 650 xaxis=(yaxis.Cross(zaxis)).Unit();
651
c9726e00 652 pMu1Dimu=pMu1Lab;
653 pMu2Dimu=pMu2Lab;
654 pMu1Dimu.Boost(beta);
655 pMu2Dimu.Boost(beta);
656
657 //Debugging part -------------------------------------
658 Double_t debugProj[4]={0.,0.,0.,0.};
659 Double_t debugTarg[4]={0.,0.,0.,0.};
660 Double_t debugMu1[4]={0.,0.,0.,0.};
661 Double_t debugMu2[4]={0.,0.,0.,0.};
662 pMu1Dimu.GetXYZT(debugMu1);
663 pMu2Dimu.GetXYZT(debugMu2);
664 pProjDimu.GetXYZT(debugProj);
665 pTargDimu.GetXYZT(debugTarg);
666 if (debugProj[0]!=debugProj[0] ||debugProj[1]!=debugProj[1] || debugProj[2]!=debugProj[2] ||debugProj[3]!=debugProj[3]) return 666;
667 if (debugTarg[0]!=debugTarg[0] ||debugTarg[1]!=debugTarg[1] || debugTarg[2]!=debugTarg[2] ||debugTarg[3]!=debugTarg[3]) return 666;
668 if (debugMu1[0]!=debugMu1[0] ||debugMu1[1]!=debugMu1[1] || debugMu1[2]!=debugMu1[2] ||debugMu1[3]!=debugMu1[3]) return 666;
669 if (debugMu2[0]!=debugMu2[0] ||debugMu2[1]!=debugMu2[1] || debugMu2[2]!=debugMu2[2] ||debugMu2[3]!=debugMu2[3]) return 666;
670 //----------------------------------------------------
9fe33e0b 671
672 Double_t phi=0.;
673 if(charge1 > 0) {
c9726e00 674 phi = TMath::ATan2((pMu1Dimu.Vect()).Dot(yaxis),(pMu1Dimu.Vect()).Dot(xaxis));
9fe33e0b 675 } else {
c9726e00 676 phi = TMath::ATan2((pMu2Dimu.Vect()).Dot(yaxis),(pMu2Dimu.Vect()).Dot(xaxis));
9fe33e0b 677 }
678 return phi;
679}
680
681//________________________________________________________________________
682void AliAnalysisTaskMuonTreeBuilder::Terminate(Option_t *)
683{
684// Terminate
685}
686
687#endif