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