PWGCF/Correlations/Base/AliAnalysisTaskCFTree.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 // Analysis task to produce trees of lightweight events
  17 // evgeny.kryshen@cern.ch
  18
  19 // aliroot
  20 #include "AliAnalysisTaskCFTree.h"
  21 #include "AliCFParticle.h"
  22 #include "AliAnalysisManager.h"
  23 #include "AliInputEventHandler.h"
  24 #include "AliESDEvent.h"
  25 #include "AliAODEvent.h"
  26 #include "AliMCEvent.h"
  27 #include "AliHeader.h"
  28 #include "AliESDHeader.h"
  29 #include "AliAODHeader.h"
  30 #include "AliAODMCHeader.h"
  31 #include "AliGenEventHeader.h"
  32 #include "AliGenCocktailEventHeader.h"
  33 #include "AliAODTrack.h"
  34 #include "AliESDtrack.h"
  35 #include "AliMCParticle.h"
  36 #include "AliAODMCParticle.h"
  37 #include "AliExternalTrackParam.h"
  38 #include "AliCentrality.h"
  39 #include "AliAnalysisFilter.h"
  40 #include "AliPIDResponse.h"
  41 #include "AliTPCPIDResponse.h"
  42 #include "AliAODTracklets.h"
  43 // root
  44 #include "TMath.h"
  45 #include "TFile.h"
  46 #include "TList.h"
  47 #include "TH1I.h"
  48 #include "TChain.h"
  49 #include "TTree.h"
  50 #include "TClonesArray.h"
  51 ClassImp(AliAnalysisTaskCFTree)
  52
  53 //-----------------------------------------------------------------------------
  54 AliAnalysisTaskCFTree::AliAnalysisTaskCFTree(const char* name) :
  55   AliAnalysisTask(name,""),
  56   fDebug(0),
  57   fMode(0),
  58   fIsAOD(1),
  59   fInputHandler(0x0),
  60   fMcHandler(0x0),
  61   fTrackFilter(0x0),
  62   fHybridConstrainedMask(0),
  63   fTPConlyConstrainedMask(0),
  64   fPIDResponse(0x0),
  65   fListOfHistos(0x0),
  66   fEventStatistics(0x0),
  67   fTree(0x0),
  68   fParticles(0x0),
  69   fTracklets(0x0),
  70   fMuons(0x0),
  71   fField(0),
  72   fCentrality(0),
  73   fZvtx(0),
  74   fRunNumber(0),
  75   fPeriod(0),
  76   fOrbit(),
  77   fBc(),
  78   fSelectMask(0),
  79   fSelectBit(AliVEvent::kMB),
  80   fZVertexCut(10.),
  81   fnTracksVertex(1),
  82   fCentralityMethod("V0M"),
  83   fTrackFilterBit(0xffffffff),
  84   fTrackEtaCut(1.0),
  85   fPtMin(0.15),
  86   fSharedClusterCut(0.4),
  87   fCrossedRowsCut(100),
  88   fFoundFractionCut(0.8),
  89   fDphiCut(1.e9),
  90   fStoreTracks(0),
  91   fStoreTracklets(0),
  92   fStoreMuons(0)
  93 {
  94   Info("AliAnalysisTaskCFTree","Calling Constructor");
  95   DefineInput(0,TChain::Class());
  96   DefineOutput(0,TList::Class());
  97   DefineOutput(1,TTree::Class());
  98 }
  99 //-----------------------------------------------------------------------------
 100
 101
 102 //-----------------------------------------------------------------------------
 103 void AliAnalysisTaskCFTree::ConnectInputData(Option_t* /*option*/){
 104   fInputHandler = dynamic_cast<AliInputEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 105   fMcHandler    = dynamic_cast<AliInputEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 106   if (fInputHandler) fInputHandler->SetNeedField();
 107 }
 108 //-----------------------------------------------------------------------------
 109
 110
 111 //-----------------------------------------------------------------------------
 112 void AliAnalysisTaskCFTree::CreateOutputObjects(){
 113   fListOfHistos = new TList();
 114   fListOfHistos->SetOwner();
 115   fEventStatistics = new TH1I("fEventStatistics","",10,0,10);
 116   fEventStatistics->SetBit(TH1::kCanRebin);
 117
 118   fListOfHistos->Add(fEventStatistics);
 119
 120   if (fStoreTracks)    fParticles = new TClonesArray("AliCFParticle",2000);
 121   if (fStoreTracklets) fTracklets = new TClonesArray("AliCFParticle",100);
 122   if (fStoreMuons)     fMuons     = new TClonesArray("AliCFParticle",100);
 123   // create file-resident tree
 124   TDirectory *owd = gDirectory;
 125   OpenFile(1);
 126   fTree = new TTree("events","events");
 127   owd->cd();
 128   fTree->Branch("cent",&fCentrality);
 129   fTree->Branch("zvtx",&fZvtx);
 130   fTree->Branch("field",&fField);
 131   fTree->Branch("run",&fRunNumber);
 132   fTree->Branch("period",&fPeriod);
 133   fTree->Branch("orbit",&fOrbit);
 134   fTree->Branch("bc",&fBc);
 135   fTree->Branch("mask",&fSelectMask);
 136   if (fStoreTracks)    fTree->Branch("particles",&fParticles);
 137   if (fStoreTracklets) fTree->Branch("tracklets",&fTracklets);
 138   if (fStoreMuons)     fTree->Branch("muons",&fMuons);
 139   PostData(0,fListOfHistos);
 140   PostData(1,fTree);
 141 }
 142 //-----------------------------------------------------------------------------
 143
 144
 145 //-----------------------------------------------------------------------------
 146 void AliAnalysisTaskCFTree::Exec(Option_t *){
 147   if (fParticles) fParticles->Clear();
 148   if (fTracklets) fTracklets->Clear();
 149   if (fMuons)     fMuons->Clear();
 150
 151   fEventStatistics->Fill("before cuts",1);
 152
 153   AliVEvent* event = fInputHandler->GetEvent();
 154   if (!event) return;
 155   fEventStatistics->Fill("after event check",1);
 156
 157   fPIDResponse = fInputHandler->GetPIDResponse();
 158
 159   fSelectMask = fInputHandler->IsEventSelected();
 160   if (!(fSelectMask & fSelectBit)) return;
 161   fEventStatistics->Fill("after trigger selection",1);
 162   fRunNumber  = event->GetRunNumber();
 163   fPeriod     = event->GetPeriodNumber();
 164   fOrbit      = event->GetOrbitNumber();
 165   fBc         = event->GetBunchCrossNumber();
 166   fField      = event->GetMagneticField();
 167   fCentrality = event->GetCentrality()->GetCentralityPercentile(fCentralityMethod);
 168
 169   if (fMode==0){ // Data analysis
 170     const AliVVertex* vertex  = event->GetPrimaryVertex();
 171     if (!vertex) return;
 172     fZvtx  = vertex->GetZ();
 173     TString name(vertex->GetName());
 174     if (name.CompareTo("PrimaryVertex") && name.CompareTo("SPDVertex")) return; // Reject TPC only vertex
 175     fEventStatistics->Fill("after rejection of TPC only vertex",1);
 176     if (TMath::Abs(fZvtx) >= fZVertexCut || vertex->GetNContributors()<fnTracksVertex)  return;
 177     fEventStatistics->Fill("after vertex selection",1);
 178
 179     if (fParticles) {
 180       for (Int_t ipart=0;ipart<event->GetNumberOfTracks();ipart++){
 181         AliVTrack* track = (AliVTrack*) event->GetTrack(ipart);
 182         if (!track) continue;
 183         UInt_t mask = GetFilterMap(track);
 184         if (!(mask & fTrackFilterBit)) continue;
 185
 186         if (track->InheritsFrom("AliAODTrack")) AddTrack(track,mask,0);
 187         else if (track->InheritsFrom("AliESDtrack")) {
 188           if (mask)                           AddTrack(track,mask,1);
 189           if (mask & fHybridConstrainedMask)  AddTrack(track,mask,2);
 190           if (mask & fTPConlyConstrainedMask) AddTrack(track,mask,3);
 191         }
 192       }
 193     }
 194     if (fIsAOD){
 195       AliAODEvent* aod = (AliAODEvent*) event;
 196
 197       if (fStoreTracklets){
 198         AliAODTracklets* tracklets = aod->GetTracklets();
 199         Int_t nTracklets = tracklets->GetNumberOfTracklets();
 200         for (Int_t i=0;i<nTracklets;i++){
 201           Float_t phi   = tracklets->GetPhi(i);
 202           Float_t theta = tracklets->GetTheta(i);
 203           Float_t dphi  = tracklets->GetDeltaPhi(i);
 204           if (TMath::Abs(dphi)>fDphiCut) continue;
 205           new ((*fTracklets)[fTracklets->GetEntriesFast()]) AliCFParticle(dphi,-TMath::Log(TMath::Tan(theta/2)),phi,0,0);
 206         }
 207       }
 208
 209       if (fStoreMuons){
 210         for (Int_t iTrack = 0; iTrack < aod->GetNTracks(); iTrack++) {
 211           AliAODTrack* track = aod->GetTrack(iTrack);
 212           if (!track->IsMuonTrack()) continue;
 213           Float_t pt     = track->Pt();
 214           Float_t eta    = track->Eta();
 215           Float_t phi    = track->Phi();
 216           Short_t charge = track->Charge();
 217           Float_t dca    = track->DCA();
 218           Float_t chi2   = track->Chi2perNDF();
 219           Float_t rabs   = track->GetRAtAbsorberEnd();
 220           Int_t   mask   = track->GetMatchTrigger();
 221           if (rabs < 17.6 || rabs > 89.5) continue;
 222           if (eta < -4 || eta > -2.5) continue;
 223           AliCFParticle* part = new ((*fMuons)[fMuons->GetEntriesFast()]) AliCFParticle(pt,eta,phi,charge,mask,3);
 224           part->SetAt(dca,0);
 225           part->SetAt(chi2,1);
 226           part->SetAt(rabs,2);
 227         }
 228       }
 229     }
 230   }
 231
 232   else { // MC analysis
 233     AliMCEvent* mcEvent = fMcHandler ? fMcHandler->MCEvent() : 0;
 234     TClonesArray* mcTracks = dynamic_cast<TClonesArray*>(event->FindListObject(AliAODMCParticle::StdBranchName()));
 235     if (!mcEvent && !mcTracks) { printf("No mc object found\n"); return; }
 236     fEventStatistics->Fill("after mc objects check",1);
 237
 238     Int_t nPrimGen = 0;
 239     Int_t nProduced = 0;
 240     if (mcEvent) { // ESD
 241       AliHeader* header = (AliHeader*) mcEvent->Header();
 242       AliGenCocktailEventHeader* cocktailHeader = dynamic_cast<AliGenCocktailEventHeader*> (header->GenEventHeader());
 243       AliGenEventHeader* mcHeader = dynamic_cast<AliGenEventHeader*> (cocktailHeader ? cocktailHeader->GetHeaders()->First() : header->GenEventHeader());
 244       if (!mcHeader) { printf("mc header not found\n"); return; }
 245       nProduced = mcEvent->GetNumberOfTracks();
 246       nPrimGen = mcHeader->NProduced();
 247       fZvtx = mcEvent->GetPrimaryVertex()->GetZ();
 248     } else { // AOD
 249       AliAODMCHeader* mcHeader = (AliAODMCHeader*) event->FindListObject(AliAODMCHeader::StdBranchName());
 250       if (!mcHeader) { printf("AliAODMCHeader not found\n"); return; }
 251       if (mcHeader->GetCocktailHeaders()) {
 252         nProduced = mcTracks->GetEntriesFast();
 253         AliGenEventHeader* header0 =  mcHeader->GetCocktailHeader(0);
 254         if (!header0) { printf("first header expected but not found\n"); return; }
 255         nPrimGen = header0->NProduced();
 256       } else  nPrimGen = nProduced;
 257       fZvtx = mcHeader->GetVtxZ();
 258     }
 259     fEventStatistics->Fill("after mc header check",1);
 260
 261     if (TMath::Abs(fZvtx)>fZVertexCut) return;
 262     fEventStatistics->Fill("after MC vertex cut",1);
 263
 264     const AliVVertex* vertex = event->GetPrimaryVertex();
 265     if (!vertex) return;
 266     fEventStatistics->Fill("after check for vertex object",1);
 267     TString name(vertex->GetName());
 268     if (name.CompareTo("PrimaryVertex") && name.CompareTo("SPDVertex")) fZvtx=1000; // Reject TPC only vertex
 269     fEventStatistics->Fill("after rejection of TPC only vertex",1);
 270     if (vertex->GetNContributors()<fnTracksVertex)  fZvtx=-1000;
 271     fEventStatistics->Fill("after check for vertex contributors",1);
 272
 273     if (fParticles) {
 274       UInt_t* masks = new UInt_t[nProduced];
 275       for (Int_t i=0;i<nProduced;i++) masks[i]=0;
 276
 277       // Loop over reconstructed tracks to set masks
 278       for (Int_t ipart=0;ipart<event->GetNumberOfTracks();ipart++){
 279         AliVTrack* part = (AliVTrack*) event->GetTrack(ipart);
 280         if (!part) continue;
 281         Int_t label = TMath::Abs(part->GetLabel());
 282         if (label>=nProduced) continue;
 283         masks[label] |= GetFilterMap(part);
 284       }
 285
 286       // Loop over mc tracks to be stored
 287       for (Int_t ipart=0;ipart<nProduced;ipart++){
 288         AliVParticle* part = 0;
 289         Bool_t isPrimary = 0;
 290         if (mcTracks) { // AOD analysis
 291           AliAODMCParticle* particle = (AliAODMCParticle*) mcTracks->At(ipart);
 292           if (!particle) continue;
 293           // skip injected signals
 294           AliAODMCParticle* mother = particle;
 295           while (!mother->IsPrimary()) mother = (AliAODMCParticle*) mcTracks->At(mother->GetMother());
 296           if (mother->GetLabel()>=nPrimGen) continue;
 297           part = particle;
 298           // check for primary
 299           isPrimary = particle->IsPhysicalPrimary();
 300         } else { // ESD analysis
 301           AliMCParticle* particle = (AliMCParticle*) mcEvent->GetTrack(ipart);
 302           if (!particle) continue;
 303           // skip injected signals
 304           AliMCParticle* mother = particle;
 305           while (mother->GetMother()>=0) mother = (AliMCParticle*) mcEvent->GetTrack(mother->GetMother());
 306           if (mother->GetLabel()>=nPrimGen) continue;
 307           part = particle;
 308           // check for primary
 309           isPrimary = mcEvent->IsPhysicalPrimary(ipart);
 310         }
 311
 312         // store only primaries and all reconstructed (non-zero mask)
 313         Int_t mask = masks[ipart];
 314         if (isPrimary) mask |= (1 << 30);
 315         if (!mask) continue;
 316         AddTrack(part,mask);
 317       }
 318       delete[] masks;
 319     }
 320
 321     if (fIsAOD){
 322       AliAODEvent* aod = (AliAODEvent*) event;
 323
 324       if (fTracklets) {
 325         AliAODTracklets* tracklets = aod->GetTracklets();
 326         Int_t nTracklets = tracklets->GetNumberOfTracklets();
 327         for (Int_t i=0;i<nTracklets;i++){
 328           Float_t phi   = tracklets->GetPhi(i);
 329           Float_t theta = tracklets->GetTheta(i);
 330           Float_t dphi  = tracklets->GetDeltaPhi(i);
 331           if (TMath::Abs(dphi)>fDphiCut) continue;
 332           AliCFParticle* tracklet = new ((*fTracklets)[fTracklets->GetEntriesFast()]) AliCFParticle(dphi,-TMath::Log(TMath::Tan(theta/2)),phi,0,0,4);
 333           Int_t label1 = tracklets->GetLabel(i,0);
 334           Int_t label2 = tracklets->GetLabel(i,1);
 335           if (label1!=label2) continue;
 336           if (!mcTracks) continue;
 337           AliAODMCParticle* particle = (AliAODMCParticle*) mcTracks->At(label1);
 338           if (!particle) continue;
 339           Short_t charge = particle->Charge();
 340           Float_t ptMC   = particle->Pt();
 341           Float_t etaMC  = particle->Eta();
 342           Float_t phiMC  = particle->Phi();
 343           Float_t pdg    = particle->PdgCode();
 344           tracklet->SetCharge(charge);
 345           tracklet->SetAt(ptMC,0);
 346           tracklet->SetAt(etaMC,1);
 347           tracklet->SetAt(phiMC,2);
 348           tracklet->SetAt(pdg,3);
 349         }
 350       }
 351
 352       if (fMuons){
 353         for (Int_t iTrack = 0; iTrack < aod->GetNTracks(); iTrack++) {
 354           AliAODTrack* track = aod->GetTrack(iTrack);
 355           if (!track->IsMuonTrack()) continue;
 356           Float_t pt     = track->Pt();
 357           Float_t eta    = track->Eta();
 358           Float_t phi    = track->Phi();
 359           Short_t charge = track->Charge();
 360           Float_t dca    = track->DCA();
 361           Float_t chi2   = track->Chi2perNDF();
 362           Float_t rabs   = track->GetRAtAbsorberEnd();
 363           Int_t   mask   = track->GetMatchTrigger();
 364           if (rabs < 17.6 || rabs > 89.5) continue;
 365           if (eta < -4 || eta > -2.5) continue;
 366           AliCFParticle* part = new ((*fMuons)[fMuons->GetEntriesFast()]) AliCFParticle(pt,eta,phi,charge,mask,3);
 367           part->SetAt(dca,0);
 368           part->SetAt(chi2,1);
 369           part->SetAt(rabs,2);
 370         }
 371       }
 372     }
 373   }
 374   fTree->Fill();
 375   PostData(0,fListOfHistos);
 376   PostData(1,fTree);
 377 }
 378 //-----------------------------------------------------------------------------
 379
 380
 381 //-----------------------------------------------------------------------------
 382 UInt_t AliAnalysisTaskCFTree::GetFilterMap(AliVParticle* particle){
 383   UInt_t mask = 0;
 384   if (particle->InheritsFrom("AliAODTrack")) {
 385     AliAODTrack* part = (AliAODTrack*) particle;
 386     mask = part->GetFilterMap();
 387     Double_t nCrossedRaws      = part->GetTPCNCrossedRows();
 388     Double_t nFindableClusters = part->GetTPCNclsF();
 389     Double_t nSharedClusters   = part->GetTPCnclsS();
 390     Double_t nClusters         = part->GetTPCncls();
 391     Bool_t itsRefit            = part->GetStatus() & AliVTrack::kITSrefit;
 392     if (nCrossedRaws/nFindableClusters > fFoundFractionCut) mask |= (1 << 26);
 393     if (nCrossedRaws>fCrossedRowsCut)                       mask |= (1 << 27);
 394     if (itsRefit)                                           mask |= (1 << 28);
 395     if (nSharedClusters/nClusters<=fSharedClusterCut)       mask |= (1 << 29);
 396     if (part->GetLabel()<0)                                 mask |= (1 << 31);
 397   } else if (particle->InheritsFrom("AliESDtrack")){
 398     AliESDtrack* part = (AliESDtrack*) particle;
 399     if (!fTrackFilter) AliFatal("Track filter undefined");
 400     mask |= fTrackFilter->IsSelected(part);
 401   }
 402
 403   return mask;
 404 }
 405 //-----------------------------------------------------------------------------
 406
 407
 408 //-----------------------------------------------------------------------------
 409 AliCFParticle* AliAnalysisTaskCFTree::AddTrack(AliVParticle* part, UInt_t mask, UInt_t flag){
 410
 411   // skip neutral mc particles
 412   Char_t charge = part->Charge();
 413   if (charge==0) return NULL;
 414
 415   // set pt,eta,phi
 416   Float_t pt=0,eta=0,phi=0;
 417   if (flag==0 || flag==1){ // AOD, MC or Global ESD tracks
 418     pt  = part->Pt();
 419     eta = part->Eta();
 420     phi = part->Phi();
 421     if  (flag==1) mask &= (~fHybridConstrainedMask) & (~fTPConlyConstrainedMask);
 422   }
 423   else if (flag==2) { // Hybrid constrained tracks (ESD)
 424     AliESDtrack* track = (AliESDtrack*) part;
 425     const AliExternalTrackParam* param = track->GetConstrainedParam();
 426     pt  = param->Pt();
 427     eta = param->Eta();
 428     phi = param->Phi();
 429     mask &= fHybridConstrainedMask;
 430   }
 431   else if (flag==3) { // TPC only constrained tracks (ESD)
 432     AliESDtrack* track = (AliESDtrack*) part;
 433     AliESDtrack tpcTrack;
 434     if (!track->FillTPCOnlyTrack(tpcTrack)) return NULL;
 435     AliExternalTrackParam param;
 436     const AliESDVertex* vtxSPD = ((AliESDEvent*) fInputHandler->GetEvent())->GetPrimaryVertexSPD();
 437     if (!tpcTrack.RelateToVertexTPC(vtxSPD,fField,1.e30,&param)) return NULL;
 438     pt  = param.Pt();
 439     eta = param.Eta();
 440     phi = param.Phi();
 441     mask &= fTPConlyConstrainedMask;
 442   }
 443
 444   // kinematic cuts
 445   if (pt < fPtMin || TMath::Abs(eta) > fTrackEtaCut) return NULL;
 446
 447
 448   AliCFParticle* cfpart = new ((*fParticles)[fParticles->GetEntriesFast()]) AliCFParticle(pt,eta,phi,charge,mask,3);
 449
 450   AliVTrack* track = dynamic_cast<AliVTrack*> (part);
 451   if (!track) return cfpart;
 452   Float_t ncl  = track->GetTPCsignalN();
 453   Float_t dedx = track->GetTPCsignalTunedOnData(); if (dedx<=0) dedx = track->GetTPCsignal();
 454   Float_t beta = -99;
 455   if (track->GetStatus()&AliESDtrack::kTOFpid){
 456     Double_t tof[5];
 457     track->GetIntegratedTimes(tof);
 458     beta = tof[0]/track->GetTOFsignal();
 459   }
 460   cfpart->SetAt(ncl,0);
 461   cfpart->SetAt(dedx,1);
 462   cfpart->SetAt(beta,2);
 463   return cfpart;
 464 }