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