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