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