]> 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
Bug in poisition corrected
[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);
625 if (mix) mix->MixRemaining(die);
626 }
627 PostData(1, &fListHistos);
628 PostData(2, &fListCF);
629}
630
631//_________________________________________________________________________________
7329335b 632void AliAnalysisTaskMultiDielectronTG::CheckGhostPairs(vector<AliDielectronSingleTG*> e1)
633{
634 ////// Check whether the pairs are in ghost pairs
635 ///// ghost means that one of the pairs is highly associated but reconstructed as a true particle
636
637
b1f673b5 638 bool reject = false;
639 if(e1.size()>1){
640 for(int i1=0; i1<(int)e1.size(); i1++){
641 reject = false;
642 for(int i2=i1+1; i2<(int)e1.size(); i2++){
ebdc847a 643 if( fabs(e1[i1]->Phi() - e1[i2]->Phi())<0.01 &&
644 fabs(e1[i1]->Eta() - e1[i2]->Eta())<0.005
645 ){
b1f673b5 646 reject = true;
647 e1[i2]->SetGstFlag(0);
648 }
649 }
650 if(reject==true)e1[i1]->SetGstFlag(0);
651 }
652 }
653}
654
ebdc847a 655//_________________________________________________________________________________
656Bool_t AliAnalysisTaskMultiDielectronTG::CheckGhost(vector<AliDielectronSingleTG*> e1, vector<AliDielectronSingleTG*> e2)
657{
658 ////// To be sure whether there are no ghost pairs in h event mixing
659
660 if(e1.size()>0 && e2.size()>0){
661 for(int i1=0; i1<(int)e1.size(); i1++){
662 for(int i2=0; i2<(int)e2.size(); i2++){
663 if( fabs(e1[i1]->Phi() - e2[i2]->Phi())<0.01 &&
664 fabs(e1[i1]->Eta() - e2[i2]->Eta())<0.005
665 ){
666 return true;
667 }
668 }
669 }
670 }
671 return false;
672}
673
b1f673b5 674//_________________________________________________________________________________
7329335b 675void AliAnalysisTaskMultiDielectronTG::RandomizePool(vector<AliDielectronSingleTG*> e1, vector<AliDielectronSingleTG*> e2)
676{
b1f673b5 677
7329335b 678 ///// Randomize the pool constent to cancel the filling scheme of the single tracks
679 ////
680
681
b1f673b5 682 int size1 = e1.size();
7329335b 683 int usedindex[1000];
b1f673b5 684 for(int i=0;i<1000;i++){
7329335b 685 usedindex[i] = -1;
b1f673b5 686 }
687 for(int i=0;i<size1;i++){
7329335b 688 usedindex[i] = 0;
b1f673b5 689 }
690
691 for(int i=0;i<size1;i++){
692 int j = (int)(gRandom->Uniform(0,size1));
7329335b 693 while(usedindex[j]==1){
b1f673b5 694 j = (int)(gRandom->Uniform(0,size1));
695 }
696 if( (e1[j]->GetGstFlag()==1) &&
697 (e1[j]->GetConvFlag()==1)
698 ){
7329335b 699 fVep.push_back(e1[j]);
b1f673b5 700 }
7329335b 701 usedindex[j] = 1;
b1f673b5 702 }
703
704
705 int size2 = e2.size();
706 for(int i=0;i<1000;i++){
7329335b 707 usedindex[i] = -1;
b1f673b5 708 }
709 for(int i=0;i<size2;i++){
7329335b 710 usedindex[i] = 0;
b1f673b5 711 }
712
713 for(int i=0;i<size2;i++){
714 int j = (int)(gRandom->Uniform(0,size2));
7329335b 715 while(usedindex[j]==1){
b1f673b5 716 j = (int)(gRandom->Uniform(0,size2));
717 }
718 if( (e2[j]->GetGstFlag()==1) &&
719 (e2[j]->GetConvFlag()==1)
720 ){
7329335b 721 fVem.push_back(e2[j]);
b1f673b5 722 }
7329335b 723 usedindex[j] = 1;
b1f673b5 724 }
725}
726
727
728//_________________________________________________________________________________
7329335b 729void AliAnalysisTaskMultiDielectronTG::CalcPair(vector<AliDielectronSingleTG*> ve1,
730 vector<AliDielectronSingleTG*> ve2, AliDielectron *die, Int_t idie)
731{
732
733 //
734 // main routine for the pair procesing
735 //
736
b1f673b5 737
738 for(int i1=0; i1<(int)ve1.size(); i1++){
739 for(int i2=0; i2<(int)ve2.size(); i2++){
805bd069 740 FillPair(ve1[i1], ve2[i2], 0, die, idie);
b1f673b5 741 }
742 }
743
744 for(int i1=0; i1<(int)ve1.size(); i1++){
745 for(int i2=i1+1; i2<(int)ve1.size(); i2++){
805bd069 746 FillPair(ve1[i1], ve1[i2], 1, die, idie);
b1f673b5 747 }
748 }
749
750 for(int i1=0; i1<(int)ve2.size(); i1++){
751 for(int i2=i1+1; i2<(int)ve2.size(); i2++){
805bd069 752 FillPair(ve2[i1], ve2[i2], 2, die, idie );
b1f673b5 753 }
754 }
755
756
7329335b 757 double dwcent = 100.0/fgkNCENT;
758 double dwiz = 20.0/fgkNZBIN;
759 double dwrp = acos(-1.0)/fgkNRPBIN;
b1f673b5 760
7329335b 761 int icent = (int)(fgValues[AliDielectronVarManager::kCentrality]/dwcent);
762 int izbin = (int)((fgValues[AliDielectronVarManager::kZvPrim]+10)/dwiz);
763 int irp = (int)((fgValues[AliDielectronVarManager::kV0ACrpH2])/dwrp);
b1f673b5 764
765 if(icent<0) icent=0;
7329335b 766 if(icent>=fgkNCENT) icent=fgkNCENT-1;
b1f673b5 767 if(izbin<0) izbin=0;
7329335b 768 if(izbin>=fgkNZBIN) izbin=fgkNZBIN-1;
b1f673b5 769 if(irp<0) irp=0;
7329335b 770 if(irp>=fgkNRPBIN) irp=fgkNRPBIN-1;
b1f673b5 771
772 int nmixed;
773 if(ve1.size()>0) {
774 //
775 // Now mixed event for +- pairs
776 //
777 nmixed = 0;
7329335b 778 for(int ibuf=0;(nmixed<fgkNMix);ibuf++) {
b1f673b5 779 int ntry = 0;
ebdc847a 780 while((fBGRejUnlike && CheckGhost(ve1, fvem[ibuf][idie][izbin][icent][irp])) && ntry<fgkMAXTRY) {
7329335b 781 ReshuffleBuffer(fvem[ibuf][idie][izbin][icent][irp],fpoolm[idie][izbin][icent][irp]);
b1f673b5 782 ntry++;
783 }
784 for(int i1=0; i1<(int)ve1.size(); i1++){
7329335b 785 for(int i2=0; i2<(int)fvem[ibuf][idie][izbin][icent][irp].size(); i2++){
805bd069 786 FillPair(ve1[i1],fvem[ibuf][idie][izbin][icent][irp][i2], 3, die, idie);
b1f673b5 787 }
788 }
789 ++nmixed;
790 }//for(ibuf)
791 }
792 if(ve2.size()>0) {
793 //
805bd069 794 // Now mixed event for -+ pairs
b1f673b5 795 //
796 nmixed = 0;
7329335b 797 for(int ibuf=0;(nmixed<fgkNMix);ibuf++) {
b1f673b5 798 int ntry = 0;
ebdc847a 799 while((fBGRejUnlike && CheckGhost(ve2, fvep[ibuf][idie][izbin][icent][irp])) && ntry<fgkMAXTRY) {
7329335b 800 ReshuffleBuffer(fvep[ibuf][idie][izbin][icent][irp],fpoolp[idie][izbin][icent][irp]);
b1f673b5 801 ntry++;
802 }
803 for(int i1=0; i1<(int)ve2.size(); i1++){
7329335b 804 for(int i2=0; i2<(int)fvep[ibuf][idie][izbin][icent][irp].size(); i2++){
805bd069 805 FillPair(fvep[ibuf][idie][izbin][icent][irp][i2],ve2[i1],4, die, idie);
b1f673b5 806 }
807 }
808 ++nmixed;
809 }//for(ibuf)
810 }
811
812 if(ve1.size()>0) {
813 //
814 // Now mixed event for ++ pairs
815 //
816 nmixed = 0;
7329335b 817 for(int ibuf=0;(nmixed<fgkNMix);ibuf++) {
b1f673b5 818 int ntry = 0;
ebdc847a 819 while((fBGRejLike && CheckGhost(ve1, fvep[ibuf][idie][izbin][icent][irp])) && ntry<fgkMAXTRY) {
7329335b 820 ReshuffleBuffer(fvep[ibuf][idie][izbin][icent][irp],fpoolp[idie][izbin][icent][irp]);
b1f673b5 821 ntry++;
822 }
823 for(int i1=0; i1<(int)ve1.size(); i1++){
7329335b 824 for(int i2=0;i2<(int)fvep[ibuf][idie][izbin][icent][irp].size(); i2++){
805bd069 825 FillPair(ve1[i1],fvep[ibuf][idie][izbin][icent][irp][i2], 5, die, idie);
b1f673b5 826 }
827 }
828 ++nmixed;
829 }//for(ibuf)
830 }
831
832 if(ve2.size()>0) {
833 //
834 // Now mixed event for +- pairs
835 //
836 nmixed = 0;
7329335b 837 for(int ibuf=0;(nmixed<fgkNMix);ibuf++) {
b1f673b5 838 int ntry = 0;
ebdc847a 839 while((fBGRejLike && CheckGhost(ve2, fvem[ibuf][idie][izbin][icent][irp])) && ntry<fgkMAXTRY) {
7329335b 840 ReshuffleBuffer(fvem[ibuf][idie][izbin][icent][irp],fpoolm[idie][izbin][icent][irp]);
b1f673b5 841 ntry++;
842 }
843 for(int i1=0; i1<(int)ve2.size(); i1++){
7329335b 844 for(int i2=0; i2<(int)fvem[ibuf][idie][izbin][icent][irp].size(); i2++){
805bd069 845 FillPair(ve2[i1],fvem[ibuf][idie][izbin][icent][irp][i2],6, die, idie);
b1f673b5 846 }
847 }
848 ++nmixed;
849 }//for(ibuf)
850 }
851
852}
853
854
855//_________________________________________________________________________________
7329335b 856void AliAnalysisTaskMultiDielectronTG::FillPair(AliDielectronSingleTG *iep,
805bd069 857 AliDielectronSingleTG *iem, int type, AliDielectron *die, Int_t idie)
7329335b 858{
859
860 //
861 // main routine for filling kinematics of pairs
862 //
863
864
b1f673b5 865
7329335b 866 double dmass, dphiv, dpxpair, dpypair, dpzpair;
867 double dptpair, depair, dphipair, detapair, dcos, dpsi;
868
869 CalcVars(iep, iem, dmass, dphiv, dpxpair, dpypair, dpzpair,
870 dptpair, depair, dphipair, detapair, dcos, dpsi);
871
872
873 double dopeningangle = -9999;
874 double dv0mass = -9999;
875 double dv0pxpair = -9999;
876 double dv0pypair = -9999;
877 double dv0pzpair = -9999;
878 double dv0ptpair = -9999;
879 double dv0epair = -9999;
880 double dv0xvpair = -9999;
881 double dv0yvpair = -9999;
882 double dv0zvpair = -9999;
883 double dv0phipair = -9999;
884 double dv0etapair = -9999;
885 double dv0rpair = -9999;
886 double dpsipair = -9999;
805bd069 887 double dphivpair = -9999;
b1f673b5 888
889 ////////////////////////////
890 ///// calculate v0 ////////
891 ///////////////////////////
7329335b 892 Bool_t v0OFF=kFALSE;
b1f673b5 893 /// for the moment, this doesn't work for MC
7329335b 894 if(fdv0mixing == kFALSE && (type==3 || type==4 || type==5 || type==6)){
895 v0OFF = kTRUE;
b1f673b5 896 }
897 if(die->GetHasMC()==kTRUE && (type==3 || type==4 || type==5 || type==6)){
7329335b 898 v0OFF = kTRUE;
b1f673b5 899 }
900 if(type==0 || type==1 || type==2){
7329335b 901 v0OFF = kFALSE;
b1f673b5 902 }
903
904
7329335b 905 if(v0OFF==kFALSE){
b1f673b5 906 AliVTrack* trackob1= iep->GetTrack();
907 AliVTrack* trackob2= iem->GetTrack();
908
909 if(!trackob1 || !trackob2){
910 return;
911 }
912
913 AliDielectronPair candidate;
914 candidate.SetTracks(trackob1, (int)(11*iep->Charge()),
915 trackob2, (int)(11*iem->Charge()));
916
917 if(trackob1==trackob2){
918 AliInfo("Same Objects!!");
919 return;
920 }
805bd069 921
922 //const AliKFParticle &kfPairLeg1 = candidate.GetKFFirstDaughter();
923 //const AliKFParticle &kfPairLeg2 = candidate.GetKFSecondDaughter();
924
b1f673b5 925 const AliKFParticle &kfPair = candidate.GetKFParticle();
805bd069 926
927 /*
7329335b 928 dv0mass = candidate.M();
929 dv0pxpair = candidate.Px();
930 dv0pypair = candidate.Py();
931 dv0pzpair = candidate.Pz();
932 dv0ptpair = candidate.Pt();
933 dv0epair = candidate.E();
934 dv0xvpair = candidate.Xv();
935 dv0yvpair = candidate.Yv();
936 dv0zvpair = candidate.Zv();
937 dv0phipair = candidate.Phi();
7329335b 938 dv0etapair = candidate.Eta();
805bd069 939 */
940 dv0mass = kfPair.GetMass();
941 dv0pxpair = kfPair.GetPx();
942 dv0pypair = kfPair.GetPy();
943 dv0pzpair = kfPair.GetPz();
944 dv0ptpair = kfPair.GetPt();
945 dv0epair = kfPair.GetE();
946 dv0xvpair = kfPair.GetX();
947 dv0yvpair = kfPair.GetY();
948 dv0zvpair = kfPair.GetZ();
949 dv0phipair = kfPair.GetPhi();
950 dv0etapair = kfPair.GetEta();
7329335b 951 dv0rpair = kfPair.GetR();
805bd069 952
953 dopeningangle = candidate.OpeningAngle();
7329335b 954 dpsipair = candidate.PsiPair(fbz);
805bd069 955 dphivpair = candidate.PhivPair(fbz);
956
957
b1f673b5 958 }
959
960 Double_t values[AliDielectronVarManager::kNMaxValues];
961 TString className1;
962 TString className2;
7329335b 963 className1.Form("MyPair_%s",kPairClassNamesTG[type]);
964 className2.Form("MyPairV0_%s",kPairClassNamesTG[type]);
b1f673b5 965
966 AliDielectronHistos *fHistos = die->GetHistoManager();
967 Bool_t pairClass1=fHistos->GetHistogramList()->FindObject(className1.Data())!=0x0;
968 Bool_t pairClass2=fHistos->GetHistogramList()->FindObject(className2.Data())!=0x0;
969
ebdc847a 970 if (pairClass1 && PairTrackcut(dphiv, dcos, dmass, idie)==true){
b1f673b5 971 ///import pair variables to values!!
7329335b 972 values[AliDielectronVarManager::kPx] = dpxpair;
973 values[AliDielectronVarManager::kPy] = dpypair;
974 values[AliDielectronVarManager::kPz] = dpzpair;
975 values[AliDielectronVarManager::kPt] = dptpair;
976 values[AliDielectronVarManager::kXv] = dv0xvpair;
977 values[AliDielectronVarManager::kYv] = dv0yvpair;
978 values[AliDielectronVarManager::kZv] = dv0zvpair;
979 values[AliDielectronVarManager::kR] = dv0rpair;
980 values[AliDielectronVarManager::kE] = depair;
981 values[AliDielectronVarManager::kEta] = detapair;
982 values[AliDielectronVarManager::kM] = dmass;
983 values[AliDielectronVarManager::kPsiPair] = dphiv;
805bd069 984 values[AliDielectronVarManager::kPhivPair] = dphiv;
7329335b 985 values[AliDielectronVarManager::kPhi] = dphipair;
986 values[AliDielectronVarManager::kOpeningAngle] = dcos;
ebdc847a 987 values[AliDielectronVarManager::kCosPointingAngle] = TMath::Abs(TMath::ATan2(TMath::Sin(iep->Phi()-iem->Phi()),TMath::Cos(iep->Phi()-iem->Phi())));
b1f673b5 988 fHistos->FillClass(className1, AliDielectronVarManager::kNMaxValues, values);
989 }
990
991
ebdc847a 992 if (pairClass2 && PairTrackcut(dphiv, dopeningangle, dv0mass, idie)==true){
7329335b 993 values[AliDielectronVarManager::kPx] = dv0pxpair;
994 values[AliDielectronVarManager::kPy] = dv0pypair;
995 values[AliDielectronVarManager::kPz] = dv0pzpair;
996 values[AliDielectronVarManager::kPt] = dv0ptpair;
997 values[AliDielectronVarManager::kXv] = dv0xvpair;
998 values[AliDielectronVarManager::kYv] = dv0yvpair;
999 values[AliDielectronVarManager::kZv] = dv0zvpair;
1000 values[AliDielectronVarManager::kR] = dv0rpair;
1001 values[AliDielectronVarManager::kE] = dv0epair;
1002 values[AliDielectronVarManager::kEta] = dv0etapair;
1003 values[AliDielectronVarManager::kM] = dv0mass;
1004 values[AliDielectronVarManager::kPsiPair] = dpsipair;
805bd069 1005 values[AliDielectronVarManager::kPhivPair] = dphivpair;
7329335b 1006 values[AliDielectronVarManager::kPhi] = dv0phipair;
1007 values[AliDielectronVarManager::kOpeningAngle] = dopeningangle;
ebdc847a 1008 values[AliDielectronVarManager::kCosPointingAngle] = TMath::Abs(TMath::ATan2(TMath::Sin(iep->Phi()-iem->Phi()),TMath::Cos(iep->Phi()-iem->Phi())));
b1f673b5 1009 fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
1010 }
1011
1012
1013
1014}
1015
1016//_________________________________________________________________________________
ebdc847a 1017bool AliAnalysisTaskMultiDielectronTG::PairTrackcut(double phiv, double op, double mass, int idie)
7329335b 1018{
1019
1020 //
1021 // pair-by-pair cuts
1022 //
b1f673b5 1023
1024
1025 bool pairOK = true;
1026
805bd069 1027 if(fRejectPairFlag[idie] == 1 || fRejectPairFlag[idie] == 2 ||
1028 fRejectPairFlag[idie] == 3 || fRejectPairFlag[idie] == 4 ){
1029 if(fRejectPairFlag[idie] == 2 || fRejectPairFlag[idie] == 4 ){
ebdc847a 1030 if(fbz>0 && (phiv>fdconvphiv && mass < fdconvMee) ){
805bd069 1031 pairOK = false;
ebdc847a 1032 }else if(fbz<0 && phiv<acos(-1.0)-fdconvphiv && mass < fdconvMee){
805bd069 1033 pairOK = false;
1034 }
1035 }else if(fRejectPairFlag[idie] == 1 || fRejectPairFlag[idie] == 3){
1036 if(op<fdop){
1037 pairOK = false;
1038 }
1039 }
b1f673b5 1040 }
1041
1042 return pairOK;
1043
1044}
1045
1046
1047//_________________________________________________________________________________
7329335b 1048void AliAnalysisTaskMultiDielectronTG::CalcVars(AliDielectronSingleTG *iep, AliDielectronSingleTG *iem,
b1f673b5 1049 double &mass, double &phiv, double &px, double &py, double&pz,
1050 double &pt, double &e, double &phi,
7329335b 1051 double &eta, double &cos, double &psi)
1052{
1053
1054
1055 //
1056 // standalone calculator for the pair variables
1057 //
b1f673b5 1058
1059 px = iep->Px()+iem->Px();
1060 py = iep->Py()+iem->Py();
1061 pz = iep->Pz()+iem->Pz();
1062 pt = sqrt(px*px+py*py);
7329335b 1063 double dppair = sqrt(pt*pt+pz*pz);
b1f673b5 1064 static const double me=0.0005109989;
1065 e = sqrt(me*me+iep->Px()*iep->Px()+iep->Py()*iep->Py()+iep->Pz()*iep->Pz())
1066 + sqrt(me*me+iem->Px()*iem->Px()+iem->Py()*iem->Py()+iem->Pz()*iem->Pz());
1067
1068 mass = e*e-px*px-py*py-pz*pz;
2ad98fc5 1069 if(mass>=0){
b1f673b5 1070 mass = sqrt(mass);
1071 }
1072
1073
1074 phi = atan2(py, px);
7329335b 1075 eta = -0.5*TMath::Log((dppair+pz)/(dppair-pz));
b1f673b5 1076 double p1 = sqrt(pow(iep->Px(),2)+pow(iep->Py(),2)+pow(iep->Pz(),2));
1077 double p2 = sqrt(pow(iem->Px(),2)+pow(iem->Py(),2)+pow(iem->Pz(),2));
1078 cos = acos((iep->Px()*iem->Px()+iep->Py()*iem->Py()+iep->Pz()*iem->Pz())/(p1*p2));
1079
1080 double dtheta = iep->Theta()-iem->Theta();
1081 psi = asin(dtheta/cos);
1082
1083
1084 //unit vector of (pep+pem)
7329335b 1085 double pl = dppair;
b1f673b5 1086 double ux = px/pl;
1087 double uy = py/pl;
1088 double uz = pz/pl;
1089 double ax = uy/sqrt(ux*ux+uy*uy);
1090 double ay = -ux/sqrt(ux*ux+uy*uy);
1091
1092 //momentum of e+ and e- in (ax,ay,az) axis. Note that az=0 by
1093 //definition.
1094 //double ptep = iep->Px()*ax + iep->Py()*ay;
1095 //double ptem = iem->Px()*ax + iem->Py()*ay;
1096
1097 double pxep = iep->Px();
1098 double pyep = iep->Py();
1099 double pzep = iep->Pz();
1100 double pxem = iem->Px();
1101 double pyem = iem->Py();
1102 double pzem = iem->Pz();
1103
1104
1105 //vector product of pep X pem
1106 double vpx = pyep*pzem - pzep*pyem;
1107 double vpy = pzep*pxem - pxep*pzem;
1108 double vpz = pxep*pyem - pyep*pxem;
1109 double vp = sqrt(vpx*vpx+vpy*vpy+vpz*vpz);
1110 //double thev = acos(vpz/vp);
1111
1112 //unit vector of pep X pem
1113 double vx = vpx/vp;
1114 double vy = vpy/vp;
1115 double vz = vpz/vp;
1116
1117 //The third axis defined by vector product (ux,uy,uz)X(vx,vy,vz)
1118 double wx = uy*vz - uz*vy;
1119 double wy = uz*vx - ux*vz;
1120 double wz = ux*vy - uy*vx;
1121 double wl = sqrt(wx*wx+wy*wy+wz*wz);
1122 // by construction, (wx,wy,wz) must be a unit vector.
1123 if(fabs(wl - 1.0) > 0.00001) std::cout << "Calculation error in W vector"<<std::endl;
1124 // measure angle between (wx,wy,wz) and (ax,ay,0). The angle between them
1125 // should be small if the pair is conversion
1126 //
1127 double cosPhiV = wx*ax + wy*ay;
1128 phiv = acos(cosPhiV);
1129
1130}
1131
1132//_________________________________________________________________________________
7329335b 1133void AliAnalysisTaskMultiDielectronTG::ReshuffleBuffer(vector<AliDielectronSingleTG*> ve, deque<AliDielectronSingleTG*> pool)
1134{
b1f673b5 1135 //If there is not enough electron in the pool, give up
7329335b 1136 //
1137 // ReshuffleBuffer for th event mixing
1138 //
b1f673b5 1139
1140 unsigned int ne = ve.size();
1141 unsigned int poolsize = pool.size();
7329335b 1142 int used[fgkMAXPOOL];
1143 for(int i=0;i<(int)fgkMAXPOOL;i++){
b1f673b5 1144 used[i]=0;
1145 }
1146
1147 if(poolsize < ne) {
1148 std::cout <<" pool size="<<poolsize<<" ne"<<ne<<std::endl;
1149 return;
1150 }
1151 for(unsigned int ie=0; ie < ne; ie++) {
1152 int j = rand()%poolsize;
1153 while(used[j]==1){
1154 j = rand()%poolsize;
1155 }
1156 ve[ie] = pool[j];
1157 used[j]=1;
1158 }
1159
1160}
805bd069 1161
1162//_________________________________________________________________________________
1163void AliAnalysisTaskMultiDielectronTG::RejectPairs(vector<AliDielectronSingleTG*> e1, vector<AliDielectronSingleTG*> e2, Int_t idie){
1164
1165 ////////////////////////////////////
1166 ///// to reject pairs from track arrays
1167 ///// conversions, ghost ..
1168 ///////////////////////////////////
1169 if(e1.size()>0 && e2.size()>0){
1170 for(int i1=0; i1<(int)e1.size(); i1++){
1171 for(int i2=0; i2<(int)e2.size(); i2++){
1172 if(fRejectPairFlag[idie]==1){
1173 Double_t cosine = GetOpeningAngle(e1[i1], e2[i2]);
1174 if(cosine<fdop){
1175 e1[i1]->SetConvFlag(0);
1176 e2[i2]->SetConvFlag(0);
1177 }
1178 }else if(fRejectPairFlag[idie]==2){
1179 Double_t phiv = GetPhiv(e1[i1], e2[i2]);
ebdc847a 1180 Double_t mee = GetMass(e1[i1], e2[i2]);
1181 if(fbz>0 && ( phiv>fdconvphiv && mee < fdconvMee) ){
805bd069 1182 e1[i1]->SetConvFlag(0);
1183 e2[i2]->SetConvFlag(0);
ebdc847a 1184 }else if(fbz<0 && phiv<acos(-1.0)-fdconvphiv && mee < fdconvMee){
805bd069 1185 e1[i1]->SetConvFlag(0);
1186 e2[i2]->SetConvFlag(0);
1187 }
1188 }
1189 }
1190 }
1191 }
805bd069 1192 if(e1.size()>0){
1193 for(int i1=0; i1<(int)e1.size(); i1++){
1194 for(int i2=i1+1; i2<(int)e1.size(); i2++){
1195 if(fRejectPairFlag[idie]==1){
1196 Double_t cosine = GetOpeningAngle(e1[i1], e1[i2]);
1197 if(cosine<fdop){
1198 e1[i1]->SetConvFlag(0);
1199 e1[i2]->SetConvFlag(0);
1200 }
1201 }else if(fRejectPairFlag[idie]==2){
1202 Double_t phiv = GetPhiv(e1[i1], e1[i2]);
ebdc847a 1203 Double_t mee = GetMass(e1[i1], e1[i2]);
1204 if(fbz>0 && phiv>fdconvphiv && mee < fdconvMee){
805bd069 1205 e1[i1]->SetConvFlag(0);
1206 e1[i2]->SetConvFlag(0);
ebdc847a 1207 }else if(fbz<0 && phiv<acos(-1.0)-fdconvphiv && mee < fdconvMee){
805bd069 1208 e1[i1]->SetConvFlag(0);
1209 e1[i2]->SetConvFlag(0);
1210 }
1211 }
1212 }
1213 }
1214 }
1215
1216 if(e2.size()>0){
1217 for(int i1=0; i1<(int)e2.size(); i1++){
1218 for(int i2=i1+1; i2<(int)e2.size(); i2++){
1219 if(fRejectPairFlag[idie]==1){
1220 Double_t cosine = GetOpeningAngle(e2[i1], e2[i2]);
1221 if(cosine<fdop){
1222 e2[i1]->SetConvFlag(0);
1223 e2[i2]->SetConvFlag(0);
1224 }
1225 }else if(fRejectPairFlag[idie]==2){
1226 Double_t phiv = GetPhiv(e2[i1], e2[i2]);
ebdc847a 1227 Double_t mee = GetMass(e2[i1], e2[i2]);
1228 if(fbz>0 && phiv>fdconvphiv && mee < fdconvMee){
805bd069 1229 e2[i1]->SetConvFlag(0);
1230 e2[i2]->SetConvFlag(0);
ebdc847a 1231 }else if(fbz<0 && phiv<acos(-1.0)-fdconvphiv && mee < fdconvMee){
805bd069 1232 e2[i1]->SetConvFlag(0);
1233 e2[i2]->SetConvFlag(0);
1234 }
1235 }
1236 }
1237 }
1238 }
1239}
1240
1241
1242//_________________________________________________________________________________
1243Double_t AliAnalysisTaskMultiDielectronTG::GetOpeningAngle(AliDielectronSingleTG* e1, AliDielectronSingleTG* e2){
1244
1245 //////////////////////
1246 //////// calculate pairs and get opening angle
1247 //////////////////////
1248
1249 double dmass, dphiv, dpxpair, dpypair, dpzpair;
1250 double dptpair, depair, dphipair, detapair, dcos, dpsi;
1251
1252 CalcVars(e1, e2, dmass, dphiv, dpxpair, dpypair, dpzpair,
1253 dptpair, depair, dphipair, detapair, dcos, dpsi);
1254
1255 return dcos;
1256}
1257
1258//_________________________________________________________________________________
1259Double_t AliAnalysisTaskMultiDielectronTG::GetPhiv(AliDielectronSingleTG* e1, AliDielectronSingleTG* e2){
1260
1261 //////////////////////
ebdc847a 1262 //////// calculate pairs and get phiv
805bd069 1263 //////////////////////
1264
1265 double dmass, dphiv, dpxpair, dpypair, dpzpair;
1266 double dptpair, depair, dphipair, detapair, dcos, dpsi;
1267
1268 CalcVars(e1, e2, dmass, dphiv, dpxpair, dpypair, dpzpair,
1269 dptpair, depair, dphipair, detapair, dcos, dpsi);
1270
1271 return dphiv;
1272}
ebdc847a 1273
1274//_________________________________________________________________________________
1275Double_t AliAnalysisTaskMultiDielectronTG::GetMass(AliDielectronSingleTG* e1, AliDielectronSingleTG* e2){
1276
1277 //////////////////////
1278 //////// calculate pairs and get mass
1279 //////////////////////
1280
1281 double dmass, dphiv, dpxpair, dpypair, dpzpair;
1282 double dptpair, depair, dphipair, detapair, dcos, dpsi;
1283
1284 CalcVars(e1, e2, dmass, dphiv, dpxpair, dpypair, dpzpair,
1285 dptpair, depair, dphipair, detapair, dcos, dpsi);
1286
1287 return dmass;
1288}