Place the config and root file at the right place
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliHFEreducedEventCreatorAOD.cxx
1 /**************************************************************************
2 * Copyright(c) 1998-1999, 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 // Debug tree task
17 // the tree is represented as reduced events
18 // 
19 // Authors:
20 //   M.Fasel <M.Fasel@gsi.de>
21 //
22 #include <TArrayI.h>
23 #include <TBits.h>
24 #include <TFile.h>
25 #include <TString.h>
26 #include <TTree.h>
27 #include <TMath.h>
28
29 #include "AliAnalysisManager.h"
30 #include "AliAnalysisUtils.h"
31 #include "AliAODEvent.h"
32 #include "AliAODMCHeader.h"
33 #include "AliAODMCParticle.h"
34 #include "AliAODPid.h"
35 #include "AliAODTrack.h"
36 #include "AliAODTracklets.h"
37 #include "AliAODVertex.h"
38 #include "AliCentrality.h"
39 #include "AliInputEventHandler.h"
40 #include "AliHFEcuts.h"
41 #include "AliHFEextraCuts.h"
42 #include "AliHFEpidTPC.h"
43 #include "AliHFEreducedEvent.h"
44 #include "AliHFEreducedTrack.h"
45 #include "AliHFEreducedMCParticle.h"
46 #include "AliHFEsignalCuts.h"
47 #include "AliLog.h"
48 #include "AliMCEvent.h"
49 #include "AliPIDResponse.h"
50 #include "AliVEvent.h"
51 #include "AliVTrack.h"
52 #include "AliVCluster.h"
53 #include "AliVVZERO.h"
54 #include "AliVZDC.h"
55 #include "TTreeStream.h"
56
57 #include "AliHFEreducedEventCreatorAOD.h"
58
59 ClassImp(AliHFEreducedEventCreatorAOD)
60
61 AliHFEreducedEventCreatorAOD::AliHFEreducedEventCreatorAOD():
62   AliAnalysisTaskSE(),
63   fAODMCHeader(NULL),
64   fAODArrayMCInfo(NULL),
65   fHFEtree(NULL),
66   fAnalysisUtils(NULL),
67   fHFEevent(NULL),
68   fTrackCuts(NULL),
69   fExtraCuts(NULL),
70   fSignalCuts(NULL),
71   fTPCpid(NULL),
72   fEventNumber(0),
73   fNclustersTPC(70),
74   fNclustersTPCPID(0),
75   fNclustersITS(2),
76   fNbOfTOFSigma(-1.0),
77   fRemoveFirstEvent(kFALSE)
78 {
79   //
80   // Default constructor
81   //
82 }
83
84 AliHFEreducedEventCreatorAOD::AliHFEreducedEventCreatorAOD(const char *name):
85   AliAnalysisTaskSE(name),
86   fAODMCHeader(NULL),
87   fAODArrayMCInfo(NULL),
88   fHFEtree(NULL),
89   fAnalysisUtils(NULL),
90   fHFEevent(NULL),
91   fTrackCuts(NULL),
92   fExtraCuts(NULL),
93   fSignalCuts(NULL),
94   fTPCpid(NULL),
95   fEventNumber(0),
96   fNclustersTPC(70),
97   fNclustersTPCPID(0),
98   fNclustersITS(2),
99   fNbOfTOFSigma(-1.0),
100   fRemoveFirstEvent(kFALSE)
101 {
102   //
103   // Default constructor
104   //
105   fTPCpid = new AliHFEpidTPC("QAtpcPID");
106   fAnalysisUtils = new AliAnalysisUtils;
107   DefineOutput(1, TTree::Class());
108 }
109
110 AliHFEreducedEventCreatorAOD::~AliHFEreducedEventCreatorAOD(){
111   //
112   // Default destructor
113   //
114     if(fTPCpid) delete fTPCpid;
115     if(fAnalysisUtils) delete fAnalysisUtils;
116     if(fHFEevent) delete fHFEevent;
117 }
118
119 void AliHFEreducedEventCreatorAOD::UserCreateOutputObjects(){
120   //
121   // Create debug tree, signal cuts and track cuts
122   //
123
124   //printf("test\n");
125
126  // printf("testa\n");
127   fSignalCuts = new AliHFEsignalCuts("HFEsignalCuts", "HFE MC Signal definition");
128   //printf("testb\n");
129   
130   fTrackCuts = new AliHFEcuts("fTrackCuts", "Basic HFE track cuts");
131   fTrackCuts->CreateStandardCuts();
132   fTrackCuts->SetAOD();
133   // Track cuts
134   fTrackCuts->SetMinNClustersTPC(fNclustersTPC);
135   fTrackCuts->SetMinRatioTPCclusters(0);
136   fTrackCuts->SetTPCmodes(AliHFEextraCuts::kFound, AliHFEextraCuts::kFoundOverFindable); 
137   fTrackCuts->SetMinNClustersTPCPID(fNclustersTPCPID);
138   fTrackCuts->SetMinNClustersITS(fNclustersITS);
139   // Event cuts
140   fTrackCuts->SetUseMixedVertex(kTRUE);
141   fTrackCuts->SetVertexRange(10.);
142   //printf("testa\n");
143   fTrackCuts->Initialize();
144   //printf("testb\n");
145
146   fExtraCuts = new AliHFEextraCuts("hfeExtraCuts","HFE Extra Cuts");
147
148   fHFEevent = new AliHFEreducedEvent;
149   OpenFile(1);
150   fHFEtree = new TTree("HFEtree", "HFE event tree");
151   fHFEtree->Branch("HFEevent", "AliHFEreducedEvent", fHFEevent,128000,0);
152   PostData(1, fHFEtree);
153 }
154
155 void AliHFEreducedEventCreatorAOD::UserExec(Option_t *){
156   //
157   // User Exec: Fill debug Tree
158   // 
159
160   // Get PID response
161   AliPIDResponse *pid = NULL;
162   AliInputEventHandler *handler = dynamic_cast<AliInputEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
163   if(handler){
164     //    printf("testb\n");
165     pid = handler->GetPIDResponse();
166   } else {
167     AliError("No Handler");
168   }
169   if(!pid){
170     //   printf("testc\n");
171     AliError("No PID response");
172     return;
173   }
174   if(!fInputEvent) {
175     AliError("No Input event");
176     return;
177   }
178
179   if(fRemoveFirstEvent){
180     if(fAnalysisUtils->IsFirstEventInChunk(fInputEvent)) return;
181   }
182
183   AliAODTrack copyTrack;
184
185   // MC info
186   Bool_t mcthere = kTRUE;
187   AliAODEvent *aodE = dynamic_cast<AliAODEvent *>(fInputEvent);
188   if(!aodE){
189     //  printf("testd\n");
190     AliError("No AOD Event");
191     return;
192   }
193   fAODMCHeader = dynamic_cast<AliAODMCHeader *>(fInputEvent->FindListObject(AliAODMCHeader::StdBranchName()));
194   if(!fAODMCHeader){ 
195       mcthere = kFALSE;
196       //   return;
197   }
198   fAODArrayMCInfo = dynamic_cast<TClonesArray *>(fInputEvent->FindListObject(AliAODMCParticle::StdBranchName()));
199   if(!fAODArrayMCInfo){ 
200       mcthere = kFALSE;
201       //  return;
202   }
203   else {
204     fSignalCuts->SetMCAODInfo(fAODArrayMCInfo);
205     fTrackCuts->SetMCEvent(aodE);
206   }
207
208   // Set for track cuts
209   if(mcthere){
210     // Check Monte-Carlo events for AODs
211     AliDebug(1, Form("Monte-Carlo Event: %p\n", fMCEvent));
212     if(fMCEvent){
213       AliDebug(1, Form("Available! Number of particles: %d\n", fMCEvent->GetNumberOfTracks()));
214     }
215   }
216   fTrackCuts->SetRecEvent(fInputEvent);
217
218   if(!fTrackCuts->CheckEventCuts("fCutsEvRec", fInputEvent)){
219     AliDebug(1, "Event rejected by the event cuts\n");
220     return;
221   }
222   if(!fExtraCuts){
223     fExtraCuts = new AliHFEextraCuts("hfeExtraCuts","HFE Extra Cuts");
224   }
225   fExtraCuts->SetRecEventInfo(fInputEvent);
226
227  
228   // Make Reduced Event 
229   //AliHFEreducedEvent hfeevent;
230   fHFEevent->~AliHFEreducedEvent();
231   new(fHFEevent)AliHFEreducedEvent();
232
233   // Get run number
234   fHFEevent->SetRunNumber(fInputEvent->GetRunNumber());
235
236   // Derive trigger 
237   AliDebug(1, "Get triggers\n");
238   UInt_t trigger = fInputHandler->IsEventSelected();
239   if(trigger & AliVEvent::kMB) fHFEevent->SetMBTrigger();
240   if(trigger & AliVEvent::kCentral) fHFEevent->SetCentralTrigger();
241   if(trigger & AliVEvent::kSemiCentral) fHFEevent->SetCentralTrigger();
242   if(trigger & AliVEvent::kEMCEJE) fHFEevent->SetEMCALTrigger();
243
244   // Get Primary Vertex
245   AliDebug(1, "Get Primary Vertex\n");
246   const AliVVertex *vertex = fInputEvent->GetPrimaryVertex();
247   Double_t vtx[3];
248   Double_t vcov[6];
249   Int_t ncontrib = -1;
250   if(vertex) {
251     AliDebug(1, "Found vertex\n");
252     vertex->GetXYZ(vtx);
253     ncontrib = vertex->GetNContributors();
254     vertex->GetCovarianceMatrix(vcov);
255   }
256   fHFEevent->SetVX(vtx[0]);
257   fHFEevent->SetVY(vtx[1]);
258   fHFEevent->SetVZ(vtx[2]);
259   fHFEevent->SetNContribVertex(ncontrib);
260   fHFEevent->SetVertexResolution(TMath::Sqrt(TMath::Abs(vcov[5])));
261   // Get Primary Vertex from SPD
262   const AliVVertex *vertexSPD = aodE->GetPrimaryVertexSPD();
263   if(vertexSPD){
264     AliDebug(1, "Found SPD vertex\n");
265     memset(vtx, 0, sizeof(Double_t) *3);
266     vertexSPD->GetXYZ(vtx);
267     fHFEevent->SetVXSPD(vtx[0]);
268     fHFEevent->SetVYSPD(vtx[1]);
269     fHFEevent->SetVZSPD(vtx[2]);
270     fHFEevent->SetNContribVertexSPD(vertexSPD->GetNContributors());
271     memset(vcov, 0, sizeof(Double_t)*6);
272     vertexSPD->GetCovarianceMatrix(vcov);
273     AliDebug(1, Form("Covariance Matrix vcov[5] %f\n",vcov[5]));
274     fHFEevent->SetVertexResolutionSPD(TMath::Sqrt(TMath::Abs(vcov[5])));
275   }
276
277   // Get centrality
278   AliDebug(1, "Centrality\n");
279   AliCentrality *hicent = fInputEvent->GetCentrality();
280   if(hicent) fHFEevent->SetCentrality(
281     hicent->GetCentralityPercentile("V0M"),
282     hicent->GetCentralityPercentile("V0A"),
283     hicent->GetCentralityPercentile("V0C"),
284     hicent->GetCentralityPercentile("TKL"),
285     hicent->GetCentralityPercentile("TRK"),
286     hicent->GetCentralityPercentile("ZNA"),
287     hicent->GetCentralityPercentile("ZNC"),
288     hicent->GetCentralityPercentile("CL0"),
289     hicent->GetCentralityPercentile("CL1"),
290     hicent->GetCentralityPercentile("CND")
291   );
292
293   // Get VZERO Information
294   AliDebug(1, "VZERO info\n");
295   AliVVZERO *vzeroinfo = fInputEvent->GetVZEROData();
296   if(vzeroinfo) fHFEevent->SetV0Multiplicity(vzeroinfo->GetMTotV0A(), vzeroinfo->GetMTotV0C());
297
298   // Get ZDC Information
299   AliDebug(1, "ZDC info\n");
300   AliVZDC *zdcinfo = fInputEvent->GetZDCData();
301   if(zdcinfo) fHFEevent->SetZDCEnergy(zdcinfo->GetZNAEnergy(), zdcinfo->GetZNCEnergy(), zdcinfo->GetZPAEnergy(), zdcinfo->GetZPCEnergy()); 
302   
303   // Set SPD multiplicity
304   AliDebug(1, "SPD multiplicity\n");
305   AliAODTracklets *tls = aodE->GetTracklets();
306   if(tls) fHFEevent->SetSPDMultiplicity(tls->GetNumberOfTracklets());
307
308   // Look for kink mother
309   AliDebug(1, "Vertices\n");
310   Int_t numberofvertices = aodE->GetNumberOfVertices();
311   Double_t listofmotherkink[numberofvertices];
312   Int_t numberofmotherkink = 0;
313   //printf("Number of vertices %d\n",numberofvertices);
314   for(Int_t ivertex=0; ivertex < numberofvertices; ivertex++) {
315     AliAODVertex *aodvertex = aodE->GetVertex(ivertex);
316     if(!aodvertex) continue;
317     //printf("Type %d\n",aodvertex->GetType());
318     if(aodvertex->GetType()==AliAODVertex::kKink) {
319       //printf("Find one kink\n");
320       AliAODTrack *mother = (AliAODTrack *) aodvertex->GetParent();
321       if(!mother) continue;
322       Int_t idmother = mother->GetID();
323       listofmotherkink[numberofmotherkink] = idmother;
324       //printf("ID %d\n",idmother);
325       numberofmotherkink++;
326     }
327   }
328   //printf("Number of kink mother in the events %d\n",numberofmotherkink);
329
330   //
331   // Loop on MC tracks only
332   //
333   AliAODMCParticle *mctrack(NULL);
334   // Monte-Carlo info
335   Int_t source(5);
336   if(mcthere){
337     AliDebug(1, "Loop MC tracks\n");
338     for(Int_t itrack = 0; itrack < fAODArrayMCInfo->GetEntriesFast(); itrack++) {
339       mctrack = (AliAODMCParticle *)(fAODArrayMCInfo->At(itrack));
340       if(!mctrack) continue;
341       AliHFEreducedMCParticle hfemcpart;
342       if(fTrackCuts->CheckParticleCuts(static_cast<UInt_t>(AliHFEcuts::kStepMCGenerated), mctrack)) continue;        
343       hfemcpart.SetSignal();
344       // Kinematics
345       hfemcpart.SetSignedPt(mctrack->Pt(), mctrack->Charge() > 0.);
346       hfemcpart.SetP(mctrack->P());
347       hfemcpart.SetEta(mctrack->Eta());
348       hfemcpart.SetPhi(mctrack->Phi());
349       hfemcpart.SetPdg(mctrack->GetPdgCode());
350       
351       // Get Production Vertex in radial direction
352       hfemcpart.SetProductionVertex(mctrack->Xv(),mctrack->Yv(),mctrack->Zv());
353
354       // Get Mother PDG code of the particle
355       Int_t motherlabel = TMath::Abs(mctrack->GetMother());
356       if(motherlabel >= 0 && motherlabel < fAODArrayMCInfo->GetEntriesFast()){
357         AliAODMCParticle *mother = dynamic_cast<AliAODMCParticle *>(fAODArrayMCInfo->At(motherlabel));
358         if(mother) hfemcpart.SetMotherPdg(mother->GetPdgCode());
359       }
360       
361       // derive source
362       source = 5;
363       if(fSignalCuts->IsCharmElectron(mctrack)) source = 0;
364       else if(fSignalCuts->IsBeautyElectron(mctrack)) source = 1;
365       else if(fSignalCuts->IsGammaElectron(mctrack)) source = 2;
366       else if(fSignalCuts->IsNonHFElectron(mctrack)) source = 3;
367       else if(TMath::Abs(mctrack->GetPdgCode()) == 11) source = 4;
368       else source = 5;
369       hfemcpart.SetSource(source);
370
371       fHFEevent->AddMCParticle(&hfemcpart);
372     }
373   }
374   
375   //
376   // Loop on reconstructed tracks
377   //
378   TArrayI arraytrack(fInputEvent->GetNumberOfTracks());
379   Int_t counterdc=0;
380   
381   AliAODTrack *track = 0x0;
382   AliDebug(1, "Loop reconstructed tracks\n");
383   for(Int_t itrack = 0; itrack < fInputEvent->GetNumberOfTracks(); itrack++){
384     // Run track loop
385     track = dynamic_cast<AliAODTrack *>(fInputEvent->GetTrack(itrack));
386     if(!track) continue;
387     // Cut track (Only basic track cuts)
388     // printf("testv\n");
389     if(!fTrackCuts->CheckParticleCuts(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepRecKineITSTPC, track)) continue;
390     //
391     //printf("testu\n");
392
393     // Kinematics
394     AliHFEreducedTrack hfetrack;
395     hfetrack.SetSignedPt(track->Pt(), track->Charge() > 0);
396     hfetrack.SetP(track->P());
397     hfetrack.SetEta(track->Eta());
398     hfetrack.SetPhi(track->Phi());
399     hfetrack.SetTPCmomentum(track->GetDetPid() ? track->GetDetPid()->GetTPCmomentum() : track->P());
400
401     // Track ID
402     hfetrack.SetTrackID(track->GetID());
403
404     // status
405     ULong_t status = track->GetStatus();
406     if((status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) hfetrack.SetITSrefit();
407     if((status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) hfetrack.SetTPCrefit();
408     if((status & AliVTrack::kTOFpid) == AliVTrack::kTOFpid) hfetrack.SetTOFpid();
409     //if((status & AliVTrack::kTOFmismatch) == AliVTrack::kTOFmismatch) hfetrack.SetTOFmismatch();
410     if(IsTOFmismatch(track, pid)) hfetrack.SetTOFmismatch(); // New version suggested by Pietro Antonioli
411     if(track->IsEMCAL()) hfetrack.SetEMCALpid();
412     // Filter
413     for(Int_t k=0; k<20; k++) {
414       Int_t u = 1<<k;     
415       if((track->TestFilterBit(u))) {
416               hfetrack.SetFilterBit(k);
417       } 
418     }
419
420     if(mcthere){
421       // Fill Monte-Carlo Information
422       Int_t label = TMath::Abs(track->GetLabel());
423       if(label < fAODArrayMCInfo->GetEntriesFast())
424         mctrack = dynamic_cast<AliAODMCParticle *>(fAODArrayMCInfo->At(label));
425       if(mctrack){ 
426         AliDebug(2, "Associated MC particle found\n");
427         if(fTrackCuts->CheckParticleCuts(static_cast<UInt_t>(AliHFEcuts::kStepMCGenerated), mctrack)) hfetrack.SetMCSignal();
428         // Kinematics
429         hfetrack.SetMCSignedPt(mctrack->Pt(),mctrack->Charge() > 0.);
430         hfetrack.SetMCP(mctrack->P());
431         hfetrack.SetMCEta(mctrack->Eta());
432         hfetrack.SetMCPhi(mctrack->Phi());
433         hfetrack.SetMCPDG(mctrack->GetPdgCode());
434       
435         // Get Production Vertex in radial direction
436         hfetrack.SetMCProdVtx(mctrack->Xv(),mctrack->Yv(),mctrack->Zv());
437       
438         // Get Mother PDG code of the particle
439         Int_t motherlabel = TMath::Abs(mctrack->GetMother());
440         if(motherlabel >= 0 && motherlabel < fAODArrayMCInfo->GetEntriesFast()){
441           AliAODMCParticle *mother = dynamic_cast<AliAODMCParticle *>(fAODArrayMCInfo->At(motherlabel));
442           if(mother) hfetrack.SetMCMotherPdg(mother->GetPdgCode());
443         }
444       
445         // derive source
446         source = 5;
447         if(fSignalCuts->IsCharmElectron(mctrack)) source = 0;
448         else if(fSignalCuts->IsBeautyElectron(mctrack)) source = 1;
449         else if(fSignalCuts->IsGammaElectron(mctrack)) source = 2;
450         else if(fSignalCuts->IsNonHFElectron(mctrack)) source = 3;
451         else if(TMath::Abs(mctrack->GetPdgCode()) == 11) source = 4;
452         else source = 5;
453         hfetrack.SetMCSource(source); 
454       } else {
455         AliDebug(2, "Associated MC particle not found\n");
456       }
457     }
458
459     // HFE DCA
460     Float_t dcaxy = -999.,
461             dcaz = -999.;
462     fExtraCuts->GetImpactParameters((AliVTrack *)track,dcaxy,dcaz);
463     hfetrack.SetDCA(dcaxy, dcaz);
464     Double_t hfeImpactParam(-999.), hfeImpactParamResol(-999.);
465     fExtraCuts->GetHFEImpactParameters((AliVTrack *)track,hfeImpactParam,hfeImpactParamResol);
466     hfetrack.SetHFEImpactParam(hfeImpactParam,hfeImpactParamResol);
467
468     // Different number of clusters definitions
469     Int_t nclustersITS(track->GetITSNcls()),
470           nclustersTPC(track->GetTPCNcls()),
471           nclustersTPCall(track->GetTPCClusterMap().CountBits()),
472           nclustersTPCshared(0);
473     UChar_t nfindableTPC = track->GetTPCNclsF();
474     const TBits &sharedTPC = track->GetTPCSharedMap();
475     for(Int_t ibit = 0; ibit < 160; ibit++) if(sharedTPC.TestBitNumber(ibit)) nclustersTPCshared++;
476     hfetrack.SetITSnclusters(nclustersITS);
477     hfetrack.SetTPCnclusters(track->GetTPCNcls());
478     hfetrack.SetTRDnclusters(track->GetTRDncls());
479     hfetrack.SetTPCnclustersPID(track->GetTPCsignalN());
480     hfetrack.SetTPCcrossedRows(track->GetTPCNCrossedRows());
481     hfetrack.SetTPCnclustersAll(nclustersTPCall);
482     hfetrack.SetTPCsharedClusters(nclustersTPCshared);
483     hfetrack.SetTPCclusterRatio(nfindableTPC ? static_cast<Float_t>(nclustersTPC)/static_cast<Float_t>(nfindableTPC) : 0);
484     hfetrack.SetTPCclusterRatioAll(nfindableTPC ? static_cast<Float_t>(nclustersTPCall)/static_cast<Float_t>(nfindableTPC) : 0);
485     UChar_t itsPixel = track->GetITSClusterMap();
486     for(int ily = 0; ily < 6; ily++) 
487             if(TESTBIT(itsPixel, ily)) hfetrack.SetITScluster(ily);
488    
489     // TRD related quantities (Yvonne)
490     hfetrack.SetTRDntrackletsPID(track->GetTRDntrackletsPID());
491     hfetrack.SetTRDnslices(track->GetDetPid()->GetTRDnSlices());
492     hfetrack.SetTRDchi2(track->GetTRDchi2());
493     AliAODPid* aodpid= track->GetDetPid();
494     Double_t* arraytrdsignals;
495     arraytrdsignals=aodpid->GetTRDslices();
496     Int_t nslicetemp=0;
497     for(Int_t iplane = 0; iplane < 6; iplane++){
498             nslicetemp=0;
499             for(Int_t b=(iplane*8);b<((iplane*8)+8);b++){
500               if(track->GetTRDntrackletsPID()>0){
501                       if(arraytrdsignals[b]>0.001) nslicetemp++;
502               }
503             }
504             if(nslicetemp > 0) hfetrack.SetTRDstatus(iplane);
505     }
506
507
508     // Kink
509     if(fExtraCuts->IsKinkDaughter(track)) hfetrack.SetIsKinkDaughter();
510
511     // kink mother
512     for(Int_t kinkmother = 0; kinkmother < numberofmotherkink; kinkmother++) {
513       if(track->GetID() == listofmotherkink[kinkmother]) {
514               hfetrack.SetIsKinkMother();
515               continue;
516       }
517     }
518     
519     
520     // Double counted
521     Int_t id(track->GetID());
522     for(Int_t l=0; l < counterdc; l++){
523       Int_t iTrack2 = arraytrack.At(l);
524       if(iTrack2==id){
525          hfetrack.SetDoubleCounted();
526          break;
527       }
528     }
529     // Add the id at this place
530     arraytrack.AddAt(id,counterdc);
531     counterdc++;
532
533     // PID
534     hfetrack.SetTPCdEdx(track->GetDetPid() ? track->GetDetPid()->GetTPCsignal() : 0.);
535     hfetrack.SetTPCsigmaEl(pid->NumberOfSigmasTPC(track, AliPID::kElectron));
536     hfetrack.SetTOFsigmaEl(pid->NumberOfSigmasTOF(track, AliPID::kElectron));
537     hfetrack.SetTOFmismatchProbability(pid->GetTOFMismatchProbability(track));
538     hfetrack.SetITSsigmaEl(pid->NumberOfSigmasITS(track, AliPID::kElectron));
539     // Eta correction
540     copyTrack.~AliAODTrack();
541     new(&copyTrack) AliAODTrack(*track);
542     if(fTPCpid->HasCentralityCorrection()) fTPCpid->ApplyCentralityCorrection(&copyTrack, static_cast<Double_t>(ncontrib),AliHFEpidObject::kAODanalysis);
543     if(fTPCpid->HasEtaCorrection()) fTPCpid->ApplyEtaCorrection(&copyTrack, AliHFEpidObject::kAODanalysis);
544     hfetrack.SetTPCsigmaElCorrected(pid->NumberOfSigmasTPC(&copyTrack, AliPID::kElectron));
545     hfetrack.SetTPCdEdxCorrected(copyTrack.GetDetPid() ? copyTrack.GetDetPid()->GetTPCsignal() : 0.);
546     if(track->IsEMCAL()){
547       // EMCAL cluster
548       Double_t emcalEnergyOverP = -1.,
549                showershape[4] = {0.,0.,0.,0.};
550       hfetrack.SetEMCALSigmaEl(pid->NumberOfSigmasEMCAL(track, AliPID::kElectron, emcalEnergyOverP, &showershape[0]));
551       hfetrack.SetEMCALEoverP(emcalEnergyOverP);
552       hfetrack.SetEMCALShowerShape(showershape);
553     }
554
555     // If TOF cut
556     if(fNbOfTOFSigma>0.0){
557       AliDebug(1, "TOF cut\n");
558       if(!((status & AliVTrack::kTOFpid) == AliVTrack::kTOFpid)) continue; 
559       if(TMath::Abs(pid->NumberOfSigmasTOF(track, AliPID::kElectron))> fNbOfTOFSigma) continue;
560     }
561
562     // Track finished, add NOW to the Event
563     fHFEevent->AddTrack(&hfetrack);
564     //printf("after\n");
565   }
566   
567   // Fill the debug tree
568   AliInfo(Form("Number of tracks: %d\n", fHFEevent->GetNumberOfTracks()));
569   AliInfo(Form("Number of MC particles: %d\n", fHFEevent->GetNumberOfMCParticles()));
570   fHFEtree->Fill();
571
572   fEventNumber++;
573   PostData(1, fHFEtree);
574 }
575
576 void AliHFEreducedEventCreatorAOD::Terminate(Option_t *){
577   //
578   // Terminate
579   //
580   AliInfo("terminating...\n");
581
582 }
583
584 Bool_t AliHFEreducedEventCreatorAOD::IsTOFmismatch(const AliVTrack *const track, const AliPIDResponse *const pid) const {
585   //
586   // is TOF mismatch
587   //
588   Double_t probs[AliPID::kSPECIESC];
589   AliPIDResponse::EDetPidStatus status = pid->ComputeTOFProbability(track, AliPID::kSPECIESC, probs);
590   return status == AliPIDResponse::kDetMismatch;
591 }
592