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