]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGDQ/dielectron/AliAnalysisTaskMultiDielectronTG.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fixing phys. selection in event mixing framework
[u/mrichter/AliRoot.git] / PWGDQ / dielectron / AliAnalysisTaskMultiDielectronTG.cxx
CommitLineData
b1f673b5 1/*************************************************************************
2* Copyright(c) 1998-2009, 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// //
18// Basic Analysis Task //
19// //
20///////////////////////////////////////////////////////////////////////////
21
22#include <TChain.h>
23#include <TH1D.h>
24
25#include <AliCFContainer.h>
26#include <AliInputEventHandler.h>
27#include <AliESDInputHandler.h>
28#include <AliAODInputHandler.h>
29#include <AliAnalysisManager.h>
30#include <AliVEvent.h>
31#include <AliTriggerAnalysis.h>
32#include <AliPIDResponse.h>
33#include <AliTPCPIDResponse.h>
34
35#include "AliDielectron.h"
36#include "AliDielectronHistos.h"
37#include "AliDielectronCF.h"
38#include "AliDielectronMC.h"
39#include "AliDielectronMixingHandler.h"
40#include "AliAnalysisTaskMultiDielectronTG.h"
41
42#include "AliESDtrack.h"
43#include "AliESDtrackCuts.h"
44
45#include "AliLog.h"
46
47#include <vector>
48#include <deque>
49#include <cstdlib>
50#include <TRandom.h>
51
7329335b 52using std::vector;
53using std::deque;
54
55const char* kPairClassNamesTG[7] = {
b1f673b5 56 "unlike",
57 "like_pp",
58 "like_ee",
59 "mixunlike_pe",
60 "mixunlike_ep",
61 "mixlike_pp",
62 "mixlike_ee"
63};
64
65
66ClassImp(AliDielectronSingleTG)
67ClassImp(AliAnalysisTaskMultiDielectronTG)
68
69//_________________________________________________________________________________
70AliAnalysisTaskMultiDielectronTG::AliAnalysisTaskMultiDielectronTG() :
71 AliAnalysisTaskSE(),
72 fListDielectron(),
73 fListHistos(),
74 fListCF(),
7329335b 75 fQAElectron(),
b1f673b5 76 fSelectPhysics(kFALSE),
77 fTriggerMask(AliVEvent::kMB),
78 fExcludeTriggerMask(0),
79 fTriggerOnV0AND(kFALSE),
80 fRejectPileup(kFALSE),
81 fTriggerLogic(kAny),
82 fTriggerAnalysis(0x0),
83 fEventFilter(0x0),
84 fCutsMother(0x0),
85 fEventStat(0x0),
86 fEvent(0x0),
87 fdEdXvsPt(0x0),
88 fdEdXnSigmaElecvsPt(0x0),
89 fdEdXvsPtTOF(0x0),
90 fdEdXnSigmaElecvsPtTOF(0x0),
91 fTOFbetavsPt(0x0),
92 fTOFnSigmaElecvsPt(0x0),
7329335b 93 fNCrossedRowsTPC(0x0),
94 fChi2ClusTPC(0x0),
95 fRatioCrossClusTPC(0x0),
96 fVem(0x0),
97 fVep(0x0),
98 fVemtmp(0x0),
99 fVeptmp(0x0),
100 fdconvphiv(acos(-1.0)),
805bd069 101 fdop(0),
7329335b 102 fbz(0),
620b34a5 103 fdv0mixing(kTRUE)
b1f673b5 104{
105 //
106 // Constructor
107 //
108}
109
110//_________________________________________________________________________________
111AliAnalysisTaskMultiDielectronTG::AliAnalysisTaskMultiDielectronTG(const char *name) :
112 AliAnalysisTaskSE(name),
113 fListDielectron(),
114 fListHistos(),
115 fListCF(),
7329335b 116 fQAElectron(),
b1f673b5 117 fSelectPhysics(kFALSE),
118 fTriggerMask(AliVEvent::kMB),
119 fExcludeTriggerMask(0),
120 fTriggerOnV0AND(kFALSE),
121 fRejectPileup(kFALSE),
122 fTriggerLogic(kAny),
123 fTriggerAnalysis(0x0),
124 fEventFilter(0x0),
125 fCutsMother(0x0),
126 fEventStat(0x0),
127 fEvent(0x0),
128 fdEdXvsPt(0x0),
129 fdEdXnSigmaElecvsPt(0x0),
130 fdEdXvsPtTOF(0x0),
131 fdEdXnSigmaElecvsPtTOF(0x0),
132 fTOFbetavsPt(0x0),
133 fTOFnSigmaElecvsPt(0x0),
7329335b 134 fNCrossedRowsTPC(0x0),
135 fChi2ClusTPC(0x0),
136 fRatioCrossClusTPC(0x0),
137 fVem(0x0),
138 fVep(0x0),
139 fVemtmp(0x0),
140 fVeptmp(0x0),
141 fdconvphiv(acos(-1.0)),
805bd069 142 fdop(0),
7329335b 143 fbz(0),
620b34a5 144 fdv0mixing(kTRUE)
b1f673b5 145{
146 //
147 // Constructor
148 //
149 DefineInput(0,TChain::Class());
150 DefineOutput(1, TList::Class());
151 DefineOutput(2, TList::Class());
152 DefineOutput(3, TH1D::Class());
153 fListHistos.SetName("Dielectron_Histos_Multi");
154 fListCF.SetName("Dielectron_CF_Multi");
155 fListDielectron.SetOwner();
156 fListHistos.SetOwner();
157 fListCF.SetOwner();
158
159 ///////////////
7329335b 160 for(int i=0;i<fgkNDIE; i++){
161 for(int j=0;j<fgkNZBIN;j++){
162 for(int k=0;k<fgkNCENT;k++){
163 for(int l=0; l<fgkNRPBIN; l++){
164 fibuf[i][j][k][l] = 0;
165 fpoolp[i][j][k][l].clear();
166 fpoolm[i][j][k][l].clear();
167 for(int ib=0;ib<fgkNBUF; ib++){
168 fvep[ib][i][j][k][l].clear();
169 fvem[ib][i][j][k][l].clear();
b1f673b5 170 }
171 }
172 }
173 }
174 }
175
176
177
178}
179
180//_________________________________________________________________________________
181AliAnalysisTaskMultiDielectronTG::~AliAnalysisTaskMultiDielectronTG()
182{
183 //
184 // Destructor
185 //
186
187 //histograms and CF are owned by the dielectron framework.
188 //however they are streamed to file, so in the first place the
189 //lists need to be owner...
190 fListHistos.SetOwner(kFALSE);
191 fListCF.SetOwner(kFALSE);
192
7329335b 193 for(int i=0;i<fgkNDIE; i++){
194 for(int j=0;j<fgkNZBIN;j++){
195 for(int k=0;k<fgkNCENT;k++){
196 for(int l=0; l<fgkNRPBIN; l++){
197 fibuf[i][j][k][l] = 0;
198 fpoolp[i][j][k][l].clear();
199 fpoolm[i][j][k][l].clear();
200 for(int ib=0;ib<fgkNBUF; ib++){
201 fvep[ib][i][j][k][l].clear();
202 fvem[ib][i][j][k][l].clear();
b1f673b5 203 }
204 }
205 }
206 }
207 }
208}
209//_________________________________________________________________________________
210void AliAnalysisTaskMultiDielectronTG::UserCreateOutputObjects()
211{
212 //
213 // Add all histogram manager histogram lists to the output TList
214 //
215
216 if (!fListHistos.IsEmpty()||!fListCF.IsEmpty()) return; //already initialised
217
218// AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
219// Bool_t isESD=man->GetInputEventHandler()->IsA()==AliESDInputHandler::Class();
220// Bool_t isAOD=man->GetInputEventHandler()->IsA()==AliAODInputHandler::Class();
221
222 TIter nextDie(&fListDielectron);
223 AliDielectron *die=0;
224 while ( (die=static_cast<AliDielectron*>(nextDie())) ){
225 die->Init();
226 if (die->GetHistogramList()) fListHistos.Add(const_cast<THashList*>(die->GetHistogramList()));
227 if (die->GetCFManagerPair()) fListCF.Add(const_cast<AliCFContainer*>(die->GetCFManagerPair()->GetContainer()));
228 }
229
230 Int_t cuts=fListDielectron.GetEntries();
231 Int_t nbins=kNbinsEvent+2*cuts;
232 if (!fEventStat){
233 fEventStat=new TH1D("hEventStat","Event statistics",nbins,0,nbins);
234 fEventStat->GetXaxis()->SetBinLabel(1,"Before Phys. Sel.");
235 fEventStat->GetXaxis()->SetBinLabel(2,"After Phys. Sel.");
236
237 //default names
238 fEventStat->GetXaxis()->SetBinLabel(3,"Bin3 not used");
239 fEventStat->GetXaxis()->SetBinLabel(4,"Bin4 not used");
240 fEventStat->GetXaxis()->SetBinLabel(5,"Bin5 not used");
241
242 if(fTriggerOnV0AND) fEventStat->GetXaxis()->SetBinLabel(3,"V0and triggers");
243 if (fEventFilter) fEventStat->GetXaxis()->SetBinLabel(4,"After Event Filter");
244 if (fRejectPileup) fEventStat->GetXaxis()->SetBinLabel(5,"After Pileup rejection");
245
246 for (Int_t i=0; i<cuts; ++i){
247 fEventStat->GetXaxis()->SetBinLabel((kNbinsEvent+1)+2*i,Form("#splitline{1 candidate}{%s}",fListDielectron.At(i)->GetName()));
248 fEventStat->GetXaxis()->SetBinLabel((kNbinsEvent+2)+2*i,Form("#splitline{With >1 candidate}{%s}",fListDielectron.At(i)->GetName()));
249 }
250 }
251
252 if (!fTriggerAnalysis) fTriggerAnalysis=new AliTriggerAnalysis;
253 fTriggerAnalysis->EnableHistograms();
254 fTriggerAnalysis->SetAnalyzeMC(AliDielectronMC::Instance()->HasMC());
255
256
257
258 int nbinx=400;
7329335b 259 float maxx=20;
260 float minx=0.2;
261 float binw = (TMath::Log(maxx)-TMath::Log(minx))/nbinx;
b1f673b5 262 double xbin[401];
263 for(int ii=0;ii<nbinx+1;ii++){
7329335b 264 xbin[ii] = TMath::Exp(TMath::Log(minx) + 0.5*binw+binw*ii);
b1f673b5 265 }
266
267
7329335b 268 fQAElectron = new TList();
269 fQAElectron->SetName("QAElectron");
270 fQAElectron->SetOwner();
b1f673b5 271
272
273 fEvent = new TH1D("Event","centrality", 100,0,100);
7329335b 274 fQAElectron->Add(fEvent);
b1f673b5 275 fdEdXvsPt = new TH2D("dEdXvsPt","dE/dX vs. PT of TPC", nbinx, xbin, 2000,0,200);
7329335b 276 fQAElectron->Add(fdEdXvsPt);
b1f673b5 277 fdEdXnSigmaElecvsPt = new TH2D("fdEdXnSigmaElecvsPt"," dE/dX normalized to electron vs. pT of TPC",
278 nbinx, xbin, 2000, -10, 10);
7329335b 279 fQAElectron->Add(fdEdXnSigmaElecvsPt);
b1f673b5 280
281 fdEdXvsPtTOF = new TH2D("dEdXvsPtTOF","dE/dX vs. PT of TPC", nbinx, xbin, 2000,0,200);
7329335b 282 fQAElectron->Add(fdEdXvsPtTOF);
b1f673b5 283 fdEdXnSigmaElecvsPtTOF = new TH2D("fdEdXnSigmaElecvsPtTOF"," dE/dX normalized to electron vs. pT of TPC",
284 nbinx, xbin, 2000, -10, 10);
7329335b 285 fQAElectron->Add(fdEdXnSigmaElecvsPtTOF);
b1f673b5 286
287
288
289 fTOFbetavsPt = new TH2D("fTOFbetavsPt","TOF beta vs. p", 400, 0, 20, 1200, 0, 1.2);
7329335b 290 fQAElectron->Add(fTOFbetavsPt);
b1f673b5 291 fTOFnSigmaElecvsPt = new TH2D("fTOFnSigmaElecvsPt","TOF nsigma for electron", 400, 0, 20, 2000, -10, 10);
7329335b 292 fQAElectron->Add(fTOFnSigmaElecvsPt);
b1f673b5 293
7329335b 294 fNCrossedRowsTPC = new TH2F("fNCrossedRowsTPC", "TPC nCrossed Rows vs. pT", 200, 0, 20, 200, 0, 200);
295 fQAElectron->Add(fNCrossedRowsTPC);
296 fChi2ClusTPC = new TH2F("fChi2ClusTPC", "hChi2ClusTPC vs. pT", 200, 0, 20, 200, 0, 10);
297 fQAElectron->Add(fChi2ClusTPC);
b1f673b5 298
7329335b 299 fRatioCrossClusTPC = new TH2F("fRatioCrossClusTPC", "hRatioCrossClusTPC vs. pT", 200, 0, 20, 200, 0, 10);
300 fQAElectron->Add(fRatioCrossClusTPC);
b1f673b5 301
7329335b 302 fListHistos.Add(fQAElectron);
b1f673b5 303
304 fListHistos.SetOwner();
305
306 PostData(1, &fListHistos);
307 PostData(2, &fListCF);
308 PostData(3, fEventStat);
309
310 fCutsMother = new AliESDtrackCuts;
311 fCutsMother->SetDCAToVertex2D(kTRUE);
312 fCutsMother->SetMaxDCAToVertexZ(3.0);
313 fCutsMother->SetMaxDCAToVertexXY(1.0);
314 fCutsMother->SetPtRange( 0.05 , 200.0);
315 fCutsMother->SetEtaRange( -0.84 , 0.84 );
316 fCutsMother->SetAcceptKinkDaughters(kFALSE);
317 fCutsMother->SetRequireITSRefit(kTRUE);
318 fCutsMother->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kFirst);
319 fCutsMother->SetMinNClustersITS(3);
320 fCutsMother->SetRequireTPCRefit(kTRUE);
321
322
805bd069 323 AliInfo(Form("PairCutType = %d %d %d %d %d",
324 fRejectPairFlag[0],
325 fRejectPairFlag[1],
326 fRejectPairFlag[2],
327 fRejectPairFlag[3],
328 fRejectPairFlag[4]));
329
b1f673b5 330}
331
332//_________________________________________________________________________________
333void AliAnalysisTaskMultiDielectronTG::UserExec(Option_t *)
334{
335 //
336 // Main loop. Called for every event
337 //
338
339 if (fListHistos.IsEmpty()&&fListCF.IsEmpty()) return;
340
341 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
342 Bool_t isESD=man->GetInputEventHandler()->IsA()==AliESDInputHandler::Class();
343 Bool_t isAOD=man->GetInputEventHandler()->IsA()==AliAODInputHandler::Class();
344
345 AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
346 if (!inputHandler) return;
347
348// AliPIDResponse *pidRes=inputHandler->GetPIDResponse();
349 if ( inputHandler->GetPIDResponse() ){
350 // for the 2.76 pass2 MC private train. Together with a sigma shift of -0.169
351// pidRes->GetTPCResponse().SetSigma(4.637e-3,2.41332105409873257e+04);
352 AliDielectronVarManager::SetPIDResponse( inputHandler->GetPIDResponse() );
353 } else {
354 AliFatal("This task needs the PID response attached to the input event handler!");
355 }
356
357 // Was event selected ?
358 ULong64_t isSelected = AliVEvent::kAny;
359 Bool_t isRejected = kFALSE;
360 if( fSelectPhysics && inputHandler){
361 if((isESD && inputHandler->GetEventSelection()) || isAOD){
362 isSelected = inputHandler->IsEventSelected();
363 if (fExcludeTriggerMask && (isSelected&fExcludeTriggerMask)) isRejected=kTRUE;
364 if (fTriggerLogic==kAny) isSelected&=fTriggerMask;
365 else if (fTriggerLogic==kExact) isSelected=((isSelected&fTriggerMask)==fTriggerMask);
366 }
367 }
368
369
370 //Before physics selection
371 fEventStat->Fill(kAllEvents);
372 if (isSelected==0||isRejected) {
373 PostData(3,fEventStat);
374 return;
375 }
376 //after physics selection
377 fEventStat->Fill(kSelectedEvents);
378
379 //V0and
380 if(fTriggerOnV0AND){
805bd069 381 if(isESD){if (!fTriggerAnalysis->IsOfflineTriggerFired(static_cast<AliESDEvent*>(InputEvent()), AliTriggerAnalysis::kV0AND))
b1f673b5 382 return;}
805bd069 383 if(isAOD){if(!((static_cast<AliAODEvent*>(InputEvent()))->GetVZEROData()->GetV0ADecision() == AliVVZERO::kV0BB &&
384 (static_cast<AliAODEvent*>(InputEvent()))->GetVZEROData()->GetV0CDecision() == AliVVZERO::kV0BB) )
385 return;}
386 }
b1f673b5 387
388
389 fEventStat->Fill(kV0andEvents);
390
391 //Fill Event histograms before the event filter
392 TIter nextDie(&fListDielectron);
393 AliDielectron *die=0;
394 Double_t values[AliDielectronVarManager::kNMaxValues]={0};
395 Double_t valuesMC[AliDielectronVarManager::kNMaxValues]={0};
396 AliDielectronVarManager::SetEvent(InputEvent());
397 AliDielectronVarManager::Fill(InputEvent(),values);
398 AliDielectronVarManager::Fill(InputEvent(),valuesMC);
399 Bool_t hasMC=AliDielectronMC::Instance()->HasMC();
400 if (hasMC) {
401 if (AliDielectronMC::Instance()->ConnectMCEvent())
402 AliDielectronVarManager::Fill(AliDielectronMC::Instance()->GetMCEvent(),valuesMC);
403 }
404
405 while ( (die=static_cast<AliDielectron*>(nextDie())) ){
406 AliDielectronHistos *h=die->GetHistoManager();
407 if (h){
408 if (h->GetHistogramList()->FindObject("Event_noCuts"))
409 h->FillClass("Event_noCuts",AliDielectronVarManager::kNMaxValues,values);
410 if (hasMC && h->GetHistogramList()->FindObject("MCEvent_noCuts"))
411 h->FillClass("Event_noCuts",AliDielectronVarManager::kNMaxValues,valuesMC);
412 }
413 }
414 nextDie.Reset();
415
416 //event filter
417 if (fEventFilter) {
418 if (!fEventFilter->IsSelected(InputEvent())) return;
419 }
420 fEventStat->Fill(kFilteredEvents);
421
422 //pileup
423 if (fRejectPileup){
424 if (InputEvent()->IsPileupFromSPD(3,0.8,3.,2.,5.)) return;
425 }
426 fEventStat->Fill(kPileupEvents);
427
428 //bz for AliKF
7329335b 429 fbz = InputEvent()->GetMagneticField();
430 AliKFParticle::SetField( fbz );
b1f673b5 431 AliDielectronPID::SetCorrVal((Double_t)InputEvent()->GetRunNumber());
432
433 //Process event in all AliDielectron instances
434 // TIter nextDie(&fListDielectron);
435 // AliDielectron *die=0;
436 Int_t idie=0;
437 while ( (die=static_cast<AliDielectron*>(nextDie())) ){
805bd069 438 //AliInfo(Form(" **************** die->Process(InputEvent()) : %d **************************", idie));
b1f673b5 439 die->SetDontClearArrays(kTRUE);
440 die->Process(InputEvent());
441 if (die->HasCandidates()){
442 Int_t ncandidates=die->GetPairArray(1)->GetEntriesFast();
443 if (ncandidates==1) fEventStat->Fill((kNbinsEvent)+2*idie);
444 else if (ncandidates>1) fEventStat->Fill((kNbinsEvent+1)+2*idie);
445 }
446
447
448 AliDielectronVarManager::Fill(InputEvent(), fgValues);
449 for (Int_t ii=0; ii<2; ++ii){
450 TObjArray *obj = (TObjArray*)die->GetTrackArray(ii);
451 Int_t ntracks=obj->GetEntriesFast();
452 //AliInfo(Form(" ************** # of tracks = %d", ntracks));
453 for (Int_t itrack=0; itrack<ntracks; ++itrack){
454
455 ////////////////////////////////////////////////////////////////////
456 AliDielectronVarManager::Fill(obj->UncheckedAt(itrack), fgValues);
457 ////////////////////////////////////////////////////////////////////
458
459 if(fgValues[AliDielectronVarManager::kCharge]>0){
7329335b 460 fVeptmp.push_back(new AliDielectronSingleTG(1,
b1f673b5 461 fgValues[AliDielectronVarManager::kCentrality],
462 fgValues[AliDielectronVarManager::kXv],
463 fgValues[AliDielectronVarManager::kYv],
464 fgValues[AliDielectronVarManager::kZv],
465 fgValues[AliDielectronVarManager::kPx],
466 fgValues[AliDielectronVarManager::kPy],
467 fgValues[AliDielectronVarManager::kPz],
468 fgValues[AliDielectronVarManager::kPt],
469 fgValues[AliDielectronVarManager::kEta],
470 fgValues[AliDielectronVarManager::kPhi],
471 fgValues[AliDielectronVarManager::kTheta],
472 1, 1, static_cast<AliVTrack*>(obj->UncheckedAt(itrack)))
473 );
474 }else if(fgValues[AliDielectronVarManager::kCharge]<0){
7329335b 475 fVemtmp.push_back(new AliDielectronSingleTG(-1,
b1f673b5 476 fgValues[AliDielectronVarManager::kCentrality],
477 fgValues[AliDielectronVarManager::kXv],
478 fgValues[AliDielectronVarManager::kYv],
479 fgValues[AliDielectronVarManager::kZv],
480 fgValues[AliDielectronVarManager::kPx],
481 fgValues[AliDielectronVarManager::kPy],
482 fgValues[AliDielectronVarManager::kPz],
483 fgValues[AliDielectronVarManager::kPt],
484 fgValues[AliDielectronVarManager::kEta],
485 fgValues[AliDielectronVarManager::kPhi],
486 fgValues[AliDielectronVarManager::kTheta],
487 1, 1, static_cast<AliVTrack*>(obj->UncheckedAt(itrack)))
488 );
489 }
490 }
491 }
7329335b 492 //AliInfo(Form("size of e and p = %d %d", (int)fVeptmp.size(), (int)fVemtmp.size()));
493
b1f673b5 494
7329335b 495 CheckGhostPairs(fVeptmp);
496 CheckGhostPairs(fVemtmp);
805bd069 497 if(fRejectPairFlag[idie]==1 || fRejectPairFlag[idie]==2){
498 RejectPairs(fVeptmp, fVemtmp, idie);
499 }
7329335b 500 RandomizePool(fVeptmp, fVemtmp);
501 CalcPair(fVep, fVem, die, idie);
b1f673b5 502
7329335b 503 // AliInfo(Form("size of e and p (after) = %d %d", (int)fVep.size(), (int)fVem.size()));
b1f673b5 504
7329335b 505 double dwcent = 100.0/fgkNCENT;
506 double dwiz = 20.0/fgkNZBIN;
507 double dwrp = acos(-1.0)/fgkNRPBIN;
b1f673b5 508
7329335b 509 int icent = (int)(fgValues[AliDielectronVarManager::kCentrality]/dwcent);
510 int izbin = (int)((fgValues[AliDielectronVarManager::kZvPrim]+10)/dwiz);
511 int irp = (int)((fgValues[AliDielectronVarManager::kV0ACrpH2])/dwrp);
b1f673b5 512
513 if(icent<0) icent=0;
7329335b 514 if(icent>=fgkNCENT) icent=fgkNCENT-1;
b1f673b5 515 if(izbin<0) izbin=0;
7329335b 516 if(izbin>=fgkNZBIN) izbin=fgkNZBIN-1;
b1f673b5 517 if(irp<0) irp=0;
7329335b 518 if(irp>=fgkNRPBIN) irp=fgkNRPBIN-1;
b1f673b5 519
7329335b 520 fvep[fibuf[idie][izbin][icent][irp]][idie][izbin][icent][irp].clear();
521 for(int iep = 0; iep<(int)fVep.size();iep++) {
522 fvep[fibuf[idie][izbin][icent][irp]][idie][izbin][icent][irp].push_back(fVep[iep]);
523 fpoolp[idie][izbin][icent][irp].push_back(fVep[iep]);
524 if(fpoolp[idie][izbin][icent][irp].size()>fgkMAXPOOL) {
525 fpoolp[idie][izbin][icent][irp].pop_front();
b1f673b5 526 }
527 }
7329335b 528 fvem[fibuf[idie][izbin][icent][irp]][idie][izbin][icent][irp].clear();
529 for(int iem = 0; iem<(int)fVem.size();iem++) {
530 fvem[fibuf[idie][izbin][icent][irp]][idie][izbin][icent][irp].push_back(fVem[iem]);
531 fpoolm[idie][izbin][icent][irp].push_back(fVem[iem]);
532 if(fpoolm[idie][izbin][icent][irp].size()>fgkMAXPOOL) {
533 fpoolm[idie][izbin][icent][irp].pop_front();
b1f673b5 534 }
535 }
536
537
7329335b 538 fibuf[idie][izbin][icent][irp]++;
539 if(fibuf[idie][izbin][icent][irp]>= fgkNBUF) fibuf[idie][izbin][icent][irp]=0;
b1f673b5 540
541
7329335b 542 fVeptmp.clear();
543 fVemtmp.clear();
544 fVep.clear();
545 fVem.clear();
b1f673b5 546
547
548 ++idie;
549 }
550
551
552 AliESDEvent *fESD = dynamic_cast<AliESDEvent*>(InputEvent());
553 fEvent->Fill(values[AliDielectronVarManager::kCentrality]);
554 for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
555 AliESDtrack* track = fESD->GetTrack(iTracks);
556 if (!track) {
557 Printf("ERROR: Could not receive track %d", iTracks);
558 continue;
559 }
560 if(!fCutsMother->AcceptTrack(track)) continue;
561 fdEdXvsPt->Fill(track->GetTPCmomentum(), track->GetTPCsignal());
562 fdEdXnSigmaElecvsPt->Fill(track->GetTPCmomentum(),
563 AliDielectronVarManager::GetPIDResponse()->NumberOfSigmasTPC(track,
564 AliPID::kElectron)
565 -AliDielectronPID::GetCorrVal());
566 /// for beta caliculaton
567 Double_t l = track->GetIntegratedLength(); // cm
568 Double_t t = track->GetTOFsignal();
569 Double_t t0 = AliDielectronVarManager::GetPIDResponse()->GetTOFResponse().GetTimeZero(); // ps
570 Double_t beta = 0;
571 if( (l < 360. || l > 800.) || (t <= 0.) || (t0 >999990.0) ) {
572 beta=-9999;
573 }
574 else {
575 t -= t0; // subtract the T0
576 l *= 0.01; // cm ->m
577 t *= 1e-12; //ps -> s
578
579 Double_t v = l / t;
580 beta = v / TMath::C();
581 }
582
583 fTOFbetavsPt->Fill(track->GetTPCmomentum(), beta);
584 fTOFnSigmaElecvsPt->Fill(track->GetTPCmomentum(),
585 AliDielectronVarManager::GetPIDResponse()->NumberOfSigmasTOF(track,
586 AliPID::kElectron));
587 ////rough PID is required
588 if( fabs(AliDielectronVarManager::GetPIDResponse()->NumberOfSigmasTOF(track, AliPID::kElectron))<3){
589 fdEdXvsPtTOF->Fill(track->GetTPCmomentum(), track->GetTPCsignal());
590 fdEdXnSigmaElecvsPtTOF->Fill(track->GetTPCmomentum(),
591 AliDielectronVarManager::GetPIDResponse()->NumberOfSigmasTPC(track,
592 AliPID::kElectron)
593 -AliDielectronPID::GetCorrVal());
594
595
596 if(track->GetTPCsignal()>70 && track->GetTPCsignal()<90){
7329335b 597 fNCrossedRowsTPC->Fill(track->GetTPCmomentum(), track->GetTPCCrossedRows());
598 //fChi2ClusTPC->Fill(track->GetTPCmomentum(), track->GetTPCchi2()/track->GetTPCNcls());
599 fChi2ClusTPC->Fill(track->GetTPCmomentum(), track->GetTPCchi2());
600 fRatioCrossClusTPC->Fill(track->GetTPCmomentum(), track->GetTPCCrossedRows()/track->GetTPCNclsF());
b1f673b5 601 }
602 }
603 }
604
605 PostData(1, &fListHistos);
606 PostData(2, &fListCF);
607 PostData(3,fEventStat);
608}
609
610//_________________________________________________________________________________
611void AliAnalysisTaskMultiDielectronTG::FinishTaskOutput()
612{
613 //
614 // Write debug tree
615 //
616 TIter nextDie(&fListDielectron);
617 AliDielectron *die=0;
618 while ( (die=static_cast<AliDielectron*>(nextDie())) ){
619 die->SaveDebugTree();
620 AliDielectronMixingHandler *mix=die->GetMixingHandler();
621// printf("\n\n\n===============\ncall mix in Terminate: %p (%p)\n=================\n\n",mix,die);
622 if (mix) mix->MixRemaining(die);
623 }
624 PostData(1, &fListHistos);
625 PostData(2, &fListCF);
626}
627
628//_________________________________________________________________________________
7329335b 629void AliAnalysisTaskMultiDielectronTG::CheckGhostPairs(vector<AliDielectronSingleTG*> e1)
630{
631 ////// Check whether the pairs are in ghost pairs
632 ///// ghost means that one of the pairs is highly associated but reconstructed as a true particle
633
634
b1f673b5 635 bool reject = false;
636 if(e1.size()>1){
637 for(int i1=0; i1<(int)e1.size(); i1++){
638 reject = false;
639 for(int i2=i1+1; i2<(int)e1.size(); i2++){
620b34a5 640 if( fabs(e1[i1]->Phi() - e1[i2]->Phi())<0.01 ){
b1f673b5 641 reject = true;
642 e1[i2]->SetGstFlag(0);
643 }
644 }
645 if(reject==true)e1[i1]->SetGstFlag(0);
646 }
647 }
648}
649
650//_________________________________________________________________________________
7329335b 651void AliAnalysisTaskMultiDielectronTG::RandomizePool(vector<AliDielectronSingleTG*> e1, vector<AliDielectronSingleTG*> e2)
652{
b1f673b5 653
7329335b 654 ///// Randomize the pool constent to cancel the filling scheme of the single tracks
655 ////
656
657
b1f673b5 658 int size1 = e1.size();
7329335b 659 int usedindex[1000];
b1f673b5 660 for(int i=0;i<1000;i++){
7329335b 661 usedindex[i] = -1;
b1f673b5 662 }
663 for(int i=0;i<size1;i++){
7329335b 664 usedindex[i] = 0;
b1f673b5 665 }
666
667 for(int i=0;i<size1;i++){
668 int j = (int)(gRandom->Uniform(0,size1));
7329335b 669 while(usedindex[j]==1){
b1f673b5 670 j = (int)(gRandom->Uniform(0,size1));
671 }
672 if( (e1[j]->GetGstFlag()==1) &&
673 (e1[j]->GetConvFlag()==1)
674 ){
7329335b 675 fVep.push_back(e1[j]);
b1f673b5 676 }
7329335b 677 usedindex[j] = 1;
b1f673b5 678 }
679
680
681 int size2 = e2.size();
682 for(int i=0;i<1000;i++){
7329335b 683 usedindex[i] = -1;
b1f673b5 684 }
685 for(int i=0;i<size2;i++){
7329335b 686 usedindex[i] = 0;
b1f673b5 687 }
688
689 for(int i=0;i<size2;i++){
690 int j = (int)(gRandom->Uniform(0,size2));
7329335b 691 while(usedindex[j]==1){
b1f673b5 692 j = (int)(gRandom->Uniform(0,size2));
693 }
694 if( (e2[j]->GetGstFlag()==1) &&
695 (e2[j]->GetConvFlag()==1)
696 ){
7329335b 697 fVem.push_back(e2[j]);
b1f673b5 698 }
7329335b 699 usedindex[j] = 1;
b1f673b5 700 }
701}
702
703
704//_________________________________________________________________________________
7329335b 705void AliAnalysisTaskMultiDielectronTG::CalcPair(vector<AliDielectronSingleTG*> ve1,
706 vector<AliDielectronSingleTG*> ve2, AliDielectron *die, Int_t idie)
707{
708
709 //
710 // main routine for the pair procesing
711 //
712
b1f673b5 713
714 for(int i1=0; i1<(int)ve1.size(); i1++){
715 for(int i2=0; i2<(int)ve2.size(); i2++){
805bd069 716 FillPair(ve1[i1], ve2[i2], 0, die, idie);
b1f673b5 717 }
718 }
719
720 for(int i1=0; i1<(int)ve1.size(); i1++){
721 for(int i2=i1+1; i2<(int)ve1.size(); i2++){
805bd069 722 FillPair(ve1[i1], ve1[i2], 1, die, idie);
b1f673b5 723 }
724 }
725
726 for(int i1=0; i1<(int)ve2.size(); i1++){
727 for(int i2=i1+1; i2<(int)ve2.size(); i2++){
805bd069 728 FillPair(ve2[i1], ve2[i2], 2, die, idie );
b1f673b5 729 }
730 }
731
732
7329335b 733 double dwcent = 100.0/fgkNCENT;
734 double dwiz = 20.0/fgkNZBIN;
735 double dwrp = acos(-1.0)/fgkNRPBIN;
b1f673b5 736
7329335b 737 int icent = (int)(fgValues[AliDielectronVarManager::kCentrality]/dwcent);
738 int izbin = (int)((fgValues[AliDielectronVarManager::kZvPrim]+10)/dwiz);
739 int irp = (int)((fgValues[AliDielectronVarManager::kV0ACrpH2])/dwrp);
b1f673b5 740
741 if(icent<0) icent=0;
7329335b 742 if(icent>=fgkNCENT) icent=fgkNCENT-1;
b1f673b5 743 if(izbin<0) izbin=0;
7329335b 744 if(izbin>=fgkNZBIN) izbin=fgkNZBIN-1;
b1f673b5 745 if(irp<0) irp=0;
7329335b 746 if(irp>=fgkNRPBIN) irp=fgkNRPBIN-1;
b1f673b5 747
748 int nmixed;
749 if(ve1.size()>0) {
750 //
751 // Now mixed event for +- pairs
752 //
753 nmixed = 0;
7329335b 754 for(int ibuf=0;(nmixed<fgkNMix);ibuf++) {
b1f673b5 755 int ntry = 0;
620b34a5 756 while(ntry<fgkMAXTRY) {
7329335b 757 ReshuffleBuffer(fvem[ibuf][idie][izbin][icent][irp],fpoolm[idie][izbin][icent][irp]);
b1f673b5 758 ntry++;
759 }
760 for(int i1=0; i1<(int)ve1.size(); i1++){
7329335b 761 for(int i2=0; i2<(int)fvem[ibuf][idie][izbin][icent][irp].size(); i2++){
805bd069 762 FillPair(ve1[i1],fvem[ibuf][idie][izbin][icent][irp][i2], 3, die, idie);
b1f673b5 763 }
764 }
765 ++nmixed;
766 }//for(ibuf)
767 }
768 if(ve2.size()>0) {
769 //
805bd069 770 // Now mixed event for -+ pairs
b1f673b5 771 //
772 nmixed = 0;
7329335b 773 for(int ibuf=0;(nmixed<fgkNMix);ibuf++) {
b1f673b5 774 int ntry = 0;
620b34a5 775 while(ntry<fgkMAXTRY) {
7329335b 776 ReshuffleBuffer(fvep[ibuf][idie][izbin][icent][irp],fpoolp[idie][izbin][icent][irp]);
b1f673b5 777 ntry++;
778 }
779 for(int i1=0; i1<(int)ve2.size(); i1++){
7329335b 780 for(int i2=0; i2<(int)fvep[ibuf][idie][izbin][icent][irp].size(); i2++){
805bd069 781 FillPair(fvep[ibuf][idie][izbin][icent][irp][i2],ve2[i1],4, die, idie);
b1f673b5 782 }
783 }
784 ++nmixed;
785 }//for(ibuf)
786 }
787
788 if(ve1.size()>0) {
789 //
790 // Now mixed event for ++ pairs
791 //
792 nmixed = 0;
7329335b 793 for(int ibuf=0;(nmixed<fgkNMix);ibuf++) {
b1f673b5 794 int ntry = 0;
620b34a5 795 while(ntry<fgkMAXTRY) {
7329335b 796 ReshuffleBuffer(fvep[ibuf][idie][izbin][icent][irp],fpoolp[idie][izbin][icent][irp]);
b1f673b5 797 ntry++;
798 }
799 for(int i1=0; i1<(int)ve1.size(); i1++){
7329335b 800 for(int i2=0;i2<(int)fvep[ibuf][idie][izbin][icent][irp].size(); i2++){
805bd069 801 FillPair(ve1[i1],fvep[ibuf][idie][izbin][icent][irp][i2], 5, die, idie);
b1f673b5 802 }
803 }
804 ++nmixed;
805 }//for(ibuf)
806 }
807
808 if(ve2.size()>0) {
809 //
810 // Now mixed event for +- pairs
811 //
812 nmixed = 0;
7329335b 813 for(int ibuf=0;(nmixed<fgkNMix);ibuf++) {
b1f673b5 814 int ntry = 0;
620b34a5 815 while(ntry<fgkMAXTRY) {
7329335b 816 ReshuffleBuffer(fvem[ibuf][idie][izbin][icent][irp],fpoolm[idie][izbin][icent][irp]);
b1f673b5 817 ntry++;
818 }
819 for(int i1=0; i1<(int)ve2.size(); i1++){
7329335b 820 for(int i2=0; i2<(int)fvem[ibuf][idie][izbin][icent][irp].size(); i2++){
805bd069 821 FillPair(ve2[i1],fvem[ibuf][idie][izbin][icent][irp][i2],6, die, idie);
b1f673b5 822 }
823 }
824 ++nmixed;
825 }//for(ibuf)
826 }
827
828}
829
830
831//_________________________________________________________________________________
7329335b 832void AliAnalysisTaskMultiDielectronTG::FillPair(AliDielectronSingleTG *iep,
805bd069 833 AliDielectronSingleTG *iem, int type, AliDielectron *die, Int_t idie)
7329335b 834{
835
836 //
837 // main routine for filling kinematics of pairs
838 //
839
840
b1f673b5 841
7329335b 842 double dmass, dphiv, dpxpair, dpypair, dpzpair;
843 double dptpair, depair, dphipair, detapair, dcos, dpsi;
844
845 CalcVars(iep, iem, dmass, dphiv, dpxpair, dpypair, dpzpair,
846 dptpair, depair, dphipair, detapair, dcos, dpsi);
847
848
849 double dopeningangle = -9999;
850 double dv0mass = -9999;
851 double dv0pxpair = -9999;
852 double dv0pypair = -9999;
853 double dv0pzpair = -9999;
854 double dv0ptpair = -9999;
855 double dv0epair = -9999;
856 double dv0xvpair = -9999;
857 double dv0yvpair = -9999;
858 double dv0zvpair = -9999;
859 double dv0phipair = -9999;
860 double dv0etapair = -9999;
861 double dv0rpair = -9999;
862 double dpsipair = -9999;
805bd069 863 double dphivpair = -9999;
b1f673b5 864
865 ////////////////////////////
866 ///// calculate v0 ////////
867 ///////////////////////////
7329335b 868 Bool_t v0OFF=kFALSE;
b1f673b5 869 /// for the moment, this doesn't work for MC
7329335b 870 if(fdv0mixing == kFALSE && (type==3 || type==4 || type==5 || type==6)){
871 v0OFF = kTRUE;
b1f673b5 872 }
873 if(die->GetHasMC()==kTRUE && (type==3 || type==4 || type==5 || type==6)){
7329335b 874 v0OFF = kTRUE;
b1f673b5 875 }
876 if(type==0 || type==1 || type==2){
7329335b 877 v0OFF = kFALSE;
b1f673b5 878 }
879
880
7329335b 881 if(v0OFF==kFALSE){
b1f673b5 882 AliVTrack* trackob1= iep->GetTrack();
883 AliVTrack* trackob2= iem->GetTrack();
884
885 if(!trackob1 || !trackob2){
886 return;
887 }
888
889 AliDielectronPair candidate;
890 candidate.SetTracks(trackob1, (int)(11*iep->Charge()),
891 trackob2, (int)(11*iem->Charge()));
892
893 if(trackob1==trackob2){
894 AliInfo("Same Objects!!");
895 return;
896 }
805bd069 897
898 //const AliKFParticle &kfPairLeg1 = candidate.GetKFFirstDaughter();
899 //const AliKFParticle &kfPairLeg2 = candidate.GetKFSecondDaughter();
900
b1f673b5 901 const AliKFParticle &kfPair = candidate.GetKFParticle();
805bd069 902
903 /*
7329335b 904 dv0mass = candidate.M();
905 dv0pxpair = candidate.Px();
906 dv0pypair = candidate.Py();
907 dv0pzpair = candidate.Pz();
908 dv0ptpair = candidate.Pt();
909 dv0epair = candidate.E();
910 dv0xvpair = candidate.Xv();
911 dv0yvpair = candidate.Yv();
912 dv0zvpair = candidate.Zv();
913 dv0phipair = candidate.Phi();
7329335b 914 dv0etapair = candidate.Eta();
805bd069 915 */
916 dv0mass = kfPair.GetMass();
917 dv0pxpair = kfPair.GetPx();
918 dv0pypair = kfPair.GetPy();
919 dv0pzpair = kfPair.GetPz();
920 dv0ptpair = kfPair.GetPt();
921 dv0epair = kfPair.GetE();
922 dv0xvpair = kfPair.GetX();
923 dv0yvpair = kfPair.GetY();
924 dv0zvpair = kfPair.GetZ();
925 dv0phipair = kfPair.GetPhi();
926 dv0etapair = kfPair.GetEta();
7329335b 927 dv0rpair = kfPair.GetR();
805bd069 928
929 dopeningangle = candidate.OpeningAngle();
7329335b 930 dpsipair = candidate.PsiPair(fbz);
805bd069 931 dphivpair = candidate.PhivPair(fbz);
932
933
b1f673b5 934 }
935
936 Double_t values[AliDielectronVarManager::kNMaxValues];
937 TString className1;
938 TString className2;
7329335b 939 className1.Form("MyPair_%s",kPairClassNamesTG[type]);
940 className2.Form("MyPairV0_%s",kPairClassNamesTG[type]);
b1f673b5 941
942 AliDielectronHistos *fHistos = die->GetHistoManager();
943 Bool_t pairClass1=fHistos->GetHistogramList()->FindObject(className1.Data())!=0x0;
944 Bool_t pairClass2=fHistos->GetHistogramList()->FindObject(className2.Data())!=0x0;
945
805bd069 946 if (pairClass1 && PairTrackcut(dphiv, dcos, idie)==true){
b1f673b5 947 ///import pair variables to values!!
7329335b 948 values[AliDielectronVarManager::kPx] = dpxpair;
949 values[AliDielectronVarManager::kPy] = dpypair;
950 values[AliDielectronVarManager::kPz] = dpzpair;
951 values[AliDielectronVarManager::kPt] = dptpair;
952 values[AliDielectronVarManager::kXv] = dv0xvpair;
953 values[AliDielectronVarManager::kYv] = dv0yvpair;
954 values[AliDielectronVarManager::kZv] = dv0zvpair;
955 values[AliDielectronVarManager::kR] = dv0rpair;
956 values[AliDielectronVarManager::kE] = depair;
957 values[AliDielectronVarManager::kEta] = detapair;
958 values[AliDielectronVarManager::kM] = dmass;
959 values[AliDielectronVarManager::kPsiPair] = dphiv;
805bd069 960 values[AliDielectronVarManager::kPhivPair] = dphiv;
7329335b 961 values[AliDielectronVarManager::kPhi] = dphipair;
962 values[AliDielectronVarManager::kOpeningAngle] = dcos;
b1f673b5 963 fHistos->FillClass(className1, AliDielectronVarManager::kNMaxValues, values);
964 }
965
966
805bd069 967 if (pairClass2 && PairTrackcut(dphiv, dopeningangle, idie)==true){
7329335b 968 values[AliDielectronVarManager::kPx] = dv0pxpair;
969 values[AliDielectronVarManager::kPy] = dv0pypair;
970 values[AliDielectronVarManager::kPz] = dv0pzpair;
971 values[AliDielectronVarManager::kPt] = dv0ptpair;
972 values[AliDielectronVarManager::kXv] = dv0xvpair;
973 values[AliDielectronVarManager::kYv] = dv0yvpair;
974 values[AliDielectronVarManager::kZv] = dv0zvpair;
975 values[AliDielectronVarManager::kR] = dv0rpair;
976 values[AliDielectronVarManager::kE] = dv0epair;
977 values[AliDielectronVarManager::kEta] = dv0etapair;
978 values[AliDielectronVarManager::kM] = dv0mass;
979 values[AliDielectronVarManager::kPsiPair] = dpsipair;
805bd069 980 values[AliDielectronVarManager::kPhivPair] = dphivpair;
7329335b 981 values[AliDielectronVarManager::kPhi] = dv0phipair;
982 values[AliDielectronVarManager::kOpeningAngle] = dopeningangle;
b1f673b5 983 fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
984 }
985
986
987
988}
989
990//_________________________________________________________________________________
805bd069 991bool AliAnalysisTaskMultiDielectronTG::PairTrackcut(double phiv, double op, int idie)
7329335b 992{
993
994 //
995 // pair-by-pair cuts
996 //
b1f673b5 997
998
999 bool pairOK = true;
1000
805bd069 1001 if(fRejectPairFlag[idie] == 1 || fRejectPairFlag[idie] == 2 ||
1002 fRejectPairFlag[idie] == 3 || fRejectPairFlag[idie] == 4 ){
1003 if(fRejectPairFlag[idie] == 2 || fRejectPairFlag[idie] == 4 ){
1004 if(fbz>0 && phiv>fdconvphiv){
1005 pairOK = false;
1006 }else if(fbz<0 && phiv<acos(-1.0)-fdconvphiv){
1007 pairOK = false;
1008 }
1009 }else if(fRejectPairFlag[idie] == 1 || fRejectPairFlag[idie] == 3){
1010 if(op<fdop){
1011 pairOK = false;
1012 }
1013 }
b1f673b5 1014 }
1015
1016 return pairOK;
1017
1018}
1019
1020
1021//_________________________________________________________________________________
7329335b 1022void AliAnalysisTaskMultiDielectronTG::CalcVars(AliDielectronSingleTG *iep, AliDielectronSingleTG *iem,
b1f673b5 1023 double &mass, double &phiv, double &px, double &py, double&pz,
1024 double &pt, double &e, double &phi,
7329335b 1025 double &eta, double &cos, double &psi)
1026{
1027
1028
1029 //
1030 // standalone calculator for the pair variables
1031 //
b1f673b5 1032
1033 px = iep->Px()+iem->Px();
1034 py = iep->Py()+iem->Py();
1035 pz = iep->Pz()+iem->Pz();
1036 pt = sqrt(px*px+py*py);
7329335b 1037 double dppair = sqrt(pt*pt+pz*pz);
b1f673b5 1038 static const double me=0.0005109989;
1039 e = sqrt(me*me+iep->Px()*iep->Px()+iep->Py()*iep->Py()+iep->Pz()*iep->Pz())
1040 + sqrt(me*me+iem->Px()*iem->Px()+iem->Py()*iem->Py()+iem->Pz()*iem->Pz());
1041
1042 mass = e*e-px*px-py*py-pz*pz;
1043 if(mass<0){
1044 mass = mass;
1045 }else{
1046 mass = sqrt(mass);
1047 }
1048
1049
1050 phi = atan2(py, px);
7329335b 1051 eta = -0.5*TMath::Log((dppair+pz)/(dppair-pz));
b1f673b5 1052 double p1 = sqrt(pow(iep->Px(),2)+pow(iep->Py(),2)+pow(iep->Pz(),2));
1053 double p2 = sqrt(pow(iem->Px(),2)+pow(iem->Py(),2)+pow(iem->Pz(),2));
1054 cos = acos((iep->Px()*iem->Px()+iep->Py()*iem->Py()+iep->Pz()*iem->Pz())/(p1*p2));
1055
1056 double dtheta = iep->Theta()-iem->Theta();
1057 psi = asin(dtheta/cos);
1058
1059
1060 //unit vector of (pep+pem)
7329335b 1061 double pl = dppair;
b1f673b5 1062 double ux = px/pl;
1063 double uy = py/pl;
1064 double uz = pz/pl;
1065 double ax = uy/sqrt(ux*ux+uy*uy);
1066 double ay = -ux/sqrt(ux*ux+uy*uy);
1067
1068 //momentum of e+ and e- in (ax,ay,az) axis. Note that az=0 by
1069 //definition.
1070 //double ptep = iep->Px()*ax + iep->Py()*ay;
1071 //double ptem = iem->Px()*ax + iem->Py()*ay;
1072
1073 double pxep = iep->Px();
1074 double pyep = iep->Py();
1075 double pzep = iep->Pz();
1076 double pxem = iem->Px();
1077 double pyem = iem->Py();
1078 double pzem = iem->Pz();
1079
1080
1081 //vector product of pep X pem
1082 double vpx = pyep*pzem - pzep*pyem;
1083 double vpy = pzep*pxem - pxep*pzem;
1084 double vpz = pxep*pyem - pyep*pxem;
1085 double vp = sqrt(vpx*vpx+vpy*vpy+vpz*vpz);
1086 //double thev = acos(vpz/vp);
1087
1088 //unit vector of pep X pem
1089 double vx = vpx/vp;
1090 double vy = vpy/vp;
1091 double vz = vpz/vp;
1092
1093 //The third axis defined by vector product (ux,uy,uz)X(vx,vy,vz)
1094 double wx = uy*vz - uz*vy;
1095 double wy = uz*vx - ux*vz;
1096 double wz = ux*vy - uy*vx;
1097 double wl = sqrt(wx*wx+wy*wy+wz*wz);
1098 // by construction, (wx,wy,wz) must be a unit vector.
1099 if(fabs(wl - 1.0) > 0.00001) std::cout << "Calculation error in W vector"<<std::endl;
1100 // measure angle between (wx,wy,wz) and (ax,ay,0). The angle between them
1101 // should be small if the pair is conversion
1102 //
1103 double cosPhiV = wx*ax + wy*ay;
1104 phiv = acos(cosPhiV);
1105
1106}
1107
1108//_________________________________________________________________________________
7329335b 1109void AliAnalysisTaskMultiDielectronTG::ReshuffleBuffer(vector<AliDielectronSingleTG*> ve, deque<AliDielectronSingleTG*> pool)
1110{
b1f673b5 1111 //If there is not enough electron in the pool, give up
7329335b 1112 //
1113 // ReshuffleBuffer for th event mixing
1114 //
b1f673b5 1115
1116 unsigned int ne = ve.size();
1117 unsigned int poolsize = pool.size();
7329335b 1118 int used[fgkMAXPOOL];
1119 for(int i=0;i<(int)fgkMAXPOOL;i++){
b1f673b5 1120 used[i]=0;
1121 }
1122
1123 if(poolsize < ne) {
1124 std::cout <<" pool size="<<poolsize<<" ne"<<ne<<std::endl;
1125 return;
1126 }
1127 for(unsigned int ie=0; ie < ne; ie++) {
1128 int j = rand()%poolsize;
1129 while(used[j]==1){
1130 j = rand()%poolsize;
1131 }
1132 ve[ie] = pool[j];
1133 used[j]=1;
1134 }
1135
1136}
805bd069 1137
1138//_________________________________________________________________________________
1139void AliAnalysisTaskMultiDielectronTG::RejectPairs(vector<AliDielectronSingleTG*> e1, vector<AliDielectronSingleTG*> e2, Int_t idie){
1140
1141 ////////////////////////////////////
1142 ///// to reject pairs from track arrays
1143 ///// conversions, ghost ..
1144 ///////////////////////////////////
1145 if(e1.size()>0 && e2.size()>0){
1146 for(int i1=0; i1<(int)e1.size(); i1++){
1147 for(int i2=0; i2<(int)e2.size(); i2++){
1148 if(fRejectPairFlag[idie]==1){
1149 Double_t cosine = GetOpeningAngle(e1[i1], e2[i2]);
1150 if(cosine<fdop){
1151 e1[i1]->SetConvFlag(0);
1152 e2[i2]->SetConvFlag(0);
1153 }
1154 }else if(fRejectPairFlag[idie]==2){
1155 Double_t phiv = GetPhiv(e1[i1], e2[i2]);
1156 if(fbz>0 && phiv>fdconvphiv){
1157 e1[i1]->SetConvFlag(0);
1158 e2[i2]->SetConvFlag(0);
1159 }else if(fbz<0 && phiv<acos(-1.0)-fdconvphiv){
1160 e1[i1]->SetConvFlag(0);
1161 e2[i2]->SetConvFlag(0);
1162 }
1163 }
1164 }
1165 }
1166 }
620b34a5 1167
805bd069 1168 if(e1.size()>0){
1169 for(int i1=0; i1<(int)e1.size(); i1++){
1170 for(int i2=i1+1; i2<(int)e1.size(); i2++){
1171 if(fRejectPairFlag[idie]==1){
1172 Double_t cosine = GetOpeningAngle(e1[i1], e1[i2]);
1173 if(cosine<fdop){
1174 e1[i1]->SetConvFlag(0);
1175 e1[i2]->SetConvFlag(0);
1176 }
1177 }else if(fRejectPairFlag[idie]==2){
1178 Double_t phiv = GetPhiv(e1[i1], e1[i2]);
1179 if(fbz>0 && phiv>fdconvphiv){
1180 e1[i1]->SetConvFlag(0);
1181 e1[i2]->SetConvFlag(0);
1182 }else if(fbz<0 && phiv<acos(-1.0)-fdconvphiv){
1183 e1[i1]->SetConvFlag(0);
1184 e1[i2]->SetConvFlag(0);
1185 }
1186 }
1187 }
1188 }
1189 }
1190
1191 if(e2.size()>0){
1192 for(int i1=0; i1<(int)e2.size(); i1++){
1193 for(int i2=i1+1; i2<(int)e2.size(); i2++){
1194 if(fRejectPairFlag[idie]==1){
1195 Double_t cosine = GetOpeningAngle(e2[i1], e2[i2]);
1196 if(cosine<fdop){
1197 e2[i1]->SetConvFlag(0);
1198 e2[i2]->SetConvFlag(0);
1199 }
1200 }else if(fRejectPairFlag[idie]==2){
1201 Double_t phiv = GetPhiv(e2[i1], e2[i2]);
1202 if(fbz>0 && phiv>fdconvphiv){
1203 e2[i1]->SetConvFlag(0);
1204 e2[i2]->SetConvFlag(0);
1205 }else if(fbz<0 && phiv<acos(-1.0)-fdconvphiv){
1206 e2[i1]->SetConvFlag(0);
1207 e2[i2]->SetConvFlag(0);
1208 }
1209 }
1210 }
1211 }
1212 }
1213}
1214
1215
1216//_________________________________________________________________________________
1217Double_t AliAnalysisTaskMultiDielectronTG::GetOpeningAngle(AliDielectronSingleTG* e1, AliDielectronSingleTG* e2){
1218
1219 //////////////////////
1220 //////// calculate pairs and get opening angle
1221 //////////////////////
1222
1223 double dmass, dphiv, dpxpair, dpypair, dpzpair;
1224 double dptpair, depair, dphipair, detapair, dcos, dpsi;
1225
1226 CalcVars(e1, e2, dmass, dphiv, dpxpair, dpypair, dpzpair,
1227 dptpair, depair, dphipair, detapair, dcos, dpsi);
1228
1229 return dcos;
1230}
1231
1232//_________________________________________________________________________________
1233Double_t AliAnalysisTaskMultiDielectronTG::GetPhiv(AliDielectronSingleTG* e1, AliDielectronSingleTG* e2){
1234
1235 //////////////////////
1236 //////// calculate pairs and get opening angle
1237 //////////////////////
1238
1239 double dmass, dphiv, dpxpair, dpypair, dpzpair;
1240 double dptpair, depair, dphipair, detapair, dcos, dpsi;
1241
1242 CalcVars(e1, e2, dmass, dphiv, dpxpair, dpypair, dpzpair,
1243 dptpair, depair, dphipair, detapair, dcos, dpsi);
1244
1245 return dphiv;
1246}