]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG/muon/AliCFMuonResUpsilon.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fix Coverity 24835
[u/mrichter/AliRoot.git] / PWG / 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++) {
f15c1f69 595 AliAODTrack *mu1 = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(j));
596 if(!mu1) AliFatal("Not a standard AOD");
014044ba 597 if(!mu1->IsMuonTrack() || !(mu1->Y()>-4 && mu1->Y()<-2.5)) continue;
f15c1f69 598 AliAODTrack *mu2 = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(jj));
599 if(!mu2) AliFatal("Not a standard AOD");
014044ba 600 if(!mu2->IsMuonTrack() || !(mu2->Y()>-4 && mu2->Y()<-2.5)) continue;
601
602 Double_t trigCondition=0;
603 if (mu1->GetMatchTrigger()>=mu2->GetMatchTrigger()) trigCondition = mu1->GetMatchTrigger()+0.1*mu2->GetMatchTrigger();
604 else trigCondition = mu2->GetMatchTrigger()+0.1*mu1->GetMatchTrigger();
605
606 // mu1
607 Double_t zr1 = mu1->Charge();
608 Double_t pxr1 = mu1->Px();
609 Double_t pyr1 = mu1->Py();
610 Double_t pzr1 = mu1->Pz();
611 Double_t ptr1 = TMath::Sqrt(pxr1*pxr1+pyr1*pyr1);
612 Double_t er1 = mu1->E();
613 Double_t momr1 = mu1->P();
614 Double_t vzr1 = mu1->Zv();
615 Double_t dcar1 = mu1->DCA();
616 Double_t rabsr1 = mu1->GetRAtAbsorberEnd();
6579f87c 617 Double_t thetaabsr1 = 180.*(1.-TMath::ATan(rabsr1/505.)/TMath::Pi()); // [deg]
618 Double_t thetaabsr12 = (thetaabsr1*TMath::Pi()/180.)/2.;
619 Double_t tanthetaabsr1 = TMath::Tan(thetaabsr12);
620 Double_t etar1 = -TMath::Log(tanthetaabsr1);
014044ba 621
622 // mu2
623 Double_t zr2 = mu2->Charge();
624 Double_t pxr2 = mu2->Px();
625 Double_t pyr2 = mu2->Py();
626 Double_t pzr2 = mu2->Pz();
627 Double_t ptr2 = TMath::Sqrt(pxr2*pxr2+pyr2*pyr2);
628 Double_t er2 = mu2->E();
629 Double_t momr2 = mu2->P();
630 Double_t vzr2 = mu2->Zv();
631 Double_t dcar2 = mu2->DCA();
632 Double_t rabsr2 = mu1->GetRAtAbsorberEnd();
6579f87c 633 Double_t thetaabsr2 = 180.*(1.-TMath::ATan(rabsr2/505.)/TMath::Pi()); // [deg]
634 Double_t thetaabsr22 = (thetaabsr2*TMath::Pi()/180.)/2.;
635 Double_t tanthetaabsr2 = TMath::Tan(thetaabsr22);
636 Double_t etar2 = -TMath::Log(tanthetaabsr2);
014044ba 637
638 // for MC
639 if(fReadMCInfo) {
640 //mu1
641 rabsr1 = 0;
642 Double_t thetar1 = (180./TMath::Pi())*TMath::ATan2(TMath::Sqrt(pxr1*pxr1+pyr1*pyr1),pzr1);
6579f87c 643 thetaabsr1 = thetar1;
014044ba 644 etar1 = mu1->Eta();
645 //mu2
646 rabsr2 = 0;
647 Double_t thetar2 = (180./TMath::Pi())*TMath::ATan2(TMath::Sqrt(pxr2*pxr2+pyr2*pyr2),pzr2);
6579f87c 648 thetaabsr2 = thetar2;
014044ba 649 etar2 = mu2->Eta();
650 }
651
652 if(TMath::Abs(etar1) > 8 || TMath::Abs(etar2) > 8) continue;
653
654 // dimuon
655 Double_t ptrec = TMath::Sqrt((pxr1+pxr2)*(pxr1+pxr2)+(pyr1+pyr2)*(pyr1+pyr2));
656 Double_t raprec= Rap((er1+er2),(pzr1+pzr2));
657 Double_t imassrec = Imass(er1,pxr1,pyr1,pzr1,er2,pxr2,pyr2,pzr2);
658 if(!imassrec) continue;
659
660 containerInput[1] = raprec ;
661 containerInput[2] = ptrec ;
662 containerInput[3] = imassrec ;
663 containerInput[4] = trigCondition+0.05 ;
664 containerInput[5] = ptr1 ;
665 containerInput[5] = ptr2 ;
666 containerInput[6] = momr1;
667 containerInput[6] = momr2;
668 containerInput[7] = trigside;
669 containerInput[8] = rabsr1;
670 containerInput[8] = rabsr2;
671 containerInput[9] = zr1 + zr2;
672 containerInput[10] = etar1;
673 containerInput[10] = etar2;
6579f87c 674 containerInput[11] = thetaabsr1;
675 containerInput[11] = thetaabsr2;
014044ba 676 containerInput[12] = vzr1;
677 containerInput[12] = vzr2;
678 containerInput[13] = dcar1;
679 containerInput[13] = dcar2;
680
681 if (trigfired==1 && trigside==1) fCFManager->GetParticleContainer()->Fill(containerInput,4);
682 if (trigfired==0 && trigside==1) fCFManager->GetParticleContainer()->Fill(containerInput,1);
683 if(fReadMCInfo) fCFManager->GetParticleContainer()->Fill(containerInput,1); // Rec container
684
685 } // second mu loop
686 } // first mu loop
687 } // end AOD-based ANALYSIS
688
689 // end User analysis loop
690}
691//________________________________________________________________________
692Double_t AliCFMuonResUpsilon::Imass(Float_t e1, Float_t px1, Float_t py1, Float_t pz1,
693 Float_t e2, Float_t px2, Float_t py2, Float_t pz2) const
694{
695// invariant mass calculation
696 Float_t imass = (e1+e2)*(e1+e2)-((px1+px2)*(px1+px2)+(py1+py2)*(py1+py2)+(pz1+pz2)*(pz1+pz2));
697 if(imass<0) return 0;
698 Float_t imassrec = TMath::Sqrt((e1+e2)*(e1+e2)-((px1+px2)*(px1+px2)+(py1+py2)*(py1+py2)+(pz1+pz2)*(pz1+pz2)));
699 return imassrec;
700}
701//________________________________________________________________________
702Double_t AliCFMuonResUpsilon::Rap(Float_t e, Float_t pz) const
703{
704// calculate rapidity
705 Float_t rap;
706 if(e!=pz){
707 rap = 0.5*TMath::Log((e+pz)/(e-pz));
708 return rap;
709 }
710 else{
711 rap = -200;
712 return rap;
713 }
714}
715//________________________________________________________________________
716Double_t AliCFMuonResUpsilon::CostCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
717Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2,
718Double_t Energy)
719{
6579f87c 720 // CS angle
721 TLorentzVector pMu1CM, pMu2CM, pProjCM, pTargCM, pDimuCM; // In the CM. frame
722 TLorentzVector pMu1Dimu, pMu2Dimu, pProjDimu, pTargDimu; // In the dimuon rest frame
014044ba 723 TVector3 beta,zaxisCS;
724 Double_t mp=0.93827231;
725 //
726 // --- Fill the Lorentz vector for projectile and target in the CM frame
727 //
6579f87c 728 pProjCM.SetPxPyPzE(0.,0.,-Energy,TMath::Sqrt(Energy*Energy+mp*mp));
729 pTargCM.SetPxPyPzE(0.,0.,Energy,TMath::Sqrt(Energy*Energy+mp*mp));
014044ba 730 //
731 // --- Get the muons parameters in the CM frame
732 //
6579f87c 733 pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
734 pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
014044ba 735 //
736 // --- Obtain the dimuon parameters in the CM frame
737 //
6579f87c 738 pDimuCM=pMu1CM+pMu2CM;
014044ba 739 //
740 // --- Translate the dimuon parameters in the dimuon rest frame
741 //
6579f87c 742 beta=(-1./pDimuCM.E())*pDimuCM.Vect();
743 pMu1Dimu=pMu1CM;
744 pMu2Dimu=pMu2CM;
745 pProjDimu=pProjCM;
746 pTargDimu=pTargCM;
747 pMu1Dimu.Boost(beta);
748 pMu2Dimu.Boost(beta);
749 pProjDimu.Boost(beta);
750 pTargDimu.Boost(beta);
014044ba 751 //
752 // --- Determine the z axis for the CS angle
753 //
6579f87c 754 zaxisCS=(((pProjDimu.Vect()).Unit())-((pTargDimu.Vect()).Unit())).Unit();
014044ba 755 //
756 // --- Determine the CS angle (angle between mu+ and the z axis defined above)
757 Double_t cost;
758
6579f87c 759 if(charge1>0) cost = zaxisCS.Dot((pMu1Dimu.Vect()).Unit());
760 else cost = zaxisCS.Dot((pMu2Dimu.Vect()).Unit());
014044ba 761 return cost;
762}
763//________________________________________________________________________
764
765//________________________________________________________________________
766Double_t AliCFMuonResUpsilon::CostHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
767Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
768{
6579f87c 769 // Helicity
770 TLorentzVector pMu1CM, pMu2CM, pDimuCM; // In the CM frame
771 TLorentzVector pMu1Dimu, pMu2Dimu; // In the dimuon rest frame
014044ba 772 TVector3 beta,zaxisCS;
773 //
774 // --- Get the muons parameters in the CM frame
775 //
6579f87c 776 pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
777 pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
014044ba 778 //
779 // --- Obtain the dimuon parameters in the CM frame
780 //
6579f87c 781 pDimuCM=pMu1CM+pMu2CM;
014044ba 782 //
783 // --- Translate the muon parameters in the dimuon rest frame
784 //
6579f87c 785 beta=(-1./pDimuCM.E())*pDimuCM.Vect();
786 pMu1Dimu=pMu1CM;
787 pMu2Dimu=pMu2CM;
788 pMu1Dimu.Boost(beta);
789 pMu2Dimu.Boost(beta);
014044ba 790 //
791 // --- Determine the z axis for the calculation of the polarization angle (i.e. the direction of the dimuon in the CM system)
792 //
793 TVector3 zaxis;
6579f87c 794 zaxis=(pDimuCM.Vect()).Unit();
014044ba 795
796 // --- Calculation of the polarization angle (angle between mu+ and the z axis defined above)
797 Double_t cost;
6579f87c 798 if(charge1>0) cost = zaxis.Dot((pMu1Dimu.Vect()).Unit());
799 else cost = zaxis.Dot((pMu2Dimu.Vect()).Unit());
014044ba 800
801 return cost;
802}
803//________________________________________________________________________
804
805//________________________________________________________________________
806Double_t AliCFMuonResUpsilon::PhiCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
807Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2,
808Double_t Energy)
809{
6579f87c 810 // CS phi
811 TLorentzVector pMu1CM, pMu2CM, pProjCM, pTargCM, pDimuCM; // In the CM frame
812 TLorentzVector pMu1Dimu, pMu2Dimu, pProjDimu, pTargDimu; // In the dimuon rest frame
014044ba 813 TVector3 beta,yaxisCS, xaxisCS, zaxisCS;
814 Double_t mp=0.93827231;
815 //
816 // --- Fill the Lorentz vector for projectile and target in the CM frame
817 //
6579f87c 818 pProjCM.SetPxPyPzE(0.,0.,-Energy,TMath::Sqrt(Energy*Energy+mp*mp));
819 pTargCM.SetPxPyPzE(0.,0.,Energy,TMath::Sqrt(Energy*Energy+mp*mp));
014044ba 820 //
821 // --- Get the muons parameters in the CM frame
822 //
6579f87c 823 pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
824 pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
014044ba 825 //
826 // --- Obtain the dimuon parameters in the CM frame
827 //
6579f87c 828 pDimuCM=pMu1CM+pMu2CM;
014044ba 829 //
830 // --- Translate the dimuon parameters in the dimuon rest frame
831 //
6579f87c 832 beta=(-1./pDimuCM.E())*pDimuCM.Vect();
833 pMu1Dimu=pMu1CM;
834 pMu2Dimu=pMu2CM;
835 pProjDimu=pProjCM;
836 pTargDimu=pTargCM;
837 pMu1Dimu.Boost(beta);
838 pMu2Dimu.Boost(beta);
839 pProjDimu.Boost(beta);
840 pTargDimu.Boost(beta);
014044ba 841 //
842 // --- Determine the z axis for the CS angle
843 //
6579f87c 844 zaxisCS=(((pProjDimu.Vect()).Unit())-((pTargDimu.Vect()).Unit())).Unit();
845 yaxisCS=(((pProjDimu.Vect()).Unit()).Cross((pTargDimu.Vect()).Unit())).Unit();
014044ba 846 xaxisCS=(yaxisCS.Cross(zaxisCS)).Unit();
847
848 Double_t phi;
6579f87c 849 if(charge1>0) phi = TMath::ATan2((pMu1Dimu.Vect()).Dot(yaxisCS),((pMu1Dimu.Vect()).Dot(xaxisCS)));
850 else phi = TMath::ATan2((pMu2Dimu.Vect()).Dot(yaxisCS),((pMu2Dimu.Vect()).Dot(xaxisCS)));
014044ba 851
852 return phi;
853}
854//________________________________________________________________________
855
856//________________________________________________________________________
857Double_t AliCFMuonResUpsilon::PhiHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
858Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2, Double_t Energy)
859{
6579f87c 860 // Helicity phi
861 TLorentzVector pMu1CM, pMu2CM, pProjCM, pTargCM, pDimuCM; // In the CM frame
862 TLorentzVector pMu1Dimu, pMu2Dimu, pProjDimu, pTargDimu; // In the dimuon rest frame
014044ba 863 TVector3 beta,xaxis,yaxis,zaxis;
864 Double_t mp=0.93827231;
865
866 //
867 // --- Get the muons parameters in the CM frame
868 //
6579f87c 869 pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
870 pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
014044ba 871 //
872 // --- Obtain the dimuon parameters in the CM frame
873 //
6579f87c 874 pDimuCM=pMu1CM+pMu2CM;
014044ba 875 //
876 // --- Translate the muon parameters in the dimuon rest frame
877 //
6579f87c 878 zaxis=(pDimuCM.Vect()).Unit();
014044ba 879
6579f87c 880 beta=(-1./pDimuCM.E())*pDimuCM.Vect();
014044ba 881
6579f87c 882 pProjCM.SetPxPyPzE(0.,0.,-Energy,TMath::Sqrt(Energy*Energy+mp*mp));
883 pTargCM.SetPxPyPzE(0.,0.,Energy,TMath::Sqrt(Energy*Energy+mp*mp));
014044ba 884
6579f87c 885 pProjDimu=pProjCM;
886 pTargDimu=pTargCM;
014044ba 887
6579f87c 888 pProjDimu.Boost(beta);
889 pTargDimu.Boost(beta);
014044ba 890
6579f87c 891 yaxis=((pProjDimu.Vect()).Cross(pTargDimu.Vect())).Unit();
014044ba 892 xaxis=(yaxis.Cross(zaxis)).Unit();
893
6579f87c 894 pMu1Dimu=pMu1CM;
895 pMu2Dimu=pMu2CM;
896 pMu1Dimu.Boost(beta);
897 pMu2Dimu.Boost(beta);
014044ba 898
899 Double_t phi;
6579f87c 900 if(charge1>0) phi = TMath::ATan2((pMu1Dimu.Vect()).Dot(yaxis),(pMu1Dimu.Vect()).Dot(xaxis));
901 else phi = TMath::ATan2((pMu2Dimu.Vect()).Dot(yaxis),(pMu2Dimu.Vect()).Dot(xaxis));
014044ba 902
903 return phi;
904}
905
906//________________________________________________________________________
907void AliCFMuonResUpsilon::Terminate(Option_t *)
908{
909 // draw result of the Invariant mass MC and ESD
910
e199f3f0 911/*
014044ba 912 TH1D *h1 = dynamic_cast<TH1D*>(GetOutputData(1));
913 AliCFContainer *cont = dynamic_cast<AliCFContainer*> (GetOutputData(2));
914
915 TH1D *kmass = cont->ShowProjection(3,0);
916 TH1D *rmass = cont->ShowProjection(3,1);
917 TH1D *mmass = cont->ShowProjection(3,4);
918
919 TCanvas *c1 = new TCanvas("AliCFMuonResUpsilon","UPSILON Container",0,0,800,800);
920 c1->Divide(2,2);
921 c1->cd(1);
922 kmass->Draw("HIST");
923 c1->cd(2);
924 rmass->Draw("HIST");
925 c1->cd(3);
926 mmass->Draw("HIST");
927 c1->cd(4);
928 h1->Draw("HIST");
e199f3f0 929 */
014044ba 930}
931//________________________________________________________________________
932
933#endif