PWG/FLOW/Tasks/AliAnalysisTaskFlowK0Candidates.cxx

   1 /*************************************************************************
   2 * Copyright(c) 1998-2008, 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 /////////////////////////////////////////////////////
  17 // AliAnalysisTaskFlowK0Candidates:
  18 // Analysis task to select K0 candidates for flow analysis.
  19 // Uses one AliESDtrackCuts object for both daughters and
  20 // QA histograms to monitor the reconstruction.
  21 // Author: Carlos Perez (cperez@cern.ch)
  22 //////////////////////////////////////////////////////
  23
  24 #include "TChain.h"
  25 #include "TList.h"
  26 #include "TH1D.h"
  27 #include "TVector3.h"
  28
  29 #include "AliAnalysisTaskSE.h"
  30 #include "AliAnalysisManager.h"
  31
  32 #include "AliESDEvent.h"
  33 #include "AliESDInputHandler.h"
  34 #include "AliESDtrack.h"
  35 #include "AliESDtrackCuts.h"
  36
  37 #include "AliAODEvent.h"
  38 #include "AliAODInputHandler.h"
  39 #include "AliAODTrack.h"
  40
  41 #include "AliCFManager.h"
  42
  43 #include "TMath.h"
  44 #include "TObjArray.h"
  45 #include "AliFlowCandidateTrack.h"
  46 #include "AliFlowEventCuts.h"
  47
  48 #include "AliAnalysisTaskFlowK0Candidates.h"
  49
  50
  51 ClassImp(AliAnalysisTaskFlowK0Candidates)
  52
  53 //_____________________________________________________________________________
  54 AliAnalysisTaskFlowK0Candidates::AliAnalysisTaskFlowK0Candidates() :
  55   AliAnalysisTaskSE(),
  56   fCutsEvent(NULL),
  57   fCuts(NULL),
  58   fQAList(NULL),
  59   fMassMin(0),
  60   fMassMax(0),
  61   fDLcut(0),
  62   fEvent(NULL),
  63   fMulti(NULL)
  64 {
  65   //
  66   for(int i=0; i!=4; ++i) {
  67     fMass[i] = NULL;
  68     fDCA[i]  = NULL;
  69     fDL[i]   = NULL;
  70     fCTP[i]  = NULL;
  71     fd0d0[i] = NULL;
  72     fPhi[i]  = NULL;
  73     fEta[i]  = NULL;
  74     fPt[i]   = NULL;
  75     fAPhi[i] = NULL;
  76     fAEta[i] = NULL;
  77     fAPt[i]  = NULL;
  78     fBPhi[i] = NULL;
  79     fBEta[i] = NULL;
  80     fBPt[i]  = NULL;
  81   }
  82 }
  83
  84 //_____________________________________________________________________________
  85 AliAnalysisTaskFlowK0Candidates::AliAnalysisTaskFlowK0Candidates(const char *name, AliFlowEventCuts *cutsEvent, AliESDtrackCuts *cuts, Double_t MassMin, Double_t MassMax) :
  86   AliAnalysisTaskSE(name),
  87   fCutsEvent(cutsEvent),
  88   fCuts(cuts),
  89   fQAList(NULL),
  90   fMassMin(MassMin),
  91   fMassMax(MassMax),
  92   fDLcut(0),
  93   fEvent(NULL),
  94   fMulti(NULL)
  95 {
  96   //
  97   for(int i=0; i!=4; ++i) {
  98     fMass[i] = NULL;
  99     fDCA[i]  = NULL;
 100     fDL[i]   = NULL;
 101     fCTP[i]  = NULL;
 102     fd0d0[i] = NULL;
 103     fPhi[i]  = NULL;
 104     fEta[i]  = NULL;
 105     fPt[i]   = NULL;
 106     fAPhi[i] = NULL;
 107     fAEta[i] = NULL;
 108     fAPt[i]  = NULL;
 109     fBPhi[i] = NULL;
 110     fBEta[i] = NULL;
 111     fBPt[i]  = NULL;
 112   }
 113   DefineInput(  0, TChain::Class() );
 114   DefineOutput( 1, TObjArray::Class() );
 115   DefineOutput( 2, TList::Class() );
 116 }
 117
 118 //_____________________________________________________________________________
 119 AliAnalysisTaskFlowK0Candidates::~AliAnalysisTaskFlowK0Candidates()
 120 {
 121   //
 122   if(fQAList) delete fQAList;
 123 }
 124
 125 //_____________________________________________________________________________
 126 void AliAnalysisTaskFlowK0Candidates::UserCreateOutputObjects()
 127 {
 128   //
 129   fQAList = new TList();
 130   fQAList->SetOwner();
 131   AddQAEvents();
 132   AddQACandidates();
 133   PostData( 2, fQAList);
 134 }
 135 //_____________________________________________________________________________
 136 void AliAnalysisTaskFlowK0Candidates::AddQAEvents()
 137 {
 138   //
 139   TList *tQAEvents = new TList();
 140   tQAEvents->SetName("Events");
 141   tQAEvents->SetOwner();
 142   fEvent = new TH1D("Event", "Number of Events",    2, 0, 2);
 143   tQAEvents->Add(fEvent);
 144   fMulti = new TH1D("Multiplicity", "Multiplicity", 180, 0, 10000);
 145   tQAEvents->Add(fMulti);
 146   fQAList->Add(tQAEvents);
 147 }
 148 //_____________________________________________________________________________
 149 void AliAnalysisTaskFlowK0Candidates::AddQACandidates()
 150 {
 151   //
 152   TList *tQACandidates[4];
 153   TList *tQADaughters[4];
 154   for(int i=0; i!=4; ++i) {
 155     tQACandidates[i] = new TList();
 156     tQACandidates[i]->SetOwner();
 157     tQACandidates[i]->SetName(Form("Candidates%d",i));
 158     fMass[i] = new TH1D( Form("Mass%i",i), "Mass;M_{#pi#pi} [GeV];Counts per MeV", 180, 0.41, 0.59); tQACandidates[i]->Add( fMass[i] );
 159     fDCA[i]  = new TH1D( Form("DCA%i" ,i), "DCA;[cm];Counts per 10 um",            180, 0.00, 0.18); tQACandidates[i]->Add( fDCA[i]  );
 160     fDL[i]   = new TH1D( Form("DL%i"  ,i), "Decay Length;[cm];Counts per 0.1 mm",  180, 0.00, 1.80); tQACandidates[i]->Add( fDL[i]   );
 161     fCTP[i]  = new TH1D( Form("CTP%i" ,i), "Cos#theta_{p}",                        180,-1.10, 1.10); tQACandidates[i]->Add( fCTP[i]  );
 162     fd0d0[i] = new TH1D( Form("d0d0%i",i), "d_{0}xd_{0};[cm^{2}];Cnts 0.01 mm^{2}",180,-0.009,0.009);tQACandidates[i]->Add( fd0d0[i] );
 163     fPhi[i]  = new TH1D( Form("Phi%i" ,i), "Phi;[rad];Counts per degree",     180,0,TMath::TwoPi()); tQACandidates[i]->Add( fPhi[i]  );
 164     fEta[i]  = new TH1D( Form("Eta%i" ,i), "Eta;;Counts per 0.04",                 180,-3.60, 3.60); tQACandidates[i]->Add( fEta[i]  );
 165     fPt[i]   = new TH1D( Form("Pt%i"  ,i), "Pt;[GeV];Counts per 0.1 GeV",          180, 0.00,18.00); tQACandidates[i]->Add( fPt[i]   );
 166     tQADaughters[i] = new TList();
 167     tQADaughters[i]->SetOwner();
 168     tQADaughters[i]->SetName(Form("Daughters%d",i));
 169     fAPhi[i] = new TH1D( Form("PhiBef%i",i), "Phi prePropagation;[rad];Cnts per degree",180,0,TMath::TwoPi()); tQADaughters[i]->Add( fAPhi[i] );
 170     fAEta[i] = new TH1D( Form("EtaBef%i",i), "Eta prePropagation;;Counts per 0.04"     ,180,-3.6,3.6); tQADaughters[i]->Add( fAEta[i] );
 171     fAPt[i]  = new TH1D( Form("PtBef%i" ,i), "Pt prePropagation;[GeV];Counts per 0.1 GeV",180, 0, 18); tQADaughters[i]->Add( fAPt[i]  );
 172     fBPhi[i] = new TH1D( Form("PhiAft%i",i), "Phi posPropagation;[rad];Cnts per degree",180,0,TMath::TwoPi()); tQADaughters[i]->Add( fBPhi[i] );
 173     fBPhi[i] = new TH1D( Form("EtaAft%i",i), "Eta posPropagation;;Counts per 0.04"     ,180,-3.6,3.6); tQADaughters[i]->Add( fBPhi[i] );
 174     fBPt[i]  = new TH1D( Form("PtAft%i" ,i), "Pt posPropagation;[GeV];Counts per 0.1 GeV",180, 0, 18); tQADaughters[i]->Add( fBPt[i]  );
 175     tQACandidates[i]->Add(tQADaughters[i]);
 176     fQAList->Add(tQACandidates[i]);
 177   }
 178 }
 179 //_____________________________________________________________________________
 180 void AliAnalysisTaskFlowK0Candidates::NotifyRun()
 181 {
 182   //
 183 //  AliESDEvent *fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 184 //  if(!fESD) return;
 185 }
 186
 187 //_____________________________________________________________________________
 188 void AliAnalysisTaskFlowK0Candidates::UserExec(Option_t *)
 189 {
 190   //
 191   AliESDEvent *fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 192   AliAODEvent *fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 193   if( (!fESD)&&(!fAOD) )
 194     return;
 195 //  printf("\nEvent found");
 196   fEvent->Fill( 0 );
 197   if(fCutsEvent)
 198     if( !(fCutsEvent->IsSelected(InputEvent())) )
 199       return;
 200   fEvent->Fill( 1 );
 201   if(fESD)
 202     ReadFromESD(fESD);
 203   else if(fAOD)
 204     ReadFromAOD(fAOD);
 205 }
 206 //_____________________________________________________________________________
 207 void AliAnalysisTaskFlowK0Candidates::ReadFromESD(const AliESDEvent *fESD)
 208 {
 209   //
 210   TObjArray *gPOIselection =  new TObjArray();
 211   gPOIselection->SetOwner();
 212   int nTracks = fESD->GetNumberOfTracks();
 213   fMulti->Fill( nTracks );
 214   for(int i=0; i!=nTracks; ++i) {
 215     AliESDtrack *ioT = fESD->GetTrack(i);
 216     if(fCuts)
 217       if( !(fCuts->IsSelected(ioT)) )
 218         continue;
 219     // printf("\n first particle OK...");
 220     for(int j=i+1; j!=nTracks; ++j) {
 221       AliESDtrack *joT = fESD->GetTrack(j);
 222       if( (ioT->Charge()*joT->Charge()) > 0 )
 223         continue;
 224       if(fCuts)
 225         if( !(fCuts->IsSelected(joT)) )
 226           continue;
 227       // printf("\n second also...");
 228       AliESDtrack *iT = new AliESDtrack(*ioT);
 229       AliESDtrack *jT = new AliESDtrack(*joT);
 230       // getting distance of closest approach
 231       double DCA = iT->PropagateToDCA(jT,fESD->GetMagneticField());
 232       // printf(" propagated...");
 233       // getting decay length
 234       TVector3 vp, vv, vl;
 235       vp = TVector3( (iT->Xv()+jT->Xv())/2, (iT->Yv()+jT->Yv())/2, (iT->Zv()+jT->Zv())/2 ); // candidate position
 236       vv = TVector3( fESD->GetPrimaryVertex()->GetX(), fESD->GetPrimaryVertex()->GetY(), fESD->GetPrimaryVertex()->GetZ() ); // vertex position
 237       vl = vp - vv; // flight line
 238       double gDL;
 239       gDL  = vl.Mag();
 240       // getting cos pointing angle
 241       double gCTP;
 242       TVector3 vi, vj, vs;
 243       vi = TVector3( iT->Px(), iT->Py(), iT->Pz() );
 244       vj = TVector3( jT->Px(), jT->Py(), jT->Pz() );
 245       vs = vi + vj;
 246       gCTP = TMath::Cos( vl.Angle(vs) ); // Marta says: "Carlos, it is faster if you calculate this by hand!"
 247       double gD0D0;
 248       gD0D0 = iT->GetD(vv.X(),vv.Y(),fESD->GetMagneticField())*jT->GetD(vv.X(),vv.Y(),fESD->GetMagneticField());
 249       // getting invariant mass
 250       double sum12 = iT->P()*iT->P()+jT->P()*jT->P();
 251       double pro12 = iT->P()*iT->P()*jT->P()*jT->P();
 252       double gInvMass = 0;
 253       gInvMass += TMath::Power(0.13957018,2);
 254       gInvMass += sqrt( TMath::Power(0.13957018,4) + TMath::Power(0.13957018,2)*(sum12) + pro12 );
 255       gInvMass -= ( iT->Px()*jT->Px()+iT->Py()*jT->Py()+iT->Pz()*jT->Pz() );
 256       gInvMass *= 2;
 257       gInvMass = sqrt(gInvMass);
 258       // filtering
 259       int iLevel = 3;
 260       if( gCTP<0.8 )
 261         iLevel = 2;
 262       if( gDL<fDLcut ) // 0.5
 263         iLevel = 1;
 264       if( DCA>0.05 )
 265         iLevel = 0;
 266       // printf(" candidate at level %d...",iLevel);
 267       for(int h=0; h!=iLevel+1; ++h) {
 268         // printf(" %d",h);
 269         // candidate
 270         fDCA[h]->Fill( DCA );
 271         fDL[h]->Fill( gDL );
 272         fCTP[h]->Fill( gCTP );
 273         fd0d0[h]->Fill( gD0D0 );
 274         fMass[h]->Fill( gInvMass );
 275         fPhi[h]->Fill( vs.Phi()+TMath::Pi() );
 276         fEta[h]->Fill( vs.Eta() );
 277         fPt[h]->Fill( vs.Pt() );
 278         // daughters
 279         fAPhi[h]->Fill( iT->Phi() );
 280         fAEta[h]->Fill( iT->Eta() );
 281         fAPt[h]->Fill(  iT->Pt()  );
 282         fAPhi[h]->Fill( jT->Phi() );
 283         fAEta[h]->Fill( jT->Eta() );
 284         fAPt[h]->Fill(  jT->Pt()  );
 285         fBPhi[h]->Fill( ioT->Phi());
 286         fBPhi[h]->Fill( ioT->Eta());
 287         fBPt[h]->Fill(  ioT->Pt() );
 288         fBPhi[h]->Fill( joT->Phi());
 289         fBPhi[h]->Fill( joT->Eta());
 290         fBPt[h]->Fill(  joT->Pt() );
 291       }
 292       // building candidate
 293       if( (gInvMass>fMassMin) && (gInvMass<fMassMax) && (iLevel==3) ) {
 294         AliFlowCandidateTrack *sTrack = new AliFlowCandidateTrack();
 295         sTrack->SetMass( gInvMass );
 296         sTrack->SetPt( vs.Pt() );
 297         sTrack->SetPhi( vs.Phi()+TMath::Pi() );
 298         sTrack->SetEta( vs.Eta() );
 299         sTrack->SetCharge( 0 );
 300         sTrack->AddDaughter( iT->GetID() );
 301         sTrack->AddDaughter( jT->GetID() );
 302         gPOIselection->AddLast( sTrack );
 303       }
 304       delete iT;
 305       delete jT;
 306     } // endfor j
 307   } // endfor i
 308   PostData( 1, gPOIselection );
 309   PostData( 2, fQAList);
 310 }
 311
 312 //_____________________________________________________________________________
 313 void AliAnalysisTaskFlowK0Candidates::ReadFromAOD(const AliAODEvent *fAOD)
 314 {
 315   //
 316   TObjArray *gPOIselection =  new TObjArray();
 317   gPOIselection->SetOwner();
 318   int nTracks = fAOD->GetNumberOfTracks();
 319   fMulti->Fill( nTracks );
 320   for(int i=0; i!=fAOD->GetNumberOfV0s(); ++i) {
 321     AliAODv0 *iV0 = (AliAODv0*) fAOD->GetV0(i);
 322     if ( iV0->GetOnFlyStatus() ) continue;
 323     if ( (iV0->ChargeProng(0)*iV0->ChargeProng(1))>0 ) continue;
 324     // getting distance of closest approach
 325     double DCA = iV0->DcaV0Daughters();
 326     // getting decay length
 327     double gDL;
 328     double vv[3];
 329     vv[0] = fAOD->GetPrimaryVertex()->GetX();
 330     vv[1] = fAOD->GetPrimaryVertex()->GetY();
 331     vv[2] = fAOD->GetPrimaryVertex()->GetZ();
 332     gDL = iV0->DecayLengthV0(vv);
 333     // cos pointing angle
 334     double gCTP;
 335     gCTP = iV0->CosPointingAngle(vv);
 336     double gD0D0;
 337     gD0D0 = iV0->Prodd0d0();
 338     // getting invariant mass
 339     double gInvMass = iV0->MassK0Short();
 340     // filtering
 341     int iLevel = 3;
 342     if( gCTP<0.8 )
 343       iLevel = 2;
 344     if( gDL<fDLcut ) // 0.5
 345       iLevel = 1;
 346     if( DCA>0.05 )
 347       iLevel = 0;
 348     for(int h=0; h!=iLevel+1; ++h) {
 349       // candidate
 350       fDCA[h]->Fill( DCA );
 351       fDL[h]->Fill( gDL );
 352       fCTP[h]->Fill( gCTP );
 353       fd0d0[h]->Fill( gD0D0 );
 354       fMass[h]->Fill( gInvMass );
 355       fPhi[h]->Fill( iV0->Phi() );
 356       fEta[h]->Fill( iV0->Eta() );
 357       fPt[h]->Fill( iV0->Pt() );
 358       // daughters
 359       fAPhi[h]->Fill( ( (AliAODTrack*) iV0->GetDaughter(0) )->Phi() );
 360       fAEta[h]->Fill( ( (AliAODTrack*) iV0->GetDaughter(0) )->Eta() );
 361       fAPt[h]->Fill(  ( (AliAODTrack*) iV0->GetDaughter(0) )->Pt()  );
 362       fAPhi[h]->Fill( ( (AliAODTrack*) iV0->GetDaughter(1) )->Phi() );
 363       fAEta[h]->Fill( ( (AliAODTrack*) iV0->GetDaughter(1) )->Eta() );
 364       fAPt[h]->Fill(  ( (AliAODTrack*) iV0->GetDaughter(1) )->Pt()  );
 365       fBPhi[h]->Fill( iV0->PhiProng(0) );
 366       fBPhi[h]->Fill( iV0->EtaProng(0) );
 367       fBPt[h]->Fill(  iV0->PtProng(0)  );
 368       fBPhi[h]->Fill( iV0->PhiProng(1) );
 369       fBPhi[h]->Fill( iV0->EtaProng(1) );
 370       fBPt[h]->Fill(  iV0->PtProng(1)  );
 371     }
 372     AliFlowCandidateTrack *sTrack = new AliFlowCandidateTrack();
 373     sTrack->SetMass( gInvMass );
 374     sTrack->SetPt( iV0->Pt() );
 375     sTrack->SetPhi( iV0->Phi() );
 376     sTrack->SetEta( iV0->Eta() );
 377     sTrack->SetCharge( 0 );
 378     sTrack->AddDaughter( iV0->GetPosID() );
 379     sTrack->AddDaughter( iV0->GetNegID() );
 380     gPOIselection->AddLast( sTrack );
 381   }
 382   PostData( 1, gPOIselection );
 383   PostData( 2, fQAList);
 384 }
 385
 386 //_____________________________________________________________________________
 387 void AliAnalysisTaskFlowK0Candidates::Terminate(Option_t *)
 388 {
 389   //
 390 }
 391