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 "AliAODTrack.h"
  28 #include "AliESDtrack.h"
  29 #include "AliExternalTrackParam.h"
  30 #include "AliCentrality.h"
  31 #include "AliAnalysisFilter.h"
  32 #include "AliVMultiplicity.h"
  33 #include "AliAnalysisUtils.h"
  34 // root
  35 #include "TMath.h"
  36 #include "TFile.h"
  37 #include "TList.h"
  38 #include "TH1I.h"
  39 #include "TChain.h"
  40 #include "TTree.h"
  41 #include "TClonesArray.h"
  42 ClassImp(AliAnalysisTaskCFTree)
  43
  44 //-----------------------------------------------------------------------------
  45 AliAnalysisTaskCFTree::AliAnalysisTaskCFTree(const char* name) :
  46   AliAnalysisTaskSE(name),
  47   fTrackFilter(0x0),
  48   fHybridConstrainedMask(0),
  49   fTPConlyConstrainedMask(0),
  50   fUtils(0x0),
  51   fListOfHistos(0x0),
  52   fEventStatistics(0x0),
  53   fTree(0x0),
  54   fTracks(0x0),
  55   fTracklets(0x0),
  56   fMuons(0x0),
  57   fMcParticles(0x0),
  58   fField(0),
  59   fCentrality(),
  60   fVtxZ(0),
  61   fVtxTPConly(0),
  62   fVtxContributors(0),
  63   fPeriod(0),
  64   fOrbit(),
  65   fBc(),
  66   fSelectMask(0),
  67   fIsPileupSPD(0),
  68   fIsPileupMV(0),
  69   fSelectBit(AliVEvent::kAny),
  70   fZVertexCut(10.),
  71   fTrackFilterBit(0xffffffff),
  72   fTrackEtaCut(1.0),
  73   fPtMin(0.15),
  74   fSharedClusterCut(0.4),
  75   fCrossedRowsCut(100),
  76   fFoundFractionCut(0.8),
  77   fDphiCut(1.e9),
  78   fStoreTracks(0),
  79   fStoreTracklets(0),
  80   fStoreMuons(0),
  81   fStoreMcTracks(0),
  82   fStoreMcTracklets(0),
  83   fStoreMcMuons(0),
  84   fStorePidInfo(0)
  85 {
  86   Info("AliAnalysisTaskCFTree","Calling Constructor");
  87   DefineInput(0,TChain::Class());
  88   DefineOutput(1,TList::Class());
  89   DefineOutput(2,TTree::Class());
  90 }
  91 //-----------------------------------------------------------------------------
  92
  93
  94 //-----------------------------------------------------------------------------
  95 void AliAnalysisTaskCFTree::UserCreateOutputObjects(){
  96   fListOfHistos = new TList();
  97   fListOfHistos->SetOwner();
  98   fEventStatistics = new TH1I("fEventStatistics","",10,0,10);
  99   fEventStatistics->SetBit(TH1::kCanRebin);
 100
 101   fListOfHistos->Add(fEventStatistics);
 102
 103   if (fStoreTracks)    fTracks      = new TClonesArray("AliCFParticle",2000);
 104   if (fStoreTracklets) fTracklets   = new TClonesArray("AliCFParticle",2000);
 105   if (fStoreMuons)     fMuons       = new TClonesArray("AliCFParticle",2000);
 106   if (fStoreMcTracks)  fMcParticles = new TClonesArray("AliCFParticle",2000);
 107   // create file-resident tree
 108   TDirectory *owd = gDirectory;
 109   OpenFile(1);
 110   fTree = new TTree("events","events");
 111   owd->cd();
 112   fTree->Branch("cent",&fCentrality,"fCentrality[6]/F");
 113   fTree->Branch("vtxz",&fVtxZ);
 114   fTree->Branch("vtxTPConly",&fVtxTPConly);
 115   fTree->Branch("vtxContributors",&fVtxContributors);
 116   fTree->Branch("field",&fField);
 117   fTree->Branch("run",&fCurrentRunNumber);
 118   fTree->Branch("period",&fPeriod);
 119   fTree->Branch("orbit",&fOrbit);
 120   fTree->Branch("bc",&fBc);
 121   fTree->Branch("mask",&fSelectMask);
 122   fTree->Branch("pileupspd",&fIsPileupSPD);
 123   fTree->Branch("pileupmv",&fIsPileupMV);
 124   if (fTracks)      fTree->Branch("tracks",&fTracks);
 125   if (fTracklets)   fTree->Branch("tracklets",&fTracklets);
 126   if (fMuons)       fTree->Branch("muons",&fMuons);
 127   if (fMcParticles) fTree->Branch("mcparticles",&fMcParticles);
 128
 129   fUtils = new AliAnalysisUtils();
 130
 131   PostData(1,fListOfHistos);
 132   PostData(2,fTree);
 133 }
 134 //-----------------------------------------------------------------------------
 135
 136
 137 //-----------------------------------------------------------------------------
 138 void AliAnalysisTaskCFTree::UserExec(Option_t *){
 139   fEventStatistics->Fill("before cuts",1);
 140
 141   if (!fInputEvent) return;
 142   fEventStatistics->Fill("after event check",1);
 143
 144   // TODO
 145   //  TString classes = fInputEvent->GetFiredTriggerClasses();
 146   //  fClassesFired = 0;
 147   //  fClassesFired |= (classes.Contains("CINT7-S-") << 0);
 148   //  fClassesFired |= (classes.Contains("CINT7-B-") << 0);
 149   //  fClassesFired |= (classes.Contains("CSHM8-S-") << 1);
 150   //  fClassesFired |= (classes.Contains("CSHM8-B-") << 1);
 151   //  fClassesFired |= (classes.Contains("CMSL7-S-") << 2);
 152   //  fClassesFired |= (classes.Contains("CMSL7-B-") << 2);
 153   //  fClassesFired |= (classes.Contains("CMSH7-S-") << 3);
 154   //  fClassesFired |= (classes.Contains("CMSH7-B-") << 3);
 155   //  fClassesFired |= (classes.Contains("CMUL7-S-") << 4);
 156   //  fClassesFired |= (classes.Contains("CMUL7-B-") << 4);
 157   //  fClassesFired |= (classes.Contains("CMLL7-S-") << 5);
 158   //  fClassesFired |= (classes.Contains("CMLL7-B-") << 5);
 159   //  fClassesFired |= (classes.Contains("CINT8-S-") << 6);
 160   //  fClassesFired |= (classes.Contains("CSPI7-S-") << 7);
 161   //  fClassesFired |= (classes.Contains("CSPI8-S-") << 8);
 162   //  fClassesFired |= (classes.Contains("CMSL8-S-") << 9);
 163   //  fClassesFired |= (classes.Contains("CMSH8-S-") <<10);
 164   //  fClassesFired |= (classes.Contains("CMUL8-S-") <<11);
 165   //  fClassesFired |= (classes.Contains("CMLL8-S-") <<12);
 166   //  fClassesFired |= (classes.Contains("C0MUL-SA") <<13);
 167   //  fClassesFired |= (classes.Contains("C0MUL-SC") <<14);
 168   //  if (!fClassesFired) return;
 169   //  fEventStatistics->Fill("after trigger check",1);
 170
 171   fSelectMask = fInputHandler->IsEventSelected();
 172   if (!(fSelectMask & fSelectBit)) return;
 173
 174   fEventStatistics->Fill("after physics selection",1);
 175
 176   fPeriod        = fInputEvent->GetPeriodNumber();
 177   fOrbit         = fInputEvent->GetOrbitNumber();
 178   fBc            = fInputEvent->GetBunchCrossNumber();
 179   fField         = fInputEvent->GetMagneticField();
 180   fCentrality[0] = fInputEvent->GetCentrality()->GetCentralityPercentile("V0M");
 181   fCentrality[1] = fInputEvent->GetCentrality()->GetCentralityPercentile("V0A");
 182   fCentrality[2] = fInputEvent->GetCentrality()->GetCentralityPercentile("V0C");
 183   fCentrality[3] = fInputEvent->GetCentrality()->GetCentralityPercentile("CL1");
 184   fCentrality[4] = fInputEvent->GetCentrality()->GetCentralityPercentile("ZNA");
 185   fCentrality[5] = fInputEvent->GetCentrality()->GetCentralityPercentile("ZNC");
 186   fIsPileupSPD   = fUtils->IsPileUpSPD(fInputEvent);
 187   fIsPileupMV    = fUtils->IsPileUpMV(fInputEvent);
 188 //  fNofITSClusters[i] = esd->GetMultiplicity()->GetNumberOfITSClusters(i);
 189
 190   const AliVVertex* vertex  = fInputEvent->GetPrimaryVertex();
 191   fVtxZ  = vertex->GetZ();
 192   fVtxTPConly = TString(vertex->GetName()).CompareTo("PrimaryVertex") && TString(vertex->GetName()).CompareTo("SPDVertex");
 193   fVtxContributors = vertex->GetNContributors();
 194   if (TMath::Abs(fVtxZ) >= fZVertexCut)  return;
 195   fEventStatistics->Fill("after vertex cut",1);
 196
 197   if (fTracks) {
 198     fTracks->Clear();
 199     for (Int_t ipart=0;ipart<fInputEvent->GetNumberOfTracks();ipart++){
 200       AliVTrack* track = (AliVTrack*) fInputEvent->GetTrack(ipart);
 201       if (!track) continue;
 202       UInt_t mask = GetFilterMap(track);
 203       if (!(mask & fTrackFilterBit)) continue;
 204
 205       if (track->InheritsFrom("AliAODTrack")) AddTrack(track,mask,0);
 206       else if (track->InheritsFrom("AliESDtrack")) {
 207         if (mask)                           AddTrack(track,mask,1);
 208         if (mask & fHybridConstrainedMask)  AddTrack(track,mask,2);
 209         if (mask & fTPConlyConstrainedMask) AddTrack(track,mask,3);
 210       }
 211     }
 212   }
 213
 214   if (fTracklets){
 215     fTracklets->Clear();
 216     AliVMultiplicity* mult = fInputEvent->GetMultiplicity();
 217     Int_t nTracklets = mult->GetNumberOfTracklets();
 218     for (Int_t i=0;i<nTracklets;i++){
 219       Float_t phi   = mult->GetPhi(i);
 220       Float_t eta   = -TMath::Log(TMath::Tan(mult->GetTheta(i)/2));
 221       Float_t dphi  = mult->GetDeltaPhi(i);
 222       if (TMath::Abs(dphi)>fDphiCut) continue;
 223       AliCFParticle* tracklet = new ((*fTracklets)[fTracklets->GetEntriesFast()]) AliCFParticle(dphi,eta,phi,0,0,fStoreMcTracklets?4:0);
 224       if (!fStoreMcTracklets || !fMCEvent) continue;
 225       Int_t label1 = mult->GetLabel(i,0);
 226       Int_t label2 = mult->GetLabel(i,1);
 227       if (label1!=label2) continue;
 228       AliVParticle* particle = fMCEvent->GetTrack(label1);
 229       if (!particle) continue;
 230       Short_t charge = particle->Charge();
 231       Float_t ptMC   = particle->Pt();
 232       Float_t etaMC  = particle->Eta();
 233       Float_t phiMC  = particle->Phi();
 234       Float_t pdg    = particle->PdgCode();
 235       tracklet->SetCharge(charge);
 236       tracklet->SetAt(ptMC,0);
 237       tracklet->SetAt(etaMC,1);
 238       tracklet->SetAt(phiMC,2);
 239       tracklet->SetAt(pdg,3);
 240     }
 241   }
 242
 243   AliAODEvent* aod = dynamic_cast<AliAODEvent*> (fInputEvent);
 244   if (fMuons && aod){ // aod only
 245     fMuons->Clear();
 246     for (Int_t iTrack = 0; iTrack < aod->GetNumberOfTracks(); iTrack++) {
 247       AliAODTrack* track = dynamic_cast<AliAODTrack*>(aod->GetTrack(iTrack));
 248       if(!track) {
 249         AliWarning("Not a standard AOD");
 250         continue;
 251       }
 252       if (!track->IsMuonTrack()) continue;
 253       Float_t pt     = track->Pt();
 254       Float_t eta    = track->Eta();
 255       Float_t phi    = track->Phi();
 256       Short_t charge = track->Charge();
 257       Float_t dca    = track->DCA();
 258       Float_t chi2   = track->Chi2perNDF();
 259       Float_t rabs   = track->GetRAtAbsorberEnd();
 260       Int_t   mask   = track->GetMatchTrigger();
 261       if (rabs < 17.6 || rabs > 89.5) continue;
 262       if (eta < -4 || eta > -2.5) continue;
 263       AliCFParticle* part = new ((*fMuons)[fMuons->GetEntriesFast()]) AliCFParticle(pt,eta,phi,charge,mask,fStoreMcMuons?11:3);
 264       part->SetAt(dca,0);
 265       part->SetAt(chi2,1);
 266       part->SetAt(rabs,2);
 267       if (!fStoreMcMuons || !fMCEvent) continue;
 268       Int_t label = TMath::Abs(track->GetLabel());
 269       AliVParticle* mcpart = fMCEvent->GetTrack(label);
 270       if (!mcpart) continue;
 271       Int_t mcpdg = mcpart->PdgCode();
 272       Float_t mcpt  = mcpart->Pt();
 273       Float_t mceta = mcpart->Eta();
 274       Float_t mcphi = mcpart->Phi();
 275       part->SetAt(mcpt,3);
 276       part->SetAt(mceta,4);
 277       part->SetAt(mcphi,5);
 278       part->SetAt(mcpdg,6);
 279       part->SetAt(0,10);
 280       Bool_t isPrimary = fMCEvent->IsPhysicalPrimary(label);
 281       if (isPrimary) continue;
 282       label = mcpart->GetMother();
 283       while (!isPrimary && label>=0) {
 284         mcpart = (AliVParticle*) fMCEvent->GetTrack(label);
 285         if (!mcpart) continue;
 286         label = mcpart->GetMother();
 287         isPrimary = fMCEvent->IsPhysicalPrimary(label);
 288       }
 289       if (!mcpart) continue;
 290       Float_t mcprimarypt  = mcpart->Pt();
 291       Float_t mcprimaryeta = mcpart->Eta();
 292       Float_t mcprimaryphi = mcpart->Phi();
 293       Int_t   mcprimarypdg = mcpart->PdgCode();
 294       part->SetAt(mcprimarypt,7);
 295       part->SetAt(mcprimaryeta,8);
 296       part->SetAt(mcprimaryphi,9);
 297       part->SetAt(mcprimarypdg,10);
 298     }
 299   }
 300
 301   if (fMcParticles && fMCEvent) {
 302     fMcParticles->Clear();
 303     TString gen;
 304     for (Int_t ipart=0;ipart<fMCEvent->GetNumberOfTracks();ipart++){
 305       AliVParticle* part = fMCEvent->GetTrack(ipart);
 306       Bool_t isCocktail = fMCEvent->GetCocktailGenerator(ipart,gen);
 307       if (isCocktail && !gen.Contains("Pythia") && !gen.Contains("Hijing") && !gen.Contains("AMPT") && !gen.Contains("DPMJET")) continue;
 308       Float_t pt     = part->Pt();
 309       Float_t eta    = part->Eta();
 310       Float_t phi    = part->Phi();
 311       Char_t  charge = part->Charge();
 312       Int_t   mask   = part->PdgCode();
 313       // TODO
 314       //  Bool_t isPrimary = fMCEvent->IsPhysicalPrimary(ipart);
 315       // if (pt < fMcPtMin) continue;
 316       // if (TMath::Abs(eta) > fMcTrackEtaCut) continue;
 317       new ((*fMcParticles)[fMcParticles->GetEntriesFast()]) AliCFParticle(pt,eta,phi,charge,mask);
 318     }
 319   }
 320
 321   fTree->Fill();
 322   PostData(1,fListOfHistos);
 323   PostData(2,fTree);
 324 }
 325 //-----------------------------------------------------------------------------
 326
 327
 328 //-----------------------------------------------------------------------------
 329 UInt_t AliAnalysisTaskCFTree::GetFilterMap(AliVTrack* track){
 330   UInt_t mask = 0;
 331   if (track->InheritsFrom("AliAODTrack")) {
 332     AliAODTrack* part = (AliAODTrack*) track;
 333     mask = part->GetFilterMap();
 334     Double_t nCrossedRaws      = part->GetTPCNCrossedRows();
 335     Double_t nFindableClusters = part->GetTPCNclsF();
 336     Double_t nSharedClusters   = part->GetTPCnclsS();
 337     Double_t nClusters         = part->GetTPCncls();
 338     Bool_t itsRefit            = part->GetStatus() & AliVTrack::kITSrefit;
 339     if (nCrossedRaws/nFindableClusters > fFoundFractionCut) mask |= (1 << 26);
 340     if (nCrossedRaws>fCrossedRowsCut)                       mask |= (1 << 27);
 341     if (itsRefit)                                           mask |= (1 << 28);
 342     if (nSharedClusters/nClusters<=fSharedClusterCut)       mask |= (1 << 29);
 343     if (part->GetLabel()<0)                                 mask |= (1 << 31);
 344   } else if (track->InheritsFrom("AliESDtrack")){
 345     AliESDtrack* part = (AliESDtrack*) track;
 346     if (!fTrackFilter) AliFatal("Track filter undefined");
 347     mask |= fTrackFilter->IsSelected(part);
 348   }
 349
 350   return mask;
 351 }
 352 //-----------------------------------------------------------------------------
 353
 354
 355 //-----------------------------------------------------------------------------
 356 AliCFParticle* AliAnalysisTaskCFTree::AddTrack(AliVTrack* track, UInt_t mask, UInt_t flag){
 357
 358   // skip neutral mc trackicles
 359   Char_t charge = track->Charge();
 360   if (charge==0) return NULL;
 361
 362   // set pt,eta,phi
 363   Float_t pt=0,eta=0,phi=0;
 364   if (flag==0 || flag==1){ // AOD or Global ESD tracks
 365     pt  = track->Pt();
 366     eta = track->Eta();
 367     phi = track->Phi();
 368     if  (flag==1) mask &= (~fHybridConstrainedMask) & (~fTPConlyConstrainedMask);
 369   }
 370   else if (flag==2) { // Hybrid constrained tracks (ESD)
 371     AliESDtrack* part = (AliESDtrack*) track;
 372     const AliExternalTrackParam* param = part->GetConstrainedParam();
 373     pt  = param->Pt();
 374     eta = param->Eta();
 375     phi = param->Phi();
 376     mask &= fHybridConstrainedMask;
 377   }
 378   else if (flag==3) { // TPC only constrained tracks (ESD)
 379     AliESDtrack* part = (AliESDtrack*) track;
 380     AliESDtrack tpcTrack;
 381     if (!part->FillTPCOnlyTrack(tpcTrack)) return NULL;
 382     AliExternalTrackParam param;
 383     const AliESDVertex* vtxSPD = ((AliESDEvent*) fInputEvent)->GetPrimaryVertexSPD();
 384     if (!tpcTrack.RelateToVertexTPC(vtxSPD,fField,1.e30,&param)) return NULL;
 385     pt  = param.Pt();
 386     eta = param.Eta();
 387     phi = param.Phi();
 388     mask &= fTPConlyConstrainedMask;
 389   }
 390
 391   // kinematic cuts
 392   if (pt < fPtMin || TMath::Abs(eta) > fTrackEtaCut) return NULL;
 393
 394   AliCFParticle* cftrack = new ((*fTracks)[fTracks->GetEntriesFast()]) AliCFParticle(pt,eta,phi,charge,mask,fStorePidInfo ? 3: 0);
 395
 396   if (!fStorePidInfo) return cftrack;
 397   Float_t ncl  = track->GetTPCsignalN();
 398   Float_t dedx = track->GetTPCsignalTunedOnData(); if (dedx<=0) dedx = track->GetTPCsignal();
 399   Float_t beta = -99;
 400   if (track->GetStatus()&AliESDtrack::kTOFpid){
 401     Double_t tof[5];
 402     track->GetIntegratedTimes(tof);
 403     beta = tof[0]/track->GetTOFsignal();
 404   }
 405   cftrack->SetAt(ncl,0);
 406   cftrack->SetAt(dedx,1);
 407   cftrack->SetAt(beta,2);
 408   return cftrack;
 409 }