Fixing coding rule violations (Sang-Un)
[u/mrichter/AliRoot.git] / PWG3 / muon / AliCFMuonResUpsilon.cxx
CommitLineData
014044ba 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16//-----------------------------------------------------------------------
17// Example of task (running locally, on AliEn and CAF),
18// which provides standard way of calculating acceptance and efficiency
19// between different steps of the procedure.
20// The ouptut of the task is a AliCFContainer from which the efficiencies
21// can be calculated
22//-----------------------------------------------------------------------
23// Author : R. Vernet, Consorzio Cometa - Catania (it)
24//-----------------------------------------------------------------------
25// Modification done by X. Lopez - LPC Clermont (fr)
26//-----------------------------------------------------------------------
27// Modification done by S. Ahn - LPC Clermont (fr), Konkuk University (kr)
28//-----------------------------------------------------------------------
29
30
31
32#ifndef ALICFMUONRESUPSILON_CXX
33#define ALICFMUONRESUPSILON_CXX
34
35#include "TH1.h"
36#include "TParticle.h"
37#include "TChain.h"
38#include "TLorentzVector.h"
39#include "TCanvas.h"
40
41#include "AliHeader.h"
42#include "AliStack.h"
43#include "AliMCEvent.h"
44#include "AliAnalysisManager.h"
45#include "AliLog.h"
46
47#include "AliVEvent.h"
48#include "AliPhysicsSelection.h"
49
50#include "AliCFMuonResUpsilon.h"
51#include "AliCFManager.h"
52#include "AliCFCutBase.h"
53#include "AliCFContainer.h"
54
55#include "AliESDEvent.h"
56#include "AliESDHeader.h"
57#include "AliESDtrack.h"
58#include "AliESDMuonTrack.h"
59#include "AliESDtrack.h"
60#include "AliESDInputHandler.h"
61
62#include "AliAODEvent.h"
63#include "AliAODInputHandler.h"
64#include "AliAODMCParticle.h"
65
66ClassImp(AliCFMuonResUpsilon)
67
68//__________________________________________________________________________
69AliCFMuonResUpsilon::AliCFMuonResUpsilon() :
70 AliAnalysisTaskSE(""),
e199f3f0 71 fReadAODData(kFALSE),
72 fReadMCInfo(kFALSE),
014044ba 73 fCFManager(0x0),
6579f87c 74 fnevts(0x0),
014044ba 75 fIsPhysSelMB(kFALSE),
76 fIsPhysSelMUON(kFALSE),
77 fQAHistList(0x0),
78 fPDG(0),
79 fPtMin(0),
80 fPtMax(0),
81 fYMin(0),
82 fYMax(0),
83 fTrigClassMuon(""),
84 fTrigClassInteraction(""),
85 fDistinguishTrigClass(kTRUE)
86{
87
88 //Default ctor
89
90}
91//___________________________________________________________________________
92AliCFMuonResUpsilon::AliCFMuonResUpsilon(const Char_t* name) :
93 AliAnalysisTaskSE(name),
e199f3f0 94 fReadAODData(kFALSE),
95 fReadMCInfo(kFALSE),
014044ba 96 fCFManager(0x0),
6579f87c 97 fnevts(0x0),
014044ba 98 fIsPhysSelMB(kFALSE),
99 fIsPhysSelMUON(kFALSE),
100 fQAHistList(0x0),
101 fPDG(0),
102 fPtMin(0),
103 fPtMax(0),
104 fYMin(0),
105 fYMax(0),
106 fTrigClassMuon(""),
107 fTrigClassInteraction(""),
108 fDistinguishTrigClass(kTRUE)
109{
110
111 // Constructor. Initialization of Inputs and Outputs
112
113 Info("AliCFMuonResUpsilon","Calling Constructor");
114
6579f87c 115 fnevts = new TH1D("fnevts","fnevts",5,0,5); // nevent CINT1B,CMUS1B + PhysSel
014044ba 116
117 TString nameside[3]={"AC","B","E"};
118
119 SetTrigClassMuonName();
120 SetTrigClassInteracName();
121 SetTrigClassSideName(nameside);
122
123 DefineInput(0, TChain::Class());
124 DefineOutput(1,TH1D::Class());
125 DefineOutput(2,AliCFContainer::Class());
126}
127
128//___________________________________________________________________________
129AliCFMuonResUpsilon& AliCFMuonResUpsilon::operator=(const AliCFMuonResUpsilon& c)
130{
131
132 // Assignment operator
133
134 if (this!=&c) {
135 AliAnalysisTaskSE::operator=(c);
136 fReadAODData = c.fReadAODData;
e199f3f0 137 fReadMCInfo = c.fReadMCInfo;
014044ba 138 fCFManager = c.fCFManager;
139 fQAHistList = c.fQAHistList;
6579f87c 140 fnevts = c.fnevts;
014044ba 141 fPDG = c.fPDG;
142 fPtMin = c.fPtMin;
143 fPtMax = c.fPtMax;
144 fYMin = c.fYMin;
145 fYMax = c.fYMax;
146 }
147 return *this;
148}
149
150//___________________________________________________________________________
151AliCFMuonResUpsilon::AliCFMuonResUpsilon(const AliCFMuonResUpsilon& c) :
152 AliAnalysisTaskSE(c),
153 fReadAODData(c.fReadAODData),
154 fReadMCInfo(c.fReadMCInfo),
155 fCFManager(c.fCFManager),
6579f87c 156 fnevts(c.fnevts),
014044ba 157 fIsPhysSelMB(kFALSE),
158 fIsPhysSelMUON(kFALSE),
159 fQAHistList(c.fQAHistList),
160 fPDG(c.fPDG),
161 fPtMin(c.fPtMin),
162 fPtMax(c.fPtMax),
163 fYMin(c.fYMin),
164 fYMax(c.fYMax),
165 fTrigClassMuon(c.fTrigClassMuon),
166 fTrigClassInteraction(c.fTrigClassInteraction),
167 fDistinguishTrigClass(c.fDistinguishTrigClass)
168{
169
170 // Copy Constructor
171
172}
173
174//___________________________________________________________________________
175AliCFMuonResUpsilon::~AliCFMuonResUpsilon() {
176
177 //destructor
178
179 Info("~AliCFMuonResUpsilon","Calling Destructor");
180 if (fCFManager) delete fCFManager ;
6579f87c 181 if (fnevts) delete fnevts ;
014044ba 182 if (fQAHistList) {fQAHistList->Clear(); delete fQAHistList;}
183}
184//___________________________________________________________________________
185void AliCFMuonResUpsilon::UserCreateOutputObjects() {
186 // UserCreateOutputObjects
187
6579f87c 188 //TH1D *fnevts = new TH1D("fnevts","fnevts",5,0,5); // nevent CINT1B,CMUS1B + PhysSel
014044ba 189
6579f87c 190 PostData(1,fnevts) ;
014044ba 191 PostData(2,fCFManager->GetParticleContainer()) ;
192}
193
194//_________________________________________________
195void AliCFMuonResUpsilon::UserExec(Option_t *)
196{
197
198 // Main loop function
199
200 // variables
201 fIsPhysSelMB=kFALSE; // physics selection : MB
202 fIsPhysSelMUON=kFALSE; // physics selection : MUON
203
204 Double_t containerInput[14] = {0,} ;
205
6579f87c 206 fnevts->Fill(0.5);
014044ba 207 containerInput[0] = 0.5;
208
209 if(!fReadAODData) { // ESD-based ANALYSIS
210
211 Bool_t flag=kFALSE;
212
213 if(fReadMCInfo) {
214 if (!fMCEvent) {
215 Error("UserExec","NO MC EVENT FOUND!");
216 return;
217 }
218
219 // MC part ----------------------------------------------------------------------------------
220 fCFManager->SetMCEventInfo(fMCEvent);
221 AliStack *stack = fMCEvent->Stack();
222 Int_t npart=fMCEvent->GetNumberOfTracks();
223
224 //printf("ESD: npart=%d\n", npart);
225
226 for (Int_t ipart=0; ipart<npart; ipart++) {
227
228 AliMCParticle *mcPart = (AliMCParticle*)fMCEvent->GetTrack(ipart);
229
230 // Mother kinematics
231 TParticle *part = mcPart->Particle();
232 Double_t e = part->Energy();
233 Double_t pz = part->Pz();
234 Double_t rapmc = Rap(e,pz);
235
236 // Selection of the resonance
237 if(!fCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,mcPart)) continue;
238 if(part->Pt()<0 || part->Pt()>1000) continue;
239 if(rapmc<-4 || rapmc>-2.4) continue;
240
241 // Decays kinematics
242 Int_t p0 = part->GetDaughter(0);
243 TParticle *part0 = stack->Particle(p0);
244 AliMCParticle *mcpart0 = new AliMCParticle(part0);
245 Int_t pdg0 = part0->GetPdgCode();
246
247 Int_t p1 = part->GetDaughter(1);
248 TParticle *part1 = stack->Particle(p1);
249 AliMCParticle *mcpart1 = new AliMCParticle(part1);
250 Int_t pdg1 = part1->GetPdgCode();
251
252 Double_t e0 = part0->Energy();
253 Double_t pz0 = part0->Pz();
254 Double_t py0 = part0->Py();
255 Double_t px0 = part0->Px();
256 Double_t charge0 = part0->GetPDG()->Charge()/3;
257 Double_t pt0 = TMath::Sqrt(px0*px0+py0*py0);
258 Double_t theta0 = (180./TMath::Pi())*TMath::ATan2(TMath::Sqrt(px0*px0+py0*py0),pz0);
259 Double_t eta0 = part0->Eta();
260 Double_t vz0 = part0->Vz();
261 Double_t mom0 = part0->P();
262
263 if(!fCFManager->CheckParticleCuts(AliCFManager::kPartAccCuts,mcpart0)) continue;
264 if(pt0<0.5) continue;
265 if(theta0<171. || theta0>178.) continue;
266
267 Double_t e1 = part1->Energy();
268 Double_t pz1 = part1->Pz();
269 Double_t py1 = part1->Py();
270 Double_t px1 = part1->Px();
271 Double_t charge1 = part1->GetPDG()->Charge()/3;
272 Double_t pt1 = TMath::Sqrt(px1*px1+py1*py1);
273 Double_t theta1 = (180./TMath::Pi())*TMath::ATan2(TMath::Sqrt(px1*px1+py1*py1),pz1);
274 Double_t eta1 = part1->Eta();
275 Double_t vz1 = part1->Vz();
276 Double_t mom1 = part1->P();
277
278 if(!fCFManager->CheckParticleCuts(AliCFManager::kPartAccCuts,mcpart1)) continue;
279 if(pt1<0.5) continue;
280 if(theta1<171. || theta1>178.) continue;
281
282 if(pdg0==-13 || pdg1==-13) {
283
284 Double_t ptmc = TMath::Sqrt((px0+px1)*(px0+px1)+(py0+py1)*(py0+py1));
285 Double_t imassmc = Imass(e0,px0,py0,pz0,e1,px1,py1,pz1);
286 if(!imassmc) continue;
287
288 containerInput[1] = rapmc ;
289 containerInput[2] = ptmc ;
290 containerInput[3] = imassmc ;
291 containerInput[4] = 1. ;
292 containerInput[5] = pt0 ;
293 containerInput[5] = pt1 ;
294 containerInput[6] = mom0 ;
295 containerInput[6] = mom1 ;
296 containerInput[7] = 1. ;
297 containerInput[8] = 0. ;
298 containerInput[8] = 0. ;
299 containerInput[9] = charge0+charge1 ;
300 containerInput[10] = eta0;
301 containerInput[10] = eta1;
302 containerInput[11] = theta0;
303 containerInput[11] = theta1;
304 containerInput[12] = vz0;
305 containerInput[12] = vz1;
306 containerInput[13] = 0;
307 containerInput[13] = 0;
308
309 // fill the container at the first step
310 fCFManager->GetParticleContainer()->Fill(containerInput,0); // MC container
311 flag=kTRUE;
312 } // second mu loop
313 } // first mu loop
314 } // end fReadMCInfo
315
316 if(fReadMCInfo && !flag) return;
317 // RECO : ESD part ----------------------------------------------------------------------------------
318
319 // ESD handler
320 AliESDEvent *fESD = 0x0;
321
322 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
323 if( ! esdH) {
324 AliError("Cannot get input event handler");
325 return;
326 }
327 fESD = esdH->GetEvent();
328
329
330 Int_t ntrk=fESD->GetNumberOfMuonTracks();
331
332 Int_t trigfired=-1;
333 Int_t trigside=-1;
334
335 if(!fReadMCInfo) {
336 fIsPhysSelMB=(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kMB));
337 fIsPhysSelMUON=(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kMUON));
338
339 if(fDistinguishTrigClass) {
340 TString trigclass = fESD->GetFiredTriggerClasses();
341 if(trigclass.Contains(fTrigClassMuon)) trigfired = 1;
342 else if(trigclass.Contains(fTrigClassInteraction)) trigfired = 0;
343 for(Int_t i=0;i<3;i++) {
344 if(trigfired==1 && trigclass.Contains(fTrigClassSide[i])) {trigside=i;}
345 if(trigfired==0 && trigclass.Contains(fTrigClassSide[i])) {trigside=i;}
346 }
347 }
348
349 if(trigside==1) { // Beam-Beam event
350 if(trigfired) { // CMUS1-B
351 containerInput[0]=3.5;
6579f87c 352 fnevts->Fill(3.5);
014044ba 353 if(fIsPhysSelMUON) {
354 containerInput[0]=4.5;
6579f87c 355 fnevts->Fill(4.5); // PhysSel
014044ba 356 }
357 } else { // CINT1-B
358 containerInput[0]=1.5;
6579f87c 359 fnevts->Fill(1.5);
014044ba 360 if(fIsPhysSelMB) {
361 containerInput[0]=2.5;
6579f87c 362 fnevts->Fill(2.5); // PhysSel
014044ba 363 }
364 }
365 }
366 }
367 else trigside = 0; // for MC
368
369 for(Int_t j=0; j<ntrk-1; j++) {
370 for(Int_t jj=j+1; jj<ntrk; jj++) {
371 AliESDMuonTrack *mu1 = new AliESDMuonTrack(*(fESD->GetMuonTrack(j)));
372 AliESDMuonTrack *mu2 = new AliESDMuonTrack(*(fESD->GetMuonTrack(jj)));
373
374 Double_t trigCondition = 0;
375 if (mu1->GetMatchTrigger()>=mu2->GetMatchTrigger()) trigCondition = mu1->GetMatchTrigger()+0.1*mu2->GetMatchTrigger();
376 else trigCondition = mu2->GetMatchTrigger()+0.1*mu1->GetMatchTrigger();
377
378 // mu1
379 Double_t zr1 = mu1->Charge();
380 Double_t pxr1 = mu1->Px();
381 Double_t pyr1 = mu1->Py();
382 Double_t pzr1 = mu1->Pz();
383 Double_t ptr1 = TMath::Sqrt(pxr1*pxr1+pyr1*pyr1);
384 Double_t er1 = mu1->E();
385 Double_t momr1 = mu1->P();
386 Double_t vzr1 = mu1->GetZ();
387 Double_t dcar1 = mu1->GetDCA();
388 Double_t rabsr1 = mu1->GetRAtAbsorberEnd();
6579f87c 389 Double_t thetaabsr1 = 180.*(1.-TMath::ATan(rabsr1/505.)/TMath::Pi()); // [deg]
390 Double_t thetaabsr12 = (thetaabsr1*TMath::Pi()/180.)/2.;
391 Double_t tanthetaabsr1 = TMath::Tan(thetaabsr12);
392 Double_t etar1 = -TMath::Log(tanthetaabsr1);
014044ba 393
394 if(Rap(er1,pzr1)<-4 || Rap(er1,pzr1)>-2.5) continue;
395
396 // mu2
397 Double_t zr2 = mu2->Charge();
398 Double_t pxr2 = mu2->Px();
399 Double_t pyr2 = mu2->Py();
400 Double_t pzr2 = mu2->Pz();
401 Double_t ptr2 = TMath::Sqrt(pxr2*pxr2+pyr2*pyr2);
402 Double_t er2 = mu2->E();
403 Double_t momr2 = mu2->P();
404 Double_t vzr2 = mu2->GetZ();
405 Double_t dcar2 = mu2->GetDCA();
406 Double_t rabsr2 = mu2->GetRAtAbsorberEnd();
6579f87c 407 Double_t thetaabsr2 = 180.*(1.-TMath::ATan(rabsr2/505.)/TMath::Pi()); // [deg]
408 Double_t thetaabsr22 = (thetaabsr2*TMath::Pi()/180.)/2.;
409 Double_t tanthetaabsr2 = TMath::Tan(thetaabsr22);
410 Double_t etar2 = -TMath::Log(tanthetaabsr2);
014044ba 411
412 if(Rap(er2,pzr2)<-4 || Rap(er2,pzr2)>-2.5) continue;
413
414 // for MC
415 if(fReadMCInfo) {
416 //mu1
417 rabsr1 = 0;
418 Double_t thetar1 = (180./TMath::Pi())*TMath::ATan2(TMath::Sqrt(pxr1*pxr1+pyr1*pyr1),pzr1);
6579f87c 419 thetaabsr1 = thetar1;
014044ba 420 etar1 = mu1->Eta();
421 //mu2
422 rabsr2 = 0;
423 Double_t thetar2 = (180./TMath::Pi())*TMath::ATan2(TMath::Sqrt(pxr2*pxr2+pyr2*pyr2),pzr2);
6579f87c 424 thetaabsr2 = thetar2;
014044ba 425 etar2 = mu2->Eta();
426 }
427
428 if(TMath::Abs(etar1) > 8 || TMath::Abs(etar2) > 8) continue;
429
430 // dimuon
431 Double_t ptrec = TMath::Sqrt((pxr1+pxr2)*(pxr1+pxr2)+(pyr1+pyr2)*(pyr1+pyr2));
432 Double_t raprec= Rap((er1+er2),(pzr1+pzr2));
433 Double_t imassrec = Imass(er1,pxr1,pyr1,pzr1,er2,pxr2,pyr2,pzr2);
434 if(!imassrec) continue;
435
436 containerInput[1] = raprec ;
437 containerInput[2] = ptrec ;
438 containerInput[3] = imassrec ;
439 containerInput[4] = trigCondition+0.05 ;
440 containerInput[5] = ptr1 ;
441 containerInput[5] = ptr2 ;
442 containerInput[6] = momr1;
443 containerInput[6] = momr2;
444 containerInput[7] = trigside;
445 containerInput[8] = rabsr1;
446 containerInput[8] = rabsr2;
447 containerInput[9] = zr1 + zr2;
448 containerInput[10] = etar1;
449 containerInput[10] = etar2;
6579f87c 450 containerInput[11] = thetaabsr1;
451 containerInput[11] = thetaabsr2;
014044ba 452 containerInput[12] = vzr1;
453 containerInput[12] = vzr2;
454 containerInput[13] = dcar1;
455 containerInput[13] = dcar2;
456
457 if (trigfired==1 && trigside==1) fCFManager->GetParticleContainer()->Fill(containerInput,4);
458 if (trigfired==0 && trigside==1) fCFManager->GetParticleContainer()->Fill(containerInput,1);
459 if(fReadMCInfo) fCFManager->GetParticleContainer()->Fill(containerInput,1); // Rec container
460 } // second mu loop
461 } // first mu loop
462 } // end ESD-based ANALYSIS
463
464 else { // AOD-based ANALYSIS
465
466 Bool_t flag=kTRUE;
467
468 // AOD handler
469 AliAODEvent *fAOD = 0x0;
470
471 AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
014044ba 472 if( ! aodH) {
473 AliError("Cannot get input event handler");
474 return;
475 }
e199f3f0 476 fAOD = aodH->GetEvent();
014044ba 477
478 // MC part -----------------------------------------------------------------------------------
479
480 if(fReadMCInfo) {
481 TClonesArray *mcArr = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
a29eedf0 482 if( ! mcArr) {
483 AliError("Cannot get MC innf in AOD MC branch");
484 return;
485 }
014044ba 486 Int_t npart=mcArr->GetEntries();
487
488 for(Int_t i=0; i<npart; i++) {
489 AliAODMCParticle *mctrack = (AliAODMCParticle*) mcArr->At(i);
490 // select resonance
491 if(mctrack->GetPdgCode()!=fPDG) continue;
492 // cuts on resonance
493 if(!(mctrack->Pt()>0 && mctrack->Pt()<1000) || !(Rap(mctrack->E(),mctrack->Pz())>-4 && Rap(mctrack->E(),mctrack->Pz())<-2.4)) continue;
494 Int_t daug0 = mctrack->GetDaughter(0);
495 Int_t daug1 = mctrack->GetDaughter(1);
496 // daughter1
497 AliAODMCParticle *mcdaug0 = (AliAODMCParticle*) mcArr->At(daug0);
498 Double_t pt0 = mcdaug0->Pt();
499 Double_t mom0 = mcdaug0->P();
500 Double_t theta0 = (180./TMath::Pi())*mcdaug0->Theta();
501 Double_t eta0 = mcdaug0->Eta();
502 Int_t charge0 = (Int_t) mcdaug0->Charge()/3;
503 Double_t vz0 = mcdaug0->Zv();
504 if(!(pt0>0.5) || !(theta0>171. && theta0<178.)) continue;
505 // daughter2
506 AliAODMCParticle *mcdaug1 = (AliAODMCParticle*) mcArr->At(daug1);
507 Double_t pt1 = mcdaug1->Pt();
508 Double_t mom1 = mcdaug1->P();
509 Double_t theta1 = (180./TMath::Pi())*mcdaug1->Theta();
510 Double_t eta1 = mcdaug1->Eta();
511 Int_t charge1 = (Int_t) mcdaug1->Charge()/3;
512 Double_t vz1 = mcdaug1->Zv();
513 if(!(pt1>0.5) || !(theta1>171. && theta1<178.)) continue;
514
515 if (TMath::Abs(mcdaug0->GetPdgCode())!=13 || TMath::Abs(mcdaug1->GetPdgCode())!=13) continue;
516
517 Double_t rapmc = Rap(mctrack->E(), mctrack->Pz());
518 Double_t ptmc = mctrack->Pt();
519 Double_t imassmc = mctrack->M();
520 if(!imassmc) continue;
521
522 containerInput[1] = rapmc ;
523 containerInput[2] = ptmc ;
524 containerInput[3] = imassmc ;
525 containerInput[4] = 1. ;
526 containerInput[5] = pt0 ;
527 containerInput[5] = pt1 ;
528 containerInput[6] = mom0 ;
529 containerInput[6] = mom1 ;
530 containerInput[7] = 1. ;
531 containerInput[8] = 0. ;
532 containerInput[8] = 0. ;
533 containerInput[9] = charge0+charge1 ;
534 containerInput[10] = eta0;
535 containerInput[10] = eta1;
536 containerInput[11] = theta0;
537 containerInput[11] = theta1;
538 containerInput[12] = vz0;
539 containerInput[12] = vz1;
540 containerInput[13] = 0;
541 containerInput[13] = 0;
542
543 fCFManager->GetParticleContainer()->Fill(containerInput,0); // MC container
544 flag=kTRUE;
545 }
546 }
547
548if(fReadMCInfo && !flag) return;
549
550 // RECO : AOD part ----------------------------------------------------------------------------------
551
552
553 Int_t trigfired=-1;
554 Int_t trigside=-1;
555
556 Int_t ntrk = fAOD->GetNumberOfTracks();
557
558 if(!fReadMCInfo) {
559 fIsPhysSelMB=(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kMB));
560 fIsPhysSelMUON=(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kMUON));
561
562 if(fDistinguishTrigClass) {
563 TString trigclass = fAOD->GetFiredTriggerClasses();
564 if(trigclass.Contains(fTrigClassMuon)) trigfired = 1;
565 else if(trigclass.Contains(fTrigClassInteraction)) trigfired = 0;
566 for(Int_t i=0;i<3;i++) {
567 if(trigfired==1 && trigclass.Contains(fTrigClassSide[i].Data())) {trigside=i;}
568 if(trigfired==0 && trigclass.Contains(fTrigClassSide[i].Data())) {trigside=i;}
569 }
570 }
571
572 if(trigside==1) { // Beam-Beam event
573 if(trigfired) { // CMUS1-B
574 containerInput[0]=3.5;
6579f87c 575 fnevts->Fill(3.5);
014044ba 576 if(fIsPhysSelMUON) {
577 containerInput[0]=4.5;
6579f87c 578 fnevts->Fill(4.5); // PhysSel
014044ba 579 }
580 } else { // CINT1-B
581 containerInput[0]=1.5;
6579f87c 582 fnevts->Fill(1.5);
014044ba 583 if(fIsPhysSelMB) {
584 containerInput[0]=2.5;
6579f87c 585 fnevts->Fill(2.5); // PhysSel
014044ba 586 }
587 }
588 }
589 }
590 else trigside = 0; // for MC
591
592
593 for(Int_t j=0; j<ntrk; j++) {
594 for(Int_t jj=j+1; jj<ntrk; jj++) {
595 AliAODTrack *mu1 = fAOD->GetTrack(j);
596 if(!mu1->IsMuonTrack() || !(mu1->Y()>-4 && mu1->Y()<-2.5)) continue;
597 AliAODTrack *mu2 = fAOD->GetTrack(jj);
598 if(!mu2->IsMuonTrack() || !(mu2->Y()>-4 && mu2->Y()<-2.5)) continue;
599
600 Double_t trigCondition=0;
601 if (mu1->GetMatchTrigger()>=mu2->GetMatchTrigger()) trigCondition = mu1->GetMatchTrigger()+0.1*mu2->GetMatchTrigger();
602 else trigCondition = mu2->GetMatchTrigger()+0.1*mu1->GetMatchTrigger();
603
604 // mu1
605 Double_t zr1 = mu1->Charge();
606 Double_t pxr1 = mu1->Px();
607 Double_t pyr1 = mu1->Py();
608 Double_t pzr1 = mu1->Pz();
609 Double_t ptr1 = TMath::Sqrt(pxr1*pxr1+pyr1*pyr1);
610 Double_t er1 = mu1->E();
611 Double_t momr1 = mu1->P();
612 Double_t vzr1 = mu1->Zv();
613 Double_t dcar1 = mu1->DCA();
614 Double_t rabsr1 = mu1->GetRAtAbsorberEnd();
6579f87c 615 Double_t thetaabsr1 = 180.*(1.-TMath::ATan(rabsr1/505.)/TMath::Pi()); // [deg]
616 Double_t thetaabsr12 = (thetaabsr1*TMath::Pi()/180.)/2.;
617 Double_t tanthetaabsr1 = TMath::Tan(thetaabsr12);
618 Double_t etar1 = -TMath::Log(tanthetaabsr1);
014044ba 619
620 // mu2
621 Double_t zr2 = mu2->Charge();
622 Double_t pxr2 = mu2->Px();
623 Double_t pyr2 = mu2->Py();
624 Double_t pzr2 = mu2->Pz();
625 Double_t ptr2 = TMath::Sqrt(pxr2*pxr2+pyr2*pyr2);
626 Double_t er2 = mu2->E();
627 Double_t momr2 = mu2->P();
628 Double_t vzr2 = mu2->Zv();
629 Double_t dcar2 = mu2->DCA();
630 Double_t rabsr2 = mu1->GetRAtAbsorberEnd();
6579f87c 631 Double_t thetaabsr2 = 180.*(1.-TMath::ATan(rabsr2/505.)/TMath::Pi()); // [deg]
632 Double_t thetaabsr22 = (thetaabsr2*TMath::Pi()/180.)/2.;
633 Double_t tanthetaabsr2 = TMath::Tan(thetaabsr22);
634 Double_t etar2 = -TMath::Log(tanthetaabsr2);
014044ba 635
636 // for MC
637 if(fReadMCInfo) {
638 //mu1
639 rabsr1 = 0;
640 Double_t thetar1 = (180./TMath::Pi())*TMath::ATan2(TMath::Sqrt(pxr1*pxr1+pyr1*pyr1),pzr1);
6579f87c 641 thetaabsr1 = thetar1;
014044ba 642 etar1 = mu1->Eta();
643 //mu2
644 rabsr2 = 0;
645 Double_t thetar2 = (180./TMath::Pi())*TMath::ATan2(TMath::Sqrt(pxr2*pxr2+pyr2*pyr2),pzr2);
6579f87c 646 thetaabsr2 = thetar2;
014044ba 647 etar2 = mu2->Eta();
648 }
649
650 if(TMath::Abs(etar1) > 8 || TMath::Abs(etar2) > 8) continue;
651
652 // dimuon
653 Double_t ptrec = TMath::Sqrt((pxr1+pxr2)*(pxr1+pxr2)+(pyr1+pyr2)*(pyr1+pyr2));
654 Double_t raprec= Rap((er1+er2),(pzr1+pzr2));
655 Double_t imassrec = Imass(er1,pxr1,pyr1,pzr1,er2,pxr2,pyr2,pzr2);
656 if(!imassrec) continue;
657
658 containerInput[1] = raprec ;
659 containerInput[2] = ptrec ;
660 containerInput[3] = imassrec ;
661 containerInput[4] = trigCondition+0.05 ;
662 containerInput[5] = ptr1 ;
663 containerInput[5] = ptr2 ;
664 containerInput[6] = momr1;
665 containerInput[6] = momr2;
666 containerInput[7] = trigside;
667 containerInput[8] = rabsr1;
668 containerInput[8] = rabsr2;
669 containerInput[9] = zr1 + zr2;
670 containerInput[10] = etar1;
671 containerInput[10] = etar2;
6579f87c 672 containerInput[11] = thetaabsr1;
673 containerInput[11] = thetaabsr2;
014044ba 674 containerInput[12] = vzr1;
675 containerInput[12] = vzr2;
676 containerInput[13] = dcar1;
677 containerInput[13] = dcar2;
678
679 if (trigfired==1 && trigside==1) fCFManager->GetParticleContainer()->Fill(containerInput,4);
680 if (trigfired==0 && trigside==1) fCFManager->GetParticleContainer()->Fill(containerInput,1);
681 if(fReadMCInfo) fCFManager->GetParticleContainer()->Fill(containerInput,1); // Rec container
682
683 } // second mu loop
684 } // first mu loop
685 } // end AOD-based ANALYSIS
686
687 // end User analysis loop
688}
689//________________________________________________________________________
690Double_t AliCFMuonResUpsilon::Imass(Float_t e1, Float_t px1, Float_t py1, Float_t pz1,
691 Float_t e2, Float_t px2, Float_t py2, Float_t pz2) const
692{
693// invariant mass calculation
694 Float_t imass = (e1+e2)*(e1+e2)-((px1+px2)*(px1+px2)+(py1+py2)*(py1+py2)+(pz1+pz2)*(pz1+pz2));
695 if(imass<0) return 0;
696 Float_t imassrec = TMath::Sqrt((e1+e2)*(e1+e2)-((px1+px2)*(px1+px2)+(py1+py2)*(py1+py2)+(pz1+pz2)*(pz1+pz2)));
697 return imassrec;
698}
699//________________________________________________________________________
700Double_t AliCFMuonResUpsilon::Rap(Float_t e, Float_t pz) const
701{
702// calculate rapidity
703 Float_t rap;
704 if(e!=pz){
705 rap = 0.5*TMath::Log((e+pz)/(e-pz));
706 return rap;
707 }
708 else{
709 rap = -200;
710 return rap;
711 }
712}
713//________________________________________________________________________
714Double_t AliCFMuonResUpsilon::CostCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
715Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2,
716Double_t Energy)
717{
6579f87c 718 // CS angle
719 TLorentzVector pMu1CM, pMu2CM, pProjCM, pTargCM, pDimuCM; // In the CM. frame
720 TLorentzVector pMu1Dimu, pMu2Dimu, pProjDimu, pTargDimu; // In the dimuon rest frame
014044ba 721 TVector3 beta,zaxisCS;
722 Double_t mp=0.93827231;
723 //
724 // --- Fill the Lorentz vector for projectile and target in the CM frame
725 //
6579f87c 726 pProjCM.SetPxPyPzE(0.,0.,-Energy,TMath::Sqrt(Energy*Energy+mp*mp));
727 pTargCM.SetPxPyPzE(0.,0.,Energy,TMath::Sqrt(Energy*Energy+mp*mp));
014044ba 728 //
729 // --- Get the muons parameters in the CM frame
730 //
6579f87c 731 pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
732 pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
014044ba 733 //
734 // --- Obtain the dimuon parameters in the CM frame
735 //
6579f87c 736 pDimuCM=pMu1CM+pMu2CM;
014044ba 737 //
738 // --- Translate the dimuon parameters in the dimuon rest frame
739 //
6579f87c 740 beta=(-1./pDimuCM.E())*pDimuCM.Vect();
741 pMu1Dimu=pMu1CM;
742 pMu2Dimu=pMu2CM;
743 pProjDimu=pProjCM;
744 pTargDimu=pTargCM;
745 pMu1Dimu.Boost(beta);
746 pMu2Dimu.Boost(beta);
747 pProjDimu.Boost(beta);
748 pTargDimu.Boost(beta);
014044ba 749 //
750 // --- Determine the z axis for the CS angle
751 //
6579f87c 752 zaxisCS=(((pProjDimu.Vect()).Unit())-((pTargDimu.Vect()).Unit())).Unit();
014044ba 753 //
754 // --- Determine the CS angle (angle between mu+ and the z axis defined above)
755 Double_t cost;
756
6579f87c 757 if(charge1>0) cost = zaxisCS.Dot((pMu1Dimu.Vect()).Unit());
758 else cost = zaxisCS.Dot((pMu2Dimu.Vect()).Unit());
014044ba 759 return cost;
760}
761//________________________________________________________________________
762
763//________________________________________________________________________
764Double_t AliCFMuonResUpsilon::CostHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
765Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
766{
6579f87c 767 // Helicity
768 TLorentzVector pMu1CM, pMu2CM, pDimuCM; // In the CM frame
769 TLorentzVector pMu1Dimu, pMu2Dimu; // In the dimuon rest frame
014044ba 770 TVector3 beta,zaxisCS;
771 //
772 // --- Get the muons parameters in the CM frame
773 //
6579f87c 774 pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
775 pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
014044ba 776 //
777 // --- Obtain the dimuon parameters in the CM frame
778 //
6579f87c 779 pDimuCM=pMu1CM+pMu2CM;
014044ba 780 //
781 // --- Translate the muon parameters in the dimuon rest frame
782 //
6579f87c 783 beta=(-1./pDimuCM.E())*pDimuCM.Vect();
784 pMu1Dimu=pMu1CM;
785 pMu2Dimu=pMu2CM;
786 pMu1Dimu.Boost(beta);
787 pMu2Dimu.Boost(beta);
014044ba 788 //
789 // --- Determine the z axis for the calculation of the polarization angle (i.e. the direction of the dimuon in the CM system)
790 //
791 TVector3 zaxis;
6579f87c 792 zaxis=(pDimuCM.Vect()).Unit();
014044ba 793
794 // --- Calculation of the polarization angle (angle between mu+ and the z axis defined above)
795 Double_t cost;
6579f87c 796 if(charge1>0) cost = zaxis.Dot((pMu1Dimu.Vect()).Unit());
797 else cost = zaxis.Dot((pMu2Dimu.Vect()).Unit());
014044ba 798
799 return cost;
800}
801//________________________________________________________________________
802
803//________________________________________________________________________
804Double_t AliCFMuonResUpsilon::PhiCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
805Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2,
806Double_t Energy)
807{
6579f87c 808 // CS phi
809 TLorentzVector pMu1CM, pMu2CM, pProjCM, pTargCM, pDimuCM; // In the CM frame
810 TLorentzVector pMu1Dimu, pMu2Dimu, pProjDimu, pTargDimu; // In the dimuon rest frame
014044ba 811 TVector3 beta,yaxisCS, xaxisCS, zaxisCS;
812 Double_t mp=0.93827231;
813 //
814 // --- Fill the Lorentz vector for projectile and target in the CM frame
815 //
6579f87c 816 pProjCM.SetPxPyPzE(0.,0.,-Energy,TMath::Sqrt(Energy*Energy+mp*mp));
817 pTargCM.SetPxPyPzE(0.,0.,Energy,TMath::Sqrt(Energy*Energy+mp*mp));
014044ba 818 //
819 // --- Get the muons parameters in the CM frame
820 //
6579f87c 821 pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
822 pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
014044ba 823 //
824 // --- Obtain the dimuon parameters in the CM frame
825 //
6579f87c 826 pDimuCM=pMu1CM+pMu2CM;
014044ba 827 //
828 // --- Translate the dimuon parameters in the dimuon rest frame
829 //
6579f87c 830 beta=(-1./pDimuCM.E())*pDimuCM.Vect();
831 pMu1Dimu=pMu1CM;
832 pMu2Dimu=pMu2CM;
833 pProjDimu=pProjCM;
834 pTargDimu=pTargCM;
835 pMu1Dimu.Boost(beta);
836 pMu2Dimu.Boost(beta);
837 pProjDimu.Boost(beta);
838 pTargDimu.Boost(beta);
014044ba 839 //
840 // --- Determine the z axis for the CS angle
841 //
6579f87c 842 zaxisCS=(((pProjDimu.Vect()).Unit())-((pTargDimu.Vect()).Unit())).Unit();
843 yaxisCS=(((pProjDimu.Vect()).Unit()).Cross((pTargDimu.Vect()).Unit())).Unit();
014044ba 844 xaxisCS=(yaxisCS.Cross(zaxisCS)).Unit();
845
846 Double_t phi;
6579f87c 847 if(charge1>0) phi = TMath::ATan2((pMu1Dimu.Vect()).Dot(yaxisCS),((pMu1Dimu.Vect()).Dot(xaxisCS)));
848 else phi = TMath::ATan2((pMu2Dimu.Vect()).Dot(yaxisCS),((pMu2Dimu.Vect()).Dot(xaxisCS)));
014044ba 849
850 return phi;
851}
852//________________________________________________________________________
853
854//________________________________________________________________________
855Double_t AliCFMuonResUpsilon::PhiHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
856Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2, Double_t Energy)
857{
6579f87c 858 // Helicity phi
859 TLorentzVector pMu1CM, pMu2CM, pProjCM, pTargCM, pDimuCM; // In the CM frame
860 TLorentzVector pMu1Dimu, pMu2Dimu, pProjDimu, pTargDimu; // In the dimuon rest frame
014044ba 861 TVector3 beta,xaxis,yaxis,zaxis;
862 Double_t mp=0.93827231;
863
864 //
865 // --- Get the muons parameters in the CM frame
866 //
6579f87c 867 pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
868 pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
014044ba 869 //
870 // --- Obtain the dimuon parameters in the CM frame
871 //
6579f87c 872 pDimuCM=pMu1CM+pMu2CM;
014044ba 873 //
874 // --- Translate the muon parameters in the dimuon rest frame
875 //
6579f87c 876 zaxis=(pDimuCM.Vect()).Unit();
014044ba 877
6579f87c 878 beta=(-1./pDimuCM.E())*pDimuCM.Vect();
014044ba 879
6579f87c 880 pProjCM.SetPxPyPzE(0.,0.,-Energy,TMath::Sqrt(Energy*Energy+mp*mp));
881 pTargCM.SetPxPyPzE(0.,0.,Energy,TMath::Sqrt(Energy*Energy+mp*mp));
014044ba 882
6579f87c 883 pProjDimu=pProjCM;
884 pTargDimu=pTargCM;
014044ba 885
6579f87c 886 pProjDimu.Boost(beta);
887 pTargDimu.Boost(beta);
014044ba 888
6579f87c 889 yaxis=((pProjDimu.Vect()).Cross(pTargDimu.Vect())).Unit();
014044ba 890 xaxis=(yaxis.Cross(zaxis)).Unit();
891
6579f87c 892 pMu1Dimu=pMu1CM;
893 pMu2Dimu=pMu2CM;
894 pMu1Dimu.Boost(beta);
895 pMu2Dimu.Boost(beta);
014044ba 896
897 Double_t phi;
6579f87c 898 if(charge1>0) phi = TMath::ATan2((pMu1Dimu.Vect()).Dot(yaxis),(pMu1Dimu.Vect()).Dot(xaxis));
899 else phi = TMath::ATan2((pMu2Dimu.Vect()).Dot(yaxis),(pMu2Dimu.Vect()).Dot(xaxis));
014044ba 900
901 return phi;
902}
903
904//________________________________________________________________________
905void AliCFMuonResUpsilon::Terminate(Option_t *)
906{
907 // draw result of the Invariant mass MC and ESD
908
e199f3f0 909/*
014044ba 910 TH1D *h1 = dynamic_cast<TH1D*>(GetOutputData(1));
911 AliCFContainer *cont = dynamic_cast<AliCFContainer*> (GetOutputData(2));
912
913 TH1D *kmass = cont->ShowProjection(3,0);
914 TH1D *rmass = cont->ShowProjection(3,1);
915 TH1D *mmass = cont->ShowProjection(3,4);
916
917 TCanvas *c1 = new TCanvas("AliCFMuonResUpsilon","UPSILON Container",0,0,800,800);
918 c1->Divide(2,2);
919 c1->cd(1);
920 kmass->Draw("HIST");
921 c1->cd(2);
922 rmass->Draw("HIST");
923 c1->cd(3);
924 mmass->Draw("HIST");
925 c1->cd(4);
926 h1->Draw("HIST");
e199f3f0 927 */
014044ba 928}
929//________________________________________________________________________
930
931#endif