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