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