o updates for AOD MC (Julian)
[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()->ConnectMCEvent()){
234     ProcessMC(ev1);
235   }
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
248     if ((ev1&&fEventFilter.IsSelected(ev1)!=selectedMask) ||
249         (ev2&&fEventFilter.IsSelected(ev2)!=selectedMask)) return;
250   
251 //   AliDielectronVarManager::SetEvent(ev1); // why a second time???
252
253   //fill track arrays for the first event
254   if (ev1){
255     FillTrackArrays(ev1);
256     if (((fPreFilterAllSigns)||(fPreFilterUnlikeOnly)) && ( fPairPreFilter.GetCuts()->GetEntries()>0 )) PairPreFilter(0, 1, fTracks[0], fTracks[1]);
257   }
258
259
260   //fill track arrays for the second event
261   if (ev2) {
262     FillTrackArrays(ev2,1);
263     if (((fPreFilterAllSigns)||(fPreFilterUnlikeOnly)) && ( fPairPreFilter.GetCuts()->GetEntries()>0 )) PairPreFilter(2, 3, fTracks[2], fTracks[3]);
264   }
265
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   
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       }
277     }
278
279     //track rotation
280     if (fTrackRotator) {
281       fTrackRotator->SetEvent(ev1);
282       FillPairArrayTR();
283     }
284   }
285
286   //fill debug tree if a manager is attached
287   if (fDebugTree) FillDebugTree();
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
299   if (!fDontClearArrays) ClearArrays();
300   AliDielectronVarManager::SetTPCEventPlane(0x0);
301   delete cevplane;
302 }
303
304 //________________________________________________________________
305 void AliDielectron::ProcessMC(AliVEvent *ev1)
306 {
307   //
308   // Process the MC data
309   //
310
311   AliDielectronMC *dieMC=AliDielectronMC::Instance();
312
313   if (fHistos) FillHistogramsMC(dieMC->GetMCEvent(), ev1);
314
315   if(!fSignalsMC) return;
316   //loop over all MC data and Fill the CF container if it exist
317   if (!fCfManagerPair) return;
318   fCfManagerPair->SetPdgMother(fPdgMother);
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;
346   }
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;
407 }
408
409 //________________________________________________________________
410 void 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 }
431
432
433 //________________________________________________________________
434 void AliDielectron::FillHistogramsMC(const AliMCEvent *ev, AliVEvent *ev1)
435 {
436   //
437   // Fill Histogram information for MCEvents
438   //
439
440   Double_t values[AliDielectronVarManager::kNMaxValues]={0.};
441   // Fill event information
442   AliDielectronVarManager::Fill(ev1, values);    // ESD/AOD information
443   AliDielectronVarManager::Fill(ev, values);     // MC truth info
444   if (fHistos->GetHistogramList()->FindObject("MCEvent"))
445     fHistos->FillClass("MCEvent", AliDielectronVarManager::kNMaxValues, values);
446 }
447
448
449 //________________________________________________________________
450 void AliDielectron::FillHistograms(const AliVEvent *ev, Bool_t pairInfoOnly)
451 {
452   //
453   // Fill Histogram information for tracks and pairs
454   //
455   
456   TString  className,className2;
457   Double_t values[AliDielectronVarManager::kNMaxValues]={0.};
458   //Fill event information
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
467     if (fHistos->GetHistogramList()->FindObject("Event"))
468 //       fHistos->FillClass("Event", AliDielectronVarManager::kNMaxValues, AliDielectronVarManager::GetData());
469       fHistos->FillClass("Event", AliDielectronVarManager::kNMaxValues, values); //check event plane stuff and replace with above...
470   }
471   
472   //Fill track information, separately for the track array candidates
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       }
482     }
483   }
484
485   //Fill Pair information, separately for all pair candidate arrays and the legs
486   TObjArray arrLegs(100);
487   for (Int_t i=0; i<10; ++i){
488     className.Form("Pair_%s",fgkPairClassNames[i]);
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;
493     Int_t ntracks=PairArray(i)->GetEntriesFast();
494     for (Int_t ipair=0; ipair<ntracks; ++ipair){
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       }
518     }
519     if (legClass) arrLegs.Clear();
520   }
521   
522 }
523 //________________________________________________________________
524 void AliDielectron::FillHistogramsPair(AliDielectronPair *pair,Bool_t fromPreFilter/*=kFALSE*/)
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!
531   //
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();
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   }
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 }
565
566 //________________________________________________________________
567 void 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);
580     //TODO: temporary solution, perhaps think about a better implementation
581     //      This is needed to use AliESDpidCuts, which relies on the ESD event
582     //      is set as a AliESDtrack attribute... somehow ugly!
583     if (ev->IsA()==AliESDEvent::Class()){
584       AliESDtrack *track=static_cast<AliESDtrack*>(particle);
585       track->SetESDEvent(static_cast<AliESDEvent*>(ev)); //only in trunk...
586     }
587     //apply track cuts
588     if (fTrackFilter.IsSelected(particle)!=selectedMask) continue;
589
590     //fill selected particle into the corresponding track arrays
591     Short_t charge=particle->Charge();
592     if (charge>0)      fTracks[eventNr*2].Add(particle);
593     else if (charge<0) fTracks[eventNr*2+1].Add(particle);
594   }
595 }
596
597 //________________________________________________________________
598 void AliDielectron::EventPlanePreFilter(Int_t arr1, Int_t arr2, TObjArray arrTracks1, TObjArray arrTracks2, const AliVEvent *ev, AliEventplane *cevplane)
599 {
600   //
601   // Prefilter tracks and tracks from pairs
602   // Needed for rejection in the Q-Vector of the event plane
603   // remove contribution of all tracks to the Q-vector that are in invariant mass window 
604   //
605   AliEventplane *evplane = const_cast<AliVEvent *>(ev)->GetEventplane();
606 //   AliEventplane *evplane = ev->GetEventplane();
607   if(!evplane) return;
608   
609   // do not change these vectors qref
610   TVector2 * const qref  = evplane->GetQVector();
611   if(!qref) return;
612   // random subevents
613   TVector2 *qrsub1 = evplane->GetQsub1();
614   TVector2 *qrsub2 = evplane->GetQsub2();
615
616   // copy references
617   TVector2 *qcorr  = new TVector2(*qref);
618   TVector2 *qcsub1 = 0x0;
619   TVector2 *qcsub2 = 0x0;
620   //  printf("qrsub1 %p %f \n",qrsub1,qrsub1->X());
621
622
623   // eta gap ?
624   Bool_t etagap = kFALSE;
625   for (Int_t iCut=0; iCut<fEventPlanePreFilter.GetCuts()->GetEntries();++iCut) {
626     TString cutName=fEventPlanePreFilter.GetCuts()->At(iCut)->GetName();
627     if(cutName.Contains("eta") || cutName.Contains("Eta"))  etagap=kTRUE;
628   }
629
630   // subevent separation
631   if(fLikeSignSubEvents || etagap) {
632     qcsub1 = new TVector2(*qcorr);
633     qcsub2 = new TVector2(*qcorr);
634
635     Int_t ntracks=ev->GetNumberOfTracks();
636     
637     // track removals
638     for (Int_t itrack=0; itrack<ntracks; ++itrack){
639       AliVParticle *particle=ev->GetTrack(itrack);
640       AliVTrack *track= static_cast<AliVTrack*>(particle);
641       if (!track) continue;
642
643       Double_t cQX     = evplane->GetQContributionX(track);
644       Double_t cQY     = evplane->GetQContributionY(track);
645       
646       // by charge sub1+ sub2-
647       if(fLikeSignSubEvents) {
648         Short_t charge=track->Charge();
649         if (charge<0) qcsub1->Set(qcsub1->X()-cQX, qcsub1->Y()-cQY);
650         if (charge>0) qcsub2->Set(qcsub2->X()-cQX, qcsub2->Y()-cQY);
651       }
652       // by eta sub1+ sub2-
653       if(etagap) {
654         Double_t eta=track->Eta();
655         if (eta<0.0) qcsub1->Set(qcsub1->X()-cQX, qcsub1->Y()-cQY);
656         if (eta>0.0) qcsub2->Set(qcsub2->X()-cQX, qcsub2->Y()-cQY);
657       }
658     }
659   }
660   else {
661     // by a random
662     qcsub1 = new TVector2(*qrsub1);
663     qcsub2 = new TVector2(*qrsub2);
664   }
665   
666   // apply cuts, e.g. etagap 
667   if(fEventPlanePreFilter.GetCuts()->GetEntries()) {
668     UInt_t selectedMask=(1<<fEventPlanePreFilter.GetCuts()->GetEntries())-1;
669     Int_t ntracks=ev->GetNumberOfTracks();
670     for (Int_t itrack=0; itrack<ntracks; ++itrack){
671       AliVParticle *particle=ev->GetTrack(itrack);
672       AliVTrack *track= static_cast<AliVTrack*>(particle);
673       if (!track) continue;
674       
675       //event plane cuts
676       UInt_t cutMask=fEventPlanePreFilter.IsSelected(track);
677       //apply cut
678       if (cutMask==selectedMask) continue;
679
680       Double_t cQX     = 0.0;
681       Double_t cQY     = 0.0;
682       if(!etagap) {
683         cQX = evplane->GetQContributionX(track);
684         cQY = evplane->GetQContributionY(track);
685       }
686       Double_t cQXsub1 = evplane->GetQContributionXsub1(track);
687       Double_t cQYsub1 = evplane->GetQContributionYsub1(track);
688       Double_t cQXsub2 = evplane->GetQContributionXsub2(track);
689       Double_t cQYsub2 = evplane->GetQContributionYsub2(track);      
690
691       // update Q vectors
692       qcorr->Set(qcorr->X()-cQX, qcorr->Y()-cQY);
693       qcsub1->Set(qcsub1->X()-cQXsub1, qcsub1->Y()-cQYsub1);
694       qcsub2->Set(qcsub2->X()-cQXsub2, qcsub2->Y()-cQYsub2);
695     }
696   }
697
698   // POI (particle of interest) rejection
699   Int_t pairIndex=GetPairIndex(arr1,arr2);
700   
701   Int_t ntrack1=arrTracks1.GetEntriesFast();
702   Int_t ntrack2=arrTracks2.GetEntriesFast();
703   AliDielectronPair candidate;
704   
705   UInt_t selectedMask=(1<<fEventPlanePOIPreFilter.GetCuts()->GetEntries())-1;
706   for (Int_t itrack1=0; itrack1<ntrack1; ++itrack1){
707     Int_t end=ntrack2;
708     if (arr1==arr2) end=itrack1;
709     Bool_t accepted=kFALSE;
710     for (Int_t itrack2=0; itrack2<end; ++itrack2){
711       TObject *track1=arrTracks1.UncheckedAt(itrack1);
712       TObject *track2=arrTracks2.UncheckedAt(itrack2);
713       if (!track1 || !track2) continue;
714       //create the pair
715       candidate.SetTracks(static_cast<AliVTrack*>(track1), fPdgLeg1,
716                           static_cast<AliVTrack*>(track2), fPdgLeg2);
717       
718       candidate.SetType(pairIndex);
719       candidate.SetLabel(AliDielectronMC::Instance()->GetLabelMotherWithPdg(&candidate,fPdgMother));
720       
721       //event plane cuts
722       UInt_t cutMask=fEventPlanePOIPreFilter.IsSelected(&candidate);
723       //apply cut
724       if (cutMask==selectedMask) continue;
725
726       accepted=kTRUE;
727       //remove the tracks from the Track arrays
728       arrTracks2.AddAt(0x0,itrack2);
729     }
730       if ( accepted ) arrTracks1.AddAt(0x0,itrack1);
731   }
732   //compress the track arrays
733   arrTracks1.Compress();
734   arrTracks2.Compress();
735   
736   
737   //Modify the components: subtract the tracks
738   ntrack1=arrTracks1.GetEntriesFast();
739   ntrack2=arrTracks2.GetEntriesFast();
740   
741   // remove leg1 contribution
742   for (Int_t itrack=0; itrack<ntrack1; ++itrack){
743     AliVTrack *track= static_cast<AliVTrack*>(arrTracks1.UncheckedAt(itrack));
744     if (!track) continue;
745     
746     Double_t cQX     = evplane->GetQContributionX(track);
747     Double_t cQY     = evplane->GetQContributionY(track);
748     Double_t cQXsub1 = evplane->GetQContributionXsub1(track);
749     Double_t cQYsub1 = evplane->GetQContributionYsub1(track);
750     Double_t cQXsub2 = evplane->GetQContributionXsub2(track);
751     Double_t cQYsub2 = evplane->GetQContributionYsub2(track);
752     
753     // update Q vectors
754     qcorr->Set(qcorr->X()-cQX, qcorr->Y()-cQY);
755     qcsub1->Set(qcsub1->X()-cQXsub1, qcsub1->Y()-cQYsub1);
756     qcsub2->Set(qcsub2->X()-cQXsub2, qcsub2->Y()-cQYsub2);
757   }
758   // remove leg2 contribution
759   for (Int_t itrack=0; itrack<ntrack2; ++itrack){
760     AliVTrack *track= static_cast<AliVTrack*>(arrTracks2.UncheckedAt(itrack));
761     if (!track) continue;
762     
763     Double_t cQX     = evplane->GetQContributionX(track);
764     Double_t cQY     = evplane->GetQContributionY(track);
765     Double_t cQXsub1 = evplane->GetQContributionXsub1(track);
766     Double_t cQYsub1 = evplane->GetQContributionYsub1(track);
767     Double_t cQXsub2 = evplane->GetQContributionXsub2(track);
768     Double_t cQYsub2 = evplane->GetQContributionYsub2(track);
769     
770     // update Q vectors
771     qcorr->Set(qcorr->X()-cQX, qcorr->Y()-cQY);
772     qcsub1->Set(qcsub1->X()-cQXsub1, qcsub1->Y()-cQYsub1);
773     qcsub2->Set(qcsub2->X()-cQXsub2, qcsub2->Y()-cQYsub2);
774   }
775
776   //  printf("qrsub1 %p %f \t qcsub1 %p %f \n",qrsub1,qrsub1->X(),qcsub1,qcsub1->X());
777   // set AliEventplane with corrected values
778   cevplane->SetQVector(qcorr);
779   cevplane->SetQsub(qcsub1, qcsub2);
780   AliDielectronVarManager::SetTPCEventPlane(cevplane);
781 }
782
783 //________________________________________________________________
784 void AliDielectron::PairPreFilter(Int_t arr1, Int_t arr2, TObjArray &arrTracks1, TObjArray &arrTracks2)
785 {
786   //
787   // Prefilter tracks from pairs
788   // Needed for datlitz rejections
789   // remove all tracks from the Single track arrays that pass the cuts in this filter
790   //
791
792   Int_t ntrack1=arrTracks1.GetEntriesFast();
793   Int_t ntrack2=arrTracks2.GetEntriesFast();
794   AliDielectronPair candidate;
795
796   // flag arrays for track removal
797   Bool_t *bTracks1 = new Bool_t[ntrack1];
798   for (Int_t itrack1=0; itrack1<ntrack1; ++itrack1) bTracks1[itrack1]=kFALSE;
799   Bool_t *bTracks2 = new Bool_t[ntrack2];
800   for (Int_t itrack2=0; itrack2<ntrack2; ++itrack2) bTracks2[itrack2]=kFALSE;
801
802   UInt_t selectedMask=(1<<fPairPreFilter.GetCuts()->GetEntries())-1;
803   UInt_t selectedMaskPair=(1<<fPairFilter.GetCuts()->GetEntries())-1;
804
805   Int_t nRejPasses = 1; //for fPreFilterUnlikeOnly and no set flag 
806   if (fPreFilterAllSigns) nRejPasses = 3;
807
808   for (Int_t iRP=0; iRP < nRejPasses; ++iRP) {
809         Int_t arr1RP=arr1, arr2RP=arr2;
810         TObjArray *arrTracks1RP=&arrTracks1;
811         TObjArray *arrTracks2RP=&arrTracks2;
812         Bool_t *bTracks1RP = bTracks1;
813         Bool_t *bTracks2RP = bTracks2;
814         switch (iRP) {
815                 case 1: arr1RP=arr1;arr2RP=arr1;
816                                 arrTracks1RP=&arrTracks1;
817                                 arrTracks2RP=&arrTracks1;
818                                 bTracks1RP = bTracks1;
819                                 bTracks2RP = bTracks1;
820                                 break;
821                 case 2: arr1RP=arr2;arr2RP=arr2;
822                                 arrTracks1RP=&arrTracks2;
823                                 arrTracks2RP=&arrTracks2;
824                                 bTracks1RP = bTracks2;
825                                 bTracks2RP = bTracks2;
826                                 break;
827                 default: ;//nothing to do
828         }
829         Int_t ntrack1RP=(*arrTracks1RP).GetEntriesFast();
830         Int_t ntrack2RP=(*arrTracks2RP).GetEntriesFast();
831
832         Int_t pairIndex=GetPairIndex(arr1RP,arr2RP);
833
834         for (Int_t itrack1=0; itrack1<ntrack1RP; ++itrack1){
835           Int_t end=ntrack2RP;
836           if (arr1RP==arr2RP) end=itrack1;
837           for (Int_t itrack2=0; itrack2<end; ++itrack2){
838                 TObject *track1=(*arrTracks1RP).UncheckedAt(itrack1);
839                 TObject *track2=(*arrTracks2RP).UncheckedAt(itrack2);
840                 if (!track1 || !track2) continue;
841                 //create the pair
842                 candidate.SetTracks(static_cast<AliVTrack*>(track1), fPdgLeg1,
843                         static_cast<AliVTrack*>(track2), fPdgLeg2);
844
845                 candidate.SetType(pairIndex);
846                 candidate.SetLabel(AliDielectronMC::Instance()->GetLabelMotherWithPdg(&candidate,fPdgMother));
847                 //relate to the production vertex
848                 //       if (AliDielectronVarManager::GetKFVertex()) candidate.SetProductionVertex(*AliDielectronVarManager::GetKFVertex());
849
850                 //pair cuts
851                 UInt_t cutMask=fPairPreFilter.IsSelected(&candidate);
852
853                 //apply cut
854                 if (cutMask!=selectedMask) continue;
855                 if (fCfManagerPair) fCfManagerPair->Fill(selectedMaskPair+1 ,&candidate);
856                 if (fHistos) FillHistogramsPair(&candidate,kTRUE);
857                 //set flags for track removal
858                 bTracks1RP[itrack1]=kTRUE;
859                 bTracks2RP[itrack2]=kTRUE;
860           }
861         }
862   }
863
864   //remove the tracks from the Track arrays
865   for (Int_t itrack1=0; itrack1<ntrack1; ++itrack1){
866     if(bTracks1[itrack1]) arrTracks1.AddAt(0x0, itrack1);
867   }
868   for (Int_t itrack2=0; itrack2<ntrack2; ++itrack2){
869     if(bTracks2[itrack2]) arrTracks2.AddAt(0x0, itrack2);
870   }
871
872   // clean up
873   delete [] bTracks1;
874   delete [] bTracks2;
875
876   //compress the track arrays
877   arrTracks1.Compress();
878   arrTracks2.Compress();
879   
880   //apply leg cuts after the pre filter
881   if ( fPairPreFilterLegs.GetCuts()->GetEntries()>0 ) {
882     selectedMask=(1<<fPairPreFilterLegs.GetCuts()->GetEntries())-1;
883     //loop over tracks from array 1
884     for (Int_t itrack=0; itrack<arrTracks1.GetEntriesFast();++itrack){
885       //test cuts
886       UInt_t cutMask=fPairPreFilterLegs.IsSelected(arrTracks1.UncheckedAt(itrack));
887       
888       //apply cut
889       if (cutMask!=selectedMask) arrTracks1.AddAt(0x0,itrack);;
890     }
891     arrTracks1.Compress();
892     
893     //in case of like sign don't loop over second array
894     if (arr1==arr2) {
895       arrTracks2=arrTracks1;
896     } else {
897       
898       //loop over tracks from array 2
899       for (Int_t itrack=0; itrack<arrTracks2.GetEntriesFast();++itrack){
900       //test cuts
901         UInt_t cutMask=fPairPreFilterLegs.IsSelected(arrTracks2.UncheckedAt(itrack));
902       //apply cut
903         if (cutMask!=selectedMask) arrTracks2.AddAt(0x0,itrack);
904       }
905       arrTracks2.Compress();
906       
907     }
908   }
909   //For unlike-sign monitor track-cuts:
910   if (arr1!=arr2&&fHistos) {
911     TObjArray *unlikesignArray[2] = {&arrTracks1,&arrTracks2};
912     FillHistogramsTracks(unlikesignArray);
913   }
914 }
915
916 //________________________________________________________________
917 void AliDielectron::FillPairArrays(Int_t arr1, Int_t arr2)
918 {
919   //
920   // select pairs and fill pair candidate arrays
921   //
922
923   TObjArray arrTracks1=fTracks[arr1];
924   TObjArray arrTracks2=fTracks[arr2];
925
926   //process pre filter if set
927   if ((!fPreFilterAllSigns) && (!fPreFilterUnlikeOnly) && ( fPairPreFilter.GetCuts()->GetEntries()>0 ))  PairPreFilter(arr1, arr2, arrTracks1, arrTracks2);
928   
929   Int_t pairIndex=GetPairIndex(arr1,arr2);
930
931   Int_t ntrack1=arrTracks1.GetEntriesFast();
932   Int_t ntrack2=arrTracks2.GetEntriesFast();
933
934   AliDielectronPair *candidate=new AliDielectronPair;
935
936   UInt_t selectedMask=(1<<fPairFilter.GetCuts()->GetEntries())-1;
937   
938   for (Int_t itrack1=0; itrack1<ntrack1; ++itrack1){
939     Int_t end=ntrack2;
940     if (arr1==arr2) end=itrack1;
941     for (Int_t itrack2=0; itrack2<end; ++itrack2){
942       //create the pair
943       candidate->SetTracks(static_cast<AliVTrack*>(arrTracks1.UncheckedAt(itrack1)), fPdgLeg1,
944                              static_cast<AliVTrack*>(arrTracks2.UncheckedAt(itrack2)), fPdgLeg2);
945       candidate->SetType(pairIndex);
946       Int_t label=AliDielectronMC::Instance()->GetLabelMotherWithPdg(candidate,fPdgMother);
947       candidate->SetLabel(label);
948       if (label>-1) candidate->SetPdgCode(fPdgMother);
949
950       //pair cuts
951       UInt_t cutMask=fPairFilter.IsSelected(candidate);
952       
953       //CF manager for the pair
954       if (fCfManagerPair) fCfManagerPair->Fill(cutMask,candidate);
955
956       //apply cut
957       if (cutMask!=selectedMask) continue;
958
959       //add the candidate to the candidate array 
960       PairArray(pairIndex)->Add(candidate);
961       //get a new candidate
962       candidate=new AliDielectronPair;
963     }
964   }
965   //delete the surplus candidate
966   delete candidate;
967 }
968
969 //________________________________________________________________
970 void AliDielectron::FillPairArrayTR()
971 {
972   //
973   // select pairs and fill pair candidate arrays
974   //
975   UInt_t selectedMask=(1<<fPairFilter.GetCuts()->GetEntries())-1;
976   
977   while ( fTrackRotator->NextCombination() ){
978     AliDielectronPair candidate;
979     candidate.SetTracks(&fTrackRotator->GetKFTrackP(), &fTrackRotator->GetKFTrackN(),
980                         fTrackRotator->GetVTrackP(),fTrackRotator->GetVTrackN());
981     candidate.SetType(kEv1PMRot);
982     
983     //pair cuts
984     UInt_t cutMask=fPairFilter.IsSelected(&candidate);
985     
986     //CF manager for the pair
987     if (fCfManagerPair) fCfManagerPair->Fill(cutMask,&candidate);
988     
989     //apply cut
990     if (cutMask==selectedMask) {
991      if(fHistos) FillHistogramsPair(&candidate);
992      if(fStoreRotatedPairs) PairArray(kEv1PMRot)->Add(new AliDielectronPair(candidate));
993     } 
994   }
995 }
996
997 //________________________________________________________________
998 void AliDielectron::FillDebugTree()
999 {
1000   //
1001   // Fill Histogram information for tracks and pairs
1002   //
1003   
1004   //Fill Debug tree
1005   for (Int_t i=0; i<10; ++i){
1006     Int_t ntracks=PairArray(i)->GetEntriesFast();
1007     for (Int_t ipair=0; ipair<ntracks; ++ipair){
1008       fDebugTree->Fill(static_cast<AliDielectronPair*>(PairArray(i)->UncheckedAt(ipair)));
1009     }
1010   }
1011 }
1012
1013 //________________________________________________________________
1014 void AliDielectron::SaveDebugTree()
1015 {
1016   //
1017   // delete the debug tree, this will also write the tree
1018   //
1019   if (fDebugTree) fDebugTree->DeleteStreamer();
1020 }
1021
1022
1023 //__________________________________________________________________
1024 void AliDielectron::AddSignalMC(AliDielectronSignalMC* signal) {
1025   //
1026   //  Add an MC signal to the signals list
1027   //
1028   if(!fSignalsMC) {
1029     fSignalsMC = new TObjArray();
1030     fSignalsMC->SetOwner();
1031   }
1032   fSignalsMC->Add(signal);
1033 }