]>
Commit | Line | Data |
---|---|---|
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 | ||
34 | ClassImp(AliAnalysisTaskMuonTreeBuilder) | |
35 | ||
36 | //__________________________________________________________________________ | |
37 | AliAnalysisTaskMuonTreeBuilder::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 | //___________________________________________________________________________ | |
54 | AliAnalysisTaskMuonTreeBuilder::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 | //___________________________________________________________________________ | |
76 | AliAnalysisTaskMuonTreeBuilder& 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 | //___________________________________________________________________________ | |
89 | AliAnalysisTaskMuonTreeBuilder::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 | //___________________________________________________________________________ | |
108 | AliAnalysisTaskMuonTreeBuilder::~AliAnalysisTaskMuonTreeBuilder() { | |
109 | // | |
110 | //destructor | |
111 | // | |
112 | Info("~AliAnalysisTaskMuonTreeBuilder","Calling Destructor"); | |
113 | } | |
114 | ||
115 | //___________________________________________________________________________ | |
116 | void 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 | //_________________________________________________ | |
174 | void 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 | //________________________________________________________________________ | |
338 | Int_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 | //________________________________________________________________________ | |
358 | Int_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 | //________________________________________________________________________ | |
405 | Double_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 | //________________________________________________________________________ | |
415 | Double_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 | //________________________________________________________________________ | |
430 | Double_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 | //________________________________________________________________________ | |
441 | Double_t AliAnalysisTaskMuonTreeBuilder::CostCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1, | |
442 | Double_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 | //________________________________________________________________________ | |
506 | Double_t AliAnalysisTaskMuonTreeBuilder::CostHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1, | |
507 | Double_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 | //________________________________________________________________________ | |
554 | Double_t AliAnalysisTaskMuonTreeBuilder::PhiCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1, | |
555 | Double_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 | //________________________________________________________________________ | |
619 | Double_t AliAnalysisTaskMuonTreeBuilder::PhiHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1, | |
620 | Double_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 | //________________________________________________________________________ | |
682 | void AliAnalysisTaskMuonTreeBuilder::Terminate(Option_t *) | |
683 | { | |
684 | // Terminate | |
685 | } | |
686 | ||
687 | #endif |