Analysis task to study the extra SPD clusters.
[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   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 }