CORRFW/test/AliCFRsnTask.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 // Author : R. Vernet, Consorzio Cometa - Catania (it)
  18 //-----------------------------------------------------------------------
  19
  20
  21 #include "AliCFRsnTask.h"
  22 #include "TCanvas.h"
  23 #include "AliStack.h"
  24 #include "TParticle.h"
  25 #include "TH1I.h"
  26 #include "TChain.h"
  27 #include "AliMCEvent.h"
  28 #include "AliESDEvent.h"
  29 #include "AliCFManager.h"
  30 #include "AliCFCutBase.h"
  31 #include "AliCFContainer.h"
  32 #include "AliESDtrack.h"
  33 #include "AliLog.h"
  34 #include "AliRsnDaughter.h"
  35 #include "AliCFPair.h"
  36 #include "AliRsnParticle.h"
  37
  38 //__________________________________________________________________________
  39 AliCFRsnTask::AliCFRsnTask() :
  40   AliAnalysisTaskSE(),
  41   fRsnPDG(0),
  42   fCFManager(0x0),
  43   fHistEventsProcessed(0x0)
  44 {
  45   //
  46   //Default ctor
  47   //
  48 }
  49 //___________________________________________________________________________
  50 AliCFRsnTask::AliCFRsnTask(const Char_t* name) :
  51   AliAnalysisTaskSE(name),
  52   fRsnPDG(0),
  53   fCFManager(0x0),
  54   fHistEventsProcessed(0x0)
  55 {
  56   //
  57   // Constructor. Initialization of Inputs and Outputs
  58   //
  59   Info("AliCFRsnTask","Calling Constructor");
  60   /*
  61     DefineInput(0) and DefineOutput(0)
  62     are taken care of by AliAnalysisTaskSE constructor
  63   */
  64   DefineOutput(1,TH1I::Class());
  65   DefineOutput(2,AliCFContainer::Class());
  66   //   DefineOutput(3,TList::Class());
  67 }
  68
  69 //___________________________________________________________________________
  70 AliCFRsnTask& AliCFRsnTask::operator=(const AliCFRsnTask& c)
  71 {
  72   //
  73   // Assignment operator
  74   //
  75   if (this!=&c) {
  76     AliAnalysisTaskSE::operator=(c) ;
  77     fRsnPDG     = c.fRsnPDG;
  78     fCFManager  = c.fCFManager;
  79     fHistEventsProcessed = c.fHistEventsProcessed;
  80   }
  81   return *this;
  82 }
  83
  84 //___________________________________________________________________________
  85 AliCFRsnTask::AliCFRsnTask(const AliCFRsnTask& c) :
  86   AliAnalysisTaskSE(c),
  87   fRsnPDG(c.fRsnPDG),
  88   fCFManager(c.fCFManager),
  89   fHistEventsProcessed(c.fHistEventsProcessed)
  90 {
  91   //
  92   // Copy Constructor
  93   //
  94 }
  95
  96 //___________________________________________________________________________
  97 AliCFRsnTask::~AliCFRsnTask() {
  98   //
  99   //destructor
 100   //
 101   Info("~AliCFRsnTask","Calling Destructor");
 102   if (fCFManager)           delete fCFManager ;
 103   if (fHistEventsProcessed) delete fHistEventsProcessed ;
 104 }
 105
 106 //_________________________________________________
 107 void AliCFRsnTask::UserExec(Option_t *)
 108 {
 109   //
 110   // Main loop function
 111   //
 112   Info("UserExec","") ;
 113
 114   AliESDEvent* fESD = dynamic_cast<AliESDEvent*>(fInputEvent);
 115   if (!fESD) {
 116     Error("UserExec","NO ESD FOUND!");
 117     return;
 118   }
 119
 120   if (!fMCEvent) Error("UserExec","NO MC INFO FOUND!");
 121   fCFManager->SetEventInfo(fMCEvent);
 122
 123   AliStack*   stack   = fMCEvent->Stack();
 124
 125   // MC-event selection
 126   Double_t containerInput[2] ;
 127
 128   //loop on the MC event
 129   Info("UserExec","Looping on MC event");
 130   for (Int_t ipart=0; ipart<stack->GetNprimary(); ipart++) {
 131     AliMCParticle *mcPart  = fMCEvent->GetTrack(ipart);
 132
 133     //check the MC-level cuts
 134     if (!fCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,mcPart)) continue;
 135     containerInput[0] = mcPart->Pt();
 136     containerInput[1] = mcPart->Y() ;
 137     //fill the container for Gen-level selection
 138     fCFManager->GetParticleContainer()->Fill(containerInput,kStepGenerated);
 139
 140     //check the Acceptance-level cuts
 141     if (!fCFManager->CheckParticleCuts(AliCFManager::kPartAccCuts,mcPart)) continue;
 142     //fill the container for Acceptance-level selection
 143     fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructible);
 144   }
 145
 146
 147   //Now go to rec level
 148   Info("UserExec","Looping on ESD event");
 149
 150   //SET THE ESD AS EVENT INFO IN RECONSTRUCTION CUTS
 151   TObjArray* fCutsReco = fCFManager->GetParticleCutsList(AliCFManager::kPartRecCuts);
 152   TObjArray* fCutsPID = fCFManager->GetParticleCutsList(AliCFManager::kPartSelCuts);
 153   TObjArrayIter iter1(fCutsReco);
 154   TObjArrayIter iter2(fCutsPID);
 155   AliCFCutBase *cut = 0;
 156   while ( (cut = (AliCFCutBase*)iter1.Next()) ) {
 157     cut->SetEvtInfo(fESD);
 158   }
 159   while ( (cut = (AliCFCutBase*)iter2.Next()) ) {
 160     cut->SetEvtInfo(fESD);
 161   }
 162
 163   // Loop on negative tracks
 164   for (Int_t iTrack1 = 0; iTrack1<fESD->GetNumberOfTracks(); iTrack1++) {
 165     AliESDtrack* esdTrack1 = fESD->GetTrack(iTrack1);
 166     //track1 is negative
 167     if (esdTrack1->Charge()>=0) continue;
 168     Int_t esdLabel1 = esdTrack1->GetLabel();
 169     if (esdLabel1<0) continue;
 170
 171     //Loop on positive tracks
 172     for (Int_t iTrack2 = 0; iTrack2<fESD->GetNumberOfTracks(); iTrack2++) {
 173       AliESDtrack* esdTrack2 = fESD->GetTrack(iTrack2);
 174       //track2 is positive
 175       if (esdTrack2->Charge()<=0) continue;
 176       Int_t esdLabel2 = esdTrack2->GetLabel();
 177       if (esdLabel2<0) continue;
 178
 179       //Create Resonance daughter objects
 180       AliRsnDaughter* tmp1 = new AliRsnDaughter(esdTrack1);
 181       AliRsnDaughter* tmp2 = new AliRsnDaughter(esdTrack2);
 182       AliRsnDaughter track1(*tmp1);
 183       AliRsnDaughter track2(*tmp2);
 184       delete tmp1;
 185       delete tmp2;
 186
 187       //Set MC information to resonance daughters
 188       TParticle *part1 = stack->Particle(esdLabel1);
 189       track1.InitParticle(part1);
 190       track1.GetParticle()->SetPDG(part1->GetPdgCode());
 191
 192       Int_t mother1 = part1->GetFirstMother();
 193       track1.GetParticle()->SetMother(mother1);
 194       if (mother1 >= 0) {
 195         TParticle *mum = stack->Particle(mother1);
 196         track1.GetParticle()->SetMotherPDG(mum->GetPdgCode());
 197       }
 198
 199       TParticle *part2 = stack->Particle(esdLabel2);
 200       track2.InitParticle(part2);
 201       track2.GetParticle()->SetPDG(part2->GetPdgCode());
 202
 203       Int_t mother2 = part2->GetFirstMother();
 204       track2.GetParticle()->SetMother(mother2);
 205       if (mother2 >= 0) {
 206         TParticle *mum = stack->Particle(mother2);
 207         track2.GetParticle()->SetMotherPDG(mum->GetPdgCode());
 208       }
 209
 210       //make a mother resonance from the 2 candidate daughters
 211       AliRsnDaughter rsn = AliRsnDaughter::Sum(track1,track2) ;
 212       AliCFPair pair(esdTrack1,esdTrack2); // This object is used for cuts (to be replaced)
 213
 214       //check if true resonance
 215       if (rsn.GetParticle()->PDG() != fRsnPDG) continue;
 216       if (!fCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,&pair)) continue;
 217
 218       //check if associated MC resonance passes the cuts
 219       Int_t motherLabel=rsn.Label();
 220       if (motherLabel<0) continue ;
 221       AliMCParticle* mcRsn = fMCEvent->GetTrack(motherLabel);
 222       if (!mcRsn) continue;
 223       if (!fCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,mcRsn)) continue;
 224
 225       //fill the container
 226       containerInput[0] = rsn.Pt() ;
 227       containerInput[1] = GetRapidity(rsn.E(),rsn.Pz());
 228       fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed) ;
 229
 230       if (!fCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,&pair)) continue ;
 231       fCFManager->GetParticleContainer()->Fill(containerInput,kStepSelected);
 232     }
 233   }
 234
 235   fHistEventsProcessed->Fill(0);
 236   /* PostData(0) is taken care of by AliAnalysisTaskSE */
 237   PostData(1,fHistEventsProcessed) ;
 238   PostData(2,fCFManager->GetParticleContainer()) ;
 239
 240   //   TList * list = new TList();
 241   //   fCFManager->AddQAHistosToList(list);
 242   //   PostData(2,list) ;
 243 }
 244
 245
 246 //___________________________________________________________________________
 247 void AliCFRsnTask::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
 255   Double_t max1 = fCFManager->GetParticleContainer()->ShowProjection(0,0)->GetMaximum();
 256   Double_t max2 = fCFManager->GetParticleContainer()->ShowProjection(1,0)->GetMaximum();
 257
 258   fCFManager->GetParticleContainer()->ShowProjection(0,0)->GetYaxis()->SetRangeUser(0,max1*1.2);
 259   fCFManager->GetParticleContainer()->ShowProjection(0,1)->GetYaxis()->SetRangeUser(0,max1*1.2);
 260   fCFManager->GetParticleContainer()->ShowProjection(0,2)->GetYaxis()->SetRangeUser(0,max1*1.2);
 261   fCFManager->GetParticleContainer()->ShowProjection(0,3)->GetYaxis()->SetRangeUser(0,max1*1.2);
 262
 263   fCFManager->GetParticleContainer()->ShowProjection(1,0)->GetYaxis()->SetRangeUser(0,max2*1.2);
 264   fCFManager->GetParticleContainer()->ShowProjection(1,1)->GetYaxis()->SetRangeUser(0,max2*1.2);
 265   fCFManager->GetParticleContainer()->ShowProjection(1,2)->GetYaxis()->SetRangeUser(0,max2*1.2);
 266   fCFManager->GetParticleContainer()->ShowProjection(1,3)->GetYaxis()->SetRangeUser(0,max2*1.2);
 267
 268   fCFManager->GetParticleContainer()->ShowProjection(0,0)->SetMarkerStyle(23) ;
 269   fCFManager->GetParticleContainer()->ShowProjection(0,1)->SetMarkerStyle(24) ;
 270   fCFManager->GetParticleContainer()->ShowProjection(0,2)->SetMarkerStyle(25) ;
 271   fCFManager->GetParticleContainer()->ShowProjection(0,3)->SetMarkerStyle(26) ;
 272
 273   fCFManager->GetParticleContainer()->ShowProjection(1,0)->SetMarkerStyle(23) ;
 274   fCFManager->GetParticleContainer()->ShowProjection(1,1)->SetMarkerStyle(24) ;
 275   fCFManager->GetParticleContainer()->ShowProjection(1,2)->SetMarkerStyle(25) ;
 276   fCFManager->GetParticleContainer()->ShowProjection(1,3)->SetMarkerStyle(26) ;
 277
 278   TCanvas * c =new TCanvas("c","",1400,800);
 279   c->Divide(4,2);
 280
 281   c->cd(1);
 282   fCFManager->GetParticleContainer()->ShowProjection(0,0)->Draw("p");
 283   c->cd(2);
 284   fCFManager->GetParticleContainer()->ShowProjection(0,1)->Draw("p");
 285   c->cd(3);
 286   fCFManager->GetParticleContainer()->ShowProjection(0,2)->Draw("p");
 287   c->cd(4);
 288   fCFManager->GetParticleContainer()->ShowProjection(0,3)->Draw("p");
 289   c->cd(5);
 290   fCFManager->GetParticleContainer()->ShowProjection(1,0)->Draw("p");
 291   c->cd(6);
 292   fCFManager->GetParticleContainer()->ShowProjection(1,1)->Draw("p");
 293   c->cd(7);
 294   fCFManager->GetParticleContainer()->ShowProjection(1,2)->Draw("p");
 295   c->cd(8);
 296   fCFManager->GetParticleContainer()->ShowProjection(1,3)->Draw("p");
 297
 298   c->SaveAs("plots.eps");
 299
 300   delete fHistEventsProcessed ;
 301 }
 302
 303 //___________________________________________________________________________
 304 void AliCFRsnTask::UserCreateOutputObjects() {
 305   //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
 306   //TO BE SET BEFORE THE EXECUTION OF THE TASK
 307   //
 308   Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
 309
 310   //slot #1
 311   OpenFile(1);
 312   fHistEventsProcessed = new TH1I("fHistEventsProcessed","",1,0,1) ;
 313 }
 314
 315 //___________________________________________________________________________
 316 Double_t AliCFRsnTask::GetRapidity(Double_t energy, Double_t pz) {
 317   if (energy == pz || energy == -pz) {
 318     printf("GetRapidity : ERROR : rapidity for 4-vector with E = Pz -- infinite result");
 319     return 999;
 320   }
 321   if (energy < pz) {
 322     printf("GetRapidity : ERROR : rapidity for 4-vector with E = Pz -- infinite result");
 323     return 999;
 324   }
 325   Double_t y = 0.5 * log((energy + pz) / (energy - pz));
 326   return y;
 327 }
 328
 329