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