PWG4/JetTasks/AliPWG4CosmicCandidates.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 // Analysis task looking for cosmic candidates embedded in a p-p event
  18 // Task checks if particles are back-to-back in eta and phi
  19 //
  20 //
  21 // Author : Marta Verweij - UU - marta.verweij@cern.ch
  22 //-----------------------------------------------------------------------
  23
  24 #include "TVector3.h"
  25 #include <iostream>
  26 #include "TH1.h"
  27 #include "TH2.h"
  28 #include "TH3.h"
  29 #include "TList.h"
  30 #include "TChain.h"
  31
  32 #include "AliAnalysisTask.h"
  33 #include "AliAnalysisManager.h"
  34
  35 #include "AliESDEvent.h"
  36 #include "AliESDInputHandler.h"
  37
  38 #include "AliESDtrack.h"
  39 #include "AliESDtrackCuts.h"
  40 #include "AliExternalTrackParam.h"
  41
  42 #include "AliLog.h"
  43
  44 #include "AliPWG4CosmicCandidates.h"
  45
  46 //using namespace std; //required for resolving the 'cout' symbol
  47 using namespace std;
  48
  49 ClassImp(AliPWG4CosmicCandidates)
  50
  51 //________________________________________________________________________
  52 AliPWG4CosmicCandidates::AliPWG4CosmicCandidates()
  53 : AliAnalysisTaskSE(),
  54   fTrackCuts(0),
  55   fPtMin(5.),
  56   fMaxCosmicAngle(0.002),
  57   fNEventAll(0),
  58   fNEventSel(0),
  59   fPtSignedCosmicCandidates(0),
  60   fDeltaPtCosmicCandidates(0),
  61   fDeltaPhiSumEta(0),
  62   fDCAZCosmicCandidates(0),
  63   fDCARCosmicCandidates(0),
  64   fTheta(0),
  65   fThetaZoom(0),
  66   fThetaPt1Pt2(0),
  67   fThetaPt1Pt2Signed(0),
  68   fDeltaPhiSumEtaPt1(0),
  69   fDeltaPhiSumEtaPt2(0),
  70   fThetaDCAZ1DCAZ2(0),
  71   fRisol(0),
  72   fRisolTheta(0),
  73   fHistListCosmics(0)
  74 {
  75   //
  76   // Default constructor
  77   //
  78 }
  79
  80 //________________________________________________________________________
  81 AliPWG4CosmicCandidates::AliPWG4CosmicCandidates(const char *name)
  82   : AliAnalysisTaskSE(name),
  83     fTrackCuts(0),
  84     fPtMin(5.),
  85     fMaxCosmicAngle(0.002),
  86     fNEventAll(0),
  87     fNEventSel(0),
  88     fPtSignedCosmicCandidates(0),
  89     fDeltaPtCosmicCandidates(0),
  90     fDeltaPhiSumEta(0),
  91     fDCAZCosmicCandidates(0),
  92     fDCARCosmicCandidates(0),
  93     fTheta(0),
  94     fThetaZoom(0),
  95     fThetaPt1Pt2(0),
  96     fThetaPt1Pt2Signed(0),
  97     fDeltaPhiSumEtaPt1(0),
  98     fDeltaPhiSumEtaPt2(0),
  99     fThetaDCAZ1DCAZ2(0),
 100     fRisol(0),
 101     fRisolTheta(0),
 102     fHistListCosmics(0)
 103 {
 104   // Constructor. Initialization of Inputs and Outputs
 105
 106   // Define input and output slots here
 107   // Input slot #0 works with a TChain
 108   //  DefineInput(0, TChain::Class());
 109   // Output slot #0 id reserved by the base class for AOD
 110   // Output slot #1 writes into a TList
 111   DefineOutput(1, TList::Class());
 112   // Output slot #2 writes into a AliESDtrackCuts
 113   DefineOutput(2, AliESDtrackCuts::Class());
 114 }
 115 //________________________________________________________________________
 116 AliPWG4CosmicCandidates::AliPWG4CosmicCandidates(const AliPWG4CosmicCandidates &res)
 117   : AliAnalysisTaskSE(res),
 118     fTrackCuts(0),
 119     fPtMin(5.),
 120     fMaxCosmicAngle(0.002),
 121     fNEventAll(0),
 122     fNEventSel(0),
 123     fPtSignedCosmicCandidates(0),
 124     fDeltaPtCosmicCandidates(0),
 125     fDeltaPhiSumEta(0),
 126     fDCAZCosmicCandidates(0),
 127     fDCARCosmicCandidates(0),
 128     fTheta(0),
 129     fThetaZoom(0),
 130     fThetaPt1Pt2(0),
 131     fThetaPt1Pt2Signed(0),
 132     fDeltaPhiSumEtaPt1(0),
 133     fDeltaPhiSumEtaPt2(0),
 134     fThetaDCAZ1DCAZ2(0),
 135     fRisol(0),
 136     fRisolTheta(0),
 137     fHistListCosmics(0)
 138 {
 139   //
 140   // Dummy copy constructor
 141   //
 142 }
 143
 144 //________________________________________________________________________
 145 AliPWG4CosmicCandidates& AliPWG4CosmicCandidates::operator=(const AliPWG4CosmicCandidates& /*trclass*/)
 146 {
 147   //
 148   // Dummy assignment operator
 149   //
 150   return *this;
 151 }
 152
 153 //________________________________________________________________________
 154 void AliPWG4CosmicCandidates::LocalInit()
 155 {
 156   //
 157   // Only called once at beginning
 158   //
 159   PostData(2,fTrackCuts);
 160 }
 161
 162 //________________________________________________________________________
 163 void AliPWG4CosmicCandidates::UserCreateOutputObjects()
 164 {
 165   //
 166   // Create output objects
 167   // Called once
 168   AliDebug(2,Form(">> AliPWG4CosmicCandidates::UserCreateOutputObjects \n"));
 169
 170   Bool_t oldStatus = TH1::AddDirectoryStatus();
 171   TH1::AddDirectory(kFALSE);
 172
 173   OpenFile(1);
 174   fHistListCosmics = new TList();
 175
 176   fNEventAll = new TH1F("fNEventAll","NEventAll",1,-0.5,0.5);
 177   fHistListCosmics->Add(fNEventAll);
 178   fNEventSel = new TH1F("fNEventSel","NEvent Selected for analysis",1,-0.5,0.5);
 179   fHistListCosmics->Add(fNEventSel);
 180
 181   Float_t fgkPtMin=0.;
 182   Float_t fgkPtMax=100.;
 183   Int_t fgkNPtBins= (int)(fgkPtMax-fgkPtMin);
 184
 185   Int_t fgkNPhiBins=18;
 186   Float_t kMinPhi = -0.5*TMath::Pi();
 187   Float_t kMaxPhi = 3./2.*TMath::Pi();
 188
 189   Int_t fgkNThetaBins=fgkNPhiBins*8;
 190   Float_t kMinTheta = -0.5*TMath::Pi();
 191   Float_t kMaxTheta = 3./2.*TMath::Pi();
 192
 193   Int_t fgkNDCARBins=80;
 194   Float_t fgkDCARMin = -0.2;
 195   Float_t fgkDCARMax = 0.2;
 196   Float_t *binsDCAR=new Float_t[fgkNDCARBins+1];
 197   for(Int_t i=0; i<=fgkNDCARBins; i++) binsDCAR[i]=(Float_t)fgkDCARMin + (fgkDCARMax-fgkDCARMin)/fgkNDCARBins*(Float_t)i ;
 198
 199   Int_t fgkNDCAZBins=80;
 200   Float_t fgkDCAZMin = -2.;
 201   Float_t fgkDCAZMax = 2.;
 202   Float_t *binsDCAZ=new Float_t[fgkNDCAZBins+1];
 203   for(Int_t i=0; i<=fgkNDCAZBins; i++) binsDCAZ[i]=(Float_t)fgkDCAZMin + (fgkDCAZMax-fgkDCAZMin)/fgkNDCAZBins*(Float_t)i ;
 204
 205   fPtSignedCosmicCandidates = new TH1F("fPtSignedCosmicCandidates","fPtSignedCosmicCandidates",2*(int)(fgkPtMax-fgkPtMin), -1.*fgkPtMax, fgkPtMax);
 206   fHistListCosmics->Add(fPtSignedCosmicCandidates);
 207
 208   fDeltaPtCosmicCandidates = new TH1F("fDeltaPtCosmicCandidates","fDeltaPtCosmicCandidates",fgkNPtBins, -50., 50.);
 209   fHistListCosmics->Add(fDeltaPtCosmicCandidates);
 210
 211   fDeltaPhiSumEta = new TH2F("fDeltaPhiSumEta","fDeltaPhiSumEta",fgkNPhiBins*4,kMinPhi,kMaxPhi,80, -2.,2.);
 212   fHistListCosmics->Add(fDeltaPhiSumEta);
 213
 214   fDCAZCosmicCandidates = new TH2F("fDCAZCosmicCandidates","fDCAZCosmicCandidates",fgkNDCAZBins,binsDCAZ,fgkNDCAZBins,binsDCAZ);
 215   fHistListCosmics->Add(fDCAZCosmicCandidates);
 216
 217   fDCARCosmicCandidates = new TH2F("fDCARCosmicCandidates","fDCARCosmicCandidates",fgkNDCARBins,binsDCAR,fgkNDCARBins,binsDCAR);
 218   fHistListCosmics->Add(fDCARCosmicCandidates);
 219
 220   fTheta = new TH1F("fTheta","fTheta",fgkNThetaBins,kMinTheta,kMaxTheta);
 221   fHistListCosmics->Add(fTheta);
 222
 223   fThetaZoom = new TH1F("fThetaZoom","fThetaZoom",100,TMath::Pi()-1.,TMath::Pi()+1.);
 224   fHistListCosmics->Add(fThetaZoom);
 225
 226   fThetaPt1Pt2 = new TH3F("fThetaPt1Pt2","fThetaPt1Pt2",fgkNThetaBins,kMinTheta,kMaxTheta,(int)(fgkPtMax-fgkPtMin),fgkPtMin,fgkPtMax,(int)(fgkPtMax-fgkPtMin),fgkPtMin,fgkPtMax);
 227   fHistListCosmics->Add(fThetaPt1Pt2);
 228
 229   fThetaPt1Pt2Signed = new TH3F("fThetaPt1Pt2Signed","fThetaPt1Pt2Signed",fgkNThetaBins,kMinTheta,kMaxTheta,4*(int)(fgkPtMax-fgkPtMin),-1.*fgkPtMax,fgkPtMax,4*(int)(fgkPtMax-fgkPtMin),-1.*fgkPtMax,fgkPtMax);
 230   fHistListCosmics->Add(fThetaPt1Pt2Signed);
 231
 232   fDeltaPhiSumEtaPt1 = new TH3F("fDeltaPhiSumEtaPt1","fDeltaPhiSumEtaPt1",fgkNThetaBins,kMinTheta,kMaxTheta,80, -2.,2.,(int)(fgkPtMax-fgkPtMin),fgkPtMin,fgkPtMax);
 233   fHistListCosmics->Add(fDeltaPhiSumEtaPt1);
 234
 235   fDeltaPhiSumEtaPt2 = new TH3F("fDeltaPhiSumEtaPt2","fDeltaPhiSumEtaPt2",fgkNThetaBins,kMinTheta,kMaxTheta,80, -2.,2.,(int)(fgkPtMax-fgkPtMin),fgkPtMin,fgkPtMax);
 236   fHistListCosmics->Add(fDeltaPhiSumEtaPt2);
 237
 238   fThetaDCAZ1DCAZ2 = new TH3F("fThetaDCAZ1DCAZ2","fThetaDCAZ1DCAZ2",fgkNThetaBins,kMinTheta,kMaxTheta,fgkNDCAZBins,-2.,2.,fgkNDCAZBins,-2.,2.);
 239   fHistListCosmics->Add(fThetaDCAZ1DCAZ2);
 240
 241   fRisol = new TH1F("fRisol","fRisol",100,0.,10.);
 242   fHistListCosmics->Add(fRisol);
 243
 244   fRisolTheta = new TH2F("fRisolTheta","fRisolTheta",100,0.,10.,fgkNThetaBins,kMinTheta,kMaxTheta);
 245   fHistListCosmics->Add(fRisolTheta);
 246
 247   TH1::AddDirectory(oldStatus);
 248
 249   if(binsDCAR) delete [] binsDCAR;
 250   if(binsDCAZ) delete [] binsDCAZ;
 251
 252 }
 253
 254 //________________________________________________________________________
 255 void AliPWG4CosmicCandidates::UserExec(Option_t *)
 256 {
 257   // Main loop
 258   // Called for each event
 259   //
 260
 261   // All events without selection
 262   fNEventAll->Fill(0.);
 263
 264   if (!fInputEvent) {
 265     AliDebug(2,Form("ERROR: fESD not available"));
 266     cout << "ERROR: fESD not available" << endl;
 267     PostData(1, fHistListCosmics);
 268     return;
 269   }
 270
 271   const AliVVertex *vtx = fInputEvent->GetPrimaryVertex();
 272   // Need vertex cut
 273   if (!vtx || vtx->GetNContributors() < 2) {
 274     // Post output data
 275     PostData(1, fHistListCosmics);
 276     return;
 277   }
 278
 279   //  AliDebug(2,Form("Vertex title %s, status %d, nCont %d\n",vtx->GetTitle(), vtx->GetStatus(), vtx->GetNContributors()));
 280   double primVtx[3];
 281   vtx->GetXYZ(primVtx);
 282   if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
 283     // Post output data
 284     PostData(1, fHistListCosmics);
 285     return;
 286   }
 287   if(!fInputEvent->GetNumberOfTracks() || fInputEvent->GetNumberOfTracks()<2){
 288     // Post output data
 289     PostData(1, fHistListCosmics);
 290     return;
 291   }
 292   Int_t nTracks = fInputEvent->GetNumberOfTracks();
 293
 294   if(!fTrackCuts) {
 295    // Post output data
 296     PostData(1, fHistListCosmics);
 297     return;
 298   }
 299
 300   fNEventSel->Fill(0.);
 301
 302   Float_t dcaR[2] = {0.,0.};
 303   Float_t dcaZ[2] = {0.,0.};
 304
 305   for (Int_t iTrack1 = 0; iTrack1 < nTracks; iTrack1++) {
 306
 307     AliESDtrack* track1 = (AliESDtrack*)fInputEvent->GetTrack(iTrack1);
 308     if (!track1)  continue;
 309     if(!(fTrackCuts->AcceptTrack(track1))) { continue; }
 310     if(track1->Pt()<fPtMin) continue;
 311     //Start 2nd track loop to look for correlations
 312     for (Int_t iTrack2 = iTrack1+1; iTrack2 < nTracks; iTrack2++) {
 313       AliESDtrack *track2 = (AliESDtrack*)fInputEvent->GetTrack(iTrack2);
 314       if(!track2) continue;
 315       if(!(fTrackCuts->AcceptTrack(track2))) { continue; }
 316
 317       //Check if back-to-back
 318       Double_t mom1[3],mom2[3];
 319       track1->GetPxPyPz(mom1);
 320       track2->GetPxPyPz(mom2);
 321       //     Double_t cosTheta = (mom1[0]*mom2[0]+mom1[1]*mom2[1]+mom1[2]*mom2[2])/( TMath::Sqrt(mom1[0]*mom1[0]+mom1[1]*mom1[1]+mom1[2]*mom1[2])*TMath::Sqrt(mom2[0]*mom2[0]+mom2[1]*mom2[1]+mom2[2]*mom2[2]) );
 322       TVector3 momv1(mom1[0],mom1[1],mom1[2]);
 323       TVector3 momv2(mom2[0],mom2[1],mom2[2]);
 324       //Double_t theta = momv1.Angle(momv2);
 325       Double_t theta = momv1.Phi()-momv2.Phi();
 326       if(theta<-0.5*TMath::Pi()) theta+=2.*TMath::Pi();
 327
 328       fDeltaPtCosmicCandidates->Fill(track1->Pt()-track2->Pt());
 329       Float_t deltaPhi = track1->Phi()-track2->Phi();
 330       if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
 331       fDeltaPhiSumEta->Fill(deltaPhi,track1->Eta()+track2->Eta());
 332
 333       track1->GetImpactParameters(dcaR[0],dcaZ[0]);
 334       track2->GetImpactParameters(dcaR[1],dcaZ[1]);
 335
 336       if(track2->Pt()<0.5) continue;
 337       Double_t rIsol = TMath::Sqrt( deltaPhi*deltaPhi+(track1->Eta()-track2->Eta())*(track1->Eta()-track2->Eta()) );
 338       fRisol->Fill(rIsol); //Fill R histogram
 339       if(track2->Pt()<fPtMin) continue;
 340
 341         fTheta->Fill(theta);
 342         fThetaZoom->Fill(theta);
 343         fThetaPt1Pt2->Fill(theta,track1->Pt(),track2->Pt());
 344         fThetaPt1Pt2Signed->Fill(theta,track1->GetSign()*track1->Pt(),track2->GetSign()*track2->Pt());
 345         fDeltaPhiSumEtaPt1->Fill(deltaPhi,track1->Eta()+track2->Eta(),track1->Pt());
 346         fDeltaPhiSumEtaPt2->Fill(deltaPhi,track1->Eta()+track2->Eta(),track2->Pt());
 347         fThetaDCAZ1DCAZ2->Fill(theta,dcaZ[0],dcaZ[1]);
 348         fRisolTheta->Fill(rIsol,theta);
 349         if(TMath::Abs(TMath::Pi()-theta)<fMaxCosmicAngle) {
 350           fDCAZCosmicCandidates->Fill(dcaZ[0],dcaZ[1]);
 351           fDCARCosmicCandidates->Fill(dcaR[0],dcaR[1]);
 352           fPtSignedCosmicCandidates->Fill(track1->GetSign()*track1->Pt());
 353         }
 354
 355     } // track2 loop
 356
 357   } // track1 loop
 358
 359   PostData(1,fHistListCosmics);
 360 }
 361
 362 //________________________________________________________________________
 363 void AliPWG4CosmicCandidates::Terminate(Option_t *)
 364 {
 365   //
 366   // Called once at the end of the query
 367   //
 368
 369 }