-add cut qa class
[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 <AliKFParticle.h>
51
52 #include <AliESDInputHandler.h>
53 #include <AliAnalysisManager.h>
54 #include <AliEPSelectionTask.h>
55 #include <AliEventplane.h>
56 #include <AliVEvent.h>
57 #include <AliVParticle.h>
58 #include <AliVTrack.h>
59 #include "AliDielectronPair.h"
60 #include "AliDielectronHistos.h"
61 #include "AliDielectronCF.h"
62 #include "AliDielectronMC.h"
63 #include "AliDielectronVarManager.h"
64 #include "AliDielectronTrackRotator.h"
65 #include "AliDielectronDebugTree.h"
66 #include "AliDielectronSignalMC.h"
67 #include "AliDielectronMixingHandler.h"
68 #include "AliDielectronV0Cuts.h"
69
70 #include "AliDielectron.h"
71
72 ClassImp(AliDielectron)
73
74 const char* AliDielectron::fgkTrackClassNames[4] = {
75   "ev1+",
76   "ev1-",
77   "ev2+",
78   "ev2-"
79 };
80
81 const char* AliDielectron::fgkPairClassNames[11] = {
82   "ev1+_ev1+",
83   "ev1+_ev1-",
84   "ev1-_ev1-",
85   "ev1+_ev2+",
86   "ev1-_ev2+",
87   "ev2+_ev2+",
88   "ev1+_ev2-",
89   "ev1-_ev2-",
90   "ev2+_ev2-",
91   "ev2-_ev2-",
92   "ev1+_ev1-_TR"
93 };
94
95 //________________________________________________________________
96 AliDielectron::AliDielectron() :
97   TNamed("AliDielectron","AliDielectron"),
98   fCutQA(kFALSE),
99   fQAmonitor(0x0),
100   fEventFilter("EventFilter"),
101   fTrackFilter("TrackFilter"),
102   fPairPreFilter("PairPreFilter"),
103   fPairPreFilterLegs("PairPreFilterLegs"),
104   fPairFilter("PairFilter"),
105   fEventPlanePreFilter("EventPlanePreFilter"),
106   fEventPlanePOIPreFilter("EventPlanePOIPreFilter"),
107   fPdgMother(443),
108   fPdgLeg1(11),
109   fPdgLeg2(11),
110   fSignalsMC(0x0),
111   fNoPairing(kFALSE),
112   fUseKF(kTRUE),
113   fHistoArray(0x0),
114   fHistos(0x0),
115   fPairCandidates(new TObjArray(11)),
116   fCfManagerPair(0x0),
117   fTrackRotator(0x0),
118   fDebugTree(0x0),
119   fMixing(0x0),
120   fPreFilterEventPlane(kFALSE),
121   fLikeSignSubEvents(kFALSE),
122   fPreFilterUnlikeOnly(kFALSE),
123   fPreFilterAllSigns(kFALSE),
124   fHasMC(kFALSE),
125   fStoreRotatedPairs(kFALSE),
126   fDontClearArrays(kFALSE),
127   fEstimatorFilename(""),
128   fTRDpidCorrectionFilename(""),
129   fVZEROCalibrationFilename(""),
130   fVZERORecenteringFilename("")
131 {
132   //
133   // Default constructor
134   //
135
136 }
137
138 //________________________________________________________________
139 AliDielectron::AliDielectron(const char* name, const char* title) :
140   TNamed(name,title),
141   fCutQA(kFALSE),
142   fQAmonitor(0x0),
143   fEventFilter("EventFilter"),
144   fTrackFilter("TrackFilter"),
145   fPairPreFilter("PairPreFilter"),
146   fPairPreFilterLegs("PairPreFilterLegs"),
147   fPairFilter("PairFilter"),
148   fEventPlanePreFilter("EventPlanePreFilter"),
149   fEventPlanePOIPreFilter("EventPlanePOIPreFilter"),
150   fPdgMother(443),
151   fPdgLeg1(11),
152   fPdgLeg2(11),
153   fSignalsMC(0x0),
154   fNoPairing(kFALSE),
155   fUseKF(kTRUE),
156   fHistoArray(0x0),
157   fHistos(0x0),
158   fPairCandidates(new TObjArray(11)),
159   fCfManagerPair(0x0),
160   fTrackRotator(0x0),
161   fDebugTree(0x0),
162   fMixing(0x0),
163   fPreFilterEventPlane(kFALSE),
164   fLikeSignSubEvents(kFALSE),
165   fPreFilterUnlikeOnly(kFALSE),
166   fPreFilterAllSigns(kFALSE),
167   fHasMC(kFALSE),
168   fStoreRotatedPairs(kFALSE),
169   fDontClearArrays(kFALSE),
170   fEstimatorFilename(""),
171   fTRDpidCorrectionFilename(""),
172   fVZEROCalibrationFilename(""),
173   fVZERORecenteringFilename("")
174 {
175   //
176   // Named constructor
177   //
178   
179 }
180
181 //________________________________________________________________
182 AliDielectron::~AliDielectron()
183 {
184   //
185   // Default destructor
186   //
187   if (fQAmonitor) delete fQAmonitor;
188   if (fHistoArray) delete fHistoArray;
189   if (fHistos) delete fHistos;
190   if (fPairCandidates) delete fPairCandidates;
191   if (fDebugTree) delete fDebugTree;
192   if (fMixing) delete fMixing;
193   if (fSignalsMC) delete fSignalsMC;
194   if (fCfManagerPair) delete fCfManagerPair;
195 }
196
197 //________________________________________________________________
198 void AliDielectron::Init()
199 {
200   //
201   // Initialise objects
202   //
203
204   if(GetHasMC()) AliDielectronMC::Instance()->SetHasMC(GetHasMC());
205    
206   InitPairCandidateArrays();
207    
208   if (fCfManagerPair) {
209     fCfManagerPair->SetSignalsMC(fSignalsMC);
210     fCfManagerPair->InitialiseContainer(fPairFilter);
211   }
212   if (fTrackRotator)  {
213     fTrackRotator->SetTrackArrays(&fTracks[0],&fTracks[1]);
214     fTrackRotator->SetPdgLegs(fPdgLeg1,fPdgLeg2);
215   }
216   if (fDebugTree) fDebugTree->SetDielectron(this);
217   if(fEstimatorFilename.Contains(".root")) AliDielectronVarManager::InitEstimatorAvg(fEstimatorFilename.Data());
218   if(fTRDpidCorrectionFilename.Contains(".root")) AliDielectronVarManager::InitTRDpidEffHistograms(fTRDpidCorrectionFilename.Data());
219   if(fVZEROCalibrationFilename.Contains(".root")) AliDielectronVarManager::SetVZEROCalibrationFile(fVZEROCalibrationFilename.Data());
220   if(fVZERORecenteringFilename.Contains(".root")) AliDielectronVarManager::SetVZERORecenteringFile(fVZERORecenteringFilename.Data());
221   
222   if (fMixing) fMixing->Init(this);
223   if (fHistoArray) {
224     fHistoArray->SetSignalsMC(fSignalsMC);
225     fHistoArray->Init();
226   }
227
228   if (fCutQA) {
229     fQAmonitor = new AliDielectronCutQA(Form("QAcuts_%s",GetName()),"QAcuts");
230     fQAmonitor->AddTrackFilter(&fTrackFilter);
231     fQAmonitor->AddPairFilter(&fPairFilter);
232     fQAmonitor->AddEventFilter(&fEventFilter);
233     fQAmonitor->Init();
234   }
235 }
236
237 //________________________________________________________________
238 void AliDielectron::Process(AliVEvent *ev1, AliVEvent *ev2)
239 {
240   //
241   // Process the events
242   //
243
244   //at least first event is needed!
245   if (!ev1){
246     AliError("At least first event must be set!");
247     return;
248   }
249
250   // modify event numbers in MC so that we can identify new events 
251   // in AliDielectronV0Cuts (not neeeded for collision data)
252   if(GetHasMC()) {
253     ev1->SetBunchCrossNumber(1);
254     ev1->SetOrbitNumber(1);
255     ev1->SetPeriodNumber(1);
256   }
257
258   // set event
259   AliDielectronVarManager::SetEvent(ev1);
260   if (fMixing){
261     //set mixing bin to event data
262     Int_t bin=fMixing->FindBin(AliDielectronVarManager::GetData());
263     AliDielectronVarManager::SetValue(AliDielectronVarManager::kMixingBin,bin);
264   }
265
266   //in case we have MC load the MC event and process the MC particles
267   // why do not apply the event cuts first ????
268   if (AliDielectronMC::Instance()->ConnectMCEvent()){
269     ProcessMC(ev1);
270   }
271
272   //if candidate array doesn't exist, create it
273   if (!fPairCandidates->UncheckedAt(0)) {
274     InitPairCandidateArrays();
275   } else {
276     ClearArrays();
277   }
278
279   //mask used to require that all cuts are fulfilled
280   UInt_t selectedMask=(1<<fEventFilter.GetCuts()->GetEntries())-1;
281
282   //apply event cuts
283   UInt_t cutmask = fEventFilter.IsSelected(ev1);
284   if(fCutQA) fQAmonitor->FillAll(ev1);
285   if(fCutQA) fQAmonitor->Fill(cutmask,ev1);
286   if ((ev1&&cutmask!=selectedMask) ||
287       (ev2&&fEventFilter.IsSelected(ev2)!=selectedMask)) return;
288
289   //fill track arrays for the first event
290   if (ev1){
291     FillTrackArrays(ev1);
292     if (((fPreFilterAllSigns)||(fPreFilterUnlikeOnly)) && ( fPairPreFilter.GetCuts()->GetEntries()>0 )) PairPreFilter(0, 1, fTracks[0], fTracks[1]);
293   }
294
295
296   //fill track arrays for the second event
297   if (ev2) {
298     FillTrackArrays(ev2,1);
299     if (((fPreFilterAllSigns)||(fPreFilterUnlikeOnly)) && ( fPairPreFilter.GetCuts()->GetEntries()>0 )) PairPreFilter(2, 3, fTracks[2], fTracks[3]);
300   }
301
302   // TPC event plane correction
303   if (ev1 && fPreFilterEventPlane && ( fEventPlanePreFilter.GetCuts()->GetEntries()>0 || fEventPlanePOIPreFilter.GetCuts()->GetEntries()>0)) 
304     EventPlanePreFilter(0, 1, fTracks[0], fTracks[1], ev1);
305
306   if (!fNoPairing){
307     // create pairs and fill pair candidate arrays
308     for (Int_t itrackArr1=0; itrackArr1<4; ++itrackArr1){
309       for (Int_t itrackArr2=itrackArr1; itrackArr2<4; ++itrackArr2){
310         FillPairArrays(itrackArr1, itrackArr2);
311       }
312     }
313
314     //track rotation
315     if (fTrackRotator) {
316       fTrackRotator->SetEvent(ev1);
317       FillPairArrayTR();
318     }
319   }
320
321   //fill debug tree if a manager is attached
322   if (fDebugTree) FillDebugTree();
323
324   //process event mixing
325   if (fMixing) {
326     fMixing->Fill(ev1,this);
327     //     FillHistograms(0x0,kTRUE);
328   }
329
330   //in case there is a histogram manager, fill the QA histograms
331   if (fHistos && fSignalsMC) FillMCHistograms(ev1);
332   if (fHistos) FillHistograms(ev1);
333
334
335   // clear arrays
336   if (!fDontClearArrays) ClearArrays();
337
338   // reset TPC EP and unique identifiers for v0 cut class 
339   AliDielectronVarManager::SetTPCEventPlane(0x0);
340   if(GetHasMC()) { // only for MC needed
341     for (Int_t iCut=0; iCut<fTrackFilter.GetCuts()->GetEntries();++iCut) {
342       if ( fTrackFilter.GetCuts()->At(iCut)->IsA() == AliDielectronV0Cuts::Class() )
343         ((AliDielectronV0Cuts*)fTrackFilter.GetCuts()->At(iCut))->ResetUniqueEventNumbers();
344     }
345   }
346
347 }
348
349 //________________________________________________________________
350 void AliDielectron::ProcessMC(AliVEvent *ev1)
351 {
352   //
353   // Process the MC data
354   //
355
356   AliDielectronMC *dieMC=AliDielectronMC::Instance();
357
358   if (fHistos) FillHistogramsMC(dieMC->GetMCEvent(), ev1);
359
360   if(!fSignalsMC) return;
361   //loop over all MC data and Fill the CF container if it exist
362   if (!fCfManagerPair && !fHistoArray) return;
363   if(fCfManagerPair) fCfManagerPair->SetPdgMother(fPdgMother);
364
365   Bool_t bFillCF = (fCfManagerPair ? fCfManagerPair->GetStepForMCtruth()  : kFALSE);
366   Bool_t bFillHF = (fHistoArray    ? fHistoArray->GetStepForMCGenerated() : kFALSE);
367   if(!bFillCF && !bFillHF) return;
368
369   // signals to be studied
370   Int_t nSignals = fSignalsMC->GetEntries();
371
372   // initialize 2D arrays of labels for particles from each MC signal
373   Int_t** labels1;      // labels for particles satisfying branch 1
374   Int_t** labels2;      // labels for particles satisfying branch 2
375   Int_t** labels12;     // labels for particles satisfying both branches
376   labels1 = new Int_t*[nSignals];
377   labels2 = new Int_t*[nSignals];
378   labels12 = new Int_t*[nSignals];
379   Int_t* indexes1=new Int_t[nSignals];
380   Int_t* indexes2=new Int_t[nSignals];
381   Int_t* indexes12=new Int_t[nSignals];
382   for(Int_t isig=0;isig<nSignals;++isig) {
383     *(labels1+isig) = new Int_t[dieMC->GetNMCTracks()];
384     *(labels2+isig) = new Int_t[dieMC->GetNMCTracks()];
385     *(labels12+isig) = new Int_t[dieMC->GetNMCTracks()];
386     for(Int_t ip=0; ip<dieMC->GetNMCTracks();++ip) {
387       labels1[isig][ip] = -1;
388       labels2[isig][ip] = -1;
389       labels12[isig][ip] = -1;
390     }
391     indexes1[isig]=0;
392     indexes2[isig]=0;
393     indexes12[isig]=0;
394   }
395
396   Bool_t truth1=kFALSE;
397   Bool_t truth2=kFALSE;
398   // loop over the MC tracks
399   for(Int_t ipart=0; ipart<dieMC->GetNMCTracks(); ++ipart) {
400     for(Int_t isig=0; isig<nSignals; ++isig) {       // loop over signals
401       // Proceed only if this signal is required in the pure MC step
402       // NOTE: Some signals can be satisfied by many particles and this leads to high
403       //       computation times (e.g. secondary electrons from the GEANT transport). Be aware of this!!
404       if(!((AliDielectronSignalMC*)fSignalsMC->At(isig))->GetFillPureMCStep()) continue;
405
406       truth1 = dieMC->IsMCTruth(ipart, (AliDielectronSignalMC*)fSignalsMC->At(isig), 1);
407       truth2 = dieMC->IsMCTruth(ipart, (AliDielectronSignalMC*)fSignalsMC->At(isig), 2);
408
409       // particles satisfying both branches are treated separately to avoid double counting during pairing
410       if(truth1 && truth2) {
411         labels12[isig][indexes12[isig]] = ipart;
412         ++indexes12[isig];
413       }
414       else {
415         if(truth1) {
416           labels1[isig][indexes1[isig]] = ipart;
417           ++indexes1[isig];
418         }
419         if(truth2) {
420           labels2[isig][indexes2[isig]] = ipart;
421           ++indexes2[isig];
422         }
423       }
424     }
425   }  // end loop over MC particles
426
427   // Do the pairing and fill the CF container with pure MC info
428   for(Int_t isig=0; isig<nSignals; ++isig) {
429     // mix the particles which satisfy only one of the signal branches
430     for(Int_t i1=0;i1<indexes1[isig];++i1) {
431       for(Int_t i2=0;i2<indexes2[isig];++i2) {
432         if(bFillCF) fCfManagerPair->FillMC(labels1[isig][i1], labels2[isig][i2], isig);
433         if(bFillHF) fHistoArray->Fill(labels1[isig][i1], labels2[isig][i2], isig);
434       }
435     }
436     // mix the particles which satisfy both branches
437     for(Int_t i1=0;i1<indexes12[isig];++i1) {
438       for(Int_t i2=0; i2<i1; ++i2) {
439         if(bFillCF) fCfManagerPair->FillMC(labels12[isig][i1], labels12[isig][i2], isig);
440         if(bFillHF) fHistoArray->Fill(labels12[isig][i1], labels12[isig][i2], isig);
441       }
442     }
443   }    // end loop over signals
444
445   // release the memory
446   for(Int_t isig=0;isig<nSignals;++isig) {
447     delete [] *(labels1+isig);
448     delete [] *(labels2+isig);
449     delete [] *(labels12+isig);
450   }
451   delete [] labels1;
452   delete [] labels2;
453   delete [] labels12;
454   delete [] indexes1;
455   delete [] indexes2;
456   delete [] indexes12;
457 }
458
459 //________________________________________________________________
460 void AliDielectron::FillHistogramsTracks(TObjArray **tracks)
461 {
462   //
463   // Fill Histogram information for tracks after prefilter
464   // ignore mixed events - for prefilter, only single tracks +/- are relevant 
465   //
466   
467   TString  className,className2;
468   Double_t values[AliDielectronVarManager::kNMaxValues];
469   
470   //Fill track information, separately for the track array candidates
471   for (Int_t i=0; i<2; ++i){
472     className.Form("Pre_%s",fgkTrackClassNames[i]);
473     if (!fHistos->GetHistogramList()->FindObject(className.Data())) continue;
474     Int_t ntracks=tracks[i]->GetEntriesFast();
475     for (Int_t itrack=0; itrack<ntracks; ++itrack){
476       AliDielectronVarManager::Fill(tracks[i]->UncheckedAt(itrack), values);
477       fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
478     }
479   }
480 }
481
482
483 //________________________________________________________________
484 void AliDielectron::FillHistogramsMC(const AliMCEvent *ev, AliVEvent *ev1)
485 {
486   //
487   // Fill Histogram information for MCEvents
488   //
489
490   Double_t values[AliDielectronVarManager::kNMaxValues]={0.};
491   // Fill event information
492   AliDielectronVarManager::Fill(ev1, values);    // ESD/AOD information
493   AliDielectronVarManager::Fill(ev, values);     // MC truth info
494   if (fHistos->GetHistogramList()->FindObject("MCEvent"))
495     fHistos->FillClass("MCEvent", AliDielectronVarManager::kNMaxValues, values);
496 }
497
498
499 //________________________________________________________________
500 void AliDielectron::FillHistograms(const AliVEvent *ev, Bool_t pairInfoOnly)
501 {
502   //
503   // Fill Histogram information for tracks and pairs
504   //
505   
506   TString  className,className2;
507   Double_t values[AliDielectronVarManager::kNMaxValues]={0.};
508
509   //Fill event information
510   if (ev){
511     if (fHistos->GetHistogramList()->FindObject("Event")) {
512       fHistos->FillClass("Event", AliDielectronVarManager::kNMaxValues, AliDielectronVarManager::GetData());
513     }
514   }
515
516   //Fill track information, separately for the track array candidates
517   if (!pairInfoOnly){
518     className2.Form("Track_%s",fgkPairClassNames[1]);  // unlike sign, SE only
519     for (Int_t i=0; i<4; ++i){
520       className.Form("Track_%s",fgkTrackClassNames[i]);
521       Bool_t mergedtrkClass=fHistos->GetHistogramList()->FindObject(className2.Data())!=0x0;
522       Bool_t trkClass=fHistos->GetHistogramList()->FindObject(className.Data())!=0x0;
523       if (!trkClass && !mergedtrkClass) continue;
524       Int_t ntracks=fTracks[i].GetEntriesFast();
525       for (Int_t itrack=0; itrack<ntracks; ++itrack){
526         AliDielectronVarManager::Fill(fTracks[i].UncheckedAt(itrack), values);
527         if(trkClass)
528           fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
529         if(mergedtrkClass && i<2)
530           fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values); //only ev1
531       }
532     }
533   }
534
535   //Fill Pair information, separately for all pair candidate arrays and the legs
536   TObjArray arrLegs(100);
537   for (Int_t i=0; i<10; ++i){
538     className.Form("Pair_%s",fgkPairClassNames[i]);
539     className2.Form("Track_Legs_%s",fgkPairClassNames[i]);
540     Bool_t pairClass=fHistos->GetHistogramList()->FindObject(className.Data())!=0x0;
541     Bool_t legClass=fHistos->GetHistogramList()->FindObject(className2.Data())!=0x0;
542     if (!pairClass&&!legClass) continue;
543     Int_t ntracks=PairArray(i)->GetEntriesFast();
544     for (Int_t ipair=0; ipair<ntracks; ++ipair){
545       AliDielectronPair *pair=static_cast<AliDielectronPair*>(PairArray(i)->UncheckedAt(ipair));
546       
547       //fill pair information
548       if (pairClass){
549         AliDielectronVarManager::Fill(pair, values);
550         fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
551       }
552
553       //fill leg information, don't fill the information twice
554       if (legClass){
555         AliVParticle *d1=pair->GetFirstDaughter();
556         AliVParticle *d2=pair->GetSecondDaughter();
557         if (!arrLegs.FindObject(d1)){
558           AliDielectronVarManager::Fill(d1, values);
559           fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
560           arrLegs.Add(d1);
561         }
562         if (!arrLegs.FindObject(d2)){
563           AliDielectronVarManager::Fill(d2, values);
564           fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
565           arrLegs.Add(d2);
566         }
567       }
568     }
569     if (legClass) arrLegs.Clear();
570   }
571   
572 }
573 //________________________________________________________________
574 void AliDielectron::FillHistogramsPair(AliDielectronPair *pair,Bool_t fromPreFilter/*=kFALSE*/)
575 {
576   //
577   // Fill Histogram information for pairs and the track in the pair
578   // NOTE: in this funtion the leg information may be filled multiple
579   //       times. This funtion is used in the track rotation pairing
580   //       and those legs are not saved!
581   //
582   TString  className,className2;
583   Double_t values[AliDielectronVarManager::kNMaxValues];
584   
585   //Fill Pair information, separately for all pair candidate arrays and the legs
586   TObjArray arrLegs(100);
587   const Int_t type=pair->GetType();
588   if (fromPreFilter) {
589     className.Form("RejPair_%s",fgkPairClassNames[type]);
590     className2.Form("RejTrack_%s",fgkPairClassNames[type]);
591   } else {
592     className.Form("Pair_%s",fgkPairClassNames[type]);
593     className2.Form("Track_Legs_%s",fgkPairClassNames[type]);
594   }
595   
596   Bool_t pairClass=fHistos->GetHistogramList()->FindObject(className.Data())!=0x0;
597   Bool_t legClass=fHistos->GetHistogramList()->FindObject(className2.Data())!=0x0;
598   
599   //fill pair information
600   if (pairClass){
601     AliDielectronVarManager::Fill(pair, values);
602     fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
603   }
604
605   if (legClass){
606     AliVParticle *d1=pair->GetFirstDaughter();
607     AliDielectronVarManager::Fill(d1, values);
608     fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
609     
610     AliVParticle *d2=pair->GetSecondDaughter();
611     AliDielectronVarManager::Fill(d2, values);
612     fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
613   }
614 }
615
616 //________________________________________________________________
617 void AliDielectron::FillTrackArrays(AliVEvent * const ev, Int_t eventNr)
618 {
619   //
620   // select tracks and fill track candidate arrays
621   // eventNr = 0: First  event, use track arrays 0 and 1
622   // eventNr = 1: Second event, use track arrays 2 and 3
623   //
624
625   Int_t ntracks=ev->GetNumberOfTracks();
626
627   UInt_t selectedMask=(1<<fTrackFilter.GetCuts()->GetEntries())-1;
628   for (Int_t itrack=0; itrack<ntracks; ++itrack){
629     //get particle
630     AliVParticle *particle=ev->GetTrack(itrack);
631
632     //apply track cuts
633     UInt_t cutmask=fTrackFilter.IsSelected(particle);
634     //fill cut QA
635     if(fCutQA) fQAmonitor->FillAll(particle);
636     if(fCutQA) fQAmonitor->Fill(cutmask,particle);
637
638     if (cutmask!=selectedMask) continue;
639
640     //fill selected particle into the corresponding track arrays
641     Short_t charge=particle->Charge();
642     if (charge>0)      fTracks[eventNr*2].Add(particle);
643     else if (charge<0) fTracks[eventNr*2+1].Add(particle);
644
645   
646   }
647 }
648
649 //________________________________________________________________
650 void AliDielectron::EventPlanePreFilter(Int_t arr1, Int_t arr2, TObjArray arrTracks1, TObjArray arrTracks2, const AliVEvent *ev)
651 {
652   //
653   // Prefilter tracks and tracks from pairs
654   // Needed for rejection in the Q-Vector of the event plane
655   // remove contribution of all tracks to the Q-vector that are in invariant mass window 
656   //
657
658   AliEventplane *evplane = const_cast<AliVEvent *>(ev)->GetEventplane();
659   if(!evplane) { // nanoAODs , here we do NOT have sub event reaction planes
660     //  if(1) {
661     // get the EPselectionTask for recalculation of weighting factors
662     AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
663     AliEPSelectionTask *eptask = dynamic_cast<AliEPSelectionTask *>(man->GetTask("EventplaneSelection"));
664     if(!eptask) return;
665
666     // track mapping
667     TMap mapRemovedTracks;
668
669     // eta gap ?
670     Bool_t etagap = kFALSE;
671     for (Int_t iCut=0; iCut<fEventPlanePreFilter.GetCuts()->GetEntries();++iCut) {
672       TString cutName=fEventPlanePreFilter.GetCuts()->At(iCut)->GetName();
673       if(cutName.Contains("eta") || cutName.Contains("Eta"))  etagap=kTRUE;
674     }
675
676     Double_t cQX=0., cQY=0.;
677     // apply cuts to the tracks, e.g. etagap
678     if(fEventPlanePreFilter.GetCuts()->GetEntries()) {
679       UInt_t selectedMask=(1<<fEventPlanePreFilter.GetCuts()->GetEntries())-1;
680       Int_t ntracks=ev->GetNumberOfTracks();
681       for (Int_t itrack=0; itrack<ntracks; ++itrack){
682         AliVParticle *particle=ev->GetTrack(itrack);
683         AliVTrack *track= static_cast<AliVTrack*>(particle);
684         if (!track) continue;
685         //event plane cuts
686         UInt_t cutMask=fEventPlanePreFilter.IsSelected(track);
687         //apply cut
688         if (cutMask==selectedMask) continue;
689
690         mapRemovedTracks.Add(track,track);
691         cQX += (eptask->GetWeight(track) * TMath::Cos(2*track->Phi()));
692         cQY += (eptask->GetWeight(track) * TMath::Sin(2*track->Phi()));
693       }
694     }
695
696     // POI (particle of interest) rejection
697     Int_t pairIndex=GetPairIndex(arr1,arr2);
698
699     Int_t ntrack1=arrTracks1.GetEntriesFast();
700     Int_t ntrack2=arrTracks2.GetEntriesFast();
701     AliDielectronPair candidate;
702     candidate.SetKFUsage(fUseKF);
703
704     UInt_t selectedMask=(1<<fEventPlanePOIPreFilter.GetCuts()->GetEntries())-1;
705     for (Int_t itrack1=0; itrack1<ntrack1; ++itrack1){
706       Int_t end=ntrack2;
707       if (arr1==arr2) end=itrack1;
708       Bool_t accepted=kFALSE;
709       for (Int_t itrack2=0; itrack2<end; ++itrack2){
710         TObject *track1=arrTracks1.UncheckedAt(itrack1);
711         TObject *track2=arrTracks2.UncheckedAt(itrack2);
712         if (!track1 || !track2) continue;
713         //create the pair
714         candidate.SetTracks(static_cast<AliVTrack*>(track1), fPdgLeg1,
715                             static_cast<AliVTrack*>(track2), fPdgLeg2);
716         candidate.SetType(pairIndex);
717         candidate.SetLabel(AliDielectronMC::Instance()->GetLabelMotherWithPdg(&candidate,fPdgMother));
718
719         //event plane pair cuts
720         UInt_t cutMask=fEventPlanePOIPreFilter.IsSelected(&candidate);
721         //apply cut
722         if (cutMask==selectedMask) continue;
723
724         accepted=kTRUE;
725         //remove the tracks from the Track arrays
726         arrTracks2.AddAt(0x0,itrack2);
727       }
728       if ( accepted ) arrTracks1.AddAt(0x0,itrack1);
729     }
730     //compress the track arrays
731     arrTracks1.Compress();
732     arrTracks2.Compress();
733
734     //Modify the components: subtract the tracks
735     ntrack1=arrTracks1.GetEntriesFast();
736     ntrack2=arrTracks2.GetEntriesFast();
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       // track contribution was already removed
742       if (mapRemovedTracks.FindObject(track)) continue;
743       else mapRemovedTracks.Add(track,track);
744
745       cQX += (eptask->GetWeight(track) * TMath::Cos(2*track->Phi()));
746       cQY += (eptask->GetWeight(track) * TMath::Sin(2*track->Phi()));
747     }
748     // remove leg2 contribution
749     for (Int_t itrack=0; itrack<ntrack2; ++itrack){
750       AliVTrack *track= static_cast<AliVTrack*>(arrTracks2.UncheckedAt(itrack));
751       if (!track) continue;
752       // track contribution was already removed
753       if (mapRemovedTracks.FindObject(track)) continue;
754       else mapRemovedTracks.Add(track,track);
755
756       cQX += (eptask->GetWeight(track) * TMath::Cos(2*track->Phi()));
757       cQY += (eptask->GetWeight(track) * TMath::Sin(2*track->Phi()));
758     }
759
760     // build a corrected alieventplane using the values from the var manager
761     // these uncorrected values are filled using the stored magnitude and angle  in the header
762     TVector2 qcorr;
763     qcorr.Set(AliDielectronVarManager::GetValue(AliDielectronVarManager::kTPCxH2uc)-cQX,
764               AliDielectronVarManager::GetValue(AliDielectronVarManager::kTPCyH2uc)-cQY);
765     // fill alieventplane
766     AliEventplane cevplane;
767     cevplane.SetQVector(&qcorr);
768     AliDielectronVarManager::SetTPCEventPlane(&cevplane);
769     cevplane.SetQVector(0);
770     return;
771   } //end: nanoAODs
772   else
773     {
774     // this is done in case of ESDs or AODs
775     Bool_t isESD=(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()->IsA()==AliESDInputHandler::Class());
776     // copy event plane object
777     AliEventplane cevplane(*evplane);
778     //    Int_t nMaxID=cevplane->GetQContributionXArray()->GetSize();
779
780     TVector2 *qcorr  = cevplane.GetQVector();
781     if(!qcorr) return;
782     TVector2 *qcsub1 = 0x0;
783     TVector2 *qcsub2 = 0x0;
784
785     // eta gap ?
786     Bool_t etagap = kFALSE;
787     for (Int_t iCut=0; iCut<fEventPlanePreFilter.GetCuts()->GetEntries();++iCut) {
788       TString cutName=fEventPlanePreFilter.GetCuts()->At(iCut)->GetName();
789       if(cutName.Contains("eta") || cutName.Contains("Eta"))  etagap=kTRUE;
790     }
791
792     // subevent configuration for eta gap or LS (default is rndm)
793     if(fLikeSignSubEvents && etagap) {
794       // start with the full Qvector/event in both sub events
795       qcsub1 = new TVector2(*qcorr);
796       qcsub2 = new TVector2(*qcorr);
797       cevplane.SetQsub(qcsub1,qcsub2);
798
799       Int_t ntracks=ev->GetNumberOfTracks();
800       // track removals
801       for (Int_t itrack=0; itrack<ntracks; ++itrack){
802         AliVParticle *particle=ev->GetTrack(itrack);
803         AliVTrack *track= static_cast<AliVTrack*>(particle);
804         if (!track) continue;
805         if (track->GetID()>=0 && !isESD) continue;
806         Int_t tmpID = isESD ? track->GetID() : track->GetID()*-1 - 1;
807
808         // set contributions to zero
809         // charge sub1+ sub2-
810         if(fLikeSignSubEvents) {
811           Short_t charge=track->Charge();
812           if (charge<0) {
813             cevplane.GetQContributionXArraysub1()->SetAt(0.0, tmpID);
814             cevplane.GetQContributionYArraysub1()->SetAt(0.0, tmpID);
815           }
816           if (charge>0) {
817             cevplane.GetQContributionXArraysub2()->SetAt(0.0, tmpID);
818             cevplane.GetQContributionYArraysub2()->SetAt(0.0, tmpID);
819           }
820         }
821         // eta sub1+ sub2-
822         if(etagap) {
823           Double_t eta=track->Eta();
824           if (eta<0.0) {
825             cevplane.GetQContributionXArraysub1()->SetAt(0.0, tmpID);
826             cevplane.GetQContributionYArraysub1()->SetAt(0.0, tmpID);
827           }
828           if (eta>0.0) {
829             cevplane.GetQContributionXArraysub2()->SetAt(0.0, tmpID);
830             cevplane.GetQContributionYArraysub2()->SetAt(0.0, tmpID);
831           }
832         }
833       } // end: loop over tracks
834     } // end: sub event configuration
835
836     // apply cuts, e.g. etagap
837     if(fEventPlanePreFilter.GetCuts()->GetEntries()) {
838       UInt_t selectedMask=(1<<fEventPlanePreFilter.GetCuts()->GetEntries())-1;
839       Int_t ntracks=ev->GetNumberOfTracks();
840       for (Int_t itrack=0; itrack<ntracks; ++itrack){
841         AliVParticle *particle=ev->GetTrack(itrack);
842         AliVTrack *track= static_cast<AliVTrack*>(particle);
843         if (!track) continue;
844         if (track->GetID()>=0 && !isESD) continue;
845         Int_t tmpID = isESD ? track->GetID() : track->GetID()*-1 - 1;
846
847         //event plane cuts
848         UInt_t cutMask=fEventPlanePreFilter.IsSelected(track);
849         //apply cut
850         if (cutMask==selectedMask) continue;
851
852         // set contributions to zero
853         cevplane.GetQContributionXArray()->SetAt(0.0, tmpID);
854         cevplane.GetQContributionYArray()->SetAt(0.0, tmpID);
855         cevplane.GetQContributionXArraysub1()->SetAt(0.0, tmpID);
856         cevplane.GetQContributionYArraysub1()->SetAt(0.0, tmpID);
857         cevplane.GetQContributionXArraysub2()->SetAt(0.0, tmpID);
858         cevplane.GetQContributionYArraysub2()->SetAt(0.0, tmpID);
859       }
860     } // end: track cuts
861
862     // POI (particle of interest) rejection
863     Int_t pairIndex=GetPairIndex(arr1,arr2);
864     Int_t ntrack1=arrTracks1.GetEntriesFast();
865     Int_t ntrack2=arrTracks2.GetEntriesFast();
866     AliDielectronPair candidate;
867     candidate.SetKFUsage(fUseKF);
868
869     UInt_t selectedMask=(1<<fEventPlanePOIPreFilter.GetCuts()->GetEntries())-1;
870     for (Int_t itrack1=0; itrack1<ntrack1; ++itrack1){
871       Int_t end=ntrack2;
872       if (arr1==arr2) end=itrack1;
873       Bool_t accepted=kFALSE;
874       for (Int_t itrack2=0; itrack2<end; ++itrack2){
875         TObject *track1=arrTracks1.UncheckedAt(itrack1);
876         TObject *track2=arrTracks2.UncheckedAt(itrack2);
877         if (!track1 || !track2) continue;
878         //create the pair
879         candidate.SetTracks(static_cast<AliVTrack*>(track1), fPdgLeg1,
880                             static_cast<AliVTrack*>(track2), fPdgLeg2);
881
882         candidate.SetType(pairIndex);
883         candidate.SetLabel(AliDielectronMC::Instance()->GetLabelMotherWithPdg(&candidate,fPdgMother));
884
885         //event plane cuts
886         UInt_t cutMask=fEventPlanePOIPreFilter.IsSelected(&candidate);
887         //apply cut
888         if (cutMask==selectedMask) continue;
889
890         accepted=kTRUE;
891         //remove the tracks from the Track arrays
892         arrTracks2.AddAt(0x0,itrack2);
893       }
894       if ( accepted ) arrTracks1.AddAt(0x0,itrack1);
895     }
896     //compress the track arrays
897     arrTracks1.Compress();
898     arrTracks2.Compress();
899
900     //Modify the components: subtract the tracks
901     ntrack1=arrTracks1.GetEntriesFast();
902     ntrack2=arrTracks2.GetEntriesFast();
903     // remove leg1 contribution
904     for (Int_t itrack=0; itrack<ntrack1; ++itrack){
905       AliVTrack *track= static_cast<AliVTrack*>(arrTracks1.UncheckedAt(itrack));
906       if (!track) continue;
907       if (track->GetID()>=0 && !isESD) continue;
908       Int_t tmpID = isESD ? track->GetID() : track->GetID()*-1 - 1;
909       // set contributions to zero
910       cevplane.GetQContributionXArray()->SetAt(0.0, tmpID);
911       cevplane.GetQContributionYArray()->SetAt(0.0, tmpID);
912       cevplane.GetQContributionXArraysub1()->SetAt(0.0, tmpID);
913       cevplane.GetQContributionYArraysub1()->SetAt(0.0, tmpID);
914       cevplane.GetQContributionXArraysub2()->SetAt(0.0, tmpID);
915       cevplane.GetQContributionYArraysub2()->SetAt(0.0, tmpID);
916     }
917     // remove leg2 contribution
918     for (Int_t itrack=0; itrack<ntrack2; ++itrack){
919       AliVTrack *track= static_cast<AliVTrack*>(arrTracks2.UncheckedAt(itrack));
920       if (!track) continue;
921       if (track->GetID()>=0 && !isESD) continue;
922       Int_t tmpID = isESD ? track->GetID() : track->GetID()*-1 - 1;
923       // set contributions to zero
924       cevplane.GetQContributionXArray()->SetAt(0.0, tmpID);
925       cevplane.GetQContributionYArray()->SetAt(0.0, tmpID);
926       cevplane.GetQContributionXArraysub1()->SetAt(0.0, tmpID);
927       cevplane.GetQContributionYArraysub1()->SetAt(0.0, tmpID);
928       cevplane.GetQContributionXArraysub2()->SetAt(0.0, tmpID);
929       cevplane.GetQContributionYArraysub2()->SetAt(0.0, tmpID);
930     }
931
932     // set corrected AliEventplane and fill variables with corrected values
933     AliDielectronVarManager::SetTPCEventPlane(&cevplane);
934     delete qcsub1;
935     delete qcsub2;
936   } // end: ESD or AOD case
937
938 }
939 //________________________________________________________________
940 void AliDielectron::PairPreFilter(Int_t arr1, Int_t arr2, TObjArray &arrTracks1, TObjArray &arrTracks2)
941 {
942   //
943   // Prefilter tracks from pairs
944   // Needed for datlitz rejections
945   // remove all tracks from the Single track arrays that pass the cuts in this filter
946   //
947
948   Int_t ntrack1=arrTracks1.GetEntriesFast();
949   Int_t ntrack2=arrTracks2.GetEntriesFast();
950   AliDielectronPair candidate;
951   candidate.SetKFUsage(fUseKF);
952   // flag arrays for track removal
953   Bool_t *bTracks1 = new Bool_t[ntrack1];
954   for (Int_t itrack1=0; itrack1<ntrack1; ++itrack1) bTracks1[itrack1]=kFALSE;
955   Bool_t *bTracks2 = new Bool_t[ntrack2];
956   for (Int_t itrack2=0; itrack2<ntrack2; ++itrack2) bTracks2[itrack2]=kFALSE;
957
958   UInt_t selectedMask=(1<<fPairPreFilter.GetCuts()->GetEntries())-1;
959   UInt_t selectedMaskPair=(1<<fPairFilter.GetCuts()->GetEntries())-1;
960
961   Int_t nRejPasses = 1; //for fPreFilterUnlikeOnly and no set flag 
962   if (fPreFilterAllSigns) nRejPasses = 3;
963
964   for (Int_t iRP=0; iRP < nRejPasses; ++iRP) {
965         Int_t arr1RP=arr1, arr2RP=arr2;
966         TObjArray *arrTracks1RP=&arrTracks1;
967         TObjArray *arrTracks2RP=&arrTracks2;
968         Bool_t *bTracks1RP = bTracks1;
969         Bool_t *bTracks2RP = bTracks2;
970         switch (iRP) {
971                 case 1: arr1RP=arr1;arr2RP=arr1;
972                                 arrTracks1RP=&arrTracks1;
973                                 arrTracks2RP=&arrTracks1;
974                                 bTracks1RP = bTracks1;
975                                 bTracks2RP = bTracks1;
976                                 break;
977                 case 2: arr1RP=arr2;arr2RP=arr2;
978                                 arrTracks1RP=&arrTracks2;
979                                 arrTracks2RP=&arrTracks2;
980                                 bTracks1RP = bTracks2;
981                                 bTracks2RP = bTracks2;
982                                 break;
983                 default: ;//nothing to do
984         }
985         Int_t ntrack1RP=(*arrTracks1RP).GetEntriesFast();
986         Int_t ntrack2RP=(*arrTracks2RP).GetEntriesFast();
987
988         Int_t pairIndex=GetPairIndex(arr1RP,arr2RP);
989
990         for (Int_t itrack1=0; itrack1<ntrack1RP; ++itrack1){
991           Int_t end=ntrack2RP;
992           if (arr1RP==arr2RP) end=itrack1;
993           for (Int_t itrack2=0; itrack2<end; ++itrack2){
994                 TObject *track1=(*arrTracks1RP).UncheckedAt(itrack1);
995                 TObject *track2=(*arrTracks2RP).UncheckedAt(itrack2);
996                 if (!track1 || !track2) continue;
997                 //create the pair
998                 candidate.SetTracks(static_cast<AliVTrack*>(track1), fPdgLeg1,
999                         static_cast<AliVTrack*>(track2), fPdgLeg2);
1000
1001                 candidate.SetType(pairIndex);
1002                 candidate.SetLabel(AliDielectronMC::Instance()->GetLabelMotherWithPdg(&candidate,fPdgMother));
1003                 //relate to the production vertex
1004                 //       if (AliDielectronVarManager::GetKFVertex()) candidate.SetProductionVertex(*AliDielectronVarManager::GetKFVertex());
1005
1006                 //pair cuts
1007                 UInt_t cutMask=fPairPreFilter.IsSelected(&candidate);
1008
1009                 //apply cut
1010                 if (cutMask!=selectedMask) continue;
1011                 if (fCfManagerPair) fCfManagerPair->Fill(selectedMaskPair+1 ,&candidate);
1012                 if (fHistos) FillHistogramsPair(&candidate,kTRUE);
1013                 //set flags for track removal
1014                 bTracks1RP[itrack1]=kTRUE;
1015                 bTracks2RP[itrack2]=kTRUE;
1016           }
1017         }
1018   }
1019
1020   //remove the tracks from the Track arrays
1021   for (Int_t itrack1=0; itrack1<ntrack1; ++itrack1){
1022     if(bTracks1[itrack1]) arrTracks1.AddAt(0x0, itrack1);
1023   }
1024   for (Int_t itrack2=0; itrack2<ntrack2; ++itrack2){
1025     if(bTracks2[itrack2]) arrTracks2.AddAt(0x0, itrack2);
1026   }
1027
1028   // clean up
1029   delete [] bTracks1;
1030   delete [] bTracks2;
1031
1032   //compress the track arrays
1033   arrTracks1.Compress();
1034   arrTracks2.Compress();
1035   
1036   //apply leg cuts after the pre filter
1037   if ( fPairPreFilterLegs.GetCuts()->GetEntries()>0 ) {
1038     selectedMask=(1<<fPairPreFilterLegs.GetCuts()->GetEntries())-1;
1039     //loop over tracks from array 1
1040     for (Int_t itrack=0; itrack<arrTracks1.GetEntriesFast();++itrack){
1041       //test cuts
1042       UInt_t cutMask=fPairPreFilterLegs.IsSelected(arrTracks1.UncheckedAt(itrack));
1043       
1044       //apply cut
1045       if (cutMask!=selectedMask) arrTracks1.AddAt(0x0,itrack);;
1046     }
1047     arrTracks1.Compress();
1048     
1049     //in case of like sign don't loop over second array
1050     if (arr1==arr2) {
1051       arrTracks2=arrTracks1;
1052     } else {
1053       
1054       //loop over tracks from array 2
1055       for (Int_t itrack=0; itrack<arrTracks2.GetEntriesFast();++itrack){
1056       //test cuts
1057         UInt_t cutMask=fPairPreFilterLegs.IsSelected(arrTracks2.UncheckedAt(itrack));
1058       //apply cut
1059         if (cutMask!=selectedMask) arrTracks2.AddAt(0x0,itrack);
1060       }
1061       arrTracks2.Compress();
1062       
1063     }
1064   }
1065   //For unlike-sign monitor track-cuts:
1066   if (arr1!=arr2&&fHistos) {
1067     TObjArray *unlikesignArray[2] = {&arrTracks1,&arrTracks2};
1068     FillHistogramsTracks(unlikesignArray);
1069   }
1070 }
1071
1072 //________________________________________________________________
1073 void AliDielectron::FillPairArrays(Int_t arr1, Int_t arr2)
1074 {
1075   //
1076   // select pairs and fill pair candidate arrays
1077   //
1078
1079   TObjArray arrTracks1=fTracks[arr1];
1080   TObjArray arrTracks2=fTracks[arr2];
1081
1082   //process pre filter if set
1083   if ((!fPreFilterAllSigns) && (!fPreFilterUnlikeOnly) && ( fPairPreFilter.GetCuts()->GetEntries()>0 ))  PairPreFilter(arr1, arr2, arrTracks1, arrTracks2);
1084   
1085   Int_t pairIndex=GetPairIndex(arr1,arr2);
1086
1087   Int_t ntrack1=arrTracks1.GetEntriesFast();
1088   Int_t ntrack2=arrTracks2.GetEntriesFast();
1089
1090   AliDielectronPair *candidate=new AliDielectronPair;
1091   candidate->SetKFUsage(fUseKF);
1092
1093   UInt_t selectedMask=(1<<fPairFilter.GetCuts()->GetEntries())-1;
1094   
1095   for (Int_t itrack1=0; itrack1<ntrack1; ++itrack1){
1096     Int_t end=ntrack2;
1097     if (arr1==arr2) end=itrack1;
1098     for (Int_t itrack2=0; itrack2<end; ++itrack2){
1099       //create the pair
1100       candidate->SetTracks(static_cast<AliVTrack*>(arrTracks1.UncheckedAt(itrack1)), fPdgLeg1,
1101                              static_cast<AliVTrack*>(arrTracks2.UncheckedAt(itrack2)), fPdgLeg2);
1102       candidate->SetType(pairIndex);
1103       Int_t label=AliDielectronMC::Instance()->GetLabelMotherWithPdg(candidate,fPdgMother);
1104       candidate->SetLabel(label);
1105       if (label>-1) candidate->SetPdgCode(fPdgMother);
1106
1107       // check for gamma kf particle
1108       label=AliDielectronMC::Instance()->GetLabelMotherWithPdg(candidate,22);
1109       if (label>-1) {
1110         candidate->SetGammaTracks(static_cast<AliVTrack*>(arrTracks1.UncheckedAt(itrack1)), fPdgLeg1,
1111                                   static_cast<AliVTrack*>(arrTracks2.UncheckedAt(itrack2)), fPdgLeg2);
1112       // should we set the pdgmothercode and the label
1113       }
1114
1115       //pair cuts
1116       UInt_t cutMask=fPairFilter.IsSelected(candidate);
1117       
1118       //CF manager for the pair
1119       if (fCfManagerPair) fCfManagerPair->Fill(cutMask,candidate);
1120       //histogram array for the pair
1121       if (fHistoArray) fHistoArray->Fill(pairIndex,candidate);
1122       // cut qa
1123       if(pairIndex==AliDielectron::kEv1PM && fCutQA) {
1124         fQAmonitor->FillAll(candidate);
1125         fQAmonitor->Fill(cutMask,candidate);
1126       }
1127
1128       //apply cut
1129       if (cutMask!=selectedMask) continue;
1130
1131       //add the candidate to the candidate array 
1132       PairArray(pairIndex)->Add(candidate);
1133       //get a new candidate
1134       candidate=new AliDielectronPair;
1135           candidate->SetKFUsage(fUseKF);
1136     }
1137   }
1138   //delete the surplus candidate
1139   delete candidate;
1140 }
1141
1142 //________________________________________________________________
1143 void AliDielectron::FillPairArrayTR()
1144 {
1145   //
1146   // select pairs and fill pair candidate arrays
1147   //
1148   UInt_t selectedMask=(1<<fPairFilter.GetCuts()->GetEntries())-1;
1149   
1150   while ( fTrackRotator->NextCombination() ){
1151     AliDielectronPair candidate;
1152     candidate.SetKFUsage(fUseKF);
1153     candidate.SetTracks(&fTrackRotator->GetKFTrackP(), &fTrackRotator->GetKFTrackN(),
1154                         fTrackRotator->GetVTrackP(),fTrackRotator->GetVTrackN());
1155     candidate.SetType(kEv1PMRot);
1156     
1157     //pair cuts
1158     UInt_t cutMask=fPairFilter.IsSelected(&candidate);
1159     
1160     //CF manager for the pair
1161     if (fCfManagerPair) fCfManagerPair->Fill(cutMask,&candidate);
1162     //histogram array for the pair
1163     if (fHistoArray) fHistoArray->Fill((Int_t)kEv1PMRot,&candidate);
1164     
1165     //apply cut
1166     if (cutMask==selectedMask) {
1167      if(fHistos) FillHistogramsPair(&candidate);
1168      if(fStoreRotatedPairs) PairArray(kEv1PMRot)->Add(new AliDielectronPair(candidate));
1169     } 
1170   }
1171 }
1172
1173 //________________________________________________________________
1174 void AliDielectron::FillDebugTree()
1175 {
1176   //
1177   // Fill Histogram information for tracks and pairs
1178   //
1179   
1180   //Fill Debug tree
1181   for (Int_t i=0; i<10; ++i){
1182     Int_t ntracks=PairArray(i)->GetEntriesFast();
1183     for (Int_t ipair=0; ipair<ntracks; ++ipair){
1184       fDebugTree->Fill(static_cast<AliDielectronPair*>(PairArray(i)->UncheckedAt(ipair)));
1185     }
1186   }
1187 }
1188
1189 //________________________________________________________________
1190 void AliDielectron::SaveDebugTree()
1191 {
1192   //
1193   // delete the debug tree, this will also write the tree
1194   //
1195   if (fDebugTree) fDebugTree->DeleteStreamer();
1196 }
1197
1198
1199 //__________________________________________________________________
1200 void AliDielectron::AddSignalMC(AliDielectronSignalMC* signal) {
1201   //
1202   //  Add an MC signal to the signals list
1203   //
1204   if(!fSignalsMC) {
1205     fSignalsMC = new TObjArray();
1206     fSignalsMC->SetOwner();
1207   }
1208   fSignalsMC->Add(signal);
1209 }
1210 //________________________________________________________________
1211 void AliDielectron::FillMCHistograms(const AliVEvent *ev) {
1212   //
1213   // fill QA MC histograms for pairs and legs of all added mc signals
1214   //
1215
1216   if (!fSignalsMC) return;
1217   TString className,className2;
1218   Double_t values[AliDielectronVarManager::kNMaxValues]={0.};
1219   AliDielectronVarManager::Fill(ev, values); // get event informations
1220   //loop over all added mc signals
1221   for(Int_t isig=0; isig<fSignalsMC->GetEntries(); isig++) {
1222
1223     className.Form("Pair_%s",fSignalsMC->At(isig)->GetName());
1224     className2.Form("Track_Legs_%s",fSignalsMC->At(isig)->GetName());
1225     Bool_t pairClass=fHistos->GetHistogramList()->FindObject(className.Data())!=0x0;
1226     Bool_t legClass=fHistos->GetHistogramList()->FindObject(className2.Data())!=0x0;
1227     if(!pairClass && !legClass) return;
1228
1229     Int_t ntracks=PairArray(AliDielectron::kEv1PM)->GetEntriesFast(); // only SE +-
1230     for (Int_t ipair=0; ipair<ntracks; ++ipair){
1231       AliDielectronPair *pair=static_cast<AliDielectronPair*>(PairArray(AliDielectron::kEv1PM)->UncheckedAt(ipair));
1232
1233       Bool_t isMCtruth = AliDielectronMC::Instance()->IsMCTruth(pair, (AliDielectronSignalMC*)fSignalsMC->At(isig));
1234       if(isMCtruth) {
1235         //fill pair information
1236         if (pairClass){
1237           AliDielectronVarManager::Fill(pair, values);
1238           fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
1239         }
1240         //fill leg information, both + and - in the same histo
1241         if (legClass){
1242           AliDielectronVarManager::Fill(pair->GetFirstDaughter(),values);
1243           fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
1244           AliDielectronVarManager::Fill(pair->GetSecondDaughter(),values);
1245           fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
1246         }
1247       } //is signal
1248     } //loop: pairs
1249   } //loop: MCsignals
1250
1251 }