]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG3/muon/AliAnalysisTaskDimuonCFContainerBuilder.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Follow PB requested that, for the time being, TRD info is not used to update
[u/mrichter/AliRoot.git] / PWG3 / muon / AliAnalysisTaskDimuonCFContainerBuilder.cxx
CommitLineData
7dd391ec 1#ifndef ALIANALYSISTASKDIMUONCFCONTAINERBUILDER_CXX
2#define ALIANALYSISTASKDIMUONCFCONTAINERBUILDER_CXX
3
4#include "AliAnalysisTaskDimuonCFContainerBuilder.h"
7dd391ec 5#include "AliStack.h"
6#include "TParticle.h"
7#include "TString.h"
8#include "TLorentzVector.h"
7dd391ec 9#include "TH2D.h"
7dd391ec 10#include "AliAnalysisManager.h"
11#include "AliESDEvent.h"
12#include "AliAODEvent.h"
13#include "AliCFManager.h"
7dd391ec 14#include "AliCFContainer.h"
7dd391ec 15#include "AliESDMuonTrack.h"
7dd391ec 16#include "AliESDInputHandler.h"
7dd391ec 17#include "AliAODInputHandler.h"
18#include "AliAODMCParticle.h"
19
20// Analysis task for the building of a dimuon CF container
21// Also some single-muon variables are stored
8320cad2 22// All the variable ranges and binning are set in the AddTask macro
23//
24// author: L. Bianchi - Universita' & INFN Torino
7dd391ec 25
26
27ClassImp(AliAnalysisTaskDimuonCFContainerBuilder)
28
29//__________________________________________________________________________
30AliAnalysisTaskDimuonCFContainerBuilder::AliAnalysisTaskDimuonCFContainerBuilder() :
31 fReadAODData(0),
32 fReadMCInfo(kTRUE),
33 fIsAccProduction(kTRUE),
34 fCFManager(0x0),
35 fQAHistList(0x0),
36 fNevt(0),
37 fBeamEnergy(3500.),
38 fOutput(0x0),
39 fCutOnzVtxSPD(kFALSE),
40 fCutOnNContributors(kFALSE),
41 fTrigClassMuon(""),
42 fTrigClassInteraction(""),
43 fDistinguishTrigClass(kFALSE)
44{
45
46 //Default constructor
47
48 Double_t chilims[2]={0.,10000.};
49 Double_t ptlims[2]={0.,100.};
50 Double_t thetalims[2]={0.,180.};
51 Double_t vtxlims[2]={-1000.,1000.};
52 SetChi2Limits(chilims);
53 SetChi2MatchLimits(chilims);
54 SetPtSingMuLimits(ptlims);
55 SetThetaSingMuLimits(thetalims);
56 SetZprimVertLimits(vtxlims);
57 SetTrigClassMuonName();
58 SetTrigClassInteracName();
59 TString namemuonside[4]={"CMUS1A-","CMUS1B-","CMUS1C-","CMUS1-E-"};
60 TString nameinteractionside[4]={"CINT1A-","CINT1B-","CINT1C-","CINT1-E"};
61 SetTrigClassMuonSideName(namemuonside);
62 SetTrigClassInteracSideName(nameinteractionside);
63}
64//___________________________________________________________________________
65AliAnalysisTaskDimuonCFContainerBuilder::AliAnalysisTaskDimuonCFContainerBuilder(const Char_t* name, Bool_t readaod, Bool_t readMC, Bool_t isaccept, Double_t beamEn) :
66 AliAnalysisTaskSE(name),
67 fReadAODData(0),
68 fReadMCInfo(kTRUE),
69 fIsAccProduction(kTRUE),
70 fCFManager(0x0),
71 fQAHistList(0x0),
72 fNevt(0),
73 fBeamEnergy(3500.),
74 fOutput(0x0),
75 fCutOnzVtxSPD(kFALSE),
76 fCutOnNContributors(kFALSE),
77 fTrigClassMuon(""),
78 fTrigClassInteraction(""),
79 fDistinguishTrigClass(kFALSE)
80{
81 //
82 // Constructor. Initialization of Inputs and Outputs
83 //
84 Info("AliAnalysisTaskDimuonCFContainerBuilder","Calling Constructor");
85
86 SetReadAODData(readaod);
87 SetReadMCinfo(readMC);
88 SetIsAccProd(isaccept);
89 SetBeamEnergy(beamEn);
90
91 Double_t chilims[2]={0.,10000.};
92 Double_t ptlims[2]={0.,100.};
93 Double_t thetalims[2]={0.,180.};
94 Double_t vtxlims[2]={-1000.,1000.};
95 SetChi2Limits(chilims);
96 SetChi2MatchLimits(chilims);
97 SetPtSingMuLimits(ptlims);
98 SetThetaSingMuLimits(thetalims);
99 SetZprimVertLimits(vtxlims);
100 SetTrigClassMuonName();
101 SetTrigClassInteracName();
102 TString namemuonside[4]={"CMUS1A-","CMUS1B-","CMUS1C-","CMUS1-E-"};
103 TString nameinteractionside[4]={"CINT1A-","CINT1B-","CINT1C-","CINT1-E"};
104 SetTrigClassMuonSideName(namemuonside);
105 SetTrigClassInteracSideName(nameinteractionside);
106
107 DefineOutput(1,TList::Class());
108 DefineOutput(2,AliCFContainer::Class());
109
110}
111
112//___________________________________________________________________________
113AliAnalysisTaskDimuonCFContainerBuilder& AliAnalysisTaskDimuonCFContainerBuilder::operator=(const AliAnalysisTaskDimuonCFContainerBuilder& c)
114{
115 //
116 // Assignment operator
117 //
118 if (this!=&c) {
119 AliAnalysisTaskSE::operator=(c) ;
120 fReadAODData = c.fReadAODData ;
121 fCFManager = c.fCFManager;
122 fQAHistList = c.fQAHistList ;
123 fNevt = c.fNevt ;
124 }
125 return *this;
126}
127
128//___________________________________________________________________________
129AliAnalysisTaskDimuonCFContainerBuilder::AliAnalysisTaskDimuonCFContainerBuilder(const AliAnalysisTaskDimuonCFContainerBuilder& c) :
130 AliAnalysisTaskSE(c),
131 fReadAODData(c.fReadAODData),
132 fReadMCInfo(c.fReadMCInfo),
133 fIsAccProduction(c.fIsAccProduction),
134 fCFManager(c.fCFManager),
135 fQAHistList(c.fQAHistList),
136 fNevt(c.fNevt),
137 fBeamEnergy(c.fBeamEnergy),
138 fOutput(c.fOutput),
139 fCutOnzVtxSPD(c.fCutOnzVtxSPD),
140 fCutOnNContributors(c.fCutOnNContributors),
141 fTrigClassMuon(c.fTrigClassMuon),
142 fTrigClassInteraction(c.fTrigClassInteraction),
143 fDistinguishTrigClass(c.fDistinguishTrigClass)
144{
145
146 // Copy Constructor
147
148}
149
150//___________________________________________________________________________
151AliAnalysisTaskDimuonCFContainerBuilder::~AliAnalysisTaskDimuonCFContainerBuilder() {
152 //
153 //destructor
154 //
155 Info("~AliAnalysisTaskDimuonCFContainerBuilder","Calling Destructor");
156 if (fCFManager) delete fCFManager ;
157 if (fQAHistList) {fQAHistList->Clear(); delete fQAHistList;}
158}
159
160//___________________________________________________________________________
161void AliAnalysisTaskDimuonCFContainerBuilder::UserCreateOutputObjects(){
8320cad2 162 // UserCreateOutputObjects
7dd391ec 163 fOutput = new TList();
164 fOutput->SetOwner();
165
166 TH1D *hnevts = new TH1D("hnevts","hnevts",1,0,1); // Stat check histos
8320cad2 167 TH1D *jpsiMult = new TH1D("jpsiMult","jpsiMult",20,0,20);
7dd391ec 168
169 TH1D *ptdimuREC = new TH1D("ptdimuREC","ptdimuREC",200,0,20); // Dimu check histos
170 TH1D *ydimuREC = new TH1D("ydimuREC","ydimuREC",200,-10.,10.);
171 TH1D *ydimuFail = new TH1D("ydimuFail","ydimuFail",10,660,670);
172 TH1D *costHEdimuREC= new TH1D("costHEdimuREC","costHEdimuREC",200,-1.,1.);
173 TH1D *costCSdimuREC= new TH1D("costCSdimuREC","costCSdimuREC",200,-1.,1.);
174 TH1D *costHEdimuFail= new TH1D("costHEdimuFail","costHEdimuFail",10,660.,670.);
175 TH1D *costCSdimuFail= new TH1D("costCSdimuFail","costCSdimuFail",10,660.,670.);
176 TH1D *phiHEdimuREC = new TH1D("phiHEdimuREC","phiHEdimuREC",100,0.,TMath::Pi());
177 TH1D *phiCSdimuREC = new TH1D("phiCSdimuREC","phiCSdimuREC",100,0.,TMath::Pi());
178 TH1D *phiHEdimuFail = new TH1D("phiHEdimuFail","phiHEdimuFail",10,660.,670.);
179 TH1D *phiCSdimuFail = new TH1D("phiCSdimuFail","phiCSdimuFail",10,660.,670.);
180 TH1D *imassTot = new TH1D("imassTot","imassTot",100,0,8);
181 TH1D *trigCond = new TH1D("trigCond","trigCond",40,0,4);
182
8320cad2 183 TH1D *emuonREC = new TH1D("emuonREC","emuonREC",200,0.,20.); // Mu check histos
7dd391ec 184 TH1D *ptmuonREC= new TH1D("ptmuonREC","ptmuonREC",200,0.,20.);
185 TH1D *ymuonREC = new TH1D("ymuonREC","ymuonREC",200,-10.,10.);
186 TH1D *hdca = new TH1D("hdca","hdca",20,0,200);
187 TH2D *hdcay = new TH2D("hdcay","hdcay",200,0,200,20,-5,0);
188
189 TH1D *zvSPD = new TH1D("zvSPD","zvSPD",100,-50,50.); // Event check histos
190 TH1D *zvSPDcut = new TH1D("zvSPDcut","zvSPDcut",100,-50,50.);
8320cad2 191 TH1D *nContrib = new TH1D("nContrib","nContrib",100,0,100);
7dd391ec 192
193
194 fOutput->Add(hnevts);
8320cad2 195 fOutput->Add(jpsiMult);
7dd391ec 196
197 fOutput->Add(ptdimuREC);
198 fOutput->Add(ydimuREC);
199 fOutput->Add(ydimuFail);
200 fOutput->Add(costHEdimuREC);
201 fOutput->Add(costCSdimuREC);
202 fOutput->Add(costHEdimuFail);
203 fOutput->Add(costCSdimuFail);
204 fOutput->Add(phiHEdimuREC);
205 fOutput->Add(phiCSdimuREC);
206 fOutput->Add(phiHEdimuFail);
207 fOutput->Add(phiCSdimuFail);
208 fOutput->Add(imassTot);
209 fOutput->Add(trigCond);
210
8320cad2 211 fOutput->Add(emuonREC);
7dd391ec 212 fOutput->Add(ptmuonREC);
213 fOutput->Add(ymuonREC);
214 fOutput->Add(hdca);
215 fOutput->Add(hdcay);
216
217 fOutput->Add(zvSPD);
218 fOutput->Add(zvSPDcut);
8320cad2 219 fOutput->Add(nContrib);
7dd391ec 220
221
222 fOutput->ls();
223
224
225}
226
227
228//_________________________________________________
229void AliAnalysisTaskDimuonCFContainerBuilder::UserExec(Option_t *)
230{
231
232
233 //Info("UserExec","");
234 fNevt++;
235 ((TH1D*)(fOutput->FindObject("hnevts")))->Fill(0.5);
236
237 Double_t containerInput[15]={666,666,666,666,666,666,666,666,666,666,666,666,666,666,666}; //y, pT, costHE, phiHE, costCS, phiCS, mass, TrigCond
238 Int_t cutAccept =1;
8320cad2 239 Int_t numbJpsis =0;
7dd391ec 240
241 if (!fReadAODData){ // ESD-based ANALYSIS
242
243 if (fReadMCInfo){
244
245 if (!fMCEvent) {
246 Error("UserExec","NO MC EVENT FOUND!");
247 //return;
248 }
249
250 // MC part ---------------------------------------------------------------------------------------------------
251
252 //fCFManager->SetEventInfo(fMCEvent); CHANGES IN NEW VERSIONS - MANUAL CUT ON PDG
253 AliStack* stack = fMCEvent->Stack();
254
255 // loop on the MC event
256 for (Int_t ipart=0; ipart<fMCEvent->GetNumberOfTracks(); ipart++) {
257
258 AliMCParticle *mcPart = (AliMCParticle*) fMCEvent->GetTrack(ipart);
259
260 TParticle *part = mcPart->Particle();
7dd391ec 261
262 // Mother kinematics
263 Double_t e = part->Energy();
264 Double_t pz = part->Pz();
265 Double_t rapmc = Rap(e,pz);
266
267 // Selection of the resonance
268 //if (!fCFManager->CheckParticleCuts(0,mcPart)) continue;
269 if (part->GetPdgCode()!=443) continue; // MANUAL CUT ON PDG
8320cad2 270 numbJpsis++;
7dd391ec 271
272 // Decays kinematics
273 Int_t p0 = part->GetDaughter(0);
b0c84a07 274 TParticle *part0 = stack->Particle(p0);
7dd391ec 275 Int_t pdg0 = part0->GetPdgCode();
276
277 Int_t p1 = part->GetDaughter(1);
b0c84a07 278 TParticle *part1 = stack->Particle(p1);
7dd391ec 279 Int_t pdg1 = part1->GetPdgCode();
280
281 Double_t e0 = part0->Energy();
282 Double_t pz0 = part0->Pz();
283 Double_t py0 = part0->Py();
284 Double_t px0 = part0->Px();
285 Double_t charge0 = part0->GetPDG()->Charge()/3;
286 Double_t theta0 = (180./TMath::Pi())*TMath::ATan2(TMath::Sqrt(px0*px0+py0*py0),pz0);
287 Double_t pt0 = TMath::Sqrt(px0*px0+py0*py0);
8320cad2 288 Double_t mom0 = part0->P();
7dd391ec 289
290 Double_t e1 = part1->Energy();
291 Double_t pz1 = part1->Pz();
292 Double_t py1 = part1->Py();
293 Double_t px1 = part1->Px();
294 //Double_t charge1 = part1->GetPDG()->Charge()/3;
295 Double_t theta1 = (180./TMath::Pi())*TMath::ATan2(TMath::Sqrt(px1*px1+py1*py1),pz1);
296 Double_t pt1 = TMath::Sqrt(px1*px1+py1*py1);
8320cad2 297 Double_t mom1 = part1->P();
7dd391ec 298
299
300 if(pdg0==13 || pdg1==13) {
301
302 Double_t ptmc = TMath::Sqrt((px0+px1)*(px0+px1)+(py0+py1)*(py0+py1));
303 Double_t imassmc = Imass(e0,px0,py0,pz0,e1,px1,py1,pz1);
304
305 Double_t costCS = CostCS(px0,py0,pz0,e0,charge0,px1,py1,pz1,e1);
306 Double_t costHE = CostHE(px0,py0,pz0,e0,charge0,px1,py1,pz1,e1);
307 Double_t phiCS = PhiCS(px0,py0,pz0,e0,charge0,px1,py1,pz1,e1);
308 Double_t phiHE = PhiHE(px0,py0,pz0,e0,charge0,px1,py1,pz1,e1);
309
310 containerInput[0] = rapmc ;
311 containerInput[1] = ptmc;
312 containerInput[2] = costHE;
313 containerInput[3] = TMath::Abs(phiHE);
314 containerInput[4] = costCS;
315 containerInput[5] = TMath::Abs(phiCS);
316 containerInput[6] = imassmc;
317 containerInput[7] = 1.; //for generated no trigger condition
318 if (pt0<pt1) {
319 containerInput[8]=pt0;
320 containerInput[9]=pt1;
321 } else {
322 containerInput[8]=pt1;
323 containerInput[9]=pt0;
324 }
325 if (theta0<theta1) {
326 containerInput[10]=theta0;
327 containerInput[11]=theta1;
328 } else {
329 containerInput[10]=theta1;
330 containerInput[11]=theta0;
331 }
8320cad2 332 if (mom0<mom1) {
333 containerInput[12]=mom0;
334 containerInput[13]=mom1;
7dd391ec 335 } else {
8320cad2 336 containerInput[12]=mom1;
337 containerInput[13]=mom0;
7dd391ec 338 }
339 containerInput[14]=1.;
340
341 // fill the container at the first step
342 fCFManager->GetParticleContainer()->Fill(containerInput,0);
343
344 if(fIsAccProduction){
345 // acceptance cuts on single mu
346 if(theta0<fThetaSingMuCut[0]||theta0>fThetaSingMuCut[1]||theta1<fThetaSingMuCut[0]||theta1>fThetaSingMuCut[1]) cutAccept=0;
347 if(pt0<fPtSingMuCut[0] || pt0>fPtSingMuCut[1] || pt1<fPtSingMuCut[0] || pt1>fPtSingMuCut[1]) cutAccept=0;
348 if (cutAccept ==1) fCFManager->GetParticleContainer()->Fill(containerInput,2);
349 }
350 }
351 }
352 } //fReadMCInfo
353
354
355 // ESD part ---------------------------------------------------------------------------------------------------
356
357 AliESDEvent *fESD;
358 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
359 fESD = esdH->GetEvent();
360 Int_t mult1 = fESD->GetNumberOfMuonTracks() ;
361
362 Int_t trigfired=-1;
363 Int_t trigside=-1;
364 if(fDistinguishTrigClass){
365 //aod->GetHeader()->SetFiredTriggerClasses("CMU");
366 TString trigclass = fESD->GetFiredTriggerClasses();
367 printf("TrigClass: %s\n",trigclass.Data());
368 if(trigclass.Contains(fTrigClassMuon)) trigfired = 1;
369 else if (trigclass.Contains(fTrigClassInteraction)) trigfired = 0;
370 printf("Muon container: %d\n",trigfired);
371 for(Int_t i=0;i<4;i++){
372 if(trigfired==1 && trigclass.Contains(fTrigClassMuonSide[i])) {trigside=i;printf("Muon Side: %d\n",trigside);}
373 if(trigfired==0 && trigclass.Contains(fTrigClassInteractionSide[i])) {trigside=i;printf("Interaction Side: %d\n",trigside);}
374 }
375 }
376
377 ((TH1D*)(fOutput->FindObject("zvSPD")))->Fill(fESD->GetPrimaryVertexSPD()->GetZ());
8320cad2 378 ((TH1D*)(fOutput->FindObject("nContrib")))->Fill(fESD->GetPrimaryVertexSPD()->GetNContributors());
7dd391ec 379 if (!fCutOnzVtxSPD || (fCutOnzVtxSPD && (fESD->GetPrimaryVertexSPD()->GetZ()>fzPrimVertexSPD[0] && fESD->GetPrimaryVertexSPD()->GetZ()<fzPrimVertexSPD[1]))){
380 ((TH1D*)(fOutput->FindObject("zvSPDcut")))->Fill(fESD->GetPrimaryVertexSPD()->GetZ());
381 if (!fCutOnNContributors || (fCutOnNContributors && (fESD->GetPrimaryVertexSPD()->GetNContributors()>fNContributors[0] && fESD->GetPrimaryVertexSPD()->GetNContributors()<fNContributors[1]))){
382 for (Int_t j = 0; j<mult1; j++) {
383
384 AliESDMuonTrack* mu1 = new AliESDMuonTrack(*(fESD->GetMuonTrack(j)));
385 if (!mu1->ContainTrackerData()) continue;
386 if (mu1->GetChi2()<fChi2Track[0] || mu1->GetChi2()>fChi2Track[1]) continue;
387 if (mu1->GetChi2MatchTrigger()<fChi2MatchTrig[0] || mu1->GetChi2MatchTrigger()>fChi2MatchTrig[1]) continue;
388 Double_t zr1 = mu1->Charge();
389 Double_t pxr1 = mu1->Px();
390 Double_t pyr1 = mu1->Py();
391 Double_t pzr1 = mu1->Pz();
392 Double_t ptmu1 = TMath::Sqrt(pxr1*pxr1+pyr1*pyr1);
393 Double_t er1 = mu1->E();
8320cad2 394 Double_t mom1 = mu1->P();
7dd391ec 395 Double_t theta1 = (180./TMath::Pi())*mu1->Theta();
396 Double_t rapiditymu1 = Rap(er1,pzr1);
8320cad2 397 ((TH1D*)(fOutput->FindObject("emuonREC")))->Fill(er1);
7dd391ec 398 ((TH1D*)(fOutput->FindObject("ptmuonREC")))->Fill(ptmu1);
399 ((TH1D*)(fOutput->FindObject("ymuonREC")))->Fill(rapiditymu1);
400
401 if(zr1<0){
402
403 for (Int_t jj = 0; jj<mult1; jj++) {
404
405 AliESDMuonTrack* mu2 = new AliESDMuonTrack(*(fESD->GetMuonTrack(jj)));
406 if (!mu2->ContainTrackerData()) continue;
407 if (mu2->GetChi2()<fChi2Track[0] || mu2->GetChi2()>fChi2Track[1]) continue;
408 if (mu2->GetChi2MatchTrigger()<fChi2MatchTrig[0] || mu2->GetChi2MatchTrigger()>fChi2MatchTrig[1]) continue;
409 Double_t zr2 = mu2->Charge();
410
411 if(zr2>0){
8320cad2 412 Double_t trigCondition=0;
413 if (mu1->GetMatchTrigger()>=mu2->GetMatchTrigger()) trigCondition = mu1->GetMatchTrigger()+0.1*mu2->GetMatchTrigger();
414 else trigCondition = mu2->GetMatchTrigger()+0.1*mu1->GetMatchTrigger();
415 ((TH1D*)(fOutput->FindObject("trigCond")))->Fill(trigCondition);
7dd391ec 416 Double_t pxr2 = mu2->Px();
417 Double_t pyr2 = mu2->Py();
418 Double_t pzr2 = mu2->Pz();
419 Double_t ptmu2 = TMath::Sqrt(pxr2*pxr2+pyr2*pyr2);
420 Double_t er2 = mu2->E();
8320cad2 421 Double_t mom2 = mu2->P();
7dd391ec 422 Double_t theta2 = (180./TMath::Pi())*mu2->Theta();
423
424 Double_t ptrec = TMath::Sqrt((pxr1+pxr2)*(pxr1+pxr2)+(pyr1+pyr2)*(pyr1+pyr2));
425 ((TH1D*)(fOutput->FindObject("ptdimuREC")))->Fill(ptrec);
426 Double_t raprec= Rap((er1+er2),(pzr1+pzr2));
427 ((TH1D*)(fOutput->FindObject("ydimuREC")))->Fill(raprec);
428 Double_t imassrec = Imass(er1,pxr1,pyr1,pzr1,er2,pxr2,pyr2,pzr2);
429 ((TH1D*)(fOutput->FindObject("imassTot")))->Fill(imassrec);
430
431 if (imassrec>12.) continue;
432
433 Double_t costCSrec = CostCS(pxr1,pyr1,pzr1,er1,zr1,pxr2,pyr2,pzr2,er2);
434 ((TH1D*)(fOutput->FindObject("costCSdimuREC")))->Fill(costCSrec);
8320cad2 435 if((Int_t)costCSrec==666) ((TH1D*)(fOutput->FindObject("costCSdimuFail")))->Fill(costCSrec);
7dd391ec 436 Double_t costHErec = CostHE(pxr1,pyr1,pzr1,er1,zr1,pxr2,pyr2,pzr2,er2);
437 ((TH1D*)(fOutput->FindObject("costHEdimuREC")))->Fill(costHErec);
8320cad2 438 if((Int_t)costHErec==666) ((TH1D*)(fOutput->FindObject("costHEdimuFail")))->Fill(costHErec);
7dd391ec 439 Double_t phiCSrec = PhiCS(pxr1,pyr1,pzr1,er1,zr1,pxr2,pyr2,pzr2,er2);
440 ((TH1D*)(fOutput->FindObject("phiCSdimuREC")))->Fill(phiCSrec);
8320cad2 441 if((Int_t)phiCSrec==666) ((TH1D*)(fOutput->FindObject("phiCSdimuFail")))->Fill(phiCSrec);
7dd391ec 442 Double_t phiHErec = PhiHE(pxr1,pyr1,pzr1,er1,zr1,pxr2,pyr2,pzr2,er2);
443 ((TH1D*)(fOutput->FindObject("phiHEdimuREC")))->Fill(phiHErec);
8320cad2 444 if((Int_t)phiHErec==666) ((TH1D*)(fOutput->FindObject("phiHEdimuFail")))->Fill(phiHErec);
7dd391ec 445
446 containerInput[0] = raprec ;
447 containerInput[1] = ptrec;
448 containerInput[2] = costHErec;
449 containerInput[3] = TMath::Abs(phiHErec);
450 containerInput[4] = costCSrec;
451 containerInput[5] = TMath::Abs(phiCSrec);
452 containerInput[6] = imassrec;
8320cad2 453 containerInput[7] = trigCondition+0.05;
7dd391ec 454 if (ptmu1<ptmu2) {
455 containerInput[8]=ptmu1;
456 containerInput[9]=ptmu2;
457 } else {
458 containerInput[8]=ptmu2;
459 containerInput[9]=ptmu1;
460 }
461 if (theta1<theta2) {
462 containerInput[10]=theta1;
463 containerInput[11]=theta2;
464 } else {
465 containerInput[10]=theta2;
466 containerInput[11]=theta1;
467 }
8320cad2 468 if (mom1<mom2) {
469 containerInput[12]=mom1;
470 containerInput[13]=mom2;
7dd391ec 471 } else {
8320cad2 472 containerInput[12]=mom2;
473 containerInput[13]=mom1;
7dd391ec 474 }
475 if (fDistinguishTrigClass && trigside==0) containerInput[14]=0.;
476 else if (fDistinguishTrigClass && trigside==1) containerInput[14]=1.;
477 else if (fDistinguishTrigClass && trigside==2) containerInput[14]=2.;
478 else if (fDistinguishTrigClass && trigside==3) containerInput[14]=3.;
479 else containerInput[14]=0.;
480
481 if(fDistinguishTrigClass){
482 if (trigfired==1) fCFManager->GetParticleContainer()->Fill(containerInput,5);
483 else if (trigfired==0) fCFManager->GetParticleContainer()->Fill(containerInput,1);
484 } else {
485 fCFManager->GetParticleContainer()->Fill(containerInput,1);
486 if(fIsAccProduction){
487 if(cutAccept==1) fCFManager->GetParticleContainer()->Fill(containerInput,3);
488 }
489 }
490 } // mu+ Selection
491
492 } // second mu Loop
493 } // mu- Selection
494 }
495 }
496 }
497 } else { // AOD-based ANALYSIS
498
499 AliAODEvent *aod;
500 AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
501 aod = aodH->GetEvent();
502 Int_t ntracks=aod->GetNumberOfTracks();
503
504 // MC part ---------------------------------------------------------------------------------------------------
505
506 if (fReadMCInfo){
507 TClonesArray *mcarray = dynamic_cast<TClonesArray*> (aod->FindListObject(AliAODMCParticle::StdBranchName())); //array of MC particles in this event
508 for(int ii=0;ii<mcarray->GetEntries();ii++){
509 AliAODMCParticle *mctrack = (AliAODMCParticle*) mcarray->At(ii);
510 if(mctrack->GetPdgCode()!=13) continue;
8320cad2 511 Int_t numbMother = mctrack->GetMother();
512 if (numbMother==-1) continue;
513 AliAODMCParticle *mother = (AliAODMCParticle*) mcarray->At(numbMother);
7dd391ec 514 if (mother->GetPdgCode()!=443) continue;
8320cad2 515 numbJpsis++;
7dd391ec 516 Int_t daught0 = mother->GetDaughter(0);
517 Int_t daught1 = mother->GetDaughter(1);
518 AliAODMCParticle *mcDaughter0 = (AliAODMCParticle*) mcarray->At(daught0);
519 Double_t pxmc0 = mcDaughter0->Px();
520 Double_t pymc0 = mcDaughter0->Py();
521 Double_t pzmc0 = mcDaughter0->Pz();
522 Double_t ptmc0 = mcDaughter0->Pt();
8320cad2 523 Double_t mommc0 = mcDaughter0->P();
524 Double_t emc0 = mcDaughter0->E();
7dd391ec 525 Double_t thetamc0 = (180./TMath::Pi())*mcDaughter0->Theta();
526 Int_t charge0 = (Int_t) mcDaughter0->Charge()/3;
527 AliAODMCParticle *mcDaughter1 = (AliAODMCParticle*) mcarray->At(daught1);
528 Double_t pxmc1 = mcDaughter1->Px();
529 Double_t pymc1 = mcDaughter1->Py();
530 Double_t pzmc1 = mcDaughter1->Pz();
531 Double_t ptmc1 = mcDaughter1->Pt();
8320cad2 532 Double_t mommc1 = mcDaughter1->P();
533 Double_t emc1 = mcDaughter1->E();
7dd391ec 534 Double_t thetamc1 = (180./TMath::Pi())*mcDaughter1->Theta();
535 Int_t charge1 = (Int_t) mcDaughter1->Charge()/3;
536
537 if (charge0==charge1 || TMath::Abs(mcDaughter0->GetPdgCode())!=13 || TMath::Abs(mcDaughter1->GetPdgCode())!=13) continue;
7dd391ec 538 Double_t rapJpsi = Rap(mother->E(),mother->Pz());
539 Double_t ptJpsi = mother->Pt();
8320cad2 540 Double_t costHE = CostHE(pxmc0,pymc0,pzmc0,emc0,(Double_t)charge0,pxmc1,pymc1,pzmc1,emc1);
541 Double_t phiHE = PhiHE(pxmc0,pymc0,pzmc0,emc0,(Double_t)charge0,pxmc1,pymc1,pzmc1,emc1);
542 Double_t costCS = CostCS(pxmc0,pymc0,pzmc0,emc0,(Double_t)charge0,pxmc1,pymc1,pzmc1,emc1);
543 Double_t phiCS = PhiCS(pxmc0,pymc0,pzmc0,emc0,(Double_t)charge0,pxmc1,pymc1,pzmc1,emc1);
7dd391ec 544 Double_t massJpsi = mother->M();
545 containerInput[0] = rapJpsi;
546 containerInput[1] = ptJpsi;
547 containerInput[2] = costHE;
548 containerInput[3] = TMath::Abs(phiHE);
549 containerInput[4] = costCS;
550 containerInput[5] = TMath::Abs(phiCS);
551 containerInput[6] = massJpsi;
552 containerInput[7] = 1.; //for generated no trigger condition
553 if (ptmc0<ptmc1) {
554 containerInput[8]=ptmc0;
555 containerInput[9]=ptmc1;
556 } else {
557 containerInput[8]=ptmc1;
558 containerInput[9]=ptmc0;
559 }
560 if (thetamc0<thetamc1) {
561 containerInput[10]=thetamc0;
562 containerInput[11]=thetamc1;
563 } else {
564 containerInput[10]=thetamc1;
565 containerInput[11]=thetamc0;
566 }
8320cad2 567 if (mommc0<mommc1) {
568 containerInput[12]=mommc0;
569 containerInput[13]=mommc1;
7dd391ec 570 } else {
8320cad2 571 containerInput[12]=mommc1;
572 containerInput[13]=mommc0;
7dd391ec 573 }
574 containerInput[14]=1.;
575
8320cad2 576 if((Int_t)rapJpsi!=666 && (Int_t)costHE!=666 && (Int_t)phiHE!=666 && (Int_t)costCS!=666 && (Int_t)phiCS!=666){
7dd391ec 577 fCFManager->GetParticleContainer()->Fill(containerInput,0);
578 if(fIsAccProduction){
579 // acceptance cuts on single mu
580 if(thetamc0<fThetaSingMuCut[0]||thetamc0>fThetaSingMuCut[1]||thetamc1<fThetaSingMuCut[0]||thetamc1>fThetaSingMuCut[1]) cutAccept=0;
581 if(ptmc0<fPtSingMuCut[0] || ptmc0>fPtSingMuCut[1] || ptmc1<fPtSingMuCut[0] || ptmc1>fPtSingMuCut[1]) cutAccept=0;
582 if (cutAccept ==1) fCFManager->GetParticleContainer()->Fill(containerInput,2);
583 }
584
585 }
586
587 }
588 }
589
590
591 // AOD part ---------------------------------------------------------------------------------------------------
592
593 Int_t trigfired=-1;
594 Int_t trigside=-1;
595 if(fDistinguishTrigClass){
596 TString trigclass = aod->GetFiredTriggerClasses();
597 printf("TrigClass: %s\n",trigclass.Data());
598 if(trigclass.Contains(fTrigClassMuon)) trigfired = 1;
599 else if (trigclass.Contains(fTrigClassInteraction)) trigfired = 0;
600 printf("Muon container: %d\n",trigfired);
601 for(Int_t i=0;i<4;i++){
602 if(trigfired==1 && trigclass.Contains(fTrigClassMuonSide[i])) {trigside=i;printf("Muon Side: %d\n",trigside);}
603 if(trigfired==0 && trigclass.Contains(fTrigClassInteractionSide[i])) {trigside=i;printf("Interaction Side: %d\n",trigside);}
604 }
605 }
606
607 ((TH1D*)(fOutput->FindObject("zvSPD")))->Fill(aod->GetPrimaryVertex()->GetZ());
608 Int_t ncontr = aod->GetPrimaryVertex()->GetNContributors();
8320cad2 609 ((TH1D*)(fOutput->FindObject("nContrib")))->Fill(ncontr);
7dd391ec 610 if (!fCutOnzVtxSPD || (fCutOnzVtxSPD && (aod->GetPrimaryVertex()->GetZ()>fzPrimVertexSPD[0] && aod->GetPrimaryVertex()->GetZ()<fzPrimVertexSPD[1]))){ //NOT REALLY SPD VERTEX
611 ((TH1D*)(fOutput->FindObject("zvSPDcut")))->Fill(aod->GetPrimaryVertex()->GetZ());
612 if (!fCutOnNContributors || (fCutOnNContributors && (aod->GetPrimaryVertex()->GetNContributors()>fNContributors[0] && aod->GetPrimaryVertex()->GetNContributors()<fNContributors[1]))){
613 for (Int_t j = 0; j<ntracks; j++) {
614 AliAODTrack *mu1 = aod->GetTrack(j);
615 if(!mu1->IsMuonTrack()) continue;
616 if (mu1->Chi2perNDF()<fChi2Track[0] || mu1->Chi2perNDF()>fChi2Track[1]) continue;
617 if (mu1->GetChi2MatchTrigger()<fChi2MatchTrig[0] || mu1->GetChi2MatchTrigger()>fChi2MatchTrig[1]) continue;
618 Double_t chargemu1 = mu1->Charge();
619 Double_t pxmu1 = mu1->Px();
620 Double_t pymu1 = mu1->Py();
621 Double_t pzmu1 = mu1->Pz();
622 Double_t emu1 = mu1->E();
623 Double_t pmu1 = mu1->P();
624 Double_t ptmu1 = mu1->Pt();
625 Double_t rapiditymu1 = Rap(emu1,pzmu1);
626 Double_t thetamu1 = (180./TMath::Pi())*mu1->Theta();
8320cad2 627 Double_t rdcamu1 = mu1->DCA();
628 ((TH1D*)(fOutput->FindObject("emuonREC")))->Fill(emu1);
7dd391ec 629 ((TH1D*)(fOutput->FindObject("ptmuonREC")))->Fill(ptmu1);
630 ((TH1D*)(fOutput->FindObject("ymuonREC")))->Fill(rapiditymu1);
631 if(chargemu1<0){
632 for (Int_t jj = 0; jj<ntracks; jj++) {
633 AliAODTrack *mu2 = aod->GetTrack(jj);
634 if(!mu2->IsMuonTrack()) continue;
635 if (mu2->Chi2perNDF()<fChi2Track[0] || mu2->Chi2perNDF()>fChi2Track[1]) continue;
636 if (mu2->GetChi2MatchTrigger()<fChi2MatchTrig[0] || mu2->GetChi2MatchTrigger()>fChi2MatchTrig[1]) continue;
637 Double_t chargemu2 = mu2->Charge();
638 Double_t pxmu2 = mu2->Px();
639 Double_t pymu2 = mu2->Py();
640 Double_t pzmu2 = mu2->Pz();
641 Double_t emu2 = mu2->E();
642 Double_t pmu2 = mu2->P();
643 Double_t ptmu2 = mu2->Pt();
644 //Double_t rapiditymu2 = Rap(emu2,pzmu2);
645 Double_t thetamu2 = (180./TMath::Pi())*mu2->Theta();
8320cad2 646 Double_t rdcamu2 = mu2->DCA();
7dd391ec 647 if(chargemu2>0){
648 if (mu1->GetMatchTrigger()>0) printf("Mu1: charge: %f, match: %d\n",chargemu1,1);
649 else printf("Mu1: charge: %f, match: %d\n",chargemu1,0);
650 if (mu2->GetMatchTrigger()>0) printf("Mu2: charge: %f, match: %d\n",chargemu2,1);
651 else printf("Mu2: charge: %f, match: %d\n",chargemu2,0);
8320cad2 652 ((TH1D*)(fOutput->FindObject("hdca")))->Fill(rdcamu2);
653 ((TH1D*)(fOutput->FindObject("hdca")))->Fill(rdcamu1);
654 Double_t trigCondition=0;
655 if (mu1->GetMatchTrigger()>=mu2->GetMatchTrigger()) trigCondition = mu1->GetMatchTrigger()+0.1*mu2->GetMatchTrigger();
656 else trigCondition = mu2->GetMatchTrigger()+0.1*mu1->GetMatchTrigger();
7dd391ec 657 containerInput[0] = (Double_t) Rap((emu1+emu2),(pzmu1+pzmu2));
658 ((TH1D*)(fOutput->FindObject("ydimuREC")))->Fill(containerInput[0]);
8320cad2 659 ((TH2D*)(fOutput->FindObject("hdcay")))->Fill(TMath::Max(rdcamu1,rdcamu2),containerInput[0]);
7dd391ec 660 //((TH2D*)(fOutput->FindObject("hdcay")))->Fill(Rdcamu2,containerInput[0]);
661 containerInput[1] = (Double_t) TMath::Sqrt((pxmu1+pxmu2)*(pxmu1+pxmu2)+(pymu1+pymu2)*(pymu1+pymu2));
662 ((TH1D*)(fOutput->FindObject("ptdimuREC")))->Fill(containerInput[1]);
663 containerInput[2] = CostHE(pxmu1,pymu1,pzmu1,emu1,chargemu1,pxmu2,pymu2,pzmu2,emu2);
664 ((TH1D*)(fOutput->FindObject("costHEdimuREC")))->Fill(containerInput[2]);
665 if(containerInput[2]==666) ((TH1D*)(fOutput->FindObject("costHEdimuFail")))->Fill(containerInput[2]);
666 containerInput[3] = TMath::Abs(PhiHE(pxmu1,pymu1,pzmu1,emu1,chargemu1,pxmu2,pymu2,pzmu2,emu2));
667 ((TH1D*)(fOutput->FindObject("phiHEdimuREC")))->Fill(containerInput[3]);
668 if(containerInput[3]==666) ((TH1D*)(fOutput->FindObject("phiHEdimuFail")))->Fill(containerInput[3]);
669 containerInput[4] = CostCS(pxmu1,pymu1,pzmu1,emu1,chargemu1,pxmu2,pymu2,pzmu2,emu2);
670 ((TH1D*)(fOutput->FindObject("costCSdimuREC")))->Fill(containerInput[4]);
671 if(containerInput[4]==666) ((TH1D*)(fOutput->FindObject("costCSdimuFail")))->Fill(containerInput[4]);
672 containerInput[5] = TMath::Abs(PhiCS(pxmu1,pymu1,pzmu1,emu1,chargemu1,pxmu2,pymu2,pzmu2,emu2));
673 ((TH1D*)(fOutput->FindObject("phiCSdimuREC")))->Fill(containerInput[5]);
674 if(containerInput[5]==666) ((TH1D*)(fOutput->FindObject("phiCSdimuFail")))->Fill(containerInput[5]);
675 containerInput[6] = Imass(emu1,pxmu1,pymu1,pzmu1,emu2,pxmu2,pymu2,pzmu2);
676 ((TH1D*)(fOutput->FindObject("imassTot")))->Fill(containerInput[6]);
8320cad2 677 containerInput[7] = trigCondition+0.05;
7dd391ec 678 ((TH1D*)(fOutput->FindObject("trigCond")))->Fill(containerInput[7]);
679 if (ptmu1<ptmu2) {
680 containerInput[8]=ptmu1;
681 containerInput[9]=ptmu2;
682 } else {
683 containerInput[8]=ptmu2;
684 containerInput[9]=ptmu1;
685 }
686 if (thetamu1<thetamu2) {
687 containerInput[10]=thetamu1;
688 containerInput[11]=thetamu2;
689 } else {
690 containerInput[10]=thetamu2;
691 containerInput[11]=thetamu1;
692 }
693 if (pmu1<pmu2) {
694 containerInput[12]=pmu1;
695 containerInput[13]=pmu2;
696 } else {
697 containerInput[12]=pmu2;
698 containerInput[13]=pmu1;
699 }
700 if (fDistinguishTrigClass && trigside==0) containerInput[14]=0.;
701 else if (fDistinguishTrigClass && trigside==1) containerInput[14]=1.;
702 else if (fDistinguishTrigClass && trigside==2) containerInput[14]=2.;
703 else if (fDistinguishTrigClass && trigside==3) containerInput[14]=3.;
704 else containerInput[14]=0.;
705
706 if(containerInput[2]!=666 && containerInput[3]!=666 && containerInput[4]!=666 && containerInput[5]!=666){
707 if(fDistinguishTrigClass){
708 if (trigfired==1) fCFManager->GetParticleContainer()->Fill(containerInput,5);
709 else if (trigfired==0) fCFManager->GetParticleContainer()->Fill(containerInput,1);
710 } else {
711 fCFManager->GetParticleContainer()->Fill(containerInput,1);
712 if(fIsAccProduction){
713 if (cutAccept ==1) fCFManager->GetParticleContainer()->Fill(containerInput,3);
714 }
715 }
716 } else if (containerInput[0]==666) ((TH1D*)(fOutput->FindObject("ydimuFail")))->Fill(containerInput[0]);
717
718 }
719
720 }
721
722 }
723
724 }
725 }
726 }
727
728
729 }
730
731
732// ----------
8320cad2 733 if (fReadMCInfo) ((TH1D*)(fOutput->FindObject("jpsiMult")))->Fill(numbJpsis);
7dd391ec 734
735 PostData(1,fOutput) ;
736 PostData(2,fCFManager->GetParticleContainer()) ;
737
738
739}
740
741//________________________________________________________________________
742Double_t AliAnalysisTaskDimuonCFContainerBuilder::Imass(Double_t e1, Double_t px1, Double_t py1, Double_t pz1,
8320cad2 743 Double_t e2, Double_t px2, Double_t py2, Double_t pz2) const
7dd391ec 744{
745// invariant mass calculation
746 Double_t imassrec = TMath::Sqrt((e1+e2)*(e1+e2)-((px1+px2)*(px1+px2)+
747 (py1+py2)*(py1+py2)+(pz1+pz2)*(pz1+pz2)));
748 return imassrec;
749}
750
751//________________________________________________________________________
8320cad2 752Double_t AliAnalysisTaskDimuonCFContainerBuilder::Rap(Double_t e, Double_t pz) const
7dd391ec 753{
754// calculate rapidity
755 Double_t rap;
756 if(e>TMath::Abs(pz)){ // in order to avoid problems with AODs
757 rap = 0.5*TMath::Log((e+pz)/(e-pz));
758 return rap;
759 }
760 else{
761 rap = 666.;
762 return rap;
763 }
764}
765
766//________________________________________________________________________
767Double_t AliAnalysisTaskDimuonCFContainerBuilder::CostCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
768Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
769{
770// Cosine of the theta decay angle (mu+) in the Collins-Soper frame
771
8320cad2 772 TLorentzVector pMu1CM, pMu2CM, pProjCM, pTargCM, pDimuCM; // In the CM. frame
773 TLorentzVector pMu1Dimu, pMu2Dimu, pProjDimu, pTargDimu; // In the dimuon rest frame
7dd391ec 774 TVector3 beta,zaxisCS;
775 Double_t mp=0.93827231;
776 //
777 // --- Fill the Lorentz vector for projectile and target in the CM frame
778 //
8320cad2 779 pProjCM.SetPxPyPzE(0.,0.,-fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp));
780 pTargCM.SetPxPyPzE(0.,0.,fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp));
7dd391ec 781 //
782 // --- Get the muons parameters in the CM frame
783 //
8320cad2 784 pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
785 pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
7dd391ec 786 //
787 // --- Obtain the dimuon parameters in the CM frame
788 //
8320cad2 789 pDimuCM=pMu1CM+pMu2CM;
7dd391ec 790 //
791 // --- Translate the dimuon parameters in the dimuon rest frame
792 //
8320cad2 793 beta=(-1./pDimuCM.E())*pDimuCM.Vect();
7dd391ec 794 if(beta.Mag()>=1) return 666.;
8320cad2 795 pMu1Dimu=pMu1CM;
796 pMu2Dimu=pMu2CM;
797 pProjDimu=pProjCM;
798 pTargDimu=pTargCM;
799 pMu1Dimu.Boost(beta);
800 pMu2Dimu.Boost(beta);
801 pProjDimu.Boost(beta);
802 pTargDimu.Boost(beta);
7dd391ec 803
804 //Debugging part -------------------------------------
805 Double_t debugProj[4]={0.,0.,0.,0.};
806 Double_t debugTarg[4]={0.,0.,0.,0.};
807 Double_t debugMu1[4]={0.,0.,0.,0.};
808 Double_t debugMu2[4]={0.,0.,0.,0.};
8320cad2 809 pMu1Dimu.GetXYZT(debugMu1);
810 pMu2Dimu.GetXYZT(debugMu2);
811 pProjDimu.GetXYZT(debugProj);
812 pTargDimu.GetXYZT(debugTarg);
7dd391ec 813 if (debugProj[0]!=debugProj[0] ||debugProj[1]!=debugProj[1] || debugProj[2]!=debugProj[2] ||debugProj[3]!=debugProj[3]) return 666;
814 if (debugTarg[0]!=debugTarg[0] ||debugTarg[1]!=debugTarg[1] || debugTarg[2]!=debugTarg[2] ||debugTarg[3]!=debugTarg[3]) return 666;
815 if (debugMu1[0]!=debugMu1[0] ||debugMu1[1]!=debugMu1[1] || debugMu1[2]!=debugMu1[2] ||debugMu1[3]!=debugMu1[3]) return 666;
816 if (debugMu2[0]!=debugMu2[0] ||debugMu2[1]!=debugMu2[1] || debugMu2[2]!=debugMu2[2] ||debugMu2[3]!=debugMu2[3]) return 666;
817 //----------------------------------------------------
818
819 // --- Determine the z axis for the CS angle
8320cad2 820 zaxisCS=(((pProjDimu.Vect()).Unit())-((pTargDimu.Vect()).Unit())).Unit();
7dd391ec 821
822 // --- Determine the CS angle (angle between mu+ and the z axis defined above)
823 Double_t cost;
824
8320cad2 825 if(charge1>0) {cost = zaxisCS.Dot((pMu1Dimu.Vect()).Unit());}
826 else {cost = zaxisCS.Dot((pMu2Dimu.Vect()).Unit());}
7dd391ec 827
828 return cost;
829}
830
831//________________________________________________________________________
832Double_t AliAnalysisTaskDimuonCFContainerBuilder::CostHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
833Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
834{
835// Cosine of the theta decay angle (mu+) in the Helicity frame
836
8320cad2 837 TLorentzVector pMu1CM, pMu2CM, pDimuCM; // In the CM frame
838 TLorentzVector pMu1Dimu, pMu2Dimu; // In the dimuon rest frame
7dd391ec 839 TVector3 beta,zaxisCS;
840 //
841 // --- Get the muons parameters in the CM frame
842 //
8320cad2 843 pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
844 pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
7dd391ec 845 //
846 // --- Obtain the dimuon parameters in the CM frame
847 //
8320cad2 848 pDimuCM=pMu1CM+pMu2CM;
7dd391ec 849 //
850 // --- Translate the muon parameters in the dimuon rest frame
851 //
8320cad2 852 beta=(-1./pDimuCM.E())*pDimuCM.Vect();
7dd391ec 853 if(beta.Mag()>=1) return 666.;
8320cad2 854 pMu1Dimu=pMu1CM;
855 pMu2Dimu=pMu2CM;
856 pMu1Dimu.Boost(beta);
857 pMu2Dimu.Boost(beta);
7dd391ec 858
859 //Debugging part -------------------------------------
860 Double_t debugMu1[4]={0.,0.,0.,0.};
861 Double_t debugMu2[4]={0.,0.,0.,0.};
8320cad2 862 pMu1Dimu.GetXYZT(debugMu1);
863 pMu2Dimu.GetXYZT(debugMu2);
7dd391ec 864 if (debugMu1[0]!=debugMu1[0] ||debugMu1[1]!=debugMu1[1] || debugMu1[2]!=debugMu1[2] ||debugMu1[3]!=debugMu1[3]) return 666;
865 if (debugMu2[0]!=debugMu2[0] ||debugMu2[1]!=debugMu2[1] || debugMu2[2]!=debugMu2[2] ||debugMu2[3]!=debugMu2[3]) return 666;
866 //----------------------------------------------------
867
868 // --- Determine the z axis for the calculation of the polarization angle (i.e. the direction of the dimuon in the CM system)
869 TVector3 zaxis;
8320cad2 870 zaxis=(pDimuCM.Vect()).Unit();
7dd391ec 871
872 // --- Calculation of the polarization angle (angle between mu+ and the z axis defined above)
873 Double_t cost;
8320cad2 874 if(charge1>0) {cost = zaxis.Dot((pMu1Dimu.Vect()).Unit());}
875 else {cost = zaxis.Dot((pMu2Dimu.Vect()).Unit());}
7dd391ec 876 return cost;
877}
878
879//________________________________________________________________________
880Double_t AliAnalysisTaskDimuonCFContainerBuilder::PhiCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
881Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
882{
883// Phi decay angle (mu+) in the Collins-Soper frame
884
8320cad2 885 TLorentzVector pMu1CM, pMu2CM, pProjCM, pTargCM, pDimuCM; // In the CM frame
886 TLorentzVector pMu1Dimu, pMu2Dimu, pProjDimu, pTargDimu; // In the dimuon rest frame
7dd391ec 887 TVector3 beta,yaxisCS, xaxisCS, zaxisCS;
888 Double_t mp=0.93827231;
889
890 // --- Fill the Lorentz vector for projectile and target in the CM frame
8320cad2 891 pProjCM.SetPxPyPzE(0.,0.,-fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp));
892 pTargCM.SetPxPyPzE(0.,0.,fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp));
7dd391ec 893
894 // --- Get the muons parameters in the CM frame
8320cad2 895 pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
896 pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
7dd391ec 897
898 // --- Obtain the dimuon parameters in the CM frame
8320cad2 899 pDimuCM=pMu1CM+pMu2CM;
7dd391ec 900
901 // --- Translate the dimuon parameters in the dimuon rest frame
8320cad2 902 beta=(-1./pDimuCM.E())*pDimuCM.Vect();
7dd391ec 903 if(beta.Mag()>=1) return 666.;
8320cad2 904 pMu1Dimu=pMu1CM;
905 pMu2Dimu=pMu2CM;
906 pProjDimu=pProjCM;
907 pTargDimu=pTargCM;
908 pMu1Dimu.Boost(beta);
909 pMu2Dimu.Boost(beta);
910 pProjDimu.Boost(beta);
911 pTargDimu.Boost(beta);
7dd391ec 912
913 //Debugging part -------------------------------------
914 Double_t debugProj[4]={0.,0.,0.,0.};
915 Double_t debugTarg[4]={0.,0.,0.,0.};
916 Double_t debugMu1[4]={0.,0.,0.,0.};
917 Double_t debugMu2[4]={0.,0.,0.,0.};
8320cad2 918 pMu1Dimu.GetXYZT(debugMu1);
919 pMu2Dimu.GetXYZT(debugMu2);
920 pProjDimu.GetXYZT(debugProj);
921 pTargDimu.GetXYZT(debugTarg);
7dd391ec 922 if (debugProj[0]!=debugProj[0] ||debugProj[1]!=debugProj[1] || debugProj[2]!=debugProj[2] ||debugProj[3]!=debugProj[3]) return 666;
923 if (debugTarg[0]!=debugTarg[0] ||debugTarg[1]!=debugTarg[1] || debugTarg[2]!=debugTarg[2] ||debugTarg[3]!=debugTarg[3]) return 666;
924 if (debugMu1[0]!=debugMu1[0] ||debugMu1[1]!=debugMu1[1] || debugMu1[2]!=debugMu1[2] ||debugMu1[3]!=debugMu1[3]) return 666;
925 if (debugMu2[0]!=debugMu2[0] ||debugMu2[1]!=debugMu2[1] || debugMu2[2]!=debugMu2[2] ||debugMu2[3]!=debugMu2[3]) return 666;
926 //----------------------------------------------------
927
928 // --- Determine the z axis for the CS angle
8320cad2 929 zaxisCS=(((pProjDimu.Vect()).Unit())-((pTargDimu.Vect()).Unit())).Unit();
930 yaxisCS=(((pProjDimu.Vect()).Unit()).Cross((pTargDimu.Vect()).Unit())).Unit();
7dd391ec 931 xaxisCS=(yaxisCS.Cross(zaxisCS)).Unit();
932
933 Double_t phi=0.;
934 if(charge1>0) {
8320cad2 935 phi = TMath::ATan2((pMu1Dimu.Vect()).Dot(yaxisCS),((pMu1Dimu.Vect()).Dot(xaxisCS)));
7dd391ec 936 } else {
8320cad2 937 phi = TMath::ATan2((pMu2Dimu.Vect()).Dot(yaxisCS),((pMu2Dimu.Vect()).Dot(xaxisCS)));
7dd391ec 938 }
939 if (phi>TMath::Pi()) phi=phi-TMath::Pi();
940
941 return phi;
942}
943
944//________________________________________________________________________
945Double_t AliAnalysisTaskDimuonCFContainerBuilder::PhiHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
946Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
947{
948// Phi decay angle (mu+) in the Helicity frame
8320cad2 949 TLorentzVector pMu1Lab, pMu2Lab, pProjLab, pTargLab, pDimuLab; // In the lab. frame
950 TLorentzVector pMu1Dimu, pMu2Dimu, pProjDimu, pTargDimu; // In the dimuon rest frame
7dd391ec 951 TVector3 beta,xaxis, yaxis,zaxis;
952 Double_t mp=0.93827231;
953
954 // --- Get the muons parameters in the LAB frame
8320cad2 955 pMu1Lab.SetPxPyPzE(px1,py1,pz1,e1);
956 pMu2Lab.SetPxPyPzE(px2,py2,pz2,e2);
7dd391ec 957
958 // --- Obtain the dimuon parameters in the LAB frame
8320cad2 959 pDimuLab=pMu1Lab+pMu2Lab;
960 zaxis=(pDimuLab.Vect()).Unit();
7dd391ec 961
962 // --- Translate the muon parameters in the dimuon rest frame
8320cad2 963 beta=(-1./pDimuLab.E())*pDimuLab.Vect();
7dd391ec 964 if(beta.Mag()>=1.) return 666.;
965
8320cad2 966 pProjLab.SetPxPyPzE(0.,0.,-fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp)); // proiettile
967 pTargLab.SetPxPyPzE(0.,0.,fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp)); // bersaglio
7dd391ec 968
8320cad2 969 pProjDimu=pProjLab;
970 pTargDimu=pTargLab;
7dd391ec 971
8320cad2 972 pProjDimu.Boost(beta);
973 pTargDimu.Boost(beta);
7dd391ec 974
8320cad2 975 yaxis=((pProjDimu.Vect()).Cross(pTargDimu.Vect())).Unit();
7dd391ec 976 xaxis=(yaxis.Cross(zaxis)).Unit();
977
8320cad2 978 pMu1Dimu=pMu1Lab;
979 pMu2Dimu=pMu2Lab;
980 pMu1Dimu.Boost(beta);
981 pMu2Dimu.Boost(beta);
7dd391ec 982
983 //Debugging part -------------------------------------
984 Double_t debugProj[4]={0.,0.,0.,0.};
985 Double_t debugTarg[4]={0.,0.,0.,0.};
986 Double_t debugMu1[4]={0.,0.,0.,0.};
987 Double_t debugMu2[4]={0.,0.,0.,0.};
8320cad2 988 pMu1Dimu.GetXYZT(debugMu1);
989 pMu2Dimu.GetXYZT(debugMu2);
990 pProjDimu.GetXYZT(debugProj);
991 pTargDimu.GetXYZT(debugTarg);
7dd391ec 992 if (debugProj[0]!=debugProj[0] ||debugProj[1]!=debugProj[1] || debugProj[2]!=debugProj[2] ||debugProj[3]!=debugProj[3]) return 666;
993 if (debugTarg[0]!=debugTarg[0] ||debugTarg[1]!=debugTarg[1] || debugTarg[2]!=debugTarg[2] ||debugTarg[3]!=debugTarg[3]) return 666;
994 if (debugMu1[0]!=debugMu1[0] ||debugMu1[1]!=debugMu1[1] || debugMu1[2]!=debugMu1[2] ||debugMu1[3]!=debugMu1[3]) return 666;
995 if (debugMu2[0]!=debugMu2[0] ||debugMu2[1]!=debugMu2[1] || debugMu2[2]!=debugMu2[2] ||debugMu2[3]!=debugMu2[3]) return 666;
996 //----------------------------------------------------
997
998 Double_t phi=0.;
999 if(charge1 > 0) {
8320cad2 1000 phi = TMath::ATan2((pMu1Dimu.Vect()).Dot(yaxis),(pMu1Dimu.Vect()).Dot(xaxis));
7dd391ec 1001 } else {
8320cad2 1002 phi = TMath::ATan2((pMu2Dimu.Vect()).Dot(yaxis),(pMu2Dimu.Vect()).Dot(xaxis));
7dd391ec 1003 }
1004 return phi;
1005}
1006
1007//________________________________________________________________________
1008void AliAnalysisTaskDimuonCFContainerBuilder::Terminate(Option_t *)
1009{
1010
1011}
1012
1013#endif