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