PWGHF/hfe/AliHFEdebugTreeTaskAOD.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 //
  18 // Authors:
  19 //   R.Bailhache <R.Bailhache@gsi.de>
  20 //
  21 #include <TBits.h>
  22 #include <TString.h>
  23 #include <TArrayI.h>
  24
  25 #include "AliAnalysisManager.h"
  26 #include "AliCentrality.h"
  27 #include "AliAODTrack.h"
  28 #include "AliAODEvent.h"
  29 #include "AliHFEcuts.h"
  30 #include "AliHFEextraCuts.h"
  31 #include "AliInputEventHandler.h"
  32 #include "AliLog.h"
  33 #include "AliPIDResponse.h"
  34 #include "AliVEvent.h"
  35 #include "AliHFEpidTPC.h"
  36 #include "TTreeStream.h"
  37 #include "AliESDtrack.h"
  38 #include "TClonesArray.h"
  39 #include "AliAODMCHeader.h"
  40 #include "AliHFEsignalCuts.h"
  41 #include "AliAODMCParticle.h"
  42 #include "AliVTrack.h"
  43 #include "AliAODVertex.h"
  44
  45 #include "AliHFEdebugTreeTaskAOD.h"
  46
  47 ClassImp(AliHFEdebugTreeTaskAOD)
  48
  49 AliHFEdebugTreeTaskAOD::AliHFEdebugTreeTaskAOD():
  50   AliAnalysisTaskSE(),
  51   fAODMCHeader(NULL),
  52   fAODArrayMCInfo(NULL),
  53   fTrackCuts(NULL),
  54   fExtraCuts(NULL),
  55   fSignalCuts(NULL),
  56   fTPCpid(NULL),
  57   fEventNumber(0),
  58   fNclustersTPC(70),
  59   fNclustersTPCPID(0),
  60   fNclustersITS(2),
  61   fFilename("HFEtree.root"),
  62   fDebugTree(NULL)
  63 {
  64
  65 }
  66
  67 AliHFEdebugTreeTaskAOD::AliHFEdebugTreeTaskAOD(const char *name):
  68   AliAnalysisTaskSE(name),
  69   fAODMCHeader(NULL),
  70   fAODArrayMCInfo(NULL),
  71   fTrackCuts(NULL),
  72   fExtraCuts(NULL),
  73   fSignalCuts(NULL),
  74   fTPCpid(NULL),
  75   fEventNumber(0),
  76   fNclustersTPC(70),
  77   fNclustersTPCPID(0),
  78   fNclustersITS(2),
  79   fFilename("HFEtree.root"),
  80   fDebugTree(NULL)
  81 {
  82   fTPCpid = new AliHFEpidTPC("QAtpcPID");
  83 }
  84
  85 AliHFEdebugTreeTaskAOD::~AliHFEdebugTreeTaskAOD(){
  86
  87     if(fDebugTree) delete fDebugTree;
  88     if(fTPCpid) delete fTPCpid;
  89 }
  90
  91 void AliHFEdebugTreeTaskAOD::UserCreateOutputObjects(){
  92   //
  93   // Create debug tree, signal cuts and track cuts
  94   //
  95
  96   //printf("test\n");
  97   fDebugTree = new TTreeSRedirector(fFilename.Data());
  98
  99   //printf("testa\n");
 100   fSignalCuts = new AliHFEsignalCuts("HFEsignalCuts", "HFE MC Signal definition");
 101   //printf("testb\n");
 102
 103   fTrackCuts = new AliHFEcuts("fTrackCuts", "Basic HFE track cuts");
 104   fTrackCuts->CreateStandardCuts();
 105   fTrackCuts->SetAOD();
 106   // Track cuts
 107   fTrackCuts->SetMinNClustersTPC(fNclustersTPC);
 108   fTrackCuts->SetMinRatioTPCclusters(0);
 109   fTrackCuts->SetTPCmodes(AliHFEextraCuts::kFound, AliHFEextraCuts::kFoundOverFindable);
 110   fTrackCuts->SetMinNClustersTPCPID(fNclustersTPCPID);
 111   fTrackCuts->SetMinNClustersITS(fNclustersITS);
 112   // Event cuts
 113   fTrackCuts->SetUseMixedVertex(kTRUE);
 114   fTrackCuts->SetVertexRange(10.);
 115   //printf("testa\n");
 116   fTrackCuts->Initialize();
 117   //printf("testb\n");
 118
 119   fExtraCuts = new AliHFEextraCuts("hfeExtraCuts","HFE Extra Cuts");
 120
 121 }
 122
 123 void AliHFEdebugTreeTaskAOD::UserExec(Option_t *){
 124   //
 125   // User Exec: Fill debug Tree
 126   //
 127
 128   // Get PID response
 129   AliPIDResponse *pid = NULL;
 130   AliInputEventHandler *handler = dynamic_cast<AliInputEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 131   if(handler){
 132     //printf("testb\n");
 133     pid = handler->GetPIDResponse();
 134   } else {
 135     AliError("No Handler");
 136   }
 137   if(!pid){
 138     printf("testc\n");
 139     AliError("No PID response");
 140     return;
 141   }
 142   if(!fInputEvent) {
 143     AliError("No Input event");
 144     return;
 145   }
 146
 147   // MC info
 148   Bool_t mcthere = kTRUE;
 149   AliAODEvent *aodE = dynamic_cast<AliAODEvent *>(fInputEvent);
 150   if(!aodE){
 151     AliError("No AOD Event");
 152     return;
 153   }
 154   fAODMCHeader = dynamic_cast<AliAODMCHeader *>(fInputEvent->FindListObject(AliAODMCHeader::StdBranchName()));
 155   if(!fAODMCHeader){
 156     mcthere = kFALSE;
 157     return;
 158   }
 159   fAODArrayMCInfo = dynamic_cast<TClonesArray *>(fInputEvent->FindListObject(AliAODMCParticle::StdBranchName()));
 160   if(!fAODArrayMCInfo){
 161     mcthere = kFALSE;
 162     return;
 163   }
 164   else {
 165     fSignalCuts->SetMCAODInfo(fAODArrayMCInfo);
 166     fTrackCuts->SetMCEvent(aodE);
 167   }
 168
 169   // Set for track cuts
 170   fTrackCuts->SetRecEvent(fInputEvent);
 171
 172   if(!fTrackCuts->CheckEventCuts("fCutsEvRec", fInputEvent)){
 173     AliDebug(1, "Event rejected by the event cuts\n");
 174     return;
 175   }
 176   if(!fExtraCuts){
 177     fExtraCuts = new AliHFEextraCuts("hfeExtraCuts","HFE Extra Cuts");
 178   }
 179   fExtraCuts->SetRecEventInfo(fInputEvent);
 180
 181   // Get run number
 182   Int_t run = fInputEvent->GetRunNumber();
 183
 184   // Derive trigger
 185   UInt_t trigger = fInputHandler->IsEventSelected();
 186   Bool_t isMBTrigger = trigger & AliVEvent::kMB;
 187   Bool_t isCentralTrigger = trigger & AliVEvent::kCentral;
 188   Bool_t isSemicentralTrigger = trigger & AliVEvent::kSemiCentral;
 189   Bool_t isEMCALTrigger = trigger & AliVEvent::kEMCEJE;
 190
 191   // Get Primary Vertex
 192   const AliVVertex *vertex = fInputEvent->GetPrimaryVertex();
 193   Double_t vtx[3];
 194   vertex->GetXYZ(vtx);
 195   Double_t ncontrib = fInputEvent->GetPrimaryVertex()->GetNContributors();
 196
 197   // Get centrality
 198   Float_t centrality = -1.;
 199   AliCentrality *hicent = fInputEvent->GetCentrality();
 200   centrality = hicent->GetCentralityPercentile("V0M");
 201
 202   // Look for kink mother
 203   Int_t numberofvertices = aodE->GetNumberOfVertices();
 204   Double_t listofmotherkink[numberofvertices];
 205   Int_t numberofmotherkink = 0;
 206   //printf("Number of vertices %d\n",numberofvertices);
 207   for(Int_t ivertex=0; ivertex < numberofvertices; ivertex++) {
 208     AliAODVertex *aodvertex = aodE->GetVertex(ivertex);
 209     if(!aodvertex) continue;
 210     //printf("Type %d\n",aodvertex->GetType());
 211     if(aodvertex->GetType()==AliAODVertex::kKink) {
 212       //printf("Find one kink\n");
 213       AliAODTrack *mother = (AliAODTrack *) aodvertex->GetParent();
 214       if(!mother) continue;
 215       Int_t idmother = mother->GetID();
 216       listofmotherkink[numberofmotherkink] = idmother;
 217       //printf("ID %d\n",idmother);
 218       numberofmotherkink++;
 219     }
 220   }
 221   //printf("Number of kink mother in the events %d\n",numberofmotherkink);
 222
 223   // Common variables
 224   Double_t charge, eta, phi, momentum, momentumTPC, transversemomentum;
 225
 226   //
 227   // Loop on reconstructed tracks
 228   //
 229
 230   TArrayI *arraytrack = new TArrayI(fInputEvent->GetNumberOfTracks());
 231
 232   AliAODTrack *track = 0x0;
 233   AliAODMCParticle *mctrack = NULL;
 234   for(Int_t itrack = 0; itrack < fInputEvent->GetNumberOfTracks(); itrack++){
 235     // fill the tree
 236     track = dynamic_cast<AliAODTrack *>(fInputEvent->GetTrack(itrack));
 237     if(!track) continue;
 238     // Cut track (Only basic track cuts)
 239     //printf("testv\n");
 240     if(!fTrackCuts->CheckParticleCuts(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepRecKineITSTPC, track)) continue;
 241     //
 242     //printf("testu\n");
 243     Double_t nSigmaTOF = pid->NumberOfSigmasTOF(track, AliPID::kElectron);
 244     Double_t nSigmaTPC = pid->NumberOfSigmasTPC(track, AliPID::kElectron);
 245     Double_t tPCdEdx = track->GetDetPid() ? track->GetDetPid()->GetTPCsignal() : 0.;
 246
 247     // Kinematics
 248     charge = track->Charge() > 0 ? 1. : -1.;
 249     eta = track->Eta();
 250     phi = track->Phi();
 251     momentum = track->P() * charge;
 252     transversemomentum = track->Pt() * charge;
 253     momentumTPC = track->GetDetPid() ? track->GetDetPid()->GetTPCmomentum() : track->P();
 254
 255     // status
 256     ULong_t status = track->GetStatus();
 257     Int_t itsrefit=0;
 258     if((status & AliESDtrack::kITSrefit) == AliESDtrack::kITSrefit) itsrefit = 1;
 259     Int_t tpcrefit=0;
 260     if((status & AliESDtrack::kTPCrefit) == AliESDtrack::kTPCrefit) tpcrefit = 1;
 261
 262     // ITS number of clusters
 263     UChar_t nclustersITS = track->GetITSNcls();
 264     // TPC number of clusters (different definitions)
 265     UChar_t nclustersTPCfit = track->GetTPCNcls();
 266     UChar_t nclustersTPCall = 0;
 267     const TBits &clusterTPC = track->GetTPCClusterMap();
 268     nclustersTPCall = clusterTPC.CountBits();
 269     UChar_t nclustersTPCPID = track->GetTPCsignalN();
 270     UChar_t nfindableTPC =  track->GetTPCNclsF();
 271     Double_t clusterRatioTPCfit = 0.0;
 272     if((static_cast<Double_t>(nfindableTPC))>0.0) clusterRatioTPCfit = static_cast<Double_t>(nclustersTPCfit)/static_cast<Double_t>(nfindableTPC);
 273     Double_t clusterRatioTPCall = 0.0;
 274     if((static_cast<Double_t>(nfindableTPC))>0.0) clusterRatioTPCall = static_cast<Double_t>(nclustersTPCall)/static_cast<Double_t>(nfindableTPC);
 275     UChar_t nclustersTPCshared = 0;
 276     Float_t ncrossedRowsTPC = track->GetTPCNCrossedRows();
 277     const TBits &sharedTPC = track->GetTPCSharedMap();
 278     for(Int_t ibit = 0; ibit < 160; ibit++) if(sharedTPC.TestBitNumber(ibit)) nclustersTPCshared++;
 279     // TRD clusters and tracklets
 280     UChar_t nclustersTRD = track->GetTRDncls();
 281     UChar_t ntrackletsTRDPID = track->GetTRDntrackletsPID();
 282     Int_t   nslicesTRD = track->GetNumberOfTRDslices();
 283     Int_t   chi2TRD = track->GetTRDchi2();
 284     // ITS and TRD acceptance maps
 285     UChar_t itsPixel = track->GetITSClusterMap();
 286     Bool_t statusL0 = kFALSE;
 287     if(TESTBIT(itsPixel, 0)) statusL0 = kTRUE;
 288     Bool_t statusL1 = kFALSE;
 289     if(TESTBIT(itsPixel, 1)) statusL1 = kTRUE;
 290
 291     // HFE DCA
 292     Float_t dcaxy = -999.;
 293     Float_t dcaz = -999.;
 294     fExtraCuts->GetImpactParameters((AliVTrack *)track,dcaxy,dcaz);
 295
 296     // Kink
 297     Int_t kink = 0;
 298     if(fExtraCuts->IsKinkDaughter(track)) kink = 1;
 299
 300     // kink mother
 301     Int_t kinkmotherpass = 0;
 302     for(Int_t kinkmother = 0; kinkmother < numberofmotherkink; kinkmother++) {
 303       if(track->GetID() == listofmotherkink[kinkmother]) {
 304         kinkmotherpass = 1;
 305         continue;
 306       }
 307     }
 308
 309     // ID track to see if negative
 310     Int_t id = track->GetID();
 311
 312     // Double counted
 313     Int_t doublec = 0;
 314     for(Int_t l=0; l < itrack; l++){
 315       Int_t iTrack2 = arraytrack->At(l);
 316       if(iTrack2==id) doublec=1;
 317     }
 318     //printf("Doublec %d\n",doublec);
 319
 320     // Add the id at this place
 321     arraytrack->AddAt(id,itrack);
 322
 323     // Filter
 324     Int_t filter[20];
 325     for(Int_t k=0; k<20; k++) {
 326       filter[k]=0;
 327       Int_t u = 1<<k;
 328       if((track->TestFilterBit(u))) {
 329         filter[k] = 1;
 330       }
 331     }
 332     Int_t filter0 = filter[0];
 333     Int_t filter1 = filter[1];
 334     Int_t filter2 = filter[2];
 335     Int_t filter3 = filter[3];
 336     Int_t filter4 = filter[4];
 337     Int_t filter5 = filter[5];
 338     Int_t filter6 = filter[6];
 339     Int_t filter7 = filter[7];
 340     Int_t filter8 = filter[8];
 341     Int_t filter9 = filter[9];
 342     Int_t filter10 = filter[10];
 343     Int_t filter11 = filter[11];
 344     Int_t filter12 = filter[12];
 345     Int_t filter13 = filter[13];
 346     Int_t filter14 = filter[14];
 347     Int_t filter15 = filter[15];
 348     Int_t filter16 = filter[16];
 349     Int_t filter17 = filter[17];
 350     Int_t filter18 = filter[18];
 351     Int_t filter19 = filter[19];
 352
 353     Int_t eventnb = fEventNumber;
 354
 355     //printf("track\n");
 356
 357     // Monte-Carlo info
 358     Double_t eR,vx,vy,vz;
 359     Double_t chargemc, etamc, phimc, momentummc, transversemomentummc;
 360     Int_t source,pdg,signal;
 361     signal = 0;
 362     if(mcthere){
 363       Int_t label = TMath::Abs(track->GetLabel());
 364       if(label && label < fAODArrayMCInfo->GetEntriesFast())
 365         mctrack = dynamic_cast<AliAODMCParticle *>(fAODArrayMCInfo->At(label));
 366       if(!mctrack) continue;
 367       if(fTrackCuts->CheckParticleCuts(static_cast<UInt_t>(AliHFEcuts::kStepMCGenerated), mctrack)) signal = 1;
 368       // Kinematics
 369       chargemc = mctrack->Charge() > 0. ? 1. : -1.;
 370       momentummc = mctrack->P() * charge;
 371       transversemomentummc = mctrack->Pt() * charge;
 372       etamc = mctrack->Eta();
 373       phimc = mctrack->Phi();
 374       pdg = mctrack->GetPdgCode();
 375
 376       // Get Production Vertex in radial direction
 377       vx = mctrack->Xv();
 378       vy = mctrack->Yv();
 379       vz = mctrack->Zv();
 380       eR = TMath::Sqrt(vx*vx+vy*vy);
 381
 382       // Get Mother PDG code of the particle
 383       Int_t motherPdg = 0;
 384       Int_t motherlabel = TMath::Abs(mctrack->GetMother());
 385       if(motherlabel >= 0 && motherlabel < fAODArrayMCInfo->GetEntriesFast()){
 386         AliAODMCParticle *mother = dynamic_cast<AliAODMCParticle *>(fAODArrayMCInfo->At(motherlabel));
 387         if(mother) motherPdg = mother->GetPdgCode();
 388       }
 389
 390       // derive source
 391       source = 5;
 392       if(fSignalCuts->IsCharmElectron(mctrack)) source = 0;
 393       else if(fSignalCuts->IsBeautyElectron(mctrack)) source = 1;
 394       else if(fSignalCuts->IsGammaElectron(mctrack)) source = 2;
 395       else if(fSignalCuts->IsNonHFElectron(mctrack)) source = 3;
 396       else if(TMath::Abs(pdg) == 11) source = 4;
 397       else source = 5;
 398
 399     }
 400
 401
 402     // Fill Tree
 403     //printf("Fill\n");
 404     (*fDebugTree) << "PIDdebug"
 405                   << "centrality="          << centrality
 406                   << "MBtrigger="           << isMBTrigger
 407                   << "CentralTrigger="      << isCentralTrigger
 408                   << "SemicentralTrigger="  << isSemicentralTrigger
 409                   << "EMCALtrigger="        << isEMCALTrigger
 410                   << "run="                 << run
 411                   << "eventnb="             << eventnb
 412                   << "vx="                  << vtx[0]
 413                   << "vy="                  << vtx[1]
 414                   << "vz="                  << vtx[2]
 415                   << "ncontrib="            << ncontrib
 416                   << "id="                  << id
 417                   << "dc="                  << doublec
 418                   << "p="                   << momentum
 419                   << "ptpc="                << momentumTPC
 420                   << "pt="                  << transversemomentum
 421                   << "eta="                 << eta
 422                   << "phi="                 << phi
 423                   << "itsrefit="            << itsrefit
 424                   << "tpcrefit="            << tpcrefit
 425                   << "nclustersTPC="        << nclustersTPCfit
 426                   << "nclustersTPCall="     << nclustersTPCall
 427                   << "nclustersTPCPID="     << nclustersTPCPID
 428                   << "nclustersTPCshared="  << nclustersTPCshared
 429                   << "ncrossedRowsTPC="     << ncrossedRowsTPC
 430                   << "clusterRatioTPC="     << clusterRatioTPCfit
 431                   << "clusterRatioTPCall="  << clusterRatioTPCall
 432                   << "nclustersITS="        << nclustersITS
 433                   << "nclustersTRD="        << nclustersTRD
 434                   << "ntrackletsTRD="       << ntrackletsTRDPID
 435                   << "nslicesTRD="          << nslicesTRD
 436                   << "chi2TRD="             << chi2TRD
 437                   << "statusITS0="          << statusL0
 438                   << "statusITS1="          << statusL1
 439                   << "TOFsigmaEl="          << nSigmaTOF
 440                   << "TPCsigmaEl="          << nSigmaTPC
 441                   << "TPCdEdx="             << tPCdEdx
 442                   << "dcaR="                << dcaxy
 443                   << "dcaZ="                << dcaz
 444                   << "kinkdaughter="        << kink
 445                   << "kinkmother="          << kinkmotherpass
 446                   << "nbofmotherkink="      << numberofmotherkink
 447                   << "filter0="             << filter0
 448                   << "filter1="             << filter1
 449                   << "filter2="             << filter2
 450                   << "filter3="             << filter3
 451                   << "filter4="             << filter4
 452                   << "filter5="             << filter5
 453                   << "filter6="             << filter6
 454                   << "filter7="             << filter7
 455                   << "filter8="             << filter8
 456                   << "filter9="             << filter9
 457                   << "filter10="            << filter10
 458                   << "filter11="            << filter11
 459                   << "filter12="            << filter12
 460                   << "filter13="            << filter13
 461                   << "filter14="            << filter14
 462                   << "filter15="            << filter15
 463                   << "filter16="            << filter16
 464                   << "filter17="            << filter17
 465                   << "filter18="            << filter18
 466                   << "filter19="            << filter19
 467                   << "mcp="                 << momentummc
 468                   << "mcpt="                << transversemomentummc
 469                   << "mceta="               << etamc
 470                   << "mcphi="               << phimc
 471                   << "mcpdg="               << pdg
 472                   << "source="              << source
 473                   << "px="                  << vx
 474                   << "py="                  << vy
 475                   << "pz="                  << vz
 476                   << "eR="                  << eR
 477                   << "mccharge="            << chargemc
 478                   << "signal="              << signal
 479                   << "\n";
 480
 481
 482     //printf("after\n");
 483
 484   }
 485
 486   arraytrack->~TArrayI();
 487   fEventNumber++;
 488
 489
 490 }
 491
 492
 493 void AliHFEdebugTreeTaskAOD::SetFileName(const char *filename){ fFilename = filename; }
 494 void AliHFEdebugTreeTaskAOD::Terminate(Option_t *){
 495
 496   if(fDebugTree) delete fDebugTree;
 497   fDebugTree=0x0;
 498
 499 }
 500