]> 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
Protection against re-initialization of histograms if not already
[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;
b2a297fa 188 if (fHistos) delete fHistos;
189 if (fPairCandidates) delete fPairCandidates;
572b0139 190 if (fDebugTree) delete fDebugTree;
5720c765 191 if (fMixing) delete fMixing;
ba15fdfb 192 if (fSignalsMC) delete fSignalsMC;
5720c765 193 if (fCfManagerPair) delete fCfManagerPair;
d4619886 194 if (fHistoArray) delete fHistoArray;
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
5720c765 683
4533e78e 684 Double_t cQX=0., cQY=0.;
685 // apply cuts to the tracks, e.g. etagap
686 if(fEventPlanePreFilter.GetCuts()->GetEntries()) {
687 UInt_t selectedMask=(1<<fEventPlanePreFilter.GetCuts()->GetEntries())-1;
688 Int_t ntracks=ev->GetNumberOfTracks();
689 for (Int_t itrack=0; itrack<ntracks; ++itrack){
690 AliVParticle *particle=ev->GetTrack(itrack);
691 AliVTrack *track= static_cast<AliVTrack*>(particle);
692 if (!track) continue;
693 //event plane cuts
694 UInt_t cutMask=fEventPlanePreFilter.IsSelected(track);
695 //apply cut
696 if (cutMask==selectedMask) continue;
697
698 mapRemovedTracks.Add(track,track);
699 cQX += (eptask->GetWeight(track) * TMath::Cos(2*track->Phi()));
700 cQY += (eptask->GetWeight(track) * TMath::Sin(2*track->Phi()));
701 }
702 }
5720c765 703
4533e78e 704 // POI (particle of interest) rejection
705 Int_t pairIndex=GetPairIndex(arr1,arr2);
706
707 Int_t ntrack1=arrTracks1.GetEntriesFast();
708 Int_t ntrack2=arrTracks2.GetEntriesFast();
709 AliDielectronPair candidate;
710 candidate.SetKFUsage(fUseKF);
711
712 UInt_t selectedMask=(1<<fEventPlanePOIPreFilter.GetCuts()->GetEntries())-1;
713 for (Int_t itrack1=0; itrack1<ntrack1; ++itrack1){
714 Int_t end=ntrack2;
715 if (arr1==arr2) end=itrack1;
716 Bool_t accepted=kFALSE;
717 for (Int_t itrack2=0; itrack2<end; ++itrack2){
718 TObject *track1=arrTracks1.UncheckedAt(itrack1);
719 TObject *track2=arrTracks2.UncheckedAt(itrack2);
720 if (!track1 || !track2) continue;
721 //create the pair
722 candidate.SetTracks(static_cast<AliVTrack*>(track1), fPdgLeg1,
723 static_cast<AliVTrack*>(track2), fPdgLeg2);
724 candidate.SetType(pairIndex);
725 candidate.SetLabel(AliDielectronMC::Instance()->GetLabelMotherWithPdg(&candidate,fPdgMother));
726
727 //event plane pair cuts
728 UInt_t cutMask=fEventPlanePOIPreFilter.IsSelected(&candidate);
729 //apply cut
730 if (cutMask==selectedMask) continue;
731
732 accepted=kTRUE;
733 //remove the tracks from the Track arrays
734 arrTracks2.AddAt(0x0,itrack2);
735 }
736 if ( accepted ) arrTracks1.AddAt(0x0,itrack1);
737 }
738 //compress the track arrays
739 arrTracks1.Compress();
740 arrTracks2.Compress();
741
742 //Modify the components: subtract the tracks
743 ntrack1=arrTracks1.GetEntriesFast();
744 ntrack2=arrTracks2.GetEntriesFast();
745 // remove leg1 contribution
746 for (Int_t itrack=0; itrack<ntrack1; ++itrack){
747 AliVTrack *track= static_cast<AliVTrack*>(arrTracks1.UncheckedAt(itrack));
5720c765 748 if (!track) continue;
4533e78e 749 // track contribution was already removed
750 if (mapRemovedTracks.FindObject(track)) continue;
751 else mapRemovedTracks.Add(track,track);
5720c765 752
4533e78e 753 cQX += (eptask->GetWeight(track) * TMath::Cos(2*track->Phi()));
754 cQY += (eptask->GetWeight(track) * TMath::Sin(2*track->Phi()));
5720c765 755 }
4533e78e 756 // remove leg2 contribution
757 for (Int_t itrack=0; itrack<ntrack2; ++itrack){
758 AliVTrack *track= static_cast<AliVTrack*>(arrTracks2.UncheckedAt(itrack));
5720c765 759 if (!track) continue;
4533e78e 760 // track contribution was already removed
761 if (mapRemovedTracks.FindObject(track)) continue;
762 else mapRemovedTracks.Add(track,track);
763
764 cQX += (eptask->GetWeight(track) * TMath::Cos(2*track->Phi()));
765 cQY += (eptask->GetWeight(track) * TMath::Sin(2*track->Phi()));
5720c765 766 }
5720c765 767
4533e78e 768 // build a corrected alieventplane using the values from the var manager
769 // these uncorrected values are filled using the stored magnitude and angle in the header
770 TVector2 qcorr;
771 qcorr.Set(AliDielectronVarManager::GetValue(AliDielectronVarManager::kTPCxH2uc)-cQX,
772 AliDielectronVarManager::GetValue(AliDielectronVarManager::kTPCyH2uc)-cQY);
773 // fill alieventplane
774 AliEventplane cevplane;
775 cevplane.SetQVector(&qcorr);
776 AliDielectronVarManager::SetTPCEventPlane(&cevplane);
777 cevplane.SetQVector(0);
778 return;
779 } //end: nanoAODs
780 else
781 {
782 // this is done in case of ESDs or AODs
783 Bool_t isESD=(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()->IsA()==AliESDInputHandler::Class());
784 // copy event plane object
785 AliEventplane cevplane(*evplane);
786 // Int_t nMaxID=cevplane->GetQContributionXArray()->GetSize();
787
788 TVector2 *qcorr = cevplane.GetQVector();
789 if(!qcorr) return;
790 TVector2 *qcsub1 = 0x0;
791 TVector2 *qcsub2 = 0x0;
792
793 // eta gap ?
794 Bool_t etagap = kFALSE;
795 for (Int_t iCut=0; iCut<fEventPlanePreFilter.GetCuts()->GetEntries();++iCut) {
796 TString cutName=fEventPlanePreFilter.GetCuts()->At(iCut)->GetName();
797 if(cutName.Contains("eta") || cutName.Contains("Eta")) etagap=kTRUE;
61d106d3 798 }
4533e78e 799
800 // subevent configuration for eta gap or LS (default is rndm)
801 if(fLikeSignSubEvents && etagap) {
802 // start with the full Qvector/event in both sub events
803 qcsub1 = new TVector2(*qcorr);
804 qcsub2 = new TVector2(*qcorr);
805 cevplane.SetQsub(qcsub1,qcsub2);
806
807 Int_t ntracks=ev->GetNumberOfTracks();
808 // track removals
809 for (Int_t itrack=0; itrack<ntracks; ++itrack){
810 AliVParticle *particle=ev->GetTrack(itrack);
811 AliVTrack *track= static_cast<AliVTrack*>(particle);
812 if (!track) continue;
813 if (track->GetID()>=0 && !isESD) continue;
814 Int_t tmpID = isESD ? track->GetID() : track->GetID()*-1 - 1;
815
816 // set contributions to zero
817 // charge sub1+ sub2-
818 if(fLikeSignSubEvents) {
819 Short_t charge=track->Charge();
820 if (charge<0) {
821 cevplane.GetQContributionXArraysub1()->SetAt(0.0, tmpID);
822 cevplane.GetQContributionYArraysub1()->SetAt(0.0, tmpID);
823 }
824 if (charge>0) {
825 cevplane.GetQContributionXArraysub2()->SetAt(0.0, tmpID);
826 cevplane.GetQContributionYArraysub2()->SetAt(0.0, tmpID);
827 }
828 }
829 // eta sub1+ sub2-
830 if(etagap) {
831 Double_t eta=track->Eta();
832 if (eta<0.0) {
833 cevplane.GetQContributionXArraysub1()->SetAt(0.0, tmpID);
834 cevplane.GetQContributionYArraysub1()->SetAt(0.0, tmpID);
835 }
836 if (eta>0.0) {
837 cevplane.GetQContributionXArraysub2()->SetAt(0.0, tmpID);
838 cevplane.GetQContributionYArraysub2()->SetAt(0.0, tmpID);
839 }
840 }
841 } // end: loop over tracks
842 } // end: sub event configuration
843
844 // apply cuts, e.g. etagap
845 if(fEventPlanePreFilter.GetCuts()->GetEntries()) {
846 UInt_t selectedMask=(1<<fEventPlanePreFilter.GetCuts()->GetEntries())-1;
847 Int_t ntracks=ev->GetNumberOfTracks();
848 for (Int_t itrack=0; itrack<ntracks; ++itrack){
849 AliVParticle *particle=ev->GetTrack(itrack);
850 AliVTrack *track= static_cast<AliVTrack*>(particle);
851 if (!track) continue;
852 if (track->GetID()>=0 && !isESD) continue;
853 Int_t tmpID = isESD ? track->GetID() : track->GetID()*-1 - 1;
854
855 //event plane cuts
856 UInt_t cutMask=fEventPlanePreFilter.IsSelected(track);
857 //apply cut
858 if (cutMask==selectedMask) continue;
859
860 // set contributions to zero
861 cevplane.GetQContributionXArray()->SetAt(0.0, tmpID);
862 cevplane.GetQContributionYArray()->SetAt(0.0, tmpID);
863 cevplane.GetQContributionXArraysub1()->SetAt(0.0, tmpID);
864 cevplane.GetQContributionYArraysub1()->SetAt(0.0, tmpID);
865 cevplane.GetQContributionXArraysub2()->SetAt(0.0, tmpID);
866 cevplane.GetQContributionYArraysub2()->SetAt(0.0, tmpID);
867 }
868 } // end: track cuts
869
870 // POI (particle of interest) rejection
871 Int_t pairIndex=GetPairIndex(arr1,arr2);
872 Int_t ntrack1=arrTracks1.GetEntriesFast();
873 Int_t ntrack2=arrTracks2.GetEntriesFast();
874 AliDielectronPair candidate;
875 candidate.SetKFUsage(fUseKF);
876
877 UInt_t selectedMask=(1<<fEventPlanePOIPreFilter.GetCuts()->GetEntries())-1;
878 for (Int_t itrack1=0; itrack1<ntrack1; ++itrack1){
879 Int_t end=ntrack2;
880 if (arr1==arr2) end=itrack1;
881 Bool_t accepted=kFALSE;
882 for (Int_t itrack2=0; itrack2<end; ++itrack2){
883 TObject *track1=arrTracks1.UncheckedAt(itrack1);
884 TObject *track2=arrTracks2.UncheckedAt(itrack2);
885 if (!track1 || !track2) continue;
886 //create the pair
887 candidate.SetTracks(static_cast<AliVTrack*>(track1), fPdgLeg1,
888 static_cast<AliVTrack*>(track2), fPdgLeg2);
889
890 candidate.SetType(pairIndex);
891 candidate.SetLabel(AliDielectronMC::Instance()->GetLabelMotherWithPdg(&candidate,fPdgMother));
892
893 //event plane cuts
894 UInt_t cutMask=fEventPlanePOIPreFilter.IsSelected(&candidate);
895 //apply cut
896 if (cutMask==selectedMask) continue;
897
898 accepted=kTRUE;
899 //remove the tracks from the Track arrays
900 arrTracks2.AddAt(0x0,itrack2);
901 }
5720c765 902 if ( accepted ) arrTracks1.AddAt(0x0,itrack1);
4533e78e 903 }
904 //compress the track arrays
905 arrTracks1.Compress();
906 arrTracks2.Compress();
907
908 //Modify the components: subtract the tracks
909 ntrack1=arrTracks1.GetEntriesFast();
910 ntrack2=arrTracks2.GetEntriesFast();
911 // remove leg1 contribution
912 for (Int_t itrack=0; itrack<ntrack1; ++itrack){
913 AliVTrack *track= static_cast<AliVTrack*>(arrTracks1.UncheckedAt(itrack));
914 if (!track) continue;
915 if (track->GetID()>=0 && !isESD) continue;
916 Int_t tmpID = isESD ? track->GetID() : track->GetID()*-1 - 1;
917 // set contributions to zero
918 cevplane.GetQContributionXArray()->SetAt(0.0, tmpID);
919 cevplane.GetQContributionYArray()->SetAt(0.0, tmpID);
920 cevplane.GetQContributionXArraysub1()->SetAt(0.0, tmpID);
921 cevplane.GetQContributionYArraysub1()->SetAt(0.0, tmpID);
922 cevplane.GetQContributionXArraysub2()->SetAt(0.0, tmpID);
923 cevplane.GetQContributionYArraysub2()->SetAt(0.0, tmpID);
924 }
925 // remove leg2 contribution
926 for (Int_t itrack=0; itrack<ntrack2; ++itrack){
927 AliVTrack *track= static_cast<AliVTrack*>(arrTracks2.UncheckedAt(itrack));
928 if (!track) continue;
929 if (track->GetID()>=0 && !isESD) continue;
930 Int_t tmpID = isESD ? track->GetID() : track->GetID()*-1 - 1;
931 // set contributions to zero
932 cevplane.GetQContributionXArray()->SetAt(0.0, tmpID);
933 cevplane.GetQContributionYArray()->SetAt(0.0, tmpID);
934 cevplane.GetQContributionXArraysub1()->SetAt(0.0, tmpID);
935 cevplane.GetQContributionYArraysub1()->SetAt(0.0, tmpID);
936 cevplane.GetQContributionXArraysub2()->SetAt(0.0, tmpID);
937 cevplane.GetQContributionYArraysub2()->SetAt(0.0, tmpID);
938 }
554e40f8 939
4533e78e 940 // set corrected AliEventplane and fill variables with corrected values
941 AliDielectronVarManager::SetTPCEventPlane(&cevplane);
942 delete qcsub1;
943 delete qcsub2;
944 } // end: ESD or AOD case
554e40f8 945
4533e78e 946}
5720c765 947//________________________________________________________________
948void AliDielectron::PairPreFilter(Int_t arr1, Int_t arr2, TObjArray &arrTracks1, TObjArray &arrTracks2)
949{
950 //
951 // Prefilter tracks from pairs
952 // Needed for datlitz rejections
953 // remove all tracks from the Single track arrays that pass the cuts in this filter
954 //
955
956 Int_t ntrack1=arrTracks1.GetEntriesFast();
957 Int_t ntrack2=arrTracks2.GetEntriesFast();
958 AliDielectronPair candidate;
08b801a6 959 candidate.SetKFUsage(fUseKF);
5720c765 960 // flag arrays for track removal
961 Bool_t *bTracks1 = new Bool_t[ntrack1];
962 for (Int_t itrack1=0; itrack1<ntrack1; ++itrack1) bTracks1[itrack1]=kFALSE;
963 Bool_t *bTracks2 = new Bool_t[ntrack2];
964 for (Int_t itrack2=0; itrack2<ntrack2; ++itrack2) bTracks2[itrack2]=kFALSE;
965
966 UInt_t selectedMask=(1<<fPairPreFilter.GetCuts()->GetEntries())-1;
967 UInt_t selectedMaskPair=(1<<fPairFilter.GetCuts()->GetEntries())-1;
968
969 Int_t nRejPasses = 1; //for fPreFilterUnlikeOnly and no set flag
970 if (fPreFilterAllSigns) nRejPasses = 3;
971
972 for (Int_t iRP=0; iRP < nRejPasses; ++iRP) {
973 Int_t arr1RP=arr1, arr2RP=arr2;
974 TObjArray *arrTracks1RP=&arrTracks1;
975 TObjArray *arrTracks2RP=&arrTracks2;
976 Bool_t *bTracks1RP = bTracks1;
977 Bool_t *bTracks2RP = bTracks2;
978 switch (iRP) {
979 case 1: arr1RP=arr1;arr2RP=arr1;
980 arrTracks1RP=&arrTracks1;
981 arrTracks2RP=&arrTracks1;
982 bTracks1RP = bTracks1;
983 bTracks2RP = bTracks1;
984 break;
985 case 2: arr1RP=arr2;arr2RP=arr2;
986 arrTracks1RP=&arrTracks2;
987 arrTracks2RP=&arrTracks2;
988 bTracks1RP = bTracks2;
989 bTracks2RP = bTracks2;
990 break;
991 default: ;//nothing to do
992 }
993 Int_t ntrack1RP=(*arrTracks1RP).GetEntriesFast();
994 Int_t ntrack2RP=(*arrTracks2RP).GetEntriesFast();
995
996 Int_t pairIndex=GetPairIndex(arr1RP,arr2RP);
997
998 for (Int_t itrack1=0; itrack1<ntrack1RP; ++itrack1){
999 Int_t end=ntrack2RP;
1000 if (arr1RP==arr2RP) end=itrack1;
1001 for (Int_t itrack2=0; itrack2<end; ++itrack2){
1002 TObject *track1=(*arrTracks1RP).UncheckedAt(itrack1);
1003 TObject *track2=(*arrTracks2RP).UncheckedAt(itrack2);
1004 if (!track1 || !track2) continue;
1005 //create the pair
1006 candidate.SetTracks(static_cast<AliVTrack*>(track1), fPdgLeg1,
1007 static_cast<AliVTrack*>(track2), fPdgLeg2);
1008
1009 candidate.SetType(pairIndex);
1010 candidate.SetLabel(AliDielectronMC::Instance()->GetLabelMotherWithPdg(&candidate,fPdgMother));
1011 //relate to the production vertex
1012 // if (AliDielectronVarManager::GetKFVertex()) candidate.SetProductionVertex(*AliDielectronVarManager::GetKFVertex());
1013
1014 //pair cuts
1015 UInt_t cutMask=fPairPreFilter.IsSelected(&candidate);
1016
1017 //apply cut
1018 if (cutMask!=selectedMask) continue;
1019 if (fCfManagerPair) fCfManagerPair->Fill(selectedMaskPair+1 ,&candidate);
1020 if (fHistos) FillHistogramsPair(&candidate,kTRUE);
1021 //set flags for track removal
1022 bTracks1RP[itrack1]=kTRUE;
1023 bTracks2RP[itrack2]=kTRUE;
1024 }
1025 }
1026 }
1027
1028 //remove the tracks from the Track arrays
1029 for (Int_t itrack1=0; itrack1<ntrack1; ++itrack1){
1030 if(bTracks1[itrack1]) arrTracks1.AddAt(0x0, itrack1);
1031 }
1032 for (Int_t itrack2=0; itrack2<ntrack2; ++itrack2){
1033 if(bTracks2[itrack2]) arrTracks2.AddAt(0x0, itrack2);
1034 }
1035
1036 // clean up
1037 delete [] bTracks1;
1038 delete [] bTracks2;
1039
1040 //compress the track arrays
61d106d3 1041 arrTracks1.Compress();
1042 arrTracks2.Compress();
1043
1044 //apply leg cuts after the pre filter
1045 if ( fPairPreFilterLegs.GetCuts()->GetEntries()>0 ) {
1046 selectedMask=(1<<fPairPreFilterLegs.GetCuts()->GetEntries())-1;
1047 //loop over tracks from array 1
1048 for (Int_t itrack=0; itrack<arrTracks1.GetEntriesFast();++itrack){
1049 //test cuts
1050 UInt_t cutMask=fPairPreFilterLegs.IsSelected(arrTracks1.UncheckedAt(itrack));
554e40f8 1051
61d106d3 1052 //apply cut
d4619886 1053 if (cutMask!=selectedMask) arrTracks1.AddAt(0x0,itrack);
61d106d3 1054 }
1055 arrTracks1.Compress();
1056
1057 //in case of like sign don't loop over second array
1058 if (arr1==arr2) {
1059 arrTracks2=arrTracks1;
1060 } else {
1061
1062 //loop over tracks from array 2
1063 for (Int_t itrack=0; itrack<arrTracks2.GetEntriesFast();++itrack){
1064 //test cuts
1065 UInt_t cutMask=fPairPreFilterLegs.IsSelected(arrTracks2.UncheckedAt(itrack));
1066 //apply cut
1067 if (cutMask!=selectedMask) arrTracks2.AddAt(0x0,itrack);
1068 }
1069 arrTracks2.Compress();
1070
1071 }
1072 }
554e40f8 1073 //For unlike-sign monitor track-cuts:
ba15fdfb 1074 if (arr1!=arr2&&fHistos) {
554e40f8 1075 TObjArray *unlikesignArray[2] = {&arrTracks1,&arrTracks2};
1076 FillHistogramsTracks(unlikesignArray);
1077 }
61d106d3 1078}
1079
b2a297fa 1080//________________________________________________________________
2a14a7b1 1081void AliDielectron::FillPairArrays(Int_t arr1, Int_t arr2)
1082{
b2a297fa 1083 //
1084 // select pairs and fill pair candidate arrays
1085 //
61d106d3 1086
1087 TObjArray arrTracks1=fTracks[arr1];
1088 TObjArray arrTracks2=fTracks[arr2];
1089
1090 //process pre filter if set
ba15fdfb 1091 if ((!fPreFilterAllSigns) && (!fPreFilterUnlikeOnly) && ( fPairPreFilter.GetCuts()->GetEntries()>0 )) PairPreFilter(arr1, arr2, arrTracks1, arrTracks2);
61d106d3 1092
b2a297fa 1093 Int_t pairIndex=GetPairIndex(arr1,arr2);
1094
61d106d3 1095 Int_t ntrack1=arrTracks1.GetEntriesFast();
1096 Int_t ntrack2=arrTracks2.GetEntriesFast();
b2a297fa 1097
1098 AliDielectronPair *candidate=new AliDielectronPair;
08b801a6 1099 candidate->SetKFUsage(fUseKF);
b2a297fa 1100
1101 UInt_t selectedMask=(1<<fPairFilter.GetCuts()->GetEntries())-1;
1102
1103 for (Int_t itrack1=0; itrack1<ntrack1; ++itrack1){
1104 Int_t end=ntrack2;
1105 if (arr1==arr2) end=itrack1;
1106 for (Int_t itrack2=0; itrack2<end; ++itrack2){
8df8e382 1107 //create the pair
61d106d3 1108 candidate->SetTracks(static_cast<AliVTrack*>(arrTracks1.UncheckedAt(itrack1)), fPdgLeg1,
5720c765 1109 static_cast<AliVTrack*>(arrTracks2.UncheckedAt(itrack2)), fPdgLeg2);
b2a297fa 1110 candidate->SetType(pairIndex);
e4339752 1111 Int_t label=AliDielectronMC::Instance()->GetLabelMotherWithPdg(candidate,fPdgMother);
1112 candidate->SetLabel(label);
1113 if (label>-1) candidate->SetPdgCode(fPdgMother);
b2a297fa 1114
4fae8ef9 1115 // check for gamma kf particle
1116 label=AliDielectronMC::Instance()->GetLabelMotherWithPdg(candidate,22);
1117 if (label>-1) {
1118 candidate->SetGammaTracks(static_cast<AliVTrack*>(arrTracks1.UncheckedAt(itrack1)), fPdgLeg1,
1119 static_cast<AliVTrack*>(arrTracks2.UncheckedAt(itrack2)), fPdgLeg2);
1120 // should we set the pdgmothercode and the label
1121 }
1122
b2a297fa 1123 //pair cuts
1124 UInt_t cutMask=fPairFilter.IsSelected(candidate);
d4619886 1125
b2a297fa 1126 //CF manager for the pair
1127 if (fCfManagerPair) fCfManagerPair->Fill(cutMask,candidate);
d4619886 1128
d8bb1abe 1129 // cut qa
dd718b9c 1130 if(pairIndex==kEv1PM && fCutQA) {
d8bb1abe 1131 fQAmonitor->FillAll(candidate);
1132 fQAmonitor->Fill(cutMask,candidate);
1133 }
b2a297fa 1134
1135 //apply cut
1136 if (cutMask!=selectedMask) continue;
1137
d4619886 1138 //histogram array for the pair
1139 if (fHistoArray) fHistoArray->Fill(pairIndex,candidate);
1140
b2a297fa 1141 //add the candidate to the candidate array
1142 PairArray(pairIndex)->Add(candidate);
1143 //get a new candidate
1144 candidate=new AliDielectronPair;
08b801a6 1145 candidate->SetKFUsage(fUseKF);
b2a297fa 1146 }
1147 }
1148 //delete the surplus candidate
1149 delete candidate;
1150}
1151
2a14a7b1 1152//________________________________________________________________
1153void AliDielectron::FillPairArrayTR()
1154{
1155 //
1156 // select pairs and fill pair candidate arrays
1157 //
1158 UInt_t selectedMask=(1<<fPairFilter.GetCuts()->GetEntries())-1;
1159
1160 while ( fTrackRotator->NextCombination() ){
1161 AliDielectronPair candidate;
08b801a6 1162 candidate.SetKFUsage(fUseKF);
1201a1a9 1163 candidate.SetTracks(&fTrackRotator->GetKFTrackP(), &fTrackRotator->GetKFTrackN(),
1164 fTrackRotator->GetVTrackP(),fTrackRotator->GetVTrackN());
ffbede40 1165 candidate.SetType(kEv1PMRot);
2a14a7b1 1166
1167 //pair cuts
1168 UInt_t cutMask=fPairFilter.IsSelected(&candidate);
1169
1170 //CF manager for the pair
1171 if (fCfManagerPair) fCfManagerPair->Fill(cutMask,&candidate);
1172
1173 //apply cut
5720c765 1174 if (cutMask==selectedMask) {
d4619886 1175
1176 //histogram array for the pair
1177 if (fHistoArray) fHistoArray->Fill((Int_t)kEv1PMRot,&candidate);
1178
1179 if(fHistos) FillHistogramsPair(&candidate);
1180 if(fStoreRotatedPairs) PairArray(kEv1PMRot)->Add(new AliDielectronPair(candidate));
1181 }
2a14a7b1 1182 }
1183}
1184
572b0139 1185//________________________________________________________________
1186void AliDielectron::FillDebugTree()
1187{
1188 //
1189 // Fill Histogram information for tracks and pairs
1190 //
1191
1192 //Fill Debug tree
1193 for (Int_t i=0; i<10; ++i){
1194 Int_t ntracks=PairArray(i)->GetEntriesFast();
1195 for (Int_t ipair=0; ipair<ntracks; ++ipair){
1196 fDebugTree->Fill(static_cast<AliDielectronPair*>(PairArray(i)->UncheckedAt(ipair)));
1197 }
1198 }
1199}
b2a297fa 1200
572b0139 1201//________________________________________________________________
1202void AliDielectron::SaveDebugTree()
1203{
1204 //
1205 // delete the debug tree, this will also write the tree
1206 //
1207 if (fDebugTree) fDebugTree->DeleteStreamer();
1208}
b2a297fa 1209
ba15fdfb 1210
1211//__________________________________________________________________
1212void AliDielectron::AddSignalMC(AliDielectronSignalMC* signal) {
1213 //
1214 // Add an MC signal to the signals list
1215 //
1216 if(!fSignalsMC) {
1217 fSignalsMC = new TObjArray();
1218 fSignalsMC->SetOwner();
1219 }
1220 fSignalsMC->Add(signal);
1221}
833da63d 1222
1223//________________________________________________________________
1224void AliDielectron::FillMCHistograms(Int_t label1, Int_t label2, Int_t nSignal) {
1225 //
1226 // fill QA MC TRUTH histograms for pairs and legs of all added mc signals
1227 //
1228
1229 TString className,className2,className3;
1230 className.Form("Pair_%s_MCtruth",fSignalsMC->At(nSignal)->GetName());
1231 className2.Form("Track_Legs_%s_MCtruth",fSignalsMC->At(nSignal)->GetName());
1232 className3.Form("Track_%s_%s_MCtruth",fgkPairClassNames[1],fSignalsMC->At(nSignal)->GetName());
1233 Bool_t pairClass=fHistos->GetHistogramList()->FindObject(className.Data())!=0x0;
1234 Bool_t legClass=fHistos->GetHistogramList()->FindObject(className2.Data())!=0x0;
1235 Bool_t trkClass=fHistos->GetHistogramList()->FindObject(className3.Data())!=0x0;
bc3cacd4 1236 // printf("fill signal %d: pair %d legs %d trk %d \n",nSignal,pairClass,legClass,trkClass);
833da63d 1237 if(!pairClass && !legClass && !trkClass) return;
1238
1239 // printf("leg labels: %d-%d \n",label1,label2);
1240 AliVParticle* part1 = AliDielectronMC::Instance()->GetMCTrackFromMCEvent(label1);
1241 AliVParticle* part2 = AliDielectronMC::Instance()->GetMCTrackFromMCEvent(label2);
1242 if(!part1 && !part2) return;
1243 if(part1&&part2) {
1244 // fill only unlike sign (and only SE)
833da63d 1245 if(part1->Charge()*part2->Charge()>=0) return;
1246 }
1247
1248
1249 AliDielectronMC* dieMC = AliDielectronMC::Instance();
1250
1251 Int_t mLabel1 = dieMC->GetMothersLabel(label1); // should work for both ESD and AOD
1252 Int_t mLabel2 = dieMC->GetMothersLabel(label2);
1253
1254 // check the same mother option
1255 AliDielectronSignalMC* sigMC = (AliDielectronSignalMC*)fSignalsMC->At(nSignal);
1256 if(sigMC->GetMothersRelation()==AliDielectronSignalMC::kSame && mLabel1!=mLabel2) return;
1257 if(sigMC->GetMothersRelation()==AliDielectronSignalMC::kDifferent && mLabel1==mLabel2) return;
1258
bc3cacd4 1259 // fill event values
833da63d 1260 Double_t values[AliDielectronVarManager::kNMaxValues];
1261 AliDielectronVarManager::Fill(dieMC->GetMCEvent(), values); // get event informations
bc3cacd4 1262
1263 // fill the leg variables
1264 // printf("leg:%d trk:%d part1:%p part2:%p \n",legClass,trkClass,part1,part2);
833da63d 1265 if (legClass || trkClass) {
1266 if(part1) AliDielectronVarManager::Fill(part1,values);
1267 if(part1 && trkClass) fHistos->FillClass(className3, AliDielectronVarManager::kNMaxValues, values);
1268 if(part1 && part2 && legClass) fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
1269 if(part2) AliDielectronVarManager::Fill(part2,values);
1270 if(part2 && trkClass) fHistos->FillClass(className3, AliDielectronVarManager::kNMaxValues, values);
1271 if(part1 && part2 && legClass) fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
1272 }
1273
1274 //fill pair information
1275 if (pairClass && part1 && part2) {
1276 AliDielectronVarManager::FillVarMCParticle2(part1,part2,values);
1277 fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
1278 }
1279
1280}
1281
f0351f9c 1282//________________________________________________________________
1283void AliDielectron::FillMCHistograms(const AliVEvent *ev) {
1284 //
1285 // fill QA MC histograms for pairs and legs of all added mc signals
1286 //
1287
4c06df48 1288 if (!fSignalsMC) return;
833da63d 1289 TString className,className2,className3;
f0351f9c 1290 Double_t values[AliDielectronVarManager::kNMaxValues]={0.};
1291 AliDielectronVarManager::Fill(ev, values); // get event informations
1292 //loop over all added mc signals
1293 for(Int_t isig=0; isig<fSignalsMC->GetEntries(); isig++) {
1294
bc3cacd4 1295 //check if and what to fill
f0351f9c 1296 className.Form("Pair_%s",fSignalsMC->At(isig)->GetName());
1297 className2.Form("Track_Legs_%s",fSignalsMC->At(isig)->GetName());
536a997a 1298 className3.Form("Track_%s_%s",fgkPairClassNames[1],fSignalsMC->At(isig)->GetName()); // unlike sign, SE only
f0351f9c 1299 Bool_t pairClass=fHistos->GetHistogramList()->FindObject(className.Data())!=0x0;
1300 Bool_t legClass=fHistos->GetHistogramList()->FindObject(className2.Data())!=0x0;
536a997a 1301 Bool_t mergedtrkClass=fHistos->GetHistogramList()->FindObject(className3.Data())!=0x0;
1302 if(!pairClass && !legClass && !mergedtrkClass) continue;
1303
1304 // fill pair and/or their leg variables
1305 if(pairClass || legClass) {
1306 Int_t npairs=PairArray(AliDielectron::kEv1PM)->GetEntriesFast(); // only SE +-
1307 for (Int_t ipair=0; ipair<npairs; ++ipair){
1308 AliDielectronPair *pair=static_cast<AliDielectronPair*>(PairArray(AliDielectron::kEv1PM)->UncheckedAt(ipair));
1309
1310 Bool_t isMCtruth = AliDielectronMC::Instance()->IsMCTruth(pair, (AliDielectronSignalMC*)fSignalsMC->At(isig));
1311 if(isMCtruth) {
1312 //fill pair information
1313 if (pairClass){
1314 AliDielectronVarManager::Fill(pair, values);
1315 fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
1316 }
1317 //fill leg information, both + and - in the same histo
1318 if (legClass){
1319 AliDielectronVarManager::Fill(pair->GetFirstDaughter(),values);
1320 fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
1321 AliDielectronVarManager::Fill(pair->GetSecondDaughter(),values);
1322 fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
1323 }
1324 } //is signal
1325 } //loop: pairs
1326 }
833da63d 1327
1328 // fill single tracks of signals
bc3cacd4 1329 if(!mergedtrkClass) continue;
1330 // loop over SE track arrays
1331 for (Int_t i=0; i<2; ++i){
1332 Int_t ntracks=fTracks[i].GetEntriesFast();
1333 for (Int_t itrack=0; itrack<ntracks; ++itrack){
1334 Int_t label=((AliVParticle*)fTracks[i].UncheckedAt(itrack))->GetLabel();
1335 Bool_t isMCtruth1 = AliDielectronMC::Instance()->IsMCTruth(label, (AliDielectronSignalMC*)fSignalsMC->At(isig), 1);
1336 Bool_t isMCtruth2 = AliDielectronMC::Instance()->IsMCTruth(label, (AliDielectronSignalMC*)fSignalsMC->At(isig), 2);
1337 // skip if track does not correspond to the signal
1338 if(!isMCtruth1 && !isMCtruth2) continue;
1339 AliDielectronVarManager::Fill(fTracks[i].UncheckedAt(itrack), values);
1340 fHistos->FillClass(className3, AliDielectronVarManager::kNMaxValues, values);
536a997a 1341 } //loop: tracks
bc3cacd4 1342 } //loop: arrays
833da63d 1343
f0351f9c 1344 } //loop: MCsignals
1345
1346}