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