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