New TTask based method to do Digits To clusters. Works with files of multiple
[u/mrichter/AliRoot.git] / MUON / AliMUONproto.cxx
CommitLineData
a9e2aefa 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$Log$
2ab0c725 17Revision 1.3 2000/10/11 09:19:12 egangler
18Corrected some bugs - it should compile now
19
4bc2de9c 20Revision 1.2 2000/06/15 07:58:49 morsch
21Code from MUON-dev joined
22
a9e2aefa 23Revision 1.1.2.1 2000/04/18 09:11:15 morsch
24Implementation of MUON Chamber Prototype Class
25Both read digits from raw data or use the Monte-Carlo.
26Rachid GUERNANE, IPN Lyon guernane@ipnl.in2p3.fr
27
28*/
29
30/*
31Implementation of MUON Chamber Prototype Class
32Both read digits from raw data or use the Monte-Carlo.
331-Feb-2000 Rachid GUERNANE, IPN Lyon guernane@ipnl.in2p3.fr
34*/
35
36
37////////////////////////////////////////////////
38// Manager and hits classes for set:PROTO //
39////////////////////////////////////////////////
40
41#include <TTUBE.h>
42#include <TBRIK.h>
43#include <TRotMatrix.h>
44#include <TNode.h>
45#include <TTree.h>
46#include <TRandom.h>
47#include <TObject.h>
48#include <TVector.h>
49#include <TObjArray.h>
50#include <TMinuit.h>
51#include <TParticle.h>
52#include <TROOT.h>
53#include <TFile.h>
54#include <TNtuple.h>
55#include <TCanvas.h>
56#include <TPad.h>
57#include <TDirectory.h>
58#include <TObjectTable.h>
59#include <AliPDG.h>
60#include <TArrayI.h>
61#include <AliDetector.h>
62
4bc2de9c 63#include "AliMUONChamber.h"
a9e2aefa 64#include "AliMUONproto.h"
4bc2de9c 65#include "AliMUONHit.h"
a9e2aefa 66#include "TTUBE.h"
67#include "AliMUONClusterFinder.h"
68#include "AliRun.h"
69#include "AliMC.h"
70#include "iostream.h"
71#include "AliCallf77.h"
72#include "AliConst.h"
4bc2de9c 73//#include "chainalice2.h"
74#include "AliMUONSegmentationV0.h"
75//#include "AliMUONSegResV11.h"
a9e2aefa 76
77ClassImp(AliMUONproto)
78
79//___________________________________________
80AliMUONproto::AliMUONproto()
81 : AliMUON()
82{
83 cout << "\n Calling AliMUONproto constructor..." << endl;
84
85 //
86 //
87 //
88}
89
90//___________________________________________
91AliMUONproto::AliMUONproto(const char *name, const char *title)
92 : AliMUON(name,title)
93// AliMUON defines everything, but the chamber for NCH=1
94{
95
96// z-Positions of Chambers
97 const Float_t zch[1] = {975.};
98//
99// inner diameter
100 const Float_t dmi[1] = {0.};
101//
102// outer diameter
103 const Float_t dma[1] = {29.2};
104//
105//
106// Default Parameters for ALICE2 prototype
107
108//
109 (*fChambers)[0] = new AliMUONChamber();
110 AliMUONChamber* chamber = (AliMUONChamber*) (*fChambers)[0];
111 chamber->SetGid(0);
4bc2de9c 112 chamber->SetZ(zch[0]);
a9e2aefa 113//
114 chamber->InitGeo(zch[0]);
115 chamber->SetRInner(dmi[0]/2);
116 chamber->SetROuter(dma[0]/2);
117
118 for (int i = 0; i <= 99; i++) {
119 fThreshold[i] = 0.;
120 }
121
122}
123
4bc2de9c 124#ifdef WE_FORGET_FOR_THE_MOMENT
a9e2aefa 125//___________________________________________
126void AliMUONproto::GetRawDigits(Int_t evnb, Int_t *lptr, Int_t ilen) {
127
128 Int_t ip = 0;
129 Int_t equip = 0;
130 Int_t nochnl;
131 Int_t loop;
132 Int_t val;
133 Int_t itype;
134 Int_t id;
135 Int_t serial;
136 Int_t equiplen;
137
138 Int_t digits[5];
139
140 AliMUON* MUON = (AliMUON*)gAlice->GetModule("MUON");
4bc2de9c 141 AliMUONSegmentationV0* seg = (AliMUONSegmentationV0*) Chamber(0).SegmentationModel(1);
a9e2aefa 142
143 Int_t tracks[10];
144 Int_t charges[10];
145
146 for (Int_t i = 0; i < 10; i++) {
147 tracks[i] = 0;
148 charges[i] = 0;
149 }
150
151 Int_t ich = 0;
152
153 nwtype:
154
155 itype = lptr[ip++];
156 id = lptr[ip++];
157 equiplen = lptr[ip++];
158
159 if (equiplen < 0 ) return;
160
161 if (itype != (int)(0XCACCA008)) {
162 ip += equiplen;
163 if (ip < ilen-2) goto nwtype;
164 }
165 else {
166 serial = id >> 16;
167 equip = id & 0x1;
168 if ((serial == 190) && (equip == 1)) {
169 for (loop = 0; loop < equiplen; loop++) {
170 nochnl = (lptr[ip] & 0x7ff000 ) >> 12;
171 val = lptr[ip] & 0x3ff;
172 // fill digits from raw data according to cathode connexions
173 if (group[nochnl][2][1]!=0)
174 digits[0] = group[nochnl][2][1] - 12;
175 else if (group[nochnl][1][1]!=0)
176 digits[0] = group[nochnl][1][1] - 12;
177 else
178 digits[0] = group[nochnl][0][1] - 12;
179 digits[1] = group[nochnl][0][0];
180 if (digits[0] != seg->Ix(digits[0], digits[1]))
181 printf("Pb pour ix=%d,iy=%d\n", digits[0], digits[1]);
182 digits[2] = val;
183 digits[3] = 0;
184 digits[4] = 0;
185 if (digits[2] >= fThreshold[nochnl])
186 MUON->AddDigits(ich, tracks, charges, digits);
187
188 ip++;
189 }
190 }
191 else
192 ip += equiplen;
193
194 if (ip < ilen-2) goto nwtype;
195 }
196
197 gAlice->TreeD()->Fill();
198 MUON->ResetDigits();
199
200 gAlice->TreeD()->Fill();
201 MUON->ResetDigits();
202
203 char hname[30];
204 sprintf(hname, "TreeD%d", evnb);
2ab0c725 205 gAlice->TreeD()->Write(hname,TObject::kOverwrite);
a9e2aefa 206 // reset tree
207 gAlice->TreeD()->Reset();
208
209}
4bc2de9c 210#endif
a9e2aefa 211
212//___________________________________________
213void AliMUONproto::SetPadSize(Int_t id, Int_t isec, Float_t p1, Float_t p2)
214{
215 Int_t i=2*(id-1);
216 ((AliMUONChamber*)(*fChambers)[i])->SetPadSize(isec,p1,p2);
217}
218
219//___________________________________________
220void AliMUONproto::SetChargeSlope(Int_t id, Float_t p1)
221{
222 Int_t i=2*(id-1);
223 ((AliMUONChamber*) (*fChambers)[i])->SetChargeSlope(p1);
224}
225
226//___________________________________________
227void AliMUONproto::SetChargeSpread(Int_t id, Float_t p1, Float_t p2)
228{
229 Int_t i=2*(id-1);
230 ((AliMUONChamber*) (*fChambers)[i])->SetChargeSpread(p1,p2);
231}
232
233//___________________________________________
234void AliMUONproto::SetMaxAdc(Int_t id, Float_t p1)
235{
236 Int_t i=2*(id-1);
237 ((AliMUONChamber*) (*fChambers)[i])->SetMaxAdc(p1);
238}
239
240//___________________________________________
241void AliMUONproto::CreateGeometry()
242{
243 Int_t *idtmed = fIdtmed->GetArray()-1099;
244//
245 Float_t zpos;
246 Float_t bpar[3];
247 Int_t idGas=idtmed[1100];
248
249 AliMUONChamber *iChamber;
250//********************************************************************
251// Prototype ALICE2 **
252//********************************************************************
253 iChamber=(AliMUONChamber*) (*fChambers)[0];
254
4bc2de9c 255 zpos=iChamber->Z();
a9e2aefa 256//
257 const Float_t X_POS = 11*.975/2; //half x
258 const Float_t Y_POS = 18*.55/2; //half y
259 const Float_t Z_POS = 0.325;
260
261 bpar[0] = X_POS;
262 bpar[1] = Y_POS;
263 bpar[2] = Z_POS;
264
265 gMC->Gsvolu("C01G", "BOX", idGas, bpar, 3);
266 gMC->Gspos("C01G", 1, "ALIC", -bpar[0], bpar[1], zpos, 0, "ONLY");
267}
268
269//___________________________________________
270void AliMUONproto::CreateMaterials()
271{
272 // *** DEFINITION OF AVAILABLE MUON MATERIALS ***
273 //
274 // Ar-CO2 gas
275 Float_t ag1[3] = { 39.95,12.01,16. };
276 Float_t zg1[3] = { 18.,6.,8. };
277 Float_t wg1[3] = { .8,.0667,.13333 };
278 Float_t dg1 = .001821;
279
280 Float_t epsil = .001; // Tracking precision,
281 Float_t tmaxfd = -20.; // Maximum angle due to field deflection
282 Float_t stmin = -.8;
283 Int_t ISXFLD = gAlice->Field()->Integ();
284 Float_t SXMGMX = gAlice->Field()->Max();
285
286 //
287 // --- Define the various materials for GEANT ---
288 AliMixture(22, "ArCO2 80%$", ag1, zg1, dg1, 3, wg1);
289
290 //
291 // Ar/CO2
292 AliMedium(1, "ARG_CO2 ", 22, 1, ISXFLD, SXMGMX, tmaxfd, fMaxStepGas,
293 fMaxDestepAlu, epsil, stmin);
294 // Air
295 //AliMedium(1, "AIR_CH_US ", 15, 1, ISXFLD, SXMGMX, tmaxfd, -1., -.3, epsil, stmin);
296}
297
298//___________________________________________
299void AliMUONproto::Init()
300{
301 printf("\n\n\n Start Init for Prototype ALICE2 - CPC chamber type\n\n\n");
302
303 //
304 // Initialize Tracking Chambers
305 //
306
307 ( (AliMUONChamber*) (*fChambers)[0])->Init();
308
309 //
310 // Set the chamber (sensitive region) GEANT identifier
311 AliMC* gMC = AliMC::GetMC();
312 ((AliMUONChamber*)(*fChambers)[0])->SetGid(gMC->VolId("C01G"));
313
314 printf("\n\n\n Finished Init for Prototype ALICE2 - CPC chamber type\n\n\n");
315
316}
317
318//___________________________________________
319void AliMUONproto::StepManager()
320{
321 Int_t copy, id;
322 static Int_t idvol;
323 static Int_t vol[2];
324 Int_t ipart;
325 TLorentzVector pos;
326 TLorentzVector mom;
327 Float_t theta,phi;
328 Float_t destep, step;
329
330 static Float_t eloss, eloss2, xhit, yhit, tlength;
331 const Float_t big=1.e10;
332
333 // modifs perso
334 static Float_t hits[14];
335
336 TClonesArray &lhits = *fHits;
337
338 //
339 // Set maximum step size for gas
340 // numed=gMC->GetMedium();
341 //
342 // Only charged tracks
343 if( !(gMC->TrackCharge()) ) return;
344 //
345 // Only gas gap inside chamber
346 // Tag chambers and record hits when track enters
347 idvol=-1;
348 id=gMC->CurrentVolID(copy);
349
4bc2de9c 350 for (Int_t i=1; i<=kNCH; i++) {
a9e2aefa 351 if(id==((AliMUONChamber*)(*fChambers)[i-1])->GetGid()){
352 vol[0]=i;
353 idvol=i-1;
354 }
355 }
356 if (idvol == -1) return;
357 //
358 // Get current particle id (ipart), track position (pos) and momentum (mom)
359 gMC->TrackPosition(pos);
360 gMC->TrackMomentum(mom);
361
362 ipart = gMC->TrackPid();
363 //Int_t ipart1 = gMC->IdFromPDG(ipart);
364 //printf("ich, ipart %d %d \n",vol[0],ipart1);
365
366 //
367 // momentum loss and steplength in last step
368 destep = gMC->Edep();
369 step = gMC->TrackStep();
370
371 //
372 // record hits when track enters ...
373 if( gMC->IsTrackEntering()) {
374 gMC->SetMaxStep(fMaxStepGas);
375 Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
376 Double_t rt = TMath::Sqrt(tc);
377 Double_t pmom = TMath::Sqrt(tc+mom[2]*mom[2]);
378 Double_t tx=mom[0]/pmom;
379 Double_t ty=mom[1]/pmom;
380 Double_t tz=mom[2]/pmom;
381 Double_t s=((AliMUONChamber*)(*fChambers)[idvol])
4bc2de9c 382 ->ResponseModel()
a9e2aefa 383 ->Pitch()/tz;
384 theta = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
385 phi = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
386 hits[0] = Float_t(ipart); // Geant3 particle type
387 hits[1] = pos[0]+s*tx; // X-position for hit
388 hits[2] = pos[1]+s*ty; // Y-position for hit
389 hits[3] = pos[2]+s*tz; // Z-position for hit
390 hits[4] = theta; // theta angle of incidence
391 hits[5] = phi; // phi angle of incidence
4bc2de9c 392 hits[8] = (Float_t) fNPadHits; // first padhit
a9e2aefa 393 hits[9] = -1; // last pad hit
394
395 // modifs perso
396 hits[10] = mom[3]; // hit momentum P
397 hits[11] = mom[0]; // Px/P
398 hits[12] = mom[1]; // Py/P
399 hits[13] = mom[2]; // Pz/P
400 // fin modifs perso
401
402 // phi angle of incidence
403 tlength = 0;
404 eloss = 0;
405 eloss2 = 0;
406 xhit = pos[0];
407 yhit = pos[1];
408 // Only if not trigger chamber
409 if(idvol<10) {
410 //
411 // Initialize hit position (cursor) in the segmentation model
412 ((AliMUONChamber*) (*fChambers)[idvol])
413 ->SigGenInit(pos[0], pos[1], pos[2]);
414 } else {
415 //geant3->Gpcxyz();
416 //printf("In the Trigger Chamber #%d\n",idvol-9);
417 }
418 }
419 eloss2+=destep;
420
421 //
422 // Calculate the charge induced on a pad (disintegration) in case
423 //
424 // Mip left chamber ...
425 if( gMC->IsTrackExiting() || gMC->IsTrackStop() || gMC->IsTrackDisappeared()){
426 gMC->SetMaxStep(big);
427 eloss += destep;
428 tlength += step;
429
430 // Only if not trigger chamber
431 if(idvol<10) {
432 if (eloss > 0)
4bc2de9c 433 MakePadHits(0.5*(xhit+pos[0]),0.5*(yhit+pos[1]),eloss,0.0,idvol);
a9e2aefa 434 }
435
436 hits[6]=tlength;
437 hits[7]=eloss2;
4bc2de9c 438 if (fNPadHits > (Int_t)hits[8]) {
a9e2aefa 439 hits[8]= hits[8]+1;
4bc2de9c 440 hits[9]= (Float_t) fNPadHits;
a9e2aefa 441 }
442
443 new(lhits[fNhits++])
444 AliMUONHit(fIshunt,gAlice->CurrentTrack(),vol,hits);
445 eloss = 0;
446 //
447 // Check additional signal generation conditions
448 // defined by the segmentation
449 // model (boundary crossing conditions)
450 } else if
451 (((AliMUONChamber*) (*fChambers)[idvol])
452 ->SigGenCond(pos[0], pos[1], pos[2]))
453 {
454 ((AliMUONChamber*) (*fChambers)[idvol])
455 ->SigGenInit(pos[0], pos[1], pos[2]);
456// printf("\n-> MakePadHits, reason special %d",ipart);
457 if (eloss > 0)
4bc2de9c 458 MakePadHits(0.5*(xhit+pos[0]),0.5*(yhit+pos[1]),eloss,0.0,idvol);
a9e2aefa 459 xhit = pos[0];
460 yhit = pos[1];
461 eloss = destep;
462 tlength += step ;
463 //
464 // nothing special happened, add up energy loss
465 } else {
466 eloss += destep;
467 tlength += step ;
468 }
469}
470
471//___________________________________________
472void AliMUONproto::BuildGeometry()
473{
474 TNode *Node;
475 TNode *Top;
476
477 const int kColorMUON = kBlue;
478 //
479 Top=gAlice->GetGeometry()->GetNode("alice");
480 //
481 //
482 //
483 float dx, dy, dz, zpos;
484
485 const Float_t cz[1]={975.};
486
487 zpos=cz[0];
488
489 dx=11*.975/2;
490 dy=18*.55/2;
491 dz=0.325;
492
493 new TBRIK("C_MUON101","Mother volume for Proto.","void",dx*2,dy*2,dz*2);
494 TBRIK* PROTO = new TBRIK("PROT", "Proto. Sens. Region","void",dx,dy,dz);
495 Top->cd();
496 Node = new TNode("MUON101","ChamberNode","C_MUON101",0,0,zpos,"");
497 Node->SetLineColor(kColorMUON);
498 Node->SetVisibility(0);
499 fNodes->Add(Node);
500 Node->cd();
501 Node = new TNode("MUON201", "Proto. Sens. Region Node", PROTO, -dx, dy, dz);
502 Node->SetLineColor(kColorMUON);
503
504}
505
506//_____________________________________________________________________________
507void AliMUONproto::FindClusters(Int_t nev, Int_t last_entry)
508{
509
510//
511// Loop on chambers and on cathode planes
512//
513 for (Int_t icat = 0; icat < 2; icat++) {
514 gAlice->ResetDigits();
515 gAlice->TreeD()->GetEvent(last_entry+icat);
516
4bc2de9c 517 for (Int_t ich = 0; ich < kNTrackingCh; ich++) {
a9e2aefa 518 AliMUONChamber* iChamber=(AliMUONChamber*) (*fChambers)[ich];
519 TClonesArray *MUONdigits = this->DigitsAddress(ich);
520 if (MUONdigits == 0) continue;
521 //
522 // Get ready the current chamber stuff
523 //
524
4bc2de9c 525 AliMUONResponse* response = iChamber->ResponseModel();
526 AliMUONSegmentation* seg = iChamber->SegmentationModel(icat+1);
527 AliMUONClusterFinder* rec = iChamber->ReconstructionModel();
a9e2aefa 528
529 if (seg) {
530 rec->SetSegmentation(seg);
531 rec->SetResponse(response);
532 rec->SetDigits(MUONdigits);
533 rec->SetChamber(ich);
534 rec->FindRawClusters();
535 }
536
537 TClonesArray *fRch;
538 fRch=RawClustAddress(ich);
539 fRch->Sort();
540 // it seems to work
541
542 } // for ich
543 // fill the tree
544 gAlice->TreeR()->Fill();
545
546 ResetRawClusters();
547 } // for icat
548
549 char hname[30];
550 sprintf(hname,"TreeR%d",nev);
551 gAlice->TreeR()->Write(hname);
552 gAlice->TreeR()->Reset();
553}
554
4bc2de9c 555#ifdef WE_FORGRET_THIS_SHIT
a9e2aefa 556//_____________________________________________________________________________
557void AliMUONproto::SetThreshold()
558{
559
560 ifstream inputFile("/home/alice/guernane/aliroot/pro/MUON/crped190.dat", ios::in);
561
562 if (inputFile.fail()) {
563 cout << "Error opening file" << endl;
564 exit(2);
565 }
566
567 char buff[32];
568 Int_t Serial;
569 Int_t Ntrigger;
570 Int_t Nchannel;
571 Int_t i1;
572 Int_t i2;
573
574 inputFile >> buff;
575
576 inputFile >> Serial;
577 inputFile >> Ntrigger;
578 inputFile >> Nchannel;
579 inputFile >> i1;
580 inputFile >> i2;
581
582 Float_t ped0[Nchannel];
583 Float_t sig0[Nchannel];
584 Float_t ped1[Nchannel];
585 Float_t sig1[Nchannel];
586 Int_t ichannel;
587
588 for (Int_t i = 0; i < Nchannel-1; i++) {
589 ped0[i] = 0;
590 sig0[i] = 0;
591 ped1[i] = 0;
592 sig1[i] = 0;
593 }
594
595 for (Int_t i = 0; i < Nchannel-1; i++) {
596 inputFile >> ichannel;
597 inputFile >> ped0[i];
598 inputFile >> sig0[i];
599 inputFile >> ped1[i];
600 inputFile >> sig1[i];
601 fThreshold[i] = fNsigma*sig1[i];
602 }
603
604 inputFile.close();
605}
4bc2de9c 606#endif