o update dielectron package
[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   fEstimatorFilename(""),
121   fTRDpidCorrectionFilename("")
122 {
123   //
124   // Default constructor
125   //
126
127 }
128
129 //________________________________________________________________
130 AliDielectron::AliDielectron(const char* name, const char* title) :
131   TNamed(name,title),
132   fEventFilter("EventFilter"),
133   fTrackFilter("TrackFilter"),
134   fPairPreFilter("PairPreFilter"),
135   fPairPreFilterLegs("PairPreFilterLegs"),
136   fPairFilter("PairFilter"),
137   fEventPlanePreFilter("EventPlanePreFilter"),
138   fEventPlanePOIPreFilter("EventPlanePOIPreFilter"),
139   fPdgMother(443),
140   fPdgLeg1(11),
141   fPdgLeg2(11),
142   fSignalsMC(0x0),
143   fNoPairing(kFALSE),
144   fHistos(0x0),
145   fPairCandidates(new TObjArray(11)),
146   fCfManagerPair(0x0),
147   fTrackRotator(0x0),
148   fDebugTree(0x0),
149   fMixing(0x0),
150   fPreFilterEventPlane(kFALSE),
151   fLikeSignSubEvents(kFALSE),
152   fPreFilterUnlikeOnly(kFALSE),
153   fPreFilterAllSigns(kFALSE),
154   fHasMC(kFALSE),
155   fStoreRotatedPairs(kFALSE),
156   fEstimatorFilename(""),
157   fTRDpidCorrectionFilename("")
158 {
159   //
160   // Named constructor
161   //
162   
163 }
164
165 //________________________________________________________________
166 AliDielectron::~AliDielectron()
167 {
168   //
169   // Default destructor
170   //
171   if (fHistos) delete fHistos;
172   if (fPairCandidates) delete fPairCandidates;
173   if (fDebugTree) delete fDebugTree;
174   if (fMixing) delete fMixing;
175   if (fSignalsMC) delete fSignalsMC;
176   if (fCfManagerPair) delete fCfManagerPair;
177 }
178
179 //________________________________________________________________
180 void AliDielectron::Init()
181 {
182   //
183   // Initialise objects
184   //
185
186   if(GetHasMC()) AliDielectronMC::Instance()->SetHasMC(GetHasMC());
187    
188   InitPairCandidateArrays();
189    
190   if (fCfManagerPair) {
191     fCfManagerPair->SetSignalsMC(fSignalsMC);
192     fCfManagerPair->InitialiseContainer(fPairFilter);
193   }
194   if (fTrackRotator)  {
195     fTrackRotator->SetTrackArrays(&fTracks[0],&fTracks[1]);
196     fTrackRotator->SetPdgLegs(fPdgLeg1,fPdgLeg2);
197   }
198   if (fDebugTree) fDebugTree->SetDielectron(this);
199   if(fEstimatorFilename.Contains(".root")) AliDielectronVarManager::InitEstimatorAvg(fEstimatorFilename.Data());
200   if(fTRDpidCorrectionFilename.Contains(".root")) AliDielectronVarManager::InitTRDpidEffHistograms(fTRDpidCorrectionFilename.Data());
201
202   if (fMixing) fMixing->Init(this);
203
204
205 //________________________________________________________________
206 void AliDielectron::Process(AliVEvent *ev1, AliVEvent *ev2)
207 {
208   //
209   // Process the events
210   //
211
212   //at least first event is needed!
213   if (!ev1){
214     AliError("At least first event must be set!");
215     return;
216   }
217   AliDielectronVarManager::SetEvent(ev1);
218
219   //in case we have MC load the MC event and process the MC particles
220   if (AliDielectronMC::Instance()->HasMC()) {
221     if (!AliDielectronMC::Instance()->ConnectMCEvent()){
222       AliError("Could not properly connect the MC event, skipping this event!");
223       return;
224     }
225     ProcessMC(ev1);
226   }
227   
228   //if candidate array doesn't exist, create it
229   if (!fPairCandidates->UncheckedAt(0)) {
230     InitPairCandidateArrays();
231   } else {
232     ClearArrays();
233   }
234
235   //mask used to require that all cuts are fulfilled
236   UInt_t selectedMask=(1<<fEventFilter.GetCuts()->GetEntries())-1;
237
238   //apply event cuts
239     if ((ev1&&fEventFilter.IsSelected(ev1)!=selectedMask) ||
240         (ev2&&fEventFilter.IsSelected(ev2)!=selectedMask)) return;
241   
242   AliDielectronVarManager::SetEvent(ev1);
243
244   //fill track arrays for the first event
245   if (ev1){
246     FillTrackArrays(ev1);
247     if (((fPreFilterAllSigns)||(fPreFilterUnlikeOnly)) && ( fPairPreFilter.GetCuts()->GetEntries()>0 )) PairPreFilter(0, 1, fTracks[0], fTracks[1]);
248   }
249
250
251   //fill track arrays for the second event
252   if (ev2) {
253     FillTrackArrays(ev2,1);
254     if (((fPreFilterAllSigns)||(fPreFilterUnlikeOnly)) && ( fPairPreFilter.GetCuts()->GetEntries()>0 )) PairPreFilter(2, 3, fTracks[2], fTracks[3]);
255   }
256
257   // TPC event plane correction
258   AliEventplane *cevplane = new AliEventplane();
259   if (ev1 && fPreFilterEventPlane && ( fEventPlanePreFilter.GetCuts()->GetEntries()>0 || fEventPlanePOIPreFilter.GetCuts()->GetEntries()>0)) 
260     EventPlanePreFilter(0, 1, fTracks[0], fTracks[1], ev1, cevplane);
261   
262   if (!fNoPairing){
263     // create pairs and fill pair candidate arrays
264     for (Int_t itrackArr1=0; itrackArr1<4; ++itrackArr1){
265       for (Int_t itrackArr2=itrackArr1; itrackArr2<4; ++itrackArr2){
266         FillPairArrays(itrackArr1, itrackArr2);
267       }
268     }
269
270     //track rotation
271     if (fTrackRotator) {
272       fTrackRotator->SetEvent(ev1);
273       FillPairArrayTR();
274     }
275   }
276
277   //fill debug tree if a manager is attached
278   if (fDebugTree) FillDebugTree();
279
280   //process event mixing
281   if (fMixing) {
282     fMixing->Fill(ev1,this);
283 //     FillHistograms(0x0,kTRUE);
284   }
285
286   //in case there is a histogram manager, fill the QA histograms
287   if (fHistos) FillHistograms(ev1);
288
289   // clear arrays
290 //   ClearArrays();
291   AliDielectronVarManager::SetTPCEventPlane(0x0);
292   delete cevplane;
293 }
294
295 //________________________________________________________________
296 void AliDielectron::ProcessMC(AliVEvent *ev1)
297 {
298   //
299   // Process the MC data
300   //
301
302   AliDielectronMC *dieMC=AliDielectronMC::Instance();
303
304   if (fHistos) FillHistogramsMC(dieMC->GetMCEvent(), ev1);
305
306   if(!fSignalsMC) return;
307   //loop over all MC data and Fill the CF container if it exist
308   if (!fCfManagerPair) return;
309   fCfManagerPair->SetPdgMother(fPdgMother);
310   if(!fCfManagerPair->GetStepForMCtruth()) return;
311
312   // signals to be studied
313   Int_t nSignals = fSignalsMC->GetEntries();
314
315   // initialize 2D arrays of labels for particles from each MC signal
316   Int_t** labels1;      // labels for particles satisfying branch 1
317   Int_t** labels2;      // labels for particles satisfying branch 2
318   Int_t** labels12;     // labels for particles satisfying both branches
319   labels1 = new Int_t*[nSignals];
320   labels2 = new Int_t*[nSignals];
321   labels12 = new Int_t*[nSignals];
322   Int_t* indexes1=new Int_t[nSignals];
323   Int_t* indexes2=new Int_t[nSignals];
324   Int_t* indexes12=new Int_t[nSignals];
325   for(Int_t isig=0;isig<nSignals;++isig) {
326     *(labels1+isig) = new Int_t[dieMC->GetNMCTracks()];
327     *(labels2+isig) = new Int_t[dieMC->GetNMCTracks()];
328     *(labels12+isig) = new Int_t[dieMC->GetNMCTracks()];
329     for(Int_t ip=0; ip<dieMC->GetNMCTracks();++ip) {
330       labels1[isig][ip] = -1;
331       labels2[isig][ip] = -1;
332       labels12[isig][ip] = -1;
333     }
334     indexes1[isig]=0;
335     indexes2[isig]=0;
336     indexes12[isig]=0;
337   }
338
339   Bool_t truth1=kFALSE;
340   Bool_t truth2=kFALSE;
341   // loop over the MC tracks
342   for(Int_t ipart=0; ipart<dieMC->GetNMCTracks(); ++ipart) {
343     for(Int_t isig=0; isig<nSignals; ++isig) {       // loop over signals
344       // Proceed only if this signal is required in the pure MC step
345       // NOTE: Some signals can be satisfied by many particles and this leads to high
346       //       computation times (e.g. secondary electrons from the GEANT transport). Be aware of this!!
347       if(!((AliDielectronSignalMC*)fSignalsMC->At(isig))->GetFillPureMCStep()) continue;
348
349       truth1 = dieMC->IsMCTruth(ipart, (AliDielectronSignalMC*)fSignalsMC->At(isig), 1);
350       truth2 = dieMC->IsMCTruth(ipart, (AliDielectronSignalMC*)fSignalsMC->At(isig), 2);
351
352       // particles satisfying both branches are treated separately to avoid double counting during pairing
353       if(truth1 && truth2) {
354         labels12[isig][indexes12[isig]] = ipart;
355         ++indexes12[isig];
356       }
357       else {
358         if(truth1) {
359           labels1[isig][indexes1[isig]] = ipart;
360           ++indexes1[isig];
361         }
362         if(truth2) {
363           labels2[isig][indexes2[isig]] = ipart;
364           ++indexes2[isig];
365         }
366       }
367     }
368   }  // end loop over MC particles
369
370   // Do the pairing and fill the CF container with pure MC info
371   for(Int_t isig=0; isig<nSignals; ++isig) {
372     // mix the particles which satisfy only one of the signal branches
373     for(Int_t i1=0;i1<indexes1[isig];++i1) {
374       for(Int_t i2=0;i2<indexes2[isig];++i2) {
375         fCfManagerPair->FillMC(labels1[isig][i1], labels2[isig][i2], isig);
376       }
377     }
378     // mix the particles which satisfy both branches
379     for(Int_t i1=0;i1<indexes12[isig];++i1) {
380       for(Int_t i2=0; i2<i1; ++i2) {
381         fCfManagerPair->FillMC(labels12[isig][i1], labels12[isig][i2], isig);
382       }
383     }
384   }    // end loop over signals
385
386   // release the memory
387   for(Int_t isig=0;isig<nSignals;++isig) {
388     delete [] *(labels1+isig);
389     delete [] *(labels2+isig);
390     delete [] *(labels12+isig);
391   }
392   delete [] labels1;
393   delete [] labels2;
394   delete [] labels12;
395   delete [] indexes1;
396   delete [] indexes2;
397   delete [] indexes12;
398 }
399
400 //________________________________________________________________
401 void AliDielectron::FillHistogramsTracks(TObjArray **tracks)
402 {
403   //
404   // Fill Histogram information for tracks after prefilter
405   // ignore mixed events - for prefilter, only single tracks +/- are relevant 
406   //
407   
408   TString  className,className2;
409   Double_t values[AliDielectronVarManager::kNMaxValues];
410   
411   //Fill track information, separately for the track array candidates
412   for (Int_t i=0; i<2; ++i){
413     className.Form("Pre_%s",fgkTrackClassNames[i]);
414     if (!fHistos->GetHistogramList()->FindObject(className.Data())) continue;
415     Int_t ntracks=tracks[i]->GetEntriesFast();
416     for (Int_t itrack=0; itrack<ntracks; ++itrack){
417       AliDielectronVarManager::Fill(tracks[i]->UncheckedAt(itrack), values);
418       fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
419     }
420   }
421 }
422
423
424 //________________________________________________________________
425 void AliDielectron::FillHistogramsMC(const AliMCEvent *ev, AliVEvent *ev1)
426 {
427   //
428   // Fill Histogram information for MCEvents
429   //
430
431   Double_t values[AliDielectronVarManager::kNMaxValues]={0.};
432   // Fill event information
433   AliDielectronVarManager::Fill(ev1, values);    // ESD/AOD information
434   AliDielectronVarManager::Fill(ev, values);     // MC truth info
435   if (fHistos->GetHistogramList()->FindObject("MCEvent"))
436     fHistos->FillClass("MCEvent", AliDielectronVarManager::kNMaxValues, values);
437 }
438
439
440 //________________________________________________________________
441 void AliDielectron::FillHistograms(const AliVEvent *ev, Bool_t pairInfoOnly)
442 {
443   //
444   // Fill Histogram information for tracks and pairs
445   //
446   
447   TString  className,className2;
448   Double_t values[AliDielectronVarManager::kNMaxValues]={0.};
449   //Fill event information
450   if (ev){
451     AliDielectronVarManager::Fill(ev, values);
452     if (fHistos->GetHistogramList()->FindObject("Event"))
453       fHistos->FillClass("Event", AliDielectronVarManager::kNMaxValues, values);
454   }
455   
456   //Fill track information, separately for the track array candidates
457   if (!pairInfoOnly){
458     for (Int_t i=0; i<4; ++i){
459       className.Form("Track_%s",fgkTrackClassNames[i]);
460       if (!fHistos->GetHistogramList()->FindObject(className.Data())) continue;
461       Int_t ntracks=fTracks[i].GetEntriesFast();
462       for (Int_t itrack=0; itrack<ntracks; ++itrack){
463         AliDielectronVarManager::Fill(fTracks[i].UncheckedAt(itrack), values);
464         fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
465       }
466     }
467   }
468
469   //Fill Pair information, separately for all pair candidate arrays and the legs
470   TObjArray arrLegs(100);
471   for (Int_t i=0; i<10; ++i){
472     className.Form("Pair_%s",fgkPairClassNames[i]);
473     className2.Form("Track_Legs_%s",fgkPairClassNames[i]);
474     Bool_t pairClass=fHistos->GetHistogramList()->FindObject(className.Data())!=0x0;
475     Bool_t legClass=fHistos->GetHistogramList()->FindObject(className2.Data())!=0x0;
476     if (!pairClass&&!legClass) continue;
477     Int_t ntracks=PairArray(i)->GetEntriesFast();
478     for (Int_t ipair=0; ipair<ntracks; ++ipair){
479       AliDielectronPair *pair=static_cast<AliDielectronPair*>(PairArray(i)->UncheckedAt(ipair));
480       
481       //fill pair information
482       if (pairClass){
483         AliDielectronVarManager::Fill(pair, values);
484         fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
485       }
486
487       //fill leg information, don't fill the information twice
488       if (legClass){
489         AliVParticle *d1=pair->GetFirstDaughter();
490         AliVParticle *d2=pair->GetSecondDaughter();
491         if (!arrLegs.FindObject(d1)){
492           AliDielectronVarManager::Fill(d1, values);
493           fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
494           arrLegs.Add(d1);
495         }
496         if (!arrLegs.FindObject(d2)){
497           AliDielectronVarManager::Fill(d2, values);
498           fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
499           arrLegs.Add(d2);
500         }
501       }
502     }
503     if (legClass) arrLegs.Clear();
504   }
505   
506 }
507 //________________________________________________________________
508 void AliDielectron::FillHistogramsPair(AliDielectronPair *pair,Bool_t fromPreFilter/*=kFALSE*/)
509 {
510   //
511   // Fill Histogram information for pairs and the track in the pair
512   // NOTE: in this funtion the leg information may be filled multiple
513   //       times. This funtion is used in the track rotation pairing
514   //       and those legs are not saved!
515   //
516   TString  className,className2;
517   Double_t values[AliDielectronVarManager::kNMaxValues];
518   
519   //Fill Pair information, separately for all pair candidate arrays and the legs
520   TObjArray arrLegs(100);
521   const Int_t type=pair->GetType();
522   if (fromPreFilter) {
523     className.Form("RejPair_%s",fgkPairClassNames[type]);
524     className2.Form("RejTrack_%s",fgkPairClassNames[type]);
525   } else {
526     className.Form("Pair_%s",fgkPairClassNames[type]);
527     className2.Form("Track_Legs_%s",fgkPairClassNames[type]);
528   }
529   
530   Bool_t pairClass=fHistos->GetHistogramList()->FindObject(className.Data())!=0x0;
531   Bool_t legClass=fHistos->GetHistogramList()->FindObject(className2.Data())!=0x0;
532   
533   //fill pair information
534   if (pairClass){
535     AliDielectronVarManager::Fill(pair, values);
536     fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
537   }
538
539   if (legClass){
540     AliVParticle *d1=pair->GetFirstDaughter();
541     AliDielectronVarManager::Fill(d1, values);
542     fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
543     
544     AliVParticle *d2=pair->GetSecondDaughter();
545     AliDielectronVarManager::Fill(d2, values);
546     fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
547   }
548 }
549
550 //________________________________________________________________
551 void AliDielectron::FillTrackArrays(AliVEvent * const ev, Int_t eventNr)
552 {
553   //
554   // select tracks and fill track candidate arrays
555   // eventNr = 0: First  event, use track arrays 0 and 1
556   // eventNr = 1: Second event, use track arrays 2 and 3
557   //
558   
559   Int_t ntracks=ev->GetNumberOfTracks();
560   UInt_t selectedMask=(1<<fTrackFilter.GetCuts()->GetEntries())-1;
561   for (Int_t itrack=0; itrack<ntracks; ++itrack){
562     //get particle
563     AliVParticle *particle=ev->GetTrack(itrack);
564     //TODO: temporary solution, perhaps think about a better implementation
565     //      This is needed to use AliESDpidCuts, which relies on the ESD event
566     //      is set as a AliESDtrack attribute... somehow ugly!
567     if (ev->IsA()==AliESDEvent::Class()){
568       AliESDtrack *track=static_cast<AliESDtrack*>(particle);
569       track->SetESDEvent(static_cast<AliESDEvent*>(ev)); //only in trunk...
570     }
571     
572     //apply track cuts
573     if (fTrackFilter.IsSelected(particle)!=selectedMask) continue;
574
575     //fill selected particle into the corresponding track arrays
576     Short_t charge=particle->Charge();
577     if (charge>0)      fTracks[eventNr*2].Add(particle);
578     else if (charge<0) fTracks[eventNr*2+1].Add(particle);
579   }
580 }
581
582 //________________________________________________________________
583 void AliDielectron::EventPlanePreFilter(Int_t arr1, Int_t arr2, TObjArray arrTracks1, TObjArray arrTracks2, const AliVEvent *ev, AliEventplane *cevplane)
584 {
585   //
586   // Prefilter tracks and tracks from pairs
587   // Needed for rejection in the Q-Vector of the event plane
588   // remove contribution of all tracks to the Q-vector that are in invariant mass window 
589   //
590   AliEventplane *evplane = const_cast<AliVEvent *>(ev)->GetEventplane();
591   if(!evplane) return;
592   
593   // do not change these vectors qref
594   TVector2 *qref  = evplane->GetQVector();
595   if(!qref) return;
596   // random subevents
597   TVector2 *qrsub1 = evplane->GetQsub1();
598   TVector2 *qrsub2 = evplane->GetQsub2();
599
600   // copy references
601   TVector2 *qstd   = dynamic_cast<TVector2 *>(qref->Clone());
602   TVector2 *qssub1 = dynamic_cast<TVector2 *>(qrsub1->Clone());
603   TVector2 *qssub2 = dynamic_cast<TVector2 *>(qrsub2->Clone());
604
605   // subevents by charge separation sub1+ sub2-
606   if(fLikeSignSubEvents) {
607     qssub1 = dynamic_cast<TVector2 *>(qstd->Clone());
608     qssub2 = dynamic_cast<TVector2 *>(qstd->Clone());
609
610     Int_t ntracks=ev->GetNumberOfTracks();
611     
612     // track removals
613     for (Int_t itrack=0; itrack<ntracks; ++itrack){
614       AliVParticle *particle=ev->GetTrack(itrack);
615       AliVTrack *track= static_cast<AliVTrack*>(particle);
616       if (!track) continue;
617
618       Double_t cQX     = evplane->GetQContributionX(track);
619       Double_t cQY     = evplane->GetQContributionY(track);
620       
621       Short_t charge=track->Charge();
622       if (charge<0) qssub1->Set(qssub1->X()-cQX, qssub1->Y()-cQY);
623       if (charge>0) qssub2->Set(qssub2->X()-cQX, qssub2->Y()-cQY);
624     }
625   }
626
627   // subevents eta division sub1+ sub2-
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")) {
632       etagap=kTRUE;
633
634       qssub1 = dynamic_cast<TVector2 *>(qstd->Clone());
635       qssub2 = dynamic_cast<TVector2 *>(qstd->Clone());
636       
637       Int_t ntracks=ev->GetNumberOfTracks();
638       
639       // track removals
640       for (Int_t itrack=0; itrack<ntracks; ++itrack){
641         AliVParticle *particle=ev->GetTrack(itrack);
642         AliVTrack *track= static_cast<AliVTrack*>(particle);
643         if (!track) continue;
644         
645         Double_t cQX     = evplane->GetQContributionX(track);
646         Double_t cQY     = evplane->GetQContributionY(track);
647       
648         Double_t eta=track->Eta();
649         if (eta<0.0) qssub1->Set(qssub1->X()-cQX, qssub1->Y()-cQY);
650         if (eta>0.0) qssub2->Set(qssub2->X()-cQX, qssub2->Y()-cQY);
651       }
652     }
653   }
654   
655   // apply cuts, e.g. etagap 
656   if(fEventPlanePreFilter.GetCuts()->GetEntries()) {
657     UInt_t selectedMask=(1<<fEventPlanePreFilter.GetCuts()->GetEntries())-1;
658     Int_t ntracks=ev->GetNumberOfTracks();
659     for (Int_t itrack=0; itrack<ntracks; ++itrack){
660       AliVParticle *particle=ev->GetTrack(itrack);
661       AliVTrack *track= static_cast<AliVTrack*>(particle);
662       if (!track) continue;
663       
664       //event plane cuts
665       UInt_t cutMask=fEventPlanePreFilter.IsSelected(track);
666       //apply cut
667       if (cutMask==selectedMask) continue;
668
669       Double_t cQX     = 0.0;
670       Double_t cQY     = 0.0;
671       if(!etagap) {
672         cQX = evplane->GetQContributionX(track);
673         cQY = evplane->GetQContributionY(track);
674       }
675       Double_t cQXsub1 = evplane->GetQContributionXsub1(track);
676       Double_t cQYsub1 = evplane->GetQContributionYsub1(track);
677       Double_t cQXsub2 = evplane->GetQContributionXsub2(track);
678       Double_t cQYsub2 = evplane->GetQContributionYsub2(track);      
679
680       // update Q vectors
681       qstd->Set(qstd->X()-cQX, qstd->Y()-cQY);
682       qssub1->Set(qssub1->X()-cQXsub1, qssub1->Y()-cQYsub1);
683       qssub2->Set(qssub2->X()-cQXsub2, qssub2->Y()-cQYsub2);
684     }
685   }
686
687   if(!qstd || !qssub1 || !qssub2) return;
688
689   TVector2 *qcorr  = dynamic_cast<TVector2 *>(qstd->Clone());
690   TVector2 *qcsub1 = dynamic_cast<TVector2 *>(qssub1->Clone());
691   TVector2 *qcsub2 = dynamic_cast<TVector2 *>(qssub2->Clone());
692
693
694   // POI (particle of interest) rejection
695   Int_t pairIndex=GetPairIndex(arr1,arr2);
696   
697   Int_t ntrack1=arrTracks1.GetEntriesFast();
698   Int_t ntrack2=arrTracks2.GetEntriesFast();
699   AliDielectronPair candidate;
700   
701   UInt_t selectedMask=(1<<fEventPlanePOIPreFilter.GetCuts()->GetEntries())-1;
702   for (Int_t itrack1=0; itrack1<ntrack1; ++itrack1){
703     Int_t end=ntrack2;
704     if (arr1==arr2) end=itrack1;
705     Bool_t accepted=kFALSE;
706     for (Int_t itrack2=0; itrack2<end; ++itrack2){
707       TObject *track1=arrTracks1.UncheckedAt(itrack1);
708       TObject *track2=arrTracks2.UncheckedAt(itrack2);
709       if (!track1 || !track2) continue;
710       //create the pair
711       candidate.SetTracks(static_cast<AliVTrack*>(track1), fPdgLeg1,
712                           static_cast<AliVTrack*>(track2), fPdgLeg2);
713       
714       candidate.SetType(pairIndex);
715       candidate.SetLabel(AliDielectronMC::Instance()->GetLabelMotherWithPdg(&candidate,fPdgMother));
716       
717       //event plane cuts
718       UInt_t cutMask=fEventPlanePOIPreFilter.IsSelected(&candidate);
719       //apply cut
720       if (cutMask==selectedMask) continue;
721
722       accepted=kTRUE;
723       //remove the tracks from the Track arrays
724       arrTracks2.AddAt(0x0,itrack2);
725     }
726       if ( accepted ) arrTracks1.AddAt(0x0,itrack1);
727   }
728   //compress the track arrays
729   arrTracks1.Compress();
730   arrTracks2.Compress();
731   
732   
733   //Modify the components: subtract the tracks
734   ntrack1=arrTracks1.GetEntriesFast();
735   ntrack2=arrTracks2.GetEntriesFast();
736   
737   // remove leg1 contribution
738   for (Int_t itrack=0; itrack<ntrack1; ++itrack){
739     AliVTrack *track= static_cast<AliVTrack*>(arrTracks1.UncheckedAt(itrack));
740     if (!track) continue;
741     
742     Double_t cQX     = evplane->GetQContributionX(track);
743     Double_t cQY     = evplane->GetQContributionY(track);
744     Double_t cQXsub1 = evplane->GetQContributionXsub1(track);
745     Double_t cQYsub1 = evplane->GetQContributionYsub1(track);
746     Double_t cQXsub2 = evplane->GetQContributionXsub2(track);
747     Double_t cQYsub2 = evplane->GetQContributionYsub2(track);
748     
749     // update Q vectors
750     qcorr->Set(qcorr->X()-cQX, qcorr->Y()-cQY);
751     qcsub1->Set(qcsub1->X()-cQXsub1, qcsub1->Y()-cQYsub1);
752     qcsub2->Set(qcsub2->X()-cQXsub2, qcsub2->Y()-cQYsub2);
753   }
754   // remove leg2 contribution
755   for (Int_t itrack=0; itrack<ntrack2; ++itrack){
756     AliVTrack *track= static_cast<AliVTrack*>(arrTracks2.UncheckedAt(itrack));
757     if (!track) continue;
758     
759     Double_t cQX     = evplane->GetQContributionX(track);
760     Double_t cQY     = evplane->GetQContributionY(track);
761     Double_t cQXsub1 = evplane->GetQContributionXsub1(track);
762     Double_t cQYsub1 = evplane->GetQContributionYsub1(track);
763     Double_t cQXsub2 = evplane->GetQContributionXsub2(track);
764     Double_t cQYsub2 = evplane->GetQContributionYsub2(track);
765     
766     // update Q vectors
767     qcorr->Set(qcorr->X()-cQX, qcorr->Y()-cQY);
768     qcsub1->Set(qcsub1->X()-cQXsub1, qcsub1->Y()-cQYsub1);
769     qcsub2->Set(qcsub2->X()-cQXsub2, qcsub2->Y()-cQYsub2);
770   }
771
772   // set AliEventplane with corrected values
773   cevplane->SetQVector(qcorr);
774   cevplane->SetQsub(qcsub1, qcsub2);
775   AliDielectronVarManager::SetTPCEventPlane(cevplane);
776 }
777
778 //________________________________________________________________
779 void 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;
790
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
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));
882       
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   }
904   //For unlike-sign monitor track-cuts:
905   if (arr1!=arr2&&fHistos) {
906     TObjArray *unlikesignArray[2] = {&arrTracks1,&arrTracks2};
907     FillHistogramsTracks(unlikesignArray);
908   }
909 }
910
911 //________________________________________________________________
912 void AliDielectron::FillPairArrays(Int_t arr1, Int_t arr2)
913 {
914   //
915   // select pairs and fill pair candidate arrays
916   //
917
918   TObjArray arrTracks1=fTracks[arr1];
919   TObjArray arrTracks2=fTracks[arr2];
920
921   //process pre filter if set
922   if ((!fPreFilterAllSigns) && (!fPreFilterUnlikeOnly) && ( fPairPreFilter.GetCuts()->GetEntries()>0 ))  PairPreFilter(arr1, arr2, arrTracks1, arrTracks2);
923   
924   Int_t pairIndex=GetPairIndex(arr1,arr2);
925
926   Int_t ntrack1=arrTracks1.GetEntriesFast();
927   Int_t ntrack2=arrTracks2.GetEntriesFast();
928
929   AliDielectronPair *candidate=new AliDielectronPair;
930
931   UInt_t selectedMask=(1<<fPairFilter.GetCuts()->GetEntries())-1;
932   
933   for (Int_t itrack1=0; itrack1<ntrack1; ++itrack1){
934     Int_t end=ntrack2;
935     if (arr1==arr2) end=itrack1;
936     for (Int_t itrack2=0; itrack2<end; ++itrack2){
937       //create the pair
938       candidate->SetTracks(static_cast<AliVTrack*>(arrTracks1.UncheckedAt(itrack1)), fPdgLeg1,
939                              static_cast<AliVTrack*>(arrTracks2.UncheckedAt(itrack2)), fPdgLeg2);
940       candidate->SetType(pairIndex);
941       candidate->SetLabel(AliDielectronMC::Instance()->GetLabelMotherWithPdg(candidate,fPdgMother));
942
943       //pair cuts
944       UInt_t cutMask=fPairFilter.IsSelected(candidate);
945       
946       //CF manager for the pair
947       if (fCfManagerPair) fCfManagerPair->Fill(cutMask,candidate);
948
949       //apply cut
950       if (cutMask!=selectedMask) continue;
951
952       //add the candidate to the candidate array 
953       PairArray(pairIndex)->Add(candidate);
954       //get a new candidate
955       candidate=new AliDielectronPair;
956     }
957   }
958   //delete the surplus candidate
959   delete candidate;
960 }
961
962 //________________________________________________________________
963 void AliDielectron::FillPairArrayTR()
964 {
965   //
966   // select pairs and fill pair candidate arrays
967   //
968   UInt_t selectedMask=(1<<fPairFilter.GetCuts()->GetEntries())-1;
969   
970   while ( fTrackRotator->NextCombination() ){
971     AliDielectronPair candidate;
972     candidate.SetTracks(&fTrackRotator->GetKFTrackP(), &fTrackRotator->GetKFTrackN(),
973                         fTrackRotator->GetVTrackP(),fTrackRotator->GetVTrackN());
974     candidate.SetType(kEv1PMRot);
975     
976     //pair cuts
977     UInt_t cutMask=fPairFilter.IsSelected(&candidate);
978     
979     //CF manager for the pair
980     if (fCfManagerPair) fCfManagerPair->Fill(cutMask,&candidate);
981     
982     //apply cut
983     if (cutMask==selectedMask) {
984      if(fHistos) FillHistogramsPair(&candidate);
985      if(fStoreRotatedPairs) PairArray(10)->Add(new AliDielectronPair(candidate));  
986     } 
987   }
988 }
989
990 //________________________________________________________________
991 void AliDielectron::FillDebugTree()
992 {
993   //
994   // Fill Histogram information for tracks and pairs
995   //
996   
997   //Fill Debug tree
998   for (Int_t i=0; i<10; ++i){
999     Int_t ntracks=PairArray(i)->GetEntriesFast();
1000     for (Int_t ipair=0; ipair<ntracks; ++ipair){
1001       fDebugTree->Fill(static_cast<AliDielectronPair*>(PairArray(i)->UncheckedAt(ipair)));
1002     }
1003   }
1004 }
1005
1006 //________________________________________________________________
1007 void AliDielectron::SaveDebugTree()
1008 {
1009   //
1010   // delete the debug tree, this will also write the tree
1011   //
1012   if (fDebugTree) fDebugTree->DeleteStreamer();
1013 }
1014
1015
1016 //__________________________________________________________________
1017 void AliDielectron::AddSignalMC(AliDielectronSignalMC* signal) {
1018   //
1019   //  Add an MC signal to the signals list
1020   //
1021   if(!fSignalsMC) {
1022     fSignalsMC = new TObjArray();
1023     fSignalsMC->SetOwner();
1024   }
1025   fSignalsMC->Add(signal);
1026 }