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