1 /*************************************************************************
2 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 ///////////////////////////////////////////////////////////////////////////
18 // Basic Analysis Task //
20 ///////////////////////////////////////////////////////////////////////////
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>
35 #include "AliDielectron.h"
36 #include "AliDielectronHistos.h"
37 #include "AliDielectronCF.h"
38 #include "AliDielectronMC.h"
39 #include "AliDielectronMixingHandler.h"
40 #include "AliAnalysisTaskMultiDielectronTG.h"
42 #include "AliESDtrack.h"
43 #include "AliESDtrackCuts.h"
55 const char* kPairClassNamesTG[7] = {
66 ClassImp(AliDielectronSingleTG)
67 ClassImp(AliAnalysisTaskMultiDielectronTG)
69 //_________________________________________________________________________________
70 AliAnalysisTaskMultiDielectronTG::AliAnalysisTaskMultiDielectronTG() :
76 fSelectPhysics(kFALSE),
77 fTriggerMask(AliVEvent::kMB),
78 fExcludeTriggerMask(0),
79 fTriggerOnV0AND(kFALSE),
80 fRejectPileup(kFALSE),
82 fTriggerAnalysis(0x0),
88 fdEdXnSigmaElecvsPt(0x0),
90 fdEdXnSigmaElecvsPtTOF(0x0),
92 fTOFnSigmaElecvsPt(0x0),
93 fNCrossedRowsTPC(0x0),
95 fRatioCrossClusTPC(0x0),
100 fdconvphiv(acos(-1.0)),
109 //_________________________________________________________________________________
110 AliAnalysisTaskMultiDielectronTG::AliAnalysisTaskMultiDielectronTG(const char *name) :
111 AliAnalysisTaskSE(name),
116 fSelectPhysics(kFALSE),
117 fTriggerMask(AliVEvent::kMB),
118 fExcludeTriggerMask(0),
119 fTriggerOnV0AND(kFALSE),
120 fRejectPileup(kFALSE),
122 fTriggerAnalysis(0x0),
128 fdEdXnSigmaElecvsPt(0x0),
130 fdEdXnSigmaElecvsPtTOF(0x0),
132 fTOFnSigmaElecvsPt(0x0),
133 fNCrossedRowsTPC(0x0),
135 fRatioCrossClusTPC(0x0),
140 fdconvphiv(acos(-1.0)),
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();
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();
178 //_________________________________________________________________________________
179 AliAnalysisTaskMultiDielectronTG::~AliAnalysisTaskMultiDielectronTG()
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);
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();
207 //_________________________________________________________________________________
208 void AliAnalysisTaskMultiDielectronTG::UserCreateOutputObjects()
211 // Add all histogram manager histogram lists to the output TList
214 if (!fListHistos.IsEmpty()||!fListCF.IsEmpty()) return; //already initialised
216 // AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
217 // Bool_t isESD=man->GetInputEventHandler()->IsA()==AliESDInputHandler::Class();
218 // Bool_t isAOD=man->GetInputEventHandler()->IsA()==AliAODInputHandler::Class();
220 TIter nextDie(&fListDielectron);
221 AliDielectron *die=0;
222 while ( (die=static_cast<AliDielectron*>(nextDie())) ){
224 if (die->GetHistogramList()) fListHistos.Add(const_cast<THashList*>(die->GetHistogramList()));
225 if (die->GetCFManagerPair()) fListCF.Add(const_cast<AliCFContainer*>(die->GetCFManagerPair()->GetContainer()));
228 Int_t cuts=fListDielectron.GetEntries();
229 Int_t nbins=kNbinsEvent+2*cuts;
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.");
236 fEventStat->GetXaxis()->SetBinLabel(3,"Bin3 not used");
237 fEventStat->GetXaxis()->SetBinLabel(4,"Bin4 not used");
238 fEventStat->GetXaxis()->SetBinLabel(5,"Bin5 not used");
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");
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()));
250 if (!fTriggerAnalysis) fTriggerAnalysis=new AliTriggerAnalysis;
251 fTriggerAnalysis->EnableHistograms();
252 fTriggerAnalysis->SetAnalyzeMC(AliDielectronMC::Instance()->HasMC());
259 float binw = (TMath::Log(maxx)-TMath::Log(minx))/nbinx;
261 for(int ii=0;ii<nbinx+1;ii++){
262 xbin[ii] = TMath::Exp(TMath::Log(minx) + 0.5*binw+binw*ii);
266 fQAElectron = new TList();
267 fQAElectron->SetName("QAElectron");
268 fQAElectron->SetOwner();
271 fEvent = new TH1D("Event","centrality", 100,0,100);
272 fQAElectron->Add(fEvent);
273 fdEdXvsPt = new TH2D("dEdXvsPt","dE/dX vs. PT of TPC", nbinx, xbin, 2000,0,200);
274 fQAElectron->Add(fdEdXvsPt);
275 fdEdXnSigmaElecvsPt = new TH2D("fdEdXnSigmaElecvsPt"," dE/dX normalized to electron vs. pT of TPC",
276 nbinx, xbin, 2000, -10, 10);
277 fQAElectron->Add(fdEdXnSigmaElecvsPt);
279 fdEdXvsPtTOF = new TH2D("dEdXvsPtTOF","dE/dX vs. PT of TPC", nbinx, xbin, 2000,0,200);
280 fQAElectron->Add(fdEdXvsPtTOF);
281 fdEdXnSigmaElecvsPtTOF = new TH2D("fdEdXnSigmaElecvsPtTOF"," dE/dX normalized to electron vs. pT of TPC",
282 nbinx, xbin, 2000, -10, 10);
283 fQAElectron->Add(fdEdXnSigmaElecvsPtTOF);
287 fTOFbetavsPt = new TH2D("fTOFbetavsPt","TOF beta vs. p", 400, 0, 20, 1200, 0, 1.2);
288 fQAElectron->Add(fTOFbetavsPt);
289 fTOFnSigmaElecvsPt = new TH2D("fTOFnSigmaElecvsPt","TOF nsigma for electron", 400, 0, 20, 2000, -10, 10);
290 fQAElectron->Add(fTOFnSigmaElecvsPt);
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);
297 fRatioCrossClusTPC = new TH2F("fRatioCrossClusTPC", "hRatioCrossClusTPC vs. pT", 200, 0, 20, 200, 0, 10);
298 fQAElectron->Add(fRatioCrossClusTPC);
300 fListHistos.Add(fQAElectron);
302 fListHistos.SetOwner();
304 PostData(1, &fListHistos);
305 PostData(2, &fListCF);
306 PostData(3, fEventStat);
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);
323 //_________________________________________________________________________________
324 void AliAnalysisTaskMultiDielectronTG::UserExec(Option_t *)
327 // Main loop. Called for every event
330 if (fListHistos.IsEmpty()&&fListCF.IsEmpty()) return;
332 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
333 Bool_t isESD=man->GetInputEventHandler()->IsA()==AliESDInputHandler::Class();
334 Bool_t isAOD=man->GetInputEventHandler()->IsA()==AliAODInputHandler::Class();
336 AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
337 if (!inputHandler) return;
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() );
345 AliFatal("This task needs the PID response attached to the input event handler!");
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);
361 //Before physics selection
362 fEventStat->Fill(kAllEvents);
363 if (isSelected==0||isRejected) {
364 PostData(3,fEventStat);
367 //after physics selection
368 fEventStat->Fill(kSelectedEvents);
372 if(isESD){if (!fTriggerAnalysis->IsOfflineTriggerFired(static_cast<AliESDEvent*>(InputEvent()), AliTriggerAnalysis::kV0AND))
374 if(isAOD){if(!((static_cast<AliAODEvent*>(InputEvent()))->GetVZEROData()->GetV0ADecision() == AliVVZERO::kV0BB &&
375 (static_cast<AliAODEvent*>(InputEvent()))->GetVZEROData()->GetV0CDecision() == AliVVZERO::kV0BB) )
380 fEventStat->Fill(kV0andEvents);
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();
392 if (AliDielectronMC::Instance()->ConnectMCEvent())
393 AliDielectronVarManager::Fill(AliDielectronMC::Instance()->GetMCEvent(),valuesMC);
396 while ( (die=static_cast<AliDielectron*>(nextDie())) ){
397 AliDielectronHistos *h=die->GetHistoManager();
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);
409 if (!fEventFilter->IsSelected(InputEvent())) return;
411 fEventStat->Fill(kFilteredEvents);
415 if (InputEvent()->IsPileupFromSPD(3,0.8,3.,2.,5.)) return;
417 fEventStat->Fill(kPileupEvents);
420 fbz = InputEvent()->GetMagneticField();
421 AliKFParticle::SetField( fbz );
422 AliDielectronPID::SetCorrVal((Double_t)InputEvent()->GetRunNumber());
424 //Process event in all AliDielectron instances
425 // TIter nextDie(&fListDielectron);
426 // AliDielectron *die=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);
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){
446 ////////////////////////////////////////////////////////////////////
447 AliDielectronVarManager::Fill(obj->UncheckedAt(itrack), fgValues);
448 ////////////////////////////////////////////////////////////////////
450 if(fgValues[AliDielectronVarManager::kCharge]>0){
451 fVeptmp.push_back(new AliDielectronSingleTG(1,
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)))
465 }else if(fgValues[AliDielectronVarManager::kCharge]<0){
466 fVemtmp.push_back(new AliDielectronSingleTG(-1,
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)))
483 //AliInfo(Form("size of e and p = %d %d", (int)fVeptmp.size(), (int)fVemtmp.size()));
486 CheckGhostPairs(fVeptmp);
487 CheckGhostPairs(fVemtmp);
488 RandomizePool(fVeptmp, fVemtmp);
489 CalcPair(fVep, fVem, die, idie);
491 // AliInfo(Form("size of e and p (after) = %d %d", (int)fVep.size(), (int)fVem.size()));
493 double dwcent = 100.0/fgkNCENT;
494 double dwiz = 20.0/fgkNZBIN;
495 double dwrp = acos(-1.0)/fgkNRPBIN;
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);
502 if(icent>=fgkNCENT) icent=fgkNCENT-1;
504 if(izbin>=fgkNZBIN) izbin=fgkNZBIN-1;
506 if(irp>=fgkNRPBIN) irp=fgkNRPBIN-1;
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();
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();
526 fibuf[idie][izbin][icent][irp]++;
527 if(fibuf[idie][izbin][icent][irp]>= fgkNBUF) fibuf[idie][izbin][icent][irp]=0;
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);
545 Printf("ERROR: Could not receive track %d", iTracks);
548 if(!fCutsMother->AcceptTrack(track)) continue;
549 fdEdXvsPt->Fill(track->GetTPCmomentum(), track->GetTPCsignal());
550 fdEdXnSigmaElecvsPt->Fill(track->GetTPCmomentum(),
551 AliDielectronVarManager::GetPIDResponse()->NumberOfSigmasTPC(track,
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
559 if( (l < 360. || l > 800.) || (t <= 0.) || (t0 >999990.0) ) {
563 t -= t0; // subtract the T0
565 t *= 1e-12; //ps -> s
568 beta = v / TMath::C();
571 fTOFbetavsPt->Fill(track->GetTPCmomentum(), beta);
572 fTOFnSigmaElecvsPt->Fill(track->GetTPCmomentum(),
573 AliDielectronVarManager::GetPIDResponse()->NumberOfSigmasTOF(track,
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,
581 -AliDielectronPID::GetCorrVal());
584 if(track->GetTPCsignal()>70 && track->GetTPCsignal()<90){
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());
593 PostData(1, &fListHistos);
594 PostData(2, &fListCF);
595 PostData(3,fEventStat);
598 //_________________________________________________________________________________
599 void AliAnalysisTaskMultiDielectronTG::FinishTaskOutput()
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);
612 PostData(1, &fListHistos);
613 PostData(2, &fListCF);
616 //_________________________________________________________________________________
617 void AliAnalysisTaskMultiDielectronTG::CheckGhostPairs(vector<AliDielectronSingleTG*> e1)
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
625 for(int i1=0; i1<(int)e1.size(); i1++){
627 for(int i2=i1+1; i2<(int)e1.size(); i2++){
628 if( fabs(e1[i1]->Phi() - e1[i2]->Phi())<0.01 ){
630 e1[i2]->SetGstFlag(0);
633 if(reject==true)e1[i1]->SetGstFlag(0);
638 //_________________________________________________________________________________
639 void AliAnalysisTaskMultiDielectronTG::RandomizePool(vector<AliDielectronSingleTG*> e1, vector<AliDielectronSingleTG*> e2)
642 ///// Randomize the pool constent to cancel the filling scheme of the single tracks
646 int size1 = e1.size();
648 for(int i=0;i<1000;i++){
651 for(int i=0;i<size1;i++){
655 for(int i=0;i<size1;i++){
656 int j = (int)(gRandom->Uniform(0,size1));
657 while(usedindex[j]==1){
658 j = (int)(gRandom->Uniform(0,size1));
660 if( (e1[j]->GetGstFlag()==1) &&
661 (e1[j]->GetConvFlag()==1)
663 fVep.push_back(e1[j]);
669 int size2 = e2.size();
670 for(int i=0;i<1000;i++){
673 for(int i=0;i<size2;i++){
677 for(int i=0;i<size2;i++){
678 int j = (int)(gRandom->Uniform(0,size2));
679 while(usedindex[j]==1){
680 j = (int)(gRandom->Uniform(0,size2));
682 if( (e2[j]->GetGstFlag()==1) &&
683 (e2[j]->GetConvFlag()==1)
685 fVem.push_back(e2[j]);
692 //_________________________________________________________________________________
693 void AliAnalysisTaskMultiDielectronTG::CalcPair(vector<AliDielectronSingleTG*> ve1,
694 vector<AliDielectronSingleTG*> ve2, AliDielectron *die, Int_t idie)
698 // main routine for the pair procesing
702 for(int i1=0; i1<(int)ve1.size(); i1++){
703 for(int i2=0; i2<(int)ve2.size(); i2++){
704 FillPair(ve1[i1], ve2[i2], 0, die);
708 for(int i1=0; i1<(int)ve1.size(); i1++){
709 for(int i2=i1+1; i2<(int)ve1.size(); i2++){
710 FillPair(ve1[i1], ve1[i2], 1, die);
714 for(int i1=0; i1<(int)ve2.size(); i1++){
715 for(int i2=i1+1; i2<(int)ve2.size(); i2++){
716 FillPair(ve2[i1], ve2[i2], 2, die);
721 double dwcent = 100.0/fgkNCENT;
722 double dwiz = 20.0/fgkNZBIN;
723 double dwrp = acos(-1.0)/fgkNRPBIN;
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);
730 if(icent>=fgkNCENT) icent=fgkNCENT-1;
732 if(izbin>=fgkNZBIN) izbin=fgkNZBIN-1;
734 if(irp>=fgkNRPBIN) irp=fgkNRPBIN-1;
739 // Now mixed event for +- pairs
742 for(int ibuf=0;(nmixed<fgkNMix);ibuf++) {
744 while(ntry<fgkMAXTRY) {
745 ReshuffleBuffer(fvem[ibuf][idie][izbin][icent][irp],fpoolm[idie][izbin][icent][irp]);
748 for(int i1=0; i1<(int)ve1.size(); i1++){
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);
758 // Now mixed event for +- pairs
761 for(int ibuf=0;(nmixed<fgkNMix);ibuf++) {
763 while(ntry<fgkMAXTRY) {
764 ReshuffleBuffer(fvep[ibuf][idie][izbin][icent][irp],fpoolp[idie][izbin][icent][irp]);
767 for(int i1=0; i1<(int)ve2.size(); i1++){
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);
778 // Now mixed event for ++ pairs
781 for(int ibuf=0;(nmixed<fgkNMix);ibuf++) {
783 while(ntry<fgkMAXTRY) {
784 ReshuffleBuffer(fvep[ibuf][idie][izbin][icent][irp],fpoolp[idie][izbin][icent][irp]);
787 for(int i1=0; i1<(int)ve1.size(); i1++){
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);
798 // Now mixed event for +- pairs
801 for(int ibuf=0;(nmixed<fgkNMix);ibuf++) {
803 while(ntry<fgkMAXTRY) {
804 ReshuffleBuffer(fvem[ibuf][idie][izbin][icent][irp],fpoolm[idie][izbin][icent][irp]);
807 for(int i1=0; i1<(int)ve2.size(); i1++){
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);
819 //_________________________________________________________________________________
820 void AliAnalysisTaskMultiDielectronTG::FillPair(AliDielectronSingleTG *iep,
821 AliDielectronSingleTG *iem, int type, AliDielectron *die)
825 // main routine for filling kinematics of pairs
830 double dmass, dphiv, dpxpair, dpypair, dpzpair;
831 double dptpair, depair, dphipair, detapair, dcos, dpsi;
833 CalcVars(iep, iem, dmass, dphiv, dpxpair, dpypair, dpzpair,
834 dptpair, depair, dphipair, detapair, dcos, dpsi);
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;
852 ////////////////////////////
853 ///// calculate v0 ////////
854 ///////////////////////////
856 /// for the moment, this doesn't work for MC
857 if(fdv0mixing == kFALSE && (type==3 || type==4 || type==5 || type==6)){
860 if(die->GetHasMC()==kTRUE && (type==3 || type==4 || type==5 || type==6)){
863 if(type==0 || type==1 || type==2){
869 AliVTrack* trackob1= iep->GetTrack();
870 AliVTrack* trackob2= iem->GetTrack();
872 if(!trackob1 || !trackob2){
876 AliDielectronPair candidate;
877 candidate.SetTracks(trackob1, (int)(11*iep->Charge()),
878 trackob2, (int)(11*iem->Charge()));
880 if(trackob1==trackob2){
881 AliInfo("Same Objects!!");
884 const AliKFParticle &kfPair = candidate.GetKFParticle();
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();
900 dpsipair = candidate.PsiPair(fbz);
903 Double_t values[AliDielectronVarManager::kNMaxValues];
906 className1.Form("MyPair_%s",kPairClassNamesTG[type]);
907 className2.Form("MyPairV0_%s",kPairClassNamesTG[type]);
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;
913 if (pairClass1 && PairTrackcut(dphiv)==true){
914 ///import pair variables to values!!
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;
929 fHistos->FillClass(className1, AliDielectronVarManager::kNMaxValues, values);
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;
948 fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
955 //_________________________________________________________________________________
956 bool AliAnalysisTaskMultiDielectronTG::PairTrackcut(double phiv)
966 //var is phiv for the moment
967 if(fbz>0 && phiv>fdconvphiv){
969 }else if(fbz<0 && phiv<acos(-1.0)-fdconvphiv){
978 //_________________________________________________________________________________
979 void AliAnalysisTaskMultiDielectronTG::CalcVars(AliDielectronSingleTG *iep, AliDielectronSingleTG *iem,
980 double &mass, double &phiv, double &px, double &py, double&pz,
981 double &pt, double &e, double &phi,
982 double &eta, double &cos, double &psi)
987 // standalone calculator for the pair variables
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);
994 double dppair = sqrt(pt*pt+pz*pz);
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());
999 mass = e*e-px*px-py*py-pz*pz;
1007 phi = atan2(py, px);
1008 eta = -0.5*TMath::Log((dppair+pz)/(dppair-pz));
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));
1013 double dtheta = iep->Theta()-iem->Theta();
1014 psi = asin(dtheta/cos);
1017 //unit vector of (pep+pem)
1022 double ax = uy/sqrt(ux*ux+uy*uy);
1023 double ay = -ux/sqrt(ux*ux+uy*uy);
1025 //momentum of e+ and e- in (ax,ay,az) axis. Note that az=0 by
1027 //double ptep = iep->Px()*ax + iep->Py()*ay;
1028 //double ptem = iem->Px()*ax + iem->Py()*ay;
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();
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);
1045 //unit vector of pep X pem
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
1060 double cosPhiV = wx*ax + wy*ay;
1061 phiv = acos(cosPhiV);
1065 //_________________________________________________________________________________
1066 void AliAnalysisTaskMultiDielectronTG::ReshuffleBuffer(vector<AliDielectronSingleTG*> ve, deque<AliDielectronSingleTG*> pool)
1068 //If there is not enough electron in the pool, give up
1070 // ReshuffleBuffer for th event mixing
1073 unsigned int ne = ve.size();
1074 unsigned int poolsize = pool.size();
1075 int used[fgkMAXPOOL];
1076 for(int i=0;i<(int)fgkMAXPOOL;i++){
1081 std::cout <<" pool size="<<poolsize<<" ne"<<ne<<std::endl;
1084 for(unsigned int ie=0; ie < ne; ie++) {
1085 int j = rand()%poolsize;
1087 j = rand()%poolsize;