]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGDQ/dielectron/AliDielectron.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Corrected EINCLUDE, compilation with Root6
[u/mrichter/AliRoot.git] / PWGDQ / dielectron / AliDielectron.cxx
CommitLineData
b2a297fa 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// Dielectron Analysis Main class //
18// //
19/*
20Framework to perform event selectoin, single track selection and track pair
21selection.
22
23Convention for the signs of the pair in fPairCandidates:
24The names are available via the function PairClassName(Int_t i)
25
260: ev1+ ev1+ (same event like sign +)
271: ev1+ ev1- (same event unlike sign)
282: ev1- ev1- (same event like sign -)
29
303: ev1+ ev2+ (mixed event like sign +)
314: ev1- ev2+ (mixed event unlike sign -+)
326: ev1+ ev2- (mixed event unlike sign +-)
337: ev1- ev2- (mixed event like sign -)
34
355: ev2+ ev2+ (same event like sign +)
368: ev2+ ev2- (same event unlike sign)
379: ev2- ev2- (same event like sign -)
38
2a14a7b1 3910: ev1+ ev1- (same event track rotation)
b2a297fa 40
41*/
42// //
43///////////////////////////////////////////////////////////////////////////
44
45#include <TString.h>
46#include <TList.h>
47#include <TMath.h>
5720c765 48#include <TObject.h>
b2a297fa 49
5720c765 50#include <AliKFParticle.h>
b2a297fa 51
4533e78e 52#include <AliESDInputHandler.h>
53#include <AliAnalysisManager.h>
54#include <AliEPSelectionTask.h>
5720c765 55#include <AliEventplane.h>
b2a297fa 56#include <AliVEvent.h>
57#include <AliVParticle.h>
58#include <AliVTrack.h>
59#include "AliDielectronPair.h"
60#include "AliDielectronHistos.h"
61#include "AliDielectronCF.h"
62#include "AliDielectronMC.h"
63#include "AliDielectronVarManager.h"
2a14a7b1 64#include "AliDielectronTrackRotator.h"
572b0139 65#include "AliDielectronDebugTree.h"
ba15fdfb 66#include "AliDielectronSignalMC.h"
5720c765 67#include "AliDielectronMixingHandler.h"
061ca303 68#include "AliDielectronV0Cuts.h"
572b0139 69
b2a297fa 70#include "AliDielectron.h"
71
72ClassImp(AliDielectron)
73
74const char* AliDielectron::fgkTrackClassNames[4] = {
75 "ev1+",
76 "ev1-",
77 "ev2+",
78 "ev2-"
79};
80
2a14a7b1 81const char* AliDielectron::fgkPairClassNames[11] = {
b2a297fa 82 "ev1+_ev1+",
83 "ev1+_ev1-",
84 "ev1-_ev1-",
85 "ev1+_ev2+",
86 "ev1-_ev2+",
87 "ev2+_ev2+",
88 "ev1+_ev2-",
89 "ev1-_ev2-",
90 "ev2+_ev2-",
2a14a7b1 91 "ev2-_ev2-",
92 "ev1+_ev1-_TR"
b2a297fa 93};
94
95//________________________________________________________________
96AliDielectron::AliDielectron() :
97 TNamed("AliDielectron","AliDielectron"),
d8bb1abe 98 fCutQA(kFALSE),
99 fQAmonitor(0x0),
b2a297fa 100 fEventFilter("EventFilter"),
101 fTrackFilter("TrackFilter"),
61d106d3 102 fPairPreFilter("PairPreFilter"),
103 fPairPreFilterLegs("PairPreFilterLegs"),
b2a297fa 104 fPairFilter("PairFilter"),
5720c765 105 fEventPlanePreFilter("EventPlanePreFilter"),
106 fEventPlanePOIPreFilter("EventPlanePOIPreFilter"),
b2a297fa 107 fPdgMother(443),
8df8e382 108 fPdgLeg1(11),
109 fPdgLeg2(11),
ba15fdfb 110 fSignalsMC(0x0),
554e40f8 111 fNoPairing(kFALSE),
08b801a6 112 fUseKF(kTRUE),
620b34a5 113 fHistoArray(0x0),
b2a297fa 114 fHistos(0x0),
5720c765 115 fPairCandidates(new TObjArray(11)),
572b0139 116 fCfManagerPair(0x0),
2a14a7b1 117 fTrackRotator(0x0),
554e40f8 118 fDebugTree(0x0),
5720c765 119 fMixing(0x0),
120 fPreFilterEventPlane(kFALSE),
121 fLikeSignSubEvents(kFALSE),
fb7d2d99 122 fPreFilterUnlikeOnly(kFALSE),
ba15fdfb 123 fPreFilterAllSigns(kFALSE),
5720c765 124 fHasMC(kFALSE),
125 fStoreRotatedPairs(kFALSE),
2e02dba4 126 fDontClearArrays(kFALSE),
5720c765 127 fEstimatorFilename(""),
a823f01b 128 fTRDpidCorrectionFilename(""),
129 fVZEROCalibrationFilename(""),
130 fVZERORecenteringFilename("")
b2a297fa 131{
132 //
133 // Default constructor
134 //
135
136}
137
138//________________________________________________________________
139AliDielectron::AliDielectron(const char* name, const char* title) :
140 TNamed(name,title),
d8bb1abe 141 fCutQA(kFALSE),
142 fQAmonitor(0x0),
b2a297fa 143 fEventFilter("EventFilter"),
144 fTrackFilter("TrackFilter"),
61d106d3 145 fPairPreFilter("PairPreFilter"),
146 fPairPreFilterLegs("PairPreFilterLegs"),
b2a297fa 147 fPairFilter("PairFilter"),
5720c765 148 fEventPlanePreFilter("EventPlanePreFilter"),
149 fEventPlanePOIPreFilter("EventPlanePOIPreFilter"),
b2a297fa 150 fPdgMother(443),
8df8e382 151 fPdgLeg1(11),
152 fPdgLeg2(11),
ba15fdfb 153 fSignalsMC(0x0),
554e40f8 154 fNoPairing(kFALSE),
08b801a6 155 fUseKF(kTRUE),
620b34a5 156 fHistoArray(0x0),
b2a297fa 157 fHistos(0x0),
5720c765 158 fPairCandidates(new TObjArray(11)),
572b0139 159 fCfManagerPair(0x0),
2a14a7b1 160 fTrackRotator(0x0),
554e40f8 161 fDebugTree(0x0),
5720c765 162 fMixing(0x0),
163 fPreFilterEventPlane(kFALSE),
164 fLikeSignSubEvents(kFALSE),
fb7d2d99 165 fPreFilterUnlikeOnly(kFALSE),
ba15fdfb 166 fPreFilterAllSigns(kFALSE),
5720c765 167 fHasMC(kFALSE),
168 fStoreRotatedPairs(kFALSE),
2e02dba4 169 fDontClearArrays(kFALSE),
5720c765 170 fEstimatorFilename(""),
a823f01b 171 fTRDpidCorrectionFilename(""),
172 fVZEROCalibrationFilename(""),
173 fVZERORecenteringFilename("")
b2a297fa 174{
175 //
176 // Named constructor
177 //
178
179}
180
181//________________________________________________________________
182AliDielectron::~AliDielectron()
183{
184 //
185 // Default destructor
186 //
d8bb1abe 187 if (fQAmonitor) delete fQAmonitor;
620b34a5 188 if (fHistoArray) delete fHistoArray;
b2a297fa 189 if (fHistos) delete fHistos;
190 if (fPairCandidates) delete fPairCandidates;
572b0139 191 if (fDebugTree) delete fDebugTree;
5720c765 192 if (fMixing) delete fMixing;
ba15fdfb 193 if (fSignalsMC) delete fSignalsMC;
5720c765 194 if (fCfManagerPair) delete fCfManagerPair;
b2a297fa 195}
196
197//________________________________________________________________
198void AliDielectron::Init()
199{
200 //
201 // Initialise objects
202 //
fb7d2d99 203
204 if(GetHasMC()) AliDielectronMC::Instance()->SetHasMC(GetHasMC());
5720c765 205
206 InitPairCandidateArrays();
207
ba15fdfb 208 if (fCfManagerPair) {
209 fCfManagerPair->SetSignalsMC(fSignalsMC);
210 fCfManagerPair->InitialiseContainer(fPairFilter);
211 }
1201a1a9 212 if (fTrackRotator) {
213 fTrackRotator->SetTrackArrays(&fTracks[0],&fTracks[1]);
214 fTrackRotator->SetPdgLegs(fPdgLeg1,fPdgLeg2);
215 }
8df8e382 216 if (fDebugTree) fDebugTree->SetDielectron(this);
5720c765 217 if(fEstimatorFilename.Contains(".root")) AliDielectronVarManager::InitEstimatorAvg(fEstimatorFilename.Data());
218 if(fTRDpidCorrectionFilename.Contains(".root")) AliDielectronVarManager::InitTRDpidEffHistograms(fTRDpidCorrectionFilename.Data());
a823f01b 219 if(fVZEROCalibrationFilename.Contains(".root")) AliDielectronVarManager::SetVZEROCalibrationFile(fVZEROCalibrationFilename.Data());
220 if(fVZERORecenteringFilename.Contains(".root")) AliDielectronVarManager::SetVZERORecenteringFile(fVZERORecenteringFilename.Data());
221
5720c765 222 if (fMixing) fMixing->Init(this);
443a091c 223 if (fHistoArray) {
224 fHistoArray->SetSignalsMC(fSignalsMC);
225 fHistoArray->Init();
226 }
d8bb1abe 227
228 if (fCutQA) {
229 fQAmonitor = new AliDielectronCutQA(Form("QAcuts_%s",GetName()),"QAcuts");
230 fQAmonitor->AddTrackFilter(&fTrackFilter);
231 fQAmonitor->AddPairFilter(&fPairFilter);
232 fQAmonitor->AddEventFilter(&fEventFilter);
233 fQAmonitor->Init();
234 }
235}
b2a297fa 236
237//________________________________________________________________
238void AliDielectron::Process(AliVEvent *ev1, AliVEvent *ev2)
239{
240 //
241 // Process the events
242 //
243
45b2b1b8 244 //at least first event is needed!
245 if (!ev1){
246 AliError("At least first event must be set!");
247 return;
248 }
5720c765 249
061ca303 250 // modify event numbers in MC so that we can identify new events
251 // in AliDielectronV0Cuts (not neeeded for collision data)
252 if(GetHasMC()) {
253 ev1->SetBunchCrossNumber(1);
254 ev1->SetOrbitNumber(1);
255 ev1->SetPeriodNumber(1);
256 }
257
c315310a 258 // set event
259 AliDielectronVarManager::SetEvent(ev1);
260 if (fMixing){
261 //set mixing bin to event data
262 Int_t bin=fMixing->FindBin(AliDielectronVarManager::GetData());
263 AliDielectronVarManager::SetValue(AliDielectronVarManager::kMixingBin,bin);
264 }
265
b2a297fa 266 //in case we have MC load the MC event and process the MC particles
d8bb1abe 267 // why do not apply the event cuts first ????
40875e45 268 if (AliDielectronMC::Instance()->ConnectMCEvent()){
5720c765 269 ProcessMC(ev1);
572b0139 270 }
0caf5fbb 271
b2a297fa 272 //if candidate array doesn't exist, create it
273 if (!fPairCandidates->UncheckedAt(0)) {
274 InitPairCandidateArrays();
275 } else {
276 ClearArrays();
277 }
278
279 //mask used to require that all cuts are fulfilled
280 UInt_t selectedMask=(1<<fEventFilter.GetCuts()->GetEntries())-1;
281
282 //apply event cuts
d8bb1abe 283 UInt_t cutmask = fEventFilter.IsSelected(ev1);
284 if(fCutQA) fQAmonitor->FillAll(ev1);
285 if(fCutQA) fQAmonitor->Fill(cutmask,ev1);
286 if ((ev1&&cutmask!=selectedMask) ||
287 (ev2&&fEventFilter.IsSelected(ev2)!=selectedMask)) return;
5720c765 288
b2a297fa 289 //fill track arrays for the first event
554e40f8 290 if (ev1){
291 FillTrackArrays(ev1);
ba15fdfb 292 if (((fPreFilterAllSigns)||(fPreFilterUnlikeOnly)) && ( fPairPreFilter.GetCuts()->GetEntries()>0 )) PairPreFilter(0, 1, fTracks[0], fTracks[1]);
554e40f8 293 }
294
b2a297fa 295
296 //fill track arrays for the second event
554e40f8 297 if (ev2) {
298 FillTrackArrays(ev2,1);
ba15fdfb 299 if (((fPreFilterAllSigns)||(fPreFilterUnlikeOnly)) && ( fPairPreFilter.GetCuts()->GetEntries()>0 )) PairPreFilter(2, 3, fTracks[2], fTracks[3]);
554e40f8 300 }
b2a297fa 301
5720c765 302 // TPC event plane correction
4533e78e 303 if (ev1 && fPreFilterEventPlane && ( fEventPlanePreFilter.GetCuts()->GetEntries()>0 || fEventPlanePOIPreFilter.GetCuts()->GetEntries()>0))
304 EventPlanePreFilter(0, 1, fTracks[0], fTracks[1], ev1);
f0351f9c 305
554e40f8 306 if (!fNoPairing){
307 // create pairs and fill pair candidate arrays
308 for (Int_t itrackArr1=0; itrackArr1<4; ++itrackArr1){
309 for (Int_t itrackArr2=itrackArr1; itrackArr2<4; ++itrackArr2){
310 FillPairArrays(itrackArr1, itrackArr2);
311 }
b2a297fa 312 }
6551594b 313
554e40f8 314 //track rotation
1201a1a9 315 if (fTrackRotator) {
316 fTrackRotator->SetEvent(ev1);
317 FillPairArrayTR();
318 }
554e40f8 319 }
572b0139 320
321 //fill debug tree if a manager is attached
322 if (fDebugTree) FillDebugTree();
5720c765 323
324 //process event mixing
325 if (fMixing) {
326 fMixing->Fill(ev1,this);
4533e78e 327 // FillHistograms(0x0,kTRUE);
5720c765 328 }
329
330 //in case there is a histogram manager, fill the QA histograms
93d4ec8a 331 if (fHistos && fSignalsMC) FillMCHistograms(ev1);
5720c765 332 if (fHistos) FillHistograms(ev1);
333
f0351f9c 334
5720c765 335 // clear arrays
2e02dba4 336 if (!fDontClearArrays) ClearArrays();
061ca303 337
338 // reset TPC EP and unique identifiers for v0 cut class
5720c765 339 AliDielectronVarManager::SetTPCEventPlane(0x0);
061ca303 340 if(GetHasMC()) { // only for MC needed
341 for (Int_t iCut=0; iCut<fTrackFilter.GetCuts()->GetEntries();++iCut) {
342 if ( fTrackFilter.GetCuts()->At(iCut)->IsA() == AliDielectronV0Cuts::Class() )
343 ((AliDielectronV0Cuts*)fTrackFilter.GetCuts()->At(iCut))->ResetUniqueEventNumbers();
344 }
345 }
346
b2a297fa 347}
348
349//________________________________________________________________
5720c765 350void AliDielectron::ProcessMC(AliVEvent *ev1)
b2a297fa 351{
352 //
353 // Process the MC data
354 //
355
ba15fdfb 356 AliDielectronMC *dieMC=AliDielectronMC::Instance();
357
5720c765 358 if (fHistos) FillHistogramsMC(dieMC->GetMCEvent(), ev1);
ba15fdfb 359
360 if(!fSignalsMC) return;
bc3cacd4 361 //loop over all MC data and Fill the HF, CF containers and histograms if they exist
88204efa 362 if(fCfManagerPair) fCfManagerPair->SetPdgMother(fPdgMother);
363
ba15fdfb 364 // signals to be studied
365 Int_t nSignals = fSignalsMC->GetEntries();
366
833da63d 367 Bool_t bFillCF = (fCfManagerPair ? fCfManagerPair->GetStepForMCtruth() : kFALSE);
368 Bool_t bFillHF = (fHistoArray ? fHistoArray->GetStepForMCGenerated() : kFALSE);
369 Bool_t bFillHist = kFALSE;
370 for(Int_t isig=0;isig<nSignals;isig++) {
371 TString sigName = fSignalsMC->At(isig)->GetName();
372 const THashList *histlist = fHistos->GetHistogramList();
373 bFillHist |= histlist->FindObject(Form("Pair_%s_MCtruth",sigName.Data()))!=0x0;
374 bFillHist |= histlist->FindObject(Form("Track_Leg_%s_MCtruth",sigName.Data()))!=0x0;
375 bFillHist |= histlist->FindObject(Form("Track_%s_%s_MCtruth",fgkPairClassNames[1],sigName.Data()))!=0x0;
376 if(bFillHist) break;
377 }
378 // check if there is anything to fill
379 if(!bFillCF && !bFillHF && !bFillHist) return;
380
381
ba15fdfb 382 // initialize 2D arrays of labels for particles from each MC signal
383 Int_t** labels1; // labels for particles satisfying branch 1
384 Int_t** labels2; // labels for particles satisfying branch 2
385 Int_t** labels12; // labels for particles satisfying both branches
386 labels1 = new Int_t*[nSignals];
387 labels2 = new Int_t*[nSignals];
388 labels12 = new Int_t*[nSignals];
389 Int_t* indexes1=new Int_t[nSignals];
390 Int_t* indexes2=new Int_t[nSignals];
391 Int_t* indexes12=new Int_t[nSignals];
392 for(Int_t isig=0;isig<nSignals;++isig) {
393 *(labels1+isig) = new Int_t[dieMC->GetNMCTracks()];
394 *(labels2+isig) = new Int_t[dieMC->GetNMCTracks()];
395 *(labels12+isig) = new Int_t[dieMC->GetNMCTracks()];
396 for(Int_t ip=0; ip<dieMC->GetNMCTracks();++ip) {
397 labels1[isig][ip] = -1;
398 labels2[isig][ip] = -1;
399 labels12[isig][ip] = -1;
400 }
401 indexes1[isig]=0;
402 indexes2[isig]=0;
403 indexes12[isig]=0;
b2a297fa 404 }
ba15fdfb 405
406 Bool_t truth1=kFALSE;
407 Bool_t truth2=kFALSE;
408 // loop over the MC tracks
409 for(Int_t ipart=0; ipart<dieMC->GetNMCTracks(); ++ipart) {
410 for(Int_t isig=0; isig<nSignals; ++isig) { // loop over signals
411 // Proceed only if this signal is required in the pure MC step
412 // NOTE: Some signals can be satisfied by many particles and this leads to high
413 // computation times (e.g. secondary electrons from the GEANT transport). Be aware of this!!
414 if(!((AliDielectronSignalMC*)fSignalsMC->At(isig))->GetFillPureMCStep()) continue;
415
416 truth1 = dieMC->IsMCTruth(ipart, (AliDielectronSignalMC*)fSignalsMC->At(isig), 1);
417 truth2 = dieMC->IsMCTruth(ipart, (AliDielectronSignalMC*)fSignalsMC->At(isig), 2);
418
419 // particles satisfying both branches are treated separately to avoid double counting during pairing
420 if(truth1 && truth2) {
421 labels12[isig][indexes12[isig]] = ipart;
422 ++indexes12[isig];
423 }
424 else {
425 if(truth1) {
426 labels1[isig][indexes1[isig]] = ipart;
427 ++indexes1[isig];
428 }
429 if(truth2) {
430 labels2[isig][indexes2[isig]] = ipart;
431 ++indexes2[isig];
432 }
433 }
434 }
435 } // end loop over MC particles
436
437 // Do the pairing and fill the CF container with pure MC info
438 for(Int_t isig=0; isig<nSignals; ++isig) {
bc3cacd4 439 // printf("INDEXES: %d-%d both%d\n",indexes1[isig],indexes2[isig],indexes12[isig]);
ba15fdfb 440 // mix the particles which satisfy only one of the signal branches
441 for(Int_t i1=0;i1<indexes1[isig];++i1) {
bc3cacd4 442 if(!indexes2[isig]) FillMCHistograms(labels1[isig][i1], -1, isig); // (e.g. single electrons only, no pairs)
443 for(Int_t i2=0;i2<indexes2[isig];++i2) {
88204efa 444 if(bFillCF) fCfManagerPair->FillMC(labels1[isig][i1], labels2[isig][i2], isig);
445 if(bFillHF) fHistoArray->Fill(labels1[isig][i1], labels2[isig][i2], isig);
833da63d 446 FillMCHistograms(labels1[isig][i1], labels2[isig][i2], isig);
ba15fdfb 447 }
448 }
449 // mix the particles which satisfy both branches
450 for(Int_t i1=0;i1<indexes12[isig];++i1) {
451 for(Int_t i2=0; i2<i1; ++i2) {
88204efa 452 if(bFillCF) fCfManagerPair->FillMC(labels12[isig][i1], labels12[isig][i2], isig);
453 if(bFillHF) fHistoArray->Fill(labels12[isig][i1], labels12[isig][i2], isig);
833da63d 454 FillMCHistograms(labels12[isig][i1], labels12[isig][i2], isig);
ba15fdfb 455 }
456 }
457 } // end loop over signals
458
459 // release the memory
460 for(Int_t isig=0;isig<nSignals;++isig) {
461 delete [] *(labels1+isig);
462 delete [] *(labels2+isig);
463 delete [] *(labels12+isig);
464 }
465 delete [] labels1;
466 delete [] labels2;
467 delete [] labels12;
468 delete [] indexes1;
469 delete [] indexes2;
470 delete [] indexes12;
b2a297fa 471}
472
554e40f8 473//________________________________________________________________
474void AliDielectron::FillHistogramsTracks(TObjArray **tracks)
475{
476 //
477 // Fill Histogram information for tracks after prefilter
478 // ignore mixed events - for prefilter, only single tracks +/- are relevant
479 //
480
481 TString className,className2;
482 Double_t values[AliDielectronVarManager::kNMaxValues];
483
484 //Fill track information, separately for the track array candidates
485 for (Int_t i=0; i<2; ++i){
486 className.Form("Pre_%s",fgkTrackClassNames[i]);
487 if (!fHistos->GetHistogramList()->FindObject(className.Data())) continue;
488 Int_t ntracks=tracks[i]->GetEntriesFast();
489 for (Int_t itrack=0; itrack<ntracks; ++itrack){
490 AliDielectronVarManager::Fill(tracks[i]->UncheckedAt(itrack), values);
491 fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
492 }
493 }
494}
ba15fdfb 495
496
497//________________________________________________________________
5720c765 498void AliDielectron::FillHistogramsMC(const AliMCEvent *ev, AliVEvent *ev1)
ba15fdfb 499{
500 //
501 // Fill Histogram information for MCEvents
502 //
503
5720c765 504 Double_t values[AliDielectronVarManager::kNMaxValues]={0.};
ba15fdfb 505 // Fill event information
5720c765 506 AliDielectronVarManager::Fill(ev1, values); // ESD/AOD information
507 AliDielectronVarManager::Fill(ev, values); // MC truth info
ba15fdfb 508 if (fHistos->GetHistogramList()->FindObject("MCEvent"))
509 fHistos->FillClass("MCEvent", AliDielectronVarManager::kNMaxValues, values);
510}
511
512
b2a297fa 513//________________________________________________________________
c315310a 514void AliDielectron::FillHistograms(const AliVEvent *ev, Bool_t pairInfoOnly)
b2a297fa 515{
516 //
517 // Fill Histogram information for tracks and pairs
518 //
519
61d106d3 520 TString className,className2;
5720c765 521 Double_t values[AliDielectronVarManager::kNMaxValues]={0.};
d327d9cd 522
0caf5fbb 523 //Fill event information
524 if (ev){
606603d9 525 if (fHistos->GetHistogramList()->FindObject("Event")) {
0caf5fbb 526 fHistos->FillClass("Event", AliDielectronVarManager::kNMaxValues, AliDielectronVarManager::GetData());
606603d9 527 }
5720c765 528 }
0caf5fbb 529
b2a297fa 530 //Fill track information, separately for the track array candidates
5720c765 531 if (!pairInfoOnly){
4533e78e 532 className2.Form("Track_%s",fgkPairClassNames[1]); // unlike sign, SE only
5720c765 533 for (Int_t i=0; i<4; ++i){
534 className.Form("Track_%s",fgkTrackClassNames[i]);
4533e78e 535 Bool_t mergedtrkClass=fHistos->GetHistogramList()->FindObject(className2.Data())!=0x0;
536 Bool_t trkClass=fHistos->GetHistogramList()->FindObject(className.Data())!=0x0;
537 if (!trkClass && !mergedtrkClass) continue;
5720c765 538 Int_t ntracks=fTracks[i].GetEntriesFast();
539 for (Int_t itrack=0; itrack<ntracks; ++itrack){
540 AliDielectronVarManager::Fill(fTracks[i].UncheckedAt(itrack), values);
4533e78e 541 if(trkClass)
542 fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
543 if(mergedtrkClass && i<2)
544 fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values); //only ev1
5720c765 545 }
b2a297fa 546 }
547 }
548
61d106d3 549 //Fill Pair information, separately for all pair candidate arrays and the legs
550 TObjArray arrLegs(100);
b2a297fa 551 for (Int_t i=0; i<10; ++i){
552 className.Form("Pair_%s",fgkPairClassNames[i]);
61d106d3 553 className2.Form("Track_Legs_%s",fgkPairClassNames[i]);
554 Bool_t pairClass=fHistos->GetHistogramList()->FindObject(className.Data())!=0x0;
555 Bool_t legClass=fHistos->GetHistogramList()->FindObject(className2.Data())!=0x0;
556 if (!pairClass&&!legClass) continue;
b2a297fa 557 Int_t ntracks=PairArray(i)->GetEntriesFast();
558 for (Int_t ipair=0; ipair<ntracks; ++ipair){
61d106d3 559 AliDielectronPair *pair=static_cast<AliDielectronPair*>(PairArray(i)->UncheckedAt(ipair));
560
561 //fill pair information
562 if (pairClass){
563 AliDielectronVarManager::Fill(pair, values);
564 fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
565 }
566
567 //fill leg information, don't fill the information twice
568 if (legClass){
569 AliVParticle *d1=pair->GetFirstDaughter();
570 AliVParticle *d2=pair->GetSecondDaughter();
571 if (!arrLegs.FindObject(d1)){
572 AliDielectronVarManager::Fill(d1, values);
573 fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
574 arrLegs.Add(d1);
575 }
576 if (!arrLegs.FindObject(d2)){
577 AliDielectronVarManager::Fill(d2, values);
578 fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
579 arrLegs.Add(d2);
580 }
581 }
b2a297fa 582 }
61d106d3 583 if (legClass) arrLegs.Clear();
b2a297fa 584 }
585
586}
2a14a7b1 587//________________________________________________________________
ffbede40 588void AliDielectron::FillHistogramsPair(AliDielectronPair *pair,Bool_t fromPreFilter/*=kFALSE*/)
2a14a7b1 589{
590 //
591 // Fill Histogram information for pairs and the track in the pair
592 // NOTE: in this funtion the leg information may be filled multiple
593 // times. This funtion is used in the track rotation pairing
594 // and those legs are not saved!
ffbede40 595 //
2a14a7b1 596 TString className,className2;
597 Double_t values[AliDielectronVarManager::kNMaxValues];
598
599 //Fill Pair information, separately for all pair candidate arrays and the legs
600 TObjArray arrLegs(100);
601 const Int_t type=pair->GetType();
ffbede40 602 if (fromPreFilter) {
603 className.Form("RejPair_%s",fgkPairClassNames[type]);
604 className2.Form("RejTrack_%s",fgkPairClassNames[type]);
605 } else {
606 className.Form("Pair_%s",fgkPairClassNames[type]);
607 className2.Form("Track_Legs_%s",fgkPairClassNames[type]);
608 }
2a14a7b1 609
610 Bool_t pairClass=fHistos->GetHistogramList()->FindObject(className.Data())!=0x0;
611 Bool_t legClass=fHistos->GetHistogramList()->FindObject(className2.Data())!=0x0;
612
613 //fill pair information
614 if (pairClass){
615 AliDielectronVarManager::Fill(pair, values);
616 fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
617 }
618
619 if (legClass){
620 AliVParticle *d1=pair->GetFirstDaughter();
621 AliDielectronVarManager::Fill(d1, values);
622 fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
623
624 AliVParticle *d2=pair->GetSecondDaughter();
625 AliDielectronVarManager::Fill(d2, values);
626 fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
627 }
628}
b2a297fa 629
630//________________________________________________________________
631void AliDielectron::FillTrackArrays(AliVEvent * const ev, Int_t eventNr)
632{
633 //
634 // select tracks and fill track candidate arrays
635 // eventNr = 0: First event, use track arrays 0 and 1
636 // eventNr = 1: Second event, use track arrays 2 and 3
637 //
99345a64 638
7ab4dc00 639 Int_t ntracks=ev->GetNumberOfTracks();
99345a64 640
b2a297fa 641 UInt_t selectedMask=(1<<fTrackFilter.GetCuts()->GetEntries())-1;
642 for (Int_t itrack=0; itrack<ntracks; ++itrack){
643 //get particle
644 AliVParticle *particle=ev->GetTrack(itrack);
67fd1119 645
b2a297fa 646 //apply track cuts
d8bb1abe 647 UInt_t cutmask=fTrackFilter.IsSelected(particle);
648 //fill cut QA
649 if(fCutQA) fQAmonitor->FillAll(particle);
650 if(fCutQA) fQAmonitor->Fill(cutmask,particle);
651
652 if (cutmask!=selectedMask) continue;
5720c765 653
b2a297fa 654 //fill selected particle into the corresponding track arrays
655 Short_t charge=particle->Charge();
656 if (charge>0) fTracks[eventNr*2].Add(particle);
657 else if (charge<0) fTracks[eventNr*2+1].Add(particle);
d8bb1abe 658
659
b2a297fa 660 }
661}
662
61d106d3 663//________________________________________________________________
4533e78e 664void AliDielectron::EventPlanePreFilter(Int_t arr1, Int_t arr2, TObjArray arrTracks1, TObjArray arrTracks2, const AliVEvent *ev)
2a14a7b1 665{
61d106d3 666 //
5720c765 667 // Prefilter tracks and tracks from pairs
668 // Needed for rejection in the Q-Vector of the event plane
669 // remove contribution of all tracks to the Q-vector that are in invariant mass window
61d106d3 670 //
0caf5fbb 671
5720c765 672 AliEventplane *evplane = const_cast<AliVEvent *>(ev)->GetEventplane();
4533e78e 673 if(!evplane) { // nanoAODs , here we do NOT have sub event reaction planes
674 // if(1) {
675 // get the EPselectionTask for recalculation of weighting factors
676 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
677 AliEPSelectionTask *eptask = dynamic_cast<AliEPSelectionTask *>(man->GetTask("EventplaneSelection"));
678 if(!eptask) return;
679
680 // track mapping
681 TMap mapRemovedTracks;
682
683 // eta gap ?
f3f5888a 684 //Bool_t etagap = kFALSE;
685 //for (Int_t iCut=0; iCut<fEventPlanePreFilter.GetCuts()->GetEntries();++iCut) {
686 // TString cutName=fEventPlanePreFilter.GetCuts()->At(iCut)->GetName();
687 // if(cutName.Contains("eta") || cutName.Contains("Eta")) etagap=kTRUE;
688 //}
5720c765 689
4533e78e 690 Double_t cQX=0., cQY=0.;
691 // apply cuts to the tracks, e.g. etagap
692 if(fEventPlanePreFilter.GetCuts()->GetEntries()) {
693 UInt_t selectedMask=(1<<fEventPlanePreFilter.GetCuts()->GetEntries())-1;
694 Int_t ntracks=ev->GetNumberOfTracks();
695 for (Int_t itrack=0; itrack<ntracks; ++itrack){
696 AliVParticle *particle=ev->GetTrack(itrack);
697 AliVTrack *track= static_cast<AliVTrack*>(particle);
698 if (!track) continue;
699 //event plane cuts
700 UInt_t cutMask=fEventPlanePreFilter.IsSelected(track);
701 //apply cut
702 if (cutMask==selectedMask) continue;
703
704 mapRemovedTracks.Add(track,track);
705 cQX += (eptask->GetWeight(track) * TMath::Cos(2*track->Phi()));
706 cQY += (eptask->GetWeight(track) * TMath::Sin(2*track->Phi()));
707 }
708 }
5720c765 709
4533e78e 710 // POI (particle of interest) rejection
711 Int_t pairIndex=GetPairIndex(arr1,arr2);
712
713 Int_t ntrack1=arrTracks1.GetEntriesFast();
714 Int_t ntrack2=arrTracks2.GetEntriesFast();
715 AliDielectronPair candidate;
716 candidate.SetKFUsage(fUseKF);
717
718 UInt_t selectedMask=(1<<fEventPlanePOIPreFilter.GetCuts()->GetEntries())-1;
719 for (Int_t itrack1=0; itrack1<ntrack1; ++itrack1){
720 Int_t end=ntrack2;
721 if (arr1==arr2) end=itrack1;
722 Bool_t accepted=kFALSE;
723 for (Int_t itrack2=0; itrack2<end; ++itrack2){
724 TObject *track1=arrTracks1.UncheckedAt(itrack1);
725 TObject *track2=arrTracks2.UncheckedAt(itrack2);
726 if (!track1 || !track2) continue;
727 //create the pair
728 candidate.SetTracks(static_cast<AliVTrack*>(track1), fPdgLeg1,
729 static_cast<AliVTrack*>(track2), fPdgLeg2);
730 candidate.SetType(pairIndex);
731 candidate.SetLabel(AliDielectronMC::Instance()->GetLabelMotherWithPdg(&candidate,fPdgMother));
732
733 //event plane pair cuts
734 UInt_t cutMask=fEventPlanePOIPreFilter.IsSelected(&candidate);
735 //apply cut
736 if (cutMask==selectedMask) continue;
737
738 accepted=kTRUE;
739 //remove the tracks from the Track arrays
740 arrTracks2.AddAt(0x0,itrack2);
741 }
742 if ( accepted ) arrTracks1.AddAt(0x0,itrack1);
743 }
744 //compress the track arrays
745 arrTracks1.Compress();
746 arrTracks2.Compress();
747
748 //Modify the components: subtract the tracks
749 ntrack1=arrTracks1.GetEntriesFast();
750 ntrack2=arrTracks2.GetEntriesFast();
751 // remove leg1 contribution
752 for (Int_t itrack=0; itrack<ntrack1; ++itrack){
753 AliVTrack *track= static_cast<AliVTrack*>(arrTracks1.UncheckedAt(itrack));
5720c765 754 if (!track) continue;
4533e78e 755 // track contribution was already removed
756 if (mapRemovedTracks.FindObject(track)) continue;
757 else mapRemovedTracks.Add(track,track);
5720c765 758
4533e78e 759 cQX += (eptask->GetWeight(track) * TMath::Cos(2*track->Phi()));
760 cQY += (eptask->GetWeight(track) * TMath::Sin(2*track->Phi()));
5720c765 761 }
4533e78e 762 // remove leg2 contribution
763 for (Int_t itrack=0; itrack<ntrack2; ++itrack){
764 AliVTrack *track= static_cast<AliVTrack*>(arrTracks2.UncheckedAt(itrack));
5720c765 765 if (!track) continue;
4533e78e 766 // track contribution was already removed
767 if (mapRemovedTracks.FindObject(track)) continue;
768 else mapRemovedTracks.Add(track,track);
769
770 cQX += (eptask->GetWeight(track) * TMath::Cos(2*track->Phi()));
771 cQY += (eptask->GetWeight(track) * TMath::Sin(2*track->Phi()));
5720c765 772 }
5720c765 773
4533e78e 774 // build a corrected alieventplane using the values from the var manager
775 // these uncorrected values are filled using the stored magnitude and angle in the header
776 TVector2 qcorr;
777 qcorr.Set(AliDielectronVarManager::GetValue(AliDielectronVarManager::kTPCxH2uc)-cQX,
778 AliDielectronVarManager::GetValue(AliDielectronVarManager::kTPCyH2uc)-cQY);
779 // fill alieventplane
780 AliEventplane cevplane;
781 cevplane.SetQVector(&qcorr);
782 AliDielectronVarManager::SetTPCEventPlane(&cevplane);
783 cevplane.SetQVector(0);
784 return;
785 } //end: nanoAODs
786 else
787 {
788 // this is done in case of ESDs or AODs
789 Bool_t isESD=(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()->IsA()==AliESDInputHandler::Class());
790 // copy event plane object
791 AliEventplane cevplane(*evplane);
792 // Int_t nMaxID=cevplane->GetQContributionXArray()->GetSize();
793
794 TVector2 *qcorr = cevplane.GetQVector();
795 if(!qcorr) return;
796 TVector2 *qcsub1 = 0x0;
797 TVector2 *qcsub2 = 0x0;
798
799 // eta gap ?
800 Bool_t etagap = kFALSE;
801 for (Int_t iCut=0; iCut<fEventPlanePreFilter.GetCuts()->GetEntries();++iCut) {
802 TString cutName=fEventPlanePreFilter.GetCuts()->At(iCut)->GetName();
803 if(cutName.Contains("eta") || cutName.Contains("Eta")) etagap=kTRUE;
61d106d3 804 }
4533e78e 805
806 // subevent configuration for eta gap or LS (default is rndm)
807 if(fLikeSignSubEvents && etagap) {
808 // start with the full Qvector/event in both sub events
809 qcsub1 = new TVector2(*qcorr);
810 qcsub2 = new TVector2(*qcorr);
811 cevplane.SetQsub(qcsub1,qcsub2);
812
813 Int_t ntracks=ev->GetNumberOfTracks();
814 // track removals
815 for (Int_t itrack=0; itrack<ntracks; ++itrack){
816 AliVParticle *particle=ev->GetTrack(itrack);
817 AliVTrack *track= static_cast<AliVTrack*>(particle);
818 if (!track) continue;
819 if (track->GetID()>=0 && !isESD) continue;
820 Int_t tmpID = isESD ? track->GetID() : track->GetID()*-1 - 1;
821
822 // set contributions to zero
823 // charge sub1+ sub2-
824 if(fLikeSignSubEvents) {
825 Short_t charge=track->Charge();
826 if (charge<0) {
827 cevplane.GetQContributionXArraysub1()->SetAt(0.0, tmpID);
828 cevplane.GetQContributionYArraysub1()->SetAt(0.0, tmpID);
829 }
830 if (charge>0) {
831 cevplane.GetQContributionXArraysub2()->SetAt(0.0, tmpID);
832 cevplane.GetQContributionYArraysub2()->SetAt(0.0, tmpID);
833 }
834 }
835 // eta sub1+ sub2-
836 if(etagap) {
837 Double_t eta=track->Eta();
838 if (eta<0.0) {
839 cevplane.GetQContributionXArraysub1()->SetAt(0.0, tmpID);
840 cevplane.GetQContributionYArraysub1()->SetAt(0.0, tmpID);
841 }
842 if (eta>0.0) {
843 cevplane.GetQContributionXArraysub2()->SetAt(0.0, tmpID);
844 cevplane.GetQContributionYArraysub2()->SetAt(0.0, tmpID);
845 }
846 }
847 } // end: loop over tracks
848 } // end: sub event configuration
849
850 // apply cuts, e.g. etagap
851 if(fEventPlanePreFilter.GetCuts()->GetEntries()) {
852 UInt_t selectedMask=(1<<fEventPlanePreFilter.GetCuts()->GetEntries())-1;
853 Int_t ntracks=ev->GetNumberOfTracks();
854 for (Int_t itrack=0; itrack<ntracks; ++itrack){
855 AliVParticle *particle=ev->GetTrack(itrack);
856 AliVTrack *track= static_cast<AliVTrack*>(particle);
857 if (!track) continue;
858 if (track->GetID()>=0 && !isESD) continue;
859 Int_t tmpID = isESD ? track->GetID() : track->GetID()*-1 - 1;
860
861 //event plane cuts
862 UInt_t cutMask=fEventPlanePreFilter.IsSelected(track);
863 //apply cut
864 if (cutMask==selectedMask) continue;
865
866 // set contributions to zero
867 cevplane.GetQContributionXArray()->SetAt(0.0, tmpID);
868 cevplane.GetQContributionYArray()->SetAt(0.0, tmpID);
869 cevplane.GetQContributionXArraysub1()->SetAt(0.0, tmpID);
870 cevplane.GetQContributionYArraysub1()->SetAt(0.0, tmpID);
871 cevplane.GetQContributionXArraysub2()->SetAt(0.0, tmpID);
872 cevplane.GetQContributionYArraysub2()->SetAt(0.0, tmpID);
873 }
874 } // end: track cuts
875
876 // POI (particle of interest) rejection
877 Int_t pairIndex=GetPairIndex(arr1,arr2);
878 Int_t ntrack1=arrTracks1.GetEntriesFast();
879 Int_t ntrack2=arrTracks2.GetEntriesFast();
880 AliDielectronPair candidate;
881 candidate.SetKFUsage(fUseKF);
882
883 UInt_t selectedMask=(1<<fEventPlanePOIPreFilter.GetCuts()->GetEntries())-1;
884 for (Int_t itrack1=0; itrack1<ntrack1; ++itrack1){
885 Int_t end=ntrack2;
886 if (arr1==arr2) end=itrack1;
887 Bool_t accepted=kFALSE;
888 for (Int_t itrack2=0; itrack2<end; ++itrack2){
889 TObject *track1=arrTracks1.UncheckedAt(itrack1);
890 TObject *track2=arrTracks2.UncheckedAt(itrack2);
891 if (!track1 || !track2) continue;
892 //create the pair
893 candidate.SetTracks(static_cast<AliVTrack*>(track1), fPdgLeg1,
894 static_cast<AliVTrack*>(track2), fPdgLeg2);
895
896 candidate.SetType(pairIndex);
897 candidate.SetLabel(AliDielectronMC::Instance()->GetLabelMotherWithPdg(&candidate,fPdgMother));
898
899 //event plane cuts
900 UInt_t cutMask=fEventPlanePOIPreFilter.IsSelected(&candidate);
901 //apply cut
902 if (cutMask==selectedMask) continue;
903
904 accepted=kTRUE;
905 //remove the tracks from the Track arrays
906 arrTracks2.AddAt(0x0,itrack2);
907 }
5720c765 908 if ( accepted ) arrTracks1.AddAt(0x0,itrack1);
4533e78e 909 }
910 //compress the track arrays
911 arrTracks1.Compress();
912 arrTracks2.Compress();
913
914 //Modify the components: subtract the tracks
915 ntrack1=arrTracks1.GetEntriesFast();
916 ntrack2=arrTracks2.GetEntriesFast();
917 // remove leg1 contribution
918 for (Int_t itrack=0; itrack<ntrack1; ++itrack){
919 AliVTrack *track= static_cast<AliVTrack*>(arrTracks1.UncheckedAt(itrack));
920 if (!track) continue;
921 if (track->GetID()>=0 && !isESD) continue;
922 Int_t tmpID = isESD ? track->GetID() : track->GetID()*-1 - 1;
923 // set contributions to zero
924 cevplane.GetQContributionXArray()->SetAt(0.0, tmpID);
925 cevplane.GetQContributionYArray()->SetAt(0.0, tmpID);
926 cevplane.GetQContributionXArraysub1()->SetAt(0.0, tmpID);
927 cevplane.GetQContributionYArraysub1()->SetAt(0.0, tmpID);
928 cevplane.GetQContributionXArraysub2()->SetAt(0.0, tmpID);
929 cevplane.GetQContributionYArraysub2()->SetAt(0.0, tmpID);
930 }
931 // remove leg2 contribution
932 for (Int_t itrack=0; itrack<ntrack2; ++itrack){
933 AliVTrack *track= static_cast<AliVTrack*>(arrTracks2.UncheckedAt(itrack));
934 if (!track) continue;
935 if (track->GetID()>=0 && !isESD) continue;
936 Int_t tmpID = isESD ? track->GetID() : track->GetID()*-1 - 1;
937 // set contributions to zero
938 cevplane.GetQContributionXArray()->SetAt(0.0, tmpID);
939 cevplane.GetQContributionYArray()->SetAt(0.0, tmpID);
940 cevplane.GetQContributionXArraysub1()->SetAt(0.0, tmpID);
941 cevplane.GetQContributionYArraysub1()->SetAt(0.0, tmpID);
942 cevplane.GetQContributionXArraysub2()->SetAt(0.0, tmpID);
943 cevplane.GetQContributionYArraysub2()->SetAt(0.0, tmpID);
944 }
554e40f8 945
4533e78e 946 // set corrected AliEventplane and fill variables with corrected values
947 AliDielectronVarManager::SetTPCEventPlane(&cevplane);
948 delete qcsub1;
949 delete qcsub2;
950 } // end: ESD or AOD case
554e40f8 951
4533e78e 952}
5720c765 953//________________________________________________________________
954void AliDielectron::PairPreFilter(Int_t arr1, Int_t arr2, TObjArray &arrTracks1, TObjArray &arrTracks2)
955{
956 //
957 // Prefilter tracks from pairs
958 // Needed for datlitz rejections
959 // remove all tracks from the Single track arrays that pass the cuts in this filter
960 //
961
962 Int_t ntrack1=arrTracks1.GetEntriesFast();
963 Int_t ntrack2=arrTracks2.GetEntriesFast();
964 AliDielectronPair candidate;
08b801a6 965 candidate.SetKFUsage(fUseKF);
5720c765 966 // flag arrays for track removal
967 Bool_t *bTracks1 = new Bool_t[ntrack1];
968 for (Int_t itrack1=0; itrack1<ntrack1; ++itrack1) bTracks1[itrack1]=kFALSE;
969 Bool_t *bTracks2 = new Bool_t[ntrack2];
970 for (Int_t itrack2=0; itrack2<ntrack2; ++itrack2) bTracks2[itrack2]=kFALSE;
971
972 UInt_t selectedMask=(1<<fPairPreFilter.GetCuts()->GetEntries())-1;
973 UInt_t selectedMaskPair=(1<<fPairFilter.GetCuts()->GetEntries())-1;
974
975 Int_t nRejPasses = 1; //for fPreFilterUnlikeOnly and no set flag
976 if (fPreFilterAllSigns) nRejPasses = 3;
977
978 for (Int_t iRP=0; iRP < nRejPasses; ++iRP) {
979 Int_t arr1RP=arr1, arr2RP=arr2;
980 TObjArray *arrTracks1RP=&arrTracks1;
981 TObjArray *arrTracks2RP=&arrTracks2;
982 Bool_t *bTracks1RP = bTracks1;
983 Bool_t *bTracks2RP = bTracks2;
984 switch (iRP) {
985 case 1: arr1RP=arr1;arr2RP=arr1;
986 arrTracks1RP=&arrTracks1;
987 arrTracks2RP=&arrTracks1;
988 bTracks1RP = bTracks1;
989 bTracks2RP = bTracks1;
990 break;
991 case 2: arr1RP=arr2;arr2RP=arr2;
992 arrTracks1RP=&arrTracks2;
993 arrTracks2RP=&arrTracks2;
994 bTracks1RP = bTracks2;
995 bTracks2RP = bTracks2;
996 break;
997 default: ;//nothing to do
998 }
999 Int_t ntrack1RP=(*arrTracks1RP).GetEntriesFast();
1000 Int_t ntrack2RP=(*arrTracks2RP).GetEntriesFast();
1001
1002 Int_t pairIndex=GetPairIndex(arr1RP,arr2RP);
1003
1004 for (Int_t itrack1=0; itrack1<ntrack1RP; ++itrack1){
1005 Int_t end=ntrack2RP;
1006 if (arr1RP==arr2RP) end=itrack1;
1007 for (Int_t itrack2=0; itrack2<end; ++itrack2){
1008 TObject *track1=(*arrTracks1RP).UncheckedAt(itrack1);
1009 TObject *track2=(*arrTracks2RP).UncheckedAt(itrack2);
1010 if (!track1 || !track2) continue;
1011 //create the pair
1012 candidate.SetTracks(static_cast<AliVTrack*>(track1), fPdgLeg1,
1013 static_cast<AliVTrack*>(track2), fPdgLeg2);
1014
1015 candidate.SetType(pairIndex);
1016 candidate.SetLabel(AliDielectronMC::Instance()->GetLabelMotherWithPdg(&candidate,fPdgMother));
1017 //relate to the production vertex
1018 // if (AliDielectronVarManager::GetKFVertex()) candidate.SetProductionVertex(*AliDielectronVarManager::GetKFVertex());
1019
1020 //pair cuts
1021 UInt_t cutMask=fPairPreFilter.IsSelected(&candidate);
1022
1023 //apply cut
1024 if (cutMask!=selectedMask) continue;
1025 if (fCfManagerPair) fCfManagerPair->Fill(selectedMaskPair+1 ,&candidate);
1026 if (fHistos) FillHistogramsPair(&candidate,kTRUE);
1027 //set flags for track removal
1028 bTracks1RP[itrack1]=kTRUE;
1029 bTracks2RP[itrack2]=kTRUE;
1030 }
1031 }
1032 }
1033
1034 //remove the tracks from the Track arrays
1035 for (Int_t itrack1=0; itrack1<ntrack1; ++itrack1){
1036 if(bTracks1[itrack1]) arrTracks1.AddAt(0x0, itrack1);
1037 }
1038 for (Int_t itrack2=0; itrack2<ntrack2; ++itrack2){
1039 if(bTracks2[itrack2]) arrTracks2.AddAt(0x0, itrack2);
1040 }
1041
1042 // clean up
1043 delete [] bTracks1;
1044 delete [] bTracks2;
1045
1046 //compress the track arrays
61d106d3 1047 arrTracks1.Compress();
1048 arrTracks2.Compress();
1049
1050 //apply leg cuts after the pre filter
1051 if ( fPairPreFilterLegs.GetCuts()->GetEntries()>0 ) {
1052 selectedMask=(1<<fPairPreFilterLegs.GetCuts()->GetEntries())-1;
1053 //loop over tracks from array 1
1054 for (Int_t itrack=0; itrack<arrTracks1.GetEntriesFast();++itrack){
1055 //test cuts
1056 UInt_t cutMask=fPairPreFilterLegs.IsSelected(arrTracks1.UncheckedAt(itrack));
554e40f8 1057
61d106d3 1058 //apply cut
1059 if (cutMask!=selectedMask) arrTracks1.AddAt(0x0,itrack);;
1060 }
1061 arrTracks1.Compress();
1062
1063 //in case of like sign don't loop over second array
1064 if (arr1==arr2) {
1065 arrTracks2=arrTracks1;
1066 } else {
1067
1068 //loop over tracks from array 2
1069 for (Int_t itrack=0; itrack<arrTracks2.GetEntriesFast();++itrack){
1070 //test cuts
1071 UInt_t cutMask=fPairPreFilterLegs.IsSelected(arrTracks2.UncheckedAt(itrack));
1072 //apply cut
1073 if (cutMask!=selectedMask) arrTracks2.AddAt(0x0,itrack);
1074 }
1075 arrTracks2.Compress();
1076
1077 }
1078 }
554e40f8 1079 //For unlike-sign monitor track-cuts:
ba15fdfb 1080 if (arr1!=arr2&&fHistos) {
554e40f8 1081 TObjArray *unlikesignArray[2] = {&arrTracks1,&arrTracks2};
1082 FillHistogramsTracks(unlikesignArray);
1083 }
61d106d3 1084}
1085
b2a297fa 1086//________________________________________________________________
2a14a7b1 1087void AliDielectron::FillPairArrays(Int_t arr1, Int_t arr2)
1088{
b2a297fa 1089 //
1090 // select pairs and fill pair candidate arrays
1091 //
61d106d3 1092
1093 TObjArray arrTracks1=fTracks[arr1];
1094 TObjArray arrTracks2=fTracks[arr2];
1095
1096 //process pre filter if set
ba15fdfb 1097 if ((!fPreFilterAllSigns) && (!fPreFilterUnlikeOnly) && ( fPairPreFilter.GetCuts()->GetEntries()>0 )) PairPreFilter(arr1, arr2, arrTracks1, arrTracks2);
61d106d3 1098
b2a297fa 1099 Int_t pairIndex=GetPairIndex(arr1,arr2);
1100
61d106d3 1101 Int_t ntrack1=arrTracks1.GetEntriesFast();
1102 Int_t ntrack2=arrTracks2.GetEntriesFast();
b2a297fa 1103
1104 AliDielectronPair *candidate=new AliDielectronPair;
08b801a6 1105 candidate->SetKFUsage(fUseKF);
b2a297fa 1106
1107 UInt_t selectedMask=(1<<fPairFilter.GetCuts()->GetEntries())-1;
1108
1109 for (Int_t itrack1=0; itrack1<ntrack1; ++itrack1){
1110 Int_t end=ntrack2;
1111 if (arr1==arr2) end=itrack1;
1112 for (Int_t itrack2=0; itrack2<end; ++itrack2){
8df8e382 1113 //create the pair
61d106d3 1114 candidate->SetTracks(static_cast<AliVTrack*>(arrTracks1.UncheckedAt(itrack1)), fPdgLeg1,
5720c765 1115 static_cast<AliVTrack*>(arrTracks2.UncheckedAt(itrack2)), fPdgLeg2);
b2a297fa 1116 candidate->SetType(pairIndex);
e4339752 1117 Int_t label=AliDielectronMC::Instance()->GetLabelMotherWithPdg(candidate,fPdgMother);
1118 candidate->SetLabel(label);
1119 if (label>-1) candidate->SetPdgCode(fPdgMother);
b2a297fa 1120
4fae8ef9 1121 // check for gamma kf particle
1122 label=AliDielectronMC::Instance()->GetLabelMotherWithPdg(candidate,22);
1123 if (label>-1) {
1124 candidate->SetGammaTracks(static_cast<AliVTrack*>(arrTracks1.UncheckedAt(itrack1)), fPdgLeg1,
1125 static_cast<AliVTrack*>(arrTracks2.UncheckedAt(itrack2)), fPdgLeg2);
1126 // should we set the pdgmothercode and the label
1127 }
1128
b2a297fa 1129 //pair cuts
1130 UInt_t cutMask=fPairFilter.IsSelected(candidate);
1131
1132 //CF manager for the pair
1133 if (fCfManagerPair) fCfManagerPair->Fill(cutMask,candidate);
620b34a5 1134 //histogram array for the pair
1135 if (fHistoArray) fHistoArray->Fill(pairIndex,candidate);
d8bb1abe 1136 // cut qa
dd718b9c 1137 if(pairIndex==kEv1PM && fCutQA) {
d8bb1abe 1138 fQAmonitor->FillAll(candidate);
1139 fQAmonitor->Fill(cutMask,candidate);
1140 }
b2a297fa 1141
1142 //apply cut
1143 if (cutMask!=selectedMask) continue;
1144
1145 //add the candidate to the candidate array
1146 PairArray(pairIndex)->Add(candidate);
1147 //get a new candidate
1148 candidate=new AliDielectronPair;
08b801a6 1149 candidate->SetKFUsage(fUseKF);
b2a297fa 1150 }
1151 }
1152 //delete the surplus candidate
1153 delete candidate;
1154}
1155
2a14a7b1 1156//________________________________________________________________
1157void AliDielectron::FillPairArrayTR()
1158{
1159 //
1160 // select pairs and fill pair candidate arrays
1161 //
1162 UInt_t selectedMask=(1<<fPairFilter.GetCuts()->GetEntries())-1;
1163
1164 while ( fTrackRotator->NextCombination() ){
1165 AliDielectronPair candidate;
08b801a6 1166 candidate.SetKFUsage(fUseKF);
1201a1a9 1167 candidate.SetTracks(&fTrackRotator->GetKFTrackP(), &fTrackRotator->GetKFTrackN(),
1168 fTrackRotator->GetVTrackP(),fTrackRotator->GetVTrackN());
ffbede40 1169 candidate.SetType(kEv1PMRot);
2a14a7b1 1170
1171 //pair cuts
1172 UInt_t cutMask=fPairFilter.IsSelected(&candidate);
1173
1174 //CF manager for the pair
1175 if (fCfManagerPair) fCfManagerPair->Fill(cutMask,&candidate);
620b34a5 1176 //histogram array for the pair
1177 if (fHistoArray) fHistoArray->Fill((Int_t)kEv1PMRot,&candidate);
2a14a7b1 1178
1179 //apply cut
5720c765 1180 if (cutMask==selectedMask) {
1181 if(fHistos) FillHistogramsPair(&candidate);
e4339752 1182 if(fStoreRotatedPairs) PairArray(kEv1PMRot)->Add(new AliDielectronPair(candidate));
5720c765 1183 }
2a14a7b1 1184 }
1185}
1186
572b0139 1187//________________________________________________________________
1188void AliDielectron::FillDebugTree()
1189{
1190 //
1191 // Fill Histogram information for tracks and pairs
1192 //
1193
1194 //Fill Debug tree
1195 for (Int_t i=0; i<10; ++i){
1196 Int_t ntracks=PairArray(i)->GetEntriesFast();
1197 for (Int_t ipair=0; ipair<ntracks; ++ipair){
1198 fDebugTree->Fill(static_cast<AliDielectronPair*>(PairArray(i)->UncheckedAt(ipair)));
1199 }
1200 }
1201}
b2a297fa 1202
572b0139 1203//________________________________________________________________
1204void AliDielectron::SaveDebugTree()
1205{
1206 //
1207 // delete the debug tree, this will also write the tree
1208 //
1209 if (fDebugTree) fDebugTree->DeleteStreamer();
1210}
b2a297fa 1211
ba15fdfb 1212
1213//__________________________________________________________________
1214void AliDielectron::AddSignalMC(AliDielectronSignalMC* signal) {
1215 //
1216 // Add an MC signal to the signals list
1217 //
1218 if(!fSignalsMC) {
1219 fSignalsMC = new TObjArray();
1220 fSignalsMC->SetOwner();
1221 }
1222 fSignalsMC->Add(signal);
1223}
833da63d 1224
1225//________________________________________________________________
1226void AliDielectron::FillMCHistograms(Int_t label1, Int_t label2, Int_t nSignal) {
1227 //
1228 // fill QA MC TRUTH histograms for pairs and legs of all added mc signals
1229 //
1230
1231 TString className,className2,className3;
1232 className.Form("Pair_%s_MCtruth",fSignalsMC->At(nSignal)->GetName());
1233 className2.Form("Track_Legs_%s_MCtruth",fSignalsMC->At(nSignal)->GetName());
1234 className3.Form("Track_%s_%s_MCtruth",fgkPairClassNames[1],fSignalsMC->At(nSignal)->GetName());
1235 Bool_t pairClass=fHistos->GetHistogramList()->FindObject(className.Data())!=0x0;
1236 Bool_t legClass=fHistos->GetHistogramList()->FindObject(className2.Data())!=0x0;
1237 Bool_t trkClass=fHistos->GetHistogramList()->FindObject(className3.Data())!=0x0;
bc3cacd4 1238 // printf("fill signal %d: pair %d legs %d trk %d \n",nSignal,pairClass,legClass,trkClass);
833da63d 1239 if(!pairClass && !legClass && !trkClass) return;
1240
1241 // printf("leg labels: %d-%d \n",label1,label2);
1242 AliVParticle* part1 = AliDielectronMC::Instance()->GetMCTrackFromMCEvent(label1);
1243 AliVParticle* part2 = AliDielectronMC::Instance()->GetMCTrackFromMCEvent(label2);
1244 if(!part1 && !part2) return;
1245 if(part1&&part2) {
1246 // fill only unlike sign (and only SE)
833da63d 1247 if(part1->Charge()*part2->Charge()>=0) return;
1248 }
1249
1250
1251 AliDielectronMC* dieMC = AliDielectronMC::Instance();
1252
1253 Int_t mLabel1 = dieMC->GetMothersLabel(label1); // should work for both ESD and AOD
1254 Int_t mLabel2 = dieMC->GetMothersLabel(label2);
1255
1256 // check the same mother option
1257 AliDielectronSignalMC* sigMC = (AliDielectronSignalMC*)fSignalsMC->At(nSignal);
1258 if(sigMC->GetMothersRelation()==AliDielectronSignalMC::kSame && mLabel1!=mLabel2) return;
1259 if(sigMC->GetMothersRelation()==AliDielectronSignalMC::kDifferent && mLabel1==mLabel2) return;
1260
bc3cacd4 1261 // fill event values
833da63d 1262 Double_t values[AliDielectronVarManager::kNMaxValues];
1263 AliDielectronVarManager::Fill(dieMC->GetMCEvent(), values); // get event informations
bc3cacd4 1264
1265 // fill the leg variables
1266 // printf("leg:%d trk:%d part1:%p part2:%p \n",legClass,trkClass,part1,part2);
833da63d 1267 if (legClass || trkClass) {
1268 if(part1) AliDielectronVarManager::Fill(part1,values);
1269 if(part1 && trkClass) fHistos->FillClass(className3, AliDielectronVarManager::kNMaxValues, values);
1270 if(part1 && part2 && legClass) fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
1271 if(part2) AliDielectronVarManager::Fill(part2,values);
1272 if(part2 && trkClass) fHistos->FillClass(className3, AliDielectronVarManager::kNMaxValues, values);
1273 if(part1 && part2 && legClass) fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
1274 }
1275
1276 //fill pair information
1277 if (pairClass && part1 && part2) {
1278 AliDielectronVarManager::FillVarMCParticle2(part1,part2,values);
1279 fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
1280 }
1281
1282}
1283
f0351f9c 1284//________________________________________________________________
1285void AliDielectron::FillMCHistograms(const AliVEvent *ev) {
1286 //
1287 // fill QA MC histograms for pairs and legs of all added mc signals
1288 //
1289
4c06df48 1290 if (!fSignalsMC) return;
833da63d 1291 TString className,className2,className3;
f0351f9c 1292 Double_t values[AliDielectronVarManager::kNMaxValues]={0.};
1293 AliDielectronVarManager::Fill(ev, values); // get event informations
1294 //loop over all added mc signals
1295 for(Int_t isig=0; isig<fSignalsMC->GetEntries(); isig++) {
1296
bc3cacd4 1297 //check if and what to fill
f0351f9c 1298 className.Form("Pair_%s",fSignalsMC->At(isig)->GetName());
1299 className2.Form("Track_Legs_%s",fSignalsMC->At(isig)->GetName());
536a997a 1300 className3.Form("Track_%s_%s",fgkPairClassNames[1],fSignalsMC->At(isig)->GetName()); // unlike sign, SE only
f0351f9c 1301 Bool_t pairClass=fHistos->GetHistogramList()->FindObject(className.Data())!=0x0;
1302 Bool_t legClass=fHistos->GetHistogramList()->FindObject(className2.Data())!=0x0;
536a997a 1303 Bool_t mergedtrkClass=fHistos->GetHistogramList()->FindObject(className3.Data())!=0x0;
1304 if(!pairClass && !legClass && !mergedtrkClass) continue;
1305
1306 // fill pair and/or their leg variables
1307 if(pairClass || legClass) {
1308 Int_t npairs=PairArray(AliDielectron::kEv1PM)->GetEntriesFast(); // only SE +-
1309 for (Int_t ipair=0; ipair<npairs; ++ipair){
1310 AliDielectronPair *pair=static_cast<AliDielectronPair*>(PairArray(AliDielectron::kEv1PM)->UncheckedAt(ipair));
1311
1312 Bool_t isMCtruth = AliDielectronMC::Instance()->IsMCTruth(pair, (AliDielectronSignalMC*)fSignalsMC->At(isig));
1313 if(isMCtruth) {
1314 //fill pair information
1315 if (pairClass){
1316 AliDielectronVarManager::Fill(pair, values);
1317 fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
1318 }
1319 //fill leg information, both + and - in the same histo
1320 if (legClass){
1321 AliDielectronVarManager::Fill(pair->GetFirstDaughter(),values);
1322 fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
1323 AliDielectronVarManager::Fill(pair->GetSecondDaughter(),values);
1324 fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
1325 }
1326 } //is signal
1327 } //loop: pairs
1328 }
833da63d 1329
1330 // fill single tracks of signals
bc3cacd4 1331 if(!mergedtrkClass) continue;
1332 // loop over SE track arrays
1333 for (Int_t i=0; i<2; ++i){
1334 Int_t ntracks=fTracks[i].GetEntriesFast();
1335 for (Int_t itrack=0; itrack<ntracks; ++itrack){
1336 Int_t label=((AliVParticle*)fTracks[i].UncheckedAt(itrack))->GetLabel();
1337 Bool_t isMCtruth1 = AliDielectronMC::Instance()->IsMCTruth(label, (AliDielectronSignalMC*)fSignalsMC->At(isig), 1);
1338 Bool_t isMCtruth2 = AliDielectronMC::Instance()->IsMCTruth(label, (AliDielectronSignalMC*)fSignalsMC->At(isig), 2);
1339 // skip if track does not correspond to the signal
1340 if(!isMCtruth1 && !isMCtruth2) continue;
1341 AliDielectronVarManager::Fill(fTracks[i].UncheckedAt(itrack), values);
1342 fHistos->FillClass(className3, AliDielectronVarManager::kNMaxValues, values);
536a997a 1343 } //loop: tracks
bc3cacd4 1344 } //loop: arrays
833da63d 1345
f0351f9c 1346 } //loop: MCsignals
1347
1348}