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