PWGPP/forward/SPDClustTask/AliSPDClustTask.cxx

   1 /*************************************************************************
   2 * Copyright(c) 1998-2008, 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 // Class AliSPDClustTask                                            //
  18 // Analysis task to produce data and MC histos needed for tracklets      //
  19 // dNdEta extraction in multiple bins in one go                          //
  20 // Author:  ruben.shahoyan@cern.ch                                       //
  21 ///////////////////////////////////////////////////////////////////////////
  22 /*
  23   Important parameters to set:
  24   1) make sure to initialize correct geometry in UserCreateOutputObjects
  25   2) The cut on signal selection variable (delta, dphi ...) should be decided beforehand
  26 ...
  27 */
  28
  29 #include "TChain.h"
  30 #include "TTree.h"
  31 #include "TRandom.h"
  32 #include "TH1F.h"
  33 #include "TH2F.h"
  34 #include "TList.h"
  35 #include "TObjArray.h"
  36 #include "TGeoGlobalMagField.h"
  37
  38 #include "AliAnalysisManager.h"
  39
  40 #include "AliMultiplicity.h"
  41 #include "AliESDEvent.h"
  42 #include "AliESDInputHandler.h"
  43 #include "AliESDInputHandlerRP.h"
  44 #include "AliCDBPath.h"
  45 #include "AliCDBManager.h"
  46 #include "AliCDBEntry.h"
  47 #include "AliCDBStorage.h"
  48 #include "AliGeomManager.h"
  49 #include "AliMagF.h"
  50 #include "AliESDVZERO.h"
  51 #include "AliESDZDC.h"
  52 #include "AliRunLoader.h"
  53 #include "AliMCEventHandler.h"
  54 #include "AliMCEvent.h"
  55 #include "AliMCParticle.h"
  56 #include "AliStack.h"
  57 #include "AliGenEventHeader.h"
  58 #include "AliCentrality.h"
  59 #include "../ITS/AliITSRecPoint.h"
  60 #include "../ITS/AliITSgeomTGeo.h"
  61 #include "../ITS/AliITSMultReconstructor.h"
  62
  63 #include "AliLog.h"
  64
  65 #include "AliPhysicsSelection.h"
  66 #include "AliSPDClustTask.h"
  67 #include "AliITSMultReconstructor.h"
  68 #include "AliGenEventHeader.h"
  69 #include "AliGenHijingEventHeader.h"
  70 #include "AliGenDPMjetEventHeader.h"
  71 #include "AliESDtrackCuts.h"
  72
  73 ClassImp(AliSPDClustTask)
  74
  75
  76 //________________________________________________________________________
  77 /*//Default constructor
  78 AliSPDClustTask::AliSPDClustTask(const char *name)
  79   : AliAnalysisTaskSE(name),
  80 */
  81 //________________________________________________________________________
  82 AliSPDClustTask::AliSPDClustTask(const char *name)
  83   : AliAnalysisTaskSE(name),
  84 //
  85   fOutput(0),
  86   fUseMC(kFALSE),
  87   fHistos(0),
  88 //
  89   fEtaMin(-3.0),
  90   fEtaMax(3.0),
  91   fZVertexMin(-20),
  92   fZVertexMax( 20),
  93 //
  94   fScaleDTBySin2T(kFALSE),
  95   fCutOnDThetaX(kFALSE),
  96   fNStdDev(1.),
  97   fDPhiWindow(0.08),
  98   fDThetaWindow(0.025),
  99   fDPhiShift(0.0045),
 100   fPhiOverlapCut(0.005),
 101   fZetaOverlap(0.05),
 102   fRemoveOverlaps(kFALSE),
 103 //
 104   fDPhiSCut(0.06),
 105   fNStdCut(1.),
 106 //
 107   fMultReco(0),
 108   fRPTree(0),
 109   fStack(0),
 110   fMCEvent(0),
 111   //
 112   fDontMerge(kFALSE),
 113   fhPtIn(0)
 114 {
 115   // Constructor
 116
 117   DefineOutput(1, TList::Class());
 118   //
 119   SetScaleDThetaBySin2T();
 120   SetNStdDev();
 121   SetPhiWindow();
 122   SetThetaWindow();
 123   SetPhiShift();
 124   SetPhiOverlapCut();
 125   SetZetaOverlapCut();
 126   SetRemoveOverlaps();
 127   //
 128 }
 129
 130 //________________________________________________________________________
 131 AliSPDClustTask::~AliSPDClustTask()
 132 {
 133   // Destructor
 134   // histograms are in the output list and deleted when the output
 135   // list is deleted by the TSelector dtor
 136   if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {  //RRR
 137     printf("Deleteing output\n");
 138     delete fOutput;
 139     fOutput = 0;
 140   }
 141   //
 142   delete fMultReco;
 143   delete fHistos;
 144   if (fhPtIn) delete fhPtIn;
 145   //
 146 }
 147
 148 //________________________________________________________________________
 149 void AliSPDClustTask::UserCreateOutputObjects()
 150 {
 151   //
 152   //  OpenFile(1);  fDontMerge = kTRUE;
 153   //
 154
 155   fOutput = new TList();
 156   fOutput->SetOwner();
 157   //
 158   AliCDBManager *man = AliCDBManager::Instance();
 159   if (fUseMC) {
 160     Bool_t newGeom = kTRUE;
 161     man->SetDefaultStorage("alien://Folder=/alice/simulation/2008/v4-15-Release/Residual");
 162     if (newGeom) {
 163       // new geom
 164       AliCDBEntry*  obj = man->Get("GRP/Geometry/Data",130844,8);
 165       AliGeomManager::SetGeometry((TGeoManager*) obj->GetObject());
 166       if (!AliGeomManager::ApplyAlignObjsToGeom("ITS",130844,6,-1)) AliFatal("Failed to misalign geometry");
 167     }
 168     else {
 169       // old geom
 170       AliCDBEntry*  obj = man->Get("GRP/Geometry/Data",130845,7);
 171       AliGeomManager::SetGeometry((TGeoManager*) obj->GetObject());
 172       if (!AliGeomManager::ApplyAlignObjsToGeom("ITS",130845,5,-1)) AliFatal("Failed to misalign geometry");
 173     }
 174   }
 175   else {
 176     man->SetDefaultStorage("alien://Folder=/alice/data/2010/OCDB"); //man->SetRun(137366);
 177     AliCDBEntry*  obj = man->Get("GRP/Geometry/Data",137366);
 178     AliGeomManager::SetGeometry((TGeoManager*) obj->GetObject());
 179     if (!AliGeomManager::ApplyAlignObjsToGeom("ITS",137366,-1,-1)) AliFatal("Failed to misalign geometry");
 180   }
 181   //
 182   // Create histograms
 183   fHistos = BookHistos();
 184   for(Int_t i=0; i<fHistos->GetEntriesFast(); i++) {
 185     if(fHistos->At(i)){
 186       fOutput->AddLast(fHistos->At(i));
 187     }
 188   }
 189
 190   PostData(1, fOutput);
 191   //
 192 }
 193
 194 //________________________________________________________________________
 195 void AliSPDClustTask::UserExec(Option_t *)
 196 {
 197   // Main loop
 198   //
 199   //
 200   AliAnalysisManager* anMan = AliAnalysisManager::GetAnalysisManager();
 201   fRPTree = 0;
 202   AliESDInputHandlerRP *handRP = (AliESDInputHandlerRP*)anMan->GetInputEventHandler();
 203   if (!handRP) { AliError("No RP handler"); return; }
 204   //
 205   fRPTree = handRP->GetTreeR("ITS");
 206   if (!fRPTree) { AliError(" Invalid ITS cluster tree !\n"); return; }
 207   //
 208   AliESDEvent *esd  = handRP->GetEvent();
 209   if (!esd) { AliError("No AliESDEvent"); return; }
 210   //
 211   AliMCEventHandler* eventHandler = 0;
 212   fMCEvent = 0;
 213   fStack = 0;
 214   //
 215   if (fUseMC) {
 216     eventHandler = (AliMCEventHandler*)anMan->GetMCtruthEventHandler();
 217     if (!eventHandler) { AliError("Could not retrieve MC event handler"); return; }
 218     fMCEvent = eventHandler->MCEvent();
 219     if (!fMCEvent) { AliError("Could not retrieve MC event"); return; }
 220     fStack = fMCEvent->Stack();
 221     if (!fStack) { AliError("Stack not available"); return; }
 222   }
 223   //
 224   const AliESDVertex* vtxESD = esd->GetPrimaryVertexSPD();
 225   if (vtxESD->GetNContributors()<1) return;
 226   //
 227   // ------------------RS: Put here additional event selection if needed ------->>
 228   // at the moment, I am just selecting on vertex
 229   TString vtxTyp = vtxESD->GetTitle();
 230   if (vtxTyp.Contains("vertexer: Z")) {
 231     if (vtxESD->GetDispersion()>0.04) return;
 232     if (vtxESD->GetZRes()>0.25) return;
 233   }
 234   //
 235   if (vtxESD->GetZ()<fZVertexMin || vtxESD->GetZ()>fZVertexMax) return;
 236   //
 237   // pile-up rejection
 238   if (esd->IsPileupFromSPD(3,0.8)) return;
 239   //
 240   // ------------------RS: Put here additional event selection if needed -------<<
 241   //
 242   fVtx[0] = vtxESD->GetX();
 243   fVtx[1] = vtxESD->GetY();
 244   fVtx[2] = vtxESD->GetZ();
 245   //
 246   // do we need to initialize the field?
 247   AliMagF* field = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
 248   if (!field && !esd->InitMagneticField()) {printf("Failed to initialize the B field\n");return;}
 249   //
 250   AliStack *stack = NULL;
 251   if (fUseMC) {
 252     if (!fMCEvent) { AliError("Could not retrieve MC event"); return; }
 253     stack = fMCEvent->Stack();
 254     if (!stack) { AliError("ERROR: Could not get stack"); return; }
 255     for (Int_t iPart = 0; iPart < stack->GetNtrack(); ++iPart) {
 256       if (!stack->IsPhysicalPrimary(iPart)) continue;
 257       TParticle* p = stack->Particle(iPart);
 258       if (TMath::Abs(p->Y()) > 0.5) continue;
 259       if (TMath::Abs(p->GetPdgCode()) != 211) continue;
 260       TH1D *hPt = (TH1D*)fHistos->At(kHPt);
 261       hPt->Fill(p->Pt());
 262     }
 263   }
 264   //
 265   InitMultReco();
 266   fMultReco->Reconstruct(fRPTree, fVtx);
 267   //  AliMultiplicity* mlt = fMultReco->GetMultiplicity();
 268   FillHistos(stack);
 269   //
 270   delete fMultReco;
 271   fMultReco = 0;
 272   //
 273 }
 274
 275 //________________________________________________________________________
 276 void AliSPDClustTask::Terminate(Option_t *)
 277 {
 278   Printf("Terminating...");
 279 }
 280
 281 //_________________________________________________________________________
 282 void AliSPDClustTask::InitMultReco()
 283 {
 284   // create mult reconstructor
 285   if (fMultReco) delete fMultReco;
 286   fMultReco = new AliITSMultReconstructor();
 287   fMultReco->SetCreateClustersCopy(kTRUE);
 288   fMultReco->SetScaleDThetaBySin2T(fScaleDTBySin2T);
 289   fMultReco->SetNStdDev(fNStdDev);
 290   fMultReco->SetPhiWindow( fDPhiWindow );
 291   fMultReco->SetThetaWindow( fDThetaWindow );
 292   fMultReco->SetPhiShift( fDPhiShift );
 293   fMultReco->SetRemoveClustersFromOverlaps(fRemoveOverlaps);
 294   fMultReco->SetPhiOverlapCut(fPhiOverlapCut);
 295   fMultReco->SetZetaOverlapCut(fZetaOverlap);
 296   fMultReco->SetHistOn(kFALSE);
 297   //
 298 }
 299
 300 //_________________________________________________________________________
 301 TObjArray* AliSPDClustTask::BookHistos()
 302 {
 303   // book set of histos
 304   //
 305   TObjArray* histos = new TObjArray();
 306   histos->SetOwner(kFALSE);
 307   // this array contains histos stored at slots:
 308   // 0-99    : tracklet related histos
 309   // 100-199 : SPD1 clusters not used by tracklets
 310   // 200-299 : SPD2 clusters not used by tracklets
 311   // 300-399 : SPD1 clusters used by tracklets
 312   // 400-499 : SPD2 clusters used by tracklets
 313   //
 314   int nEtaBins = int((fEtaMax-fEtaMin)/0.1);
 315   if (nEtaBins<1) nEtaBins = 1;
 316   //
 317   // just an example: cluster type vs eta
 318   for (int iused=0;iused<2;iused++) {
 319     for (int spd=0;spd<2;spd++) {
 320       TH2F* h = new TH2F(Form("clType_SPD%d_%s",spd,iused ? "used":"free"),
 321                          Form("clType SPD%d %s",spd,iused ? "used":"free"),
 322                          nEtaBins, fEtaMin,fEtaMax,
 323                          15, -0.5, 14.5);
 324       histos->AddAtAndExpand(h, kHClusters+iused*200+spd*100 + kClTypevsEta);
 325       TH1F* hZ = new TH1F(Form("clZ_SPD%d_%s",spd,iused ? "used":"free"),
 326                           Form("clZ SPD%d %s",spd,iused ? "used":"free"),
 327                           200, -15, 15);
 328       histos->AddAtAndExpand(hZ, kHClusters+iused*200+spd*100 + kClZ);
 329       TH1F* hEta = new TH1F(Form("clEta_SPD%d_%s",spd,iused ? "used":"free"),
 330                             Form("clEta SPD%d %s",spd,iused ? "used":"free"),
 331                             nEtaBins, fEtaMin,fEtaMax);
 332       histos->AddAtAndExpand(hEta, kHClusters+iused*200+spd*100 + kClEta);
 333
 334       hZ = new TH1F(Form("clZ_pions_weighted_SPD%d_%s",spd,iused ? "used":"free"),
 335                     Form("clZ pions weighted SPD%d %s",spd,iused ? "used":"free"),
 336                     200, -15, 15);
 337       histos->AddAtAndExpand(hZ, kHClusters+iused*200+spd*100 + kClZPionsW);
 338       hEta = new TH1F(Form("clEta_pions_weighted_SPD%d_%s",spd,iused ? "used":"free"),
 339                       Form("clEta pions weighted SPD%d %s",spd,iused ? "used":"free"),
 340                       nEtaBins, fEtaMin,fEtaMax);
 341       histos->AddAtAndExpand(hEta, kHClusters+iused*200+spd*100 + kClEtaPionsW);
 342
 343       hZ = new TH1F(Form("clZ_pions_SPD%d_%s",spd,iused ? "used":"free"),
 344                     Form("clZ pions SPD%d %s",spd,iused ? "used":"free"),
 345                     200, -15, 15);
 346       histos->AddAtAndExpand(hZ, kHClusters+iused*200+spd*100 + kClZPions);
 347       hEta = new TH1F(Form("clEta_pions_SPD%d_%s",spd,iused ? "used":"free"),
 348                       Form("clEta pions SPD%d %s",spd,iused ? "used":"free"),
 349                       nEtaBins, fEtaMin,fEtaMax);
 350       histos->AddAtAndExpand(hEta, kHClusters+iused*200+spd*100 + kClEtaPions);
 351     }
 352   }
 353   if (fUseMC) {
 354     TH1D *hPt = new TH1D("hPt","Pions Pt spectra (MC)",fhPtIn->GetXaxis()->GetNbins(),fhPtIn->GetXaxis()->GetXbins()->GetArray());
 355     histos->AddAtAndExpand(hPt, kHPt);
 356   }
 357   //
 358   return histos;
 359 }
 360
 361 //_________________________________________________________________________
 362 void AliSPDClustTask::FillHistos(AliStack *stack)
 363 {
 364   // fill info on clusters, separately on associated to tracklets and singles
 365   const int kUsedBit = BIT(15);
 366   //
 367   // 0) --------- just in case: reset kUsedBit of clusters --->>>
 368   for (int ilr=0;ilr<2;ilr++) {
 369     for (int icl=fMultReco->GetNClustersLayer(ilr);icl--;) {  // loop on clusters of layer
 370       AliITSRecPoint *clus = fMultReco->GetRecPoint(ilr,icl);
 371       if (!clus) {
 372         AliError(Form("Failed to extract cluster %d of layer %d",icl,ilr));
 373         return;
 374       }
 375       clus->ResetBit(kUsedBit);
 376     } // loop on clusters of layer
 377   } // loop on layers
 378   //
 379   // 1) --------- mark used clusters by kUsedBit ------------>>>
 380   int ntr = fMultReco->GetNTracklets();
 381   for (int itr=ntr;itr--;) {
 382     Float_t *trc = fMultReco->GetTracklet(itr);
 383     AliITSRecPoint *cl0 = fMultReco->GetRecPoint(0,(int)trc[AliITSMultReconstructor::kClID1]); // cluster on SPD1
 384     AliITSRecPoint *cl1 = fMultReco->GetRecPoint(1,(int)trc[AliITSMultReconstructor::kClID2]); // cluster on SPD1
 385     //
 386     if (!cl0 || !cl1) {
 387       AliError(Form("Failed to extract clusters associated with tracklet %d",itr));
 388       continue;
 389     }
 390     cl0->SetBit(kUsedBit);
 391     cl1->SetBit(kUsedBit);
 392   }
 393   //
 394   // 2) --------- fill needed info for used and unuses clusters
 395   TH2F* h2d = 0;
 396   TH1F* h1d = 0;
 397   //
 398   for (int ilr=0;ilr<2;ilr++) { // loop on layers
 399     int ncl = fMultReco->GetNClustersLayer(ilr);
 400     for (int icl=ncl;icl--;) {  // loop on clusters of layer
 401       AliITSRecPoint *clus = fMultReco->GetRecPoint(ilr,icl); // extract the cluster to access clType
 402       Float_t* clInfo = fMultReco->GetClusterOfLayer(ilr,icl); // we already have in the MultReco cluster params extracted and
 403       // stored in the float array, see AliITSMultReconstructor
 404       //
 405       float phi = clInfo[AliITSMultReconstructor::kClPh];
 406       float eta = -TMath::Log(TMath::Tan(clInfo[AliITSMultReconstructor::kClTh]/2));
 407       float z   = clInfo[AliITSMultReconstructor::kClZ];
 408       //
 409       int used = clus->TestBit(kUsedBit) ? 1:0;
 410       //
 411       // Annalisa, here you should fill the info on used/unused clusters
 412       if (phi > 0 && phi < 1.6) {
 413         h2d = (TH2F*)fHistos->At(kHClusters+used*200+ilr*100 + kClTypevsEta);
 414         h2d->Fill(eta, clus->GetClusterType());
 415         h1d = (TH1F*)fHistos->At(kHClusters+used*200+ilr*100 + kClZ);
 416         h1d->Fill(z);
 417         h1d = (TH1F*)fHistos->At(kHClusters+used*200+ilr*100 + kClEta);
 418         h1d->Fill(eta);
 419         if (stack) {
 420           for(Int_t iLabel = 0; iLabel < 3; ++iLabel) {
 421             Int_t clLabel = clus->GetLabel(iLabel);
 422             if (clLabel <= 0) continue;
 423             Int_t moLabel = FindMotherParticle(stack, clLabel);
 424             if (moLabel <= 0) continue;
 425             TParticle* p = stack->Particle(moLabel);
 426             if (TMath::Abs(p->GetPdgCode()) != 211) continue;
 427             Double_t weight = PtWeight(p->Pt());
 428             h1d = (TH1F*)fHistos->At(kHClusters+used*200+ilr*100 + kClZPionsW);
 429             h1d->Fill(z,weight);
 430             h1d = (TH1F*)fHistos->At(kHClusters+used*200+ilr*100 + kClEtaPionsW);
 431             h1d->Fill(eta,weight);
 432             h1d = (TH1F*)fHistos->At(kHClusters+used*200+ilr*100 + kClZPions);
 433             h1d->Fill(z);
 434             h1d = (TH1F*)fHistos->At(kHClusters+used*200+ilr*100 + kClEtaPions);
 435             h1d->Fill(eta);
 436             break;
 437           }
 438         }
 439       }
 440     } // loop on clusters of layer
 441   } // loop on layers
 442   //
 443 }
 444
 445 //________________________________________________________________________
 446 void AliSPDClustTask::SetInput(const char *filename)
 447 {
 448   // Read pions Pt spectra
 449   // from a local input file
 450
 451   AliInfo(Form("Reading input histograms from %s",filename));
 452   if (fUseMC) {
 453     TFile *fPt = TFile::Open("spectraCombine.root");
 454     if (!fPt) { AliError("Failed to open input file"); return; }
 455     TList *lPt = (TList*)fPt->Get("output");
 456     fhPtIn = (TH1D*)lPt->FindObject("PionsPos_MB_combine_sum")->Clone("fhPtIn");
 457     fhPtIn->SetDirectory(0);
 458     fPt->Close();
 459   }
 460 }
 461
 462 //________________________________________________________________________
 463 Int_t AliSPDClustTask::FindMotherParticle(AliStack* stack, Int_t i)
 464 {
 465   if(stack->IsPhysicalPrimary(i)) return i;
 466   TParticle* p = stack->Particle(i);
 467   Int_t imo =  p->GetFirstMother();
 468
 469   //  printf("imo = %d\n",imo);
 470
 471   if(imo<=0)
 472     {
 473       AliInfo("particle is not a PhysPrim and has no mother");
 474       return imo;
 475     }
 476   return FindMotherParticle(stack, imo);
 477
 478 }
 479
 480 //________________________________________________________________________
 481 Double_t AliSPDClustTask::PtWeight(Double_t pt)
 482 {
 483   if (pt > 0.35)
 484     return 1.;
 485   else
 486     return (1.71-0.71*pt/0.35);
 487 }