Change Mult binning scheme
[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) mix->MixRemaining(die);
626   }
627   PostData(1, &fListHistos);
628   PostData(2, &fListCF);
629 }
630
631 //_________________________________________________________________________________
632 void AliAnalysisTaskMultiDielectronTG::CheckGhostPairs(vector<AliDielectronSingleTG*> e1)
633 {
634   ////// Check whether the pairs are in ghost pairs 
635   /////  ghost means that one of the pairs is highly associated but reconstructed as a true particle
636
637
638   bool reject = false;
639   if(e1.size()>1){
640     for(int i1=0; i1<(int)e1.size(); i1++){
641       reject = false;
642       for(int i2=i1+1; i2<(int)e1.size(); i2++){
643         if( fabs(e1[i1]->Phi() - e1[i2]->Phi())<0.01 && 
644             fabs(e1[i1]->Eta() - e1[i2]->Eta())<0.005
645             ){
646           reject = true;
647           e1[i2]->SetGstFlag(0);
648         }
649       }
650       if(reject==true)e1[i1]->SetGstFlag(0);
651     }
652   }
653 }
654
655 //_________________________________________________________________________________
656 Bool_t AliAnalysisTaskMultiDielectronTG::CheckGhost(vector<AliDielectronSingleTG*> e1, vector<AliDielectronSingleTG*> e2)
657 {
658   ////// To be sure whether there are no ghost pairs in h event mixing 
659
660   if(e1.size()>0 && e2.size()>0){
661     for(int i1=0; i1<(int)e1.size(); i1++){
662       for(int i2=0; i2<(int)e2.size(); i2++){
663         if( fabs(e1[i1]->Phi() - e2[i2]->Phi())<0.01 && 
664             fabs(e1[i1]->Eta() - e2[i2]->Eta())<0.005
665             ){
666           return true;
667         }
668       }
669     }
670   }
671   return false;
672 }
673
674 //_________________________________________________________________________________
675 void AliAnalysisTaskMultiDielectronTG::RandomizePool(vector<AliDielectronSingleTG*> e1, vector<AliDielectronSingleTG*> e2)
676 {
677   
678   ///// Randomize the pool constent to cancel the filling scheme of the single tracks
679   //// 
680
681
682   int size1 = e1.size();
683   int usedindex[1000];
684   for(int i=0;i<1000;i++){
685     usedindex[i] = -1;
686   }
687   for(int i=0;i<size1;i++){
688     usedindex[i] = 0;
689   }
690
691   for(int i=0;i<size1;i++){
692     int j = (int)(gRandom->Uniform(0,size1));
693     while(usedindex[j]==1){
694       j = (int)(gRandom->Uniform(0,size1));
695     }
696     if( (e1[j]->GetGstFlag()==1) &&
697         (e1[j]->GetConvFlag()==1)
698         ){
699       fVep.push_back(e1[j]);
700     }
701     usedindex[j] = 1;
702   }
703   
704
705   int size2 = e2.size();
706   for(int i=0;i<1000;i++){
707     usedindex[i] = -1;
708   }
709   for(int i=0;i<size2;i++){
710     usedindex[i] = 0;
711   }
712
713   for(int i=0;i<size2;i++){
714     int j = (int)(gRandom->Uniform(0,size2));
715     while(usedindex[j]==1){
716       j = (int)(gRandom->Uniform(0,size2));
717     }
718     if( (e2[j]->GetGstFlag()==1) &&
719         (e2[j]->GetConvFlag()==1)
720         ){
721       fVem.push_back(e2[j]);
722     }
723     usedindex[j] = 1;
724   }
725 }
726
727
728 //_________________________________________________________________________________
729 void AliAnalysisTaskMultiDielectronTG::CalcPair(vector<AliDielectronSingleTG*> ve1, 
730                                                 vector<AliDielectronSingleTG*> ve2, AliDielectron *die, Int_t idie)
731 {
732
733   //
734   // main routine for the pair procesing 
735   //
736
737
738   for(int i1=0; i1<(int)ve1.size(); i1++){
739     for(int i2=0; i2<(int)ve2.size(); i2++){
740       FillPair(ve1[i1], ve2[i2], 0, die, idie);  
741     }    
742   }
743
744   for(int i1=0; i1<(int)ve1.size(); i1++){
745     for(int i2=i1+1; i2<(int)ve1.size(); i2++){
746       FillPair(ve1[i1], ve1[i2], 1, die, idie);  
747     }    
748   }
749
750   for(int i1=0; i1<(int)ve2.size(); i1++){
751     for(int i2=i1+1; i2<(int)ve2.size(); i2++){
752       FillPair(ve2[i1], ve2[i2], 2, die, idie );  
753     }    
754   }
755
756
757   double dwcent = 100.0/fgkNCENT;
758   double dwiz = 20.0/fgkNZBIN;
759   double dwrp = acos(-1.0)/fgkNRPBIN;
760
761   int icent = (int)(fgValues[AliDielectronVarManager::kCentrality]/dwcent);
762   int izbin = (int)((fgValues[AliDielectronVarManager::kZvPrim]+10)/dwiz);
763   int irp = (int)((fgValues[AliDielectronVarManager::kV0ACrpH2])/dwrp);
764
765   if(icent<0) icent=0;
766   if(icent>=fgkNCENT) icent=fgkNCENT-1;
767   if(izbin<0) izbin=0;
768   if(izbin>=fgkNZBIN) izbin=fgkNZBIN-1;
769   if(irp<0) irp=0;
770   if(irp>=fgkNRPBIN) irp=fgkNRPBIN-1;
771
772   int nmixed;
773   if(ve1.size()>0) {
774     //
775     // Now mixed event for +- pairs
776     //
777     nmixed = 0;
778     for(int ibuf=0;(nmixed<fgkNMix);ibuf++) {
779       int ntry = 0;
780       while((fBGRejUnlike && CheckGhost(ve1, fvem[ibuf][idie][izbin][icent][irp])) &&  ntry<fgkMAXTRY) {
781         ReshuffleBuffer(fvem[ibuf][idie][izbin][icent][irp],fpoolm[idie][izbin][icent][irp]);
782         ntry++;
783       }
784       for(int i1=0; i1<(int)ve1.size(); i1++){
785         for(int i2=0; i2<(int)fvem[ibuf][idie][izbin][icent][irp].size(); i2++){
786           FillPair(ve1[i1],fvem[ibuf][idie][izbin][icent][irp][i2], 3, die, idie);
787         }
788       }
789       ++nmixed;
790     }//for(ibuf)
791   }
792   if(ve2.size()>0) {
793     //
794     // Now mixed event for -+ pairs
795     //
796     nmixed = 0;
797     for(int ibuf=0;(nmixed<fgkNMix);ibuf++) {
798       int ntry = 0;
799       while((fBGRejUnlike && CheckGhost(ve2, fvep[ibuf][idie][izbin][icent][irp])) &&  ntry<fgkMAXTRY) {
800         ReshuffleBuffer(fvep[ibuf][idie][izbin][icent][irp],fpoolp[idie][izbin][icent][irp]);
801         ntry++;
802       }
803       for(int i1=0; i1<(int)ve2.size(); i1++){
804         for(int i2=0; i2<(int)fvep[ibuf][idie][izbin][icent][irp].size(); i2++){
805           FillPair(fvep[ibuf][idie][izbin][icent][irp][i2],ve2[i1],4, die, idie);
806         }
807       }
808       ++nmixed;
809     }//for(ibuf)
810   }
811
812   if(ve1.size()>0) {
813     //
814     // Now mixed event for ++ pairs
815     //
816     nmixed = 0;
817     for(int ibuf=0;(nmixed<fgkNMix);ibuf++) {
818       int ntry = 0;
819       while((fBGRejLike && CheckGhost(ve1, fvep[ibuf][idie][izbin][icent][irp])) &&  ntry<fgkMAXTRY) {
820         ReshuffleBuffer(fvep[ibuf][idie][izbin][icent][irp],fpoolp[idie][izbin][icent][irp]);
821         ntry++;
822       }
823       for(int i1=0; i1<(int)ve1.size(); i1++){
824         for(int i2=0;i2<(int)fvep[ibuf][idie][izbin][icent][irp].size(); i2++){
825           FillPair(ve1[i1],fvep[ibuf][idie][izbin][icent][irp][i2], 5, die, idie);
826         }
827       }
828       ++nmixed;
829     }//for(ibuf)
830   }
831
832   if(ve2.size()>0) {
833     //
834     // Now mixed event for +- pairs
835     //
836     nmixed = 0;
837     for(int ibuf=0;(nmixed<fgkNMix);ibuf++) {
838       int ntry = 0;
839       while((fBGRejLike && CheckGhost(ve2, fvem[ibuf][idie][izbin][icent][irp])) &&  ntry<fgkMAXTRY) {
840         ReshuffleBuffer(fvem[ibuf][idie][izbin][icent][irp],fpoolm[idie][izbin][icent][irp]);
841         ntry++;
842       }
843       for(int i1=0; i1<(int)ve2.size(); i1++){
844         for(int i2=0; i2<(int)fvem[ibuf][idie][izbin][icent][irp].size(); i2++){
845           FillPair(ve2[i1],fvem[ibuf][idie][izbin][icent][irp][i2],6, die, idie);
846         }
847       }
848       ++nmixed;
849     }//for(ibuf)
850   }
851
852 }
853
854
855 //_________________________________________________________________________________
856 void AliAnalysisTaskMultiDielectronTG::FillPair(AliDielectronSingleTG *iep, 
857                                                 AliDielectronSingleTG *iem, int type, AliDielectron *die, Int_t idie)
858 {                      
859
860   //
861   // main routine for filling kinematics of pairs
862   //
863
864
865   
866   double dmass, dphiv, dpxpair, dpypair, dpzpair;
867   double dptpair, depair, dphipair, detapair, dcos, dpsi;
868
869   CalcVars(iep, iem, dmass, dphiv, dpxpair, dpypair, dpzpair, 
870             dptpair, depair, dphipair, detapair, dcos, dpsi);
871
872
873   double dopeningangle =  -9999;
874   double dv0mass =  -9999;
875   double dv0pxpair = -9999;
876   double dv0pypair = -9999;
877   double dv0pzpair = -9999;
878   double dv0ptpair = -9999;
879   double dv0epair = -9999;
880   double dv0xvpair = -9999;
881   double dv0yvpair = -9999;
882   double dv0zvpair = -9999;
883   double dv0phipair = -9999;
884   double dv0etapair = -9999;
885   double dv0rpair = -9999;
886   double dpsipair =  -9999;
887   double dphivpair =  -9999;
888
889   ////////////////////////////
890   ///// calculate v0 ////////
891   ///////////////////////////
892   Bool_t v0OFF=kFALSE;
893   /// for the moment, this doesn't work for MC
894   if(fdv0mixing == kFALSE && (type==3 || type==4 || type==5 || type==6)){
895     v0OFF = kTRUE;
896   }
897   if(die->GetHasMC()==kTRUE && (type==3 || type==4 || type==5 || type==6)){
898     v0OFF = kTRUE;
899   }
900   if(type==0 || type==1 || type==2){
901     v0OFF = kFALSE;
902   }
903   
904
905   if(v0OFF==kFALSE){
906     AliVTrack* trackob1= iep->GetTrack();    
907     AliVTrack* trackob2= iem->GetTrack();    
908
909     if(!trackob1 || !trackob2){
910       return; 
911     }
912
913     AliDielectronPair candidate;
914     candidate.SetTracks(trackob1, (int)(11*iep->Charge()), 
915                         trackob2, (int)(11*iem->Charge()));
916     
917     if(trackob1==trackob2){
918       AliInfo("Same Objects!!");
919       return; 
920     }
921
922     //const AliKFParticle &kfPairLeg1 = candidate.GetKFFirstDaughter();
923     //const AliKFParticle &kfPairLeg2 = candidate.GetKFSecondDaughter();
924     
925     const AliKFParticle &kfPair = candidate.GetKFParticle();
926
927     /*
928     dv0mass = candidate.M();
929     dv0pxpair = candidate.Px();
930     dv0pypair = candidate.Py();
931     dv0pzpair = candidate.Pz();
932     dv0ptpair = candidate.Pt();
933     dv0epair = candidate.E();
934     dv0xvpair = candidate.Xv();
935     dv0yvpair = candidate.Yv();
936     dv0zvpair = candidate.Zv();
937     dv0phipair = candidate.Phi();
938     dv0etapair = candidate.Eta();
939     */
940     dv0mass = kfPair.GetMass();
941     dv0pxpair = kfPair.GetPx();
942     dv0pypair = kfPair.GetPy();
943     dv0pzpair = kfPair.GetPz();
944     dv0ptpair = kfPair.GetPt();
945     dv0epair = kfPair.GetE();
946     dv0xvpair = kfPair.GetX();
947     dv0yvpair = kfPair.GetY();
948     dv0zvpair = kfPair.GetZ();
949     dv0phipair = kfPair.GetPhi();
950     dv0etapair = kfPair.GetEta();
951     dv0rpair =  kfPair.GetR();
952
953     dopeningangle = candidate.OpeningAngle();    
954     dpsipair = candidate.PsiPair(fbz);
955     dphivpair = candidate.PhivPair(fbz);
956
957
958   }
959
960   Double_t values[AliDielectronVarManager::kNMaxValues];
961   TString  className1;
962   TString  className2;
963   className1.Form("MyPair_%s",kPairClassNamesTG[type]);
964   className2.Form("MyPairV0_%s",kPairClassNamesTG[type]);
965   
966   AliDielectronHistos *fHistos = die->GetHistoManager();
967   Bool_t pairClass1=fHistos->GetHistogramList()->FindObject(className1.Data())!=0x0;
968   Bool_t pairClass2=fHistos->GetHistogramList()->FindObject(className2.Data())!=0x0;
969
970   if (pairClass1 && PairTrackcut(dphiv, dcos, dmass, idie)==true){
971     ///import pair variables to values!!
972     values[AliDielectronVarManager::kPx] = dpxpair;
973     values[AliDielectronVarManager::kPy] = dpypair;
974     values[AliDielectronVarManager::kPz] = dpzpair;
975     values[AliDielectronVarManager::kPt] = dptpair;
976     values[AliDielectronVarManager::kXv] = dv0xvpair;
977     values[AliDielectronVarManager::kYv] = dv0yvpair;
978     values[AliDielectronVarManager::kZv] = dv0zvpair;
979     values[AliDielectronVarManager::kR] = dv0rpair;
980     values[AliDielectronVarManager::kE] = depair;
981     values[AliDielectronVarManager::kEta] = detapair;
982     values[AliDielectronVarManager::kM] = dmass;
983     values[AliDielectronVarManager::kPsiPair] = dphiv;
984     values[AliDielectronVarManager::kPhivPair] = dphiv;
985     values[AliDielectronVarManager::kPhi]  = dphipair;
986     values[AliDielectronVarManager::kOpeningAngle]  = dcos;
987     values[AliDielectronVarManager::kCosPointingAngle]  = TMath::Abs(TMath::ATan2(TMath::Sin(iep->Phi()-iem->Phi()),TMath::Cos(iep->Phi()-iem->Phi())));
988     fHistos->FillClass(className1, AliDielectronVarManager::kNMaxValues, values);
989   }
990
991
992   if (pairClass2 && PairTrackcut(dphiv, dopeningangle, dv0mass, idie)==true){
993     values[AliDielectronVarManager::kPx] = dv0pxpair;
994     values[AliDielectronVarManager::kPy] = dv0pypair;
995     values[AliDielectronVarManager::kPz] = dv0pzpair;
996     values[AliDielectronVarManager::kPt] = dv0ptpair;
997     values[AliDielectronVarManager::kXv] = dv0xvpair;
998     values[AliDielectronVarManager::kYv] = dv0yvpair;
999     values[AliDielectronVarManager::kZv] = dv0zvpair;
1000     values[AliDielectronVarManager::kR] = dv0rpair;
1001     values[AliDielectronVarManager::kE] = dv0epair;
1002     values[AliDielectronVarManager::kEta] = dv0etapair;
1003     values[AliDielectronVarManager::kM] = dv0mass;
1004     values[AliDielectronVarManager::kPsiPair] = dpsipair;
1005     values[AliDielectronVarManager::kPhivPair] = dphivpair;
1006     values[AliDielectronVarManager::kPhi]  = dv0phipair;
1007     values[AliDielectronVarManager::kOpeningAngle]  = dopeningangle;
1008     values[AliDielectronVarManager::kCosPointingAngle]  = TMath::Abs(TMath::ATan2(TMath::Sin(iep->Phi()-iem->Phi()),TMath::Cos(iep->Phi()-iem->Phi())));
1009     fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
1010   }
1011
1012
1013   
1014 }
1015
1016 //_________________________________________________________________________________
1017 bool AliAnalysisTaskMultiDielectronTG::PairTrackcut(double phiv, double op, double mass, int idie)
1018 {
1019
1020   //
1021   // pair-by-pair cuts
1022   //
1023   
1024
1025   bool pairOK = true;
1026
1027   if(fRejectPairFlag[idie] == 1 || fRejectPairFlag[idie] == 2 ||
1028      fRejectPairFlag[idie] == 3 || fRejectPairFlag[idie] == 4 ){
1029     if(fRejectPairFlag[idie] == 2 || fRejectPairFlag[idie] == 4 ){
1030       if(fbz>0 && (phiv>fdconvphiv && mass < fdconvMee) ){
1031         pairOK = false;
1032       }else if(fbz<0 && phiv<acos(-1.0)-fdconvphiv  && mass < fdconvMee){
1033         pairOK = false;
1034       }
1035     }else if(fRejectPairFlag[idie] == 1 || fRejectPairFlag[idie] == 3){
1036       if(op<fdop){
1037         pairOK = false;
1038       }
1039     }
1040   }
1041
1042   return pairOK;
1043
1044 }
1045
1046
1047 //_________________________________________________________________________________
1048 void AliAnalysisTaskMultiDielectronTG::CalcVars(AliDielectronSingleTG *iep, AliDielectronSingleTG *iem, 
1049                                                 double &mass, double &phiv, double &px, double &py, double&pz,
1050                                                 double &pt, double &e, double &phi, 
1051                                                 double &eta, double &cos, double &psi)
1052 {
1053
1054
1055   //
1056   // standalone calculator for the pair variables
1057   //
1058
1059   px = iep->Px()+iem->Px();
1060   py = iep->Py()+iem->Py();
1061   pz = iep->Pz()+iem->Pz();
1062   pt = sqrt(px*px+py*py);
1063   double dppair = sqrt(pt*pt+pz*pz);
1064   static const double me=0.0005109989;
1065   e = sqrt(me*me+iep->Px()*iep->Px()+iep->Py()*iep->Py()+iep->Pz()*iep->Pz())
1066     + sqrt(me*me+iem->Px()*iem->Px()+iem->Py()*iem->Py()+iem->Pz()*iem->Pz());
1067   
1068   mass =  e*e-px*px-py*py-pz*pz;
1069   if(mass>=0){
1070     mass = sqrt(mass);
1071   }
1072    
1073   
1074   phi = atan2(py, px);
1075   eta = -0.5*TMath::Log((dppair+pz)/(dppair-pz));
1076   double p1 = sqrt(pow(iep->Px(),2)+pow(iep->Py(),2)+pow(iep->Pz(),2));
1077   double p2 = sqrt(pow(iem->Px(),2)+pow(iem->Py(),2)+pow(iem->Pz(),2));
1078   cos = acos((iep->Px()*iem->Px()+iep->Py()*iem->Py()+iep->Pz()*iem->Pz())/(p1*p2));
1079
1080   double dtheta = iep->Theta()-iem->Theta();
1081   psi = asin(dtheta/cos);
1082
1083
1084   //unit vector of (pep+pem) 
1085   double pl = dppair;
1086   double ux = px/pl;
1087   double uy = py/pl;
1088   double uz = pz/pl;
1089   double ax = uy/sqrt(ux*ux+uy*uy);
1090   double ay = -ux/sqrt(ux*ux+uy*uy); 
1091   
1092   //momentum of e+ and e- in (ax,ay,az) axis. Note that az=0 by 
1093   //definition. 
1094   //double ptep = iep->Px()*ax + iep->Py()*ay; 
1095   //double ptem = iem->Px()*ax + iem->Py()*ay; 
1096   
1097   double pxep = iep->Px();
1098   double pyep = iep->Py();
1099   double pzep = iep->Pz();
1100   double pxem = iem->Px();
1101   double pyem = iem->Py();
1102   double pzem = iem->Pz();
1103   
1104   
1105   //vector product of pep X pem 
1106   double vpx = pyep*pzem - pzep*pyem; 
1107   double vpy = pzep*pxem - pxep*pzem; 
1108   double vpz = pxep*pyem - pyep*pxem; 
1109   double vp = sqrt(vpx*vpx+vpy*vpy+vpz*vpz); 
1110   //double thev = acos(vpz/vp); 
1111   
1112   //unit vector of pep X pem 
1113   double vx = vpx/vp; 
1114   double vy = vpy/vp; 
1115   double vz = vpz/vp; 
1116
1117   //The third axis defined by vector product (ux,uy,uz)X(vx,vy,vz) 
1118   double wx = uy*vz - uz*vy; 
1119   double wy = uz*vx - ux*vz; 
1120   double wz = ux*vy - uy*vx; 
1121   double wl = sqrt(wx*wx+wy*wy+wz*wz); 
1122   // by construction, (wx,wy,wz) must be a unit vector. 
1123   if(fabs(wl - 1.0) > 0.00001) std::cout << "Calculation error in W vector"<<std::endl; 
1124   // measure angle between (wx,wy,wz) and (ax,ay,0). The angle between them 
1125   // should be small if the pair is conversion 
1126   //
1127   double cosPhiV = wx*ax + wy*ay; 
1128   phiv = acos(cosPhiV); 
1129   
1130 }
1131
1132 //_________________________________________________________________________________
1133 void AliAnalysisTaskMultiDielectronTG::ReshuffleBuffer(vector<AliDielectronSingleTG*> ve, deque<AliDielectronSingleTG*> pool)
1134 {
1135   //If there is not enough electron in the pool, give up
1136   //
1137   // ReshuffleBuffer for th event mixing 
1138   //
1139
1140   unsigned int ne = ve.size();
1141   unsigned int poolsize = pool.size();
1142   int used[fgkMAXPOOL];
1143   for(int i=0;i<(int)fgkMAXPOOL;i++){
1144     used[i]=0;
1145   }
1146
1147   if(poolsize < ne) {
1148     std::cout <<" pool size="<<poolsize<<" ne"<<ne<<std::endl;
1149     return;
1150   }
1151   for(unsigned int ie=0; ie < ne; ie++) {
1152     int j = rand()%poolsize;
1153     while(used[j]==1){
1154       j = rand()%poolsize;    
1155     }
1156     ve[ie] = pool[j];
1157     used[j]=1;
1158   }
1159
1160 }
1161
1162 //_________________________________________________________________________________
1163 void AliAnalysisTaskMultiDielectronTG::RejectPairs(vector<AliDielectronSingleTG*> e1, vector<AliDielectronSingleTG*> e2, Int_t idie){
1164
1165   ////////////////////////////////////
1166   ///// to reject pairs from track arrays 
1167   ///// conversions, ghost ..
1168   ///////////////////////////////////
1169   if(e1.size()>0 && e2.size()>0){
1170     for(int i1=0; i1<(int)e1.size(); i1++){
1171       for(int i2=0; i2<(int)e2.size(); i2++){
1172         if(fRejectPairFlag[idie]==1){
1173           Double_t cosine = GetOpeningAngle(e1[i1], e2[i2]);
1174           if(cosine<fdop){
1175             e1[i1]->SetConvFlag(0);
1176             e2[i2]->SetConvFlag(0);
1177           }
1178         }else if(fRejectPairFlag[idie]==2){
1179           Double_t phiv = GetPhiv(e1[i1], e2[i2]);
1180           Double_t mee = GetMass(e1[i1], e2[i2]);
1181           if(fbz>0 && ( phiv>fdconvphiv && mee < fdconvMee) ){
1182             e1[i1]->SetConvFlag(0);
1183             e2[i2]->SetConvFlag(0);
1184           }else if(fbz<0 && phiv<acos(-1.0)-fdconvphiv && mee < fdconvMee){
1185             e1[i1]->SetConvFlag(0);
1186             e2[i2]->SetConvFlag(0);
1187           }
1188         }
1189       }
1190     }
1191   }
1192   if(e1.size()>0){
1193     for(int i1=0; i1<(int)e1.size(); i1++){
1194       for(int i2=i1+1; i2<(int)e1.size(); i2++){
1195         if(fRejectPairFlag[idie]==1){
1196           Double_t cosine = GetOpeningAngle(e1[i1], e1[i2]);
1197           if(cosine<fdop){
1198             e1[i1]->SetConvFlag(0);
1199             e1[i2]->SetConvFlag(0);
1200           }
1201         }else if(fRejectPairFlag[idie]==2){
1202           Double_t phiv = GetPhiv(e1[i1], e1[i2]);
1203           Double_t mee = GetMass(e1[i1], e1[i2]);
1204           if(fbz>0 && phiv>fdconvphiv && mee < fdconvMee){
1205             e1[i1]->SetConvFlag(0);
1206             e1[i2]->SetConvFlag(0);
1207           }else if(fbz<0 && phiv<acos(-1.0)-fdconvphiv && mee < fdconvMee){
1208             e1[i1]->SetConvFlag(0);
1209             e1[i2]->SetConvFlag(0);
1210           }
1211         }
1212       }
1213     }
1214   }
1215
1216   if(e2.size()>0){
1217     for(int i1=0; i1<(int)e2.size(); i1++){
1218       for(int i2=i1+1; i2<(int)e2.size(); i2++){
1219         if(fRejectPairFlag[idie]==1){
1220           Double_t cosine = GetOpeningAngle(e2[i1], e2[i2]);
1221           if(cosine<fdop){
1222             e2[i1]->SetConvFlag(0);
1223             e2[i2]->SetConvFlag(0);
1224           }
1225         }else if(fRejectPairFlag[idie]==2){
1226           Double_t phiv = GetPhiv(e2[i1], e2[i2]);
1227           Double_t mee = GetMass(e2[i1], e2[i2]);
1228           if(fbz>0 && phiv>fdconvphiv && mee < fdconvMee){
1229             e2[i1]->SetConvFlag(0);
1230             e2[i2]->SetConvFlag(0);
1231           }else if(fbz<0 && phiv<acos(-1.0)-fdconvphiv && mee < fdconvMee){
1232             e2[i1]->SetConvFlag(0);
1233             e2[i2]->SetConvFlag(0);
1234           }
1235         }
1236       }
1237     }
1238   }
1239 }
1240
1241
1242 //_________________________________________________________________________________
1243 Double_t AliAnalysisTaskMultiDielectronTG::GetOpeningAngle(AliDielectronSingleTG* e1, AliDielectronSingleTG* e2){
1244
1245   //////////////////////
1246   //////// calculate pairs and get opening angle 
1247   //////////////////////
1248
1249   double dmass, dphiv, dpxpair, dpypair, dpzpair;
1250   double dptpair, depair, dphipair, detapair, dcos, dpsi;
1251
1252   CalcVars(e1, e2, dmass, dphiv, dpxpair, dpypair, dpzpair, 
1253             dptpair, depair, dphipair, detapair, dcos, dpsi);
1254
1255   return dcos;
1256 }
1257
1258 //_________________________________________________________________________________
1259 Double_t AliAnalysisTaskMultiDielectronTG::GetPhiv(AliDielectronSingleTG* e1, AliDielectronSingleTG* e2){
1260
1261   //////////////////////
1262   //////// calculate pairs and get phiv
1263   //////////////////////
1264
1265   double dmass, dphiv, dpxpair, dpypair, dpzpair;
1266   double dptpair, depair, dphipair, detapair, dcos, dpsi;
1267
1268   CalcVars(e1, e2, dmass, dphiv, dpxpair, dpypair, dpzpair, 
1269             dptpair, depair, dphipair, detapair, dcos, dpsi);
1270
1271   return dphiv;
1272 }
1273
1274 //_________________________________________________________________________________
1275 Double_t AliAnalysisTaskMultiDielectronTG::GetMass(AliDielectronSingleTG* e1, AliDielectronSingleTG* e2){
1276
1277   //////////////////////
1278   //////// calculate pairs and get mass
1279   //////////////////////
1280
1281   double dmass, dphiv, dpxpair, dpypair, dpzpair;
1282   double dptpair, depair, dphipair, detapair, dcos, dpsi;
1283   
1284   CalcVars(e1, e2, dmass, dphiv, dpxpair, dpypair, dpzpair, 
1285             dptpair, depair, dphipair, detapair, dcos, dpsi);
1286
1287   return dmass;
1288 }