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