[u/mrichter/AliRoot.git] / ITS / SSDrecpointTest.C

void SSDrecpointTest (Int_t evNumber1=0,Int_t evNumber2=0)
  //void SSDrecpointTest (Int_t evNumber1=0,Int_t evNumber2=999)
{
/////////////////////////////////////////////////////////////////////////
//   This macro is a small example of a ROOT macro
//   illustrating how to read the output of GALICE
//   and fill some histograms.
//   
//     Root > .L anal.C   //this loads the macro in memory
//     Root > anal();     //by default process first event   
//     Root > anal(2);    //process third event
//Begin_Html
/*
<img src="gif/anal.gif">
*/
//End_Html
/////////////////////////////////////////////////////////////////////////
    
// Dynamically link some shared libs

   if (gClassTable->GetID("AliRun") < 0) {
      gROOT->LoadMacro("loadlibs.C");
      loadlibs();
   } else {
      delete gAlice;
      gAlice=0;
   }

// Connect the Root Galice file containing Geometry, Kine and Hits
   TString *str = new TString("galice.root");
   TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(str->Data());
   if (!file) file = new TFile(str->Data(),"UPDATE");

// Get AliRun object from file or create it if not on file
   //   if (!gAlice) {
     gAlice = (AliRun*)file->Get("gAlice");
     if (gAlice) printf("AliRun object found on file\n");
     if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
     //}


     // -------------- Create ntuples --------------------

     //  ntuple structures:


          struct {
            Int_t lay;
            Int_t nxP;
            Int_t nxN;
            Int_t hitprim;
            Int_t partcode;
            Int_t ntrover;
            Float_t x;
            Float_t z;
            Float_t dx;
            Float_t dz;
            Float_t pmod;
          } ntuple_st;

          struct {
            Int_t lay;
            Int_t lad;
            Int_t det;
            Int_t nxP;
            Int_t nxN;
            Int_t noverlaps;
            Int_t noverprim;
            Int_t ntrover;
            Float_t qclP;
            Float_t qclN;
            Float_t qrec;
            Float_t dx;
            Float_t dz;
          } ntuple1_st;

          struct {
            Int_t nxP;
            Int_t nxN;
            Float_t x;
            Float_t z;
          } ntuple2_st;

          ntuple = new TTree("ntuple","Demo ntuple");
          ntuple->Branch("lay",&ntuple_st.lay,"lay/I");
          ntuple->Branch("nxP",&ntuple_st.nxP,"nxP/I");
          ntuple->Branch("nxN",&ntuple_st.nxN,"nxN/I");
          ntuple->Branch("hitprim",&ntuple_st.hitprim,"hitprim/I");
          ntuple->Branch("partcode",&ntuple_st.partcode,"partcode/I");
          ntuple->Branch("ntrover",&ntuple_st.ntrover,"ntrover/I");
          ntuple->Branch("x",&ntuple_st.x,"x/F");
          ntuple->Branch("z",&ntuple_st.z,"z/F");
          ntuple->Branch("dx",&ntuple_st.dx,"dx/F");
          ntuple->Branch("dz",&ntuple_st.dz,"dz/F");
          ntuple->Branch("pmod",&ntuple_st.pmod,"pmod/F");

	  ntuple1 = new TTree("ntuple1","Demo ntuple1");
	  ntuple1->Branch("lay",&ntuple1_st.lay,"lay/I");
	  ntuple1->Branch("lad",&ntuple1_st.lad,"lad/I");
	  ntuple1->Branch("det",&ntuple1_st.det,"det/I");
	  ntuple1->Branch("nxP",&ntuple1_st.nxP,"nxP/I");
	  ntuple1->Branch("nxN",&ntuple1_st.nxN,"nxN/I");
	  ntuple1->Branch("qclP",&ntuple1_st.qclP,"qclP/F");
          ntuple1->Branch("qclN",&ntuple1_st.qclN,"qclN/F");
          ntuple1->Branch("qrec",&ntuple1_st.qrec,"qrec/F");
          ntuple1->Branch("dx",&ntuple1_st.dx,"dx/F");
          ntuple1->Branch("dz",&ntuple1_st.dz,"dz/F");
	  ntuple1->Branch("noverlaps",&ntuple1_st.noverlaps,"noverlaps/I");
          ntuple1->Branch("noverprim",&ntuple1_st.noverprim,"noverprim/I");
          ntuple1->Branch("ntrover",&ntuple1_st.ntrover,"ntrover/I");

          ntuple2 = new TTree("ntuple2","Demo ntuple2");
          ntuple2->Branch("nxP",&ntuple2_st.nxP,"nxP/I");
          ntuple2->Branch("nxN",&ntuple2_st.nxN,"nxN/I");
          ntuple2->Branch("x",&ntuple2_st.x,"x/F");
          ntuple2->Branch("z",&ntuple2_st.z,"z/F");


	  // Create Histogramms

	  TH1F *NxP5 = new TH1F("NxP5","P cluster size for layer 5",20,0.,20.);
	  TH1F *NxN5 = new TH1F("NxN5","N cluster size for layer 5",20,0.,20.);
	  TH1F *NxP6 = new TH1F("NxP6","P cluster size for layer 6",20,0.,20.);
	  TH1F *NxN6 = new TH1F("NxN6","N cluster size for layer 6",20,0.,20.);

	  TH1F *Xres5 = new TH1F("Xres5","Xrec and Xgen difference (micr) for layers 5",100,-200.,200.);
	  TH1F *Xres6 = new TH1F("Xres6","Xrec and Xgen difference (micr) for layers 6",100,-200.,200.);
	  TH1F *Zres5 = new TH1F("Zres5","Zrec and Zgen difference (micr) for layers 5",100,-8000.,8000.);
	  TH1F *Zres6 = new TH1F("Zres6","Zrec and Zgen difference (micr) for layers 6",100,-8000.,8000.);
          TH1F *Path5 = new TH1F("Path5","Path length in Si",100,0.,600.);
          TH1F *Path6 = new TH1F("Path6","Path length in Si",100,0.,600.);
          TH1F *dEdX = new TH1F("dEdX","dEdX  (KeV)",100,0.,500.);
          TH2F *adcPadcN5all = new TH2F("adcPadcN5all","adcP/N correlation for lay5",100,0.,200.,100,0.,200.);
          TH2F *adcPadcN6all = new TH2F("adcPadcN6all","adcP/N correlation for lay6",100,0.,200.,100,0.,200.);
          TH2F *adcPadcN5cut = new TH2F("adcPadcN5cut","adcP/N correlation for lay5 and cut of P-N signas",100,0.,200.,100,0.,200.);
          TH2F *adcPadcN6cut = new TH2F("adcPadcN6cut","adcP/N correlation for lay6 and cut of P-N signals",100,0.,200.,100,0.,200.);


   AliITS *ITS  = (AliITS*) gAlice->GetModule("ITS");
   if (!ITS) { cout << "no ITS" << endl; return; }
   
   //AliITSgeom *aliitsgeo = ITS->GetITSgeom();
   AliITSgeom *geom = ITS->GetITSgeom();


   //Int_t cp[8]={0,0,0,0,0,0,0,0};

   //cout << "SSD" << endl;

   AliITSDetType *iDetType=ITS->DetType(2);
   AliITSsegmentationSSD *seg2=(AliITSsegmentationSSD*)iDetType->GetSegmentationModel();
   AliITSresponseSSD *res2 = (AliITSresponseSSD*)iDetType->GetResponseModel();
   //res2->SetSigmaSpread(3.,2.);
   AliITSsimulationSSD *sim2=new AliITSsimulationSSD(seg2,res2);
   ITS->SetSimulationModel(2,sim2);

   TClonesArray *dig2  = ITS->DigitsAddress(2);
   TClonesArray *recp2  = ITS->ClustersAddress(2);
   AliITSClusterFinderSSD *rec2=new AliITSClusterFinderSSD(seg2,dig2);
   ITS->SetReconstructionModel(2,rec2);

   // test
   printf("SSD dimensions %f %f %f \n",seg2->Dx(),seg2->Dz(),seg2->Dy());
   printf("SSD nstrips %d %d \n",seg2->Npz(),seg2->Npx());
   Float_t ylim = seg2->Dy()/2 - 12;
   
//
//   Loop over events
//


   Int_t Nh=0;
   Int_t Nh1=0;
   for (int nev=0; nev<= evNumber2; nev++) {
     Int_t nparticles = 0;
     nparticles = gAlice->GetEvent(nev);
     cout << "nev         " << nev <<endl;
     cout << "nparticles  " << nparticles <<endl;
     if (nev < evNumber1) continue;
     if (nparticles <= 0) return;
     
     AliITShit *itsHit;
     AliITSRecPoint *itsPnt = 0;
     AliITSRawClusterSSD *itsClu = 0;
     
     // Get Hit, Cluster & Recpoints Tree Pointers

     TTree *TH = gAlice->TreeH();
     Int_t nenthit=TH->GetEntries();
     printf("Found %d entries in the Hit tree (must be one per track per event!)\n",nenthit);

     ITS->GetTreeC(nev);
     TTree *TC=ITS->TreeC();
     Int_t nentclu=TC->GetEntries();
     printf("Found %d entries in the Cluster tree (must be one per module per event!)\n",nentclu);

     TTree *TR = gAlice->TreeR();
     Int_t nentrec=TR->GetEntries();
     printf("Found %d entries in the RecPoints tree\n",nentrec);

     // Get Pointers to Clusters & Recpoints TClonesArrays

     TClonesArray *ITSclu  = ITS->ClustersAddress(2); 
     printf ("ITSclu %p \n",ITSclu);
     TClonesArray *ITSrec  = ITS->RecPoints(); 
     printf ("ITSrec %p \n",ITSrec);

     // check recpoints

     //Int_t nbytes = 0;
     Int_t totpoints = 0;
     Int_t totclust = 0;

     // check hits
     
     Int_t nmodules=0;
     Int_t mod;
     
     ITS->InitModules(-1,nmodules); 
     ITS->FillModules(nev,0,nmodules,"","");
     
     TObjArray *fITSmodules = ITS->GetModules();
     
     Int_t first0 = geom->GetStartDet(0);  // SPD
     Int_t last0 = geom->GetLastDet(0);    // SPD
     Int_t first1 = geom->GetStartDet(1);  // SDD
     Int_t last1 = geom->GetLastDet(1);    // SDD
     Int_t first2 = geom->GetStartDet(2);  // SSD
     Int_t last2 = geom->GetLastDet(2);    // SSD

     //  For the SPD: first0 = 0, last0 = 239     (240 modules);  
     //  for the SDD: first1 = 240, last1 = 499   (260 modules);  
     //  for the SSD: first2 = 500, last2 = 2269  (1770 modules).  

     printf("det type %d first0, last0 %d %d \n",0,first0,last0);
     printf("det type %d first1, last1 %d %d \n",1,first1,last1);
     printf("det type %d first2, last2 %d %d \n",2,first2,last2);

     // module loop for the SSD
     for (mod=first2; mod<last2+1; mod++) {  // for the "ALL" option
       //for (mod=0; mod<last2-first2+1; mod++) { //for the "SSD" option

       TTree *TR = gAlice->TreeR();
       Int_t nentrec=TR->GetEntries();
       //printf("Found %d entries in the RecPoints tree\n",nentrec);
      
              //cout << "CLUSTERS: reset" << endl;
       ITS->ResetClusters();
       //cout << "CLUSTERS: get" << endl;
       TC->GetEvent(mod);
       //cout << "RECPOINTS: reset" << endl;
       ITS->ResetRecPoints();
       //cout << "RECPOINTS: get" << endl;
       //TR->GetEvent(mod+1);   // for the V3.04 AliRoot
       TR->GetEvent(mod);       // for the V3.05 AliRoot

       Int_t nrecp = ITSrec->GetEntries();
       totpoints += nrecp;
       //if (nrecp) printf("Found %d rec points for module %d\n",nrecp,mod);
       if (!nrecp) continue;
       Int_t nclusters = ITSclu->GetEntries();
       totclust += nclusters;
       //if (nclusters) printf("Found %d clusters for module %d\n",nrecc,mod);
       
       //AliITSmodule *Mod = (AliITSmodule *)fITSmodules->At(mod+first2);
       // for the "SSD" option

       AliITSmodule *Mod = (AliITSmodule *)fITSmodules->At(mod);
       // for the "ALL" option

       //       printf("Mod: %X\n",Mod);
       Int_t nhits = Mod->GetNhits();
       Float_t epart = 0;
       cout <<" module,nrecp,nclusters,nhits ="<<mod<<","<<nrecp<<","<<nclusters<<","<<nhits<< endl;

       // ---------------- cluster/hit analysis ---------------------


     Float_t pathInSSD = 300.;

       // ---- Recpoint loop
       for (Int_t pnt=0;pnt<nrecp;pnt++) {
	 itsPnt  = (AliITSRecPoint*)ITSrec->At(pnt);
	 if(!itsPnt) continue;
	 itsClu  = (AliITSRawClusterSSD*)ITSclu->At(pnt);
	 if(!itsClu) continue;

	 Int_t nxP = itsClu->fMultiplicity;
	 Int_t nxN = itsClu->fMultiplicityN;
	 Int_t ntrover = itsClu->fNtracks;
	 Float_t qclP = itsClu->fSignalP;     // in ADC
	 Float_t qclN = itsClu->fSignalN;     // in ADC
	 //Float_t dq = qclP - qclN;
	 Float_t qcut = itsClu->fQErr;        // abs(dq)/signal,
	                                      // where signal is
	                                      // max of qclP,qclN        
 	 Float_t xrec = 10000*itsPnt->GetX();
	 Float_t zrec = 10000*itsPnt->GetZ();
	 Float_t qrec = itsPnt->GetQ();      // in ADC, maximum from fSignalP/N
	 //Float_t dedx = itsPnt->GetdEdX();   // in KeV (ADC * 2.16)
	 Float_t dedx = itsPnt->fdEdX;   // in KeV (ADC * 2.16)
         Int_t ii = 0;
	 Int_t tr1 = itsPnt->GetLabel(ii);
         Int_t ii = 1;
	 Int_t tr2 = itsPnt->GetLabel(ii);
         Int_t ii = 2;
	 Int_t tr3 = itsPnt->GetLabel(ii);
	 // fill ntuple2
	     ntuple2_st.nxP = nxP;
             ntuple2_st.nxN = nxN;
	     ntuple2_st.x = xrec/1000;
             ntuple2_st.z = zrec/1000;

             if(qcut < 0.18) ntuple2->Fill();


	  Int_t noverlaps = 0;
	  Int_t noverprim = 0;
 	  Int_t flaghit = 0;
          Float_t xhit0 = 1e+7;
          Float_t yhit0 = 1e+7;
          Float_t zhit0 = 1e+7;

       // Hit loop
        for (Int_t hit=0;hit<nhits;hit++) {

	 itsHit   = (AliITShit*)Mod->GetHit(hit);

	 Int_t flagtrack = 0;
	 Int_t hitlayer = itsHit->GetLayer();
	 Int_t hitladder= itsHit->GetLadder();
	 Int_t hitdet= itsHit->GetDetector();

	 Int_t track = itsHit->GetTrack();
	 Int_t dray = 0;
	 Int_t hitstat = itsHit->GetTrackStatus();

 	  Float_t zhit = 10000*itsHit->GetZL();
	  Float_t xhit = 10000*itsHit->GetXL();
	  Float_t yhit = 10000*itsHit->GetYL();
	  Float_t ehit = 1.0e+6*itsHit->GetIonization(); // hit energy, KeV 

	   Int_t parent = itsHit->GetParticle()->GetFirstMother();
	   Int_t partcode = itsHit->GetParticle()->GetPdgCode();

   //  partcode (pdgCode): 11 - e-, 13 - mu-, 22 - gamma, 111 - pi0, 211 - i+
   //  310 - K0s, 321 - K+, 2112 - n, 2212 - p, 3122 - lambda

           Float_t pmod = itsHit->GetParticle()->P(); // the momentum at the
	                                              // vertex
	   pmod *= 1.0e+3;

	  if(hitstat == 66 && yhit < -ylim) {  // entering hit
	    xhit0 = xhit;
	    yhit0 = yhit;
	    zhit0 = zhit;
	  }

	  if(hitstat == 66) continue; // Take the not entering hits only 

	  if(xhit0 > 9e+6 || zhit0 > 9e+6 || yhit0 > 9e+6) {
	    //cout<<"default xhit0,zhit0,yhit0 ="<<xhit0<<","<<zhit0<<","<<yhit0<<endl;
	    continue;
	  }


	  // Consider the hits only with the track number equaled to one
	  // of the recpoint
	  if((track == tr1) || (track == tr2) || (track == tr3)) flagtrack = 1;

         if(flagtrack == 1) {     // the hit corresponds to the recpoint

	   flaghit = 1;

	   //Float_t px = itsHit->GetPXL(); // the momenta at this GEANT point
	   //Float_t py = itsHit->GetPYL();
	   //Float_t pz = itsHit->GetPZL();

         Int_t hitprim = 0;

	 if(partcode == 11 && pmod < 6) dray = 1; // delta ray is e-
	                                          // at p < 6 MeV/c

         if((hitstat == 68 || hitstat == 33) && dray == 0)  noverlaps=noverlaps + 1;
                                                  // overlapps for all hits but
	                                          // not for delta ray which
	                                          // also went out from the
	                                          // detector and returned
	                                          // again


	  // x,z resolution colculation
          if((hitstat == 68 || hitsat == 33) && dray == 0) {
  	     Float_t xmed = (xhit + xhit0)/2;
	     Float_t zmed = (zhit + zhit0)/2;
	     Float_t xdif = xmed - xrec;
	     Float_t zdif = zmed - zrec;

            if(parent < 0)  {
	      hitprim = 1; // hitprim=1 for the primery particles
	      noverprim += 1;
	    }
	     pathInSSD = TMath::Sqrt((xhit0-xhit)*(xhit0-xhit)+(yhit0-yhit)*(yhit0-yhit)+(zhit0-zhit)*(zhit0-zhit));

	     // fill ntuple
             ntuple_st.lay = hitlayer;
	     ntuple_st.nxP = nxP;
             ntuple_st.nxN = nxN;
	     ntuple_st.hitprim = hitprim;
             ntuple_st.partcode = partcode;
             ntuple_st.ntrover = ntrover;
	     ntuple_st.x = xrec/1000;
             ntuple_st.z = zrec/1000;
	     ntuple_st.dx = xdif;
             ntuple_st.dz = zdif;
             ntuple_st.pmod = pmod;

             //if(qcut < 0.18) ntuple->Fill();
             ntuple->Fill();

	     //if(hitlayer == 5 && qcut < 0.18) {
             
	    if(hitlayer == 5 ) {
             Xres5->Fill(xdif);
             Zres5->Fill(zdif);
             Path5->Fill(pathInSSD);
	    }
            //if(hitlayer == 6 && qcut < 0.18) {
            if(hitlayer == 6) {
             Xres6->Fill(xdif);
             Zres6->Fill(zdif);
             Path6->Fill(pathInSSD);
            }
	  } // hitstat 68/33
	 } else {       // non correspondent hit
	  xhit0 = 1e+7;
	  zhit0 = 1e+7;
	 } // end of hit-recpoint correspondence
	} // hit loop       

	if(flaghit == 1) {

	  if(noverlaps == 0) noverlaps = 1; // cluster contains one or more
	  // delta rays only

	  // fill ntuple1
	  ntuple1_st.lay = hitlayer;
	  ntuple1_st.lad = hitladder;
	  ntuple1_st.det = hitdet;
	  ntuple1_st.nxP = nxP;
	  ntuple1_st.nxN = nxN;
	  ntuple1_st.qclP = qclP*300/pathInSSD; 
	  ntuple1_st.qclN = qclN*300/pathInSSD; 
	  ntuple1_st.qrec = qrec*300/pathInSSD; 
	  ntuple1_st.dx = xdif;
	  ntuple1_st.dz = zdif;
	  noverlaps -= 1;
	  noverprim -= 1;
	  ntuple1_st.noverlaps = noverlaps;
	  ntuple1_st.noverprim = noverprim;
          ntuple1_st.ntrover = ntrover;

	  //if(qcut < 0.18) ntuple1->Fill();
	  ntuple1->Fill();

          Float_t de = dedx*300./pathInSSD;
          dEdX->Fill(de);
	    if(hitlayer == 5 ) {
             adcPadcN5all->Fill(qclP,qclN);
            }
	    if(hitlayer == 6 ) {
             adcPadcN6all->Fill(qclP,qclN);
            }
	    if(hitlayer == 5 && qcut < 0.18) {
             adcPadcN5cut->Fill(qclP,qclN);
             NxP5->Fill(nxP);
             NxN5->Fill(nxN);
            }
	    if(hitlayer == 6 && qcut < 0.18) {
             adcPadcN6cut->Fill(qclP,qclN);
             NxP6->Fill(nxP);
             NxN6->Fill(nxN);
            }
	} // flaghit = 1
       } //b.b. recpoint loop
     } //b.b. module loop
   } //b.b. evnt loop

   TFile fhistos("SSD_his.root","RECREATE");

   ntuple->Write();
   ntuple1->Write();
   ntuple2->Write();
   NxP5->Write();
   NxN5->Write();
   NxP6->Write();
   NxN6->Write();
   Xres5->Write();
   Zres5->Write();
   Xres6->Write();
   Zres6->Write();
   Path5->Write();
   Path6->Write();
   adcPadcN5all->Write();
   adcPadcN6all->Write();
   adcPadcN5cut->Write();
   adcPadcN6cut->Write();
   dEdX->Write();

   fhistos.Close();

   cout<<"!!! Histogramms and ntuples were written"<<endl;

   TCanvas *c1 = new TCanvas("c1","ITS clusters",400,10,600,700);
   c1->Divide(2,2);
   c1->cd(1);
   gPad->SetFillColor(33);
         Xres5->SetFillColor(42);
         Xres5->Draw();
   c1->cd(2);
   gPad->SetFillColor(33);
         Zres5->SetFillColor(46);
         Zres5->Draw();
   c1->cd(3);
   gPad->SetFillColor(33);
         Xres6->SetFillColor(42);
         Xres6->Draw();
   c1->cd(4);
   gPad->SetFillColor(33);
         Zres6->SetFillColor(46);
         Zres6->Draw();

   cout<<"END  test for clusters and hits "<<endl;

}
Commit	Line	Data
6b8f55ce	1	void SSDrecpointTest (Int_t evNumber1=0,Int_t evNumber2=0)
	2	//void SSDrecpointTest (Int_t evNumber1=0,Int_t evNumber2=999)
	3	{
	4	/////////////////////////////////////////////////////////////////////////
	5	// This macro is a small example of a ROOT macro
	6	// illustrating how to read the output of GALICE
	7	// and fill some histograms.
	8	//
	9	// Root > .L anal.C //this loads the macro in memory
	10	// Root > anal(); //by default process first event
	11	// Root > anal(2); //process third event
	12	//Begin_Html
	13	/*
	14	<img src="gif/anal.gif">
	15	*/
	16	//End_Html
	17	/////////////////////////////////////////////////////////////////////////
	18
	19	// Dynamically link some shared libs
	20
	21	if (gClassTable->GetID("AliRun") < 0) {
	22	gROOT->LoadMacro("loadlibs.C");
	23	loadlibs();
487a6bce	24	} else {
	25	delete gAlice;
	26	gAlice=0;
6b8f55ce	27	}
	28
	29	// Connect the Root Galice file containing Geometry, Kine and Hits
	30	TString *str = new TString("galice.root");
	31	TFile file = (TFile)gROOT->GetListOfFiles()->FindObject(str->Data());
	32	if (!file) file = new TFile(str->Data(),"UPDATE");
	33
	34	// Get AliRun object from file or create it if not on file
	35	// if (!gAlice) {
	36	gAlice = (AliRun*)file->Get("gAlice");
	37	if (gAlice) printf("AliRun object found on file\n");
	38	if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
	39	//}
	40
	41
	42	// -------------- Create ntuples --------------------
	43
	44	// ntuple structures:
	45
	46
	47	struct {
	48	Int_t lay;
	49	Int_t nxP;
	50	Int_t nxN;
	51	Int_t hitprim;
	52	Int_t partcode;
f239b2fe	53	Int_t ntrover;
6b8f55ce	54	Float_t x;
	55	Float_t z;
	56	Float_t dx;
	57	Float_t dz;
	58	Float_t pmod;
	59	} ntuple_st;
	60
	61	struct {
	62	Int_t lay;
	63	Int_t lad;
	64	Int_t det;
	65	Int_t nxP;
	66	Int_t nxN;
	67	Int_t noverlaps;
	68	Int_t noverprim;
f239b2fe	69	Int_t ntrover;
6b8f55ce	70	Float_t qclP;
	71	Float_t qclN;
	72	Float_t qrec;
	73	Float_t dx;
	74	Float_t dz;
	75	} ntuple1_st;
	76
	77	struct {
	78	Int_t nxP;
	79	Int_t nxN;
	80	Float_t x;
	81	Float_t z;
	82	} ntuple2_st;
	83
	84	ntuple = new TTree("ntuple","Demo ntuple");
	85	ntuple->Branch("lay",&ntuple_st.lay,"lay/I");
	86	ntuple->Branch("nxP",&ntuple_st.nxP,"nxP/I");
	87	ntuple->Branch("nxN",&ntuple_st.nxN,"nxN/I");
	88	ntuple->Branch("hitprim",&ntuple_st.hitprim,"hitprim/I");
	89	ntuple->Branch("partcode",&ntuple_st.partcode,"partcode/I");
f239b2fe	90	ntuple->Branch("ntrover",&ntuple_st.ntrover,"ntrover/I");
6b8f55ce	91	ntuple->Branch("x",&ntuple_st.x,"x/F");
	92	ntuple->Branch("z",&ntuple_st.z,"z/F");
	93	ntuple->Branch("dx",&ntuple_st.dx,"dx/F");
	94	ntuple->Branch("dz",&ntuple_st.dz,"dz/F");
	95	ntuple->Branch("pmod",&ntuple_st.pmod,"pmod/F");
	96
	97	ntuple1 = new TTree("ntuple1","Demo ntuple1");
	98	ntuple1->Branch("lay",&ntuple1_st.lay,"lay/I");
	99	ntuple1->Branch("lad",&ntuple1_st.lad,"lad/I");
	100	ntuple1->Branch("det",&ntuple1_st.det,"det/I");
	101	ntuple1->Branch("nxP",&ntuple1_st.nxP,"nxP/I");
	102	ntuple1->Branch("nxN",&ntuple1_st.nxN,"nxN/I");
	103	ntuple1->Branch("qclP",&ntuple1_st.qclP,"qclP/F");
	104	ntuple1->Branch("qclN",&ntuple1_st.qclN,"qclN/F");
	105	ntuple1->Branch("qrec",&ntuple1_st.qrec,"qrec/F");
	106	ntuple1->Branch("dx",&ntuple1_st.dx,"dx/F");
	107	ntuple1->Branch("dz",&ntuple1_st.dz,"dz/F");
	108	ntuple1->Branch("noverlaps",&ntuple1_st.noverlaps,"noverlaps/I");
	109	ntuple1->Branch("noverprim",&ntuple1_st.noverprim,"noverprim/I");
f239b2fe	110	ntuple1->Branch("ntrover",&ntuple1_st.ntrover,"ntrover/I");
6b8f55ce	111
	112	ntuple2 = new TTree("ntuple2","Demo ntuple2");
	113	ntuple2->Branch("nxP",&ntuple2_st.nxP,"nxP/I");
	114	ntuple2->Branch("nxN",&ntuple2_st.nxN,"nxN/I");
	115	ntuple2->Branch("x",&ntuple2_st.x,"x/F");
	116	ntuple2->Branch("z",&ntuple2_st.z,"z/F");
	117
	118
	119	// Create Histogramms
	120
	121	TH1F *NxP5 = new TH1F("NxP5","P cluster size for layer 5",20,0.,20.);
	122	TH1F *NxN5 = new TH1F("NxN5","N cluster size for layer 5",20,0.,20.);
	123	TH1F *NxP6 = new TH1F("NxP6","P cluster size for layer 6",20,0.,20.);
	124	TH1F *NxN6 = new TH1F("NxN6","N cluster size for layer 6",20,0.,20.);
	125
	126	TH1F *Xres5 = new TH1F("Xres5","Xrec and Xgen difference (micr) for layers 5",100,-200.,200.);
	127	TH1F *Xres6 = new TH1F("Xres6","Xrec and Xgen difference (micr) for layers 6",100,-200.,200.);
	128	TH1F *Zres5 = new TH1F("Zres5","Zrec and Zgen difference (micr) for layers 5",100,-8000.,8000.);
	129	TH1F *Zres6 = new TH1F("Zres6","Zrec and Zgen difference (micr) for layers 6",100,-8000.,8000.);
	130	TH1F *Path5 = new TH1F("Path5","Path length in Si",100,0.,600.);
	131	TH1F *Path6 = new TH1F("Path6","Path length in Si",100,0.,600.);
	132	TH1F *dEdX = new TH1F("dEdX","dEdX (KeV)",100,0.,500.);
	133	TH2F *adcPadcN5all = new TH2F("adcPadcN5all","adcP/N correlation for lay5",100,0.,200.,100,0.,200.);
	134	TH2F *adcPadcN6all = new TH2F("adcPadcN6all","adcP/N correlation for lay6",100,0.,200.,100,0.,200.);
	135	TH2F *adcPadcN5cut = new TH2F("adcPadcN5cut","adcP/N correlation for lay5 and cut of P-N signas",100,0.,200.,100,0.,200.);
	136	TH2F *adcPadcN6cut = new TH2F("adcPadcN6cut","adcP/N correlation for lay6 and cut of P-N signals",100,0.,200.,100,0.,200.);
	137
	138
	139	AliITS ITS = (AliITS) gAlice->GetModule("ITS");
	140	if (!ITS) { cout << "no ITS" << endl; return; }
	141
3f0d03a1	142	//AliITSgeom *aliitsgeo = ITS->GetITSgeom();
3f0d03a1	143	AliITSgeom *geom = ITS->GetITSgeom();
6b8f55ce	144
43c39b72	145
6b8f55ce	146	//Int_t cp[8]={0,0,0,0,0,0,0,0};
6b8f55ce	147
43c39b72	148	//cout << "SSD" << endl;
6b8f55ce	149
	150	AliITSDetType *iDetType=ITS->DetType(2);
	151	AliITSsegmentationSSD seg2=(AliITSsegmentationSSD)iDetType->GetSegmentationModel();
	152	AliITSresponseSSD res2 = (AliITSresponseSSD)iDetType->GetResponseModel();
	153	//res2->SetSigmaSpread(3.,2.);
	154	AliITSsimulationSSD *sim2=new AliITSsimulationSSD(seg2,res2);
	155	ITS->SetSimulationModel(2,sim2);
	156
	157	TClonesArray *dig2 = ITS->DigitsAddress(2);
	158	TClonesArray *recp2 = ITS->ClustersAddress(2);
6b8f55ce	159	AliITSClusterFinderSSD *rec2=new AliITSClusterFinderSSD(seg2,dig2);
6b8f55ce	160	ITS->SetReconstructionModel(2,rec2);
43c39b72	161
6b8f55ce	162	// test
43c39b72	163	printf("SSD dimensions %f %f %f \n",seg2->Dx(),seg2->Dz(),seg2->Dy());
6b8f55ce	164	printf("SSD nstrips %d %d \n",seg2->Npz(),seg2->Npx());
43c39b72	165	Float_t ylim = seg2->Dy()/2 - 12;
6b8f55ce	166
	167	//
	168	// Loop over events
	169	//
	170
	171
	172	Int_t Nh=0;
	173	Int_t Nh1=0;
	174	for (int nev=0; nev<= evNumber2; nev++) {
	175	Int_t nparticles = 0;
	176	nparticles = gAlice->GetEvent(nev);
	177	cout << "nev " << nev <<endl;
	178	cout << "nparticles " << nparticles <<endl;
	179	if (nev < evNumber1) continue;
	180	if (nparticles <= 0) return;
	181
	182	AliITShit *itsHit;
	183	AliITSRecPoint *itsPnt = 0;
	184	AliITSRawClusterSSD *itsClu = 0;
	185
	186	// Get Hit, Cluster & Recpoints Tree Pointers
	187
	188	TTree *TH = gAlice->TreeH();
	189	Int_t nenthit=TH->GetEntries();
	190	printf("Found %d entries in the Hit tree (must be one per track per event!)\n",nenthit);
	191
	192	ITS->GetTreeC(nev);
	193	TTree *TC=ITS->TreeC();
	194	Int_t nentclu=TC->GetEntries();
	195	printf("Found %d entries in the Cluster tree (must be one per module per event!)\n",nentclu);
	196
	197	TTree *TR = gAlice->TreeR();
	198	Int_t nentrec=TR->GetEntries();
	199	printf("Found %d entries in the RecPoints tree\n",nentrec);
	200
	201	// Get Pointers to Clusters & Recpoints TClonesArrays
	202
	203	TClonesArray *ITSclu = ITS->ClustersAddress(2);
	204	printf ("ITSclu %p \n",ITSclu);
	205	TClonesArray *ITSrec = ITS->RecPoints();
	206	printf ("ITSrec %p \n",ITSrec);
	207
	208	// check recpoints
	209
	210	//Int_t nbytes = 0;
	211	Int_t totpoints = 0;
	212	Int_t totclust = 0;
	213
	214	// check hits
	215
	216	Int_t nmodules=0;
3f0d03a1	217	Int_t mod;
6b8f55ce	218
	219	ITS->InitModules(-1,nmodules);
	220	ITS->FillModules(nev,0,nmodules,"","");
	221
	222	TObjArray *fITSmodules = ITS->GetModules();
	223
3f0d03a1	224	Int_t first0 = geom->GetStartDet(0); // SPD
	225	Int_t last0 = geom->GetLastDet(0); // SPD
	226	Int_t first1 = geom->GetStartDet(1); // SDD
	227	Int_t last1 = geom->GetLastDet(1); // SDD
	228	Int_t first2 = geom->GetStartDet(2); // SSD
	229	Int_t last2 = geom->GetLastDet(2); // SSD
	230
	231	// For the SPD: first0 = 0, last0 = 239 (240 modules);
	232	// for the SDD: first1 = 240, last1 = 499 (260 modules);
	233	// for the SSD: first2 = 500, last2 = 2269 (1770 modules).
	234
	235	printf("det type %d first0, last0 %d %d \n",0,first0,last0);
	236	printf("det type %d first1, last1 %d %d \n",1,first1,last1);
	237	printf("det type %d first2, last2 %d %d \n",2,first2,last2);
6b8f55ce	238
3f0d03a1	239	// module loop for the SSD
43c39b72	240	for (mod=first2; mod<last2+1; mod++) { // for the "ALL" option
43c39b72	241	//for (mod=0; mod<last2-first2+1; mod++) { //for the "SSD" option
6b8f55ce	242
3f0d03a1	243	TTree *TR = gAlice->TreeR();
6b8f55ce	244	Int_t nentrec=TR->GetEntries();
	245	//printf("Found %d entries in the RecPoints tree\n",nentrec);
	246
	247	//cout << "CLUSTERS: reset" << endl;
	248	ITS->ResetClusters();
	249	//cout << "CLUSTERS: get" << endl;
	250	TC->GetEvent(mod);
	251	//cout << "RECPOINTS: reset" << endl;
	252	ITS->ResetRecPoints();
	253	//cout << "RECPOINTS: get" << endl;
3f0d03a1	254	//TR->GetEvent(mod+1); // for the V3.04 AliRoot
3f0d03a1	255	TR->GetEvent(mod); // for the V3.05 AliRoot
6b8f55ce	256
	257	Int_t nrecp = ITSrec->GetEntries();
	258	totpoints += nrecp;
43c39b72	259	//if (nrecp) printf("Found %d rec points for module %d\n",nrecp,mod);
43c39b72	260	if (!nrecp) continue;
6b8f55ce	261	Int_t nclusters = ITSclu->GetEntries();
6b8f55ce	262	totclust += nclusters;
6b8f55ce	263	//if (nclusters) printf("Found %d clusters for module %d\n",nrecc,mod);
6b8f55ce	264
487a6bce	265	//AliITSmodule Mod = (AliITSmodule )fITSmodules->At(mod+first2);
0315d466	266	// for the "SSD" option
0315d466	267
487a6bce	268	AliITSmodule Mod = (AliITSmodule )fITSmodules->At(mod);
0315d466	269	// for the "ALL" option
0315d466	270
6b8f55ce	271	// printf("Mod: %X\n",Mod);
	272	Int_t nhits = Mod->GetNhits();
	273	Float_t epart = 0;
43c39b72	274	cout <<" module,nrecp,nclusters,nhits ="<<mod<<","<<nrecp<<","<<nclusters<<","<<nhits<< endl;
6b8f55ce	275
	276	// ---------------- cluster/hit analysis ---------------------
	277
	278
	279	Float_t pathInSSD = 300.;
	280
	281	// ---- Recpoint loop
	282	for (Int_t pnt=0;pnt<nrecp;pnt++) {
	283	itsPnt = (AliITSRecPoint*)ITSrec->At(pnt);
	284	if(!itsPnt) continue;
	285	itsClu = (AliITSRawClusterSSD*)ITSclu->At(pnt);
	286	if(!itsClu) continue;
	287
	288	Int_t nxP = itsClu->fMultiplicity;
	289	Int_t nxN = itsClu->fMultiplicityN;
f239b2fe	290	Int_t ntrover = itsClu->fNtracks;
6b8f55ce	291	Float_t qclP = itsClu->fSignalP; // in ADC
	292	Float_t qclN = itsClu->fSignalN; // in ADC
	293	//Float_t dq = qclP - qclN;
	294	Float_t qcut = itsClu->fQErr; // abs(dq)/signal,
	295	// where signal is
	296	// max of qclP,qclN
	297	Float_t xrec = 10000*itsPnt->GetX();
	298	Float_t zrec = 10000*itsPnt->GetZ();
	299	Float_t qrec = itsPnt->GetQ(); // in ADC, maximum from fSignalP/N
	300	//Float_t dedx = itsPnt->GetdEdX(); // in KeV (ADC * 2.16)
	301	Float_t dedx = itsPnt->fdEdX; // in KeV (ADC * 2.16)
	302	Int_t ii = 0;
	303	Int_t tr1 = itsPnt->GetLabel(ii);
	304	Int_t ii = 1;
	305	Int_t tr2 = itsPnt->GetLabel(ii);
	306	Int_t ii = 2;
	307	Int_t tr3 = itsPnt->GetLabel(ii);
6b8f55ce	308	// fill ntuple2
	309	ntuple2_st.nxP = nxP;
	310	ntuple2_st.nxN = nxN;
	311	ntuple2_st.x = xrec/1000;
	312	ntuple2_st.z = zrec/1000;
	313
	314	if(qcut < 0.18) ntuple2->Fill();
	315
	316
	317	Int_t noverlaps = 0;
	318	Int_t noverprim = 0;
	319	Int_t flaghit = 0;
	320	Float_t xhit0 = 1e+7;
	321	Float_t yhit0 = 1e+7;
	322	Float_t zhit0 = 1e+7;
	323
	324	// Hit loop
	325	for (Int_t hit=0;hit<nhits;hit++) {
	326
	327	itsHit = (AliITShit*)Mod->GetHit(hit);
	328
	329	Int_t flagtrack = 0;
	330	Int_t hitlayer = itsHit->GetLayer();
	331	Int_t hitladder= itsHit->GetLadder();
	332	Int_t hitdet= itsHit->GetDetector();
	333
	334	Int_t track = itsHit->GetTrack();
	335	Int_t dray = 0;
	336	Int_t hitstat = itsHit->GetTrackStatus();
	337
	338	Float_t zhit = 10000*itsHit->GetZL();
	339	Float_t xhit = 10000*itsHit->GetXL();
	340	Float_t yhit = 10000*itsHit->GetYL();
	341	Float_t ehit = 1.0e+6*itsHit->GetIonization(); // hit energy, KeV
	342
	343	Int_t parent = itsHit->GetParticle()->GetFirstMother();
	344	Int_t partcode = itsHit->GetParticle()->GetPdgCode();
	345
	346	// partcode (pdgCode): 11 - e-, 13 - mu-, 22 - gamma, 111 - pi0, 211 - i+
	347	// 310 - K0s, 321 - K+, 2112 - n, 2212 - p, 3122 - lambda
	348
	349	Float_t pmod = itsHit->GetParticle()->P(); // the momentum at the
	350	// vertex
	351	pmod *= 1.0e+3;
	352
43c39b72	353	if(hitstat == 66 && yhit < -ylim) { // entering hit
6b8f55ce	354	xhit0 = xhit;
	355	yhit0 = yhit;
	356	zhit0 = zhit;
	357	}
	358
	359	if(hitstat == 66) continue; // Take the not entering hits only
	360
	361	if(xhit0 > 9e+6 \|\| zhit0 > 9e+6 \|\| yhit0 > 9e+6) {
	362	//cout<<"default xhit0,zhit0,yhit0 ="<<xhit0<<","<<zhit0<<","<<yhit0<<endl;
	363	continue;
	364	}
	365
	366
	367
	368	// Consider the hits only with the track number equaled to one
	369	// of the recpoint
f239b2fe	370	if((track == tr1) \|\| (track == tr2) \|\| (track == tr3)) flagtrack = 1;
6b8f55ce	371
	372	if(flagtrack == 1) { // the hit corresponds to the recpoint
	373
	374	flaghit = 1;
	375
	376	//Float_t px = itsHit->GetPXL(); // the momenta at this GEANT point
	377	//Float_t py = itsHit->GetPYL();
	378	//Float_t pz = itsHit->GetPZL();
	379
	380	Int_t hitprim = 0;
	381
	382	if(partcode == 11 && pmod < 6) dray = 1; // delta ray is e-
	383	// at p < 6 MeV/c
	384
	385	if((hitstat == 68 \|\| hitstat == 33) && dray == 0) noverlaps=noverlaps + 1;
	386	// overlapps for all hits but
	387	// not for delta ray which
	388	// also went out from the
	389	// detector and returned
	390	// again
	391
	392
	393	// x,z resolution colculation
f239b2fe	394	if((hitstat == 68 \|\| hitsat == 33) && dray == 0) {
6b8f55ce	395	Float_t xmed = (xhit + xhit0)/2;
	396	Float_t zmed = (zhit + zhit0)/2;
	397	Float_t xdif = xmed - xrec;
	398	Float_t zdif = zmed - zrec;
	399
0315d466	400	if(parent < 0) {
	401	hitprim = 1; // hitprim=1 for the primery particles
	402	noverprim += 1;
	403	}
6b8f55ce	404	pathInSSD = TMath::Sqrt((xhit0-xhit)(xhit0-xhit)+(yhit0-yhit)(yhit0-yhit)+(zhit0-zhit)*(zhit0-zhit));
6b8f55ce	405
487a6bce	406	// fill ntuple
6b8f55ce	407	ntuple_st.lay = hitlayer;
	408	ntuple_st.nxP = nxP;
	409	ntuple_st.nxN = nxN;
	410	ntuple_st.hitprim = hitprim;
	411	ntuple_st.partcode = partcode;
f239b2fe	412	ntuple_st.ntrover = ntrover;
6b8f55ce	413	ntuple_st.x = xrec/1000;
	414	ntuple_st.z = zrec/1000;
	415	ntuple_st.dx = xdif;
	416	ntuple_st.dz = zdif;
	417	ntuple_st.pmod = pmod;
	418
0315d466	419	//if(qcut < 0.18) ntuple->Fill();
0315d466	420	ntuple->Fill();
6b8f55ce	421
0315d466	422	//if(hitlayer == 5 && qcut < 0.18) {
487a6bce	423
0315d466	424	if(hitlayer == 5 ) {
6b8f55ce	425	Xres5->Fill(xdif);
	426	Zres5->Fill(zdif);
	427	Path5->Fill(pathInSSD);
0315d466	428	}
	429	//if(hitlayer == 6 && qcut < 0.18) {
	430	if(hitlayer == 6) {
6b8f55ce	431	Xres6->Fill(xdif);
	432	Zres6->Fill(zdif);
	433	Path6->Fill(pathInSSD);
	434	}
	435	} // hitstat 68/33
	436	} else { // non correspondent hit
	437	xhit0 = 1e+7;
	438	zhit0 = 1e+7;
	439	} // end of hit-recpoint correspondence
	440	} // hit loop
	441
	442	if(flaghit == 1) {
6b8f55ce	443
	444	if(noverlaps == 0) noverlaps = 1; // cluster contains one or more
	445	// delta rays only
	446
	447	// fill ntuple1
	448	ntuple1_st.lay = hitlayer;
	449	ntuple1_st.lad = hitladder;
	450	ntuple1_st.det = hitdet;
	451	ntuple1_st.nxP = nxP;
	452	ntuple1_st.nxN = nxN;
	453	ntuple1_st.qclP = qclP*300/pathInSSD;
	454	ntuple1_st.qclN = qclN*300/pathInSSD;
	455	ntuple1_st.qrec = qrec*300/pathInSSD;
	456	ntuple1_st.dx = xdif;
	457	ntuple1_st.dz = zdif;
	458	noverlaps -= 1;
	459	noverprim -= 1;
	460	ntuple1_st.noverlaps = noverlaps;
	461	ntuple1_st.noverprim = noverprim;
f239b2fe	462	ntuple1_st.ntrover = ntrover;
6b8f55ce	463
0315d466	464	//if(qcut < 0.18) ntuple1->Fill();
0315d466	465	ntuple1->Fill();
6b8f55ce	466
	467	Float_t de = dedx*300./pathInSSD;
	468	dEdX->Fill(de);
	469	if(hitlayer == 5 ) {
	470	adcPadcN5all->Fill(qclP,qclN);
	471	}
	472	if(hitlayer == 6 ) {
	473	adcPadcN6all->Fill(qclP,qclN);
	474	}
	475	if(hitlayer == 5 && qcut < 0.18) {
	476	adcPadcN5cut->Fill(qclP,qclN);
	477	NxP5->Fill(nxP);
	478	NxN5->Fill(nxN);
	479	}
	480	if(hitlayer == 6 && qcut < 0.18) {
	481	adcPadcN6cut->Fill(qclP,qclN);
	482	NxP6->Fill(nxP);
	483	NxN6->Fill(nxN);
	484	}
	485	} // flaghit = 1
	486	} //b.b. recpoint loop
	487	} //b.b. module loop
	488	} //b.b. evnt loop
	489
	490	TFile fhistos("SSD_his.root","RECREATE");
	491
	492	ntuple->Write();
	493	ntuple1->Write();
	494	ntuple2->Write();
	495	NxP5->Write();
	496	NxN5->Write();
	497	NxP6->Write();
	498	NxN6->Write();
	499	Xres5->Write();
	500	Zres5->Write();
	501	Xres6->Write();
	502	Zres6->Write();
	503	Path5->Write();
	504	Path6->Write();
	505	adcPadcN5all->Write();
	506	adcPadcN6all->Write();
	507	adcPadcN5cut->Write();
	508	adcPadcN6cut->Write();
	509	dEdX->Write();
	510
	511	fhistos.Close();
	512
	513	cout<<"!!! Histogramms and ntuples were written"<<endl;
	514
	515	TCanvas *c1 = new TCanvas("c1","ITS clusters",400,10,600,700);
	516	c1->Divide(2,2);
	517	c1->cd(1);
	518	gPad->SetFillColor(33);
	519	Xres5->SetFillColor(42);
	520	Xres5->Draw();
	521	c1->cd(2);
	522	gPad->SetFillColor(33);
	523	Zres5->SetFillColor(46);
	524	Zres5->Draw();
	525	c1->cd(3);
	526	gPad->SetFillColor(33);
	527	Xres6->SetFillColor(42);
	528	Xres6->Draw();
	529	c1->cd(4);
530	gPad->SetFillColor(33);
531	Zres6->SetFillColor(46);
532	Zres6->Draw();
533
534	cout<<"END test for clusters and hits "<<endl;
535
536	}
537
538
539
540
541
542