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