]> 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
Merge branch 'feature-movesplit'
[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 52const char* kPairClassNamesTG[7] = {
b1f673b5 53 "unlike",
54 "like_pp",
55 "like_ee",
56 "mixunlike_pe",
57 "mixunlike_ep",
58 "mixlike_pp",
59 "mixlike_ee"
60};
61
62
63ClassImp(AliDielectronSingleTG)
64ClassImp(AliAnalysisTaskMultiDielectronTG)
65
66//_________________________________________________________________________________
67AliAnalysisTaskMultiDielectronTG::AliAnalysisTaskMultiDielectronTG() :
68 AliAnalysisTaskSE(),
69 fListDielectron(),
70 fListHistos(),
71 fListCF(),
7329335b 72 fQAElectron(),
b1f673b5 73 fSelectPhysics(kFALSE),
74 fTriggerMask(AliVEvent::kMB),
75 fExcludeTriggerMask(0),
76 fTriggerOnV0AND(kFALSE),
77 fRejectPileup(kFALSE),
78 fTriggerLogic(kAny),
79 fTriggerAnalysis(0x0),
80 fEventFilter(0x0),
81 fCutsMother(0x0),
82 fEventStat(0x0),
83 fEvent(0x0),
84 fdEdXvsPt(0x0),
85 fdEdXnSigmaElecvsPt(0x0),
86 fdEdXvsPtTOF(0x0),
87 fdEdXnSigmaElecvsPtTOF(0x0),
88 fTOFbetavsPt(0x0),
89 fTOFnSigmaElecvsPt(0x0),
7329335b 90 fNCrossedRowsTPC(0x0),
91 fChi2ClusTPC(0x0),
92 fRatioCrossClusTPC(0x0),
93 fVem(0x0),
94 fVep(0x0),
95 fVemtmp(0x0),
96 fVeptmp(0x0),
97 fdconvphiv(acos(-1.0)),
ebdc847a 98 fdconvMee(100),
805bd069 99 fdop(0),
7329335b 100 fbz(0),
ebdc847a 101 fdv0mixing(kTRUE),
102 fBGRejUnlike(kFALSE),
103 fBGRejLike(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)),
ebdc847a 142 fdconvMee(100),
805bd069 143 fdop(0),
7329335b 144 fbz(0),
ebdc847a 145 fdv0mixing(kTRUE),
146 fBGRejUnlike(kFALSE),
147 fBGRejLike(kTRUE)
b1f673b5 148{
149 //
150 // Constructor
151 //
152 DefineInput(0,TChain::Class());
153 DefineOutput(1, TList::Class());
154 DefineOutput(2, TList::Class());
155 DefineOutput(3, TH1D::Class());
156 fListHistos.SetName("Dielectron_Histos_Multi");
157 fListCF.SetName("Dielectron_CF_Multi");
158 fListDielectron.SetOwner();
159 fListHistos.SetOwner();
160 fListCF.SetOwner();
161
162 ///////////////
7329335b 163 for(int i=0;i<fgkNDIE; i++){
164 for(int j=0;j<fgkNZBIN;j++){
165 for(int k=0;k<fgkNCENT;k++){
166 for(int l=0; l<fgkNRPBIN; l++){
167 fibuf[i][j][k][l] = 0;
168 fpoolp[i][j][k][l].clear();
169 fpoolm[i][j][k][l].clear();
170 for(int ib=0;ib<fgkNBUF; ib++){
171 fvep[ib][i][j][k][l].clear();
172 fvem[ib][i][j][k][l].clear();
b1f673b5 173 }
174 }
175 }
176 }
177 }
178
179
180
181}
182
183//_________________________________________________________________________________
184AliAnalysisTaskMultiDielectronTG::~AliAnalysisTaskMultiDielectronTG()
185{
186 //
187 // Destructor
188 //
189
190 //histograms and CF are owned by the dielectron framework.
191 //however they are streamed to file, so in the first place the
192 //lists need to be owner...
193 fListHistos.SetOwner(kFALSE);
194 fListCF.SetOwner(kFALSE);
195
7329335b 196 for(int i=0;i<fgkNDIE; i++){
197 for(int j=0;j<fgkNZBIN;j++){
198 for(int k=0;k<fgkNCENT;k++){
199 for(int l=0; l<fgkNRPBIN; l++){
200 fibuf[i][j][k][l] = 0;
201 fpoolp[i][j][k][l].clear();
202 fpoolm[i][j][k][l].clear();
203 for(int ib=0;ib<fgkNBUF; ib++){
204 fvep[ib][i][j][k][l].clear();
205 fvem[ib][i][j][k][l].clear();
b1f673b5 206 }
207 }
208 }
209 }
210 }
211}
212//_________________________________________________________________________________
213void AliAnalysisTaskMultiDielectronTG::UserCreateOutputObjects()
214{
215 //
216 // Add all histogram manager histogram lists to the output TList
217 //
218
219 if (!fListHistos.IsEmpty()||!fListCF.IsEmpty()) return; //already initialised
220
221// AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
222// Bool_t isESD=man->GetInputEventHandler()->IsA()==AliESDInputHandler::Class();
223// Bool_t isAOD=man->GetInputEventHandler()->IsA()==AliAODInputHandler::Class();
224
225 TIter nextDie(&fListDielectron);
226 AliDielectron *die=0;
227 while ( (die=static_cast<AliDielectron*>(nextDie())) ){
228 die->Init();
229 if (die->GetHistogramList()) fListHistos.Add(const_cast<THashList*>(die->GetHistogramList()));
230 if (die->GetCFManagerPair()) fListCF.Add(const_cast<AliCFContainer*>(die->GetCFManagerPair()->GetContainer()));
231 }
232
233 Int_t cuts=fListDielectron.GetEntries();
234 Int_t nbins=kNbinsEvent+2*cuts;
235 if (!fEventStat){
236 fEventStat=new TH1D("hEventStat","Event statistics",nbins,0,nbins);
237 fEventStat->GetXaxis()->SetBinLabel(1,"Before Phys. Sel.");
238 fEventStat->GetXaxis()->SetBinLabel(2,"After Phys. Sel.");
239
240 //default names
241 fEventStat->GetXaxis()->SetBinLabel(3,"Bin3 not used");
242 fEventStat->GetXaxis()->SetBinLabel(4,"Bin4 not used");
243 fEventStat->GetXaxis()->SetBinLabel(5,"Bin5 not used");
244
245 if(fTriggerOnV0AND) fEventStat->GetXaxis()->SetBinLabel(3,"V0and triggers");
246 if (fEventFilter) fEventStat->GetXaxis()->SetBinLabel(4,"After Event Filter");
247 if (fRejectPileup) fEventStat->GetXaxis()->SetBinLabel(5,"After Pileup rejection");
248
249 for (Int_t i=0; i<cuts; ++i){
250 fEventStat->GetXaxis()->SetBinLabel((kNbinsEvent+1)+2*i,Form("#splitline{1 candidate}{%s}",fListDielectron.At(i)->GetName()));
251 fEventStat->GetXaxis()->SetBinLabel((kNbinsEvent+2)+2*i,Form("#splitline{With >1 candidate}{%s}",fListDielectron.At(i)->GetName()));
252 }
253 }
254
255 if (!fTriggerAnalysis) fTriggerAnalysis=new AliTriggerAnalysis;
256 fTriggerAnalysis->EnableHistograms();
257 fTriggerAnalysis->SetAnalyzeMC(AliDielectronMC::Instance()->HasMC());
258
259
260
261 int nbinx=400;
7329335b 262 float maxx=20;
263 float minx=0.2;
264 float binw = (TMath::Log(maxx)-TMath::Log(minx))/nbinx;
b1f673b5 265 double xbin[401];
266 for(int ii=0;ii<nbinx+1;ii++){
7329335b 267 xbin[ii] = TMath::Exp(TMath::Log(minx) + 0.5*binw+binw*ii);
b1f673b5 268 }
269
270
7329335b 271 fQAElectron = new TList();
272 fQAElectron->SetName("QAElectron");
273 fQAElectron->SetOwner();
b1f673b5 274
275
276 fEvent = new TH1D("Event","centrality", 100,0,100);
7329335b 277 fQAElectron->Add(fEvent);
b1f673b5 278 fdEdXvsPt = new TH2D("dEdXvsPt","dE/dX vs. PT of TPC", nbinx, xbin, 2000,0,200);
7329335b 279 fQAElectron->Add(fdEdXvsPt);
b1f673b5 280 fdEdXnSigmaElecvsPt = new TH2D("fdEdXnSigmaElecvsPt"," dE/dX normalized to electron vs. pT of TPC",
281 nbinx, xbin, 2000, -10, 10);
7329335b 282 fQAElectron->Add(fdEdXnSigmaElecvsPt);
b1f673b5 283
284 fdEdXvsPtTOF = new TH2D("dEdXvsPtTOF","dE/dX vs. PT of TPC", nbinx, xbin, 2000,0,200);
7329335b 285 fQAElectron->Add(fdEdXvsPtTOF);
b1f673b5 286 fdEdXnSigmaElecvsPtTOF = new TH2D("fdEdXnSigmaElecvsPtTOF"," dE/dX normalized to electron vs. pT of TPC",
287 nbinx, xbin, 2000, -10, 10);
7329335b 288 fQAElectron->Add(fdEdXnSigmaElecvsPtTOF);
b1f673b5 289
290
291
292 fTOFbetavsPt = new TH2D("fTOFbetavsPt","TOF beta vs. p", 400, 0, 20, 1200, 0, 1.2);
7329335b 293 fQAElectron->Add(fTOFbetavsPt);
b1f673b5 294 fTOFnSigmaElecvsPt = new TH2D("fTOFnSigmaElecvsPt","TOF nsigma for electron", 400, 0, 20, 2000, -10, 10);
7329335b 295 fQAElectron->Add(fTOFnSigmaElecvsPt);
b1f673b5 296
7329335b 297 fNCrossedRowsTPC = new TH2F("fNCrossedRowsTPC", "TPC nCrossed Rows vs. pT", 200, 0, 20, 200, 0, 200);
298 fQAElectron->Add(fNCrossedRowsTPC);
299 fChi2ClusTPC = new TH2F("fChi2ClusTPC", "hChi2ClusTPC vs. pT", 200, 0, 20, 200, 0, 10);
300 fQAElectron->Add(fChi2ClusTPC);
b1f673b5 301
7329335b 302 fRatioCrossClusTPC = new TH2F("fRatioCrossClusTPC", "hRatioCrossClusTPC vs. pT", 200, 0, 20, 200, 0, 10);
303 fQAElectron->Add(fRatioCrossClusTPC);
b1f673b5 304
7329335b 305 fListHistos.Add(fQAElectron);
b1f673b5 306
307 fListHistos.SetOwner();
308
309 PostData(1, &fListHistos);
310 PostData(2, &fListCF);
311 PostData(3, fEventStat);
312
313 fCutsMother = new AliESDtrackCuts;
314 fCutsMother->SetDCAToVertex2D(kTRUE);
315 fCutsMother->SetMaxDCAToVertexZ(3.0);
316 fCutsMother->SetMaxDCAToVertexXY(1.0);
317 fCutsMother->SetPtRange( 0.05 , 200.0);
318 fCutsMother->SetEtaRange( -0.84 , 0.84 );
319 fCutsMother->SetAcceptKinkDaughters(kFALSE);
320 fCutsMother->SetRequireITSRefit(kTRUE);
321 fCutsMother->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kFirst);
322 fCutsMother->SetMinNClustersITS(3);
323 fCutsMother->SetRequireTPCRefit(kTRUE);
324
325
805bd069 326 AliInfo(Form("PairCutType = %d %d %d %d %d",
327 fRejectPairFlag[0],
328 fRejectPairFlag[1],
329 fRejectPairFlag[2],
330 fRejectPairFlag[3],
331 fRejectPairFlag[4]));
332
b1f673b5 333}
334
335//_________________________________________________________________________________
336void AliAnalysisTaskMultiDielectronTG::UserExec(Option_t *)
337{
338 //
339 // Main loop. Called for every event
340 //
341
342 if (fListHistos.IsEmpty()&&fListCF.IsEmpty()) return;
343
344 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
345 Bool_t isESD=man->GetInputEventHandler()->IsA()==AliESDInputHandler::Class();
346 Bool_t isAOD=man->GetInputEventHandler()->IsA()==AliAODInputHandler::Class();
347
348 AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
349 if (!inputHandler) return;
350
351// AliPIDResponse *pidRes=inputHandler->GetPIDResponse();
352 if ( inputHandler->GetPIDResponse() ){
353 // for the 2.76 pass2 MC private train. Together with a sigma shift of -0.169
354// pidRes->GetTPCResponse().SetSigma(4.637e-3,2.41332105409873257e+04);
355 AliDielectronVarManager::SetPIDResponse( inputHandler->GetPIDResponse() );
356 } else {
357 AliFatal("This task needs the PID response attached to the input event handler!");
358 }
359
360 // Was event selected ?
361 ULong64_t isSelected = AliVEvent::kAny;
362 Bool_t isRejected = kFALSE;
363 if( fSelectPhysics && inputHandler){
364 if((isESD && inputHandler->GetEventSelection()) || isAOD){
365 isSelected = inputHandler->IsEventSelected();
366 if (fExcludeTriggerMask && (isSelected&fExcludeTriggerMask)) isRejected=kTRUE;
367 if (fTriggerLogic==kAny) isSelected&=fTriggerMask;
368 else if (fTriggerLogic==kExact) isSelected=((isSelected&fTriggerMask)==fTriggerMask);
369 }
370 }
371
372
373 //Before physics selection
374 fEventStat->Fill(kAllEvents);
375 if (isSelected==0||isRejected) {
376 PostData(3,fEventStat);
377 return;
378 }
379 //after physics selection
380 fEventStat->Fill(kSelectedEvents);
381
382 //V0and
383 if(fTriggerOnV0AND){
805bd069 384 if(isESD){if (!fTriggerAnalysis->IsOfflineTriggerFired(static_cast<AliESDEvent*>(InputEvent()), AliTriggerAnalysis::kV0AND))
b1f673b5 385 return;}
805bd069 386 if(isAOD){if(!((static_cast<AliAODEvent*>(InputEvent()))->GetVZEROData()->GetV0ADecision() == AliVVZERO::kV0BB &&
387 (static_cast<AliAODEvent*>(InputEvent()))->GetVZEROData()->GetV0CDecision() == AliVVZERO::kV0BB) )
388 return;}
389 }
b1f673b5 390
391
392 fEventStat->Fill(kV0andEvents);
393
394 //Fill Event histograms before the event filter
395 TIter nextDie(&fListDielectron);
396 AliDielectron *die=0;
397 Double_t values[AliDielectronVarManager::kNMaxValues]={0};
398 Double_t valuesMC[AliDielectronVarManager::kNMaxValues]={0};
399 AliDielectronVarManager::SetEvent(InputEvent());
400 AliDielectronVarManager::Fill(InputEvent(),values);
401 AliDielectronVarManager::Fill(InputEvent(),valuesMC);
402 Bool_t hasMC=AliDielectronMC::Instance()->HasMC();
403 if (hasMC) {
404 if (AliDielectronMC::Instance()->ConnectMCEvent())
405 AliDielectronVarManager::Fill(AliDielectronMC::Instance()->GetMCEvent(),valuesMC);
406 }
407
408 while ( (die=static_cast<AliDielectron*>(nextDie())) ){
409 AliDielectronHistos *h=die->GetHistoManager();
410 if (h){
411 if (h->GetHistogramList()->FindObject("Event_noCuts"))
412 h->FillClass("Event_noCuts",AliDielectronVarManager::kNMaxValues,values);
413 if (hasMC && h->GetHistogramList()->FindObject("MCEvent_noCuts"))
414 h->FillClass("Event_noCuts",AliDielectronVarManager::kNMaxValues,valuesMC);
415 }
416 }
417 nextDie.Reset();
418
419 //event filter
420 if (fEventFilter) {
421 if (!fEventFilter->IsSelected(InputEvent())) return;
422 }
423 fEventStat->Fill(kFilteredEvents);
424
425 //pileup
426 if (fRejectPileup){
427 if (InputEvent()->IsPileupFromSPD(3,0.8,3.,2.,5.)) return;
428 }
429 fEventStat->Fill(kPileupEvents);
430
431 //bz for AliKF
7329335b 432 fbz = InputEvent()->GetMagneticField();
433 AliKFParticle::SetField( fbz );
b1f673b5 434 AliDielectronPID::SetCorrVal((Double_t)InputEvent()->GetRunNumber());
435
436 //Process event in all AliDielectron instances
437 // TIter nextDie(&fListDielectron);
438 // AliDielectron *die=0;
439 Int_t idie=0;
440 while ( (die=static_cast<AliDielectron*>(nextDie())) ){
805bd069 441 //AliInfo(Form(" **************** die->Process(InputEvent()) : %d **************************", idie));
b1f673b5 442 die->SetDontClearArrays(kTRUE);
443 die->Process(InputEvent());
444 if (die->HasCandidates()){
445 Int_t ncandidates=die->GetPairArray(1)->GetEntriesFast();
446 if (ncandidates==1) fEventStat->Fill((kNbinsEvent)+2*idie);
447 else if (ncandidates>1) fEventStat->Fill((kNbinsEvent+1)+2*idie);
448 }
449
450
451 AliDielectronVarManager::Fill(InputEvent(), fgValues);
452 for (Int_t ii=0; ii<2; ++ii){
453 TObjArray *obj = (TObjArray*)die->GetTrackArray(ii);
454 Int_t ntracks=obj->GetEntriesFast();
455 //AliInfo(Form(" ************** # of tracks = %d", ntracks));
456 for (Int_t itrack=0; itrack<ntracks; ++itrack){
457
458 ////////////////////////////////////////////////////////////////////
459 AliDielectronVarManager::Fill(obj->UncheckedAt(itrack), fgValues);
460 ////////////////////////////////////////////////////////////////////
461
462 if(fgValues[AliDielectronVarManager::kCharge]>0){
7329335b 463 fVeptmp.push_back(new AliDielectronSingleTG(1,
b1f673b5 464 fgValues[AliDielectronVarManager::kCentrality],
465 fgValues[AliDielectronVarManager::kXv],
466 fgValues[AliDielectronVarManager::kYv],
467 fgValues[AliDielectronVarManager::kZv],
468 fgValues[AliDielectronVarManager::kPx],
469 fgValues[AliDielectronVarManager::kPy],
470 fgValues[AliDielectronVarManager::kPz],
471 fgValues[AliDielectronVarManager::kPt],
472 fgValues[AliDielectronVarManager::kEta],
473 fgValues[AliDielectronVarManager::kPhi],
474 fgValues[AliDielectronVarManager::kTheta],
475 1, 1, static_cast<AliVTrack*>(obj->UncheckedAt(itrack)))
476 );
477 }else if(fgValues[AliDielectronVarManager::kCharge]<0){
7329335b 478 fVemtmp.push_back(new AliDielectronSingleTG(-1,
b1f673b5 479 fgValues[AliDielectronVarManager::kCentrality],
480 fgValues[AliDielectronVarManager::kXv],
481 fgValues[AliDielectronVarManager::kYv],
482 fgValues[AliDielectronVarManager::kZv],
483 fgValues[AliDielectronVarManager::kPx],
484 fgValues[AliDielectronVarManager::kPy],
485 fgValues[AliDielectronVarManager::kPz],
486 fgValues[AliDielectronVarManager::kPt],
487 fgValues[AliDielectronVarManager::kEta],
488 fgValues[AliDielectronVarManager::kPhi],
489 fgValues[AliDielectronVarManager::kTheta],
490 1, 1, static_cast<AliVTrack*>(obj->UncheckedAt(itrack)))
491 );
492 }
493 }
494 }
7329335b 495 //AliInfo(Form("size of e and p = %d %d", (int)fVeptmp.size(), (int)fVemtmp.size()));
496
b1f673b5 497
7329335b 498 CheckGhostPairs(fVeptmp);
499 CheckGhostPairs(fVemtmp);
805bd069 500 if(fRejectPairFlag[idie]==1 || fRejectPairFlag[idie]==2){
501 RejectPairs(fVeptmp, fVemtmp, idie);
502 }
7329335b 503 RandomizePool(fVeptmp, fVemtmp);
504 CalcPair(fVep, fVem, die, idie);
b1f673b5 505
7329335b 506 // AliInfo(Form("size of e and p (after) = %d %d", (int)fVep.size(), (int)fVem.size()));
b1f673b5 507
7329335b 508 double dwcent = 100.0/fgkNCENT;
509 double dwiz = 20.0/fgkNZBIN;
510 double dwrp = acos(-1.0)/fgkNRPBIN;
b1f673b5 511
7329335b 512 int icent = (int)(fgValues[AliDielectronVarManager::kCentrality]/dwcent);
513 int izbin = (int)((fgValues[AliDielectronVarManager::kZvPrim]+10)/dwiz);
514 int irp = (int)((fgValues[AliDielectronVarManager::kV0ACrpH2])/dwrp);
b1f673b5 515
516 if(icent<0) icent=0;
7329335b 517 if(icent>=fgkNCENT) icent=fgkNCENT-1;
b1f673b5 518 if(izbin<0) izbin=0;
7329335b 519 if(izbin>=fgkNZBIN) izbin=fgkNZBIN-1;
b1f673b5 520 if(irp<0) irp=0;
7329335b 521 if(irp>=fgkNRPBIN) irp=fgkNRPBIN-1;
b1f673b5 522
7329335b 523 fvep[fibuf[idie][izbin][icent][irp]][idie][izbin][icent][irp].clear();
524 for(int iep = 0; iep<(int)fVep.size();iep++) {
525 fvep[fibuf[idie][izbin][icent][irp]][idie][izbin][icent][irp].push_back(fVep[iep]);
526 fpoolp[idie][izbin][icent][irp].push_back(fVep[iep]);
527 if(fpoolp[idie][izbin][icent][irp].size()>fgkMAXPOOL) {
528 fpoolp[idie][izbin][icent][irp].pop_front();
b1f673b5 529 }
530 }
7329335b 531 fvem[fibuf[idie][izbin][icent][irp]][idie][izbin][icent][irp].clear();
532 for(int iem = 0; iem<(int)fVem.size();iem++) {
533 fvem[fibuf[idie][izbin][icent][irp]][idie][izbin][icent][irp].push_back(fVem[iem]);
534 fpoolm[idie][izbin][icent][irp].push_back(fVem[iem]);
535 if(fpoolm[idie][izbin][icent][irp].size()>fgkMAXPOOL) {
536 fpoolm[idie][izbin][icent][irp].pop_front();
b1f673b5 537 }
538 }
539
540
7329335b 541 fibuf[idie][izbin][icent][irp]++;
542 if(fibuf[idie][izbin][icent][irp]>= fgkNBUF) fibuf[idie][izbin][icent][irp]=0;
b1f673b5 543
544
7329335b 545 fVeptmp.clear();
546 fVemtmp.clear();
547 fVep.clear();
548 fVem.clear();
b1f673b5 549
550
551 ++idie;
552 }
553
554
555 AliESDEvent *fESD = dynamic_cast<AliESDEvent*>(InputEvent());
556 fEvent->Fill(values[AliDielectronVarManager::kCentrality]);
557 for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
558 AliESDtrack* track = fESD->GetTrack(iTracks);
559 if (!track) {
560 Printf("ERROR: Could not receive track %d", iTracks);
561 continue;
562 }
563 if(!fCutsMother->AcceptTrack(track)) continue;
564 fdEdXvsPt->Fill(track->GetTPCmomentum(), track->GetTPCsignal());
565 fdEdXnSigmaElecvsPt->Fill(track->GetTPCmomentum(),
566 AliDielectronVarManager::GetPIDResponse()->NumberOfSigmasTPC(track,
567 AliPID::kElectron)
568 -AliDielectronPID::GetCorrVal());
569 /// for beta caliculaton
570 Double_t l = track->GetIntegratedLength(); // cm
571 Double_t t = track->GetTOFsignal();
572 Double_t t0 = AliDielectronVarManager::GetPIDResponse()->GetTOFResponse().GetTimeZero(); // ps
573 Double_t beta = 0;
574 if( (l < 360. || l > 800.) || (t <= 0.) || (t0 >999990.0) ) {
575 beta=-9999;
576 }
577 else {
578 t -= t0; // subtract the T0
579 l *= 0.01; // cm ->m
580 t *= 1e-12; //ps -> s
581
582 Double_t v = l / t;
583 beta = v / TMath::C();
584 }
585
586 fTOFbetavsPt->Fill(track->GetTPCmomentum(), beta);
587 fTOFnSigmaElecvsPt->Fill(track->GetTPCmomentum(),
588 AliDielectronVarManager::GetPIDResponse()->NumberOfSigmasTOF(track,
589 AliPID::kElectron));
590 ////rough PID is required
591 if( fabs(AliDielectronVarManager::GetPIDResponse()->NumberOfSigmasTOF(track, AliPID::kElectron))<3){
592 fdEdXvsPtTOF->Fill(track->GetTPCmomentum(), track->GetTPCsignal());
593 fdEdXnSigmaElecvsPtTOF->Fill(track->GetTPCmomentum(),
594 AliDielectronVarManager::GetPIDResponse()->NumberOfSigmasTPC(track,
595 AliPID::kElectron)
596 -AliDielectronPID::GetCorrVal());
597
598
599 if(track->GetTPCsignal()>70 && track->GetTPCsignal()<90){
7329335b 600 fNCrossedRowsTPC->Fill(track->GetTPCmomentum(), track->GetTPCCrossedRows());
601 //fChi2ClusTPC->Fill(track->GetTPCmomentum(), track->GetTPCchi2()/track->GetTPCNcls());
602 fChi2ClusTPC->Fill(track->GetTPCmomentum(), track->GetTPCchi2());
603 fRatioCrossClusTPC->Fill(track->GetTPCmomentum(), track->GetTPCCrossedRows()/track->GetTPCNclsF());
b1f673b5 604 }
605 }
606 }
607
608 PostData(1, &fListHistos);
609 PostData(2, &fListCF);
610 PostData(3,fEventStat);
611}
612
613//_________________________________________________________________________________
614void AliAnalysisTaskMultiDielectronTG::FinishTaskOutput()
615{
616 //
617 // Write debug tree
618 //
619 TIter nextDie(&fListDielectron);
620 AliDielectron *die=0;
621 while ( (die=static_cast<AliDielectron*>(nextDie())) ){
622 die->SaveDebugTree();
623 AliDielectronMixingHandler *mix=die->GetMixingHandler();
624// printf("\n\n\n===============\ncall mix in Terminate: %p (%p)\n=================\n\n",mix,die);
42c76bea 625 if(!mix) continue;
626 for (Int_t ipool=0; ipool<mix->GetNumberOfBins(); ++ipool){
627 mix->MixRemaining(die, ipool);
628 }
b1f673b5 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;
2ad98fc5 1072 if(mass>=0){
b1f673b5 1073 mass = sqrt(mass);
1074 }
1075
1076
1077 phi = atan2(py, px);
7329335b 1078 eta = -0.5*TMath::Log((dppair+pz)/(dppair-pz));
b1f673b5 1079 double p1 = sqrt(pow(iep->Px(),2)+pow(iep->Py(),2)+pow(iep->Pz(),2));
1080 double p2 = sqrt(pow(iem->Px(),2)+pow(iem->Py(),2)+pow(iem->Pz(),2));
1081 cos = acos((iep->Px()*iem->Px()+iep->Py()*iem->Py()+iep->Pz()*iem->Pz())/(p1*p2));
1082
1083 double dtheta = iep->Theta()-iem->Theta();
1084 psi = asin(dtheta/cos);
1085
1086
1087 //unit vector of (pep+pem)
7329335b 1088 double pl = dppair;
b1f673b5 1089 double ux = px/pl;
1090 double uy = py/pl;
1091 double uz = pz/pl;
1092 double ax = uy/sqrt(ux*ux+uy*uy);
1093 double ay = -ux/sqrt(ux*ux+uy*uy);
1094
1095 //momentum of e+ and e- in (ax,ay,az) axis. Note that az=0 by
1096 //definition.
1097 //double ptep = iep->Px()*ax + iep->Py()*ay;
1098 //double ptem = iem->Px()*ax + iem->Py()*ay;
1099
1100 double pxep = iep->Px();
1101 double pyep = iep->Py();
1102 double pzep = iep->Pz();
1103 double pxem = iem->Px();
1104 double pyem = iem->Py();
1105 double pzem = iem->Pz();
1106
1107
1108 //vector product of pep X pem
1109 double vpx = pyep*pzem - pzep*pyem;
1110 double vpy = pzep*pxem - pxep*pzem;
1111 double vpz = pxep*pyem - pyep*pxem;
1112 double vp = sqrt(vpx*vpx+vpy*vpy+vpz*vpz);
1113 //double thev = acos(vpz/vp);
1114
1115 //unit vector of pep X pem
1116 double vx = vpx/vp;
1117 double vy = vpy/vp;
1118 double vz = vpz/vp;
1119
1120 //The third axis defined by vector product (ux,uy,uz)X(vx,vy,vz)
1121 double wx = uy*vz - uz*vy;
1122 double wy = uz*vx - ux*vz;
1123 double wz = ux*vy - uy*vx;
1124 double wl = sqrt(wx*wx+wy*wy+wz*wz);
1125 // by construction, (wx,wy,wz) must be a unit vector.
1126 if(fabs(wl - 1.0) > 0.00001) std::cout << "Calculation error in W vector"<<std::endl;
1127 // measure angle between (wx,wy,wz) and (ax,ay,0). The angle between them
1128 // should be small if the pair is conversion
1129 //
1130 double cosPhiV = wx*ax + wy*ay;
1131 phiv = acos(cosPhiV);
1132
1133}
1134
1135//_________________________________________________________________________________
7329335b 1136void AliAnalysisTaskMultiDielectronTG::ReshuffleBuffer(vector<AliDielectronSingleTG*> ve, deque<AliDielectronSingleTG*> pool)
1137{
b1f673b5 1138 //If there is not enough electron in the pool, give up
7329335b 1139 //
1140 // ReshuffleBuffer for th event mixing
1141 //
b1f673b5 1142
1143 unsigned int ne = ve.size();
1144 unsigned int poolsize = pool.size();
7329335b 1145 int used[fgkMAXPOOL];
1146 for(int i=0;i<(int)fgkMAXPOOL;i++){
b1f673b5 1147 used[i]=0;
1148 }
1149
1150 if(poolsize < ne) {
1151 std::cout <<" pool size="<<poolsize<<" ne"<<ne<<std::endl;
1152 return;
1153 }
1154 for(unsigned int ie=0; ie < ne; ie++) {
1155 int j = rand()%poolsize;
1156 while(used[j]==1){
1157 j = rand()%poolsize;
1158 }
1159 ve[ie] = pool[j];
1160 used[j]=1;
1161 }
1162
1163}
805bd069 1164
1165//_________________________________________________________________________________
1166void AliAnalysisTaskMultiDielectronTG::RejectPairs(vector<AliDielectronSingleTG*> e1, vector<AliDielectronSingleTG*> e2, Int_t idie){
1167
1168 ////////////////////////////////////
1169 ///// to reject pairs from track arrays
1170 ///// conversions, ghost ..
1171 ///////////////////////////////////
1172 if(e1.size()>0 && e2.size()>0){
1173 for(int i1=0; i1<(int)e1.size(); i1++){
1174 for(int i2=0; i2<(int)e2.size(); i2++){
1175 if(fRejectPairFlag[idie]==1){
1176 Double_t cosine = GetOpeningAngle(e1[i1], e2[i2]);
1177 if(cosine<fdop){
1178 e1[i1]->SetConvFlag(0);
1179 e2[i2]->SetConvFlag(0);
1180 }
1181 }else if(fRejectPairFlag[idie]==2){
1182 Double_t phiv = GetPhiv(e1[i1], e2[i2]);
ebdc847a 1183 Double_t mee = GetMass(e1[i1], e2[i2]);
1184 if(fbz>0 && ( phiv>fdconvphiv && mee < fdconvMee) ){
805bd069 1185 e1[i1]->SetConvFlag(0);
1186 e2[i2]->SetConvFlag(0);
ebdc847a 1187 }else if(fbz<0 && phiv<acos(-1.0)-fdconvphiv && mee < fdconvMee){
805bd069 1188 e1[i1]->SetConvFlag(0);
1189 e2[i2]->SetConvFlag(0);
1190 }
1191 }
1192 }
1193 }
1194 }
805bd069 1195 if(e1.size()>0){
1196 for(int i1=0; i1<(int)e1.size(); i1++){
1197 for(int i2=i1+1; i2<(int)e1.size(); i2++){
1198 if(fRejectPairFlag[idie]==1){
1199 Double_t cosine = GetOpeningAngle(e1[i1], e1[i2]);
1200 if(cosine<fdop){
1201 e1[i1]->SetConvFlag(0);
1202 e1[i2]->SetConvFlag(0);
1203 }
1204 }else if(fRejectPairFlag[idie]==2){
1205 Double_t phiv = GetPhiv(e1[i1], e1[i2]);
ebdc847a 1206 Double_t mee = GetMass(e1[i1], e1[i2]);
1207 if(fbz>0 && phiv>fdconvphiv && mee < fdconvMee){
805bd069 1208 e1[i1]->SetConvFlag(0);
1209 e1[i2]->SetConvFlag(0);
ebdc847a 1210 }else if(fbz<0 && phiv<acos(-1.0)-fdconvphiv && mee < fdconvMee){
805bd069 1211 e1[i1]->SetConvFlag(0);
1212 e1[i2]->SetConvFlag(0);
1213 }
1214 }
1215 }
1216 }
1217 }
1218
1219 if(e2.size()>0){
1220 for(int i1=0; i1<(int)e2.size(); i1++){
1221 for(int i2=i1+1; i2<(int)e2.size(); i2++){
1222 if(fRejectPairFlag[idie]==1){
1223 Double_t cosine = GetOpeningAngle(e2[i1], e2[i2]);
1224 if(cosine<fdop){
1225 e2[i1]->SetConvFlag(0);
1226 e2[i2]->SetConvFlag(0);
1227 }
1228 }else if(fRejectPairFlag[idie]==2){
1229 Double_t phiv = GetPhiv(e2[i1], e2[i2]);
ebdc847a 1230 Double_t mee = GetMass(e2[i1], e2[i2]);
1231 if(fbz>0 && phiv>fdconvphiv && mee < fdconvMee){
805bd069 1232 e2[i1]->SetConvFlag(0);
1233 e2[i2]->SetConvFlag(0);
ebdc847a 1234 }else if(fbz<0 && phiv<acos(-1.0)-fdconvphiv && mee < fdconvMee){
805bd069 1235 e2[i1]->SetConvFlag(0);
1236 e2[i2]->SetConvFlag(0);
1237 }
1238 }
1239 }
1240 }
1241 }
1242}
1243
1244
1245//_________________________________________________________________________________
1246Double_t AliAnalysisTaskMultiDielectronTG::GetOpeningAngle(AliDielectronSingleTG* e1, AliDielectronSingleTG* e2){
1247
1248 //////////////////////
1249 //////// calculate pairs and get opening angle
1250 //////////////////////
1251
1252 double dmass, dphiv, dpxpair, dpypair, dpzpair;
1253 double dptpair, depair, dphipair, detapair, dcos, dpsi;
1254
1255 CalcVars(e1, e2, dmass, dphiv, dpxpair, dpypair, dpzpair,
1256 dptpair, depair, dphipair, detapair, dcos, dpsi);
1257
1258 return dcos;
1259}
1260
1261//_________________________________________________________________________________
1262Double_t AliAnalysisTaskMultiDielectronTG::GetPhiv(AliDielectronSingleTG* e1, AliDielectronSingleTG* e2){
1263
1264 //////////////////////
ebdc847a 1265 //////// calculate pairs and get phiv
805bd069 1266 //////////////////////
1267
1268 double dmass, dphiv, dpxpair, dpypair, dpzpair;
1269 double dptpair, depair, dphipair, detapair, dcos, dpsi;
1270
1271 CalcVars(e1, e2, dmass, dphiv, dpxpair, dpypair, dpzpair,
1272 dptpair, depair, dphipair, detapair, dcos, dpsi);
1273
1274 return dphiv;
1275}
ebdc847a 1276
1277//_________________________________________________________________________________
1278Double_t AliAnalysisTaskMultiDielectronTG::GetMass(AliDielectronSingleTG* e1, AliDielectronSingleTG* e2){
1279
1280 //////////////////////
1281 //////// calculate pairs and get mass
1282 //////////////////////
1283
1284 double dmass, dphiv, dpxpair, dpypair, dpzpair;
1285 double dptpair, depair, dphipair, detapair, dcos, dpsi;
1286
1287 CalcVars(e1, e2, dmass, dphiv, dpxpair, dpypair, dpzpair,
1288 dptpair, depair, dphipair, detapair, dcos, dpsi);
1289
1290 return dmass;
1291}