PWG3/muon/AliCFMuonSingleTask1.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 locally, on AliEn and 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 // Modification done by X. Lopez - LPC Clermont (fr)
  26 //-----------------------------------------------------------------------
  27
  28 #ifndef ALICFMUONSINGLETASK1_CXX
  29 #define ALICFMUONSINGLETASK1_CXX
  30
  31 #include "AliCFMuonSingleTask1.h"
  32 #include "AliHeader.h"
  33 #include "AliESDHeader.h"
  34 #include "AliStack.h"
  35 #include "TParticle.h"
  36 #include "TLorentzVector.h"
  37 #include "TH1I.h"
  38 #include "AliMCEvent.h"
  39 #include "AliAnalysisManager.h"
  40 #include "AliESDEvent.h"
  41 #include "AliAODEvent.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 "AliESDMuonTrack.h"
  49 #include "AliESDtrack.h"
  50 #include "AliESDInputHandler.h"
  51 #include "TCanvas.h"
  52
  53 ClassImp(AliCFMuonSingleTask1)
  54
  55 //__________________________________________________________________________
  56 AliCFMuonSingleTask1::AliCFMuonSingleTask1() :
  57   fReadAODData(0),
  58   fCFManager(0x0),
  59   fQAHistList(0x0),
  60   fHistEventsProcessed(0x0),
  61   fNevt(0),
  62   fIsMC(kFALSE)
  63 {
  64   //
  65   //Default ctor
  66   //
  67 }
  68 //___________________________________________________________________________
  69 AliCFMuonSingleTask1::AliCFMuonSingleTask1(const Char_t* name) :
  70   AliAnalysisTaskSE(name),
  71   fReadAODData(0),
  72   fCFManager(0x0),
  73   fQAHistList(0x0),
  74   fHistEventsProcessed(0x0),
  75   fNevt(0),
  76   fIsMC(kFALSE)
  77 {
  78   //
  79   // Constructor. Initialization of Inputs and Outputs
  80   //
  81   Info("AliCFMuonSingleTask1","Calling Constructor");
  82
  83   fHistEventsProcessed = new TH1I("fHistEventsProcessed","",1,0,1) ;
  84
  85   DefineOutput(1,TH1I::Class());
  86   DefineOutput(2,AliCFContainer::Class());
  87
  88 }
  89
  90 //___________________________________________________________________________
  91 AliCFMuonSingleTask1& AliCFMuonSingleTask1::operator=(const AliCFMuonSingleTask1& c)
  92 {
  93   //
  94   // Assignment operator
  95   //
  96   if (this!=&c) {
  97     AliAnalysisTaskSE::operator=(c) ;
  98     fReadAODData = c.fReadAODData ;
  99     fCFManager  = c.fCFManager;
 100     fQAHistList = c.fQAHistList ;
 101     fHistEventsProcessed = c.fHistEventsProcessed;
 102     fNevt = c.fNevt ;
 103   }
 104   return *this;
 105 }
 106
 107 //___________________________________________________________________________
 108 AliCFMuonSingleTask1::AliCFMuonSingleTask1(const AliCFMuonSingleTask1& c) :
 109   AliAnalysisTaskSE(c),
 110   fReadAODData(c.fReadAODData),
 111   fCFManager(c.fCFManager),
 112   fQAHistList(c.fQAHistList),
 113   fHistEventsProcessed(c.fHistEventsProcessed),
 114   fNevt(c.fNevt),
 115   fIsMC(kFALSE)
 116 {
 117   //
 118   // Copy Constructor
 119   //
 120 }
 121
 122 //___________________________________________________________________________
 123 AliCFMuonSingleTask1::~AliCFMuonSingleTask1() {
 124   //
 125   //destructor
 126   //
 127   Info("~AliCFMuonSingleTask1","Calling Destructor");
 128   if (fCFManager)           delete fCFManager ;
 129   if (fHistEventsProcessed) delete fHistEventsProcessed ;
 130   if (fQAHistList) {fQAHistList->Clear(); delete fQAHistList;}
 131 }
 132
 133 //_________________________________________________
 134 void AliCFMuonSingleTask1::UserExec(Option_t *)
 135 {
 136   //
 137   // Main loop function
 138   //
 139
 140   Info("UserExec"," ") ;
 141   if (!fMCEvent && fIsMC) {
 142     Error("UserExec","NO MC EVENT FOUND!");
 143     return;
 144   }
 145
 146   fNevt++;
 147   fCFManager->SetEventInfo(fMCEvent);
 148   Double_t containerInput[18] ;
 149
 150 ////////
 151 //// MC
 152 ////////
 153
 154 // loop on the MC part
 155
 156   if(fIsMC){
 157
 158       for (Int_t ipart=0; ipart<fMCEvent->GetNumberOfTracks(); ipart++) {
 159
 160           AliMCParticle *mcPart  = (AliMCParticle*) fMCEvent->GetTrack(ipart);
 161           TParticle *part = mcPart->Particle();
 162
 163           // rapidity and Pt cuts (default -4<y<-2.5 et 0<pt<20)
 164           if (!fCFManager->CheckParticleCuts(AliCFManager::kPartAccCuts,mcPart)) continue;
 165
 166           // Selection of mu
 167           Int_t pdg=part->GetPdgCode();
 168           if(TMath::Abs(pdg)==13){
 169
 170               Float_t emc = part->Energy();
 171               Float_t pzmc = part->Pz();
 172               Float_t rapmc = Rap(emc,pzmc);
 173               Float_t phimc = part->Phi();
 174               phimc = Phideg(phimc);
 175               Float_t ptmc = part->Pt();
 176               Float_t pmc = part->P();
 177               Float_t zmc = part->Vz();
 178               Float_t etamc = part->Eta();
 179               Float_t chargemc=0;
 180               if(pdg==13) chargemc=-1;
 181               if(pdg==-13) chargemc=1;
 182               Float_t thetamc;
 183               thetamc = TMath::Abs(-TMath::Pi()+TMath::ASin(ptmc/pmc)); // convention for EVGEN
 184               Float_t ymc = part->Vy();
 185               Float_t xmc = part->Vx();
 186               Float_t rmc = TMath::Sqrt(xmc*xmc+ymc*ymc);
 187
 188               containerInput[0] = etamc ;
 189               containerInput[1] = rapmc ;
 190               containerInput[2] = phimc ;
 191               containerInput[3] = ptmc ;
 192               containerInput[4] = pmc ;
 193               containerInput[14] = zmc ;
 194               containerInput[16] = chargemc ;
 195               containerInput[17] = thetamc ;
 196
 197               // 10 var calculated only for ESD .i.e set at 1 in MC step
 198               for(Int_t i=5; i<=13; i++){
 199                   containerInput[i] = 1;
 200               }
 201
 202               containerInput[7] = 3;
 203               containerInput[10] = rmc ; // radius at production
 204               containerInput[15] = 45 ;  // rabsmc
 205
 206               // fill the container at the first step
 207               fCFManager->GetParticleContainer()->Fill(containerInput,0);
 208           }
 209       }
 210   }
 211
 212
 213 ////////
 214 //// ESD
 215 ////////
 216
 217   Int_t trig1MuL = 1 ,trig1MuH = 2, trigUSL = 4, trigUSH = 8, trigLSL = 16, trigLSH = 32;
 218   Int_t trig = 0;
 219
 220   AliESDEvent *fESD;
 221   AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*>
 222       (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 223   fESD = esdH->GetEvent();
 224
 225 // trigger
 226
 227   if (fESD->GetTriggerMask() & trig1MuL) trig = 1;
 228   if (fESD->GetTriggerMask() & trig1MuH) trig = 2;
 229   if (fESD->GetTriggerMask() & trigUSL)  trig = 3;
 230   if (fESD->GetTriggerMask() & trigUSH)  trig = 4;
 231   if (fESD->GetTriggerMask() & trigLSL)  trig = 5;
 232   if (fESD->GetTriggerMask() & trigLSH)  trig = 6;
 233
 234 // vertex
 235
 236   Float_t vx = -200 , vy = -200 , vz = -200, vt=-200;
 237
 238   Double_t Ncont = 0.;
 239   AliESDVertex* vertex = (AliESDVertex*) fESD->GetVertex();
 240   if (vertex->GetNContributors()) {
 241       vz = vertex->GetZv();
 242       vy = vertex->GetYv();
 243       vx = vertex->GetXv();
 244       vt=TMath::Sqrt(vx*vx+vy*vy);
 245       Ncont=vertex->GetNContributors();
 246   }
 247
 248 // tracks
 249
 250   Int_t mult1 = fESD->GetNumberOfMuonTracks() ;
 251
 252     for (Int_t j = 0; j<mult1; j++) {
 253         AliESDMuonTrack* mu1 = new AliESDMuonTrack(*(fESD->GetMuonTrack(j)));
 254
 255         if (!mu1->ContainTrackerData()) { mu1=0; continue; }
 256
 257         Float_t charger = mu1->Charge();
 258         Float_t er = mu1->E();
 259         Float_t pzr = mu1->Pz();
 260         Float_t rapr=Rap(er,pzr);
 261         Float_t phir = mu1->Phi();
 262         phir = Phideg(phir);
 263         Float_t ptr = mu1->Pt();
 264         Float_t pr = mu1->P();
 265         Float_t hit = mu1->GetNHit();
 266         Float_t chi2 = mu1->GetChi2() / (2.0 * hit - 5);
 267         Int_t matchtrig = mu1->GetMatchTrigger();
 268         Float_t chi2match = mu1->GetChi2MatchTrigger();
 269         Float_t dcar = mu1->GetDCA();
 270         Float_t zr = mu1->GetZ();
 271         Float_t etar = mu1->Eta();
 272         Float_t rabs = mu1->GetRAtAbsorberEnd();
 273         Float_t thetar ;
 274         thetar = TMath::Abs(-TMath::Pi()+TMath::ASin(ptr/pr));
 275
 276 // rapidity and Pt cuts (default -4<y<-2.5 et 0<pt<20)
 277         if (!fCFManager->CheckParticleCuts(AliCFManager::kPartAccCuts,mu1)) continue;
 278
 279         containerInput[0] = etar ;
 280         containerInput[1] = rapr ;
 281         containerInput[2] = phir ;
 282         containerInput[3] = ptr ;
 283         containerInput[4] = pr ;
 284         containerInput[5] = hit ;
 285         containerInput[6] = chi2 ;
 286         containerInput[7] = matchtrig ;
 287         containerInput[8] = chi2match ;
 288         containerInput[9] = Ncont ;
 289         containerInput[10] = vt ;
 290         containerInput[11] = vz ;
 291         containerInput[12] = trig ;
 292         containerInput[13] = dcar ;
 293         containerInput[14] = zr ;
 294         containerInput[15] = rabs ;
 295         containerInput[16] = charger ;
 296         containerInput[17] = thetar ;
 297
 298 // fill the container at the second step (for simu, first for data)
 299         if(fIsMC){
 300             fCFManager->GetParticleContainer()->Fill(containerInput,1);
 301         }
 302         else{
 303             fCFManager->GetParticleContainer()->Fill(containerInput,0);
 304         }
 305
 306
 307     }      // mu Loop
 308
 309 //  ----------
 310   fHistEventsProcessed->Fill(0);
 311   PostData(1,fHistEventsProcessed) ;
 312   PostData(2,fCFManager->GetParticleContainer()) ;
 313
 314 }
 315 //________________________________________________________________________
 316 Float_t AliCFMuonSingleTask1::Rap(Float_t e, Float_t pz)
 317 {
 318 // calculate rapidity
 319     Float_t rap;
 320     if(e!=pz){
 321         rap = 0.5*TMath::Log((e+pz)/(e-pz));
 322         return rap;
 323     }
 324     else{
 325         rap = -200;
 326         return rap;
 327     }
 328 }
 329 //________________________________________________________________________
 330 Float_t AliCFMuonSingleTask1::Phideg(Float_t phi)
 331 {
 332 // calculate Phi in range [-180,180]
 333     Float_t phideg;
 334
 335         phideg = phi-TMath::Pi();
 336         phideg = phideg*57.296;
 337         return phideg;
 338 }
 339 //________________________________________________________________________
 340 void AliCFMuonSingleTask1::Terminate(Option_t *)
 341 {
 342   // project pt (var[3]) from the two steps MC(0) and ESD(1)
 343
 344     AliCFContainer *cont = dynamic_cast<AliCFContainer*> (GetOutputData(2));
 345
 346     TH1D *kpt = cont->ShowProjection(16,0);
 347     TH1D *rpt = cont->ShowProjection(16,1);
 348
 349     TCanvas *c1 = new TCanvas("AliCFMuonSingleTask1"," MC & ESD",10,10,510,510);
 350     c1->Divide(1,2);
 351     c1->cd(1);
 352     kpt->Draw("HIST");
 353     c1->cd(2);
 354     rpt->Draw("HIST");
 355
 356 }
 357 //________________________________________________________________________
 358
 359 #endif