PWGJE/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   Int_t fgkNPtBins2D= (int)((fgkPtMax-fgkPtMin)/4.);
 185
 186   Int_t fgkNPhiBins=18;
 187   Float_t kMinPhi = -0.5*TMath::Pi();
 188   Float_t kMaxPhi = 3./2.*TMath::Pi();
 189
 190   Int_t fgkNThetaBins=fgkNPhiBins*8;
 191   Float_t kMinTheta = -0.5*TMath::Pi();
 192   Float_t kMaxTheta = 3./2.*TMath::Pi();
 193
 194   Int_t fgkNDCARBins=40;
 195   Float_t fgkDCARMin = -0.2;
 196   Float_t fgkDCARMax = 0.2;
 197   Float_t *binsDCAR=new Float_t[fgkNDCARBins+1];
 198   for(Int_t i=0; i<=fgkNDCARBins; i++) binsDCAR[i]=(Float_t)fgkDCARMin + (fgkDCARMax-fgkDCARMin)/fgkNDCARBins*(Float_t)i ;
 199
 200   Int_t fgkNDCAZBins=40;
 201   Float_t fgkDCAZMin = -2.;
 202   Float_t fgkDCAZMax = 2.;
 203   Float_t *binsDCAZ=new Float_t[fgkNDCAZBins+1];
 204   for(Int_t i=0; i<=fgkNDCAZBins; i++) binsDCAZ[i]=(Float_t)fgkDCAZMin + (fgkDCAZMax-fgkDCAZMin)/fgkNDCAZBins*(Float_t)i ;
 205
 206   fPtSignedCosmicCandidates = new TH1F("fPtSignedCosmicCandidates","fPtSignedCosmicCandidates",2*fgkNPtBins, -1.*fgkPtMax, fgkPtMax);
 207   fHistListCosmics->Add(fPtSignedCosmicCandidates);
 208
 209   fDeltaPtCosmicCandidates = new TH1F("fDeltaPtCosmicCandidates","fDeltaPtCosmicCandidates",fgkNPtBins, -50., 50.);
 210   fHistListCosmics->Add(fDeltaPtCosmicCandidates);
 211
 212   fDeltaPhiSumEta = new TH2F("fDeltaPhiSumEta","fDeltaPhiSumEta",fgkNPhiBins*4,kMinPhi,kMaxPhi,80, -2.,2.);
 213   fHistListCosmics->Add(fDeltaPhiSumEta);
 214
 215   fDCAZCosmicCandidates = new TH2F("fDCAZCosmicCandidates","fDCAZCosmicCandidates",fgkNDCAZBins,binsDCAZ,fgkNDCAZBins,binsDCAZ);
 216   fHistListCosmics->Add(fDCAZCosmicCandidates);
 217
 218   fDCARCosmicCandidates = new TH2F("fDCARCosmicCandidates","fDCARCosmicCandidates",fgkNDCARBins,binsDCAR,fgkNDCARBins,binsDCAR);
 219   fHistListCosmics->Add(fDCARCosmicCandidates);
 220
 221   fTheta = new TH1F("fTheta","fTheta",fgkNThetaBins,kMinTheta,kMaxTheta);
 222   fHistListCosmics->Add(fTheta);
 223
 224   fThetaZoom = new TH1F("fThetaZoom","fThetaZoom",100,TMath::Pi()-1.,TMath::Pi()+1.);
 225   fHistListCosmics->Add(fThetaZoom);
 226
 227   fThetaPt1Pt2 = new TH3F("fThetaPt1Pt2","fThetaPt1Pt2",fgkNThetaBins,kMinTheta,kMaxTheta,fgkNPtBins2D,fgkPtMin,fgkPtMax,fgkNPtBins2D,fgkPtMin,fgkPtMax);
 228   fHistListCosmics->Add(fThetaPt1Pt2);
 229
 230   fThetaPt1Pt2Signed = new TH3F("fThetaPt1Pt2Signed","fThetaPt1Pt2Signed",fgkNThetaBins,kMinTheta,kMaxTheta,fgkNPtBins2D,-1.*fgkPtMax,fgkPtMax,fgkNPtBins2D,-1.*fgkPtMax,fgkPtMax);
 231   fHistListCosmics->Add(fThetaPt1Pt2Signed);
 232
 233   fDeltaPhiSumEtaPt1 = new TH3F("fDeltaPhiSumEtaPt1","fDeltaPhiSumEtaPt1",fgkNThetaBins,kMinTheta,kMaxTheta,80, -2.,2.,fgkNPtBins2D,fgkPtMin,fgkPtMax);
 234   fHistListCosmics->Add(fDeltaPhiSumEtaPt1);
 235
 236   fDeltaPhiSumEtaPt2 = new TH3F("fDeltaPhiSumEtaPt2","fDeltaPhiSumEtaPt2",fgkNThetaBins,kMinTheta,kMaxTheta,80, -2.,2.,fgkNPtBins2D,fgkPtMin,fgkPtMax);
 237   fHistListCosmics->Add(fDeltaPhiSumEtaPt2);
 238
 239   fThetaDCAZ1DCAZ2 = new TH3F("fThetaDCAZ1DCAZ2","fThetaDCAZ1DCAZ2",fgkNThetaBins,kMinTheta,kMaxTheta,fgkNDCAZBins,-2.,2.,fgkNDCAZBins,-2.,2.);
 240   fHistListCosmics->Add(fThetaDCAZ1DCAZ2);
 241
 242   fRisol = new TH1F("fRisol","fRisol",100,0.,10.);
 243   fHistListCosmics->Add(fRisol);
 244
 245   fRisolTheta = new TH2F("fRisolTheta","fRisolTheta",100,0.,10.,fgkNThetaBins,kMinTheta,kMaxTheta);
 246   fHistListCosmics->Add(fRisolTheta);
 247
 248   TH1::AddDirectory(oldStatus);
 249
 250   PostData(1, fHistListCosmics);
 251
 252   if(binsDCAR) delete [] binsDCAR;
 253   if(binsDCAZ) delete [] binsDCAZ;
 254
 255 }
 256
 257 //________________________________________________________________________
 258 void AliPWG4CosmicCandidates::UserExec(Option_t *)
 259 {
 260   // Main loop
 261   // Called for each event
 262   //
 263
 264   // All events without selection
 265   fNEventAll->Fill(0.);
 266
 267   if (!fInputEvent) {
 268     AliDebug(2,Form("ERROR: fESD not available"));
 269     cout << "ERROR: fESD not available" << endl;
 270     PostData(1, fHistListCosmics);
 271     return;
 272   }
 273
 274   const AliVVertex *vtx = fInputEvent->GetPrimaryVertex();
 275   if(!vtx){
 276     PostData(1, fHistListCosmics);
 277     return;
 278   }
 279
 280 // Need vertex cut
 281   TString vtxName(vtx->GetName());
 282   if( vtx->GetNContributors() < 2 || (vtxName.Contains("TPCVertex")) ) {
 283     // Post output data
 284     PostData(1, fHistListCosmics);
 285     return;
 286   }
 287
 288   //  AliDebug(2,Form("Vertex title %s, status %d, nCont %d\n",vtx->GetTitle(), vtx->GetStatus(), vtx->GetNContributors()));
 289   double primVtx[3];
 290   vtx->GetXYZ(primVtx);
 291   if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
 292     // Post output data
 293     PostData(1, fHistListCosmics);
 294     return;
 295   }
 296   if(!fInputEvent->GetNumberOfTracks() || fInputEvent->GetNumberOfTracks()<2){
 297     // Post output data
 298     PostData(1, fHistListCosmics);
 299     return;
 300   }
 301   Int_t nTracks = fInputEvent->GetNumberOfTracks();
 302
 303   if(!fTrackCuts) {
 304    // Post output data
 305     PostData(1, fHistListCosmics);
 306     return;
 307   }
 308
 309   fNEventSel->Fill(0.);
 310
 311   Float_t dcaR[2] = {0.,0.};
 312   Float_t dcaZ[2] = {0.,0.};
 313
 314   for (Int_t iTrack1 = 0; iTrack1 < nTracks; iTrack1++) {
 315
 316     AliESDtrack* track1 = (AliESDtrack*)fInputEvent->GetTrack(iTrack1);
 317     if (!track1)  continue;
 318     if(!(fTrackCuts->AcceptTrack(track1))) { continue; }
 319     if(track1->Pt()<fPtMin) continue;
 320     //Start 2nd track loop to look for correlations
 321     for (Int_t iTrack2 = iTrack1+1; iTrack2 < nTracks; iTrack2++) {
 322       AliESDtrack *track2 = (AliESDtrack*)fInputEvent->GetTrack(iTrack2);
 323       if(!track2) continue;
 324       if(!(fTrackCuts->AcceptTrack(track2))) { continue; }
 325
 326       //Check if back-to-back
 327       Double_t mom1[3],mom2[3];
 328       track1->GetPxPyPz(mom1);
 329       track2->GetPxPyPz(mom2);
 330       //     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]) );
 331       TVector3 momv1(mom1[0],mom1[1],mom1[2]);
 332       TVector3 momv2(mom2[0],mom2[1],mom2[2]);
 333       //Double_t theta = momv1.Angle(momv2);
 334       Double_t theta = momv1.Phi()-momv2.Phi();
 335       if(theta<-0.5*TMath::Pi()) theta+=2.*TMath::Pi();
 336
 337       fDeltaPtCosmicCandidates->Fill(track1->Pt()-track2->Pt());
 338       Float_t deltaPhi = track1->Phi()-track2->Phi();
 339       if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
 340       fDeltaPhiSumEta->Fill(deltaPhi,track1->Eta()+track2->Eta());
 341
 342       track1->GetImpactParameters(dcaR[0],dcaZ[0]);
 343       track2->GetImpactParameters(dcaR[1],dcaZ[1]);
 344
 345       if(track2->Pt()<0.5) continue;
 346       Double_t rIsol = TMath::Sqrt( deltaPhi*deltaPhi+(track1->Eta()-track2->Eta())*(track1->Eta()-track2->Eta()) );
 347       fRisol->Fill(rIsol); //Fill R histogram
 348       if(track2->Pt()<fPtMin) continue;
 349
 350         fTheta->Fill(theta);
 351         fThetaZoom->Fill(theta);
 352         fThetaPt1Pt2->Fill(theta,track1->Pt(),track2->Pt());
 353         fThetaPt1Pt2Signed->Fill(theta,track1->GetSign()*track1->Pt(),track2->GetSign()*track2->Pt());
 354         fDeltaPhiSumEtaPt1->Fill(deltaPhi,track1->Eta()+track2->Eta(),track1->Pt());
 355         fDeltaPhiSumEtaPt2->Fill(deltaPhi,track1->Eta()+track2->Eta(),track2->Pt());
 356         fThetaDCAZ1DCAZ2->Fill(theta,dcaZ[0],dcaZ[1]);
 357         fRisolTheta->Fill(rIsol,theta);
 358         if(TMath::Abs(TMath::Pi()-theta)<fMaxCosmicAngle) {
 359           fDCAZCosmicCandidates->Fill(dcaZ[0],dcaZ[1]);
 360           fDCARCosmicCandidates->Fill(dcaR[0],dcaR[1]);
 361           fPtSignedCosmicCandidates->Fill(track1->GetSign()*track1->Pt());
 362         }
 363
 364     } // track2 loop
 365
 366   } // track1 loop
 367
 368   PostData(1,fHistListCosmics);
 369 }
 370
 371 //________________________________________________________________________
 372 void AliPWG4CosmicCandidates::Terminate(Option_t *)
 373 {
 374   //
 375   // Called once at the end of the query
 376   //
 377
 378 }