]> 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
-its trk cut (julian)
[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
5720c765 52#include <AliEventplane.h>
b2a297fa 53#include <AliVEvent.h>
54#include <AliVParticle.h>
55#include <AliVTrack.h>
56#include "AliDielectronPair.h"
57#include "AliDielectronHistos.h"
58#include "AliDielectronCF.h"
59#include "AliDielectronMC.h"
60#include "AliDielectronVarManager.h"
2a14a7b1 61#include "AliDielectronTrackRotator.h"
572b0139 62#include "AliDielectronDebugTree.h"
ba15fdfb 63#include "AliDielectronSignalMC.h"
5720c765 64#include "AliDielectronMixingHandler.h"
572b0139 65
b2a297fa 66#include "AliDielectron.h"
67
68ClassImp(AliDielectron)
69
70const char* AliDielectron::fgkTrackClassNames[4] = {
71 "ev1+",
72 "ev1-",
73 "ev2+",
74 "ev2-"
75};
76
2a14a7b1 77const char* AliDielectron::fgkPairClassNames[11] = {
b2a297fa 78 "ev1+_ev1+",
79 "ev1+_ev1-",
80 "ev1-_ev1-",
81 "ev1+_ev2+",
82 "ev1-_ev2+",
83 "ev2+_ev2+",
84 "ev1+_ev2-",
85 "ev1-_ev2-",
86 "ev2+_ev2-",
2a14a7b1 87 "ev2-_ev2-",
88 "ev1+_ev1-_TR"
b2a297fa 89};
90
91//________________________________________________________________
92AliDielectron::AliDielectron() :
93 TNamed("AliDielectron","AliDielectron"),
94 fEventFilter("EventFilter"),
95 fTrackFilter("TrackFilter"),
61d106d3 96 fPairPreFilter("PairPreFilter"),
97 fPairPreFilterLegs("PairPreFilterLegs"),
b2a297fa 98 fPairFilter("PairFilter"),
5720c765 99 fEventPlanePreFilter("EventPlanePreFilter"),
100 fEventPlanePOIPreFilter("EventPlanePOIPreFilter"),
b2a297fa 101 fPdgMother(443),
8df8e382 102 fPdgLeg1(11),
103 fPdgLeg2(11),
ba15fdfb 104 fSignalsMC(0x0),
554e40f8 105 fNoPairing(kFALSE),
08b801a6 106 fUseKF(kTRUE),
b2a297fa 107 fHistos(0x0),
5720c765 108 fPairCandidates(new TObjArray(11)),
572b0139 109 fCfManagerPair(0x0),
2a14a7b1 110 fTrackRotator(0x0),
554e40f8 111 fDebugTree(0x0),
5720c765 112 fMixing(0x0),
113 fPreFilterEventPlane(kFALSE),
114 fLikeSignSubEvents(kFALSE),
fb7d2d99 115 fPreFilterUnlikeOnly(kFALSE),
ba15fdfb 116 fPreFilterAllSigns(kFALSE),
5720c765 117 fHasMC(kFALSE),
118 fStoreRotatedPairs(kFALSE),
2e02dba4 119 fDontClearArrays(kFALSE),
5720c765 120 fEstimatorFilename(""),
a823f01b 121 fTRDpidCorrectionFilename(""),
122 fVZEROCalibrationFilename(""),
123 fVZERORecenteringFilename("")
b2a297fa 124{
125 //
126 // Default constructor
127 //
128
129}
130
131//________________________________________________________________
132AliDielectron::AliDielectron(const char* name, const char* title) :
133 TNamed(name,title),
134 fEventFilter("EventFilter"),
135 fTrackFilter("TrackFilter"),
61d106d3 136 fPairPreFilter("PairPreFilter"),
137 fPairPreFilterLegs("PairPreFilterLegs"),
b2a297fa 138 fPairFilter("PairFilter"),
5720c765 139 fEventPlanePreFilter("EventPlanePreFilter"),
140 fEventPlanePOIPreFilter("EventPlanePOIPreFilter"),
b2a297fa 141 fPdgMother(443),
8df8e382 142 fPdgLeg1(11),
143 fPdgLeg2(11),
ba15fdfb 144 fSignalsMC(0x0),
554e40f8 145 fNoPairing(kFALSE),
08b801a6 146 fUseKF(kTRUE),
b2a297fa 147 fHistos(0x0),
5720c765 148 fPairCandidates(new TObjArray(11)),
572b0139 149 fCfManagerPair(0x0),
2a14a7b1 150 fTrackRotator(0x0),
554e40f8 151 fDebugTree(0x0),
5720c765 152 fMixing(0x0),
153 fPreFilterEventPlane(kFALSE),
154 fLikeSignSubEvents(kFALSE),
fb7d2d99 155 fPreFilterUnlikeOnly(kFALSE),
ba15fdfb 156 fPreFilterAllSigns(kFALSE),
5720c765 157 fHasMC(kFALSE),
158 fStoreRotatedPairs(kFALSE),
2e02dba4 159 fDontClearArrays(kFALSE),
5720c765 160 fEstimatorFilename(""),
a823f01b 161 fTRDpidCorrectionFilename(""),
162 fVZEROCalibrationFilename(""),
163 fVZERORecenteringFilename("")
b2a297fa 164{
165 //
166 // Named constructor
167 //
168
169}
170
171//________________________________________________________________
172AliDielectron::~AliDielectron()
173{
174 //
175 // Default destructor
176 //
177 if (fHistos) delete fHistos;
178 if (fPairCandidates) delete fPairCandidates;
572b0139 179 if (fDebugTree) delete fDebugTree;
5720c765 180 if (fMixing) delete fMixing;
ba15fdfb 181 if (fSignalsMC) delete fSignalsMC;
5720c765 182 if (fCfManagerPair) delete fCfManagerPair;
b2a297fa 183}
184
185//________________________________________________________________
186void AliDielectron::Init()
187{
188 //
189 // Initialise objects
190 //
fb7d2d99 191
192 if(GetHasMC()) AliDielectronMC::Instance()->SetHasMC(GetHasMC());
5720c765 193
194 InitPairCandidateArrays();
195
ba15fdfb 196 if (fCfManagerPair) {
197 fCfManagerPair->SetSignalsMC(fSignalsMC);
198 fCfManagerPair->InitialiseContainer(fPairFilter);
199 }
1201a1a9 200 if (fTrackRotator) {
201 fTrackRotator->SetTrackArrays(&fTracks[0],&fTracks[1]);
202 fTrackRotator->SetPdgLegs(fPdgLeg1,fPdgLeg2);
203 }
8df8e382 204 if (fDebugTree) fDebugTree->SetDielectron(this);
5720c765 205 if(fEstimatorFilename.Contains(".root")) AliDielectronVarManager::InitEstimatorAvg(fEstimatorFilename.Data());
206 if(fTRDpidCorrectionFilename.Contains(".root")) AliDielectronVarManager::InitTRDpidEffHistograms(fTRDpidCorrectionFilename.Data());
a823f01b 207 if(fVZEROCalibrationFilename.Contains(".root")) AliDielectronVarManager::SetVZEROCalibrationFile(fVZEROCalibrationFilename.Data());
208 if(fVZERORecenteringFilename.Contains(".root")) AliDielectronVarManager::SetVZERORecenteringFile(fVZERORecenteringFilename.Data());
209
5720c765 210 if (fMixing) fMixing->Init(this);
572b0139 211}
b2a297fa 212
213//________________________________________________________________
214void AliDielectron::Process(AliVEvent *ev1, AliVEvent *ev2)
215{
216 //
217 // Process the events
218 //
219
45b2b1b8 220 //at least first event is needed!
221 if (!ev1){
222 AliError("At least first event must be set!");
223 return;
224 }
572b0139 225 AliDielectronVarManager::SetEvent(ev1);
d327d9cd 226 if (fMixing){
227 //set mixing bin to event data
228 Int_t bin=fMixing->FindBin(AliDielectronVarManager::GetData());
229 AliDielectronVarManager::SetValue(AliDielectronVarManager::kMixingBin,bin);
230 }
5720c765 231
b2a297fa 232 //in case we have MC load the MC event and process the MC particles
40875e45 233 if (AliDielectronMC::Instance()->ConnectMCEvent()){
5720c765 234 ProcessMC(ev1);
572b0139 235 }
b2a297fa 236
237 //if candidate array doesn't exist, create it
238 if (!fPairCandidates->UncheckedAt(0)) {
239 InitPairCandidateArrays();
240 } else {
241 ClearArrays();
242 }
243
244 //mask used to require that all cuts are fulfilled
245 UInt_t selectedMask=(1<<fEventFilter.GetCuts()->GetEntries())-1;
246
247 //apply event cuts
37e9382d 248 if ((ev1&&fEventFilter.IsSelected(ev1)!=selectedMask) ||
249 (ev2&&fEventFilter.IsSelected(ev2)!=selectedMask)) return;
6551594b 250
d327d9cd 251// AliDielectronVarManager::SetEvent(ev1); // why a second time???
5720c765 252
b2a297fa 253 //fill track arrays for the first event
554e40f8 254 if (ev1){
255 FillTrackArrays(ev1);
ba15fdfb 256 if (((fPreFilterAllSigns)||(fPreFilterUnlikeOnly)) && ( fPairPreFilter.GetCuts()->GetEntries()>0 )) PairPreFilter(0, 1, fTracks[0], fTracks[1]);
554e40f8 257 }
258
b2a297fa 259
260 //fill track arrays for the second event
554e40f8 261 if (ev2) {
262 FillTrackArrays(ev2,1);
ba15fdfb 263 if (((fPreFilterAllSigns)||(fPreFilterUnlikeOnly)) && ( fPairPreFilter.GetCuts()->GetEntries()>0 )) PairPreFilter(2, 3, fTracks[2], fTracks[3]);
554e40f8 264 }
b2a297fa 265
5720c765 266 // TPC event plane correction
267 AliEventplane *cevplane = new AliEventplane();
268 if (ev1 && fPreFilterEventPlane && ( fEventPlanePreFilter.GetCuts()->GetEntries()>0 || fEventPlanePOIPreFilter.GetCuts()->GetEntries()>0))
269 EventPlanePreFilter(0, 1, fTracks[0], fTracks[1], ev1, cevplane);
270
554e40f8 271 if (!fNoPairing){
272 // create pairs and fill pair candidate arrays
273 for (Int_t itrackArr1=0; itrackArr1<4; ++itrackArr1){
274 for (Int_t itrackArr2=itrackArr1; itrackArr2<4; ++itrackArr2){
275 FillPairArrays(itrackArr1, itrackArr2);
276 }
b2a297fa 277 }
6551594b 278
554e40f8 279 //track rotation
1201a1a9 280 if (fTrackRotator) {
281 fTrackRotator->SetEvent(ev1);
282 FillPairArrayTR();
283 }
554e40f8 284 }
572b0139 285
286 //fill debug tree if a manager is attached
287 if (fDebugTree) FillDebugTree();
5720c765 288
289 //process event mixing
290 if (fMixing) {
291 fMixing->Fill(ev1,this);
292// FillHistograms(0x0,kTRUE);
293 }
294
295 //in case there is a histogram manager, fill the QA histograms
296 if (fHistos) FillHistograms(ev1);
297
298 // clear arrays
2e02dba4 299 if (!fDontClearArrays) ClearArrays();
5720c765 300 AliDielectronVarManager::SetTPCEventPlane(0x0);
301 delete cevplane;
b2a297fa 302}
303
304//________________________________________________________________
5720c765 305void AliDielectron::ProcessMC(AliVEvent *ev1)
b2a297fa 306{
307 //
308 // Process the MC data
309 //
310
ba15fdfb 311 AliDielectronMC *dieMC=AliDielectronMC::Instance();
312
5720c765 313 if (fHistos) FillHistogramsMC(dieMC->GetMCEvent(), ev1);
ba15fdfb 314
315 if(!fSignalsMC) return;
b2a297fa 316 //loop over all MC data and Fill the CF container if it exist
317 if (!fCfManagerPair) return;
318 fCfManagerPair->SetPdgMother(fPdgMother);
ba15fdfb 319 if(!fCfManagerPair->GetStepForMCtruth()) return;
320
321 // signals to be studied
322 Int_t nSignals = fSignalsMC->GetEntries();
323
324 // initialize 2D arrays of labels for particles from each MC signal
325 Int_t** labels1; // labels for particles satisfying branch 1
326 Int_t** labels2; // labels for particles satisfying branch 2
327 Int_t** labels12; // labels for particles satisfying both branches
328 labels1 = new Int_t*[nSignals];
329 labels2 = new Int_t*[nSignals];
330 labels12 = new Int_t*[nSignals];
331 Int_t* indexes1=new Int_t[nSignals];
332 Int_t* indexes2=new Int_t[nSignals];
333 Int_t* indexes12=new Int_t[nSignals];
334 for(Int_t isig=0;isig<nSignals;++isig) {
335 *(labels1+isig) = new Int_t[dieMC->GetNMCTracks()];
336 *(labels2+isig) = new Int_t[dieMC->GetNMCTracks()];
337 *(labels12+isig) = new Int_t[dieMC->GetNMCTracks()];
338 for(Int_t ip=0; ip<dieMC->GetNMCTracks();++ip) {
339 labels1[isig][ip] = -1;
340 labels2[isig][ip] = -1;
341 labels12[isig][ip] = -1;
342 }
343 indexes1[isig]=0;
344 indexes2[isig]=0;
345 indexes12[isig]=0;
b2a297fa 346 }
ba15fdfb 347
348 Bool_t truth1=kFALSE;
349 Bool_t truth2=kFALSE;
350 // loop over the MC tracks
351 for(Int_t ipart=0; ipart<dieMC->GetNMCTracks(); ++ipart) {
352 for(Int_t isig=0; isig<nSignals; ++isig) { // loop over signals
353 // Proceed only if this signal is required in the pure MC step
354 // NOTE: Some signals can be satisfied by many particles and this leads to high
355 // computation times (e.g. secondary electrons from the GEANT transport). Be aware of this!!
356 if(!((AliDielectronSignalMC*)fSignalsMC->At(isig))->GetFillPureMCStep()) continue;
357
358 truth1 = dieMC->IsMCTruth(ipart, (AliDielectronSignalMC*)fSignalsMC->At(isig), 1);
359 truth2 = dieMC->IsMCTruth(ipart, (AliDielectronSignalMC*)fSignalsMC->At(isig), 2);
360
361 // particles satisfying both branches are treated separately to avoid double counting during pairing
362 if(truth1 && truth2) {
363 labels12[isig][indexes12[isig]] = ipart;
364 ++indexes12[isig];
365 }
366 else {
367 if(truth1) {
368 labels1[isig][indexes1[isig]] = ipart;
369 ++indexes1[isig];
370 }
371 if(truth2) {
372 labels2[isig][indexes2[isig]] = ipart;
373 ++indexes2[isig];
374 }
375 }
376 }
377 } // end loop over MC particles
378
379 // Do the pairing and fill the CF container with pure MC info
380 for(Int_t isig=0; isig<nSignals; ++isig) {
381 // mix the particles which satisfy only one of the signal branches
382 for(Int_t i1=0;i1<indexes1[isig];++i1) {
383 for(Int_t i2=0;i2<indexes2[isig];++i2) {
384 fCfManagerPair->FillMC(labels1[isig][i1], labels2[isig][i2], isig);
385 }
386 }
387 // mix the particles which satisfy both branches
388 for(Int_t i1=0;i1<indexes12[isig];++i1) {
389 for(Int_t i2=0; i2<i1; ++i2) {
390 fCfManagerPair->FillMC(labels12[isig][i1], labels12[isig][i2], isig);
391 }
392 }
393 } // end loop over signals
394
395 // release the memory
396 for(Int_t isig=0;isig<nSignals;++isig) {
397 delete [] *(labels1+isig);
398 delete [] *(labels2+isig);
399 delete [] *(labels12+isig);
400 }
401 delete [] labels1;
402 delete [] labels2;
403 delete [] labels12;
404 delete [] indexes1;
405 delete [] indexes2;
406 delete [] indexes12;
b2a297fa 407}
408
554e40f8 409//________________________________________________________________
410void AliDielectron::FillHistogramsTracks(TObjArray **tracks)
411{
412 //
413 // Fill Histogram information for tracks after prefilter
414 // ignore mixed events - for prefilter, only single tracks +/- are relevant
415 //
416
417 TString className,className2;
418 Double_t values[AliDielectronVarManager::kNMaxValues];
419
420 //Fill track information, separately for the track array candidates
421 for (Int_t i=0; i<2; ++i){
422 className.Form("Pre_%s",fgkTrackClassNames[i]);
423 if (!fHistos->GetHistogramList()->FindObject(className.Data())) continue;
424 Int_t ntracks=tracks[i]->GetEntriesFast();
425 for (Int_t itrack=0; itrack<ntracks; ++itrack){
426 AliDielectronVarManager::Fill(tracks[i]->UncheckedAt(itrack), values);
427 fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
428 }
429 }
430}
ba15fdfb 431
432
433//________________________________________________________________
5720c765 434void AliDielectron::FillHistogramsMC(const AliMCEvent *ev, AliVEvent *ev1)
ba15fdfb 435{
436 //
437 // Fill Histogram information for MCEvents
438 //
439
5720c765 440 Double_t values[AliDielectronVarManager::kNMaxValues]={0.};
ba15fdfb 441 // Fill event information
5720c765 442 AliDielectronVarManager::Fill(ev1, values); // ESD/AOD information
443 AliDielectronVarManager::Fill(ev, values); // MC truth info
ba15fdfb 444 if (fHistos->GetHistogramList()->FindObject("MCEvent"))
445 fHistos->FillClass("MCEvent", AliDielectronVarManager::kNMaxValues, values);
446}
447
448
b2a297fa 449//________________________________________________________________
5720c765 450void AliDielectron::FillHistograms(const AliVEvent *ev, Bool_t pairInfoOnly)
b2a297fa 451{
452 //
453 // Fill Histogram information for tracks and pairs
454 //
455
61d106d3 456 TString className,className2;
5720c765 457 Double_t values[AliDielectronVarManager::kNMaxValues]={0.};
6551594b 458 //Fill event information
d327d9cd 459 if (ev){ //TODO: Why not use GetData() ??? See below event plane stuff!!!
460 AliDielectronVarManager::Fill(ev, values); //data should already be stored in AliDielectronVarManager from SetEvent, does EV plane correction rely on this???
461 if (fMixing){
462 //set mixing bin to event data
463 Int_t bin=fMixing->FindBin(values);
464 values[AliDielectronVarManager::kMixingBin]=bin;
465 }
466
5720c765 467 if (fHistos->GetHistogramList()->FindObject("Event"))
d327d9cd 468// fHistos->FillClass("Event", AliDielectronVarManager::kNMaxValues, AliDielectronVarManager::GetData());
469 fHistos->FillClass("Event", AliDielectronVarManager::kNMaxValues, values); //check event plane stuff and replace with above...
5720c765 470 }
6551594b 471
b2a297fa 472 //Fill track information, separately for the track array candidates
5720c765 473 if (!pairInfoOnly){
474 for (Int_t i=0; i<4; ++i){
475 className.Form("Track_%s",fgkTrackClassNames[i]);
476 if (!fHistos->GetHistogramList()->FindObject(className.Data())) continue;
477 Int_t ntracks=fTracks[i].GetEntriesFast();
478 for (Int_t itrack=0; itrack<ntracks; ++itrack){
479 AliDielectronVarManager::Fill(fTracks[i].UncheckedAt(itrack), values);
480 fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
481 }
b2a297fa 482 }
483 }
484
61d106d3 485 //Fill Pair information, separately for all pair candidate arrays and the legs
486 TObjArray arrLegs(100);
b2a297fa 487 for (Int_t i=0; i<10; ++i){
488 className.Form("Pair_%s",fgkPairClassNames[i]);
61d106d3 489 className2.Form("Track_Legs_%s",fgkPairClassNames[i]);
490 Bool_t pairClass=fHistos->GetHistogramList()->FindObject(className.Data())!=0x0;
491 Bool_t legClass=fHistos->GetHistogramList()->FindObject(className2.Data())!=0x0;
492 if (!pairClass&&!legClass) continue;
b2a297fa 493 Int_t ntracks=PairArray(i)->GetEntriesFast();
494 for (Int_t ipair=0; ipair<ntracks; ++ipair){
61d106d3 495 AliDielectronPair *pair=static_cast<AliDielectronPair*>(PairArray(i)->UncheckedAt(ipair));
496
497 //fill pair information
498 if (pairClass){
499 AliDielectronVarManager::Fill(pair, values);
500 fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
501 }
502
503 //fill leg information, don't fill the information twice
504 if (legClass){
505 AliVParticle *d1=pair->GetFirstDaughter();
506 AliVParticle *d2=pair->GetSecondDaughter();
507 if (!arrLegs.FindObject(d1)){
508 AliDielectronVarManager::Fill(d1, values);
509 fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
510 arrLegs.Add(d1);
511 }
512 if (!arrLegs.FindObject(d2)){
513 AliDielectronVarManager::Fill(d2, values);
514 fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
515 arrLegs.Add(d2);
516 }
517 }
b2a297fa 518 }
61d106d3 519 if (legClass) arrLegs.Clear();
b2a297fa 520 }
521
522}
2a14a7b1 523//________________________________________________________________
ffbede40 524void AliDielectron::FillHistogramsPair(AliDielectronPair *pair,Bool_t fromPreFilter/*=kFALSE*/)
2a14a7b1 525{
526 //
527 // Fill Histogram information for pairs and the track in the pair
528 // NOTE: in this funtion the leg information may be filled multiple
529 // times. This funtion is used in the track rotation pairing
530 // and those legs are not saved!
ffbede40 531 //
2a14a7b1 532 TString className,className2;
533 Double_t values[AliDielectronVarManager::kNMaxValues];
534
535 //Fill Pair information, separately for all pair candidate arrays and the legs
536 TObjArray arrLegs(100);
537 const Int_t type=pair->GetType();
ffbede40 538 if (fromPreFilter) {
539 className.Form("RejPair_%s",fgkPairClassNames[type]);
540 className2.Form("RejTrack_%s",fgkPairClassNames[type]);
541 } else {
542 className.Form("Pair_%s",fgkPairClassNames[type]);
543 className2.Form("Track_Legs_%s",fgkPairClassNames[type]);
544 }
2a14a7b1 545
546 Bool_t pairClass=fHistos->GetHistogramList()->FindObject(className.Data())!=0x0;
547 Bool_t legClass=fHistos->GetHistogramList()->FindObject(className2.Data())!=0x0;
548
549 //fill pair information
550 if (pairClass){
551 AliDielectronVarManager::Fill(pair, values);
552 fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
553 }
554
555 if (legClass){
556 AliVParticle *d1=pair->GetFirstDaughter();
557 AliDielectronVarManager::Fill(d1, values);
558 fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
559
560 AliVParticle *d2=pair->GetSecondDaughter();
561 AliDielectronVarManager::Fill(d2, values);
562 fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
563 }
564}
b2a297fa 565
566//________________________________________________________________
567void AliDielectron::FillTrackArrays(AliVEvent * const ev, Int_t eventNr)
568{
569 //
570 // select tracks and fill track candidate arrays
571 // eventNr = 0: First event, use track arrays 0 and 1
572 // eventNr = 1: Second event, use track arrays 2 and 3
573 //
574
575 Int_t ntracks=ev->GetNumberOfTracks();
576 UInt_t selectedMask=(1<<fTrackFilter.GetCuts()->GetEntries())-1;
577 for (Int_t itrack=0; itrack<ntracks; ++itrack){
578 //get particle
579 AliVParticle *particle=ev->GetTrack(itrack);
67fd1119 580
b2a297fa 581 //apply track cuts
582 if (fTrackFilter.IsSelected(particle)!=selectedMask) continue;
5720c765 583
b2a297fa 584 //fill selected particle into the corresponding track arrays
585 Short_t charge=particle->Charge();
586 if (charge>0) fTracks[eventNr*2].Add(particle);
587 else if (charge<0) fTracks[eventNr*2+1].Add(particle);
588 }
589}
590
61d106d3 591//________________________________________________________________
5720c765 592void AliDielectron::EventPlanePreFilter(Int_t arr1, Int_t arr2, TObjArray arrTracks1, TObjArray arrTracks2, const AliVEvent *ev, AliEventplane *cevplane)
2a14a7b1 593{
61d106d3 594 //
5720c765 595 // Prefilter tracks and tracks from pairs
596 // Needed for rejection in the Q-Vector of the event plane
597 // remove contribution of all tracks to the Q-vector that are in invariant mass window
61d106d3 598 //
5720c765 599 AliEventplane *evplane = const_cast<AliVEvent *>(ev)->GetEventplane();
b9d223bb 600// AliEventplane *evplane = ev->GetEventplane();
5720c765 601 if(!evplane) return;
602
603 // do not change these vectors qref
b9d223bb 604 TVector2 * const qref = evplane->GetQVector();
5720c765 605 if(!qref) return;
606 // random subevents
607 TVector2 *qrsub1 = evplane->GetQsub1();
608 TVector2 *qrsub2 = evplane->GetQsub2();
609
610 // copy references
b9d223bb 611 TVector2 *qcorr = new TVector2(*qref);
612 TVector2 *qcsub1 = 0x0;
613 TVector2 *qcsub2 = 0x0;
614 // printf("qrsub1 %p %f \n",qrsub1,qrsub1->X());
615
616
617 // eta gap ?
618 Bool_t etagap = kFALSE;
619 for (Int_t iCut=0; iCut<fEventPlanePreFilter.GetCuts()->GetEntries();++iCut) {
620 TString cutName=fEventPlanePreFilter.GetCuts()->At(iCut)->GetName();
621 if(cutName.Contains("eta") || cutName.Contains("Eta")) etagap=kTRUE;
622 }
5720c765 623
b9d223bb 624 // subevent separation
625 if(fLikeSignSubEvents || etagap) {
626 qcsub1 = new TVector2(*qcorr);
627 qcsub2 = new TVector2(*qcorr);
5720c765 628
629 Int_t ntracks=ev->GetNumberOfTracks();
630
631 // track removals
632 for (Int_t itrack=0; itrack<ntracks; ++itrack){
633 AliVParticle *particle=ev->GetTrack(itrack);
634 AliVTrack *track= static_cast<AliVTrack*>(particle);
635 if (!track) continue;
636
637 Double_t cQX = evplane->GetQContributionX(track);
638 Double_t cQY = evplane->GetQContributionY(track);
639
b9d223bb 640 // by charge sub1+ sub2-
641 if(fLikeSignSubEvents) {
642 Short_t charge=track->Charge();
643 if (charge<0) qcsub1->Set(qcsub1->X()-cQX, qcsub1->Y()-cQY);
644 if (charge>0) qcsub2->Set(qcsub2->X()-cQX, qcsub2->Y()-cQY);
645 }
646 // by eta sub1+ sub2-
647 if(etagap) {
5720c765 648 Double_t eta=track->Eta();
b9d223bb 649 if (eta<0.0) qcsub1->Set(qcsub1->X()-cQX, qcsub1->Y()-cQY);
650 if (eta>0.0) qcsub2->Set(qcsub2->X()-cQX, qcsub2->Y()-cQY);
5720c765 651 }
652 }
653 }
b9d223bb 654 else {
655 // by a random
656 qcsub1 = new TVector2(*qrsub1);
657 qcsub2 = new TVector2(*qrsub2);
658 }
5720c765 659
660 // apply cuts, e.g. etagap
661 if(fEventPlanePreFilter.GetCuts()->GetEntries()) {
662 UInt_t selectedMask=(1<<fEventPlanePreFilter.GetCuts()->GetEntries())-1;
663 Int_t ntracks=ev->GetNumberOfTracks();
664 for (Int_t itrack=0; itrack<ntracks; ++itrack){
665 AliVParticle *particle=ev->GetTrack(itrack);
666 AliVTrack *track= static_cast<AliVTrack*>(particle);
667 if (!track) continue;
668
669 //event plane cuts
670 UInt_t cutMask=fEventPlanePreFilter.IsSelected(track);
671 //apply cut
672 if (cutMask==selectedMask) continue;
673
674 Double_t cQX = 0.0;
675 Double_t cQY = 0.0;
676 if(!etagap) {
677 cQX = evplane->GetQContributionX(track);
678 cQY = evplane->GetQContributionY(track);
679 }
680 Double_t cQXsub1 = evplane->GetQContributionXsub1(track);
681 Double_t cQYsub1 = evplane->GetQContributionYsub1(track);
682 Double_t cQXsub2 = evplane->GetQContributionXsub2(track);
683 Double_t cQYsub2 = evplane->GetQContributionYsub2(track);
684
685 // update Q vectors
b9d223bb 686 qcorr->Set(qcorr->X()-cQX, qcorr->Y()-cQY);
687 qcsub1->Set(qcsub1->X()-cQXsub1, qcsub1->Y()-cQYsub1);
688 qcsub2->Set(qcsub2->X()-cQXsub2, qcsub2->Y()-cQYsub2);
5720c765 689 }
690 }
691
5720c765 692 // POI (particle of interest) rejection
61d106d3 693 Int_t pairIndex=GetPairIndex(arr1,arr2);
694
695 Int_t ntrack1=arrTracks1.GetEntriesFast();
696 Int_t ntrack2=arrTracks2.GetEntriesFast();
61d106d3 697 AliDielectronPair candidate;
08b801a6 698 candidate.SetKFUsage(fUseKF);
61d106d3 699
5720c765 700 UInt_t selectedMask=(1<<fEventPlanePOIPreFilter.GetCuts()->GetEntries())-1;
61d106d3 701 for (Int_t itrack1=0; itrack1<ntrack1; ++itrack1){
702 Int_t end=ntrack2;
703 if (arr1==arr2) end=itrack1;
704 Bool_t accepted=kFALSE;
705 for (Int_t itrack2=0; itrack2<end; ++itrack2){
5720c765 706 TObject *track1=arrTracks1.UncheckedAt(itrack1);
707 TObject *track2=arrTracks2.UncheckedAt(itrack2);
61d106d3 708 if (!track1 || !track2) continue;
709 //create the pair
5720c765 710 candidate.SetTracks(static_cast<AliVTrack*>(track1), fPdgLeg1,
711 static_cast<AliVTrack*>(track2), fPdgLeg2);
712
61d106d3 713 candidate.SetType(pairIndex);
714 candidate.SetLabel(AliDielectronMC::Instance()->GetLabelMotherWithPdg(&candidate,fPdgMother));
61d106d3 715
5720c765 716 //event plane cuts
717 UInt_t cutMask=fEventPlanePOIPreFilter.IsSelected(&candidate);
61d106d3 718 //apply cut
5720c765 719 if (cutMask==selectedMask) continue;
720
61d106d3 721 accepted=kTRUE;
722 //remove the tracks from the Track arrays
723 arrTracks2.AddAt(0x0,itrack2);
61d106d3 724 }
5720c765 725 if ( accepted ) arrTracks1.AddAt(0x0,itrack1);
61d106d3 726 }
727 //compress the track arrays
5720c765 728 arrTracks1.Compress();
729 arrTracks2.Compress();
730
554e40f8 731
5720c765 732 //Modify the components: subtract the tracks
733 ntrack1=arrTracks1.GetEntriesFast();
734 ntrack2=arrTracks2.GetEntriesFast();
735
736 // remove leg1 contribution
737 for (Int_t itrack=0; itrack<ntrack1; ++itrack){
738 AliVTrack *track= static_cast<AliVTrack*>(arrTracks1.UncheckedAt(itrack));
739 if (!track) continue;
740
741 Double_t cQX = evplane->GetQContributionX(track);
742 Double_t cQY = evplane->GetQContributionY(track);
743 Double_t cQXsub1 = evplane->GetQContributionXsub1(track);
744 Double_t cQYsub1 = evplane->GetQContributionYsub1(track);
745 Double_t cQXsub2 = evplane->GetQContributionXsub2(track);
746 Double_t cQYsub2 = evplane->GetQContributionYsub2(track);
747
748 // update Q vectors
749 qcorr->Set(qcorr->X()-cQX, qcorr->Y()-cQY);
750 qcsub1->Set(qcsub1->X()-cQXsub1, qcsub1->Y()-cQYsub1);
751 qcsub2->Set(qcsub2->X()-cQXsub2, qcsub2->Y()-cQYsub2);
752 }
753 // remove leg2 contribution
754 for (Int_t itrack=0; itrack<ntrack2; ++itrack){
755 AliVTrack *track= static_cast<AliVTrack*>(arrTracks2.UncheckedAt(itrack));
756 if (!track) continue;
757
758 Double_t cQX = evplane->GetQContributionX(track);
759 Double_t cQY = evplane->GetQContributionY(track);
760 Double_t cQXsub1 = evplane->GetQContributionXsub1(track);
761 Double_t cQYsub1 = evplane->GetQContributionYsub1(track);
762 Double_t cQXsub2 = evplane->GetQContributionXsub2(track);
763 Double_t cQYsub2 = evplane->GetQContributionYsub2(track);
764
765 // update Q vectors
766 qcorr->Set(qcorr->X()-cQX, qcorr->Y()-cQY);
767 qcsub1->Set(qcsub1->X()-cQXsub1, qcsub1->Y()-cQYsub1);
768 qcsub2->Set(qcsub2->X()-cQXsub2, qcsub2->Y()-cQYsub2);
769 }
554e40f8 770
b9d223bb 771 // printf("qrsub1 %p %f \t qcsub1 %p %f \n",qrsub1,qrsub1->X(),qcsub1,qcsub1->X());
5720c765 772 // set AliEventplane with corrected values
773 cevplane->SetQVector(qcorr);
774 cevplane->SetQsub(qcsub1, qcsub2);
775 AliDielectronVarManager::SetTPCEventPlane(cevplane);
776}
554e40f8 777
5720c765 778//________________________________________________________________
779void AliDielectron::PairPreFilter(Int_t arr1, Int_t arr2, TObjArray &arrTracks1, TObjArray &arrTracks2)
780{
781 //
782 // Prefilter tracks from pairs
783 // Needed for datlitz rejections
784 // remove all tracks from the Single track arrays that pass the cuts in this filter
785 //
786
787 Int_t ntrack1=arrTracks1.GetEntriesFast();
788 Int_t ntrack2=arrTracks2.GetEntriesFast();
789 AliDielectronPair candidate;
08b801a6 790 candidate.SetKFUsage(fUseKF);
5720c765 791 // flag arrays for track removal
792 Bool_t *bTracks1 = new Bool_t[ntrack1];
793 for (Int_t itrack1=0; itrack1<ntrack1; ++itrack1) bTracks1[itrack1]=kFALSE;
794 Bool_t *bTracks2 = new Bool_t[ntrack2];
795 for (Int_t itrack2=0; itrack2<ntrack2; ++itrack2) bTracks2[itrack2]=kFALSE;
796
797 UInt_t selectedMask=(1<<fPairPreFilter.GetCuts()->GetEntries())-1;
798 UInt_t selectedMaskPair=(1<<fPairFilter.GetCuts()->GetEntries())-1;
799
800 Int_t nRejPasses = 1; //for fPreFilterUnlikeOnly and no set flag
801 if (fPreFilterAllSigns) nRejPasses = 3;
802
803 for (Int_t iRP=0; iRP < nRejPasses; ++iRP) {
804 Int_t arr1RP=arr1, arr2RP=arr2;
805 TObjArray *arrTracks1RP=&arrTracks1;
806 TObjArray *arrTracks2RP=&arrTracks2;
807 Bool_t *bTracks1RP = bTracks1;
808 Bool_t *bTracks2RP = bTracks2;
809 switch (iRP) {
810 case 1: arr1RP=arr1;arr2RP=arr1;
811 arrTracks1RP=&arrTracks1;
812 arrTracks2RP=&arrTracks1;
813 bTracks1RP = bTracks1;
814 bTracks2RP = bTracks1;
815 break;
816 case 2: arr1RP=arr2;arr2RP=arr2;
817 arrTracks1RP=&arrTracks2;
818 arrTracks2RP=&arrTracks2;
819 bTracks1RP = bTracks2;
820 bTracks2RP = bTracks2;
821 break;
822 default: ;//nothing to do
823 }
824 Int_t ntrack1RP=(*arrTracks1RP).GetEntriesFast();
825 Int_t ntrack2RP=(*arrTracks2RP).GetEntriesFast();
826
827 Int_t pairIndex=GetPairIndex(arr1RP,arr2RP);
828
829 for (Int_t itrack1=0; itrack1<ntrack1RP; ++itrack1){
830 Int_t end=ntrack2RP;
831 if (arr1RP==arr2RP) end=itrack1;
832 for (Int_t itrack2=0; itrack2<end; ++itrack2){
833 TObject *track1=(*arrTracks1RP).UncheckedAt(itrack1);
834 TObject *track2=(*arrTracks2RP).UncheckedAt(itrack2);
835 if (!track1 || !track2) continue;
836 //create the pair
837 candidate.SetTracks(static_cast<AliVTrack*>(track1), fPdgLeg1,
838 static_cast<AliVTrack*>(track2), fPdgLeg2);
839
840 candidate.SetType(pairIndex);
841 candidate.SetLabel(AliDielectronMC::Instance()->GetLabelMotherWithPdg(&candidate,fPdgMother));
842 //relate to the production vertex
843 // if (AliDielectronVarManager::GetKFVertex()) candidate.SetProductionVertex(*AliDielectronVarManager::GetKFVertex());
844
845 //pair cuts
846 UInt_t cutMask=fPairPreFilter.IsSelected(&candidate);
847
848 //apply cut
849 if (cutMask!=selectedMask) continue;
850 if (fCfManagerPair) fCfManagerPair->Fill(selectedMaskPair+1 ,&candidate);
851 if (fHistos) FillHistogramsPair(&candidate,kTRUE);
852 //set flags for track removal
853 bTracks1RP[itrack1]=kTRUE;
854 bTracks2RP[itrack2]=kTRUE;
855 }
856 }
857 }
858
859 //remove the tracks from the Track arrays
860 for (Int_t itrack1=0; itrack1<ntrack1; ++itrack1){
861 if(bTracks1[itrack1]) arrTracks1.AddAt(0x0, itrack1);
862 }
863 for (Int_t itrack2=0; itrack2<ntrack2; ++itrack2){
864 if(bTracks2[itrack2]) arrTracks2.AddAt(0x0, itrack2);
865 }
866
867 // clean up
868 delete [] bTracks1;
869 delete [] bTracks2;
870
871 //compress the track arrays
61d106d3 872 arrTracks1.Compress();
873 arrTracks2.Compress();
874
875 //apply leg cuts after the pre filter
876 if ( fPairPreFilterLegs.GetCuts()->GetEntries()>0 ) {
877 selectedMask=(1<<fPairPreFilterLegs.GetCuts()->GetEntries())-1;
878 //loop over tracks from array 1
879 for (Int_t itrack=0; itrack<arrTracks1.GetEntriesFast();++itrack){
880 //test cuts
881 UInt_t cutMask=fPairPreFilterLegs.IsSelected(arrTracks1.UncheckedAt(itrack));
554e40f8 882
61d106d3 883 //apply cut
884 if (cutMask!=selectedMask) arrTracks1.AddAt(0x0,itrack);;
885 }
886 arrTracks1.Compress();
887
888 //in case of like sign don't loop over second array
889 if (arr1==arr2) {
890 arrTracks2=arrTracks1;
891 } else {
892
893 //loop over tracks from array 2
894 for (Int_t itrack=0; itrack<arrTracks2.GetEntriesFast();++itrack){
895 //test cuts
896 UInt_t cutMask=fPairPreFilterLegs.IsSelected(arrTracks2.UncheckedAt(itrack));
897 //apply cut
898 if (cutMask!=selectedMask) arrTracks2.AddAt(0x0,itrack);
899 }
900 arrTracks2.Compress();
901
902 }
903 }
554e40f8 904 //For unlike-sign monitor track-cuts:
ba15fdfb 905 if (arr1!=arr2&&fHistos) {
554e40f8 906 TObjArray *unlikesignArray[2] = {&arrTracks1,&arrTracks2};
907 FillHistogramsTracks(unlikesignArray);
908 }
61d106d3 909}
910
b2a297fa 911//________________________________________________________________
2a14a7b1 912void AliDielectron::FillPairArrays(Int_t arr1, Int_t arr2)
913{
b2a297fa 914 //
915 // select pairs and fill pair candidate arrays
916 //
61d106d3 917
918 TObjArray arrTracks1=fTracks[arr1];
919 TObjArray arrTracks2=fTracks[arr2];
920
921 //process pre filter if set
ba15fdfb 922 if ((!fPreFilterAllSigns) && (!fPreFilterUnlikeOnly) && ( fPairPreFilter.GetCuts()->GetEntries()>0 )) PairPreFilter(arr1, arr2, arrTracks1, arrTracks2);
61d106d3 923
b2a297fa 924 Int_t pairIndex=GetPairIndex(arr1,arr2);
925
61d106d3 926 Int_t ntrack1=arrTracks1.GetEntriesFast();
927 Int_t ntrack2=arrTracks2.GetEntriesFast();
b2a297fa 928
929 AliDielectronPair *candidate=new AliDielectronPair;
08b801a6 930 candidate->SetKFUsage(fUseKF);
b2a297fa 931
932 UInt_t selectedMask=(1<<fPairFilter.GetCuts()->GetEntries())-1;
933
934 for (Int_t itrack1=0; itrack1<ntrack1; ++itrack1){
935 Int_t end=ntrack2;
936 if (arr1==arr2) end=itrack1;
937 for (Int_t itrack2=0; itrack2<end; ++itrack2){
8df8e382 938 //create the pair
61d106d3 939 candidate->SetTracks(static_cast<AliVTrack*>(arrTracks1.UncheckedAt(itrack1)), fPdgLeg1,
5720c765 940 static_cast<AliVTrack*>(arrTracks2.UncheckedAt(itrack2)), fPdgLeg2);
b2a297fa 941 candidate->SetType(pairIndex);
e4339752 942 Int_t label=AliDielectronMC::Instance()->GetLabelMotherWithPdg(candidate,fPdgMother);
943 candidate->SetLabel(label);
944 if (label>-1) candidate->SetPdgCode(fPdgMother);
b2a297fa 945
946 //pair cuts
947 UInt_t cutMask=fPairFilter.IsSelected(candidate);
948
949 //CF manager for the pair
950 if (fCfManagerPair) fCfManagerPair->Fill(cutMask,candidate);
951
952 //apply cut
953 if (cutMask!=selectedMask) continue;
954
955 //add the candidate to the candidate array
956 PairArray(pairIndex)->Add(candidate);
957 //get a new candidate
958 candidate=new AliDielectronPair;
08b801a6 959 candidate->SetKFUsage(fUseKF);
b2a297fa 960 }
961 }
962 //delete the surplus candidate
963 delete candidate;
964}
965
2a14a7b1 966//________________________________________________________________
967void AliDielectron::FillPairArrayTR()
968{
969 //
970 // select pairs and fill pair candidate arrays
971 //
972 UInt_t selectedMask=(1<<fPairFilter.GetCuts()->GetEntries())-1;
973
974 while ( fTrackRotator->NextCombination() ){
975 AliDielectronPair candidate;
08b801a6 976 candidate.SetKFUsage(fUseKF);
1201a1a9 977 candidate.SetTracks(&fTrackRotator->GetKFTrackP(), &fTrackRotator->GetKFTrackN(),
978 fTrackRotator->GetVTrackP(),fTrackRotator->GetVTrackN());
ffbede40 979 candidate.SetType(kEv1PMRot);
2a14a7b1 980
981 //pair cuts
982 UInt_t cutMask=fPairFilter.IsSelected(&candidate);
983
984 //CF manager for the pair
985 if (fCfManagerPair) fCfManagerPair->Fill(cutMask,&candidate);
986
987 //apply cut
5720c765 988 if (cutMask==selectedMask) {
989 if(fHistos) FillHistogramsPair(&candidate);
e4339752 990 if(fStoreRotatedPairs) PairArray(kEv1PMRot)->Add(new AliDielectronPair(candidate));
5720c765 991 }
2a14a7b1 992 }
993}
994
572b0139 995//________________________________________________________________
996void AliDielectron::FillDebugTree()
997{
998 //
999 // Fill Histogram information for tracks and pairs
1000 //
1001
1002 //Fill Debug tree
1003 for (Int_t i=0; i<10; ++i){
1004 Int_t ntracks=PairArray(i)->GetEntriesFast();
1005 for (Int_t ipair=0; ipair<ntracks; ++ipair){
1006 fDebugTree->Fill(static_cast<AliDielectronPair*>(PairArray(i)->UncheckedAt(ipair)));
1007 }
1008 }
1009}
b2a297fa 1010
572b0139 1011//________________________________________________________________
1012void AliDielectron::SaveDebugTree()
1013{
1014 //
1015 // delete the debug tree, this will also write the tree
1016 //
1017 if (fDebugTree) fDebugTree->DeleteStreamer();
1018}
b2a297fa 1019
ba15fdfb 1020
1021//__________________________________________________________________
1022void AliDielectron::AddSignalMC(AliDielectronSignalMC* signal) {
1023 //
1024 // Add an MC signal to the signals list
1025 //
1026 if(!fSignalsMC) {
1027 fSignalsMC = new TObjArray();
1028 fSignalsMC->SetOwner();
1029 }
1030 fSignalsMC->Add(signal);
1031}