CORRFW/test/AliCFV0Task.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 //-----------------------------------------------------------------------
  17 // Example of task running on AliEn (CAF?)
  18 // which provides standard way of calculating acceptance and efficiency
  19 // between different steps of the procedure.
  20 // The ouptut of the task is a AliCFContainer from which the efficiencies
  21 // can be calculated
  22 //-----------------------------------------------------------------------
  23 // Author : R. Vernet, Consorzio Cometa - Catania (it)
  24 //-----------------------------------------------------------------------
  25
  26
  27 #ifndef ALICFV0TASK_CXX
  28 #define ALICFV0TASK_CXX
  29
  30 #include "AliCFV0Task.h"
  31 #include "TCanvas.h"
  32 #include "AliStack.h"
  33 #include "TParticle.h"
  34 #include "TH1I.h"
  35 #include "AliMCEvent.h"
  36 #include "AliCFManager.h"
  37 #include "AliCFContainer.h"
  38 #include "AliESDtrack.h"
  39 #include "AliLog.h"
  40 #include "AliESDv0.h"
  41 #include "AliV0vertexer.h"
  42 #include "AliCFPair.h"
  43 #include "AliESDEvent.h"
  44 #include "AliAODEvent.h"
  45 #include "TChain.h"
  46 #include "AliCFParticleGenCuts.h"
  47 #include "AliAODv0.h"
  48 #include "TDatabasePDG.h"
  49
  50 //__________________________________________________________________________
  51 AliCFV0Task::AliCFV0Task() :
  52   AliAnalysisTaskSE(),
  53   fRebuildV0s(0),
  54   fV0PDG(310),
  55   fCFManager(0x0),
  56   fHistEventsProcessed(0x0)
  57 {
  58   //
  59   //Default ctor
  60   //
  61 }
  62 //___________________________________________________________________________
  63 AliCFV0Task::AliCFV0Task(const Char_t* name) :
  64   AliAnalysisTaskSE(name),
  65   fRebuildV0s(0),
  66   fV0PDG(310),
  67   fCFManager(0x0),
  68   fHistEventsProcessed(0x0)
  69 {
  70   //
  71   // Constructor. Initialization of Inputs and Outputs
  72   //
  73   Info("AliCFV0Task","Calling Constructor");
  74   /*
  75     DefineInput(0) and DefineOutput(0)
  76     are taken care of by AliAnalysisTaskSE constructor
  77   */
  78   DefineOutput(1,TH1I::Class());
  79   DefineOutput(2,AliCFContainer::Class());
  80 //   DefineOutput(3,TList::Class());
  81 }
  82
  83 //___________________________________________________________________________
  84 AliCFV0Task& AliCFV0Task::operator=(const AliCFV0Task& c)
  85 {
  86   //
  87   // Assignment operator
  88   //
  89   if (this!=&c) {
  90     AliAnalysisTaskSE::operator=(c) ;
  91     fRebuildV0s = c.fRebuildV0s;
  92     fV0PDG      = c.fV0PDG;
  93     fCFManager  = c.fCFManager;
  94     fHistEventsProcessed = c.fHistEventsProcessed;
  95   }
  96   return *this;
  97 }
  98
  99 //___________________________________________________________________________
 100 AliCFV0Task::AliCFV0Task(const AliCFV0Task& c) :
 101   AliAnalysisTaskSE(c),
 102   fRebuildV0s(c.fRebuildV0s),
 103   fV0PDG(c.fV0PDG),
 104   fCFManager(c.fCFManager),
 105   fHistEventsProcessed(c.fHistEventsProcessed)
 106 {
 107   //
 108   // Copy Constructor
 109   //
 110 }
 111
 112 //___________________________________________________________________________
 113 AliCFV0Task::~AliCFV0Task() {
 114   //
 115   //destructor
 116   //
 117   Info("~AliCFV0Task","Calling Destructor");
 118   if (fCFManager)           delete fCFManager ;
 119   if (fHistEventsProcessed) delete fHistEventsProcessed ;
 120 }
 121
 122 //_________________________________________________
 123 void AliCFV0Task::UserExec(Option_t *)
 124 {
 125   //
 126   // Main loop function
 127   //
 128   Info("UserExec","") ;
 129
 130   if (!fInputEvent) {
 131     Error("UserExec","NO EVENT FOUND!");
 132     return;
 133   }
 134
 135   if (!fMCEvent) {
 136     Error("UserExec","NO MC INFO FOUND!");
 137     return;
 138   }
 139
 140   fCFManager->SetMCEventInfo (fMCEvent);
 141   fCFManager->SetRecEventInfo(fInputEvent);
 142
 143   Bool_t isESDEvent = strcmp(fInputEvent->ClassName(),"AliESDEvent") == 0 ? kTRUE : kFALSE ;
 144   Bool_t isAODEvent = strcmp(fInputEvent->ClassName(),"AliAODEvent") == 0 ? kTRUE : kFALSE ;
 145
 146   AliESDEvent* fESD = dynamic_cast<AliESDEvent*>(fInputEvent);
 147   AliAODEvent* fAOD = dynamic_cast<AliAODEvent*>(fInputEvent);
 148
 149
 150   // MC-event selection
 151   Double_t containerInput[2] ;
 152
 153   //loop on the MC event
 154   Info("UserExec","Looping on MC event");
 155   for (Int_t ipart=0; ipart<fMCEvent->GetNumberOfTracks(); ipart++) {
 156     AliMCParticle *mcPart  = (AliMCParticle*)fMCEvent->GetTrack(ipart);
 157
 158     //check the MC-level cuts
 159     if (!fCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,mcPart)) continue;
 160     containerInput[0] = mcPart->Pt();
 161     containerInput[1] = mcPart->Y() ;
 162     //fill the container for Gen-level selection
 163     fCFManager->GetParticleContainer()->Fill(containerInput,kStepGenerated);
 164
 165     //check the Acceptance-level cuts
 166     if (!fCFManager->CheckParticleCuts(AliCFManager::kPartAccCuts,mcPart)) continue;
 167     //fill the container for Acceptance-level selection
 168     fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructible);
 169   }
 170
 171
 172   //Now go to rec level
 173   Info("UserExec","Looping on %s",fInputEvent->ClassName());
 174
 175   if (isESDEvent && fRebuildV0s) RebuildV0s(fESD) ;
 176
 177   Info("UserExec","There are %d V0s in event",fInputEvent->GetNumberOfV0s());
 178
 179   AliESDv0*  esdV0 = 0x0;
 180   AliAODv0*  aodV0 = 0x0;
 181   AliCFPair* pair  = 0x0;
 182   Int_t      labMCV0 = 0;
 183
 184   for (Int_t iV0 = 0; iV0<fInputEvent->GetNumberOfV0s(); iV0++) {
 185
 186     if (isESDEvent) {
 187       esdV0 = fESD->GetV0(iV0);
 188       //check if mother reconstructed V0 can be associated to a MC V0
 189       labMCV0 = IsMcV0(esdV0);
 190       if (labMCV0 <0) continue;
 191       pair = new AliCFPair(esdV0,fESD);
 192     }
 193     else if (isAODEvent) {
 194       aodV0 = fAOD->GetV0(iV0);
 195       labMCV0 = IsMcV0(aodV0);
 196       if (labMCV0 <0) continue;
 197       pair = new AliCFPair(aodV0);
 198     }
 199     else {
 200       Error("UserExec","Error: input data file is not an ESD nor an AOD");
 201       return;
 202     }
 203
 204     pair->SetV0PDG(fV0PDG);
 205     pair->SetLabel(labMCV0);
 206     if (!fCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pair)) continue;
 207
 208     //check if associated MC v0 passes the cuts
 209     AliMCParticle* mcV0 = (AliMCParticle*)fMCEvent->GetTrack(labMCV0);
 210     if (!mcV0) continue;
 211
 212     if (!fCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,mcV0)) continue;
 213
 214     //fill the container
 215     Double_t pt, rapidity;
 216
 217     if (isESDEvent) {
 218       Double_t mom[3];
 219       esdV0->GetPxPyPz(mom[0],mom[1],mom[2]);
 220       Double_t pt2 = mom[0]*mom[0]+mom[1]*mom[1] ;
 221       pt  = TMath::Sqrt(pt2);
 222       Double_t energy  = TMath::Sqrt(pt2 + mom[2]*mom[2] + TMath::Power(TDatabasePDG::Instance()->GetParticle(fV0PDG)->Mass(),2));
 223       rapidity = GetRapidity(energy,mom[2]) ;
 224     }
 225     else {
 226       pt = aodV0->Pt();
 227       rapidity = aodV0->Y(fV0PDG);
 228     }
 229
 230     containerInput[0] = pt ;
 231     containerInput[1] = rapidity ;
 232     fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed) ;
 233     if (!fCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,pair)) continue ;
 234     fCFManager->GetParticleContainer()->Fill(containerInput,kStepSelected);
 235
 236     delete pair;
 237   }
 238
 239   fHistEventsProcessed->Fill(0);
 240   /* PostData(0) is taken care of by AliAnalysisTaskSE */
 241   PostData(1,fHistEventsProcessed) ;
 242   PostData(2,fCFManager->GetParticleContainer()) ;
 243 }
 244
 245
 246 //___________________________________________________________________________
 247 void AliCFV0Task::Terminate(Option_t*)
 248 {
 249   // The Terminate() function is the last function to be called during
 250   // a query. It always runs on the client, it can be used to present
 251   // the results graphically or save the results to file.
 252
 253   Info("Terminate","");
 254   AliAnalysisTaskSE::Terminate();
 255
 256
 257   //draw some example plots....
 258
 259   AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
 260
 261   TH1D* h00 =   cont->ShowProjection(0,0) ;
 262   TH1D* h01 =   cont->ShowProjection(0,1) ;
 263   TH1D* h02 =   cont->ShowProjection(0,2) ;
 264   TH1D* h03 =   cont->ShowProjection(0,3) ;
 265
 266   TH1D* h10 =   cont->ShowProjection(1,0) ;
 267   TH1D* h11 =   cont->ShowProjection(1,1) ;
 268   TH1D* h12 =   cont->ShowProjection(1,2) ;
 269   TH1D* h13 =   cont->ShowProjection(1,3) ;
 270
 271   Double_t max1 = h00->GetMaximum();
 272   Double_t max2 = h10->GetMaximum();
 273
 274   h00->GetYaxis()->SetRangeUser(0,max1*1.2);
 275   h01->GetYaxis()->SetRangeUser(0,max1*1.2);
 276   h02->GetYaxis()->SetRangeUser(0,max1*1.2);
 277   h03->GetYaxis()->SetRangeUser(0,max1*1.2);
 278
 279   h10->GetYaxis()->SetRangeUser(0,max2*1.2);
 280   h11->GetYaxis()->SetRangeUser(0,max2*1.2);
 281   h12->GetYaxis()->SetRangeUser(0,max2*1.2);
 282   h13->GetYaxis()->SetRangeUser(0,max2*1.2);
 283
 284   h00->SetMarkerStyle(23) ;
 285   h01->SetMarkerStyle(24) ;
 286   h02->SetMarkerStyle(25) ;
 287   h03->SetMarkerStyle(26) ;
 288
 289   h10->SetMarkerStyle(23) ;
 290   h11->SetMarkerStyle(24) ;
 291   h12->SetMarkerStyle(25) ;
 292   h13->SetMarkerStyle(26) ;
 293
 294   TCanvas * c =new TCanvas("c","",1400,800);
 295   c->Divide(4,2);
 296
 297   c->cd(1);
 298   h00->Draw("p");
 299   c->cd(2);
 300   h01->Draw("p");
 301   c->cd(3);
 302   h02->Draw("p");
 303   c->cd(4);
 304   h03->Draw("p");
 305   c->cd(5);
 306   h10->Draw("p");
 307   c->cd(6);
 308   h11->Draw("p");
 309   c->cd(7);
 310   h12->Draw("p");
 311   c->cd(8);
 312   h13->Draw("p");
 313
 314   c->SaveAs("plots.eps");
 315 }
 316
 317 //___________________________________________________________________________
 318 void AliCFV0Task::UserCreateOutputObjects() {
 319   //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
 320   //TO BE SET BEFORE THE EXECUTION OF THE TASK
 321   //
 322   Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
 323
 324   //slot #1
 325   OpenFile(1);
 326   fHistEventsProcessed = new TH1I("fHistEventsProcessed","",1,0,1) ;
 327 }
 328
 329 //___________________________________________________________________________
 330 Int_t AliCFV0Task::GetV0Label(UInt_t lab1, UInt_t lab2) const {
 331   //
 332   // returns the label of the V0, given the labels of the 2 daughter tracks
 333   // returns -1 if the V0 is fake
 334   //
 335
 336   AliStack* stack = fMCEvent->Stack();
 337   TParticle* part1 = stack->Particle(lab1) ;
 338   TParticle* part2 = stack->Particle(lab2) ;
 339
 340   Int_t part1MotherLab=part1->GetFirstMother();
 341   Int_t part2MotherLab=part2->GetFirstMother();
 342
 343   if (part1MotherLab==-1 || part2MotherLab==-1) return -1 ;
 344   if (part1MotherLab != part2MotherLab )        return -1 ;
 345
 346   if (stack->Particle(part1MotherLab)->GetPdgCode() != fV0PDG ) return -1 ;
 347
 348   switch (fV0PDG) {
 349   case kK0Short :
 350     if ( (part1->GetPdgCode()==  -211 && part2->GetPdgCode()==  211)  ||
 351          (part1->GetPdgCode()==   211 && part2->GetPdgCode()== -211)   )   return part1MotherLab ;
 352     break ;
 353   case kLambda0 :
 354     if ( (part1->GetPdgCode()==  -211 && part2->GetPdgCode()== 2212)  ||
 355          (part1->GetPdgCode()==  2212 && part2->GetPdgCode()== -211)   )   return part1MotherLab ;
 356     break ;
 357   case kLambda0Bar :
 358     if ( (part1->GetPdgCode()==   211 && part2->GetPdgCode()== -2212) ||
 359          (part1->GetPdgCode()== -2212 && part2->GetPdgCode()==   211)  )   return part1MotherLab ;
 360     break ;
 361   default :
 362     return -1;
 363     break ;
 364   }
 365
 366   return -1 ;
 367 }
 368
 369 //___________________________________________________________________________
 370 Double_t AliCFV0Task::GetRapidity(Double_t energy, Double_t pz) {
 371   //
 372   // calculates rapidity, checking energy is larger than pz, otherwise returns 999.
 373   //
 374
 375   if (energy == pz || energy == -pz) {
 376     printf("GetRapidity : ERROR : rapidity for 4-vector with E = Pz -- infinite result");
 377     return 999;
 378   }
 379   if (energy < pz) {
 380     printf("GetRapidity : ERROR : rapidity for 4-vector with E = Pz -- infinite result");
 381     return 999;
 382   }
 383   Double_t y = 0.5 * log((energy + pz) / (energy - pz));
 384   return y;
 385 }
 386
 387 //___________________________________________________________________________
 388 Int_t AliCFV0Task::IsMcV0(AliESDv0* v0) const {
 389   //
 390   // check if the passed V0 is associated to a MC one,
 391   //   and returns the corresponding geant label.
 392   // returns -1 if the V0 is fake (i.e. label<0).
 393   //
 394
 395   AliESDEvent* fESD = dynamic_cast<AliESDEvent*>(fInputEvent);
 396
 397   Int_t nindex = v0->GetNindex();
 398   Int_t pindex = v0->GetPindex();
 399   AliESDtrack *nTrack = fESD->GetTrack(nindex) ;
 400   AliESDtrack *pTrack = fESD->GetTrack(pindex) ;
 401
 402   if (!nTrack || !pTrack) return -1 ;
 403
 404   Int_t nlab  = nTrack->GetLabel() ;
 405   Int_t plab  = pTrack->GetLabel() ;
 406
 407   if (nlab <0 || plab <0) return -1 ;
 408
 409   return GetV0Label((UInt_t)nlab,(UInt_t)plab) ;
 410 }
 411
 412 //___________________________________________________________________________
 413 Int_t AliCFV0Task::IsMcV0(AliAODv0* v0) const {
 414   //
 415   // check if the passed V0 is associated to a MC one,
 416   //   and returns the corresponding geant label.
 417   // returns -1 if the V0 is fake (i.e. label<0).
 418   //
 419
 420   AliAODVertex * vtx = v0->GetSecondaryVtx();
 421
 422   AliAODTrack *nTrack = (AliAODTrack*)vtx->GetDaughter(1); //neg is filled after pos in AliAnalysisTaskStrange
 423   AliAODTrack *pTrack = (AliAODTrack*)vtx->GetDaughter(0);
 424
 425   if (!nTrack || !pTrack) return -1 ;
 426
 427   Int_t nlab  = nTrack->GetLabel() ;
 428   Int_t plab  = pTrack->GetLabel() ;
 429
 430   if (nlab <0 || plab <0) return -1 ;
 431
 432   return GetV0Label((UInt_t)nlab,(UInt_t)plab) ;
 433 }
 434
 435
 436 //___________________________________________________________________________
 437 void AliCFV0Task::RebuildV0s(AliESDEvent* fESD) {
 438
 439   fESD->ResetV0s();
 440
 441   //These are pp cuts : to change if Pb+Pb !
 442   Double_t cuts[]={33,  // max. allowed chi2
 443                    0.01,// min. allowed negative daughter's impact parameter
 444                    0.01,// min. allowed positive daughter's impact parameter
 445                    0.5,// max. allowed DCA between the daughter tracks
 446                    0.98,// max. allowed cosine of V0's pointing angle
 447                    0.2,  // min. radius of the fiducial volume
 448                    100.   // max. radius of the fiducial volume
 449   };
 450   //   AliV0vertexer* v0Vertexer = new AliV0vertexer(cuts);
 451   //   v0Vertexer->SetVertex(primVertexPos);
 452   AliV0vertexer* v0Vertexer = new AliV0vertexer();
 453   v0Vertexer->SetCuts(cuts) ;
 454   v0Vertexer->Tracks2V0vertices(fESD);
 455   delete v0Vertexer ;
 456 }
 457
 458 #endif