]> git.uio.no Git - u/mrichter/AliRoot.git/blob - PWGDQ/dielectron/AliAnalysisTaskMultiDielectronTG.cxx
Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[u/mrichter/AliRoot.git] / PWGDQ / dielectron / AliAnalysisTaskMultiDielectronTG.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 //                                                                       //
18 //                        Basic Analysis Task                            //
19 //                                                                       //
20 ///////////////////////////////////////////////////////////////////////////
21
22 #include <TChain.h>
23 #include <TH1D.h>
24
25 #include <AliCFContainer.h>
26 #include <AliInputEventHandler.h>
27 #include <AliESDInputHandler.h>
28 #include <AliAODInputHandler.h>
29 #include <AliAnalysisManager.h>
30 #include <AliVEvent.h>
31 #include <AliTriggerAnalysis.h>
32 #include <AliPIDResponse.h>
33 #include <AliTPCPIDResponse.h>
34
35 #include "AliDielectron.h"
36 #include "AliDielectronHistos.h"
37 #include "AliDielectronCF.h"
38 #include "AliDielectronMC.h"
39 #include "AliDielectronMixingHandler.h"
40 #include "AliAnalysisTaskMultiDielectronTG.h"
41
42 #include "AliESDtrack.h"
43 #include "AliESDtrackCuts.h"
44
45 #include "AliLog.h"
46
47 #include <vector>
48 #include <deque>
49 #include <cstdlib> 
50 #include <TRandom.h>
51
52 const char* kPairClassNamesTG[7] = {
53   "unlike",
54   "like_pp",
55   "like_ee",
56   "mixunlike_pe",
57   "mixunlike_ep",
58   "mixlike_pp",
59   "mixlike_ee"
60 };
61
62
63 ClassImp(AliDielectronSingleTG)
64 ClassImp(AliAnalysisTaskMultiDielectronTG)
65
66 //_________________________________________________________________________________
67 AliAnalysisTaskMultiDielectronTG::AliAnalysisTaskMultiDielectronTG() :
68   AliAnalysisTaskSE(),
69   fListDielectron(),
70   fListHistos(),
71   fListCF(),
72   fQAElectron(),
73   fSelectPhysics(kFALSE),
74   fTriggerMask(AliVEvent::kMB),
75   fExcludeTriggerMask(0),
76   fTriggerOnV0AND(kFALSE),
77   fRejectPileup(kFALSE),
78   fTriggerLogic(kAny),
79   fTriggerAnalysis(0x0),
80   fEventFilter(0x0),
81   fCutsMother(0x0),   
82   fEventStat(0x0),
83   fEvent(0x0),
84   fdEdXvsPt(0x0),
85   fdEdXnSigmaElecvsPt(0x0),
86   fdEdXvsPtTOF(0x0),
87   fdEdXnSigmaElecvsPtTOF(0x0),
88   fTOFbetavsPt(0x0),
89   fTOFnSigmaElecvsPt(0x0),
90   fNCrossedRowsTPC(0x0),
91   fChi2ClusTPC(0x0),
92   fRatioCrossClusTPC(0x0),
93   fVem(0x0),
94   fVep(0x0),
95   fVemtmp(0x0),
96   fVeptmp(0x0),
97   fdconvphiv(acos(-1.0)),
98   fdconvMee(100),
99   fdop(0),
100   fbz(0),
101   fdv0mixing(kTRUE),
102   fBGRejUnlike(kFALSE),
103   fBGRejLike(kTRUE)
104 {
105   //
106   // Constructor
107   //
108 }
109
110 //_________________________________________________________________________________
111 AliAnalysisTaskMultiDielectronTG::AliAnalysisTaskMultiDielectronTG(const char *name) :
112   AliAnalysisTaskSE(name),
113   fListDielectron(),
114   fListHistos(),
115   fListCF(),
116   fQAElectron(),
117   fSelectPhysics(kFALSE),
118   fTriggerMask(AliVEvent::kMB),
119   fExcludeTriggerMask(0),
120   fTriggerOnV0AND(kFALSE),
121   fRejectPileup(kFALSE),
122   fTriggerLogic(kAny),
123   fTriggerAnalysis(0x0),
124   fEventFilter(0x0),
125   fCutsMother(0x0),   
126   fEventStat(0x0),
127   fEvent(0x0),
128   fdEdXvsPt(0x0),
129   fdEdXnSigmaElecvsPt(0x0),
130   fdEdXvsPtTOF(0x0),
131   fdEdXnSigmaElecvsPtTOF(0x0),
132   fTOFbetavsPt(0x0),
133   fTOFnSigmaElecvsPt(0x0),
134   fNCrossedRowsTPC(0x0),
135   fChi2ClusTPC(0x0),
136   fRatioCrossClusTPC(0x0),
137   fVem(0x0),
138   fVep(0x0),
139   fVemtmp(0x0),
140   fVeptmp(0x0),
141   fdconvphiv(acos(-1.0)),
142   fdconvMee(100),
143   fdop(0),
144   fbz(0),
145   fdv0mixing(kTRUE),
146   fBGRejUnlike(kFALSE),
147   fBGRejLike(kTRUE)
148 {
149   //
150   // Constructor
151   //
152   DefineInput(0,TChain::Class());
153   DefineOutput(1, TList::Class());
154   DefineOutput(2, TList::Class());
155   DefineOutput(3, TH1D::Class());
156   fListHistos.SetName("Dielectron_Histos_Multi");
157   fListCF.SetName("Dielectron_CF_Multi");
158   fListDielectron.SetOwner();
159   fListHistos.SetOwner();
160   fListCF.SetOwner();
161
162   ///////////////
163   for(int i=0;i<fgkNDIE; i++){
164     for(int j=0;j<fgkNZBIN;j++){
165       for(int k=0;k<fgkNCENT;k++){
166         for(int l=0; l<fgkNRPBIN; l++){
167           fibuf[i][j][k][l] = 0;
168           fpoolp[i][j][k][l].clear();
169           fpoolm[i][j][k][l].clear();
170           for(int ib=0;ib<fgkNBUF; ib++){
171             fvep[ib][i][j][k][l].clear();
172             fvem[ib][i][j][k][l].clear();
173           }
174         }
175       }
176     }
177   }
178
179
180
181 }
182
183 //_________________________________________________________________________________
184 AliAnalysisTaskMultiDielectronTG::~AliAnalysisTaskMultiDielectronTG()
185 {
186   //
187   // Destructor
188   //
189
190   //histograms and CF are owned by the dielectron framework.
191   //however they are streamed to file, so in the first place the
192   //lists need to be owner...
193   fListHistos.SetOwner(kFALSE);
194   fListCF.SetOwner(kFALSE);
195   
196   for(int i=0;i<fgkNDIE; i++){
197     for(int j=0;j<fgkNZBIN;j++){
198       for(int k=0;k<fgkNCENT;k++){
199         for(int l=0; l<fgkNRPBIN; l++){
200           fibuf[i][j][k][l] = 0;
201           fpoolp[i][j][k][l].clear();
202           fpoolm[i][j][k][l].clear();
203           for(int ib=0;ib<fgkNBUF; ib++){
204             fvep[ib][i][j][k][l].clear();
205             fvem[ib][i][j][k][l].clear();
206           }
207         }
208       }
209     }
210   }
211 }
212 //_________________________________________________________________________________
213 void AliAnalysisTaskMultiDielectronTG::UserCreateOutputObjects()
214 {
215   //
216   // Add all histogram manager histogram lists to the output TList
217   //
218
219   if (!fListHistos.IsEmpty()||!fListCF.IsEmpty()) return; //already initialised
220
221 //   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
222 //   Bool_t isESD=man->GetInputEventHandler()->IsA()==AliESDInputHandler::Class();
223 //   Bool_t isAOD=man->GetInputEventHandler()->IsA()==AliAODInputHandler::Class();
224   
225   TIter nextDie(&fListDielectron);
226   AliDielectron *die=0;
227   while ( (die=static_cast<AliDielectron*>(nextDie())) ){
228     die->Init();
229     if (die->GetHistogramList()) fListHistos.Add(const_cast<THashList*>(die->GetHistogramList()));
230     if (die->GetCFManagerPair()) fListCF.Add(const_cast<AliCFContainer*>(die->GetCFManagerPair()->GetContainer()));
231   }
232
233   Int_t cuts=fListDielectron.GetEntries();
234   Int_t nbins=kNbinsEvent+2*cuts;
235   if (!fEventStat){
236     fEventStat=new TH1D("hEventStat","Event statistics",nbins,0,nbins);
237     fEventStat->GetXaxis()->SetBinLabel(1,"Before Phys. Sel.");
238     fEventStat->GetXaxis()->SetBinLabel(2,"After Phys. Sel.");
239
240     //default names
241     fEventStat->GetXaxis()->SetBinLabel(3,"Bin3 not used");
242     fEventStat->GetXaxis()->SetBinLabel(4,"Bin4 not used");
243     fEventStat->GetXaxis()->SetBinLabel(5,"Bin5 not used");
244     
245     if(fTriggerOnV0AND) fEventStat->GetXaxis()->SetBinLabel(3,"V0and triggers");
246     if (fEventFilter) fEventStat->GetXaxis()->SetBinLabel(4,"After Event Filter");
247     if (fRejectPileup) fEventStat->GetXaxis()->SetBinLabel(5,"After Pileup rejection");
248     
249     for (Int_t i=0; i<cuts; ++i){
250       fEventStat->GetXaxis()->SetBinLabel((kNbinsEvent+1)+2*i,Form("#splitline{1 candidate}{%s}",fListDielectron.At(i)->GetName()));
251       fEventStat->GetXaxis()->SetBinLabel((kNbinsEvent+2)+2*i,Form("#splitline{With >1 candidate}{%s}",fListDielectron.At(i)->GetName()));
252     }
253   }
254
255   if (!fTriggerAnalysis) fTriggerAnalysis=new AliTriggerAnalysis;
256   fTriggerAnalysis->EnableHistograms();
257   fTriggerAnalysis->SetAnalyzeMC(AliDielectronMC::Instance()->HasMC());
258   
259
260   
261   int nbinx=400;
262   float maxx=20;
263   float minx=0.2;
264   float binw = (TMath::Log(maxx)-TMath::Log(minx))/nbinx;
265   double xbin[401];
266   for(int ii=0;ii<nbinx+1;ii++){
267     xbin[ii] = TMath::Exp(TMath::Log(minx) + 0.5*binw+binw*ii);
268   }
269
270   
271   fQAElectron = new TList();
272   fQAElectron->SetName("QAElectron");
273   fQAElectron->SetOwner();
274
275
276   fEvent = new TH1D("Event","centrality",   100,0,100);
277   fQAElectron->Add(fEvent);
278   fdEdXvsPt = new TH2D("dEdXvsPt","dE/dX vs. PT of TPC", nbinx, xbin, 2000,0,200);
279   fQAElectron->Add(fdEdXvsPt);
280   fdEdXnSigmaElecvsPt = new TH2D("fdEdXnSigmaElecvsPt"," dE/dX normalized to electron vs. pT of TPC",
281                                  nbinx, xbin, 2000, -10, 10);
282   fQAElectron->Add(fdEdXnSigmaElecvsPt);
283
284   fdEdXvsPtTOF = new TH2D("dEdXvsPtTOF","dE/dX vs. PT of TPC", nbinx, xbin, 2000,0,200);
285   fQAElectron->Add(fdEdXvsPtTOF);
286   fdEdXnSigmaElecvsPtTOF = new TH2D("fdEdXnSigmaElecvsPtTOF"," dE/dX normalized to electron vs. pT of TPC",
287                                  nbinx, xbin, 2000, -10, 10);
288   fQAElectron->Add(fdEdXnSigmaElecvsPtTOF);
289
290
291
292   fTOFbetavsPt = new TH2D("fTOFbetavsPt","TOF beta vs. p", 400, 0, 20, 1200, 0, 1.2);
293   fQAElectron->Add(fTOFbetavsPt);
294   fTOFnSigmaElecvsPt = new TH2D("fTOFnSigmaElecvsPt","TOF nsigma for electron", 400, 0, 20, 2000, -10, 10);
295   fQAElectron->Add(fTOFnSigmaElecvsPt);
296
297   fNCrossedRowsTPC = new TH2F("fNCrossedRowsTPC", "TPC nCrossed Rows vs. pT", 200, 0, 20, 200, 0, 200);
298   fQAElectron->Add(fNCrossedRowsTPC);
299   fChi2ClusTPC = new TH2F("fChi2ClusTPC", "hChi2ClusTPC vs. pT", 200, 0, 20, 200, 0, 10);
300   fQAElectron->Add(fChi2ClusTPC);
301   
302   fRatioCrossClusTPC = new TH2F("fRatioCrossClusTPC", "hRatioCrossClusTPC vs. pT", 200, 0, 20, 200, 0, 10);     
303   fQAElectron->Add(fRatioCrossClusTPC);
304
305   fListHistos.Add(fQAElectron);
306
307   fListHistos.SetOwner();  
308   
309   PostData(1, &fListHistos);
310   PostData(2, &fListCF);
311   PostData(3, fEventStat);
312
313   fCutsMother = new AliESDtrackCuts;
314   fCutsMother->SetDCAToVertex2D(kTRUE);
315   fCutsMother->SetMaxDCAToVertexZ(3.0);
316   fCutsMother->SetMaxDCAToVertexXY(1.0);
317   fCutsMother->SetPtRange(  0.05 , 200.0);
318   fCutsMother->SetEtaRange( -0.84 , 0.84 );
319   fCutsMother->SetAcceptKinkDaughters(kFALSE);
320   fCutsMother->SetRequireITSRefit(kTRUE);
321   fCutsMother->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kFirst);
322   fCutsMother->SetMinNClustersITS(3);
323   fCutsMother->SetRequireTPCRefit(kTRUE);
324
325
326   AliInfo(Form("PairCutType = %d %d %d %d %d",
327                fRejectPairFlag[0],
328                fRejectPairFlag[1],
329                fRejectPairFlag[2],
330                fRejectPairFlag[3],
331                fRejectPairFlag[4]));
332
333 }
334
335 //_________________________________________________________________________________
336 void AliAnalysisTaskMultiDielectronTG::UserExec(Option_t *)
337 {
338   //
339   // Main loop. Called for every event
340   //
341
342   if (fListHistos.IsEmpty()&&fListCF.IsEmpty()) return;
343
344   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
345   Bool_t isESD=man->GetInputEventHandler()->IsA()==AliESDInputHandler::Class();
346   Bool_t isAOD=man->GetInputEventHandler()->IsA()==AliAODInputHandler::Class();
347   
348   AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
349   if (!inputHandler) return;
350   
351 //   AliPIDResponse *pidRes=inputHandler->GetPIDResponse();
352   if ( inputHandler->GetPIDResponse() ){
353     // for the 2.76 pass2 MC private train. Together with a sigma shift of -0.169
354 //    pidRes->GetTPCResponse().SetSigma(4.637e-3,2.41332105409873257e+04);
355     AliDielectronVarManager::SetPIDResponse( inputHandler->GetPIDResponse() );
356   } else {
357     AliFatal("This task needs the PID response attached to the input event handler!");
358   }
359   
360   // Was event selected ?
361   ULong64_t isSelected = AliVEvent::kAny;
362   Bool_t isRejected = kFALSE;
363   if( fSelectPhysics && inputHandler){
364     if((isESD && inputHandler->GetEventSelection()) || isAOD){
365       isSelected = inputHandler->IsEventSelected();
366       if (fExcludeTriggerMask && (isSelected&fExcludeTriggerMask)) isRejected=kTRUE;
367       if (fTriggerLogic==kAny) isSelected&=fTriggerMask;
368       else if (fTriggerLogic==kExact) isSelected=((isSelected&fTriggerMask)==fTriggerMask);
369     }
370    }
371  
372  
373   //Before physics selection
374   fEventStat->Fill(kAllEvents);
375   if (isSelected==0||isRejected) {
376     PostData(3,fEventStat);
377     return;
378   }
379   //after physics selection
380   fEventStat->Fill(kSelectedEvents);
381
382   //V0and
383   if(fTriggerOnV0AND){
384     if(isESD){if (!fTriggerAnalysis->IsOfflineTriggerFired(static_cast<AliESDEvent*>(InputEvent()), AliTriggerAnalysis::kV0AND))
385             return;}
386     if(isAOD){if(!((static_cast<AliAODEvent*>(InputEvent()))->GetVZEROData()->GetV0ADecision() == AliVVZERO::kV0BB &&
387                    (static_cast<AliAODEvent*>(InputEvent()))->GetVZEROData()->GetV0CDecision() == AliVVZERO::kV0BB) )
388         return;}
389   }
390   
391
392   fEventStat->Fill(kV0andEvents);
393
394   //Fill Event histograms before the event filter
395   TIter nextDie(&fListDielectron);
396   AliDielectron *die=0;
397   Double_t values[AliDielectronVarManager::kNMaxValues]={0};
398   Double_t valuesMC[AliDielectronVarManager::kNMaxValues]={0};
399   AliDielectronVarManager::SetEvent(InputEvent());
400   AliDielectronVarManager::Fill(InputEvent(),values);
401   AliDielectronVarManager::Fill(InputEvent(),valuesMC);
402   Bool_t hasMC=AliDielectronMC::Instance()->HasMC();
403   if (hasMC) {
404     if (AliDielectronMC::Instance()->ConnectMCEvent())
405       AliDielectronVarManager::Fill(AliDielectronMC::Instance()->GetMCEvent(),valuesMC);
406   }
407
408   while ( (die=static_cast<AliDielectron*>(nextDie())) ){
409     AliDielectronHistos *h=die->GetHistoManager();
410     if (h){
411       if (h->GetHistogramList()->FindObject("Event_noCuts"))
412         h->FillClass("Event_noCuts",AliDielectronVarManager::kNMaxValues,values);
413       if (hasMC && h->GetHistogramList()->FindObject("MCEvent_noCuts"))
414         h->FillClass("Event_noCuts",AliDielectronVarManager::kNMaxValues,valuesMC);
415     }
416   }
417   nextDie.Reset();
418   
419   //event filter
420   if (fEventFilter) {
421     if (!fEventFilter->IsSelected(InputEvent())) return;
422   }
423   fEventStat->Fill(kFilteredEvents);
424   
425   //pileup
426   if (fRejectPileup){
427     if (InputEvent()->IsPileupFromSPD(3,0.8,3.,2.,5.)) return;
428   }
429   fEventStat->Fill(kPileupEvents);
430   
431   //bz for AliKF
432   fbz = InputEvent()->GetMagneticField();
433   AliKFParticle::SetField( fbz );
434   AliDielectronPID::SetCorrVal((Double_t)InputEvent()->GetRunNumber());
435   
436   //Process event in all AliDielectron instances
437   //   TIter nextDie(&fListDielectron);
438   //   AliDielectron *die=0;
439   Int_t idie=0;
440   while ( (die=static_cast<AliDielectron*>(nextDie())) ){
441     //AliInfo(Form(" **************** die->Process(InputEvent()) : %d **************************", idie));
442     die->SetDontClearArrays(kTRUE);
443     die->Process(InputEvent());
444     if (die->HasCandidates()){
445       Int_t ncandidates=die->GetPairArray(1)->GetEntriesFast();
446       if (ncandidates==1) fEventStat->Fill((kNbinsEvent)+2*idie);
447       else if (ncandidates>1) fEventStat->Fill((kNbinsEvent+1)+2*idie);
448     }
449
450
451     AliDielectronVarManager::Fill(InputEvent(), fgValues);
452     for (Int_t ii=0; ii<2; ++ii){ 
453       TObjArray *obj = (TObjArray*)die->GetTrackArray(ii);
454       Int_t ntracks=obj->GetEntriesFast();
455       //AliInfo(Form(" ************** # of tracks = %d", ntracks));
456       for (Int_t itrack=0; itrack<ntracks; ++itrack){
457         
458         ////////////////////////////////////////////////////////////////////
459         AliDielectronVarManager::Fill(obj->UncheckedAt(itrack), fgValues);
460         ////////////////////////////////////////////////////////////////////
461
462         if(fgValues[AliDielectronVarManager::kCharge]>0){
463           fVeptmp.push_back(new  AliDielectronSingleTG(1, 
464                                                        fgValues[AliDielectronVarManager::kCentrality],
465                                                        fgValues[AliDielectronVarManager::kXv],
466                                                        fgValues[AliDielectronVarManager::kYv],
467                                                        fgValues[AliDielectronVarManager::kZv],
468                                                        fgValues[AliDielectronVarManager::kPx],
469                                                        fgValues[AliDielectronVarManager::kPy],
470                                                        fgValues[AliDielectronVarManager::kPz],
471                                                        fgValues[AliDielectronVarManager::kPt],
472                                                        fgValues[AliDielectronVarManager::kEta],
473                                                        fgValues[AliDielectronVarManager::kPhi],
474                                                        fgValues[AliDielectronVarManager::kTheta],
475                                                        1, 1, static_cast<AliVTrack*>(obj->UncheckedAt(itrack)))
476                             );
477         }else if(fgValues[AliDielectronVarManager::kCharge]<0){
478           fVemtmp.push_back(new  AliDielectronSingleTG(-1, 
479                                                        fgValues[AliDielectronVarManager::kCentrality],
480                                                        fgValues[AliDielectronVarManager::kXv],
481                                                        fgValues[AliDielectronVarManager::kYv],
482                                                        fgValues[AliDielectronVarManager::kZv],
483                                                        fgValues[AliDielectronVarManager::kPx],
484                                                        fgValues[AliDielectronVarManager::kPy],
485                                                        fgValues[AliDielectronVarManager::kPz],
486                                                        fgValues[AliDielectronVarManager::kPt],
487                                                        fgValues[AliDielectronVarManager::kEta],
488                                                        fgValues[AliDielectronVarManager::kPhi],
489                                                        fgValues[AliDielectronVarManager::kTheta],
490                                                        1, 1, static_cast<AliVTrack*>(obj->UncheckedAt(itrack)))
491                             );
492         }
493       }
494     }
495     //AliInfo(Form("size of e and p = %d %d", (int)fVeptmp.size(), (int)fVemtmp.size()));
496     
497
498     CheckGhostPairs(fVeptmp);
499     CheckGhostPairs(fVemtmp);
500     if(fRejectPairFlag[idie]==1 || fRejectPairFlag[idie]==2){
501       RejectPairs(fVeptmp, fVemtmp, idie);
502     }
503     RandomizePool(fVeptmp, fVemtmp);    
504     CalcPair(fVep, fVem, die, idie);
505
506     //    AliInfo(Form("size of e and p (after) = %d %d", (int)fVep.size(), (int)fVem.size()));
507
508     double dwcent = 100.0/fgkNCENT;
509     double dwiz = 20.0/fgkNZBIN;
510     double dwrp = acos(-1.0)/fgkNRPBIN;
511
512     int icent = (int)(fgValues[AliDielectronVarManager::kCentrality]/dwcent);
513     int izbin = (int)((fgValues[AliDielectronVarManager::kZvPrim]+10)/dwiz);
514     int irp = (int)((fgValues[AliDielectronVarManager::kV0ACrpH2])/dwrp);
515     
516     if(icent<0) icent=0;
517     if(icent>=fgkNCENT) icent=fgkNCENT-1;
518     if(izbin<0) izbin=0;
519     if(izbin>=fgkNZBIN) izbin=fgkNZBIN-1;
520     if(irp<0) irp=0;
521     if(irp>=fgkNRPBIN) irp=fgkNRPBIN-1;
522     
523     fvep[fibuf[idie][izbin][icent][irp]][idie][izbin][icent][irp].clear();
524     for(int iep = 0; iep<(int)fVep.size();iep++) {
525       fvep[fibuf[idie][izbin][icent][irp]][idie][izbin][icent][irp].push_back(fVep[iep]);
526       fpoolp[idie][izbin][icent][irp].push_back(fVep[iep]);
527       if(fpoolp[idie][izbin][icent][irp].size()>fgkMAXPOOL) {
528         fpoolp[idie][izbin][icent][irp].pop_front();
529       }
530     }
531     fvem[fibuf[idie][izbin][icent][irp]][idie][izbin][icent][irp].clear();
532     for(int iem = 0; iem<(int)fVem.size();iem++) {
533       fvem[fibuf[idie][izbin][icent][irp]][idie][izbin][icent][irp].push_back(fVem[iem]);
534       fpoolm[idie][izbin][icent][irp].push_back(fVem[iem]);
535       if(fpoolm[idie][izbin][icent][irp].size()>fgkMAXPOOL) {
536         fpoolm[idie][izbin][icent][irp].pop_front();
537       }
538     }
539
540
541     fibuf[idie][izbin][icent][irp]++;
542     if(fibuf[idie][izbin][icent][irp]>= fgkNBUF) fibuf[idie][izbin][icent][irp]=0; 
543
544
545     fVeptmp.clear();
546     fVemtmp.clear();
547     fVep.clear();
548     fVem.clear();
549
550
551     ++idie;
552   }
553
554
555   AliESDEvent *fESD = dynamic_cast<AliESDEvent*>(InputEvent());
556   fEvent->Fill(values[AliDielectronVarManager::kCentrality]);
557   for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
558     AliESDtrack* track = fESD->GetTrack(iTracks);
559     if (!track) {
560       Printf("ERROR: Could not receive track %d", iTracks);
561       continue;
562     }
563     if(!fCutsMother->AcceptTrack(track)) continue;
564     fdEdXvsPt->Fill(track->GetTPCmomentum(), track->GetTPCsignal());
565     fdEdXnSigmaElecvsPt->Fill(track->GetTPCmomentum(), 
566                               AliDielectronVarManager::GetPIDResponse()->NumberOfSigmasTPC(track,
567                                                                                       AliPID::kElectron)
568                               -AliDielectronPID::GetCorrVal());
569     /// for beta caliculaton 
570     Double_t l = track->GetIntegratedLength();  // cm
571     Double_t t = track->GetTOFsignal();
572     Double_t t0 = AliDielectronVarManager::GetPIDResponse()->GetTOFResponse().GetTimeZero(); // ps
573     Double_t beta = 0;
574     if( (l < 360. || l > 800.) || (t <= 0.) || (t0 >999990.0) ) {
575       beta=-9999;
576     }
577     else {
578       t -= t0; // subtract the T0
579       l *= 0.01;  // cm ->m
580       t *= 1e-12; //ps -> s
581     
582       Double_t v = l / t;
583       beta = v / TMath::C();
584     }
585
586     fTOFbetavsPt->Fill(track->GetTPCmomentum(), beta);
587     fTOFnSigmaElecvsPt->Fill(track->GetTPCmomentum(), 
588                              AliDielectronVarManager::GetPIDResponse()->NumberOfSigmasTOF(track,
589                                                                                      AliPID::kElectron));
590     ////rough PID is required 
591     if( fabs(AliDielectronVarManager::GetPIDResponse()->NumberOfSigmasTOF(track, AliPID::kElectron))<3){
592       fdEdXvsPtTOF->Fill(track->GetTPCmomentum(), track->GetTPCsignal());
593       fdEdXnSigmaElecvsPtTOF->Fill(track->GetTPCmomentum(), 
594                                    AliDielectronVarManager::GetPIDResponse()->NumberOfSigmasTPC(track,
595                                                                                                 AliPID::kElectron)
596                                    -AliDielectronPID::GetCorrVal());
597
598       
599       if(track->GetTPCsignal()>70 && track->GetTPCsignal()<90){
600         fNCrossedRowsTPC->Fill(track->GetTPCmomentum(), track->GetTPCCrossedRows());
601         //fChi2ClusTPC->Fill(track->GetTPCmomentum(), track->GetTPCchi2()/track->GetTPCNcls());
602         fChi2ClusTPC->Fill(track->GetTPCmomentum(), track->GetTPCchi2());
603         fRatioCrossClusTPC->Fill(track->GetTPCmomentum(), track->GetTPCCrossedRows()/track->GetTPCNclsF());
604       }
605     }
606   }
607   
608   PostData(1, &fListHistos);
609   PostData(2, &fListCF);
610   PostData(3,fEventStat);
611 }
612
613 //_________________________________________________________________________________
614 void AliAnalysisTaskMultiDielectronTG::FinishTaskOutput()
615 {
616   //
617   // Write debug tree
618   //
619   TIter nextDie(&fListDielectron);
620   AliDielectron *die=0;
621   while ( (die=static_cast<AliDielectron*>(nextDie())) ){
622     die->SaveDebugTree();
623     AliDielectronMixingHandler *mix=die->GetMixingHandler();
624 //    printf("\n\n\n===============\ncall mix in Terminate: %p (%p)\n=================\n\n",mix,die);
625     if(!mix) continue;
626     for (Int_t ipool=0; ipool<mix->GetNumberOfBins(); ++ipool){
627       mix->MixRemaining(die, ipool);
628     }
629   }
630   PostData(1, &fListHistos);
631   PostData(2, &fListCF);
632 }
633
634 //_________________________________________________________________________________
635 void AliAnalysisTaskMultiDielectronTG::CheckGhostPairs(vector<AliDielectronSingleTG*> e1)
636 {
637   ////// Check whether the pairs are in ghost pairs 
638   /////  ghost means that one of the pairs is highly associated but reconstructed as a true particle
639
640
641   bool reject = false;
642   if(e1.size()>1){
643     for(int i1=0; i1<(int)e1.size(); i1++){
644       reject = false;
645       for(int i2=i1+1; i2<(int)e1.size(); i2++){
646         if( fabs(e1[i1]->Phi() - e1[i2]->Phi())<0.01 && 
647             fabs(e1[i1]->Eta() - e1[i2]->Eta())<0.005
648             ){
649           reject = true;
650           e1[i2]->SetGstFlag(0);
651         }
652       }
653       if(reject==true)e1[i1]->SetGstFlag(0);
654     }
655   }
656 }
657
658 //_________________________________________________________________________________
659 Bool_t AliAnalysisTaskMultiDielectronTG::CheckGhost(vector<AliDielectronSingleTG*> e1, vector<AliDielectronSingleTG*> e2)
660 {
661   ////// To be sure whether there are no ghost pairs in h event mixing 
662
663   if(e1.size()>0 && e2.size()>0){
664     for(int i1=0; i1<(int)e1.size(); i1++){
665       for(int i2=0; i2<(int)e2.size(); i2++){
666         if( fabs(e1[i1]->Phi() - e2[i2]->Phi())<0.01 && 
667             fabs(e1[i1]->Eta() - e2[i2]->Eta())<0.005
668             ){
669           return true;
670         }
671       }
672     }
673   }
674   return false;
675 }
676
677 //_________________________________________________________________________________
678 void AliAnalysisTaskMultiDielectronTG::RandomizePool(vector<AliDielectronSingleTG*> e1, vector<AliDielectronSingleTG*> e2)
679 {
680   
681   ///// Randomize the pool constent to cancel the filling scheme of the single tracks
682   //// 
683
684
685   int size1 = e1.size();
686   int usedindex[1000];
687   for(int i=0;i<1000;i++){
688     usedindex[i] = -1;
689   }
690   for(int i=0;i<size1;i++){
691     usedindex[i] = 0;
692   }
693
694   for(int i=0;i<size1;i++){
695     int j = (int)(gRandom->Uniform(0,size1));
696     while(usedindex[j]==1){
697       j = (int)(gRandom->Uniform(0,size1));
698     }
699     if( (e1[j]->GetGstFlag()==1) &&
700         (e1[j]->GetConvFlag()==1)
701         ){
702       fVep.push_back(e1[j]);
703     }
704     usedindex[j] = 1;
705   }
706   
707
708   int size2 = e2.size();
709   for(int i=0;i<1000;i++){
710     usedindex[i] = -1;
711   }
712   for(int i=0;i<size2;i++){
713     usedindex[i] = 0;
714   }
715
716   for(int i=0;i<size2;i++){
717     int j = (int)(gRandom->Uniform(0,size2));
718     while(usedindex[j]==1){
719       j = (int)(gRandom->Uniform(0,size2));
720     }
721     if( (e2[j]->GetGstFlag()==1) &&
722         (e2[j]->GetConvFlag()==1)
723         ){
724       fVem.push_back(e2[j]);
725     }
726     usedindex[j] = 1;
727   }
728 }
729
730
731 //_________________________________________________________________________________
732 void AliAnalysisTaskMultiDielectronTG::CalcPair(vector<AliDielectronSingleTG*> ve1, 
733                                                 vector<AliDielectronSingleTG*> ve2, AliDielectron *die, Int_t idie)
734 {
735
736   //
737   // main routine for the pair procesing 
738   //
739
740
741   for(int i1=0; i1<(int)ve1.size(); i1++){
742     for(int i2=0; i2<(int)ve2.size(); i2++){
743       FillPair(ve1[i1], ve2[i2], 0, die, idie);  
744     }    
745   }
746
747   for(int i1=0; i1<(int)ve1.size(); i1++){
748     for(int i2=i1+1; i2<(int)ve1.size(); i2++){
749       FillPair(ve1[i1], ve1[i2], 1, die, idie);  
750     }    
751   }
752
753   for(int i1=0; i1<(int)ve2.size(); i1++){
754     for(int i2=i1+1; i2<(int)ve2.size(); i2++){
755       FillPair(ve2[i1], ve2[i2], 2, die, idie );  
756     }    
757   }
758
759
760   double dwcent = 100.0/fgkNCENT;
761   double dwiz = 20.0/fgkNZBIN;
762   double dwrp = acos(-1.0)/fgkNRPBIN;
763
764   int icent = (int)(fgValues[AliDielectronVarManager::kCentrality]/dwcent);
765   int izbin = (int)((fgValues[AliDielectronVarManager::kZvPrim]+10)/dwiz);
766   int irp = (int)((fgValues[AliDielectronVarManager::kV0ACrpH2])/dwrp);
767
768   if(icent<0) icent=0;
769   if(icent>=fgkNCENT) icent=fgkNCENT-1;
770   if(izbin<0) izbin=0;
771   if(izbin>=fgkNZBIN) izbin=fgkNZBIN-1;
772   if(irp<0) irp=0;
773   if(irp>=fgkNRPBIN) irp=fgkNRPBIN-1;
774
775   int nmixed;
776   if(ve1.size()>0) {
777     //
778     // Now mixed event for +- pairs
779     //
780     nmixed = 0;
781     for(int ibuf=0;(nmixed<fgkNMix);ibuf++) {
782       int ntry = 0;
783       while((fBGRejUnlike && CheckGhost(ve1, fvem[ibuf][idie][izbin][icent][irp])) &&  ntry<fgkMAXTRY) {
784         ReshuffleBuffer(fvem[ibuf][idie][izbin][icent][irp],fpoolm[idie][izbin][icent][irp]);
785         ntry++;
786       }
787       for(int i1=0; i1<(int)ve1.size(); i1++){
788         for(int i2=0; i2<(int)fvem[ibuf][idie][izbin][icent][irp].size(); i2++){
789           FillPair(ve1[i1],fvem[ibuf][idie][izbin][icent][irp][i2], 3, die, idie);
790         }
791       }
792       ++nmixed;
793     }//for(ibuf)
794   }
795   if(ve2.size()>0) {
796     //
797     // Now mixed event for -+ pairs
798     //
799     nmixed = 0;
800     for(int ibuf=0;(nmixed<fgkNMix);ibuf++) {
801       int ntry = 0;
802       while((fBGRejUnlike && CheckGhost(ve2, fvep[ibuf][idie][izbin][icent][irp])) &&  ntry<fgkMAXTRY) {
803         ReshuffleBuffer(fvep[ibuf][idie][izbin][icent][irp],fpoolp[idie][izbin][icent][irp]);
804         ntry++;
805       }
806       for(int i1=0; i1<(int)ve2.size(); i1++){
807         for(int i2=0; i2<(int)fvep[ibuf][idie][izbin][icent][irp].size(); i2++){
808           FillPair(fvep[ibuf][idie][izbin][icent][irp][i2],ve2[i1],4, die, idie);
809         }
810       }
811       ++nmixed;
812     }//for(ibuf)
813   }
814
815   if(ve1.size()>0) {
816     //
817     // Now mixed event for ++ pairs
818     //
819     nmixed = 0;
820     for(int ibuf=0;(nmixed<fgkNMix);ibuf++) {
821       int ntry = 0;
822       while((fBGRejLike && CheckGhost(ve1, fvep[ibuf][idie][izbin][icent][irp])) &&  ntry<fgkMAXTRY) {
823         ReshuffleBuffer(fvep[ibuf][idie][izbin][icent][irp],fpoolp[idie][izbin][icent][irp]);
824         ntry++;
825       }
826       for(int i1=0; i1<(int)ve1.size(); i1++){
827         for(int i2=0;i2<(int)fvep[ibuf][idie][izbin][icent][irp].size(); i2++){
828           FillPair(ve1[i1],fvep[ibuf][idie][izbin][icent][irp][i2], 5, die, idie);
829         }
830       }
831       ++nmixed;
832     }//for(ibuf)
833   }
834
835   if(ve2.size()>0) {
836     //
837     // Now mixed event for +- pairs
838     //
839     nmixed = 0;
840     for(int ibuf=0;(nmixed<fgkNMix);ibuf++) {
841       int ntry = 0;
842       while((fBGRejLike && CheckGhost(ve2, fvem[ibuf][idie][izbin][icent][irp])) &&  ntry<fgkMAXTRY) {
843         ReshuffleBuffer(fvem[ibuf][idie][izbin][icent][irp],fpoolm[idie][izbin][icent][irp]);
844         ntry++;
845       }
846       for(int i1=0; i1<(int)ve2.size(); i1++){
847         for(int i2=0; i2<(int)fvem[ibuf][idie][izbin][icent][irp].size(); i2++){
848           FillPair(ve2[i1],fvem[ibuf][idie][izbin][icent][irp][i2],6, die, idie);
849         }
850       }
851       ++nmixed;
852     }//for(ibuf)
853   }
854
855 }
856
857
858 //_________________________________________________________________________________
859 void AliAnalysisTaskMultiDielectronTG::FillPair(AliDielectronSingleTG *iep, 
860                                                 AliDielectronSingleTG *iem, int type, AliDielectron *die, Int_t idie)
861 {                      
862
863   //
864   // main routine for filling kinematics of pairs
865   //
866
867
868   
869   double dmass, dphiv, dpxpair, dpypair, dpzpair;
870   double dptpair, depair, dphipair, detapair, dcos, dpsi;
871
872   CalcVars(iep, iem, dmass, dphiv, dpxpair, dpypair, dpzpair, 
873             dptpair, depair, dphipair, detapair, dcos, dpsi);
874
875
876   double dopeningangle =  -9999;
877   double dv0mass =  -9999;
878   double dv0pxpair = -9999;
879   double dv0pypair = -9999;
880   double dv0pzpair = -9999;
881   double dv0ptpair = -9999;
882   double dv0epair = -9999;
883   double dv0xvpair = -9999;
884   double dv0yvpair = -9999;
885   double dv0zvpair = -9999;
886   double dv0phipair = -9999;
887   double dv0etapair = -9999;
888   double dv0rpair = -9999;
889   double dpsipair =  -9999;
890   double dphivpair =  -9999;
891
892   ////////////////////////////
893   ///// calculate v0 ////////
894   ///////////////////////////
895   Bool_t v0OFF=kFALSE;
896   /// for the moment, this doesn't work for MC
897   if(fdv0mixing == kFALSE && (type==3 || type==4 || type==5 || type==6)){
898     v0OFF = kTRUE;
899   }
900   if(die->GetHasMC()==kTRUE && (type==3 || type==4 || type==5 || type==6)){
901     v0OFF = kTRUE;
902   }
903   if(type==0 || type==1 || type==2){
904     v0OFF = kFALSE;
905   }
906   
907
908   if(v0OFF==kFALSE){
909     AliVTrack* trackob1= iep->GetTrack();    
910     AliVTrack* trackob2= iem->GetTrack();    
911
912     if(!trackob1 || !trackob2){
913       return; 
914     }
915
916     AliDielectronPair candidate;
917     candidate.SetTracks(trackob1, (int)(11*iep->Charge()), 
918                         trackob2, (int)(11*iem->Charge()));
919     
920     if(trackob1==trackob2){
921       AliInfo("Same Objects!!");
922       return; 
923     }
924
925     //const AliKFParticle &kfPairLeg1 = candidate.GetKFFirstDaughter();
926     //const AliKFParticle &kfPairLeg2 = candidate.GetKFSecondDaughter();
927     
928     const AliKFParticle &kfPair = candidate.GetKFParticle();
929
930     /*
931     dv0mass = candidate.M();
932     dv0pxpair = candidate.Px();
933     dv0pypair = candidate.Py();
934     dv0pzpair = candidate.Pz();
935     dv0ptpair = candidate.Pt();
936     dv0epair = candidate.E();
937     dv0xvpair = candidate.Xv();
938     dv0yvpair = candidate.Yv();
939     dv0zvpair = candidate.Zv();
940     dv0phipair = candidate.Phi();
941     dv0etapair = candidate.Eta();
942     */
943     dv0mass = kfPair.GetMass();
944     dv0pxpair = kfPair.GetPx();
945     dv0pypair = kfPair.GetPy();
946     dv0pzpair = kfPair.GetPz();
947     dv0ptpair = kfPair.GetPt();
948     dv0epair = kfPair.GetE();
949     dv0xvpair = kfPair.GetX();
950     dv0yvpair = kfPair.GetY();
951     dv0zvpair = kfPair.GetZ();
952     dv0phipair = kfPair.GetPhi();
953     dv0etapair = kfPair.GetEta();
954     dv0rpair =  kfPair.GetR();
955
956     dopeningangle = candidate.OpeningAngle();    
957     dpsipair = candidate.PsiPair(fbz);
958     dphivpair = candidate.PhivPair(fbz);
959
960
961   }
962
963   Double_t values[AliDielectronVarManager::kNMaxValues];
964   TString  className1;
965   TString  className2;
966   className1.Form("MyPair_%s",kPairClassNamesTG[type]);
967   className2.Form("MyPairV0_%s",kPairClassNamesTG[type]);
968   
969   AliDielectronHistos *fHistos = die->GetHistoManager();
970   Bool_t pairClass1=fHistos->GetHistogramList()->FindObject(className1.Data())!=0x0;
971   Bool_t pairClass2=fHistos->GetHistogramList()->FindObject(className2.Data())!=0x0;
972
973   if (pairClass1 && PairTrackcut(dphiv, dcos, dmass, idie)==true){
974     ///import pair variables to values!!
975     values[AliDielectronVarManager::kPx] = dpxpair;
976     values[AliDielectronVarManager::kPy] = dpypair;
977     values[AliDielectronVarManager::kPz] = dpzpair;
978     values[AliDielectronVarManager::kPt] = dptpair;
979     values[AliDielectronVarManager::kXv] = dv0xvpair;
980     values[AliDielectronVarManager::kYv] = dv0yvpair;
981     values[AliDielectronVarManager::kZv] = dv0zvpair;
982     values[AliDielectronVarManager::kR] = dv0rpair;
983     values[AliDielectronVarManager::kE] = depair;
984     values[AliDielectronVarManager::kEta] = detapair;
985     values[AliDielectronVarManager::kM] = dmass;
986     values[AliDielectronVarManager::kPsiPair] = dphiv;
987     values[AliDielectronVarManager::kPhivPair] = dphiv;
988     values[AliDielectronVarManager::kPhi]  = dphipair;
989     values[AliDielectronVarManager::kOpeningAngle]  = dcos;
990     values[AliDielectronVarManager::kCosPointingAngle]  = TMath::Abs(TMath::ATan2(TMath::Sin(iep->Phi()-iem->Phi()),TMath::Cos(iep->Phi()-iem->Phi())));
991     fHistos->FillClass(className1, AliDielectronVarManager::kNMaxValues, values);
992   }
993
994
995   if (pairClass2 && PairTrackcut(dphiv, dopeningangle, dv0mass, idie)==true){
996     values[AliDielectronVarManager::kPx] = dv0pxpair;
997     values[AliDielectronVarManager::kPy] = dv0pypair;
998     values[AliDielectronVarManager::kPz] = dv0pzpair;
999     values[AliDielectronVarManager::kPt] = dv0ptpair;
1000     values[AliDielectronVarManager::kXv] = dv0xvpair;
1001     values[AliDielectronVarManager::kYv] = dv0yvpair;
1002     values[AliDielectronVarManager::kZv] = dv0zvpair;
1003     values[AliDielectronVarManager::kR] = dv0rpair;
1004     values[AliDielectronVarManager::kE] = dv0epair;
1005     values[AliDielectronVarManager::kEta] = dv0etapair;
1006     values[AliDielectronVarManager::kM] = dv0mass;
1007     values[AliDielectronVarManager::kPsiPair] = dpsipair;
1008     values[AliDielectronVarManager::kPhivPair] = dphivpair;
1009     values[AliDielectronVarManager::kPhi]  = dv0phipair;
1010     values[AliDielectronVarManager::kOpeningAngle]  = dopeningangle;
1011     values[AliDielectronVarManager::kCosPointingAngle]  = TMath::Abs(TMath::ATan2(TMath::Sin(iep->Phi()-iem->Phi()),TMath::Cos(iep->Phi()-iem->Phi())));
1012     fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
1013   }
1014
1015
1016   
1017 }
1018
1019 //_________________________________________________________________________________
1020 bool AliAnalysisTaskMultiDielectronTG::PairTrackcut(double phiv, double op, double mass, int idie)
1021 {
1022
1023   //
1024   // pair-by-pair cuts
1025   //
1026   
1027
1028   bool pairOK = true;
1029
1030   if(fRejectPairFlag[idie] == 1 || fRejectPairFlag[idie] == 2 ||
1031      fRejectPairFlag[idie] == 3 || fRejectPairFlag[idie] == 4 ){
1032     if(fRejectPairFlag[idie] == 2 || fRejectPairFlag[idie] == 4 ){
1033       if(fbz>0 && (phiv>fdconvphiv && mass < fdconvMee) ){
1034         pairOK = false;
1035       }else if(fbz<0 && phiv<acos(-1.0)-fdconvphiv  && mass < fdconvMee){
1036         pairOK = false;
1037       }
1038     }else if(fRejectPairFlag[idie] == 1 || fRejectPairFlag[idie] == 3){
1039       if(op<fdop){
1040         pairOK = false;
1041       }
1042     }
1043   }
1044
1045   return pairOK;
1046
1047 }
1048
1049
1050 //_________________________________________________________________________________
1051 void AliAnalysisTaskMultiDielectronTG::CalcVars(AliDielectronSingleTG *iep, AliDielectronSingleTG *iem, 
1052                                                 double &mass, double &phiv, double &px, double &py, double&pz,
1053                                                 double &pt, double &e, double &phi, 
1054                                                 double &eta, double &cos, double &psi)
1055 {
1056
1057
1058   //
1059   // standalone calculator for the pair variables
1060   //
1061
1062   px = iep->Px()+iem->Px();
1063   py = iep->Py()+iem->Py();
1064   pz = iep->Pz()+iem->Pz();
1065   pt = sqrt(px*px+py*py);
1066   double dppair = sqrt(pt*pt+pz*pz);
1067   static const double me=0.0005109989;
1068   e = sqrt(me*me+iep->Px()*iep->Px()+iep->Py()*iep->Py()+iep->Pz()*iep->Pz())
1069     + sqrt(me*me+iem->Px()*iem->Px()+iem->Py()*iem->Py()+iem->Pz()*iem->Pz());
1070   
1071   mass =  e*e-px*px-py*py-pz*pz;
1072   if(mass>=0){
1073     mass = sqrt(mass);
1074   }
1075    
1076   
1077   phi = atan2(py, px);
1078   eta = -0.5*TMath::Log((dppair+pz)/(dppair-pz));
1079   double p1 = sqrt(pow(iep->Px(),2)+pow(iep->Py(),2)+pow(iep->Pz(),2));
1080   double p2 = sqrt(pow(iem->Px(),2)+pow(iem->Py(),2)+pow(iem->Pz(),2));
1081   cos = acos((iep->Px()*iem->Px()+iep->Py()*iem->Py()+iep->Pz()*iem->Pz())/(p1*p2));
1082
1083   double dtheta = iep->Theta()-iem->Theta();
1084   psi = asin(dtheta/cos);
1085
1086
1087   //unit vector of (pep+pem) 
1088   double pl = dppair;
1089   double ux = px/pl;
1090   double uy = py/pl;
1091   double uz = pz/pl;
1092   double ax = uy/sqrt(ux*ux+uy*uy);
1093   double ay = -ux/sqrt(ux*ux+uy*uy); 
1094   
1095   //momentum of e+ and e- in (ax,ay,az) axis. Note that az=0 by 
1096   //definition. 
1097   //double ptep = iep->Px()*ax + iep->Py()*ay; 
1098   //double ptem = iem->Px()*ax + iem->Py()*ay; 
1099   
1100   double pxep = iep->Px();
1101   double pyep = iep->Py();
1102   double pzep = iep->Pz();
1103   double pxem = iem->Px();
1104   double pyem = iem->Py();
1105   double pzem = iem->Pz();
1106   
1107   
1108   //vector product of pep X pem 
1109   double vpx = pyep*pzem - pzep*pyem; 
1110   double vpy = pzep*pxem - pxep*pzem; 
1111   double vpz = pxep*pyem - pyep*pxem; 
1112   double vp = sqrt(vpx*vpx+vpy*vpy+vpz*vpz); 
1113   //double thev = acos(vpz/vp); 
1114   
1115   //unit vector of pep X pem 
1116   double vx = vpx/vp; 
1117   double vy = vpy/vp; 
1118   double vz = vpz/vp; 
1119
1120   //The third axis defined by vector product (ux,uy,uz)X(vx,vy,vz) 
1121   double wx = uy*vz - uz*vy; 
1122   double wy = uz*vx - ux*vz; 
1123   double wz = ux*vy - uy*vx; 
1124   double wl = sqrt(wx*wx+wy*wy+wz*wz); 
1125   // by construction, (wx,wy,wz) must be a unit vector. 
1126   if(fabs(wl - 1.0) > 0.00001) std::cout << "Calculation error in W vector"<<std::endl; 
1127   // measure angle between (wx,wy,wz) and (ax,ay,0). The angle between them 
1128   // should be small if the pair is conversion 
1129   //
1130   double cosPhiV = wx*ax + wy*ay; 
1131   phiv = acos(cosPhiV); 
1132   
1133 }
1134
1135 //_________________________________________________________________________________
1136 void AliAnalysisTaskMultiDielectronTG::ReshuffleBuffer(vector<AliDielectronSingleTG*> ve, deque<AliDielectronSingleTG*> pool)
1137 {
1138   //If there is not enough electron in the pool, give up
1139   //
1140   // ReshuffleBuffer for th event mixing 
1141   //
1142
1143   unsigned int ne = ve.size();
1144   unsigned int poolsize = pool.size();
1145   int used[fgkMAXPOOL];
1146   for(int i=0;i<(int)fgkMAXPOOL;i++){
1147     used[i]=0;
1148   }
1149
1150   if(poolsize < ne) {
1151     std::cout <<" pool size="<<poolsize<<" ne"<<ne<<std::endl;
1152     return;
1153   }
1154   for(unsigned int ie=0; ie < ne; ie++) {
1155     int j = rand()%poolsize;
1156     while(used[j]==1){
1157       j = rand()%poolsize;    
1158     }
1159     ve[ie] = pool[j];
1160     used[j]=1;
1161   }
1162
1163 }
1164
1165 //_________________________________________________________________________________
1166 void AliAnalysisTaskMultiDielectronTG::RejectPairs(vector<AliDielectronSingleTG*> e1, vector<AliDielectronSingleTG*> e2, Int_t idie){
1167
1168   ////////////////////////////////////
1169   ///// to reject pairs from track arrays 
1170   ///// conversions, ghost ..
1171   ///////////////////////////////////
1172   if(e1.size()>0 && e2.size()>0){
1173     for(int i1=0; i1<(int)e1.size(); i1++){
1174       for(int i2=0; i2<(int)e2.size(); i2++){
1175         if(fRejectPairFlag[idie]==1){
1176           Double_t cosine = GetOpeningAngle(e1[i1], e2[i2]);
1177           if(cosine<fdop){
1178             e1[i1]->SetConvFlag(0);
1179             e2[i2]->SetConvFlag(0);
1180           }
1181         }else if(fRejectPairFlag[idie]==2){
1182           Double_t phiv = GetPhiv(e1[i1], e2[i2]);
1183           Double_t mee = GetMass(e1[i1], e2[i2]);
1184           if(fbz>0 && ( phiv>fdconvphiv && mee < fdconvMee) ){
1185             e1[i1]->SetConvFlag(0);
1186             e2[i2]->SetConvFlag(0);
1187           }else if(fbz<0 && phiv<acos(-1.0)-fdconvphiv && mee < fdconvMee){
1188             e1[i1]->SetConvFlag(0);
1189             e2[i2]->SetConvFlag(0);
1190           }
1191         }
1192       }
1193     }
1194   }
1195   if(e1.size()>0){
1196     for(int i1=0; i1<(int)e1.size(); i1++){
1197       for(int i2=i1+1; i2<(int)e1.size(); i2++){
1198         if(fRejectPairFlag[idie]==1){
1199           Double_t cosine = GetOpeningAngle(e1[i1], e1[i2]);
1200           if(cosine<fdop){
1201             e1[i1]->SetConvFlag(0);
1202             e1[i2]->SetConvFlag(0);
1203           }
1204         }else if(fRejectPairFlag[idie]==2){
1205           Double_t phiv = GetPhiv(e1[i1], e1[i2]);
1206           Double_t mee = GetMass(e1[i1], e1[i2]);
1207           if(fbz>0 && phiv>fdconvphiv && mee < fdconvMee){
1208             e1[i1]->SetConvFlag(0);
1209             e1[i2]->SetConvFlag(0);
1210           }else if(fbz<0 && phiv<acos(-1.0)-fdconvphiv && mee < fdconvMee){
1211             e1[i1]->SetConvFlag(0);
1212             e1[i2]->SetConvFlag(0);
1213           }
1214         }
1215       }
1216     }
1217   }
1218
1219   if(e2.size()>0){
1220     for(int i1=0; i1<(int)e2.size(); i1++){
1221       for(int i2=i1+1; i2<(int)e2.size(); i2++){
1222         if(fRejectPairFlag[idie]==1){
1223           Double_t cosine = GetOpeningAngle(e2[i1], e2[i2]);
1224           if(cosine<fdop){
1225             e2[i1]->SetConvFlag(0);
1226             e2[i2]->SetConvFlag(0);
1227           }
1228         }else if(fRejectPairFlag[idie]==2){
1229           Double_t phiv = GetPhiv(e2[i1], e2[i2]);
1230           Double_t mee = GetMass(e2[i1], e2[i2]);
1231           if(fbz>0 && phiv>fdconvphiv && mee < fdconvMee){
1232             e2[i1]->SetConvFlag(0);
1233             e2[i2]->SetConvFlag(0);
1234           }else if(fbz<0 && phiv<acos(-1.0)-fdconvphiv && mee < fdconvMee){
1235             e2[i1]->SetConvFlag(0);
1236             e2[i2]->SetConvFlag(0);
1237           }
1238         }
1239       }
1240     }
1241   }
1242 }
1243
1244
1245 //_________________________________________________________________________________
1246 Double_t AliAnalysisTaskMultiDielectronTG::GetOpeningAngle(AliDielectronSingleTG* e1, AliDielectronSingleTG* e2){
1247
1248   //////////////////////
1249   //////// calculate pairs and get opening angle 
1250   //////////////////////
1251
1252   double dmass, dphiv, dpxpair, dpypair, dpzpair;
1253   double dptpair, depair, dphipair, detapair, dcos, dpsi;
1254
1255   CalcVars(e1, e2, dmass, dphiv, dpxpair, dpypair, dpzpair, 
1256             dptpair, depair, dphipair, detapair, dcos, dpsi);
1257
1258   return dcos;
1259 }
1260
1261 //_________________________________________________________________________________
1262 Double_t AliAnalysisTaskMultiDielectronTG::GetPhiv(AliDielectronSingleTG* e1, AliDielectronSingleTG* e2){
1263
1264   //////////////////////
1265   //////// calculate pairs and get phiv
1266   //////////////////////
1267
1268   double dmass, dphiv, dpxpair, dpypair, dpzpair;
1269   double dptpair, depair, dphipair, detapair, dcos, dpsi;
1270
1271   CalcVars(e1, e2, dmass, dphiv, dpxpair, dpypair, dpzpair, 
1272             dptpair, depair, dphipair, detapair, dcos, dpsi);
1273
1274   return dphiv;
1275 }
1276
1277 //_________________________________________________________________________________
1278 Double_t AliAnalysisTaskMultiDielectronTG::GetMass(AliDielectronSingleTG* e1, AliDielectronSingleTG* e2){
1279
1280   //////////////////////
1281   //////// calculate pairs and get mass
1282   //////////////////////
1283
1284   double dmass, dphiv, dpxpair, dpypair, dpzpair;
1285   double dptpair, depair, dphipair, detapair, dcos, dpsi;
1286   
1287   CalcVars(e1, e2, dmass, dphiv, dpxpair, dpypair, dpzpair, 
1288             dptpair, depair, dphipair, detapair, dcos, dpsi);
1289
1290   return dmass;
1291 }