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