cc2e4af8d0aa7fc4681a0c3ae70a41fd8f047939
[u/mrichter/AliRoot.git] / 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   fhPtPionIn(0),
114   fhPtKaonIn(0),
115   fhPtProtonIn(0)
116 {
117   // Constructor
118
119   DefineOutput(1, TList::Class());
120   //
121   SetScaleDThetaBySin2T();
122   SetNStdDev();
123   SetPhiWindow();
124   SetThetaWindow();
125   SetPhiShift();
126   SetPhiOverlapCut();
127   SetZetaOverlapCut();
128   SetRemoveOverlaps();
129   //
130 }
131
132 //________________________________________________________________________
133 AliSPDClustTask::~AliSPDClustTask()
134 {
135   // Destructor
136   // histograms are in the output list and deleted when the output
137   // list is deleted by the TSelector dtor
138   if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {  //RRR
139     printf("Deleteing output\n");
140     delete fOutput;
141     fOutput = 0;
142   }
143   //
144   delete fMultReco;
145   delete fHistos;
146   if (fhPtPionIn) delete fhPtPionIn;
147   if (fhPtKaonIn) delete fhPtKaonIn;
148   if (fhPtProtonIn) delete fhPtProtonIn;
149   //
150 }
151
152 //________________________________________________________________________
153 void AliSPDClustTask::UserCreateOutputObjects() 
154 {
155   //
156   //  OpenFile(1);  fDontMerge = kTRUE;
157   //  
158
159   fOutput = new TList();
160   fOutput->SetOwner(); 
161   //
162   AliCDBManager *man = AliCDBManager::Instance();
163   if (fUseMC) {
164     Bool_t newGeom = kTRUE;
165     man->SetDefaultStorage("alien://Folder=/alice/simulation/2008/v4-15-Release/Residual");
166     if (newGeom) {
167       // new geom
168       AliCDBEntry*  obj = man->Get("GRP/Geometry/Data",130844,8);
169       AliGeomManager::SetGeometry((TGeoManager*) obj->GetObject());
170       if (!AliGeomManager::ApplyAlignObjsToGeom("ITS",130844,6,-1)) AliFatal("Failed to misalign geometry");
171     }
172     else {
173       // old geom
174       AliCDBEntry*  obj = man->Get("GRP/Geometry/Data",130845,7);
175       AliGeomManager::SetGeometry((TGeoManager*) obj->GetObject());
176       if (!AliGeomManager::ApplyAlignObjsToGeom("ITS",130845,5,-1)) AliFatal("Failed to misalign geometry");
177     }
178   }
179   else {
180     man->SetDefaultStorage("alien://Folder=/alice/data/2010/OCDB"); //man->SetRun(137366);
181     AliCDBEntry*  obj = man->Get("GRP/Geometry/Data",137366);
182     AliGeomManager::SetGeometry((TGeoManager*) obj->GetObject());
183     if (!AliGeomManager::ApplyAlignObjsToGeom("ITS",137366,-1,-1)) AliFatal("Failed to misalign geometry");
184   }
185   //
186   // Create histograms
187   fHistos = BookHistos();
188   for(Int_t i=0; i<fHistos->GetEntriesFast(); i++) {
189     if(fHistos->At(i)){
190       fOutput->AddLast(fHistos->At(i));
191     }
192   }
193
194   PostData(1, fOutput);
195   //
196 }
197
198 //________________________________________________________________________
199 void AliSPDClustTask::UserExec(Option_t *) 
200 {
201   // Main loop
202   //
203   //
204   AliAnalysisManager* anMan = AliAnalysisManager::GetAnalysisManager();
205   fRPTree = 0;
206   AliESDInputHandlerRP *handRP = (AliESDInputHandlerRP*)anMan->GetInputEventHandler(); 
207   if (!handRP) { AliError("No RP handler"); return; }
208   //
209   fRPTree = handRP->GetTreeR("ITS");
210   if (!fRPTree) { AliError(" Invalid ITS cluster tree !\n"); return; }
211   //
212   AliESDEvent *esd  = handRP->GetEvent();
213   if (!esd) { AliError("No AliESDEvent"); return; }
214   //
215   AliMCEventHandler* eventHandler = 0;
216   fMCEvent = 0;
217   fStack = 0;
218   //
219   if (fUseMC) {
220     eventHandler = (AliMCEventHandler*)anMan->GetMCtruthEventHandler();
221     if (!eventHandler) { AliError("Could not retrieve MC event handler"); return; }
222     fMCEvent = eventHandler->MCEvent();
223     if (!fMCEvent) { AliError("Could not retrieve MC event"); return; }
224     fStack = fMCEvent->Stack();
225     if (!fStack) { AliError("Stack not available"); return; }
226   }
227   //
228   const AliESDVertex* vtxESD = esd->GetPrimaryVertexSPD();
229   if (vtxESD->GetNContributors()<1) return;
230   //
231   // ------------------RS: Put here additional event selection if needed ------->>
232   // at the moment, I am just selecting on vertex
233   TString vtxTyp = vtxESD->GetTitle();
234   if (vtxTyp.Contains("vertexer: Z")) {
235     if (vtxESD->GetDispersion()>0.04) return;
236     if (vtxESD->GetZRes()>0.25) return;
237   }
238   //
239   if (vtxESD->GetZ()<fZVertexMin || vtxESD->GetZ()>fZVertexMax) return;
240   //
241   // pile-up rejection
242   if (esd->IsPileupFromSPD(3,0.8)) return;
243   //
244   // ------------------RS: Put here additional event selection if needed -------<<
245   //
246   fVtx[0] = vtxESD->GetX();
247   fVtx[1] = vtxESD->GetY();
248   fVtx[2] = vtxESD->GetZ();
249   //
250   // do we need to initialize the field?
251   AliMagF* field = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
252   if (!field && !esd->InitMagneticField()) {printf("Failed to initialize the B field\n");return;}
253   //
254   AliStack *stack = NULL;
255   if (fUseMC) {
256     if (!fMCEvent) { AliError("Could not retrieve MC event"); return; }
257     stack = fMCEvent->Stack();
258     if (!stack) { AliError("ERROR: Could not get stack"); return; }
259     for (Int_t iPart = 0; iPart < stack->GetNtrack(); ++iPart) {
260       if (!stack->IsPhysicalPrimary(iPart)) continue;
261       TParticle* p = stack->Particle(iPart);
262       if (TMath::Abs(p->Y()) > 0.5) continue;
263       if (TMath::Abs(p->GetPdgCode()) != 211) continue;
264       TH1D *hPt = (TH1D*)fHistos->At(kHPt);
265       hPt->Fill(p->Pt());
266     }
267   }
268   //
269   InitMultReco();
270   fMultReco->Reconstruct(fRPTree, fVtx);
271   //  AliMultiplicity* mlt = fMultReco->GetMultiplicity();
272   FillHistos(stack);
273   //
274   delete fMultReco; 
275   fMultReco = 0;
276   //
277 }      
278
279 //________________________________________________________________________
280 void AliSPDClustTask::Terminate(Option_t *) 
281 {
282   Printf("Terminating...");
283 }
284
285 //_________________________________________________________________________
286 void AliSPDClustTask::InitMultReco()
287 {
288   // create mult reconstructor
289   if (fMultReco) delete fMultReco;
290   fMultReco = new AliITSMultReconstructor();
291   fMultReco->SetCreateClustersCopy(kTRUE);
292   fMultReco->SetScaleDThetaBySin2T(fScaleDTBySin2T);
293   fMultReco->SetNStdDev(fNStdDev);
294   fMultReco->SetPhiWindow( fDPhiWindow );
295   fMultReco->SetThetaWindow( fDThetaWindow );
296   fMultReco->SetPhiShift( fDPhiShift );
297   fMultReco->SetRemoveClustersFromOverlaps(fRemoveOverlaps);
298   fMultReco->SetPhiOverlapCut(fPhiOverlapCut);
299   fMultReco->SetZetaOverlapCut(fZetaOverlap);
300   fMultReco->SetHistOn(kFALSE); 
301   //
302 }
303
304 //_________________________________________________________________________
305 TObjArray* AliSPDClustTask::BookHistos()
306 {
307   // book set of histos
308   //
309   TObjArray* histos = new TObjArray();
310   histos->SetOwner(kFALSE);
311   // this array contains histos stored at slots:
312   // 0-99    : tracklet related histos
313   // 100-199 : SPD1 clusters not used by tracklets
314   // 200-299 : SPD2 clusters not used by tracklets
315   // 300-399 : SPD1 clusters used by tracklets
316   // 400-499 : SPD2 clusters used by tracklets
317   //
318   int nEtaBins = int((fEtaMax-fEtaMin)/0.1);
319   if (nEtaBins<1) nEtaBins = 1;
320   //
321   // just an example: cluster type vs eta
322   for (int iused=0;iused<2;iused++) {
323     for (int spd=0;spd<2;spd++) {
324       TH2F* h = new TH2F(Form("clType_SPD%d_%s",spd,iused ? "used":"free"),
325                          Form("clType SPD%d %s",spd,iused ? "used":"free"),
326                          nEtaBins, fEtaMin,fEtaMax,
327                          15, -0.5, 14.5);
328       histos->AddAtAndExpand(h, kHClusters+iused*200+spd*100 + kClTypevsEta);
329       TH1F* hZ = new TH1F(Form("clZ_SPD%d_%s",spd,iused ? "used":"free"),
330                           Form("clZ SPD%d %s",spd,iused ? "used":"free"),
331                           200, -15, 15);
332       histos->AddAtAndExpand(hZ, kHClusters+iused*200+spd*100 + kClZ);
333       TH1F* hEta = new TH1F(Form("clEta_SPD%d_%s",spd,iused ? "used":"free"),
334                             Form("clEta SPD%d %s",spd,iused ? "used":"free"),
335                             nEtaBins, fEtaMin,fEtaMax);
336       histos->AddAtAndExpand(hEta, kHClusters+iused*200+spd*100 + kClEta);
337
338       hZ = new TH1F(Form("clZ_pions_weighted_SPD%d_%s",spd,iused ? "used":"free"),
339                     Form("clZ pions weighted SPD%d %s",spd,iused ? "used":"free"),
340                     200, -15, 15);
341       histos->AddAtAndExpand(hZ, kHClusters+iused*200+spd*100 + kClZPionsW);
342       hEta = new TH1F(Form("clEta_pions_weighted_SPD%d_%s",spd,iused ? "used":"free"),
343                       Form("clEta pions weighted SPD%d %s",spd,iused ? "used":"free"),
344                       nEtaBins, fEtaMin,fEtaMax);
345       histos->AddAtAndExpand(hEta, kHClusters+iused*200+spd*100 + kClEtaPionsW);
346
347       hZ = new TH1F(Form("clZ_pions_SPD%d_%s",spd,iused ? "used":"free"),
348                     Form("clZ pions SPD%d %s",spd,iused ? "used":"free"),
349                     200, -15, 15);
350       histos->AddAtAndExpand(hZ, kHClusters+iused*200+spd*100 + kClZPions);
351       hEta = new TH1F(Form("clEta_pions_SPD%d_%s",spd,iused ? "used":"free"),
352                       Form("clEta pions SPD%d %s",spd,iused ? "used":"free"),
353                       nEtaBins, fEtaMin,fEtaMax);
354       histos->AddAtAndExpand(hEta, kHClusters+iused*200+spd*100 + kClEtaPions);
355     }
356   }
357   if (fUseMC) {
358     TH1D *hPt = new TH1D("hPt","Pions Pt spectra (MC)",fhPtPionIn->GetXaxis()->GetNbins(),fhPtPionIn->GetXaxis()->GetXbins()->GetArray());
359     histos->AddAtAndExpand(hPt, kHPt);
360   }
361   //
362   return histos;
363 }
364
365 //_________________________________________________________________________
366 void AliSPDClustTask::FillHistos(AliStack *stack)
367 {
368   // fill info on clusters, separately on associated to tracklets and singles
369   const int kUsedBit = BIT(15);
370   //
371   // 0) --------- just in case: reset kUsedBit of clusters --->>>
372   for (int ilr=0;ilr<2;ilr++) {
373     for (int icl=fMultReco->GetNClustersLayer(ilr);icl--;) {  // loop on clusters of layer
374       AliITSRecPoint *clus = fMultReco->GetRecPoint(ilr,icl);
375       if (!clus) {
376         AliError(Form("Failed to extract cluster %d of layer %d",icl,ilr));
377         return;
378       }
379       clus->ResetBit(kUsedBit);
380     } // loop on clusters of layer
381   } // loop on layers
382   //
383   // 1) --------- mark used clusters by kUsedBit ------------>>>
384   int ntr = fMultReco->GetNTracklets();
385   for (int itr=ntr;itr--;) {
386     Float_t *trc = fMultReco->GetTracklet(itr);
387     AliITSRecPoint *cl0 = fMultReco->GetRecPoint(0,(int)trc[AliITSMultReconstructor::kClID1]); // cluster on SPD1
388     AliITSRecPoint *cl1 = fMultReco->GetRecPoint(1,(int)trc[AliITSMultReconstructor::kClID2]); // cluster on SPD1
389     //
390     if (!cl0 || !cl1) {
391       AliError(Form("Failed to extract clusters associated with tracklet %d",itr));
392       continue;
393     }
394     cl0->SetBit(kUsedBit);
395     cl1->SetBit(kUsedBit);
396   }
397   //
398   // 2) --------- fill needed info for used and unuses clusters
399   TH2F* h2d = 0;
400   TH1F* h1d = 0;
401   //
402   for (int ilr=0;ilr<2;ilr++) { // loop on layers
403     int ncl = fMultReco->GetNClustersLayer(ilr);
404     for (int icl=ncl;icl--;) {  // loop on clusters of layer
405       AliITSRecPoint *clus = fMultReco->GetRecPoint(ilr,icl); // extract the cluster to access clType
406       Float_t* clInfo = fMultReco->GetClusterOfLayer(ilr,icl); // we already have in the MultReco cluster params extracted and
407       // stored in the float array, see AliITSMultReconstructor
408       //
409       float phi = clInfo[AliITSMultReconstructor::kClPh];
410       float eta = -TMath::Log(TMath::Tan(clInfo[AliITSMultReconstructor::kClTh]/2));
411       float z   = clInfo[AliITSMultReconstructor::kClZ];
412       //
413       int used = clus->TestBit(kUsedBit) ? 1:0;
414       //
415       // Annalisa, here you should fill the info on used/unused clusters
416       if (phi > 0 && phi < 1.6) {
417         h2d = (TH2F*)fHistos->At(kHClusters+used*200+ilr*100 + kClTypevsEta);
418         h2d->Fill(eta, clus->GetClusterType());
419         h1d = (TH1F*)fHistos->At(kHClusters+used*200+ilr*100 + kClZ);
420         h1d->Fill(z);
421         h1d = (TH1F*)fHistos->At(kHClusters+used*200+ilr*100 + kClEta);
422         h1d->Fill(eta);
423         Double_t weight = 1.;
424         if (stack) {
425           for(Int_t iLabel = 0; iLabel < 3; ++iLabel) {
426             Int_t clLabel = clus->GetLabel(iLabel);
427             if (clLabel <= 0) continue;
428             Int_t moLabel = FindMotherParticle(stack, clLabel);
429             if (moLabel <= 0) continue;
430             TParticle* p = stack->Particle(moLabel);
431             if ((TMath::Abs(p->GetPdgCode()) != 211) &&
432                 (TMath::Abs(p->GetPdgCode()) != 321) &&
433                 (TMath::Abs(p->GetPdgCode()) != 2212)) continue;
434             weight = PtWeight(p->Pt(),p->GetPdgCode());
435             h1d = (TH1F*)fHistos->At(kHClusters+used*200+ilr*100 + kClZPionsW);
436             h1d->Fill(z,weight);
437             h1d = (TH1F*)fHistos->At(kHClusters+used*200+ilr*100 + kClEtaPionsW);
438             h1d->Fill(eta,weight);
439             h1d = (TH1F*)fHistos->At(kHClusters+used*200+ilr*100 + kClZPions);
440             h1d->Fill(z);
441             h1d = (TH1F*)fHistos->At(kHClusters+used*200+ilr*100 + kClEtaPions);
442             h1d->Fill(eta);
443             break;
444           }
445         }
446       }
447     } // loop on clusters of layer
448   } // loop on layers
449   //
450 }
451
452 //________________________________________________________________________
453 void AliSPDClustTask::SetInput(const char *filename)
454 {
455   // Read pions Pt spectra
456   // from a local input file
457
458   AliInfo(Form("Reading input histograms from %s",filename));
459   if (fUseMC) {
460     TFile *fPt = TFile::Open("FitOutput.root");
461     if (!fPt) { AliError("Failed to open input file"); return; }
462     fhPtPionIn = (TH1F*)fPt->Get("pionPtWeight")->Clone("fhPtPionIn");
463     fhPtPionIn->SetDirectory(0);
464     fhPtKaonIn = (TH1F*)fPt->Get("kaonPtWeight")->Clone("fhPtKaonIn");
465     fhPtKaonIn->SetDirectory(0);
466     fhPtProtonIn = (TH1F*)fPt->Get("protPtWeight")->Clone("fhPtProtonIn");
467     fhPtProtonIn->SetDirectory(0);
468     fPt->Close();
469   }
470 }
471
472 //________________________________________________________________________
473 Int_t AliSPDClustTask::FindMotherParticle(AliStack* stack, Int_t i)
474 {
475   if(stack->IsPhysicalPrimary(i)) return i;
476   TParticle* p = stack->Particle(i);
477   Int_t imo =  p->GetFirstMother();
478
479   //  printf("imo = %d\n",imo);
480
481   if(imo<=0)
482     {
483       AliInfo("particle is not a PhysPrim and has no mother");
484       return imo;
485     }
486   return FindMotherParticle(stack, imo);
487
488 }
489
490 //________________________________________________________________________
491 Double_t AliSPDClustTask::PtWeight(Double_t pt, Int_t pdgCode)
492 {
493   switch (pdgCode) {
494   case 211:
495     return fhPtPionIn->GetBinContent(fhPtPionIn->FindBin(pt));
496   case 321:
497     return fhPtKaonIn->GetBinContent(fhPtKaonIn->FindBin(pt));
498   case 2212:
499     return fhPtProtonIn->GetBinContent(fhPtProtonIn->FindBin(pt));
500   default:
501     return 1.;
502   }
503 }