PWG1/cosmic/AliAnalysisTaskCosmic.cxx

   1 #include "TChain.h"
   2 #include "TTree.h"
   3 #include "TH1F.h"
   4 #include "TH2F.h"
   5 #include "TProfile.h"
   6 #include "TList.h"
   7 #include "TTree.h"
   8 #include "TParticle.h"
   9 #include "TParticlePDG.h"
  10
  11
  12 #include "AliESDEvent.h"
  13 #include "AliESDtrack.h"
  14 #include "AliAnalysisTaskCosmic.h"
  15 #include "AliESDInputHandler.h"
  16
  17 ClassImp(AliAnalysisTaskCosmic)
  18
  19 //________________________________________________________________________
  20 AliAnalysisTaskCosmic::AliAnalysisTaskCosmic(const char *name)
  21     : AliAnalysisTaskSE(name),
  22       fHists(0),
  23       fhDZvsZ(0),
  24       fhDZvsPhi(0),
  25       fhCh1Ch2(0),
  26       fhPh1Ph2(0),
  27       fhCl1Cl2G(0),
  28       fhCl1Cl2B(0)
  29 {
  30   //
  31   // Constructor
  32   //
  33     for (Int_t i = 0; i < 6; i++) {
  34         fhPt[i]      = 0;
  35         fhTheta[i]   = 0;
  36         fhPhi[i]     = 0;
  37         fhDPhi[i]    = 0;
  38         fhDTheta[i]  = 0;
  39         fhDZ[i]      = 0;
  40         fhDX[i]      = 0;
  41         fhDY[i]      = 0;
  42         fhDPt[i]     = 0;
  43         fhD1ovPt[i]  = 0;
  44         fpDPt[i]     = 0;
  45         fhDPtovPt[i] = 0;
  46         fpDPtS[i]    = 0;
  47     }
  48
  49   DefineOutput(1,  TList::Class());
  50 }
  51
  52
  53 //________________________________________________________________________
  54 void AliAnalysisTaskCosmic::UserCreateOutputObjects()
  55 {
  56   // Create histograms
  57   // Called once
  58     TString ext[6] = {"PC", "NC", "PZ", "NZ", "_Good", "_Bad"};
  59
  60     char name[12];
  61
  62     fHists = new TList();
  63
  64     for (Int_t i = 0; i < 6; i++) {
  65         // Pt
  66         sprintf(name, "fhPt%2s", ext[i].Data());
  67         fhPt[i]   = new TH1F(name,   " pT distribution",        800, 0., 200.);
  68         fhPt[i]->SetXTitle("p_{T} [Gev]");
  69         // Phi
  70         sprintf(name, "fhPhi%2s", ext[i].Data());
  71         fhPhi[i]   = new TH1F(name,   "Phi distribution",        62,  0., 2. * TMath::Pi());
  72         fhPhi[i]->SetXTitle("#phi [rad]");
  73         // Theta
  74         sprintf(name, "fhTheta%2s", ext[i].Data());
  75         fhTheta[i]   = new TH1F(name,   "Theta distribution",    62,  0., TMath::Pi());
  76         fhTheta[i]->SetXTitle("#theta [rad]");
  77         // Delta Phi
  78         sprintf(name, "fhDPhi%2s", ext[i].Data());
  79         fhDPhi[i]   = new TH1F(name,   "DeltaPhi distribution",        320, -0.4, 0.4);
  80         fhDPhi[i]->SetXTitle("#Delta#phi [rad]");
  81         // Delta Theta
  82         sprintf(name, "fhDTheta%2s", ext[i].Data());
  83         fhDTheta[i]   = new TH1F(name,   "DeltaTheta distribution",    320, -0.4, 0.4);
  84         fhDTheta[i]->SetXTitle("#Delta#theta [rad]");
  85         // Delta Z
  86         sprintf(name, "fhDZ%2s", ext[i].Data());
  87         fhDZ[i]   = new TH1F(name,   "DeltaZ distribution",        200, -10., 10.);
  88         fhDZ[i]->SetXTitle("#DeltaZ [cm]");
  89         // Delta X
  90         sprintf(name, "fhDX%2s", ext[i].Data());
  91         fhDX[i]   = new TH1F(name,   "DeltaX distribution",        200, -10., 10.);
  92         fhDX[i]->SetXTitle("#DeltaX [cm]");
  93         // Delta Y
  94         sprintf(name, "fhDY%2s", ext[i].Data());
  95         fhDY[i]   = new TH1F(name,   "DeltaY distribution",        200, -10, 10.);
  96         fhDY[i]->SetXTitle("#DeltaY [cm]");
  97         // Delta Pt
  98         sprintf(name, "fhDPt%2s", ext[i].Data());
  99         fhDPt[i]   = new TH1F(name,   "DeltaPt distribution",        200, -20., 20.);
 100         fhDPt[i]->SetXTitle("#Delta p_{T}  [GeV]");
 101
 102         // Delta 1/Pt
 103         sprintf(name, "fhD1ovPt%2s", ext[i].Data());
 104         fhD1ovPt[i]   = new TH1F(name,   "Delta 1/Pt distribution",        200, -1., 1.);
 105         fhD1ovPt[i]->SetXTitle("#Delta 1/Pt");
 106
 107         // Delta Pt over Pt
 108         sprintf(name, "fhDPtovPt%2s", ext[i].Data());
 109         fhDPtovPt[i]   = new TH1F(name,   "DeltaPt/Pt distribution",        200, -2., 2.);
 110         fhDPtovPt[i]->SetXTitle("#DeltaPt/Pt");
 111
 112
 113
 114         // Delta Pt/ Pt vs Pt
 115         sprintf(name, "fpDPt%2s", ext[i].Data());
 116         fpDPt[i] = new TProfile(name, "#Delta Pt / Pt", 20, 0., 20., -1, 1., "S");
 117         fpDPt[i]->SetXTitle("p_{T} [GeV]");
 118         fpDPt[i]->SetYTitle("#Delta 1/p_{T} [GeV^{-1}]");
 119         // Delta Pt error
 120         sprintf(name, "fpDPtS%2s", ext[i].Data());
 121         fpDPtS[i] = new TProfile(name, "#Delta Pt / <sigma>", 20, 0., 20., 0., 10.);
 122         fpDPtS[i]->SetXTitle("p_{T}");
 123         fpDPtS[i]->SetYTitle("#Delta p_{T} / <#sigma_{p_{T}}>");
 124
 125
 126         fHists->Add(fhPt[i]);
 127         fHists->Add(fhPhi[i]);
 128         fHists->Add(fhTheta[i]);
 129         fHists->Add(fhDPhi[i]);
 130         fHists->Add(fhDPt[i]);
 131         fHists->Add(fhD1ovPt[i]);
 132         fHists->Add(fhDPtovPt[i]);
 133         fHists->Add(fhDTheta[i]);
 134         fHists->Add(fhDZ[i]);
 135         fHists->Add(fhDX[i]);
 136         fHists->Add(fhDY[i]);
 137         fHists->Add(fpDPt[i]);
 138         fHists->Add(fpDPtS[i]);
 139     }
 140
 141     fhDZvsZ      = new TH2F("fhDZvsZ",    "dz vs z", 500, -250., 250., 100, -10., 10.);
 142     fhDZvsZ->SetXTitle("z_{in} * sign(z_{in}) * sign(z_{out}) [cm]");
 143     fhDZvsZ->SetYTitle("#Delta z [cm]");
 144
 145
 146     fhDZvsPhi    = new TH2F("fhDZvsPhi", "dz vs phi", 36, 0., TMath::Pi(),  50, -2., 2.);
 147
 148     fhCh1Ch2   = new TH2F("fCh1Ch2",   "ch1 vs ch2", 8, -2., 2., 8, -2., 2.);
 149     fhPh1Ph2   = new TH2F("fPh1Ph2",   "ph1 vs ph2", 128, 0., 2. * TMath::Pi(), 128, 0., 2. * TMath::Pi());
 150     fhCl1Cl2G  = new TH2F("fCl1Cl2G",   "#cl vs #cl", 200, 0., 200., 200, 0., 200);
 151     fhCl1Cl2B  = new TH2F("fCl1Cl2B",   "#cl vs #cl", 200, 0., 200., 200, 0., 200);
 152
 153     fHists->Add(fhDZvsZ);
 154     fHists->Add(fhDZvsPhi);
 155     fHists->Add(fhCh1Ch2);
 156     fHists->Add(fhPh1Ph2);
 157     fHists->Add(fhCl1Cl2G);
 158     fHists->Add(fhCl1Cl2B);
 159     fHists->SetOwner();
 160
 161 }
 162
 163 //________________________________________________________________________
 164 void AliAnalysisTaskCosmic::UserExec(Option_t *)
 165 {
 166   // Main loop
 167   // Called for each event
 168
 169   if (!fInputEvent) {
 170     Printf("ERROR: fESD not available");
 171     return;
 172   }
 173
 174   AliESDEvent* esdE = (AliESDEvent*)fInputEvent;
 175   if (esdE->GetNumberOfTracks() != 2) return;
 176
 177   for (Int_t iTracks = 0; iTracks < esdE->GetNumberOfTracks(); iTracks++) {
 178       AliESDtrack* track = esdE->GetTrack(iTracks);
 179
 180
 181
 182 //      const AliExternalTrackParam * track = trackG->GetTPCInnerParam();
 183 //      if (!track) continue;
 184       // Pt cut
 185       Float_t pt     = track->Pt();
 186       Float_t charge = track->Charge();
 187       Float_t phi    = track->Phi();
 188       if (phi > 0. && phi  < TMath::Pi()) charge*=(-1.);
 189
 190
 191
 192       // TPC track
 193       UInt_t status = track->GetStatus();
 194       if (status&AliESDtrack::kTPCrefit ==0) continue;
 195
 196       Int_t nClustersTPC = track->GetTPCclusters(0);
 197       if (nClustersTPC < 50) continue;
 198
 199       //
 200
 201       Float_t z      = track->GetZ();
 202       Float_t x      = track->Xv();
 203       Float_t y      = track->Yv();
 204       Float_t theta  = track->Theta();
 205
 206
 207       const AliExternalTrackParam * trackOut = track->GetOuterParam();
 208       Float_t zOut = 0.;
 209       if (trackOut)zOut = trackOut->Zv();
 210
 211
 212       if (charge > 0) {
 213           fhPt[kPosC]   ->Fill(pt);
 214       } else {
 215           fhPt[kNegC]   ->Fill(pt);
 216       }
 217
 218       //      if ((TMath::Abs(esdE->GetCurrentL3()) > 1.e-3) && (pt < 1. || pt > 50.)) continue;
 219
 220
 221       //
 222       if (charge > 0) {
 223           fhPhi[kPosC]  ->Fill(phi);
 224           fhTheta[kPosC]->Fill(theta);
 225       } else {
 226           fhPhi[kNegC]  ->Fill(phi);
 227           fhTheta[kNegC]->Fill(theta);
 228       }
 229
 230
 231       if (z > 0) {
 232           fhPt[kPosZ]   ->Fill(pt);
 233           fhPhi[kPosZ]  ->Fill(phi);
 234           fhTheta[kPosZ]->Fill(theta);
 235       } else {
 236           fhPt[kNegZ]   ->Fill(pt);
 237           fhPhi[kNegZ]  ->Fill(phi);
 238           fhTheta[kNegZ]->Fill(theta);
 239       }
 240
 241
 242
 243       // Tracks coming from above
 244       if (phi > 0. &&  phi < TMath::Pi()) {
 245
 246 //        printf("Track#1 %5d %5d %13.3f %13.3f \n", Int_t(Entry()), iTracks, phi, zOut);
 247
 248           Float_t dphiMin = 999.;
 249           Float_t rMin    = 999.;
 250           Int_t   jMin    = -1;
 251           // Search for a matching track (in dphi-dtheta)
 252           for (Int_t jTracks = 0; jTracks < esdE->GetNumberOfTracks(); jTracks++) {
 253               if (jTracks == iTracks) continue;
 254
 255               AliESDtrack* track2 = esdE->GetTrack(jTracks);
 256
 257               UInt_t status = track2->GetStatus();
 258               if (status&AliESDtrack::kTPCrefit ==0) continue;
 259
 260               Int_t nClustersTPC2 = track2->GetTPCclusters(0);
 261               if (nClustersTPC2 < 50) continue;
 262               //              if ((TMath::Abs(esdE->GetCurrentL3()) > 1.e-3) && (track2->Pt() < 1. || track2->Pt() > 50.)) continue;
 263 //            const AliExternalTrackParam * track2 = trackG2->GetTPCInnerParam();
 264 //            if (!track2) continue;
 265
 266               Float_t phi2   = track2->Phi() - TMath::Pi();
 267               Float_t theta2 = TMath::Pi()   - track2->Theta();
 268               Float_t dphi   = phi2   - phi;
 269               Float_t dtheta = theta2 - theta;
 270
 271               if (dphi >  TMath::Pi()) dphi -= 2. * TMath::Pi();
 272               if (dphi < -TMath::Pi()) dphi += 2. * TMath::Pi();
 273
 274               Float_t dR = TMath::Sqrt(dphi * dphi + dtheta * dtheta);
 275
 276               if (dR < rMin) {
 277                   rMin    = dR;
 278                   dphiMin = dphi;
 279                   jMin = jTracks;
 280               }
 281
 282           } // tracks 2
 283           if (jMin != -1) {
 284               // we found a matching track candidate ...
 285               AliESDtrack* track2 = esdE->GetTrack(jMin);
 286               const AliExternalTrackParam * trackOut2 = track2->GetOuterParam();
 287               Float_t zOut2 = 0.;
 288               if (trackOut2) zOut2 = trackOut2->Zv();
 289
 290
 291               Float_t theta2 = - track2->Theta() + TMath::Pi();
 292               Float_t z2 = track2->GetZ();
 293               Float_t x2 = track2->Xv();
 294               Float_t y2 = track2->Yv();
 295               Float_t charge2 = track2->Charge();
 296               Float_t dz  = z2 - z;
 297               Float_t dx  = x2 - x;
 298               Float_t dy  = y2 - y;
 299               Float_t pt1 = track->Pt();
 300               Float_t pt2 = track2->Pt();
 301               Float_t dpt = pt2 - pt1;
 302               Float_t d1pt = 1./pt2 - 1./pt1;
 303
 304               Float_t ptm = 0.5 * (pt1 + pt2);
 305
 306               Int_t nClustersTPC2 = track2->GetTPCclusters(0);
 307
 308               if (charge > 0.) {
 309                   fhDPhi[kPosC]  ->Fill(dphiMin);
 310                   fhDTheta[kPosC]->Fill(theta2 - theta);
 311               } else {
 312                   fhDPhi[kNegC]  ->Fill(dphiMin);
 313                   fhDTheta[kNegC]->Fill(theta2 - theta);
 314               }
 315
 316               if (z > 0.) {
 317                   fhDPhi[kPosZ]  ->Fill(dphiMin);
 318                   fhDTheta[kPosZ]->Fill(theta2 - theta);
 319               } else {
 320                   fhDPhi[kNegZ]  ->Fill(dphiMin);
 321                   fhDTheta[kNegZ]->Fill(theta2 - theta);
 322               }
 323
 324               if (TMath::Abs(dpt)/ptm > 0.5) {
 325                   fhDPhi[kBad]  ->Fill(dphiMin);
 326                   fhDTheta[kBad]->Fill(theta2 - theta);
 327               } else {
 328                   fhDPhi[kGood]  ->Fill(dphiMin);
 329                   fhDTheta[kGood]->Fill(theta2 - theta);
 330               }
 331               // Good matches ...
 332 //            if (TMath::Abs(rMin < 0.04) && (charge == charge2) && TMath::Abs(dz) < 0.5 &&  TMath::Abs(dy) && TMath::Abs(dx) < 0.5 )
 333 //            if (TMath::Abs(rMin < 0.04) && (charge == charge2) && (zOut * zOut2) < 0.)
 334               if (TMath::Abs(rMin < 0.04) && (charge == charge2))
 335               {
 336
 337                   if (TMath::Abs(dpt)/ptm  < .1) {
 338                       fhCl1Cl2G->Fill(nClustersTPC, nClustersTPC2);
 339                   }
 340
 341                   if (TMath::Abs(dpt)/ptm  > .5) fhCl1Cl2B->Fill(nClustersTPC, nClustersTPC2);
 342                   fhPh1Ph2->Fill(track->Phi(), track2->Phi());
 343
 344
 345                   Double_t sigmaPt1 = TMath::Sqrt(track->GetSigma1Pt2());
 346                   Double_t sigmaPt2 = TMath::Sqrt(track->GetSigma1Pt2());
 347                   Double_t sigmaPt = 0.5 * TMath::Sqrt(sigmaPt1 * sigmaPt1 + sigmaPt2 * sigmaPt2);
 348                   if (TMath::Abs(dpt)/ptm > 0.2 && pt > 10. &&  charge > 0.)
 349 //                    printf("Track#2 %5d %5d %13.3f %13.3f %13.3f %13.3f\n", Entry(), jMin, track2->Phi(), dz, dpt, zOut2);
 350                   if (zOut != 0. && zOut2 != 0.) fhDZvsZ->Fill(TMath::Abs(zOut) * zOut/zOut2, dz);
 351                   if (zOut * zOut2 > 0. && zOut < 0) fhDZvsPhi->Fill(phi, dz);
 352
 353                   fhCh1Ch2->Fill(charge, track2->Charge());
 354
 355
 356                   if (charge > 0.) {
 357                       fhDPt[kPosC]->Fill(dpt);
 358                       fhD1ovPt[kPosC]->Fill(d1pt);
 359                       fpDPt[kPosC]->Fill(ptm, d1pt);
 360                       if (ptm > 5. && ptm < 15.) fhDPtovPt[kPosC]->Fill(dpt/ptm);
 361                       fhDZ[kPosC]->Fill(dz);
 362                       fhDX[kPosC]->Fill(dx);
 363                       fhDY[kPosC]->Fill(dy);
 364                       fpDPtS[kPosC]->Fill(ptm, 2. * TMath::Abs(1./pt1 - 1./pt2)/sigmaPt);
 365                   } else {
 366                       fhDPt[kNegC]->Fill(dpt);
 367                       fhD1ovPt[kNegC]->Fill(d1pt);
 368                       fpDPt[kNegC]->Fill(ptm, d1pt);
 369                       if (ptm > 5. && ptm < 15.) fhDPtovPt[kNegC]->Fill(dpt/ptm);
 370                       fhDZ[kNegC]->Fill(dz);
 371                       fhDX[kNegC]->Fill(dx);
 372                       fhDY[kNegC]->Fill(dy);
 373                       fpDPtS[kNegC]->Fill(ptm, 2. * TMath::Abs(1./pt1 - 1./pt2)/sigmaPt);
 374                   }
 375
 376                   if (z > 0.) {
 377                       fhDPt[kPosZ]->Fill(dpt);
 378                       fhD1ovPt[kPosZ]->Fill(d1pt);
 379                       fpDPt[kPosZ]->Fill(ptm, d1pt);
 380                       fhDZ[kPosZ]->Fill(dz);
 381                       fhDX[kPosZ]->Fill(dx);
 382                       fhDY[kPosZ]->Fill(dy);
 383                       fpDPtS[kPosZ]->Fill(ptm, 2. * TMath::Abs(1./pt1 - 1./pt2)/sigmaPt);
 384                   } else {
 385                       fhDPt[kNegZ]->Fill(dpt);
 386                       fhD1ovPt[kNegZ]->Fill(d1pt);
 387                       fpDPt[kNegZ]->Fill(ptm, d1pt);
 388                       fhDZ[kNegZ]->Fill(dz);
 389                       fhDX[kNegZ]->Fill(dx);
 390                       fhDY[kNegZ]->Fill(dy);
 391                       fpDPtS[kNegZ]->Fill(ptm, 2. * TMath::Abs(1./pt1 - 1./pt2)/sigmaPt);
 392                   }
 393
 394
 395                   if (TMath::Abs(dpt)/ptm > 0.5) {
 396                       fhDPt[kBad]->Fill(dpt);
 397                       fpDPt[kBad]->Fill(ptm, TMath::Abs(dpt)/ptm);
 398                       fhDZ[kBad]->Fill(dz);
 399                       fhDX[kBad]->Fill(dx);
 400                       fhDY[kBad]->Fill(dy);
 401                       fpDPtS[kBad]->Fill(ptm, 2. * TMath::Abs(1./pt1 - 1./pt2)/sigmaPt);
 402                   } else {
 403                       fhDPt[kGood]->Fill(dpt);
 404                       fpDPt[kGood]->Fill(ptm, TMath::Abs(dpt)/ptm);
 405                       fhDZ[kGood]->Fill(dz);
 406                       fhDX[kGood]->Fill(dx);
 407                       fhDY[kGood]->Fill(dy);
 408                       fpDPtS[kGood]->Fill(ptm, 2. * TMath::Abs(1./pt1 - 1./pt2)/sigmaPt);
 409                   }
 410
 411               } // good matches
 412           } // found possible match
 413       } // upper
 414   } // tracks 1
 415   PostData(1, fHists);
 416 }
 417
 418 //________________________________________________________________________
 419 void AliAnalysisTaskCosmic::Terminate(Option_t *)
 420 {
 421 }